diff --git a/pyphare/pyphare/pharein/diagnostics.py b/pyphare/pyphare/pharein/diagnostics.py index 0cded1fe8..fd37b128d 100644 --- a/pyphare/pyphare/pharein/diagnostics.py +++ b/pyphare/pyphare/pharein/diagnostics.py @@ -205,7 +205,13 @@ def population_in_model(population): class FluidDiagnostics_(Diagnostics): - fluid_quantities = ["density", "flux", "bulkVelocity", "momentum_tensor"] + fluid_quantities = [ + "density", + "mass_density", + "flux", + "bulkVelocity", + "momentum_tensor", + ] type = "fluid" def __init__(self, **kwargs): @@ -268,7 +274,7 @@ class FluidDiagnostics: def __init__(self, **kwargs): if kwargs["quantity"] == "pressure_tensor": if for_total_ions(**kwargs): - needed_quantities = ["density", "bulkVelocity", "momentum_tensor"] + needed_quantities = ["mass_density", "bulkVelocity", "momentum_tensor"] else: needed_quantities = ["density", "flux", "momentum_tensor"] diff --git a/pyphare/pyphare/pharesee/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy.py index 4bf6a1b64..626d3b1fd 100644 --- a/pyphare/pyphare/pharesee/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy.py @@ -1169,6 +1169,7 @@ def is_root_lvl(patch_level): "momentum_tensor_yz": "Myz", "momentum_tensor_zz": "Mzz", "density": "rho", + "mass_density": "rho", "tags": "tags", } @@ -1256,7 +1257,60 @@ def compute_hier_from(h, compute): def are_compatible_hierarchies(hierarchies): - pass + return True + + +def extract_patchdatas(hierarchies, ilvl, ipatch): + """ + returns a dict {patchdata_name:patchdata} from a list of hierarchies for patch ipatch at level ilvl + """ + patches = [h.patch_levels[ilvl].patches[ipatch] for h in hierarchies] + patch_datas = { + pdname: pd for p in patches for pdname, pd in list(p.patch_datas.items()) + } + return patch_datas + + +def new_patchdatas_from(compute, patchdatas, layout): + new_patch_datas = {} + datas = compute(patchdatas) + for data in datas: + pd = FieldData(layout, data["name"], data["data"], centering=data["centering"]) + new_patch_datas[data["name"]] = pd + return new_patch_datas + + +def new_patches_from(compute, hierarchies, ilvl): + reference_hier = hierarchies[0] + new_patches = [] + patch_nbr = len(reference_hier.patch_levels[ilvl].patches) + for ip in range(patch_nbr): + current_patch = reference_hier.patch_levels[ilvl].patches[ip] + layout = current_patch.layout + patch_datas = extract_patchdatas(hierarchies, ilvl, ip) + new_patch_datas = new_patchdatas_from(compute, patch_datas, layout) + new_patches.append(Patch(new_patch_datas, current_patch.id)) + return new_patches + + +def compute_hier_from__(compute, hierarchies): + """return a new hierarchy using the callback 'compute' on all patchdatas + of the given hierarchies + """ + if not are_compatible_hierarchies(hierarchies): + raise RuntimeError("hierarchies are not compatible") + + reference_hier = hierarchies[0] + domain_box = reference_hier.domain_box + patch_levels = {} + for ilvl in range(reference_hier.levelNbr()): + patch_levels[ilvl] = PatchLevel( + ilvl, new_patches_from(compute, hierarchies, ilvl) + ) + + assert len(reference_hier.time_hier) == 1 # only single time hierarchies now + t = list(reference_hier.time_hier.keys())[0] + return PatchHierarchy(patch_levels, domain_box, refinement_ratio, time=t) def compute_hier_from_(compute, hierarchies): diff --git a/src/core/data/ions/ions.hpp b/src/core/data/ions/ions.hpp index 9182dabde..53f4e4995 100644 --- a/src/core/data/ions/ions.hpp +++ b/src/core/data/ions/ions.hpp @@ -73,6 +73,23 @@ namespace core throw std::runtime_error("Error - cannot access density data"); } + NO_DISCARD field_type const& massDensity() const + { + if (isUsable()) + return sameMasses_ ? *rho_ : *massDensity_; + else + throw std::runtime_error("Error - cannot access density data"); + } + + + + NO_DISCARD field_type& massDensity() + { + if (isUsable()) + return sameMasses_ ? *rho_ : *massDensity_; + else + throw std::runtime_error("Error - cannot access density data"); + } NO_DISCARD vecfield_type const& velocity() const { return bulkVelocity_; } @@ -195,6 +212,10 @@ namespace core { bool usable = rho_ != nullptr and bulkVelocity_.isUsable() and momentumTensor_.isUsable(); + + // if all populations have the same mass, we don't need the massDensity_ + usable &= (sameMasses_) ? true : massDensity_ != nullptr; + for (auto const& pop : populations_) { usable = usable and pop.isUsable(); @@ -208,6 +229,10 @@ namespace core { bool settable = rho_ == nullptr and bulkVelocity_.isSettable() and momentumTensor_.isSettable(); + + // if all populations have the same mass, we don't need the massDensity_ + settable &= (sameMasses_) ? true : massDensity_ == nullptr; + for (auto const& pop : populations_) { settable = settable and pop.isSettable(); @@ -228,12 +253,15 @@ namespace core typename HybridQuantity::Scalar qty; }; - using MomentProperties = std::array; + using MomentProperties = std::vector; NO_DISCARD MomentProperties getFieldNamesAndQuantities() const { - return {{{densityName(), HybridQuantity::Scalar::rho}, - {massDensityName(), HybridQuantity::Scalar::rho}}}; + if (sameMasses_) + return {{{densityName(), HybridQuantity::Scalar::rho}}}; + else + return {{{densityName(), HybridQuantity::Scalar::rho}, + {massDensityName(), HybridQuantity::Scalar::rho}}}; } @@ -246,6 +274,7 @@ namespace core } else if (bufferName == massDensityName()) { + assert(sameMasses_ == false); massDensity_ = field; } else @@ -284,6 +313,8 @@ namespace core return ss.str(); } + NO_DISCARD bool sameMasses() const { return sameMasses_; } + private: bool allSameMass_() const diff --git a/src/diagnostic/detail/types/fluid.hpp b/src/diagnostic/detail/types/fluid.hpp index 2b20209ca..eb0d2afae 100644 --- a/src/diagnostic/detail/types/fluid.hpp +++ b/src/diagnostic/detail/types/fluid.hpp @@ -100,9 +100,9 @@ void FluidDiagnosticWriter::compute(DiagnosticProperties& diagnostic) // and the former levelGhostPartsNew has been moved to levelGhostPartsOld auto levelGhostParts = core::makeIndexRange(pop.levelGhostParticlesOld()); - interpolator(domainParts, pop_momentum_tensor, layout); - interpolator(patchGhostParts, pop_momentum_tensor, layout); - interpolator(levelGhostParts, pop_momentum_tensor, layout); + interpolator(domainParts, pop_momentum_tensor, layout, pop.mass()); + interpolator(patchGhostParts, pop_momentum_tensor, layout, pop.mass()); + interpolator(levelGhostParts, pop_momentum_tensor, layout, pop.mass()); } ions.computeFullMomentumTensor(); }; @@ -126,9 +126,9 @@ void FluidDiagnosticWriter::compute(DiagnosticProperties& diagnostic) // and the former levelGhostPartsNew has been moved to levelGhostPartsOld auto levelGhostParts = core::makeIndexRange(pop.levelGhostParticlesOld()); - interpolator(domainParts, pop_momentum_tensor, layout); - interpolator(patchGhostParts, pop_momentum_tensor, layout); - interpolator(levelGhostParts, pop_momentum_tensor, layout); + interpolator(domainParts, pop_momentum_tensor, layout, pop.mass()); + interpolator(patchGhostParts, pop_momentum_tensor, layout, pop.mass()); + interpolator(levelGhostParts, pop_momentum_tensor, layout, pop.mass()); }; if (isActiveDiag(diagnostic, tree, "momentum_tensor")) { @@ -150,7 +150,8 @@ void FluidDiagnosticWriter::createFiles(DiagnosticProperties& diagnost } std::string tree{"/ions/"}; - checkCreateFileFor_(diagnostic, fileData_, tree, "density", "bulkVelocity", "momentum_tensor"); + checkCreateFileFor_(diagnostic, fileData_, tree, "density", "mass_density", "bulkVelocity", + "momentum_tensor"); } @@ -207,6 +208,8 @@ void FluidDiagnosticWriter::getDataSetInfo(DiagnosticProperties& diagn std::string tree{"/ions/"}; if (isActiveDiag(diagnostic, tree, "density")) infoDS(ions.density(), "density", patchAttributes[lvlPatchID]["ion"]); + if (isActiveDiag(diagnostic, tree, "mass_density")) + infoDS(ions.massDensity(), "mass_density", patchAttributes[lvlPatchID]["ion"]); if (isActiveDiag(diagnostic, tree, "bulkVelocity")) infoVF(ions.velocity(), "bulkVelocity", patchAttributes[lvlPatchID]["ion"]); if (isActiveDiag(diagnostic, tree, "momentum_tensor")) @@ -274,6 +277,8 @@ void FluidDiagnosticWriter::initDataSets( std::string tree{"/ions/"}; if (isActiveDiag(diagnostic, tree, "density")) initDS(path, attr["ion"], "density", null); + if (isActiveDiag(diagnostic, tree, "mass_density")) + initDS(path, attr["ion"], "mass_density", null); if (isActiveDiag(diagnostic, tree, "bulkVelocity")) initVF(path, attr["ion"], "bulkVelocity", null); if (isActiveDiag(diagnostic, tree, "momentum_tensor")) @@ -310,9 +315,10 @@ void FluidDiagnosticWriter::write(DiagnosticProperties& diagnostic) } std::string tree{"/ions/"}; - auto& density = ions.density(); if (isActiveDiag(diagnostic, tree, "density")) - writeDS(path + "density", density); + writeDS(path + "density", ions.density()); + if (isActiveDiag(diagnostic, tree, "mass_density")) + writeDS(path + "mass_density", ions.massDensity()); if (isActiveDiag(diagnostic, tree, "bulkVelocity")) writeTF(path + "bulkVelocity", ions.velocity()); if (isActiveDiag(diagnostic, tree, "momentum_tensor")) diff --git a/tests/core/numerics/ion_updater/test_updater.cpp b/tests/core/numerics/ion_updater/test_updater.cpp index 9f757f853..0c55ff92d 100644 --- a/tests/core/numerics/ion_updater/test_updater.cpp +++ b/tests/core/numerics/ion_updater/test_updater.cpp @@ -466,7 +466,8 @@ struct IonsBuffers void setBuffers(Ions& ions) { ions.setBuffer("rho", &ionDensity); - ions.setBuffer("massDensity", &ionMassDensity); + if (!ions.sameMasses()) + ions.setBuffer("massDensity", &ionMassDensity); auto& v = ions.velocity(); v.setBuffer("bulkVel_x", &Vx); v.setBuffer("bulkVel_y", &Vy);