New Linear Solver: Curl of Curl

Add a new linear solver for curl (alpha curl E) + beta E = rhs in 3D, where E is an array of MultiFabs on edges, and alpha and beta are scalars. Currently this is just a prototype. It only supports periodic boundaries.
AMReX-Codes · Dec 22, 2023 · 874d525 · 874d525
1 parent 75571e2
commit 874d525
Show file tree

Hide file tree

Showing 19 changed files with 1,194 additions and 25 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/LinearSolvers/CMakeLists.txt b/Src/LinearSolvers/CMakeLists.txt
@@ -64,6 +64,15 @@ foreach(D IN LISTS AMReX_SPACEDIM)
           )
     endif ()
 
+    if (D EQUAL 3)
+       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,142 @@
+#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 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.
+ */
+class MLCurlCurl
+    : public MLLinOpT<Array<MultiFab,AMREX_SPACEDIM> >
+{
+public:
+    using MF        = Array<MultiFab,AMREX_SPACEDIM>;
+    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;
+
+    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_x (int amrlev, int mglev, MF& sol, MultiFab const& rhs,
+                   int redblack) const;
+    void smooth_y (int amrlev, int mglev, MF& sol, MultiFab const& rhs,
+                   int redblack) const;
+    void smooth_z (int amrlev, int mglev, MF& sol, MultiFab const& rhs,
+                   int redblack) const;
+
+private:
+
+    void applyBC (int amrlev, int mglev, MF& in, BCMode bc_mode, StateMode s_mode) const;
+    iMultiFab const& getDotMask (int amrlev, int mglev, int idim) const;
+
+    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;
+    static constexpr int m_ncomp = 1;
+};
+
+}
+
+#endif