New Linear Solver: Curl of Curl (#3682)

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.

---------

Co-authored-by: Axel Huebl <[email protected]>

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

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 ----------------------------

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

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().

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
-

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";

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

Src/LinearSolvers/MLMG/AMReX_MLCurlCurl.H

@@ -0,0 +1,127 @@
+#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 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 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;
+
+    [[nodiscard]] RT normInf (int amrlev, MF const& mf, bool local) const override;
+
+    void averageDownAndSync (Vector<MF>& sol) 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