New Linear Solver: Curl of Curl

Add a new linear solver for curl (alpha curl E) + beta E = rhs in 2D & 3D, where E is an array of 3 MultiFabs on edges, and alpha and beta are scalars. It supports periodic, homogeneous Dirichlet, and symmetry boundaries. At the symmetry boundary, the normal component changes the sign, whereas the transverse components do not.
AMReX-Codes · Jan 20, 2024 · 3a8ce86 · 3a8ce86
1 parent 9996b68
commit 3a8ce86
Show file tree

Hide file tree

Showing 23 changed files with 1,517 additions and 28 deletions.
diff --git a/Src/Base/AMReX_FabDataType.H b/Src/Base/AMReX_FabDataType.H
@@ -0,0 +1,27 @@
+#ifndef AMREX_FAB_DATA_TYPE_H_
+#define AMREX_FAB_DATA_TYPE_H_
+#include <AMReX_Config.H>
+
+#include <AMReX_TypeTraits.H>
+
+namespace amrex {
+
+template <typename T, class Enable = void> struct FabDataType {};
+//
+template <typename T>
+struct FabDataType <T, std::enable_if_t<IsMultiFabLike_v<T> > >
+{
+    using   fab_type = typename T::fab_type;
+    using value_type = typename T::value_type;
+};
+
+template <typename T>
+struct FabDataType <T, std::enable_if_t<IsMultiFabLike_v<typename T::value_type> > >
+{
+    using   fab_type = typename T::value_type::fab_type;
+    using value_type = typename T::value_type::value_type;
+};
+
+}
+
+#endif
diff --git a/Src/Base/CMakeLists.txt b/Src/Base/CMakeLists.txt
@@ -132,6 +132,7 @@ foreach(D IN LISTS AMReX_SPACEDIM)
        AMReX_Array4.H
        AMReX_MakeType.H
        AMReX_TypeTraits.H
+       AMReX_FabDataType.H
        AMReX_FabFactory.H
        AMReX_BaseFabUtility.H
        # Fortran data defined on unions of rectangles ----------------------------

diff --git a/Src/Base/Make.package b/Src/Base/Make.package
@@ -149,6 +149,7 @@ C$(AMREX_BASE)_headers += AMReX_IArrayBox.H
 
 C$(AMREX_BASE)_headers += AMReX_MakeType.H
 C$(AMREX_BASE)_headers += AMReX_TypeTraits.H
+C$(AMREX_BASE)_headers += AMReX_FabDataType.H
 
 C$(AMREX_BASE)_headers += AMReX_Array4.H
 C$(AMREX_BASE)_sources += AMReX_BaseFab.cpp

diff --git a/Src/Boundary/AMReX_FabSet.H b/Src/Boundary/AMReX_FabSet.H
@@ -3,6 +3,7 @@
 #define AMREX_FABSET_H_
 #include <AMReX_Config.H>
 
+#include <AMReX_FabDataType.H>
 #include <AMReX_MultiFab.H>
 #include <AMReX_ParallelDescriptor.H>
 #include <AMReX_BLProfiler.H>
@@ -46,8 +47,8 @@ class FabSetT
     friend class FabSetIter;
     friend class FluxRegister;
 public:
-    using value_type = typename MF::value_type;
-    using FAB = typename MF::fab_type;
+    using value_type = typename FabDataType<MF>::value_type;
+    using FAB        = typename FabDataType<MF>::fab_type;
 
     //
     //! The default constructor -- you must later call define().

diff --git a/Src/Boundary/AMReX_LO_BCTYPES.H b/Src/Boundary/AMReX_LO_BCTYPES.H
@@ -20,6 +20,7 @@
 #define AMREX_LO_INFLOW   106
 #define AMREX_LO_INHOMOG_NEUMANN 107
 #define AMREX_LO_ROBIN 108
+#define AMREX_LO_SYMMETRY 109
 #define AMREX_LO_PERIODIC 200
 #define AMREX_LO_BOGUS    1729
 
@@ -38,6 +39,7 @@ namespace amrex {
         inflow           = AMREX_LO_INFLOW,
         inhomogNeumann   = AMREX_LO_INHOMOG_NEUMANN,
         Robin            = AMREX_LO_ROBIN,
+        symmetry         = AMREX_LO_SYMMETRY,
         Periodic         = AMREX_LO_PERIODIC,
         bogus            = AMREX_LO_BOGUS
     };
@@ -48,4 +50,3 @@ namespace amrex {
 #endif
 
 #endif
-
diff --git a/Src/Boundary/AMReX_LO_BCTYPES.cpp b/Src/Boundary/AMReX_LO_BCTYPES.cpp
@@ -52,6 +52,11 @@ std::ostream& operator<< (std::ostream& os, const LinOpBCType& t)
             os << "Robin";
             break;
         }
+        case LinOpBCType::symmetry:
+        {
+            os << "symmetry";
+            break;
+        }
         case LinOpBCType::Periodic:
         {
             os << "Periodic";

diff --git a/Src/LinearSolvers/CMakeLists.txt b/Src/LinearSolvers/CMakeLists.txt
@@ -67,6 +67,15 @@ foreach(D IN LISTS AMReX_SPACEDIM)
           )
     endif ()
 
+    if (NOT D EQUAL 1)
+       target_sources(amrex_${D}d
+          PRIVATE
+          MLMG/AMReX_MLCurlCurl.H
+          MLMG/AMReX_MLCurlCurl.cpp
+          MLMG/AMReX_MLCurlCurl_K.H
+          )
+    endif ()
+
     if (AMReX_EB AND NOT D EQUAL 1)
        target_sources(amrex_${D}d
           PRIVATE

diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H b/Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H
@@ -0,0 +1,155 @@
+#ifndef AMREX_ML_CURL_CURL_H_
+#define AMREX_ML_CURL_CURL_H_
+#include <AMReX_Config.H>
+
+#include <AMReX_MLLinOp.H>
+
+namespace amrex {
+
+/**
+ * \brief curl (alpha curl E) + beta E = rhs
+ *
+ * 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: Try different restriction & interpolation strategies.
+ */
+class MLCurlCurl
+    : public MLLinOpT<Array<MultiFab,3> >
+{
+public:
+    using MF        = Array<MultiFab,3>;
+    using RT        = typename MLLinOpT<MF>::RT;
+    using BCType    = typename MLLinOpT<MF>::BCType;
+    using BCMode    = typename MLLinOpT<MF>::BCMode;
+    using StateMode = typename MLLinOpT<MF>::StateMode;
+    using Location  = typename MLLinOpT<MF>::Location;
+
+    MLCurlCurl () = default;
+    MLCurlCurl (const Vector<Geometry>& a_geom,
+                const Vector<BoxArray>& a_grids,
+                const Vector<DistributionMapping>& a_dmap,
+                const LPInfo& a_info = LPInfo());
+
+    void define (const Vector<Geometry>& a_geom,
+                 const Vector<BoxArray>& a_grids,
+                 const Vector<DistributionMapping>& a_dmap,
+                 const LPInfo& a_info = LPInfo());
+
+    void setScalars (RT a_alpha, RT a_beta) noexcept;
+
+    [[nodiscard]] std::string name () const override {
+        return std::string("curl of curl");
+    }
+
+    void setLevelBC (int amrlev, const MF* levelbcdata,
+                     const MF* robinbc_a = nullptr,
+                     const MF* robinbc_b = nullptr,
+                     const MF* robinbc_f = nullptr) override;
+
+    void restriction (int amrlev, int cmglev, MF& crse, MF& fine) const override;
+
+    void interpolation (int amrlev, int fmglev, MF& fine, const MF& crse) const override;
+
+    void interpAssign (int amrlev, int fmglev, MF& fine, MF& crse) const override;
+
+    void interpolationAmr (int famrlev, MF& fine, const MF& crse,
+                           IntVect const& nghost) const override;
+
+    void averageDownSolutionRHS (int camrlev, MF& crse_sol, MF& crse_rhs,
+                                 const MF& fine_sol, const MF& fine_rhs) override;
+
+    void apply (int amrlev, int mglev, MF& out, MF& in, BCMode bc_mode,
+                StateMode s_mode, const MLMGBndryT<MF>* bndry=nullptr) const override;
+
+    void smooth (int amrlev, int mglev, MF& sol, const MF& rhs,
+                 bool skip_fillboundary=false) const override;
+
+    void solutionResidual (int amrlev, MF& resid, MF& x, const MF& b,
+                           const MF* crse_bcdata=nullptr) override;
+
+    void correctionResidual (int amrlev, int mglev, MF& resid, MF& x,
+                             const MF& b, BCMode bc_mode,
+                             const MF* crse_bcdata=nullptr) override;
+
+    void reflux (int crse_amrlev,
+                 MF& res, const MF& crse_sol, const MF& crse_rhs,
+                 MF& fine_res, MF& fine_sol, const MF& fine_rhs) const override;
+
+    void compFlux (int amrlev, const Array<MF*,AMREX_SPACEDIM>& fluxes,
+                   MF& sol, Location loc) const override;
+
+    void compGrad (int amrlev, const Array<MF*,AMREX_SPACEDIM>& grad,
+                   MF& sol, Location loc) const override;
+
+    void applyMetricTerm (int /*amrlev*/, int /*mglev*/, MF& /*rhs*/) const override {}
+
+    void unapplyMetricTerm (int /*amrlev*/, int /*mglev*/, MF& /*rhs*/) const override {}
+
+    void prepareForSolve () override;
+
+    [[nodiscard]] bool isSingular (int /*amrlev*/) const override { return false; }
+    [[nodiscard]] bool isBottomSingular () const override { return false; }
+
+    RT xdoty (int amrlev, int mglev, const MF& x, const MF& y, bool local) const override;
+
+    std::unique_ptr<MLLinOpT<MF>> makeNLinOp (int /*grid_size*/) const override { return nullptr; }
+
+    [[nodiscard]] RT normInf (int amrlev, MF const& mf, bool local) const override;
+
+    void averageDownAndSync (Vector<MF>& sol) const override;
+
+    void avgDownResAmr (int clev, MF& cres, MF const& fres) const override;
+
+    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;
+
+    [[nodiscard]] MF makeAlias (MF const& mf) const override;
+
+    [[nodiscard]] MF makeCoarseMG (int amrlev, int mglev, IntVect const& ng) const override;
+
+    [[nodiscard]] MF makeCoarseAmr (int famrlev, IntVect const& ng) const override;
+
+// public for cuda
+
+    void smooth (int amrlev, int mglev, MF& sol, MultiFab const& rhs,
+                 int redblack) const;
+
+    void compresid (int amrlev, int mglev, MF& resid, MF const& b) const;
+
+    void applyPhysBC (int amrlev, int mglev, MultiFab& mf) const;
+
+private:
+
+    void applyBC (int amrlev, int mglev, MF& in) const;
+    void applyBC (int amrlev, int mglev, MultiFab& mf) const;
+
+    [[nodiscard]] iMultiFab const& getDotMask (int amrlev, int mglev, int idim) const;
+
+    [[nodiscard]] int getDirichlet (int amrlev, int mglev, int idim, int face) const;
+
+    RT m_alpha = std::numeric_limits<RT>::lowest();
+    RT m_beta  = std::numeric_limits<RT>::lowest();
+
+    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;
+};
+
+}
+
+#endif