em diags compatibility

Author: UCaromel

pyphare/pyphare/core/gridlayout.py

@@ -43,6 +43,18 @@
         "Pyy": "primal",
         "Pyz": "primal",
         "Pzz": "primal",
+
+        # finite volume mhd quantities are 3ple dual
+        "mhdRho": "dual",
+        "mhdVx": "dual",
+        "mhdVy": "dual",
+        "mhdVz": "dual",
+        "mhdP": "dual",
+        "mhdRhoVx": "dual",
+        "mhdRhoVy": "dual",
+        "mhdRhoVz": "dual",
+        "mhdEtot": "dual",
+
         "tags": "dual",
     },
     "y": {
@@ -78,6 +90,17 @@
         "Pyy": "primal",
         "Pyz": "primal",
         "Pzz": "primal",
+
+        "mhdRho": "dual",
+        "mhdVx": "dual",
+        "mhdVy": "dual",
+        "mhdVz": "dual",
+        "mhdP": "dual",
+        "mhdRhoVx": "dual",
+        "mhdRhoVy": "dual",
+        "mhdRhoVz": "dual",
+        "mhdEtot": "dual",
+
         "tags": "dual",
     },
     "z": {
@@ -113,6 +136,17 @@
         "Pyy": "primal",
         "Pyz": "primal",
         "Pzz": "primal",
+
+        "mhdRho": "dual",
+        "mhdVx": "dual",
+        "mhdVy": "dual",
+        "mhdVz": "dual",
+        "mhdP": "dual",
+        "mhdRhoVx": "dual",
+        "mhdRhoVy": "dual",
+        "mhdRhoVz": "dual",
+        "mhdEtot": "dual",
+
         "tags": "dual",
     },
 }

pyphare/pyphare/pharein/diagnostics.py

@@ -180,7 +180,7 @@ def _setSubTypeAttributes(self, **kwargs):  # stop pyline complaining
 
 # ------------------------------------------------------------------------------
 class MHDDiagnostics(Diagnostics):
-    mhd_quantities = ["rho", "V", "B", "P", "rhoV", "Etot"]
+    mhd_quantities = ["rho", "V", "P", "rhoV", "Etot"]
     type = "mhd"
 
     def __init__(self, **kwargs):
@@ -196,8 +196,6 @@ def _setSubTypeAttributes(self, **kwargs):
                 MHDDiagnostics.mhd_quantities
             )
             raise ValueError(error_msg.format(kwargs["quantity"]))
-        elif kwargs["quantity"] == "B":
-            self.quantity = "/EM_B"
         else:
             self.quantity = "/mhd/" + kwargs["quantity"]
 

pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py

@@ -1,19 +1,17 @@
-from dataclasses import dataclass, field
 from copy import deepcopy
-import numpy as np
-
+from dataclasses import dataclass, field
 from typing import Any, List, Tuple
 
+import numpy as np
+from pyphare.core import phare_utilities as phut
+
+from ...core.box import Box
+from ...core.gridlayout import GridLayout
+from ...core.phare_utilities import listify, refinement_ratio
 from .hierarchy import PatchHierarchy, format_timestamp
+from .patch import Patch
 from .patchdata import FieldData, ParticleData
 from .patchlevel import PatchLevel
-from .patch import Patch
-from ...core.box import Box
-from ...core.gridlayout import GridLayout
-from ...core.phare_utilities import listify
-from ...core.phare_utilities import refinement_ratio
-from pyphare.core import phare_utilities as phut
-
 
 field_qties = {
     "EM_B_x": "Bx",
@@ -36,6 +34,18 @@
     "momentum_tensor_zz": "Mzz",
     "density": "rho",
     "mass_density": "rho",
+
+    # for now mhd specific quantities
+    "rho": "mhdRho",
+    "V_x": "mhdVx",
+    "V_y": "mhdVy",
+    "V_z": "mhdVz",
+    "P": "mhdP",
+    "rhoV_x": "mhdRhoVx",
+    "rhoV_y": "mhdRhoVy",
+    "rhoV_z": "mhdRhoVz",
+    "Etot": "mhdEtot",
+
     "tags": "tags",
 }
 

pyphare/pyphare/pharesee/run/run.py

@@ -1,25 +1,16 @@
-import os
 import glob
-import numpy as np
-
-from pyphare.pharesee.hierarchy import hierarchy_from
-from pyphare.pharesee.hierarchy import ScalarField, VectorField
+import os
 
-from pyphare.pharesee.hierarchy.hierarchy_utils import compute_hier_from
-from pyphare.pharesee.hierarchy.hierarchy_utils import flat_finest_field
+import numpy as np
 from pyphare.core.phare_utilities import listify
-
 from pyphare.logger import getLogger
-from .utils import (
-    _compute_to_primal,
-    _compute_pop_pressure,
-    _compute_pressure,
-    _compute_current,
-    _compute_divB,
-    _get_rank,
-    make_interpolator,
-)
+from pyphare.pharesee.hierarchy import ScalarField, VectorField, hierarchy_from
+from pyphare.pharesee.hierarchy.hierarchy_utils import (compute_hier_from,
+                                                        flat_finest_field)
 
+from .utils import (_compute_current, _compute_divB, _compute_pop_pressure,
+                    _compute_pressure, _compute_to_primal, _get_rank,
+                    make_interpolator)
 
 logger = getLogger(__name__)
 
@@ -176,6 +167,64 @@ def GetDivB(self, time, merged=False, interp="nearest", **kwargs):
         db = compute_hier_from(_compute_divB, B)
         return ScalarField(self._get(db, time, merged, interp))
 
+    def GetMHDrho(
+        self, time, merged=False, interp="nearest", all_primal=True, **kwargs
+    ):
+        if merged:
+            all_primal = False
+        hier = self._get_hierarchy(time, "mhd_rho.h5", **kwargs)
+        if not all_primal:
+            return self._get(hier, time, merged, interp)
+
+        h = compute_hier_from(_compute_to_primal, hier, value="mhdRho")
+        return VectorField(h)
+
+    def GetMHDV(self, time, merged=False, interp="nearest", all_primal=True, **kwargs):
+        if merged:
+            all_primal = False
+        hier = self._get_hierarchy(time, "mhd_V.h5", **kwargs)
+        if not all_primal:
+            return self._get(hier, time, merged, interp)
+
+        h = compute_hier_from(_compute_to_primal, hier, x="mhdVx", y="mhdVy", z="mhdVz")
+        return VectorField(h)
+
+    def GetMHDP(
+        self, time, merged=False, interp="nearest", all_primal=True, **kwargs
+    ):
+        if merged:
+            all_primal = False
+        hier = self._get_hierarchy(time, "mhd_P.h5", **kwargs)
+        if not all_primal:
+            return self._get(hier, time, merged, interp)
+
+        h = compute_hier_from(_compute_to_primal, hier, value="mhdP")
+        return VectorField(h)
+
+    def GetMHDrhoV(
+        self, time, merged=False, interp="nearest", all_primal=True, **kwargs
+    ):
+        if merged:
+            all_primal = False
+        hier = self._get_hierarchy(time, "mhd_rhoV.h5", **kwargs)
+        if not all_primal:
+            return self._get(hier, time, merged, interp)
+
+        h = compute_hier_from(
+            _compute_to_primal, hier, x="mhdRhoVx", y="mhdRhoVy", z="mhdRhoVz"
+        )
+        return VectorField(h)
+
+    def GetMHDEtot(self, time, merged=False, interp="nearest", all_primal=True, **kwargs):
+        if merged:
+            all_primal = False
+        hier = self._get_hierarchy(time, "mhd_Etot.h5", **kwargs)
+        if not all_primal:
+            return self._get(hier, time, merged, interp)
+
+        h = compute_hier_from(_compute_to_primal, hier, value="mhdEtot")
+        return VectorField(h)
+
     def GetRanks(self, time, merged=False, interp="nearest", **kwargs):
         """
         returns a hierarchy of MPI ranks

pyphare/pyphare/pharesee/run/utils.py

@@ -1,5 +1,5 @@
-from pyphare.core.gridlayout import yee_centering
 import numpy as np
+from pyphare.core.gridlayout import yee_centering
 
 
 def _current1d(by, bz, xby, xbz):
@@ -184,7 +184,7 @@ def _dpd_to_ppp_domain_slicing(**kwargs):
 
 def _ddp_to_ppp_domain_slicing(**kwargs):
     """
-    return the slicing for (dual,primal,primal) to (primal,primal,primal)
+    return the slicing for (dual,dual,primal) to (primal,primal,primal)
     centering that is the centering of Bz on a Yee grid
     """
 
@@ -270,6 +270,31 @@ def _ppd_to_ppp_domain_slicing(**kwargs):
     else:
         raise RuntimeError("dimension not yet implemented")
 
+def _ddd_to_ppp_domain_slicing(**kwargs):
+    """
+    return the slicing for (dual,dual,dual) to (primal,primal,primal)
+    centering that is the centering of Bz on a Yee grid
+    """
+
+    nb_ghosts = kwargs["nb_ghosts"]
+    ndim = kwargs["ndim"]
+
+    inner, inner_shift_left, inner_shift_right = _inner_slices(nb_ghosts)
+
+    if ndim == 1:
+        inner_all = tuple([inner] * ndim)
+        return inner_all, (inner_shift_left, inner_shift_right)
+    elif ndim == 2:
+        inner_all = tuple([inner] * ndim)
+        return inner_all, (
+            (inner_shift_left, inner_shift_left),
+            (inner_shift_left, inner_shift_right),
+            (inner_shift_right, inner_shift_left),
+            (inner_shift_right, inner_shift_right),
+        )
+    else:
+        raise RuntimeError("dimension not yet implemented")
+
 
 slices_to_primal_ = {
     "primal_primal_primal": _ppp_to_ppp_domain_slicing,
@@ -279,6 +304,7 @@ def _ppd_to_ppp_domain_slicing(**kwargs):
     "dual_primal_primal": _dpp_to_ppp_domain_slicing,
     "primal_dual_primal": _pdp_to_ppp_domain_slicing,
     "primal_primal_dual": _ppd_to_ppp_domain_slicing,
+    "dual_dual_dual": _ddd_to_ppp_domain_slicing,
 }
 
 
@@ -468,7 +494,8 @@ def make_interpolator(data, coords, interp, domain, dl, qty, nbrGhosts):
         finest_coords = (x,)
 
     elif dim == 2:
-        from scipy.interpolate import LinearNDInterpolator, NearestNDInterpolator
+        from scipy.interpolate import (LinearNDInterpolator,
+                                       NearestNDInterpolator)
 
         if interp == "nearest":
             interpolator = NearestNDInterpolator(coords, data)

src/diagnostic/detail/types/mhd.hpp

@@ -79,7 +79,7 @@ void MHDDiagnosticWriter<H5Writer>::compute(DiagnosticProperties& diagnostic)
 
     auto& rho  = modelView.getRho();
     auto& V    = modelView.getV();
-    auto& B    = *modelView.getElectromagFields()[0];
+    auto& B    = modelView.getB();
     auto& P    = modelView.getP();
     auto& rhoV = modelView.getRhoV();
     auto& Etot = modelView.getEtot();

src/diagnostic/diagnostic_model_view.hpp

@@ -190,10 +190,11 @@ class ModelView<Hierarchy, Model, MHDIdentifier> : public BaseModelView<Hierarch
 
     NO_DISCARD VecField& getV() const { return this->model_.state.V; }
 
-    // compatibility with the electromag diagnostic writer
-    NO_DISCARD std::vector<VecField*> getElectromagFields() const
+    NO_DISCARD VecField& getB() const { return this->model_.state.B; }
+
+    NO_DISCARD VecField& getE() const
     {
-        return {&this->model_.state.B};
+        throw std::runtime_error("E not currently available in MHD diagnostics");
     }
 
     NO_DISCARD Field& getP() const { return this->model_.state.P; }

tests/functional/mhd_orszagtang/orszag_tang_merged.py

@@ -80,8 +80,8 @@ def p(x, y):
 
     ph.ElectromagDiagnostics(quantity="B", write_timestamps=timestamps)
 
-    # for quantity in ["rho", "V", "B", "P"]:
-    #     ph.MHDDiagnostics(quantity=quantity, write_timestamps=timestamps)
+    for quantity in ["rho", "V", "P"]:
+        ph.MHDDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
     return sim