Skip to content

Commit b251fc8

Browse files
JustinRayAngusJustinRayAngus
committed
add blanking for components of E.
1 parent b1c5d6e commit b251fc8

File tree

3 files changed

+54
-10
lines changed

3 files changed

+54
-10
lines changed

‎Source/FieldSolver/ImplicitSolvers/ImplicitSolver.H‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -129,6 +129,10 @@ public:
129129

130130
void ComputeJfromMassMatrices (const bool a_J_from_MM_only);
131131

132+
const amrex::Array<bool,3>& GetBlankElectricField () const {
133+
return m_blank_electric_field;
134+
}
135+
132136
protected:
133137

134138
/**
@@ -148,6 +152,11 @@ protected:
148152
*/
149153
mutable amrex::Real m_dt = 0.0;
150154

155+
/**
156+
* \brief Flags to turn off evolving select components of the electric field
157+
*/
158+
amrex::Array<bool,3> m_blank_electric_field{false,false,false};
159+
151160
/**
152161
* \brief Time-biasing parameter for fields used on RHS to advance system.
153162
*/

‎Source/FieldSolver/ImplicitSolvers/ThetaImplicitEM.cpp‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,13 @@ void ThetaImplicitEM::Define ( WarpX* const a_WarpX )
4848
m_theta>=0.5 && m_theta<=1.0,
4949
"theta parameter for theta implicit time solver must be between 0.5 and 1.0");
5050

51+
amrex::Vector<int> tmp(3);
52+
if (pp.queryarr("blank_electric_field", tmp)) {
53+
for (int dir = 0; dir < 3; ++dir) {
54+
m_blank_electric_field[dir] = (tmp[dir] != 0);
55+
}
56+
}
57+
5158
// Parse nonlinear solver parameters
5259
parseNonlinearSolverParams( pp );
5360

@@ -74,6 +81,9 @@ void ThetaImplicitEM::PrintParameters () const
7481
amrex::Print() << "----------- THETA IMPLICIT EM SOLVER PARAMETERS -----------\n";
7582
amrex::Print() << "-----------------------------------------------------------\n";
7683
amrex::Print() << "Time-bias parameter theta: " << m_theta << "\n";
84+
amrex::Print() << "Blank x-electric field: " << (m_blank_electric_field[0] ? "true":"false") << "\n";
85+
amrex::Print() << "Blank y-electric field: " << (m_blank_electric_field[1] ? "true":"false") << "\n";
86+
amrex::Print() << "Blank z-electric field: " << (m_blank_electric_field[2] ? "true":"false") << "\n";
7787
PrintBaseImplicitSolverParameters();
7888
m_nlsolver->PrintParams();
7989
amrex::Print() << "-----------------------------------------------------------\n\n";
@@ -146,6 +156,16 @@ void ThetaImplicitEM::ComputeRHS ( WarpXSolverVec& a_RHS,
146156

147157
// RHS = cvac^2*m_theta*dt*( curl(Bg^{n+theta}) - mu0*Jg^{n+1/2} )
148158
m_WarpX->ImplicitComputeRHSE( m_theta*m_dt, a_RHS);
159+
160+
// Apply blanking to electric field RHS vector
161+
for (int lev = 0; lev < m_num_amr_levels; ++lev) {
162+
for (int dir = 0; dir < 3; ++dir) {
163+
if (m_blank_electric_field[dir]) {
164+
a_RHS.getArrayVec()[lev][dir]->setVal(0._rt);
165+
}
166+
}
167+
}
168+
149169
}
150170

151171
void ThetaImplicitEM::UpdateWarpXFields ( const WarpXSolverVec& a_E,
@@ -157,6 +177,14 @@ void ThetaImplicitEM::UpdateWarpXFields ( const WarpXSolverVec& a_E,
157177
const amrex::Real theta_time = start_time + m_theta*m_dt;
158178
m_WarpX->SetElectricFieldAndApplyBCs( a_E, theta_time );
159179

180+
// Apply blanking to the electric field vector
181+
for (int lev = 0; lev < m_num_amr_levels; ++lev) {
182+
ablastr::fields::VectorField Efp = m_WarpX->m_fields.get_alldirs(FieldType::Efield_fp, lev);
183+
for (int dir = 0; dir < 3; ++dir) {
184+
if (m_blank_electric_field[dir]) { Efp[dir]->setVal(0._rt); }
185+
}
186+
}
187+
160188
// Update Bfield_fp owned by WarpX
161189
ablastr::fields::MultiLevelVectorField const& B_old = m_WarpX->m_fields.get_mr_levels_alldirs(FieldType::B_old, 0);
162190
m_WarpX->UpdateMagneticFieldAndApplyBCs( B_old, m_theta*m_dt, start_time );

‎Source/NonlinearSolvers/MatrixPC.H‎

Lines changed: 17 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -384,6 +384,8 @@ int MatrixPC<T,Ops>::Assemble (const T& a_U)
384384
Gpu::Buffer<int> nnz_actual_d({nnz_max});
385385
auto* nnz_actual_ptr = nnz_actual_d.data();
386386

387+
const amrex::Array<bool,3>& blank_Efield = m_ops->GetBlankElectricField();
388+
387389
for (int dir = 0; dir < 3; dir++) {
388390

389391
const amrex::MultiFab* BC_mask_Edir = m_ops->GetCurl2BCmask(lev,dir);
@@ -446,6 +448,7 @@ int MatrixPC<T,Ops>::Assemble (const T& a_U)
446448
if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
447449
}
448450

451+
if (!blank_Efield[dir]) {
449452
#if defined(WARPX_DIM_1D_Z)
450453
if (dir != 2) {
451454
{
@@ -484,7 +487,7 @@ int MatrixPC<T,Ops>::Assemble (const T& a_U)
484487
val += 2.0_rt*alpha * dxi[1]*dxi[1] * BC_mask_Edir_arr(i,j,k,0);
485488
} else if (dir == 2) {
486489
val += 2.0_rt*alpha * dxi[0]*dxi[0] * BC_mask_Edir_arr(i,j,k,0);
487-
} else {
490+
} else if (dir == 1) {
488491
val += 2.0_rt*alpha * ( dxi[0]*dxi[0] * BC_mask_Edir_arr(i,j,k,0)
489492
+ dxi[1]*dxi[1] * BC_mask_Edir_arr(i,j,k,2) );
490493
}
@@ -516,43 +519,45 @@ int MatrixPC<T,Ops>::Assemble (const T& a_U)
516519
// The following four blocks are for
517520
// dir = 0: d/dx(dEz/dz) at Ex(i,j) with Ex centered in x and nodal in z
518521
// dir = 2: d/dz(dEx/dx) at Ez(i,j) with Ez centered in z and nodal in x
519-
{
522+
if (!blank_Efield[tdir]) {
523+
{
520524
int cidx_g_rhs = (dir == 0 ? dof_tdir_arr(i,j-1,k,1) : dof_tdir_arr(i-1,j,k,1));
521525
amrex::Real val = alpha * dxi[0] * dxi[1] * BC_mask_Edir_arr(i,j,k,2);
522526
auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
523527
&c_indices_g_ptr[ridx_l*nnz_max],
524528
&a_ij_ptr[ridx_l*nnz_max],
525529
nnz_max, icol );
526530
if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
527-
}
528-
{
531+
}
532+
{
529533
int cidx_g_rhs = dof_tdir_arr(i,j,k,1);
530534
amrex::Real val = -alpha * dxi[0] * dxi[1] * BC_mask_Edir_arr(i,j,k,2);
531535
auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
532536
&c_indices_g_ptr[ridx_l*nnz_max],
533537
&a_ij_ptr[ridx_l*nnz_max],
534538
nnz_max, icol );
535539
if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
536-
}
537-
{
540+
}
541+
{
538542
int cidx_g_rhs = (dir == 0 ? dof_tdir_arr(i+1,j-1,k,1) : dof_tdir_arr(i-1,j+1,k,1));
539543
amrex::Real val = -alpha * dxi[0] * dxi[1] * BC_mask_Edir_arr(i,j,k,2);
540544
auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
541545
&c_indices_g_ptr[ridx_l*nnz_max],
542546
&a_ij_ptr[ridx_l*nnz_max],
543547
nnz_max, icol );
544548
if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
545-
}
546-
{
549+
}
550+
{
547551
int cidx_g_rhs = (dir == 0 ? dof_tdir_arr(i+1,j,k,1) : dof_tdir_arr(i,j+1,k,1));
548552
amrex::Real val = alpha * dxi[0] * dxi[1] * BC_mask_Edir_arr(i,j,k,2);
549553
auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
550554
&c_indices_g_ptr[ridx_l*nnz_max],
551555
&a_ij_ptr[ridx_l*nnz_max],
552556
nnz_max, icol );
553557
if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
558+
}
554559
}
555-
} else {
560+
} else if (dir==1) {
556561
for (int jdir = 0; jdir <= 2; jdir+=2) {
557562
{
558563
int cidx_g_rhs = (jdir == 0 ? dof_arr(i-1,j,k,1) : dof_arr(i,j-1,k,1));
@@ -603,6 +608,7 @@ int MatrixPC<T,Ops>::Assemble (const T& a_U)
603608
for (int ctr_dir = -1; ctr_dir <= 0; ctr_dir++) {
604609
for (int ctr_tdir = -1; ctr_tdir <= 0; ctr_tdir++) {
605610
for (int tdir = 1; tdir <= 2; tdir++) {
611+
if (blank_Efield[dvec[tdir]]) { continue; }
606612
auto iv = ic; iv[dvec[0]] += (ctr_dir+1); iv[dvec[tdir]] += ctr_tdir;
607613
auto sign = std::copysign(1,ctr_dir) * std::copysign(1,ctr_tdir);
608614
int cidx_g_rhs = dof_arrays[tdir](iv,1);
@@ -617,12 +623,13 @@ int MatrixPC<T,Ops>::Assemble (const T& a_U)
617623
}
618624
}
619625
#endif
626+
}
620627
num_nz_ptr[ridx_l] = icol;
621628
});
622629
}
623630

624631
// Add the mass matrix piece
625-
if (m_include_mass_matrices) {
632+
if (m_include_mass_matrices && !blank_Efield[dir]) {
626633

627634
auto sigma_ii_arr = (*m_bcoefs)[lev][dir]->const_array(mfi);
628635
const amrex::IntVect MM_ncomps = m_ops->GetMassMatricesPCNcomp(dir);

0 commit comments

Comments
 (0)