symmetrize fields (#1038)

docs/source/run/parameters.rst

@@ -247,6 +247,10 @@ The default is to use the explicit solver. **We strongly recommend to use the ex
     ``fields.extended_solve = true`` and ``geometry.is_periodic = false false false``.
     Only available with the predictor-corrector solver.
 
+* ``fields.do_symmetrize`` (`bool`) optional (default `0`)
+    Symmetrizes current and charge densities transversely before the field solve.
+    Each cell at (`x`, `y`) is averaged with cells at (`-x`, `y`), (`x`, `-y`) and (`-x`, `-y`).
+
 Explicit solver parameters
 ^^^^^^^^^^^^^^^^^^^^^^^^^^
 

src/Hipace.cpp

@@ -731,6 +731,12 @@ Hipace::ExplicitMGSolveBxBy (const int lev, const int which_slice)
     m_fields.LevelUpBoundary(m_3D_geom, lev, which_slice, "By",
         m_fields.m_slices_nguards, m_fields.m_poisson_nguards);
 
+    if (m_fields.m_do_symmetrize) {
+        m_fields.SymmetrizeFields(Comps[which_slice_chi]["chi"], lev, 1, 1);
+        m_fields.SymmetrizeFields(Comps[which_slice]["Sx"], lev, -1, 1);
+        m_fields.SymmetrizeFields(Comps[which_slice]["Sy"], lev, 1, -1);
+    }
+
 #ifdef AMREX_USE_LINEAR_SOLVERS
     if (m_use_amrex_mlmg) {
         // Copy chi to chi2

src/fields/Fields.H

@@ -241,6 +241,15 @@ public:
      */
     void SolvePoissonBxBy (amrex::Vector<amrex::Geometry> const& geom, const int current_N_level,
                            const int which_slice);
+    /** \brief Symmetrize fields by averaging over (x,y), symm_x*(-x,y),
+     * symm_y*(x,-y) and symm_x*symm_y*(-x,-y) where symm_x and symm_y can be 1 or -1.
+     *
+     * \param[in] field_comp component index of the filed
+     * \param[in] lev current level
+     * \param[in] symm_x type of reflection in x direction
+     * \param[in] symm_y type of reflection in y direction
+     */
+    void SymmetrizeFields (int field_comp, const int lev, const int symm_x, const int symm_y);
     /** \brief Sets the initial guess of the B field from the two previous slices
      *
      * This modifies component Bx or By of slice 1 in m_fields.m_slices
@@ -474,6 +483,8 @@ public:
     inline static amrex::IntVect m_poisson_nguards {-1, -1, -1};
     /** If the poisson solver should include the guard cells */
     bool m_extended_solve = false;
+    /** Whether the currents should be symmetrized for the field solve */
+    bool m_do_symmetrize = false;
 
 private:
     /** Vector over levels of all required fields to compute current slice */

src/fields/Fields.cpp

@@ -25,6 +25,7 @@ Fields::Fields (const int nlev)
     queryWithParser(ppf, "do_dirichlet_poisson", m_do_dirichlet_poisson);
     queryWithParser(ppf, "extended_solve", m_extended_solve);
     queryWithParser(ppf, "open_boundary", m_open_boundary);
+    queryWithParser(ppf, "do_symmetrize", m_do_symmetrize);
 }
 
 void
@@ -863,6 +864,12 @@ Fields::SolvePoissonPsiExmByEypBxEzBz (amrex::Vector<amrex::Geometry> const& geo
             m_slices_nguards, -m_slices_nguards);
         LevelUpBoundary(geom, lev, WhichSlice::This, "jy",
             m_slices_nguards, -m_slices_nguards);
+
+        if (m_do_symmetrize) {
+            SymmetrizeFields(Comps[WhichSlice::This]["rhomjz"], lev, 1, 1);
+            SymmetrizeFields(Comps[WhichSlice::This]["jx"], lev, -1, 1);
+            SymmetrizeFields(Comps[WhichSlice::This]["jy"], lev, 1, -1);
+        }
     }
 
     for (int lev=0; lev<current_N_level; ++lev) {
@@ -958,6 +965,11 @@ Fields::SolvePoissonEz (amrex::Vector<amrex::Geometry> const& geom,
         // also inside ghost cells to account for x and y derivative
         LevelUpBoundary(geom, lev, which_slice, "jx", m_slices_nguards, -m_slices_nguards);
         LevelUpBoundary(geom, lev, which_slice, "jy", m_slices_nguards, -m_slices_nguards);
+
+        if (m_do_symmetrize) {
+            SymmetrizeFields(Comps[which_slice]["jx"], lev, -1, 1);
+            SymmetrizeFields(Comps[which_slice]["jy"], lev, 1, -1);
+        }
     }
 
     for (int lev=0; lev<current_N_level; ++lev) {
@@ -1006,6 +1018,12 @@ Fields::SolvePoissonBxBy (amrex::Vector<amrex::Geometry> const& geom,
         // also inside ghost cells to account for x and y derivative
         LevelUpBoundary(geom, lev, WhichSlice::This, "jz", m_slices_nguards, -m_slices_nguards);
         // jx and jy on WhichSlice::Previous was already leveled up on previous slice
+
+        if (m_do_symmetrize) {
+            SymmetrizeFields(Comps[WhichSlice::This]["jz"], lev, 1, 1);
+            SymmetrizeFields(Comps[WhichSlice::Next]["jx"], lev, -1, 1);
+            SymmetrizeFields(Comps[WhichSlice::Next]["jy"], lev, 1, -1);
+        }
     }
 
     for (int lev=0; lev<current_N_level; ++lev) {
@@ -1048,6 +1066,42 @@ Fields::SolvePoissonBxBy (amrex::Vector<amrex::Geometry> const& geom,
     }
 }
 
+void
+Fields::SymmetrizeFields (int field_comp, const int lev, const int symm_x, const int symm_y)
+{
+    HIPACE_PROFILE("Fields::SymmetrizeFields()");
+
+    AMREX_ALWAYS_ASSERT(symm_x*symm_x == 1 && symm_y*symm_y == 1);
+
+    amrex::MultiFab& slicemf = getSlices(lev);
+
+    for ( amrex::MFIter mfi(slicemf, DfltMfiTlng); mfi.isValid(); ++mfi ) {
+        const Array2<amrex::Real> arr = slicemf.array(mfi, field_comp);
+
+        const amrex::Box full_box = mfi.growntilebox();
+
+        const int upper_x = full_box.smallEnd(0) + full_box.bigEnd(0);
+        const int upper_y = full_box.smallEnd(1) + full_box.bigEnd(1);
+
+        amrex::Box quarter_box = full_box;
+        quarter_box.setBig(0, full_box.smallEnd(0) + (full_box.length(0)+1)/2 - 1);
+        quarter_box.setBig(1, full_box.smallEnd(1) + (full_box.length(1)+1)/2 - 1);
+
+        amrex::ParallelFor(quarter_box,
+            [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
+            {
+                const amrex::Real avg = 0.25_rt*(arr(i, j) + arr(upper_x - i, j)*symm_x
+                    + arr(i, upper_y - j)*symm_y + arr(upper_x - i, upper_y - j)*symm_x*symm_y);
+
+                // Note: this may write to the same cell multiple times in the center.
+                arr(i, j) = avg;
+                arr(upper_x - i, j) = avg*symm_x;
+                arr(i, upper_y - j) = avg*symm_y;
+                arr(upper_x - i, upper_y - j) = avg*symm_x*symm_y;
+            });
+    }
+}
+
 void
 Fields::InitialBfieldGuess (const amrex::Real relative_Bfield_error,
                             const amrex::Real predcorr_B_error_tolerance, const int lev)