check for domain overlaps (#1129)

check for domain overlaps & test emergency dump is made

Author: PhilipDeegan

src/amr/level_initializer/hybrid_level_initializer.hpp

@@ -1,20 +1,24 @@
 #ifndef PHARE_HYBRID_LEVEL_INITIALIZER_HPP
 #define PHARE_HYBRID_LEVEL_INITIALIZER_HPP
 
-#include "amr/level_initializer/level_initializer.hpp"
-#include "amr/messengers/hybrid_messenger.hpp"
-#include "amr/messengers/messenger.hpp"
-#include "amr/physical_models/hybrid_model.hpp"
-#include "amr/physical_models/physical_model.hpp"
-#include "amr/resources_manager/amr_utils.hpp"
-#include "core/data/grid/gridlayout_utils.hpp"
-#include "core/data/ions/ions.hpp"
+#include "core/errors.hpp"
+#include "core/numerics/ohm/ohm.hpp"
+#include "core/utilities/mpi_utils.hpp"
 #include "core/numerics/ampere/ampere.hpp"
-#include "core/numerics/interpolator/interpolator.hpp"
 #include "core/numerics/moments/moments.hpp"
-#include "core/numerics/ohm/ohm.hpp"
+#include "core/data/grid/gridlayout_utils.hpp"
+#include "core/numerics/interpolator/interpolator.hpp"
+
+#include "amr/messengers/messenger.hpp"
+#include "amr/messengers/hybrid_messenger.hpp"
+#include "amr/resources_manager/amr_utils.hpp"
+#include "amr/physical_models/physical_model.hpp"
+#include "amr/level_initializer/level_initializer.hpp"
+
 #include "initializer/data_provider.hpp"
 
+#include <exception>
+
 namespace PHARE
 {
 namespace solver
@@ -36,7 +40,7 @@ namespace solver
         PHARE::core::Ohm<GridLayoutT> ohm_;
         PHARE::core::Ampere<GridLayoutT> ampere_;
 
-        inline bool isRootLevel(int levelNumber) const { return levelNumber == 0; }
+        inline bool isRootLevel(int const levelNumber) const { return levelNumber == 0; }
 
     public:
         explicit HybridLevelInitializer(PHARE::initializer::PHAREDict const& dict)
@@ -79,6 +83,17 @@ namespace solver
                 }
             }
 
+            try
+            {
+                amr::noDomainOverlapsOn(*hierarchy, levelNumber); // VERY BAD!
+            }
+            catch (std::exception const& ex)
+            {
+                PHARE_LOG_ERROR(ex.what());
+            }
+            if (core::mpi::any(core::Errors::instance().any()))
+                throw core::DictionaryException{}("ID", "HybridLevelInitializer::initialize");
+
             // now all particles are here, we must compute moments.
             auto& ions = hybridModel.state.ions;
             auto& rm   = *hybridModel.resourcesManager;
@@ -167,6 +182,7 @@ namespace solver
                     hybMessenger.fillElectricGhosts(E, level, 0.);
                 }
 
+
             // quantities have been computed on the level,like the moments and J
             // that we later in the code need to get on level ghost nodes via
             // space and TIME interpolation. We thus need to save current values

src/amr/resources_manager/amr_utils.hpp

@@ -16,6 +16,7 @@
 #include <SAMRAI/hier/BoxOverlap.h>
 #include <SAMRAI/hier/HierarchyNeighbors.h>
 #include <SAMRAI/geom/CartesianPatchGeometry.h>
+#include <stdexcept>
 
 namespace PHARE
 {
@@ -196,6 +197,33 @@ namespace amr
     }
 
 
+
+    NO_DISCARD auto inline getSameLevelNeighbors(SAMRAI::hier::Patch const& patch,
+                                                 SAMRAI::hier::PatchHierarchy const& hierarchy)
+    {
+        auto const lvlNbr = patch.getPatchLevelNumber();
+
+        return SAMRAI::hier::HierarchyNeighbors{hierarchy, lvlNbr, lvlNbr}.getSameLevelNeighbors(
+            patch.getBox(), lvlNbr);
+    }
+
+    void inline noDomainOverlapsOn(SAMRAI::hier::PatchHierarchy const& hierarchy, int const ilvl)
+    {
+        for (auto const& patch : *hierarchy.getPatchLevel(ilvl))
+            for (auto const& neighbox : getSameLevelNeighbors(*patch, hierarchy))
+                if (auto const overlap = patch->getBox() * neighbox; !overlap.empty())
+                    throw std::runtime_error(
+                        "CATASTROPHIC ERROR: Patch domain overlap detected on level: "
+                        + std::to_string(ilvl));
+    }
+
+    void inline noDomainOverlapsOn(SAMRAI::hier::PatchHierarchy const& hierarchy)
+    {
+        for (int iLevel = 0; iLevel < hierarchy.getNumberOfLevels(); ++iLevel)
+            noDomainOverlapsOn(hierarchy, iLevel);
+    }
+
+
     // potentially to replace with SAMRAI coarse to fine boundary stuff
     template<typename GridLayoutT> // fow now it gives us a box for only patch ghost layer
     NO_DISCARD auto makeNonLevelGhostBoxFor(SAMRAI::hier::Patch const& patch,
@@ -207,8 +235,7 @@ namespace amr
         auto const domBox              = phare_box_from<dimension>(domain);
         auto const particleGhostBox    = grow(domBox, GridLayoutT::nbrParticleGhosts());
 
-        SAMRAI::hier::HierarchyNeighbors const hier_nbrs{hierarchy, lvlNbr, lvlNbr};
-        auto const neighbors = hier_nbrs.getSameLevelNeighbors(domain, lvlNbr);
+        auto const neighbors = getSameLevelNeighbors(patch, hierarchy);
         std::vector<core::Box<int, GridLayoutT::dimension>> patchGhostLayerBoxes;
         patchGhostLayerBoxes.reserve(neighbors.size() + 1);
         patchGhostLayerBoxes.emplace_back(domBox);
@@ -245,7 +272,7 @@ namespace amr
                         int minLevel, int maxLevel, Args&&... args)
     {
         for (int iLevel = minLevel; iLevel < hierarchy.getNumberOfLevels() && iLevel <= maxLevel;
-             iLevel++)
+             ++iLevel)
         {
             visitLevel<GridLayout>(*hierarchy.getPatchLevel(iLevel), resman,
                                    std::forward<Action>(action), std::forward<Args...>(args...));

src/core/errors.hpp

@@ -3,6 +3,7 @@
 
 #include "core/def.hpp"
 
+#include <stdexcept>
 #include <string>
 #include <sstream>
 #include <iostream>
@@ -44,6 +45,13 @@ class DictionaryException : public std::exception
         return (msg = (*this)()).c_str();
     }
 
+    std::string id() const
+    {
+        if (!dict_.contains("ID"))
+            throw std::runtime_error("DictionaryException has no ID");
+        return (*this)["ID"].template to<std::string>();
+    }
+
 private:
     Dict_t dict_;
 };

src/simulator/simulator.hpp

@@ -262,9 +262,7 @@ void Simulator<opts>::diagnostics_init(initializer::PHAREDict const& dict)
         }
     }
 
-    this->allowEmergencyDumps = dict.contains("allow_emergency_dumps")
-                                    ? dict["allow_emergency_dumps"].template to<bool>()
-                                    : false;
+    this->allowEmergencyDumps = cppdict::get_value(dict, "allow_emergency_dumps", false);
 }
 
 
@@ -410,6 +408,12 @@ void Simulator<opts>::initialize()
 
         integrator_->initialize();
     }
+    catch (core::DictionaryException const& ex)
+    {
+        handle_dictionary_exception(ex);
+        error = std::string{"DICTIONARY EXCEPTION CAUGHT: "} + ex.what();
+        PHARE_LOG_ERROR(*error);
+    }
     catch (std::exception const& e)
     {
         error = std::string{"EXCEPTION CAUGHT: "} + e.what();
@@ -486,16 +490,14 @@ double Simulator<opts>::advance(double dt)
 template<auto opts>
 void Simulator<opts>::handle_dictionary_exception(core::DictionaryException const& ex)
 {
+    constexpr static std::array dump_exceptions{"Updater::updatePopulations",
+                                                "HybridLevelInitializer::initialize"};
+
     if (this->allowEmergencyDumps)
-    {
-        if (ex["ID"].template to<std::string>() == "Updater::updatePopulations")
-        {
-            for (auto level_nbr = 0; level_nbr < hierarchy_->getMaxNumberOfLevels(); ++level_nbr)
-            {
-                this->dMan->dump_level(level_nbr, currentTime_);
-            }
-        }
-    }
+        for (auto const& exception_id : dump_exceptions)
+            if (ex.id() == exception_id)
+                for (int ilvl = 0; ilvl < hierarchy_->getMaxNumberOfLevels(); ++ilvl)
+                    this->dMan->dump_level(ilvl, currentTime_);
 }
 
 template<auto opts>

tests/simulator/test_exceptions.py

@@ -3,9 +3,13 @@
 import unittest
 
 import pyphare.pharein as ph
+from pyphare.pharesee.run import Run
 from pyphare.simulator import simulator
 
 
+ph.NO_GUI()
+
+
 def setup_model(ppc=100):
     def density(*xyz):
         return 1.0
@@ -60,6 +64,10 @@ def vthz(*xyz):
         },
     )
     ph.ElectronModel(closure="isothermal", Te=0.12)
+
+    ph.FluidDiagnostics(  # NO TIMESTAMPS!
+        quantity="density", write_timestamps=[], population_name="protons"
+    )
     return model
 
 
@@ -72,7 +80,7 @@ def vthz(*xyz):
     "dl": 0.3,
     "diag_options": {
         "format": "phareh5",
-        "options": {"dir": out, "mode": "overwrite"},
+        "options": {"dir": out, "mode": "overwrite", "allow_emergency_dumps": True},
     },
 }
 
@@ -94,6 +102,9 @@ def test(self):
             simulator.exit_on_exception = False
             simulator.Simulator(simulation).run()
 
+        # the following will fail if there is no emergency diag
+        Run(out).GetN(0, "protons").plot(filename="exceptional_pop_density.png")
+
 
 if __name__ == "__main__":
     simulator.startMPI()