divergence free refluxing

BROKEN

adding the patch refine strategy for the magnetic field

3D for magnetic regriding

the path of the SAMRAI

cleanup and comments

back to NaNs

3d fixes + refactoring

back to iterating on the fine field box (in mpi the fine_box might be different from the patch box)

pr comments

initial commit, E accumulation

BROKEN

refluxing without divB

new magnetic field refinement implementing Toth and Roe (2002)

adding the patch refine strategy for the magnetic field

forcing even number of ghost cells

3D for magnetic regriding

the path of the SAMRAI

cleanup and comments

back to NaNs

3d fixes + refactoring

back to iterating on the fine field box (in mpi the fine_box might be different from the patch box)

pr comments

initial commit, E accumulation

BROKEN

adding the patch refine strategy for the magnetic field

3D for magnetic regriding

the path of the SAMRAI

cleanup and comments

back to NaNs

3d fixes + refactoring

back to iterating on the fine field box (in mpi the fine_box might be different from the patch box)

pr comments

BROKEN

refluxing without divB

rebase master

Author: UCaromel

src/amr/CMakeLists.txt

@@ -12,6 +12,7 @@ set( SOURCES_INC
      data/field/coarsening/coarsen_weighter.hpp
      data/field/coarsening/default_field_coarsener.hpp
      data/field/coarsening/magnetic_field_coarsener.hpp
+     data/field/coarsening/flux_sum_coarsener.hpp
      data/field/field_data.hpp
      data/field/field_data_factory.hpp
      data/field/field_geometry.hpp

src/amr/data/field/coarsening/flux_sum_coarsener.hpp

@@ -0,0 +1,158 @@
+#ifndef PHARE_FLUX_SUM_COARSENER
+#define PHARE_FLUX_SUM_COARSENER
+
+#include "core/data/grid/gridlayoutdefs.hpp"
+#include "core/utilities/constants.hpp"
+#include "amr/resources_manager/amr_utils.hpp"
+
+
+#include <SAMRAI/hier/Box.h>
+#include <stdexcept>
+
+namespace PHARE::amr
+{
+using core::dirX;
+using core::dirY;
+using core::dirZ;
+/** @brief This class gives an operator() that performs the coarsening of N fine nodes onto a
+ * given coarse node
+ *
+ * A MagneticFieldCoarsener object is created each time the refine() method of the
+ * FieldCoarsenOperator is called and its operator() is called for each coarse index.
+ * It is the default coarsening policy and used for any field that does not come with
+ * specific constraints (such as conserving some property in the coarsening process).
+ *
+ *
+ * This coarsening operation is defined so to conserve the magnetic flux.
+ * This is done by assigning to a magnetic field component on a coarse face, the average
+ * of the enclosed fine faces
+ *
+ */
+template<std::size_t dimension>
+class FluxSumCoarsener
+{
+public:
+    FluxSumCoarsener(std::array<core::QtyCentering, dimension> const centering,
+                     SAMRAI::hier::Box const& sourceBox, SAMRAI::hier::Box const& destinationBox,
+                     SAMRAI::hier::IntVector const& ratio)
+        : centering_{centering}
+        , sourceBox_{sourceBox}
+        , destinationBox_{destinationBox}
+
+    {
+    }
+
+    template<typename FieldT>
+    void operator()(FieldT const& fineField, FieldT& coarseField,
+                    core::Point<int, dimension> coarseIndex)
+    {
+        // For the moment we only take the case of field with the same centering
+        TBOX_ASSERT(fineField.physicalQuantity() == coarseField.physicalQuantity());
+
+        core::Point<int, dimension> fineStartIndex;
+
+        fineStartIndex[dirX] = coarseIndex[dirX] * this->ratio_;
+
+        if constexpr (dimension > 1)
+        {
+            fineStartIndex[dirY] = coarseIndex[dirY] * this->ratio_;
+            if constexpr (dimension > 2)
+            {
+                fineStartIndex[dirZ] = coarseIndex[dirZ] * this->ratio_;
+            }
+        }
+
+        fineStartIndex = AMRToLocal(fineStartIndex, sourceBox_);
+        coarseIndex    = AMRToLocal(coarseIndex, destinationBox_);
+
+        if constexpr (dimension == 1)
+        {
+            if (centering_[dirX] == core::QtyCentering::dual) // ex
+            {
+                coarseField(coarseIndex[dirX])
+                    = 0.5 * (fineField(fineStartIndex[dirX] + 1) + fineField(fineStartIndex[dirX]));
+            }
+            else if (centering_[dirX] == core::QtyCentering::primal) // ey, ez
+            {
+                coarseField(coarseIndex[dirX]) = fineField(fineStartIndex[dirX]);
+            }
+        }
+
+        if constexpr (dimension == 2)
+        {
+            if (centering_[dirX] == core::QtyCentering::dual
+                and centering_[dirY] == core::QtyCentering::primal) // ex
+            {
+                coarseField(coarseIndex[dirX], coarseIndex[dirY])
+                    = 0.5
+                      * (fineField(fineStartIndex[dirX], fineStartIndex[dirY])
+                         + fineField(fineStartIndex[dirX] + 1, fineStartIndex[dirY]));
+            }
+            else if (centering_[dirX] == core::QtyCentering::primal
+                     and centering_[dirY] == core::QtyCentering::dual) // ey
+            {
+                coarseField(coarseIndex[dirX], coarseIndex[dirY])
+                    = 0.5
+                      * (fineField(fineStartIndex[dirX], fineStartIndex[dirY])
+                         + fineField(fineStartIndex[dirX], fineStartIndex[dirY] + 1));
+            }
+            else if (centering_[dirX] == core::QtyCentering::primal
+                     and centering_[dirY] == core::QtyCentering::primal) // ez
+            {
+                coarseField(coarseIndex[dirX], coarseIndex[dirY])
+                    = fineField(fineStartIndex[dirX], fineStartIndex[dirY]);
+            }
+            else
+            {
+                throw std::runtime_error("no electric field should end up here");
+            }
+        }
+        else if constexpr (dimension == 3)
+        {
+            if (centering_[dirX] == core::QtyCentering::dual
+                and centering_[dirY] == core::QtyCentering::primal
+                and centering_[dirZ] == core::QtyCentering::primal) // ex
+            {
+                coarseField(coarseIndex[dirX], coarseIndex[dirY], coarseIndex[dirZ])
+                    = 0.5
+                      * (fineField(fineStartIndex[dirX], fineStartIndex[dirY], fineStartIndex[dirZ])
+                         + fineField(fineStartIndex[dirX] + 1, fineStartIndex[dirY],
+                                     fineStartIndex[dirZ]));
+            }
+            else if (centering_[dirX] == core::QtyCentering::primal
+                     and centering_[dirY] == core::QtyCentering::dual
+                     and centering_[dirZ] == core::QtyCentering::primal) // ey
+            {
+                coarseField(coarseIndex[dirX], coarseIndex[dirY], coarseIndex[dirZ])
+                    = 0.5
+                      * (fineField(fineStartIndex[dirX], fineStartIndex[dirY], fineStartIndex[dirZ])
+                         + fineField(fineStartIndex[dirX], fineStartIndex[dirY] + 1,
+                                     fineStartIndex[dirZ]));
+            }
+            else if (centering_[dirX] == core::QtyCentering::primal
+                     and centering_[dirY] == core::QtyCentering::primal
+                     and centering_[dirZ] == core::QtyCentering::dual) // ez
+            {
+                coarseField(coarseIndex[dirX], coarseIndex[dirY], coarseIndex[dirZ])
+                    = 0.5
+                      * (fineField(fineStartIndex[dirX], fineStartIndex[dirY], fineStartIndex[dirZ])
+                         + fineField(fineStartIndex[dirX], fineStartIndex[dirY],
+                                     fineStartIndex[dirZ] + 1));
+            }
+            else
+            {
+                throw std::runtime_error("no electric field should end up here");
+            }
+        }
+    }
+
+private:
+    std::array<core::QtyCentering, dimension> const centering_;
+    SAMRAI::hier::Box const sourceBox_;
+    SAMRAI::hier::Box const destinationBox_;
+    static int constexpr ratio_ = 2;
+};
+
+} // namespace PHARE::amr
+
+#endif

src/amr/data/field/coarsening/magnetic_field_coarsener.hpp

@@ -70,7 +70,7 @@ class MagneticFieldCoarsener
         coarseIndex    = AMRToLocal(coarseIndex, destinationBox_);
 
         // the following kinda assumes where B is, i.e. Yee layout centering
-        // as it only does faces pirmal-dual, dual-primal and dual-dual
+        // as it only does faces primal-dual, dual-primal and dual-dual
 
         if constexpr (dimension == 1)
         {

src/amr/data/field/refine/magnetic_refine_patch_strategy.hpp

@@ -67,8 +67,8 @@ class MagneticRefinePatchStrategy : public SAMRAI::xfer::RefinePatchStrategy
     {
     }
 
-    // We compute the values of the new fine magnetic faces using what was already refined, ie the
-    // values on the old coarse faces.
+    // We compute the values of the new fine magnetic faces using what was already refined, ie
+    // the values on the old coarse faces.
     void postprocessRefine(SAMRAI::hier::Patch& fine, SAMRAI::hier::Patch const& coarse,
                            SAMRAI::hier::Box const& fine_box,
                            SAMRAI::hier::IntVector const& ratio) override

src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp

@@ -10,14 +10,15 @@
 #include "core/utilities/index/index.hpp"
 #include "refiner_pool.hpp"
 #include "synchronizer_pool.hpp"
+#include "amr/data/field/coarsening/field_coarsen_operator.hpp"
 #include "amr/data/field/coarsening/default_field_coarsener.hpp"
 #include "amr/data/field/coarsening/magnetic_field_coarsener.hpp"
+#include "amr/data/field/coarsening/flux_sum_coarsener.hpp"
 #include "amr/data/field/refine/field_refiner.hpp"
 #include "amr/data/field/refine/magnetic_field_refiner.hpp"
 #include "amr/data/field/refine/electric_field_refiner.hpp"
 #include "amr/data/field/time_interpolate/field_linear_time_interpolate.hpp"
 #include "amr/data/field/refine/field_refine_operator.hpp"
-#include "amr/data/field/coarsening/field_coarsen_operator.hpp"
 #include "amr/messengers/messenger_info.hpp"
 #include "amr/messengers/hybrid_messenger_info.hpp"
 #include "amr/messengers/hybrid_messenger_strategy.hpp"
@@ -35,6 +36,8 @@
 
 #include "SAMRAI/xfer/RefineAlgorithm.h"
 #include "SAMRAI/xfer/RefineSchedule.h"
+#include "SAMRAI/xfer/CoarsenAlgorithm.h"
+#include "SAMRAI/xfer/CoarsenSchedule.h"
 #include "SAMRAI/xfer/BoxGeometryVariableFillPattern.h"
 
 
@@ -99,6 +102,7 @@ namespace amr
         using BaseCoarsenOp     = FieldCoarsenOperator<GridLayoutT, GridT, Policy>;
         using MagneticCoarsenOp = BaseCoarsenOp<MagneticFieldCoarsener<dimension>>;
         using DefaultCoarsenOp  = BaseCoarsenOp<DefaultFieldCoarsener<dimension>>;
+        using FluxSumCoarsenOp  = BaseCoarsenOp<FluxSumCoarsener<dimension>>;
 
     public:
         static inline std::string const stratName    = "HybridModel-HybridModel";
@@ -185,6 +189,42 @@ namespace amr
             BalgoNode.registerRefine(*bz_id, *bz_id, *bz_id, BfieldNodeRefineOp_,
                                      zVariableFillPattern);
 
+            // refluxing
+            // we first want to coarsen the flux sum onto the coarser level
+            auto ex_reflux_id = resourcesManager_->getID(hybridInfo->refluxElectric.xName);
+            auto ey_reflux_id = resourcesManager_->getID(hybridInfo->refluxElectric.yName);
+            auto ez_reflux_id = resourcesManager_->getID(hybridInfo->refluxElectric.zName);
+
+            auto ex_fluxsum_id = resourcesManager_->getID(hybridInfo->fluxSumElectric.xName);
+            auto ey_fluxsum_id = resourcesManager_->getID(hybridInfo->fluxSumElectric.yName);
+            auto ez_fluxsum_id = resourcesManager_->getID(hybridInfo->fluxSumElectric.zName);
+
+            if (!ex_reflux_id or !ey_reflux_id or !ez_reflux_id or !ex_fluxsum_id or !ey_fluxsum_id
+                or !ez_fluxsum_id)
+            {
+                throw std::runtime_error(
+                    "HybridHybridMessengerStrategy: missing electric field variable IDs");
+            }
+
+            RefluxAlgo.registerCoarsen(*ex_reflux_id, *ex_fluxsum_id, fluxSumCoarseningOp_,
+                                       xVariableFillPattern);
+            RefluxAlgo.registerCoarsen(*ey_reflux_id, *ey_fluxsum_id, fluxSumCoarseningOp_,
+                                       yVariableFillPattern);
+            RefluxAlgo.registerCoarsen(*ez_reflux_id, *ez_fluxsum_id, fluxSumCoarseningOp_,
+                                       zVariableFillPattern);
+
+            // we then need to refill the ghosts so that they agree with the newly refluxed cells
+
+            GhostRefluxedAlgo.registerRefine(*ex_reflux_id, *ex_reflux_id, *ex_reflux_id,
+                                             EfieldRefineOp_, xVariableFillPattern);
+
+            GhostRefluxedAlgo.registerRefine(*ey_reflux_id, *ey_reflux_id, *ey_reflux_id,
+                                             EfieldRefineOp_, yVariableFillPattern);
+
+            GhostRefluxedAlgo.registerRefine(*ez_reflux_id, *ez_reflux_id, *ez_reflux_id,
+                                             EfieldRefineOp_, zVariableFillPattern);
+
+
             registerGhostComms_(hybridInfo);
             registerInitComms(hybridInfo);
             registerSyncComms(hybridInfo);
@@ -210,6 +250,9 @@ namespace amr
             magGhostsRefineSchedules[levelNumber] = Balgo.createSchedule(
                 level, levelNumber - 1, hierarchy, &magneticRefinePatchStrategy_);
 
+            // technically not needed for finest
+            ghostRefluxedSchedules[levelNumber] = GhostRefluxedAlgo.createSchedule(level);
+
             elecSharedNodesRefiners_.registerLevel(hierarchy, level);
             currentSharedNodesRefiners_.registerLevel(hierarchy, level);
 
@@ -227,9 +270,14 @@ namespace amr
             // TODO this 'if' may not be OK if L0 is regrided
             if (levelNumber != rootLevelNumber)
             {
+                // refluxing
+                auto const& coarseLevel      = hierarchy->getPatchLevel(levelNumber - 1);
+                refluxSchedules[levelNumber] = RefluxAlgo.createSchedule(coarseLevel, level);
+
                 // those are for refinement
                 magInitRefineSchedules[levelNumber] = Balgo.createSchedule(
                     level, nullptr, levelNumber - 1, hierarchy, &magneticRefinePatchStrategy_);
+
                 electricInitRefiners_.registerLevel(hierarchy, level);
                 domainParticlesRefiners_.registerLevel(hierarchy, level);
                 lvlGhostPartOldRefiners_.registerLevel(hierarchy, level);
@@ -544,6 +592,7 @@ namespace amr
 
 
 
+
         /**
          * @brief prepareStep is the concrete implementation of the
          * HybridMessengerStrategy::prepareStep method For hybrid-Hybrid communications.
@@ -564,8 +613,7 @@ namespace amr
             for (auto& patch : level)
             {
                 auto dataOnPatch = resourcesManager_->setOnPatch(
-                    *patch, hybridModel.state.electromag, hybridModel.state.J,
-                    hybridModel.state.ions, Jold_, NiOld_, ViOld_);
+                    *patch, hybridModel.state.J, hybridModel.state.ions, Jold_, NiOld_, ViOld_);
 
                 resourcesManager_->setTime(Jold_, *patch, currentTime);
                 resourcesManager_->setTime(NiOld_, *patch, currentTime);
@@ -623,12 +671,20 @@ namespace amr
             PHARE_LOG_LINE_STR("synchronizing level " + std::to_string(levelNumber));
 
             // call coarsning schedules...
-            magnetoSynchronizers_.sync(levelNumber);
+            // magnetoSynchronizers_.sync(levelNumber);
             electroSynchronizers_.sync(levelNumber);
             densitySynchronizers_.sync(levelNumber);
             ionBulkVelSynchronizers_.sync(levelNumber);
         }
 
+
+        void reflux(int const coarserLevelNumber, int const fineLevelNumber,
+                    double const syncTime) override
+        {
+            refluxSchedules[fineLevelNumber]->coarsenData();
+            ghostRefluxedSchedules[coarserLevelNumber]->fillData(syncTime);
+        }
+
         // after coarsening, domain nodes have been updated and therefore patch ghost nodes
         // will probably stop having the exact same value as their overlapped neighbor
         // domain node we thus fill ghost nodes. note that we first fill shared border nodes
@@ -654,7 +710,7 @@ namespace amr
             // level border with next coarser model B would invalidate divB on the first
             // fine domain cell since its border face only received a fraction of the
             // induction that has occured on the shared coarse face.
-            magPatchGhostsRefineSchedules[levelNumber]->fillData(time);
+            // magPatchGhostsRefineSchedules[levelNumber]->fillData(time);
             elecGhostsRefiners_.fill(hybridModel.state.electromag.E, levelNumber, time);
             rhoGhostsRefiners_.fill(levelNumber, time);
             velGhostsRefiners_.fill(hybridModel.state.ions.velocity(), levelNumber, time);
@@ -973,7 +1029,6 @@ namespace amr
 
 
 
-
         VecFieldT Jold_{stratName + "_Jold", core::HybridQuantity::Vector::J};
         VecFieldT ViOld_{stratName + "_VBulkOld", core::HybridQuantity::Vector::V};
         FieldT NiOld_{stratName + "_NiOld", core::HybridQuantity::Scalar::rho};
@@ -1011,6 +1066,10 @@ namespace amr
         std::map<int, std::shared_ptr<SAMRAI::xfer::RefineSchedule>> magPatchGhostsRefineSchedules;
         std::map<int, std::shared_ptr<SAMRAI::xfer::RefineSchedule>> magSharedNodeRefineSchedules;
 
+        SAMRAI::xfer::CoarsenAlgorithm RefluxAlgo{SAMRAI::tbox::Dimension{dimension}};
+        SAMRAI::xfer::RefineAlgorithm GhostRefluxedAlgo;
+        std::map<int, std::shared_ptr<SAMRAI::xfer::CoarsenSchedule>> refluxSchedules;
+        std::map<int, std::shared_ptr<SAMRAI::xfer::RefineSchedule>> ghostRefluxedSchedules;
 
         //! store refiners for electric fields that need ghosts to be filled
         SharedNodeRefinerPool elecSharedNodesRefiners_{resourcesManager_};
@@ -1069,6 +1128,7 @@ namespace amr
         using CoarsenOperator_ptr = std::shared_ptr<SAMRAI::hier::CoarsenOperator>;
         CoarsenOperator_ptr fieldCoarseningOp_{std::make_shared<DefaultCoarsenOp>()};
         CoarsenOperator_ptr magneticCoarseningOp_{std::make_shared<MagneticCoarsenOp>()};
+        CoarsenOperator_ptr fluxSumCoarseningOp_{std::make_shared<FluxSumCoarsenOp>()};
 
         MagneticRefinePatchStrategy<ResourcesManagerT, FieldDataT> magneticRefinePatchStrategy_{
             *resourcesManager_};

src/amr/messengers/hybrid_messenger.hpp

@@ -190,6 +190,14 @@ namespace amr
 
         void synchronize(SAMRAI::hier::PatchLevel& level) override { strat_->synchronize(level); }
 
+
+        void reflux(int const coarserLevelNumber, int const fineLevelNumber,
+                    double const syncTime) override
+        {
+            strat_->reflux(coarserLevelNumber, fineLevelNumber, syncTime);
+        }
+
+
         void postSynchronize(IPhysicalModel& model, SAMRAI::hier::PatchLevel& level,
                              double const time) override
         {
@@ -345,7 +353,7 @@ namespace amr
         virtual ~HybridMessenger() = default;
 
     private:
-        const std::unique_ptr<stratT> strat_;
+        std::unique_ptr<stratT> const strat_;
     };
 
 

src/amr/messengers/hybrid_messenger_info.hpp

@@ -67,6 +67,10 @@ namespace amr
         std::vector<VecFieldNames> ghostCurrent;
         std::vector<VecFieldNames> ghostBulkVelocity;
 
+        // below are the descriptions of the electric field that we use in the refluxing
+        VecFieldNames refluxElectric;
+        VecFieldNames fluxSumElectric;
+
         virtual ~HybridMessengerInfo() = default;
     };