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.

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

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