EdgeFluxRegister for MHD

This implements EdgeFluxRegister for refluxing in MHD. The equation here is

  \frac{\partial B}{\partial t} + curl E = 0.

Here, B is on faces and E is on edges.

Src/AmrCore/AMReX_AmrCore.cpp

@@ -1,6 +1,5 @@
 
 #include <AMReX_AmrCore.H>
-#include <AMReX_Print.H>
 
 #ifdef AMREX_PARTICLES
 #include <AMReX_AmrParGDB.H>

Src/Base/AMReX_BaseFab.H

@@ -543,7 +543,7 @@ public:
  int numcomp = 1) noexcept;
  /**
  * \brief As above, except that the destination Box is specified,
- * but the source Box is taken to the equal to the source
+ * but the source Box is taken to the equal to the destination
  * Box, and all components of the destination BaseFab are
  * copied.
  */

Src/Boundary/AMReX_BoundaryFwd.H

@@ -11,6 +11,7 @@ class InterpBndryDataT;
 class Mask;
 class MultiMask;
 class YAFluxRegisterT;
+class EdgeFluxRegister;
 
 }
 

Src/Boundary/AMReX_EdgeFluxRegister.H

@@ -0,0 +1,98 @@
+#ifndef AMREX_EDGE_FLUX_REGISTER_H_
+#define AMREX_EDGE_FLUX_REGISTER_H_
+#include <AMReX_Config.H>
+
+#include <AMReX_iMultiFab.H>
+#include <AMReX_LayoutData.H>
+#include <AMReX_MultiFab.H>
+
+namespace amrex {
+
+/**
+ * Edfe Flux Register for Constrained Transport
+ *
+ * This Flux Register is useful for solving system like dB/dt + curl E = 0
+ * on a staggered mesh. (Here d is of course partial derivation.) B is a
+ * vector on cell faces, and E is a vector on cell edges. In 2D, E has only
+ * one component, Ez, and it is on the nodes of a 2d mesh.
+ *
+ * At the beginning of a coarse step, `reset()` is called. In MFIter for
+ * the coarse level advance, `CrseAdd` is called with coarse flux (i.e., E).
+ * The flux is not scaled. In MFIter for the fine level advance, `FineAdd`
+ * is called. After the fine level finished its time steps, `Reflux` is
+ * called to update the coarse level B on the coarse/fine boundary.
+ *
+ * Note that both CrseAdd and FineAdd are async for GPU builds. That means
+ * it's the user's responsibility to keep the FArrayBox arguments alive or
+ * call Gpu::streamSynchronize().
+ *
+ * We try to keep the interface simple by not providing overloads that
+ * specify the component index. If the user's data does not start with
+ * component 0, it can be worked around by creating alias FArrayBox and
+ * MulitiFab.
+ */
+class EdgeFluxRegister
+{
+public:
+
+ EdgeFluxRegister () = default;
+
+ EdgeFluxRegister (const BoxArray& fba, const BoxArray& cba,
+ const DistributionMapping& fdm, const DistributionMapping& cdm,
+ const Geometry& fgeom, const Geometry& cgeom,
+ int nvar = 1);
+
+ void define (const BoxArray& fba, const BoxArray& cba,
+ const DistributionMapping& fdm, const DistributionMapping& cdm,
+ const Geometry& fgeom, const Geometry& cgeom,
+ int nvar = 1);
+
+ void reset ();
+
+#if (AMREX_SPACEDIM == 3)
+
+ void CrseAdd (MFIter const& mfi, const Array<FArrayBox const*,3>& E_crse, Real dt_crse);
+ void FineAdd (MFIter const& mfi, const Array<FArrayBox const*,3>& E_fine, Real dt_fine);
+
+#else /* 2D */
+
+ void CrseAdd (MFIter const& mfi, FArrayBox const& E_crse, Real dt_crse);
+ void FineAdd (MFIter const& mfi, FArrayBox const& E_fine, Real dt_fine);
+
+#endif
+
+ void Reflux (Array<MultiFab*,AMREX_SPACEDIM> const& B_crse) const;
+
+private:
+
+ Geometry m_fine_geom;
+ Geometry m_crse_geom;
+
+ IntVect m_ratio;
+ int m_ncomp;
+
+#if (AMREX_SPACEDIM == 3)
+
+ Array<MultiFab,AMREX_SPACEDIM> m_E_crse; // on original grids
+
+ // There are AMREX_SPACEDIM*2 faces. For each face, we need to store two
+ // component. For example, at the x-faces, we need to store Ey and Ez.
+ Array<Array<MultiFab,2>,AMREX_SPACEDIM*2> m_E_fine;
+
+ // Mask on the coarse level indicating overlapping with m_E_fine
+ Array<iMultiFab,AMREX_SPACEDIM> m_fine_mask;
+
+#else
+
+ MultiFab m_E_crse;
+ Array<MultiFab,AMREX_SPACEDIM*2> m_E_fine;
+ iMultiFab m_fine_mask;
+
+#endif
+
+ LayoutData<int> m_has_cf; // Flag on the coarse level indicating c/f interface
+};
+
+}
+
+#endif