Initial version of use_damp_fields_in_z_guard option (#1304)

* Initial version of use_damp_in_z_guard option

* Fix typo in WarpX.H

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Changed use_damp_in_z_guard to use_damp_fields_in_z_guard

* Added documentatin for use_damp_fields_in_z_guard

* Clean up for use_damp_fields_in_z_guard

* Added use_damp_fields_in_z_guard to picmi interface, and moved it to the psatd input group

* For RZ spectral, damping of the fields is on by default when not periodic

* Updated documentation for use_damp_fields_in_z_guard

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

1 files changed, 88 insertions, 0 deletions

diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp
index a2789996e..1b38d798e 100644
--- a/Source/FieldSolver/WarpXPushFieldsEM.cpp
+++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp
@@ -135,6 +135,9 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */)
         PushPSATDSinglePatch( *spectral_solver_cp[lev],
              Efield_cp[lev], Bfield_cp[lev], Efield_avg_cp[lev], Bfield_avg_cp[lev], current_cp[lev], rho_cp[lev] );
     }
+    if (use_damp_fields_in_z_guard) {
+        DampFieldsInGuards( Efield_fp[lev], Bfield_fp[lev] );
+    }
 }
 #endif
 
@@ -328,6 +331,91 @@ WarpX::MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real a_dt) {
     }
 }
 
+void
+WarpX::DampFieldsInGuards(std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
+                          std::array<std::unique_ptr<amrex::MultiFab>,3>& Bfield) {
+
+    constexpr int zdir = (AMREX_SPACEDIM - 1);
+
+    for ( amrex::MFIter mfi(*Efield[0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
+    {
+
+        amrex::Array4<amrex::Real> const& Ex_arr = Efield[0]->array(mfi);
+        amrex::Array4<amrex::Real> const& Ey_arr = Efield[1]->array(mfi);
+        amrex::Array4<amrex::Real> const& Ez_arr = Efield[2]->array(mfi);
+        amrex::Array4<amrex::Real> const& Bx_arr = Bfield[0]->array(mfi);
+        amrex::Array4<amrex::Real> const& By_arr = Bfield[1]->array(mfi);
+        amrex::Array4<amrex::Real> const& Bz_arr = Bfield[2]->array(mfi);
+
+        // Get the tilebox from Efield so that it includes the guard cells.
+        amrex::Box tilebox = (*Efield[0])[mfi].box();
+        int const nz_tile = tilebox.bigEnd(zdir);
+
+        // Box for the whole simulation domain
+        amrex::Box const& domain = Geom(0).Domain();
+        int const nz_domain = domain.bigEnd(zdir);
+
+        if (tilebox.smallEnd(zdir) < 0) {
+
+            // Apply damping factor in guards cells below the lower end of the domain
+            int const nz_guard = -tilebox.smallEnd(zdir);
+
+            // Set so the box only covers the lower half of the guard cells
+            tilebox.setBig(zdir, -nz_guard/2-1);
+
+            amrex::ParallelFor(tilebox, Efield[0]->nComp(),
+            [=] AMREX_GPU_DEVICE (int i, int j, int k, int icomp)
+            {
+#if (AMREX_SPACEDIM == 3)
+                amrex::Real zcell = static_cast<amrex::Real>(k + nz_guard);
+#else
+                amrex::Real zcell = static_cast<amrex::Real>(j + nz_guard);
+#endif
+                amrex::Real phase = MathConst::pi*zcell/nz_guard;
+                amrex::Real damp_factor = std::pow(std::sin(phase), 2);
+
+                Ex_arr(i,j,k,icomp) *= damp_factor;
+                Ey_arr(i,j,k,icomp) *= damp_factor;
+                Ez_arr(i,j,k,icomp) *= damp_factor;
+                Bx_arr(i,j,k,icomp) *= damp_factor;
+                By_arr(i,j,k,icomp) *= damp_factor;
+                Bz_arr(i,j,k,icomp) *= damp_factor;
+
+            });
+
+        }
+        else if (nz_tile > nz_domain) {
+
+            // Apply damping factor in guards cells above the upper end of the domain
+            int nz_guard = nz_tile - nz_domain;
+
+            // Set so the box only covers the upper half of the guard cells
+            tilebox.setSmall(zdir, nz_domain + nz_guard/2 + 1);
+
+            amrex::ParallelFor(tilebox, Efield[0]->nComp(),
+            [=] AMREX_GPU_DEVICE (int i, int j, int k, int icomp)
+            {
+#if (AMREX_SPACEDIM == 3)
+                amrex::Real zcell = static_cast<amrex::Real>(nz_tile - k);
+#else
+                amrex::Real zcell = static_cast<amrex::Real>(nz_tile - j);
+#endif
+                amrex::Real phase = MathConst::pi*zcell/nz_guard;
+                amrex::Real damp_factor = std::pow(std::sin(phase), 2);
+
+                Ex_arr(i,j,k,icomp) *= damp_factor;
+                Ey_arr(i,j,k,icomp) *= damp_factor;
+                Ez_arr(i,j,k,icomp) *= damp_factor;
+                Bx_arr(i,j,k,icomp) *= damp_factor;
+                By_arr(i,j,k,icomp) *= damp_factor;
+                Bz_arr(i,j,k,icomp) *= damp_factor;
+
+            });
+
+        }
+    }
+}
+
 #ifdef WARPX_DIM_RZ
 // This scales the current by the inverse volume and wraps around the depostion at negative radius.
 // It is faster to apply this on the grid than to do it particle by particle.