Usable MHD solver

.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

pyphare/pyphare/pharesee/hierarchy/hierarchy.py

@@ -275,10 +275,11 @@ def global_min(self, qty, **kwargs):
             for patch in lvl.patches:
                 pd = patch.patch_datas[qty]
                 if first:
-                    m = pd.dataset[:].min()
+                    m = np.nanmin(pd.dataset[:])
                     first = False
                 else:
-                    m = min(m, pd.dataset[:].min())
+                    data_and_min = np.concatenate(([m], pd.dataset[:].flatten()))
+                    m = np.nanmin(data_and_min)
 
         return m
 
@@ -289,10 +290,11 @@ def global_max(self, qty, **kwargs):
             for patch in lvl.patches:
                 pd = patch.patch_datas[qty]
                 if first:
-                    m = pd.dataset[:].max()
+                    m = np.nanmax(pd.dataset[:])
                     first = False
                 else:
-                    m = max(m, pd.dataset[:].max())
+                    data_and_max = np.concatenate(([m], pd.dataset[:].flatten()))
+                    m = np.nanmax(data_and_max)
 
         return m
 

pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py

@@ -40,6 +40,35 @@
 }
 
 
+def nbr_ranks(hier):
+    """
+    returns the number of mpi ranks used in the given hierarchy
+    """
+    max_rank = 0
+    t0 = hier.times()[0]
+    for _, lvl in hier.levels(t0).items():
+        for patch in lvl.patches:
+            rank = patch.attrs["mpi_rank"]
+            if rank > max_rank:
+                max_rank = rank
+    return max_rank
+
+
+def patch_per_rank(hier):
+    """
+    returns the number of patch per mpi rank for each time step
+    """
+    nbranks = nbr_ranks(hier)
+    ppr = {}
+    for t in hier.times():
+        ppr[t] = {ir: 0 for ir in np.arange(nbranks + 1)}
+        for _, lvl in hier.levels(t).items():
+            for patch in lvl.patches:
+                ppr[t][patch.attrs["mpi_rank"]] += 1
+
+    return ppr
+
+
 def are_compatible_hierarchies(hierarchies):
     ref = hierarchies[0]
     same_box = True
@@ -163,15 +192,16 @@ def new_patchdatas_from(compute, patchdatas, layout, id, **kwargs):
 def new_patches_from(compute, hierarchies, ilvl, t, **kwargs):
     reference_hier = hierarchies[0]
     new_patches = []
-    patch_nbr = len(reference_hier.level(ilvl, time=t).patches)
-    for ip in range(patch_nbr):
-        current_patch = reference_hier.level(ilvl, time=t).patches[ip]
+    ref_patches = reference_hier.level(ilvl, time=t).patches
+    for ip, current_patch in enumerate(ref_patches):
         layout = current_patch.layout
         patch_datas = extract_patchdatas(hierarchies, ilvl, t, ip)
         new_patch_datas = new_patchdatas_from(
             compute, patch_datas, layout, id=current_patch.id, **kwargs
         )
-        new_patches.append(Patch(new_patch_datas, current_patch.id))
+        new_patches.append(
+            Patch(new_patch_datas, current_patch.id, attrs=current_patch.attrs)
+        )
     return new_patches