2D

WeiqunZhang · Jan 19, 2024 · 1ff845b · 1ff845b
1 parent 8618d96
commit 1ff845b
Show file tree

Hide file tree

Showing 9 changed files with 132 additions and 31 deletions.
diff --git a/Src/LinearSolvers/CMakeLists.txt b/Src/LinearSolvers/CMakeLists.txt
@@ -67,7 +67,7 @@ foreach(D IN LISTS AMReX_SPACEDIM)
           )
     endif ()
 
-    if (D EQUAL 3)
+    if (NOT D EQUAL 1)
        target_sources(amrex_${D}d
           PRIVATE
           MLMG/AMReX_MLCurlCurl.H

diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H
@@ -9,17 +9,17 @@ namespace amrex {
 /**
  * \brief curl (alpha curl E) + beta E = rhs
  *
- * Here E is an Array of MultiFabs on staggered grid, alpha is a positive
+ * Here E is an Array of 3 MultiFabs on staggered grid, alpha is a positive
  * scalar, and beta is a non-negative scalar.
  *
  * TODO: If beta is zero, the system could be singular.
- * TODO: Make sure isCellCentered in MLMG makes sense.
+ * TODO: Try different restriction & interpolation strategies.
  */
 class MLCurlCurl
-    : public MLLinOpT<Array<MultiFab,AMREX_SPACEDIM> >
+    : public MLLinOpT<Array<MultiFab,3> >
 {
 public:
-    using MF        = Array<MultiFab,AMREX_SPACEDIM>;
+    using MF        = Array<MultiFab,3>;
     using RT        = typename MLLinOpT<MF>::RT;
     using BCType    = typename MLLinOpT<MF>::BCType;
     using BCMode    = typename MLLinOpT<MF>::BCMode;
@@ -104,6 +104,10 @@ public:
 
     void avgDownResMG (int clev, MF& cres, MF const& fres) const override;
 
+    [[nodiscard]] IntVect getNGrowVectRestriction () const override {
+        return IntVect(0);
+    }
+
     void make (Vector<Vector<MF> >& mf, IntVect const& ng) const override;
 
     [[nodiscard]] MF make (int amrlev, int mglev, IntVect const& ng) const override;
@@ -134,10 +138,15 @@ private:
 
     RT m_alpha = std::numeric_limits<RT>::lowest();
     RT m_beta  = std::numeric_limits<RT>::lowest();
-    Array<IntVect,AMREX_SPACEDIM> m_etype{IntVect(0,1,1),
-                                          IntVect(1,0,1),
-                                          IntVect(1,1,0)};
-    mutable Vector<Vector<Array<std::unique_ptr<iMultiFab>,AMREX_SPACEDIM>>> m_dotmask;
+
+    Array<IntVect,3> m_etype
+#if (AMREX_SPACEDIM == 3)
+         {IntVect(0,1,1), IntVect(1,0,1), IntVect(1,1,0)};
+#else
+         {IntVect(0,1), IntVect(1,0), IntVect(1,1)};
+#endif
+
+    mutable Vector<Vector<Array<std::unique_ptr<iMultiFab>,3>>> m_dotmask;
     static constexpr int m_ncomp = 1;
 };
 

diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.cpp
@@ -44,6 +44,9 @@ void MLCurlCurl::restriction (int amrlev, int cmglev, MF& crse, MF& fine) const
     for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
         amrex::average_down_edges(fine[idim], crse[idim], ratio);
     }
+#if (AMREX_SPACEDIM == 2)
+    amrex::average_down_nodal(fine[2], crse[2], ratio);
+#endif
 }
 
 void MLCurlCurl::interpolation (int amrlev, int fmglev, MF& fine,
@@ -52,7 +55,7 @@ void MLCurlCurl::interpolation (int amrlev, int fmglev, MF& fine,
     IntVect ratio = (amrlev > 0) ? IntVect(2) : this->mg_coarsen_ratio_vec[fmglev];
     AMREX_ALWAYS_ASSERT(ratio == 2);
 
-    for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+    for (int idim = 0; idim < 3; ++idim) {
         auto const& finema = fine[idim].arrays();
         auto const& crsema = crse[idim].const_arrays();
         ParallelFor(fine[idim], [=] AMREX_GPU_DEVICE (int bno, int i, int j, int k)
@@ -166,11 +169,13 @@ void MLCurlCurl::smooth (int amrlev, int mglev, MF& sol, const MF& rhs,
     applyBC(amrlev, mglev, sol[0]);
 
     smooth(amrlev, mglev, sol, rhs[1], 1); // Ey black
+#if (AMREX_SPACEDIM == 3)
     applyBC(amrlev, mglev, sol[1]);
+#endif
 
     smooth(amrlev, mglev, sol, rhs[2], 1); // Ez black
 
-    for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+    for (int idim = 0; idim < 3; ++idim) {
         amrex::OverrideSync(sol[idim], getDotMask(amrlev,mglev,idim),
                             this->m_geom[amrlev][mglev].periodicity());
     }
@@ -338,7 +343,7 @@ Real MLCurlCurl::xdoty (int amrlev, int mglev, const MF& x, const MF& y,
                         bool local) const
 {
     auto result = Real(0.0);
-    for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+    for (int idim = 0; idim < 3; ++idim) {
         auto rtmp = MultiFab::Dot(getDotMask(amrlev,mglev,idim),
                                   x[idim], 0, y[idim], 0, 1, 0, false);
         result += rtmp;
@@ -360,7 +365,7 @@ void MLCurlCurl::averageDownAndSync (Vector<MF>& sol) const
     AMREX_ALWAYS_ASSERT(sol.size() == 1);
     const int amrlev = 0;
     const int mglev = 0;
-    for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+    for (int idim = 0; idim < 3; ++idim) {
         amrex::OverrideSync(sol[amrlev][idim], getDotMask(amrlev,mglev,idim),
                             this->m_geom[amrlev][mglev].periodicity());
     }
@@ -383,52 +388,52 @@ void MLCurlCurl::make (Vector<Vector<MF> >& mf, IntVect const& ng) const
     MLLinOpT<MF>::make(mf, ng);
 }
 
-Array<MultiFab,AMREX_SPACEDIM>
+Array<MultiFab,3>
 MLCurlCurl::make (int amrlev, int mglev, IntVect const& ng) const
 {
     MF r;
-    for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+    for (int idim = 0; idim < 3; ++idim) {
         r[idim].define(amrex::convert(this->m_grids[amrlev][mglev], m_etype[idim]),
                        this->m_dmap[amrlev][mglev], m_ncomp, ng, MFInfo(),
                        *(this->m_factory)[amrlev][mglev]);
     }
     return r;
 }
 
-Array<MultiFab,AMREX_SPACEDIM>
+Array<MultiFab,3>
 MLCurlCurl::makeAlias (MF const& mf) const
 {
     MF r;
-    for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+    for (int idim = 0; idim < 3; ++idim) {
         r[idim] = MultiFab(mf[idim], amrex::make_alias, 0, mf[idim].nComp());
     }
     return r;
 }
 
-Array<MultiFab,AMREX_SPACEDIM>
+Array<MultiFab,3>
 MLCurlCurl::makeCoarseMG (int amrlev, int mglev, IntVect const& ng) const
 {
     BoxArray cba = this->m_grids[amrlev][mglev];
     IntVect ratio = (amrlev > 0) ? IntVect(2) : this->mg_coarsen_ratio_vec[mglev];
     cba.coarsen(ratio);
 
     MF r;
-    for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+    for (int idim = 0; idim < 3; ++idim) {
         r[idim].define(amrex::convert(cba, m_etype[idim]),
                        this->m_dmap[amrlev][mglev], m_ncomp, ng);
     }
     return r;
 }
 
-Array<MultiFab,AMREX_SPACEDIM>
+Array<MultiFab,3>
 MLCurlCurl::makeCoarseAmr (int famrlev, IntVect const& ng) const
 {
     BoxArray cba = this->m_grids[famrlev][0];
     IntVect ratio(this->AMRRefRatio(famrlev-1));
     cba.coarsen(ratio);
 
     MF r;
-    for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+    for (int idim = 0; idim < 3; ++idim) {
         r[idim].define(amrex::convert(cba, m_etype[idim]),
                        this->m_dmap[famrlev][0], m_ncomp, ng);
     }
@@ -439,9 +444,9 @@ void MLCurlCurl::applyBC (int amrlev, int mglev, MF& in) const
 {
     Vector<MultiFab*> mfs{AMREX_D_DECL(in.data(),&(in[1]),&(in[2]))};
     FillBoundary(mfs, this->m_geom[amrlev][mglev].periodicity());
-    applyPhysBC(amrlev, mglev, in[0]);
-    applyPhysBC(amrlev, mglev, in[1]);
-    applyPhysBC(amrlev, mglev, in[2]);
+    for (auto& mf : in) {
+        applyPhysBC(amrlev, mglev, mf);
+    }
 }
 
 void MLCurlCurl::applyBC (int amrlev, int mglev, MultiFab& mf) const
@@ -474,7 +479,7 @@ void MLCurlCurl::applyPhysBC (int amrlev, int mglev, MultiFab& mf) const
 #endif
 
 #ifdef AMREX_USE_OMP
-#pragma omp parallel if (Gpu::notInLaunchRegion());
+#pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
     for (MFIter mfi(mf,mfi_info); mfi.isValid(); ++mfi) {
         auto const& vbx = mfi.validbox();
@@ -523,6 +528,12 @@ iMultiFab const& MLCurlCurl::getDotMask (int amrlev, int mglev, int idim) const
 
 int MLCurlCurl::getDirichlet (int amrlev, int mglev, int idim, int face) const
 {
+#if (AMREX_SPACEDIM == 2)
+    if (idim == 2) {
+        return std::numeric_limits<int>::lowest();
+    }
+#endif
+
     if (face == 0) {
         if (m_lobc[0][idim] == LinOpBCType::Dirichlet) {
             return m_geom[amrlev][mglev].Domain().smallEnd(idim);

diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl_K.H
@@ -23,6 +23,18 @@ void mlcurlcurl_adotx_x (int i, int j, int k, Array4<Real> const& Ax,
                          Real alpha, Real beta,
                          GpuArray<Real,AMREX_SPACEDIM> const& dxinv)
 {
+#if (AMREX_SPACEDIM == 2)
+    amrex::ignore_unused(ez);
+    Real dyy = dxinv[1] * dxinv[1];
+    Real dxy = dxinv[0] * dxinv[1];
+    Real ccex =  ex(i  ,j  ,k  ) * dyy * Real(2.0)
+        - dyy * (ex(i  ,j-1,k  ) +
+                 ex(i  ,j+1,k  ))
+        + dxy * (ey(i  ,j-1,k  )
+               - ey(i  ,j  ,k  )
+               - ey(i+1,j-1,k  )
+               + ey(i+1,j  ,k  ));
+#else
     Real dyy = dxinv[1] * dxinv[1];
     Real dzz = dxinv[2] * dxinv[2];
     Real dxy = dxinv[0] * dxinv[1];
@@ -40,6 +52,7 @@ void mlcurlcurl_adotx_x (int i, int j, int k, Array4<Real> const& Ax,
                - ez(i  ,j  ,k  )
                - ez(i+1,j  ,k-1)
                + ez(i+1,j  ,k  ));
+#endif
     Ax(i,j,k) = alpha * ccex + beta * ex(i,j,k);
 }
 
@@ -51,6 +64,18 @@ void mlcurlcurl_adotx_y (int i, int j, int k, Array4<Real> const& Ay,
                          Real alpha, Real beta,
                          GpuArray<Real,AMREX_SPACEDIM> const& dxinv)
 {
+#if (AMREX_SPACEDIM == 2)
+    amrex::ignore_unused(ez);
+    Real dxx = dxinv[0] * dxinv[0];
+    Real dxy = dxinv[0] * dxinv[1];
+    Real ccey =  ey(i  ,j  ,k  ) * dxx * Real(2.0)
+        - dxx * (ey(i-1,j  ,k  ) +
+                 ey(i+1,j  ,k  ))
+        + dxy * (ex(i-1,j  ,k  )
+               - ex(i  ,j  ,k  )
+               - ex(i-1,j+1,k  )
+               + ex(i  ,j+1,k  ));
+#else
     Real dxx = dxinv[0] * dxinv[0];
     Real dzz = dxinv[2] * dxinv[2];
     Real dxy = dxinv[0] * dxinv[1];
@@ -68,6 +93,7 @@ void mlcurlcurl_adotx_y (int i, int j, int k, Array4<Real> const& Ay,
                - ez(i  ,j  ,k  )
                - ez(i  ,j+1,k-1)
                + ez(i  ,j+1,k  ));
+#endif
     Ay(i,j,k) = alpha * ccey + beta * ey(i,j,k);
 }
 
@@ -79,6 +105,16 @@ void mlcurlcurl_adotx_z (int i, int j, int k, Array4<Real> const& Az,
                          Real alpha, Real beta,
                          GpuArray<Real,AMREX_SPACEDIM> const& dxinv)
 {
+#if (AMREX_SPACEDIM == 2)
+    amrex::ignore_unused(ex,ey);
+    Real dxx = dxinv[0] * dxinv[0];
+    Real dyy = dxinv[1] * dxinv[1];
+    Real ccez =  ez(i  ,j  ,k  ) * (dxx+dyy)*Real(2.0)
+        - dxx * (ez(i-1,j  ,k  ) +
+                 ez(i+1,j  ,k  ))
+        - dyy * (ez(i  ,j-1,k  ) +
+                 ez(i  ,j+1,k  ));
+#else
     Real dxx = dxinv[0] * dxinv[0];
     Real dyy = dxinv[1] * dxinv[1];
     Real dxz = dxinv[0] * dxinv[2];
@@ -96,6 +132,7 @@ void mlcurlcurl_adotx_z (int i, int j, int k, Array4<Real> const& Az,
                - ey(i  ,j  ,k  )
                - ey(i  ,j-1,k+1)
                + ey(i  ,j  ,k+1));
+#endif
     Az(i,j,k) = alpha * ccez + beta * ez(i,j,k);
 }
 
@@ -112,6 +149,19 @@ void mlcurlcurl_gsrb_x (int i, int j, int k,
 
     if ((i+j+k+redblack) % 2 != 0) { return; }
 
+#if (AMREX_SPACEDIM == 2)
+    amrex::ignore_unused(ez);
+    Real dyy = dxinv[1] * dxinv[1];
+    Real dxy = dxinv[0] * dxinv[1];
+    Real gamma = alpha * (dyy)*Real(2.0) + beta;
+    Real ccex =
+        - dyy * (ex(i  ,j-1,k  ) +
+                 ex(i  ,j+1,k  ))
+        + dxy * (ey(i  ,j-1,k  )
+               - ey(i  ,j  ,k  )
+               - ey(i+1,j-1,k  )
+               + ey(i+1,j  ,k  ));
+#else
     Real dyy = dxinv[1] * dxinv[1];
     Real dzz = dxinv[2] * dxinv[2];
     Real dxy = dxinv[0] * dxinv[1];
@@ -130,6 +180,7 @@ void mlcurlcurl_gsrb_x (int i, int j, int k,
                - ez(i  ,j  ,k  )
                - ez(i+1,j  ,k-1)
                + ez(i+1,j  ,k  ));
+#endif
     Real res = rhs(i,j,k) - (gamma*ex(i,j,k) + alpha*ccex);
     ex(i,j,k) += omega/gamma * res;
 }
@@ -147,6 +198,19 @@ void mlcurlcurl_gsrb_y (int i, int j, int k,
 
     if ((i+j+k+redblack) % 2 != 0) { return; }
 
+#if (AMREX_SPACEDIM == 2)
+    amrex::ignore_unused(ez);
+    Real dxx = dxinv[0] * dxinv[0];
+    Real dxy = dxinv[0] * dxinv[1];
+    Real gamma = alpha * (dxx)*Real(2.0) + beta;
+    Real ccey =
+        - dxx * (ey(i-1,j  ,k  ) +
+                 ey(i+1,j  ,k  ))
+        + dxy * (ex(i-1,j  ,k  )
+               - ex(i  ,j  ,k  )
+               - ex(i-1,j+1,k  )
+               + ex(i  ,j+1,k  ));
+#else
     Real dxx = dxinv[0] * dxinv[0];
     Real dzz = dxinv[2] * dxinv[2];
     Real dxy = dxinv[0] * dxinv[1];
@@ -165,6 +229,7 @@ void mlcurlcurl_gsrb_y (int i, int j, int k,
                - ez(i  ,j  ,k  )
                - ez(i  ,j+1,k-1)
                + ez(i  ,j+1,k  ));
+#endif
     Real res = rhs(i,j,k) - (gamma*ey(i,j,k) + alpha*ccey);
     ey(i,j,k) += omega/gamma * res;
 }
@@ -182,6 +247,17 @@ void mlcurlcurl_gsrb_z (int i, int j, int k,
 
     if ((i+j+k+redblack) % 2 != 0) { return; }
 
+#if (AMREX_SPACEDIM == 2)
+    amrex::ignore_unused(ex,ey);
+    Real dxx = dxinv[0] * dxinv[0];
+    Real dyy = dxinv[1] * dxinv[1];
+    Real gamma = alpha * (dxx+dyy)*Real(2.0) + beta;
+    Real ccez =
+        - dxx * (ez(i-1,j  ,k  ) +
+                 ez(i+1,j  ,k  ))
+        - dyy * (ez(i  ,j-1,k  ) +
+                 ez(i  ,j+1,k  ));
+#else
     Real dxx = dxinv[0] * dxinv[0];
     Real dyy = dxinv[1] * dxinv[1];
     Real dxz = dxinv[0] * dxinv[2];
@@ -200,6 +276,7 @@ void mlcurlcurl_gsrb_z (int i, int j, int k,
                - ey(i  ,j  ,k  )
                - ey(i  ,j-1,k+1)
                + ey(i  ,j  ,k+1));
+#endif
     Real res = rhs(i,j,k) - (gamma*ez(i,j,k) + alpha*ccez);
     ez(i,j,k) += omega/gamma * res;
 }

diff --git a/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H b/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H
@@ -617,6 +617,10 @@ protected:
 
     [[nodiscard]] bool isCellCentered () const noexcept { return m_ixtype == 0; }
 
+    [[nodiscard]] virtual IntVect getNGrowVectRestriction () const {
+        return isCellCentered() ? IntVect(0) : IntVect(1);
+    }
+
     virtual void make (Vector<Vector<MF> >& mf, IntVect const& ng) const;
 
     [[nodiscard]] virtual MF make (int amrlev, int mglev, IntVect const& ng) const;

diff --git a/Src/LinearSolvers/MLMG/AMReX_MLMG.H b/Src/LinearSolvers/MLMG/AMReX_MLMG.H
@@ -1020,7 +1020,7 @@ MLMGT<MF>::prepareForSolve (Vector<AMF*> const& a_sol, Vector<AMF const*> const&
         makeSolvable();
     }
 
-    IntVect ng = linop.isCellCentered() ? IntVect(0) : IntVect(1);
+    IntVect ng = linop.getNGrowVectRestriction();
     if (cf_strategy == CFStrategy::ghostnodes) { ng = ng_rhs; }
     if (!solve_called) {
         linop.make(res, ng);

diff --git a/Src/LinearSolvers/MLMG/Make.package b/Src/LinearSolvers/MLMG/Make.package
@@ -87,7 +87,7 @@ ifneq ($(BL_NO_FORT),TRUE)
   F90EXE_sources += AMReX_MLLinOp_nd.F90
 endif
 
-ifeq ($(DIM),3)
+ifneq ($(DIM),1)
   CEXE_headers += AMReX_MLCurlCurl.H
   CEXE_sources += AMReX_MLCurlCurl.cpp
   CEXE_headers += AMReX_MLCurlCurl_K.H