usable mhd simulator

updated + tested mhd_state, updated mhd_model

revert back my typo on hybrid_tagger

fixed typo in test_mhd_state

preparing mhd usablevecfield

new field initialisation logic

Added end of file newline in  field_user_initializer.hpp

added mhd state fixture

wrote test for mhd_state_fixtures

added mhd yee grid, modified gridlayout to handle both hybrid and mhd

up

init_functions

MHD solver

solver mhd update

more solver updates

preparing rebase

formatting

riemann reconstruction

quick commit fix

solved conflict

more solver updates

initial branch commit: new architecture for reconstruction + flux computation

improved test, fixed some bugs in godunov-fluxes

better timing directory creation (#914)

safer nu (#911)

try fallback for dl on keyerror + ruffage (#916)

convenient utils for mpi/hierarchies

keep pyattrs on compute_from_hier

nan min/max to handle possible nan ghosts (#923)

* nan min/max to handle possible nan ghosts

* flatten

rm atefact file

fixed missing template keyword for macos-12 build

more explicit unwrapping in godunov fluxes with variadic arguments of unclear size

fixed lambda capture problem

fixed lambda capture problem

added time integration

added projection functions for centering, primitive/conservetive converter

added projection functions for centering, primitive/conservetive converter

solver ref binding in a lambda issue

added template deduction guide on for macOS compile on a struct of MHDSolver

New approach for template deduction for the struct in MHDSolver

added usable ghost cells

syntaxe fixes for Werror

full boundary conditions in solver

changed ghost cells function names to be more consistant with hybrid

corrected dummyhierarchy in testmhdsolver (operation on nullptr)

MHDMock simulator setup

added pybind wrappers and first tests

successful orszag-tang test

fixed typo

fixed typo

fixed typo

bug fixes (wrong usage of prevIndex nextIndex)

Author: UCaromel

.gitignore

@@ -22,3 +22,4 @@ perf.*
 .vscode
 .phare*
 PHARE_REPORT.zip
+

pyphare/pyphare/mock_mhd_simulator/__init__.py

@@ -0,0 +1,94 @@
+import numpy as np
+
+def is_scalar(arg):
+    return not isinstance(arg, (list, tuple)) and not is_nd_array(arg)
+
+def is_nd_array(arg):
+    return isinstance(arg, np.ndarray)
+
+
+# converts scalars to array of expected size
+# converts lists to arrays
+class py_fn_wrapper:
+    def __init__(self, fn):
+        self.fn = fn
+
+    def __call__(self, *xyz):
+        args = [np.asarray(arg) for arg in xyz]
+        ret = self.fn(*args)
+        if isinstance(ret, list):
+            ret = np.asarray(ret)
+        if is_scalar(ret):
+            ret = np.full(len(args[-1]), ret)
+        return ret
+
+
+# Wrap calls to user init functions to turn C++ vectors to ndarrays,
+#  and returned ndarrays to C++ span
+class fn_wrapper(py_fn_wrapper):
+    def __init__(self, fn):
+        super().__init__(fn)
+
+    def __call__(self, *xyz):
+        from pyphare.cpp import cpp_etc_lib
+
+        # convert numpy array to C++ SubSpan
+        # couples vector init functions to C++
+        return cpp_etc_lib().makePyArrayWrapper(super().__call__(*xyz))
+
+
+def clearDict():
+    """
+    dict may contain dangling references from a previous simulation unless cleared
+    """
+    import pybindlibs.dictator as pp
+
+    pp.stop()
+
+def populateDict():
+    from .global_vars import sim as simulation
+    import pybindlibs.dictator as pp
+
+    def add_int(path, val):
+        pp.add_int(path, int(val))
+
+    def add_double(path, val):
+        pp.add_double(path, float(val))
+
+    add_string = pp.add_string
+    addInitFunction = getattr(pp, "addInitFunction{:d}".format(simulation.ndim) + "D")
+
+    add_double("time_step", simulation.timestep)
+    add_double("final_time", simulation.final_time)
+
+    add_int("nbr_cells/x", simulation.cells[0])
+    add_double("mesh_size/x", simulation.dl[0])
+    add_double("origin/x", simulation.origin[0])
+
+    if simulation.ndim > 1:
+        add_int("nbr_cells/y", simulation.cells[1])
+        add_double("mesh_size/y", simulation.dl[1])
+        add_double("origin/y", simulation.origin[1])
+
+        if simulation.ndim > 2:
+            add_int("nbr_cells/z", simulation.cells[2])
+            add_double("mesh_size/z", simulation.dl[2])
+            add_double("origin/z", simulation.origin[2])
+
+    add_double("godunov/resistivity", simulation.eta)
+    add_double("godunov/hyper_resistivity", simulation.nu)
+    add_double("godunov/heat_capacity_ratio", simulation.gamma)
+    add_string("godunov/terms", simulation.terms)
+    add_double("to_primitive/heat_capacity_ratio", simulation.gamma)
+    add_double("to_conservative/heat_capacity_ratio", simulation.gamma)
+
+    d=simulation.model.model_dict
+
+    addInitFunction("density/initializer", fn_wrapper(d["density"]))
+    addInitFunction("velocity/initializer/x_component", fn_wrapper(d["vx"]))
+    addInitFunction("velocity/initializer/y_component", fn_wrapper(d["vy"]))
+    addInitFunction("velocity/initializer/z_component", fn_wrapper(d["vz"]))
+    addInitFunction("magnetic/initializer/x_component", fn_wrapper(d["bx"]))
+    addInitFunction("magnetic/initializer/y_component", fn_wrapper(d["by"]))
+    addInitFunction("magnetic/initializer/z_component", fn_wrapper(d["bz"]))
+    addInitFunction("pressure/initializer", fn_wrapper(d["p"]))

pyphare/pyphare/mock_mhd_simulator/global_vars.py

@@ -0,0 +1 @@
+sim = None

pyphare/pyphare/mock_mhd_simulator/mhd_model.py

@@ -0,0 +1,44 @@
+from . import global_vars
+
+class MHDModel(object):
+    def defaulter(self, input, value):
+        if input is not None:
+            import inspect
+
+            params = list(inspect.signature(input).parameters.values())
+            assert len(params)
+            param_per_dim = len(params) == self.dim
+            has_vargs = params[0].kind == inspect.Parameter.VAR_POSITIONAL
+            assert param_per_dim or has_vargs
+            return input
+        if self.dim == 1:
+            return lambda x: value + x * 0
+        if self.dim == 2:
+            return lambda x, y: value
+        if self.dim == 3:
+            return lambda x, y, z: value
+
+    def __init__(self, density=None, vx=None, vy=None, vz=None, bx=None, by=None, bz=None, p=None):
+        if global_vars.sim is None:
+            raise RuntimeError("A simulation must be declared before a model")
+
+        if global_vars.sim.model is not None:
+            raise RuntimeError("A model is already created")
+
+        self.dim = global_vars.sim.ndim
+
+        density = self.defaulter(density, 1.0)
+        vx = self.defaulter(vx, 1.0)
+        vy = self.defaulter(vy, 0.0)
+        vz = self.defaulter(vz, 0.0)
+        bx = self.defaulter(bx, 1.0)
+        by = self.defaulter(by, 0.0)
+        bz = self.defaulter(bz, 0.0)
+        p = self.defaulter(p, 1.0)
+
+        self.model_dict = {}
+
+        self.model_dict.update({"density": density, "vx": vx, "vy": vy, "vz": vz, "bx": bx, "by": by, "bz": bz, "p": p})
+
+        global_vars.sim.set_model(self)
+

pyphare/pyphare/mock_mhd_simulator/simulation.py

@@ -0,0 +1,16 @@
+from . import global_vars
+
+class Simulation(object):
+    def __init__(self, **kwargs):
+        if global_vars.sim is not None:
+            raise RuntimeError("simulation is already created")
+
+        global_vars.sim = self
+
+        for k, v in kwargs.items():
+            object.__setattr__(self, k, v)
+
+        self.model = None
+
+    def set_model(self, model):
+        self.model = model

pyphare/pyphare/mock_mhd_simulator/simulator.py

@@ -0,0 +1,48 @@
+from . import clearDict
+from . import populateDict
+
+def pyMHD():
+    import importlib
+
+    return importlib.import_module("pybindlibs.pyMHD")
+
+def make_cpp_simulator(dim, interp):
+    import pybindlibs.pyMHD
+
+    make_sim = f"make_mhd_mock_simulator_{dim}_{interp}"
+    return getattr(pyMHD(), make_sim)()
+
+class MHDMockSimulator:
+    def __init__(self, simulation):
+        self.cpp_sim = None
+        self.simulation = simulation
+
+    def __del__(self):
+        self.reset()
+
+    def reset(self):
+        if self.cpp_sim is not None:
+            clearDict()
+        self.cpp_sim = None
+        return self
+
+    def initialize(self):
+        if self.cpp_sim is not None:
+            raise ValueError(
+                "Simulator already initialized: requires reset to re-initialize"
+            )
+
+        populateDict()
+        self.cpp_sim = make_cpp_simulator(
+            self.simulation.ndim,
+            self.simulation.order,
+        )
+
+    def _check_init(self):
+        if self.cpp_sim is None:
+            self.initialize()
+
+    def run(self, filename, dumpfrequency=1):
+        self._check_init()
+        self.cpp_sim.advance(filename, dumpfrequency)
+        return self

res/cmake/test.cmake

@@ -17,6 +17,7 @@ if (test AND ${PHARE_EXEC_LEVEL_MIN} GREATER 0) # 0 = no tests
   add_subdirectory(tests/core/data/ion_population)
   add_subdirectory(tests/core/data/maxwellian_particle_initializer)
   add_subdirectory(tests/core/data/particle_initializer)
+  add_subdirectory(tests/core/data/mhd_state)
   add_subdirectory(tests/core/utilities/box)
   add_subdirectory(tests/core/utilities/range)
   add_subdirectory(tests/core/utilities/index)

src/amr/messengers/mhd_messenger_info.hpp

@@ -2,6 +2,7 @@
 #ifndef PHARE_MHD_MESSENGER_INFO_HPP
 #define PHARE_MHD_MESSENGER_INFO_HPP
 
+#include "core/data/vecfield/vecfield.hpp"
 #include "messenger_info.hpp"
 
 
@@ -12,7 +13,9 @@ namespace amr
 {
     class MHDMessengerInfo : public IMessengerInfo
     {
+        using VecFieldNames = core::VecFieldNames;
     public:
+        // What is needed here ?
         virtual ~MHDMessengerInfo() = default;
     };
 

src/amr/messengers/mhd_messenger_stategy.hpp

@@ -16,60 +16,87 @@
 
 
 
-namespace PHARE
+namespace PHARE::solver
 {
-namespace solver
+template<typename GridLayoutT, typename VecFieldT, typename AMR_Types, typename Grid_t>
+class MHDModel : public IPhysicalModel<AMR_Types>
 {
-    template<typename GridLayoutT, typename VecFieldT, typename AMR_Types, typename Grid_t>
-    class MHDModel : public IPhysicalModel<AMR_Types>
-    {
-    public:
-        using patch_t   = typename AMR_Types::patch_t;
-        using level_t   = typename AMR_Types::level_t;
-        using Interface = IPhysicalModel<AMR_Types>;
+public:
+    using patch_t   = typename AMR_Types::patch_t;
+    using level_t   = typename AMR_Types::level_t;
+    using Interface = IPhysicalModel<AMR_Types>;
+
+    using field_type      = typename VecFieldT::field_type;
+    using gridlayout_type = GridLayoutT;
+
+    static const inline std::string model_name = "MHDModel";
+    static constexpr auto dimension            = gridlayout_type::dimension;
+    using resources_manager_type               = amr::ResourcesManager<gridlayout_type, Grid_t>;
+
+    core::MHDState<VecFieldT> state;
+    std::shared_ptr<resources_manager_type> resourcesManager;
+
+    virtual void initialize(level_t& level) override;
+
 
-        static const inline std::string model_name = "MHDModel";
-        static constexpr auto dimension            = GridLayoutT::dimension;
-        using resources_manager_type               = amr::ResourcesManager<GridLayoutT, Grid_t>;
+    virtual void allocate(patch_t& patch, double const allocateTime) override
+    {
+        resourcesManager->allocate(state, patch, allocateTime);
+    }
 
 
-        explicit MHDModel(std::shared_ptr<resources_manager_type> const& _resourcesManager)
-            : IPhysicalModel<AMR_Types>{model_name}
-            , resourcesManager{std::move(_resourcesManager)}
-        {
-        }
 
-        virtual void initialize(level_t& /*level*/) override {}
+    virtual void
+    fillMessengerInfo(std::unique_ptr<amr::IMessengerInfo> const& /*info*/) const override
+    {
+    }
 
+    NO_DISCARD auto setOnPatch(patch_t& patch)
+    {
+        return resourcesManager->setOnPatch(patch, *this);
+    }
+
+    explicit MHDModel(PHARE::initializer::PHAREDict const& dict,
+                      std::shared_ptr<resources_manager_type> const& _resourcesManager)
+        : IPhysicalModel<AMR_Types>{model_name}
+        , state{dict}
+        , resourcesManager{std::move(_resourcesManager)}
+    {
+    }
 
-        virtual void allocate(patch_t& patch, double const allocateTime) override
-        {
-            resourcesManager->allocate(state.B, patch, allocateTime);
-            resourcesManager->allocate(state.V, patch, allocateTime);
-        }
+    virtual ~MHDModel() override = default;
 
 
+    //-------------------------------------------------------------------------
+    //                  start the ResourcesUser interface
+    //-------------------------------------------------------------------------
 
-        virtual void
-        fillMessengerInfo(std::unique_ptr<amr::IMessengerInfo> const& /*info*/) const override
-        {
-        }
+    NO_DISCARD bool isUsable() const { return state.isUsable(); }
 
-        NO_DISCARD auto setOnPatch(patch_t& patch)
-        {
-            return resourcesManager->setOnPatch(patch, *this);
-        }
+    NO_DISCARD bool isSettable() const { return state.isSettable(); }
 
+    NO_DISCARD auto getCompileTimeResourcesViewList() const { return std::forward_as_tuple(state); }
 
-        virtual ~MHDModel() override = default;
+    NO_DISCARD auto getCompileTimeResourcesViewList() { return std::forward_as_tuple(state); }
 
-        core::MHDState<VecFieldT> state;
-        std::shared_ptr<resources_manager_type> resourcesManager;
-    };
+    //-------------------------------------------------------------------------
+    //                  ends the ResourcesUser interface
+    //-------------------------------------------------------------------------
+};
 
+template<typename GridLayoutT, typename VecFieldT, typename AMR_Types, typename Grid_t>
+void MHDModel<GridLayoutT, VecFieldT, AMR_Types, Grid_t>::initialize(level_t& level)
+{
+    for (auto& patch : level)
+    {
+        auto layout = amr::layoutFromPatch<GridLayoutT>(*patch);
+        auto _      = this->resourcesManager->setOnPatch(*patch, state);
 
+        state.initialize(layout);
+    }
+    resourcesManager->registerForRestarts(*this);
+}
 
-} // namespace solver
-} // namespace PHARE
+} // namespace PHARE::solver
 
 #endif