Interpolate jx_beam and jy_beam to level 1 (#936)

Hi-PACE · May 7, 2023 · de2affa · de2affa
1 parent 0fec6e1
commit de2affa
Show file tree

Hide file tree

Showing 3 changed files with 54 additions and 1 deletion.
diff --git a/src/Hipace.cpp b/src/Hipace.cpp
@@ -524,6 +524,15 @@ Hipace::SolveOneSlice (int islice, const int islice_local, int step)
             if (islice < m_3D_geom[lev].Domain().smallEnd(Direction::z) ||
                 islice > m_3D_geom[lev].Domain().bigEnd(Direction::z)) {
                 continue;
+            } else if (islice == m_3D_geom[lev].Domain().bigEnd(Direction::z)) {
+                // first slice of level 1 (islice goes backwards)
+                // iterpolate jx_beam and jy_beam from level 0 to level 1
+                m_fields.LevelUp(m_3D_geom, lev, WhichSlice::Previous1, "jx_beam");
+                m_fields.LevelUp(m_3D_geom, lev, WhichSlice::Previous1, "jy_beam");
+                m_fields.LevelUp(m_3D_geom, lev, WhichSlice::This, "jx_beam");
+                m_fields.LevelUp(m_3D_geom, lev, WhichSlice::This, "jy_beam");
+                m_fields.duplicate(lev, WhichSlice::This, {"jx"     , "jy"     },
+                                        WhichSlice::This, {"jx_beam", "jy_beam"});
             }
         }
 

diff --git a/src/fields/Fields.H b/src/fields/Fields.H
@@ -183,14 +183,24 @@ public:
      *
      * \param[in] geom Geometry
      * \param[in] lev current level
-     * \param[in] component which can be Psi or rho
+     * \param[in] component which can be Psi or rho etc.
      * \param[in] outer_edge start writing interpolated values at domain + outer_edge
      * \param[in] inner_edge stop writing interpolated values at domain + inner_edge
      */
     void InterpolateFromLev0toLev1 (amrex::Vector<amrex::Geometry> const& geom, const int lev,
                                     std::string component, const amrex::IntVect outer_edge,
                                     const amrex::IntVect inner_edge);
 
+    /** \brief Interpolate the full field from the coarse grid (lev-1) to the fine grid (lev).
+     *
+     * \param[in] geom Geometry
+     * \param[in] lev current level
+     * \param[in] component which can be jx_beam or jy_beam etc.
+     * \param[in] which_slice slice of the field to interpolate
+     */
+    void LevelUp (amrex::Vector<amrex::Geometry> const& geom, const int lev,
+                  const int which_slice, const std::string component);
+
     /** \brief Compute ExmBy and EypBx on the slice container from J by solving a Poisson equation
      * ExmBy and EypBx are solved in the same function because both rely on Psi.
      *

diff --git a/src/fields/Fields.cpp b/src/fields/Fields.cpp
@@ -695,6 +695,40 @@ Fields::InterpolateFromLev0toLev1 (amrex::Vector<amrex::Geometry> const& geom, c
 }
 
 
+void
+Fields::LevelUp (amrex::Vector<amrex::Geometry> const& geom, const int lev,
+                 const int which_slice, const std::string component)
+{
+    if (lev == 0) return; // only interpolate field to lev 1
+    HIPACE_PROFILE("Fields::LevelUp()");
+    constexpr int interp_order = 2;
+
+    auto field_coarse_interp = interpolated_field_xy<interp_order, amrex::MultiFab>{
+        getField(lev-1, which_slice, component), geom[lev-1]};
+    amrex::MultiFab field_fine = getField(lev, which_slice, component);
+
+    for (amrex::MFIter mfi( field_fine, DfltMfi); mfi.isValid(); ++mfi)
+    {
+        auto arr_field_coarse_interp = field_coarse_interp.array(mfi);
+        const Array2<amrex::Real> arr_field_fine = field_fine.array(mfi);
+
+        const amrex::Real dx = geom[lev].CellSize(0);
+        const amrex::Real dy = geom[lev].CellSize(1);
+        const amrex::Real offset0 = GetPosOffset(0, geom[lev], geom[lev].Domain());
+        const amrex::Real offset1 = GetPosOffset(1, geom[lev], geom[lev].Domain());
+
+        amrex::ParallelFor(field_fine[mfi].box(),
+            [=] AMREX_GPU_DEVICE (int i, int j , int) noexcept
+            {
+                // interpolate the full field
+                const amrex::Real x = i * dx + offset0;
+                const amrex::Real y = j * dy + offset1;
+                arr_field_fine(i,j) = arr_field_coarse_interp(x,y);
+            });
+    }
+}
+
+
 void
 Fields::SolvePoissonExmByAndEypBx (amrex::Vector<amrex::Geometry> const& geom, const int lev)
 {