diff --git a/pyphare/pyphare/pharein/diagnostics.py b/pyphare/pyphare/pharein/diagnostics.py index 3c999aa26..3bd164232 100644 --- a/pyphare/pyphare/pharein/diagnostics.py +++ b/pyphare/pyphare/pharein/diagnostics.py @@ -291,7 +291,7 @@ def __init__(self, **kwargs): if for_total_ions(**kwargs): needed_quantities = ["mass_density", "bulkVelocity", "momentum_tensor"] else: - needed_quantities = ["charge_density", "flux", "momentum_tensor"] + needed_quantities = ["density", "flux", "momentum_tensor"] for quantity in needed_quantities: kwargs["quantity"] = quantity diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py index 2e3b62df5..50c55325d 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py @@ -353,6 +353,8 @@ def box_to_Rectangle(self, box): return Rectangle(box.lower, *box.shape) def plot_2d_patches(self, ilvl, collections, **kwargs): + from matplotlib.patches import Rectangle + if isinstance(collections, list) and all( [isinstance(el, Box) for el in collections] ): @@ -363,19 +365,37 @@ def plot_2d_patches(self, ilvl, collections, **kwargs): level_domain_box = self.level_domain_box(ilvl) mi, ma = level_domain_box.lower.min(), level_domain_box.upper.max() - fig, ax = kwargs.get("subplot", plt.subplots(figsize=(6, 6))) + fig, ax = kwargs.get("subplot", plt.subplots(figsize=(16, 16))) + + color = 1 + i0, j0 = level_domain_box.lower + i1, j1 = level_domain_box.upper + ij = np.zeros((i1 - i0 + 1, j1 - j0 + 1)) + np.nan + ix = np.arange(i0, i1 + 1) + iy = np.arange(j0, j1 + 1) for collection in collections: - facecolor = collection.get("facecolor", "none") - edgecolor = collection.get("edgecolor", "purple") - alpha = collection.get("alpha", 1) - rects = [self.box_to_Rectangle(box) for box in collection["boxes"]] - - ax.add_collection( - PatchCollection( - rects, facecolor=facecolor, alpha=alpha, edgecolor=edgecolor - ) - ) + facecolor = collection.get("facecolor", np.nan) + for box in collection["boxes"]: + if isinstance(box, Box): + i0, j0 = box.lower + i1, j1 = box.upper + ij[i0 : i1 + 1, j0 : j1 + 1] = facecolor + else: + ij[box] = facecolor + + ax.pcolormesh(ix, iy, ij.T, edgecolors="k", cmap="jet") + ax.set_xticks(ix) + ax.set_yticks(iy) + + for patch in self.level(ilvl).patches: + box = patch.box + r = Rectangle(box.lower - 0.5, *(box.upper + 0.5)) + + r.set_edgecolor("r") + r.set_facecolor("none") + r.set_linewidth(2) + ax.add_patch(r) if "title" in kwargs: from textwrap import wrap @@ -383,15 +403,10 @@ def plot_2d_patches(self, ilvl, collections, **kwargs): xfigsize = int(fig.get_size_inches()[0] * 10) # 10 characters per inch ax.set_title("\n".join(wrap(kwargs["title"], xfigsize))) - major_ticks = np.arange(mi - 5, ma + 5 + 5, 5) - ax.set_xticks(major_ticks) - ax.set_yticks(major_ticks) - - minor_ticks = np.arange(mi - 5, ma + 5 + 5, 1) - ax.set_xticks(minor_ticks, minor=True) - ax.set_yticks(minor_ticks, minor=True) - - ax.grid(which="both") + if "xlim" in kwargs: + ax.set_xlim(kwargs["xlim"]) + if "ylim" in kwargs: + ax.set_ylim(kwargs["ylim"]) return fig diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py index 2efadfc3d..85f67f265 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py @@ -28,6 +28,7 @@ "momentum_tensor_zz": "Mzz", "density": "rho", "mass_density": "rho", + "charge_density": "rho", "tags": "tags", } diff --git a/src/amr/level_initializer/hybrid_level_initializer.hpp b/src/amr/level_initializer/hybrid_level_initializer.hpp index 0042f1548..3e112a34e 100644 --- a/src/amr/level_initializer/hybrid_level_initializer.hpp +++ b/src/amr/level_initializer/hybrid_level_initializer.hpp @@ -98,7 +98,7 @@ namespace solver core::depositParticles(ions, layout, interpolate_, core::LevelGhostDeposit{}); } - ions.computeParticleDensity(); // TODO ouam : do we need here the charge density + ions.computeChargeDensity(); ions.computeBulkVelocity(); } diff --git a/src/amr/physical_models/hybrid_model.hpp b/src/amr/physical_models/hybrid_model.hpp index c008161d2..6076836f5 100644 --- a/src/amr/physical_models/hybrid_model.hpp +++ b/src/amr/physical_models/hybrid_model.hpp @@ -154,7 +154,7 @@ void HybridModel::f hybridInfo.modelMagnetic = core::VecFieldNames{state.electromag.B}; hybridInfo.modelElectric = core::VecFieldNames{state.electromag.E}; - hybridInfo.modelIonDensity = state.ions.particleDensityName(); + hybridInfo.modelIonDensity = state.ions.chargeDensityName(); hybridInfo.modelIonBulkVelocity = core::VecFieldNames{state.ions.velocity()}; hybridInfo.modelCurrent = core::VecFieldNames{state.J}; diff --git a/src/core/data/electrons/electrons.hpp b/src/core/data/electrons/electrons.hpp index c3b573ec3..b18d0c290 100644 --- a/src/core/data/electrons/electrons.hpp +++ b/src/core/data/electrons/electrons.hpp @@ -65,7 +65,7 @@ class StandardHybridElectronFluxComputer { if (isUsable()) { - return ions_.density(); + return ions_.chargeDensity(); } else { @@ -78,7 +78,7 @@ class StandardHybridElectronFluxComputer { if (isUsable()) { - return ions_.density(); + return ions_.chargeDensity(); } else { @@ -102,7 +102,7 @@ class StandardHybridElectronFluxComputer } - void computeParticleDensity() {} // TODO ouam : which is the same as the charge particle... except the sign + void computeDensity() {} void computeBulkVelocity(GridLayout const& layout) { @@ -112,23 +112,23 @@ class StandardHybridElectronFluxComputer auto const& Vix = ions_.velocity()(Component::X); auto const& Viy = ions_.velocity()(Component::Y); auto const& Viz = ions_.velocity()(Component::Z); - auto const& Ni = ions_.density(); + auto const& Ne = ions_.chargeDensity(); auto& Vex = Ve_(Component::X); auto& Vey = Ve_(Component::Y); auto& Vez = Ve_(Component::Z); // from Ni because all components defined on primal - layout.evalOnBox(Ni, [&](auto const&... args) { + layout.evalOnBox(Ne, [&](auto const&... args) { auto arr = std::array{args...}; auto const JxOnVx = GridLayout::project(Jx, arr, GridLayout::JxToMoments()); auto const JyOnVy = GridLayout::project(Jy, arr, GridLayout::JyToMoments()); auto const JzOnVz = GridLayout::project(Jz, arr, GridLayout::JzToMoments()); - Vex(arr) = Vix(arr) - JxOnVx / Ni(arr); - Vey(arr) = Viy(arr) - JyOnVy / Ni(arr); - Vez(arr) = Viz(arr) - JzOnVz / Ni(arr); + Vex(arr) = Vix(arr) - JxOnVx / Ne(arr); + Vey(arr) = Viy(arr) - JyOnVy / Ne(arr); + Vez(arr) = Viz(arr) - JzOnVz / Ne(arr); }); } @@ -211,7 +211,7 @@ class IsothermalElectronPressureClosure if (!Pe_.isUsable()) throw std::runtime_error("Error - isothermal closure pressure not usable"); - auto const& Ne_ = ions_.density(); + auto const& Ne_ = ions_.chargeDensity(); std::transform(std::begin(Ne_), std::end(Ne_), std::begin(Pe_), [this](auto n) { return n * Te_; }); } @@ -281,7 +281,7 @@ class ElectronMomentModel - void computeParticleDensity() { fluxComput_.computeParticleDensity(); } + void computeDensity() { fluxComput_.computeDensity(); } void computeBulkVelocity(GridLayout const& layout) { fluxComput_.computeBulkVelocity(layout); } void computePressure(GridLayout const& layout) { pressureClosure_.computePressure(layout); } @@ -310,7 +310,7 @@ class Electrons : public LayoutHolder { if (isUsable()) { - momentModel_.computeParticleDensity(); + momentModel_.computeDensity(); momentModel_.computeBulkVelocity(layout); momentModel_.computePressure(layout); } diff --git a/src/core/data/ions/ion_population/ion_population.hpp b/src/core/data/ions/ion_population/ion_population.hpp index 4cfa05e38..d6c5ab4df 100644 --- a/src/core/data/ions/ion_population/ion_population.hpp +++ b/src/core/data/ions/ion_population/ion_population.hpp @@ -50,13 +50,14 @@ namespace core NO_DISCARD bool isUsable() const { - return particles_.isUsable() && particleDensity_.isUsable() && chargeDensity_.isUsable() && flux_.isUsable() - && momentumTensor_.isUsable(); + return particles_.isUsable() && particleDensity_.isUsable() && chargeDensity_.isUsable() + && flux_.isUsable() && momentumTensor_.isUsable(); } NO_DISCARD bool isSettable() const { - return particles_.isSettable() && particleDensity_.isSettable() && chargeDensity_.isSettable() && flux_.isSettable() + return particles_.isSettable() && particleDensity_.isSettable() + && chargeDensity_.isSettable() && flux_.isSettable() && momentumTensor_.isSettable(); } @@ -81,9 +82,6 @@ namespace core return particles_.levelGhostParticlesNew(); } - NO_DISCARD field_type const& density() const { return particleDensity_; } // TODO ouam : to remove - NO_DISCARD field_type& density() { return particleDensity_; } // TODO ouam : to remove - NO_DISCARD field_type const& particleDensity() const { return particleDensity_; } NO_DISCARD field_type& particleDensity() { return particleDensity_; } @@ -107,7 +105,8 @@ namespace core NO_DISCARD auto getCompileTimeResourcesViewList() { - return std::forward_as_tuple(flux_, momentumTensor_, particleDensity_, chargeDensity_, particles_); + return std::forward_as_tuple(flux_, momentumTensor_, particleDensity_, chargeDensity_, + particles_); } diff --git a/src/core/data/ions/ions.hpp b/src/core/data/ions/ions.hpp index 9f0133059..a9a536802 100644 --- a/src/core/data/ions/ions.hpp +++ b/src/core/data/ions/ions.hpp @@ -41,8 +41,7 @@ namespace core explicit Ions(PHARE::initializer::PHAREDict const& dict) - : particleDensity_{particleDensityName(), HybridQuantity::Scalar::rho} - , massDensity_{massDensityName(), HybridQuantity::Scalar::rho} + : massDensity_{massDensityName(), HybridQuantity::Scalar::rho} , chargeDensity_{chargeDensityName(), HybridQuantity::Scalar::rho} , bulkVelocity_{"bulkVel", HybridQuantity::Vector::V} , populations_{generate( @@ -50,7 +49,6 @@ namespace core return IonPopulation{dict["pop" + std::to_string(ipop)]}; }, dict["nbrPopulations"].template to())} - , sameMasses_{allSameMass_()} , momentumTensor_{"momentumTensor", HybridQuantity::Tensor::M} { } @@ -59,18 +57,8 @@ namespace core NO_DISCARD auto nbrPopulations() const { return populations_.size(); } NO_DISCARD auto size() const { return nbrPopulations(); } - - NO_DISCARD field_type const& density() const { return particleDensity_; } // TODO ouam : to remove - NO_DISCARD field_type& density() { return particleDensity_; } // TODO ouam : to remove - - NO_DISCARD field_type const& particleDensity() const { return particleDensity_; } - NO_DISCARD field_type& particleDensity() { return particleDensity_; } - - NO_DISCARD field_type const& massDensity() const - { - return sameMasses_ ? particleDensity_ : massDensity_; - } - NO_DISCARD field_type& massDensity() { return sameMasses_ ? particleDensity_ : massDensity_; } + NO_DISCARD field_type const& massDensity() const { return massDensity_; } + NO_DISCARD field_type const& massDensity() { return massDensity_; } NO_DISCARD field_type const& chargeDensity() const { return chargeDensity_; } NO_DISCARD field_type& chargeDensity() { return chargeDensity_; } @@ -78,16 +66,15 @@ namespace core NO_DISCARD vecfield_type const& velocity() const { return bulkVelocity_; } NO_DISCARD vecfield_type& velocity() { return bulkVelocity_; } - NO_DISCARD std::string static particleDensityName() { return "particleDensity"; } NO_DISCARD std::string static chargeDensityName() { return "chargeDensity"; } NO_DISCARD std::string static massDensityName() { return "massDensity"; } tensorfield_type const& momentumTensor() const { return momentumTensor_; } tensorfield_type& momentumTensor() { return momentumTensor_; } - void computeParticleDensity() + void computeChargeDensity() { - particleDensity_.zero(); + chargeDensity_.zero(); for (auto const& pop : populations_) { @@ -95,11 +82,13 @@ namespace core // nodes. This is more efficient and easier to code as we don't // have to account for the field dimensionality. - auto& popDensity = pop.particleDensity(); - std::transform(std::begin(particleDensity_), std::end(particleDensity_), std::begin(popDensity), - std::begin(particleDensity_), std::plus{}); + auto& popDensity = pop.chargeDensity(); + std::transform(std::begin(chargeDensity_), std::end(chargeDensity_), + std::begin(popDensity), std::begin(chargeDensity_), + std::plus{}); } } + void computeMassDensity() { massDensity_.zero(); @@ -121,14 +110,7 @@ namespace core void computeBulkVelocity() { - // the bulk velocity is sum(pop_mass * pop_flux) / sum(pop_mass * pop_density) - // if all populations have the same mass, this is equivalent to sum(pop_flux) / - // sum(pop_density) sum(pop_density) is particleDensity_ and already known by the time we get here. - // sum(pop_mass * pop_flux) is massDensity_ and is computed by computeMassDensity() if - // needed - if (!sameMasses_) - computeMassDensity(); - auto const& density = (sameMasses_) ? particleDensity_ : massDensity_; + computeMassDensity(); bulkVelocity_.zero(); auto& vx = bulkVelocity_.getComponent(Component::X); @@ -138,28 +120,26 @@ namespace core for (auto& pop : populations_) { // account for mass only if populations have different masses - std::function plus = std::plus{}; std::function plusMass = [&pop](Float const& v, Float const& f) { return v + f * pop.mass(); }; auto const& flux = pop.flux(); auto&& [fx, fy, fz] = flux(); - std::transform(std::begin(vx), std::end(vx), std::begin(fx), std::begin(vx), - sameMasses_ ? plus : plusMass); + plusMass); std::transform(std::begin(vy), std::end(vy), std::begin(fy), std::begin(vy), - sameMasses_ ? plus : plusMass); + plusMass); std::transform(std::begin(vz), std::end(vz), std::begin(fz), std::begin(vz), - sameMasses_ ? plus : plusMass); + plusMass); } - std::transform(std::begin(vx), std::end(vx), std::begin(density), std::begin(vx), + std::transform(std::begin(vx), std::end(vx), std::begin(massDensity_), std::begin(vx), std::divides{}); - std::transform(std::begin(vy), std::end(vy), std::begin(density), std::begin(vy), + std::transform(std::begin(vy), std::end(vy), std::begin(massDensity_), std::begin(vy), std::divides{}); - std::transform(std::begin(vz), std::end(vz), std::begin(density), std::begin(vz), + std::transform(std::begin(vz), std::end(vz), std::begin(massDensity_), std::begin(vz), std::divides{}); } @@ -192,11 +172,8 @@ namespace core // because it is for internal use only so no object will ever need to access it. NO_DISCARD bool isUsable() const { - bool usable - = particleDensity_.isUsable() and chargeDensity_.isUsable() and bulkVelocity_.isUsable() and momentumTensor_.isUsable(); - - // if all populations have the same mass, we don't need the massDensity_ - usable &= (sameMasses_) ? true : massDensity_.isUsable(); + bool usable = chargeDensity_.isUsable() and bulkVelocity_.isUsable() + and momentumTensor_.isUsable() and massDensity_.isUsable(); for (auto const& pop : populations_) { @@ -209,11 +186,8 @@ namespace core NO_DISCARD bool isSettable() const { - bool settable - = particleDensity_.isSettable() and chargeDensity_.isSettable() and bulkVelocity_.isSettable() and momentumTensor_.isSettable(); - - // if all populations have the same mass, we don't need the massDensity_ - settable &= (sameMasses_) ? true : massDensity_.isSettable(); + bool settable = massDensity_.isSettable() and chargeDensity_.isSettable() + and bulkVelocity_.isSettable() and momentumTensor_.isSettable(); for (auto const& pop : populations_) { @@ -237,7 +211,8 @@ namespace core NO_DISCARD auto getCompileTimeResourcesViewList() { - return std::forward_as_tuple(bulkVelocity_, momentumTensor_, particleDensity_, chargeDensity_, massDensity_); + return std::forward_as_tuple(bulkVelocity_, momentumTensor_, chargeDensity_, + massDensity_); } @@ -258,30 +233,12 @@ namespace core return ss.str(); } - NO_DISCARD bool sameMasses() const { return sameMasses_; } - private: - bool allSameMass_() const - { - return all(populations_, [this](auto const& pop) { // arbitrary small diff - return float_equals(pop.mass(), populations_.front().mass(), /*abs_tol=*/1e-10); - }); - } - - - - - field_type particleDensity_; field_type massDensity_; field_type chargeDensity_; vecfield_type bulkVelocity_; std::vector populations_; - - // note this is set at construction when all populations are added - // in the future if some populations are dynamically created during the simulation - // this should be updated accordingly - bool sameMasses_{false}; tensorfield_type momentumTensor_; }; } // namespace core diff --git a/src/core/numerics/interpolator/interpolator.hpp b/src/core/numerics/interpolator/interpolator.hpp index 2aa5f4a86..e8858d6a4 100644 --- a/src/core/numerics/interpolator/interpolator.hpp +++ b/src/core/numerics/interpolator/interpolator.hpp @@ -469,8 +469,9 @@ class Interpolator : private Weighter * onto the particle. */ template - inline void operator()(ParticleRange& particleRange, Field& particleDensity, Field& chargeDensity, VecField& flux, - GridLayout const& layout, double coef = 1.) + inline void operator()(ParticleRange& particleRange, Field& particleDensity, + Field& chargeDensity, VecField& flux, GridLayout const& layout, + double coef = 1.) { auto begin = particleRange.begin(); auto end = particleRange.end(); @@ -487,11 +488,11 @@ class Interpolator : private Weighter currPart->delta); particleToMesh_( - particleDensity, *currPart, [](auto const& part) { return 1.; }, startIndex_, weights_, - coef); + particleDensity, *currPart, [](auto const& part) { return 1.; }, startIndex_, + weights_, coef); particleToMesh_( - chargeDensity, *currPart, [](auto const& part) { return part.charge; }, startIndex_, weights_, - coef); + chargeDensity, *currPart, [](auto const& part) { return part.charge; }, startIndex_, + weights_, coef); particleToMesh_( xFlux, *currPart, [](auto const& part) { return part.v[0]; }, startIndex_, weights_, coef); @@ -505,8 +506,8 @@ class Interpolator : private Weighter PHARE_LOG_STOP(3, "ParticleToMesh::operator()"); } template - inline void operator()(ParticleRange&& range, Field& particleDensity, Field& chargeDensity, VecField& flux, - GridLayout const& layout, double coef = 1.) + inline void operator()(ParticleRange&& range, Field& particleDensity, Field& chargeDensity, + VecField& flux, GridLayout const& layout, double coef = 1.) { (*this)(range, particleDensity, chargeDensity, flux, layout, coef); } diff --git a/src/core/numerics/ion_updater/ion_updater.hpp b/src/core/numerics/ion_updater/ion_updater.hpp index bf5ca46e3..6f51a6147 100644 --- a/src/core/numerics/ion_updater/ion_updater.hpp +++ b/src/core/numerics/ion_updater/ion_updater.hpp @@ -107,8 +107,7 @@ void IonUpdater::updatePopulations(Ions& ions, Ele template void IonUpdater::updateIons(Ions& ions) { - ions.computeParticleDensity(); // TODO ouam : should we need here to compute the charge density ? - // ions.computeChargeDensity(); // TODO ouam : this one might be better for the diagnostics + ions.computeChargeDensity(); ions.computeBulkVelocity(); } @@ -182,7 +181,8 @@ void IonUpdater::updateAndDepositDomain_(Ions& ion auto enteredInDomain = pusher_->move(inRange, outRange, em, pop.mass(), interpolator_, layout, inGhostBox, inDomainBox); - interpolator_(enteredInDomain, pop.particleDensity(), pop.chargeDensity(), pop.flux(), layout); + interpolator_(enteredInDomain, pop.particleDensity(), pop.chargeDensity(), pop.flux(), + layout); if (copyInDomain) { @@ -270,7 +270,8 @@ void IonUpdater::updateAndDepositAll_(Ions& ions, pushAndCopyInDomain(makeIndexRange(pop.patchGhostParticles())); pushAndCopyInDomain(makeIndexRange(pop.levelGhostParticles())); - interpolator_(makeIndexRange(domainParticles), pop.particleDensity(), pop.chargeDensity(), pop.flux(), layout); + interpolator_(makeIndexRange(domainParticles), pop.particleDensity(), pop.chargeDensity(), + pop.flux(), layout); } } diff --git a/src/core/numerics/moments/moments.hpp b/src/core/numerics/moments/moments.hpp index f28a75215..87b805fa1 100644 --- a/src/core/numerics/moments/moments.hpp +++ b/src/core/numerics/moments/moments.hpp @@ -16,6 +16,7 @@ namespace core for (auto& pop : ions) { pop.particleDensity().zero(); + pop.chargeDensity().zero(); pop.flux().zero(); } } @@ -41,8 +42,8 @@ namespace core for (auto& pop : ions) { auto& particleDensity = pop.particleDensity(); - auto& chargeDensity = pop.chargeDensity(); - auto& flux = pop.flux(); + auto& chargeDensity = pop.chargeDensity(); + auto& flux = pop.flux(); if constexpr (std::is_same_v) { diff --git a/src/diagnostic/detail/types/fluid.hpp b/src/diagnostic/detail/types/fluid.hpp index da0b5d7d3..89839c5d7 100644 --- a/src/diagnostic/detail/types/fluid.hpp +++ b/src/diagnostic/detail/types/fluid.hpp @@ -146,12 +146,13 @@ void FluidDiagnosticWriter::createFiles(DiagnosticProperties& diagnost for (auto const& pop : this->h5Writer_.modelView().getIons()) { std::string tree{"/ions/pop/" + pop.name() + "/"}; - checkCreateFileFor_(diagnostic, fileData_, tree, "charge_density", "flux", "momentum_tensor"); + checkCreateFileFor_(diagnostic, fileData_, tree, "density", "charge_density", "flux", + "momentum_tensor"); } std::string tree{"/ions/"}; - checkCreateFileFor_(diagnostic, fileData_, tree, "density", "mass_density", "bulkVelocity", - "momentum_tensor"); + checkCreateFileFor_(diagnostic, fileData_, tree, "charge_density", "mass_density", + "bulkVelocity", "momentum_tensor"); } @@ -198,7 +199,9 @@ void FluidDiagnosticWriter::getDataSetInfo(DiagnosticProperties& diagn std::string tree{"/ions/pop/" + pop.name() + "/"}; auto& popAttr = patchAttributes[lvlPatchID]["fluid_" + pop.name()]; if (isActiveDiag(diagnostic, tree, "density")) - infoDS(pop.chargeDensity(), "density", popAttr); + infoDS(pop.particleDensity(), "density", popAttr); + if (isActiveDiag(diagnostic, tree, "charge_density")) + infoDS(pop.chargeDensity(), "charge_density", popAttr); if (isActiveDiag(diagnostic, tree, "flux")) infoVF(pop.flux(), "flux", popAttr); if (isActiveDiag(diagnostic, tree, "momentum_tensor")) @@ -268,6 +271,8 @@ void FluidDiagnosticWriter::initDataSets( std::string popPath(path + "pop/" + pop.name() + "/"); if (isActiveDiag(diagnostic, tree, "density")) initDS(path, attr[popId], "density", null); + if (isActiveDiag(diagnostic, tree, "charge_density")) + initDS(path, attr[popId], "charge_density", null); if (isActiveDiag(diagnostic, tree, "flux")) initVF(path, attr[popId], "flux", null); if (isActiveDiag(diagnostic, tree, "momentum_tensor")) @@ -307,7 +312,9 @@ void FluidDiagnosticWriter::write(DiagnosticProperties& diagnostic) { std::string tree{"/ions/pop/" + pop.name() + "/"}; if (isActiveDiag(diagnostic, tree, "density")) - writeDS(path + "density", pop.chargeDensity()); + writeDS(path + "density", pop.particleDensity()); + if (isActiveDiag(diagnostic, tree, "charge_density")) + writeDS(path + "charge_density", pop.chargeDensity()); if (isActiveDiag(diagnostic, tree, "flux")) writeTF(path + "flux", pop.flux()); if (isActiveDiag(diagnostic, tree, "momentum_tensor")) @@ -344,6 +351,7 @@ void FluidDiagnosticWriter::writeAttributes( for (auto& pop : h5Writer.modelView().getIons()) { std::string tree = "/ions/pop/" + pop.name() + "/"; + checkWrite(tree, "density", pop); checkWrite(tree, "charge_density", pop); checkWrite(tree, "flux", pop); checkWrite(tree, "momentum_tensor", pop); diff --git a/tests/core/data/electrons/test_electrons.cpp b/tests/core/data/electrons/test_electrons.cpp index 421140891..bf0bad9ac 100644 --- a/tests/core/data/electrons/test_electrons.cpp +++ b/tests/core/data/electrons/test_electrons.cpp @@ -126,8 +126,8 @@ class NDlayout template, InterpConst<1>>*/> struct ElectronsTest : public ::testing::Test { - static constexpr auto dim = typename TypeInfo::first_type{}(); - static constexpr auto interp = typename TypeInfo::second_type{}(); + static constexpr auto dim = typename TypeInfo::first_type{}(); + static constexpr auto interp = typename TypeInfo::second_type{}(); using GridYee = GridLayout>; @@ -150,7 +150,7 @@ struct ElectronsTest : public ::testing::Test UsableVecField J, F, Ve, Vi; UsableTensorField ionTensor, protonTensor; - GridND ionParticleDensity, ionChargeDensity, protonParticleDensity, protonChargeDensity, Pe; + GridND ionChargeDensity, ionMassDensity, protonParticleDensity, protonChargeDensity, Pe; ParticleArray_t domainParticles{layout.AMRBox()}; ParticleArray_t patchGhostParticles = domainParticles; @@ -163,13 +163,14 @@ struct ElectronsTest : public ::testing::Test template auto static _ions(Args&... args) { - auto const& [ionFlux, ionParticleDensity, ionChargeDensity, protonParticleDensity, protonChargeDensity, Vi, ionTensor, protonTensor, pack] + auto const& [ionFlux, ionChargeDensity, ionMassDensity, protonParticleDensity, + protonChargeDensity, Vi, ionTensor, protonTensor, pack] = std::forward_as_tuple(args...); IonsT ions{createDict()["ions"]}; { - auto const& [V, m, d, d_c, d_m] = ions.getCompileTimeResourcesViewList(); - d.setBuffer(&ionParticleDensity); + auto const& [V, m, d_c, d_m] = ions.getCompileTimeResourcesViewList(); d_c.setBuffer(&ionChargeDensity); + d_m.setBuffer(&ionMassDensity); Vi.set_on(V); ionTensor.set_on(m); } @@ -194,16 +195,17 @@ struct ElectronsTest : public ::testing::Test , Vi{"bulkVel", layout, HybridQuantity::Vector::V} , ionTensor{"momentumTensor", layout, HybridQuantity::Tensor::M} , protonTensor{"protons_momentumTensor", layout, HybridQuantity::Tensor::M} - , ionParticleDensity{"particleDensity", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} , ionChargeDensity{"chargeDensity", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} + layout.allocSize(HybridQuantity::Scalar::rho)} + , ionMassDensity{"chargeDensity", HybridQuantity::Scalar::rho, + layout.allocSize(HybridQuantity::Scalar::rho)} , protonParticleDensity{"protons_particleDensity", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} + layout.allocSize(HybridQuantity::Scalar::rho)} , protonChargeDensity{"protons_chargeDensity", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} + layout.allocSize(HybridQuantity::Scalar::rho)} , Pe{"Pe", HybridQuantity::Scalar::P, layout.allocSize(HybridQuantity::Scalar::P)} - , ions{_ions(F, ionParticleDensity, ionChargeDensity, protonParticleDensity, protonChargeDensity, Vi, ionTensor, protonTensor, pack)} + , ions{_ions(F, ionChargeDensity, ionMassDensity, protonParticleDensity, + protonChargeDensity, Vi, ionTensor, protonTensor, pack)} , electrons{createDict()["electrons"], ions, J} { auto&& emm = std::get<0>(electrons.getCompileTimeResourcesViewList()); @@ -242,7 +244,7 @@ struct ElectronsTest : public ::testing::Test fill(Jy, [](double x) { return std::sinh(0.3 * x); }); fill(Jz, [](double x) { return std::sinh(0.4 * x); }); - fill(ionParticleDensity, [](double x) { return std::cosh(0.1 * x); }); + fill(ionChargeDensity, [](double x) { return std::cosh(0.1 * x); }); } else if constexpr (dim == 2) { @@ -271,7 +273,7 @@ struct ElectronsTest : public ::testing::Test fill(Jy, [](double x, double y) { return std::sinh(0.3 * x) * std::sinh(0.3 * y); }); fill(Jz, [](double x, double y) { return std::sinh(0.4 * x) * std::sinh(0.4 * y); }); - fill(ionParticleDensity, + fill(ionChargeDensity, [](double x, double y) { return std::cosh(0.1 * x) * std::cosh(0.1 * y); }); } else if constexpr (dim == 3) @@ -318,7 +320,7 @@ struct ElectronsTest : public ::testing::Test return std::sinh(0.4 * x) * std::sinh(0.4 * y) * std::sinh(0.4 * z); }); - fill(ionParticleDensity, [](double x, double y, double z) { + fill(ionChargeDensity, [](double x, double y, double z) { return std::cosh(0.1 * x) * std::cosh(0.1 * y) * std::cosh(0.1 * z); }); } diff --git a/tests/core/numerics/ion_updater/test_updater.cpp b/tests/core/numerics/ion_updater/test_updater.cpp index 622a83c7f..09013e517 100644 --- a/tests/core/numerics/ion_updater/test_updater.cpp +++ b/tests/core/numerics/ion_updater/test_updater.cpp @@ -211,7 +211,6 @@ struct IonsBuffers using ParticleArray = typename PHARETypes::ParticleArray_t; using ParticleInitializerFactory = typename PHARETypes::ParticleInitializerFactory; - Grid ionParticleDensity; Grid ionChargeDensity; Grid ionMassDensity; Grid protonParticleDensity; @@ -240,20 +239,18 @@ struct IonsBuffers ParticlesPack alphaPack; IonsBuffers(GridLayout const& layout) - : ionParticleDensity{"particleDensity", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} - , ionChargeDensity{"chargeDensity", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} + : ionChargeDensity{"chargeDensity", HybridQuantity::Scalar::rho, + layout.allocSize(HybridQuantity::Scalar::rho)} , ionMassDensity{"massDensity", HybridQuantity::Scalar::rho, layout.allocSize(HybridQuantity::Scalar::rho)} , protonParticleDensity{"protons_particleDensity", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} + layout.allocSize(HybridQuantity::Scalar::rho)} , protonChargeDensity{"protons_chargeDensity", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} + layout.allocSize(HybridQuantity::Scalar::rho)} , alphaParticleDensity{"alpha_particleDensity", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} + layout.allocSize(HybridQuantity::Scalar::rho)} , alphaChargeDensity{"alpha_chargeDensity", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} + layout.allocSize(HybridQuantity::Scalar::rho)} , protonF{"protons_flux", layout, HybridQuantity::Vector::V} , alphaF{"alpha_flux", layout, HybridQuantity::Vector::V} , Vi{"bulkVel", layout, HybridQuantity::Vector::V} @@ -279,20 +276,18 @@ struct IonsBuffers IonsBuffers(IonsBuffers const& source, GridLayout const& layout) - : ionParticleDensity{"rho", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} - , ionChargeDensity{"chargeDensity", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} + : ionChargeDensity{"chargeDensity", HybridQuantity::Scalar::rho, + layout.allocSize(HybridQuantity::Scalar::rho)} , ionMassDensity{"massDensity", HybridQuantity::Scalar::rho, layout.allocSize(HybridQuantity::Scalar::rho)} , protonParticleDensity{"protons_particleDensity", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} + layout.allocSize(HybridQuantity::Scalar::rho)} , protonChargeDensity{"protons_chargeDensity", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} + layout.allocSize(HybridQuantity::Scalar::rho)} , alphaParticleDensity{"alpha_particleDensity", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} + layout.allocSize(HybridQuantity::Scalar::rho)} , alphaChargeDensity{"alpha_chargeDensity", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} + layout.allocSize(HybridQuantity::Scalar::rho)} , protonF{"protons_flux", layout, HybridQuantity::Vector::V} , alphaF{"alpha_flux", layout, HybridQuantity::Vector::V} , Vi{"bulkVel", layout, HybridQuantity::Vector::V} @@ -315,7 +310,6 @@ struct IonsBuffers &alphaLevelGhost, &alphaLevelGhostOld, &alphaLevelGhostNew} { - ionParticleDensity.copyData(source.ionParticleDensity); ionChargeDensity.copyData(source.ionChargeDensity); ionMassDensity.copyData(source.ionMassDensity); protonParticleDensity.copyData(source.protonParticleDensity); @@ -331,10 +325,9 @@ struct IonsBuffers void setBuffers(Ions& ions) { { - auto const& [V, m, d, cd, md] = ions.getCompileTimeResourcesViewList(); + auto const& [V, m, cd, md] = ions.getCompileTimeResourcesViewList(); Vi.set_on(V); M.set_on(m); - d.setBuffer(&ionParticleDensity); cd.setBuffer(&ionChargeDensity); md.setBuffer(&ionMassDensity); } @@ -536,7 +529,7 @@ struct IonUpdaterTest : public ::testing::Test } } // end 1D - } // end pop loop + } // end pop loop PHARE::core::depositParticles(ions, layout, Interpolator{}, PHARE::core::DomainDeposit{}); @@ -547,8 +540,7 @@ struct IonUpdaterTest : public ::testing::Test PHARE::core::LevelGhostDeposit{}); - ions.computeParticleDensity(); // TODO ouam : should we need here the charge density - // ions.computeChargeDensity(); // TODO ouam : should be included when ions will have a coimputeChargeDensity + ions.computeChargeDensity(); ions.computeBulkVelocity(); } // end Ctor @@ -561,17 +553,17 @@ struct IonUpdaterTest : public ::testing::Test for (auto& pop : this->ions) { - interpolate(makeIndexRange(pop.patchGhostParticles()), pop.particleDensity(), pop.chargeDensity(), pop.flux(), - layout); + interpolate(makeIndexRange(pop.patchGhostParticles()), pop.particleDensity(), + pop.chargeDensity(), pop.flux(), layout); double alpha = 0.5; - interpolate(makeIndexRange(pop.levelGhostParticlesNew()), pop.particleDensity(), pop.chargeDensity(), pop.flux(), - layout, + interpolate(makeIndexRange(pop.levelGhostParticlesNew()), pop.particleDensity(), + pop.chargeDensity(), pop.flux(), layout, /*coef = */ alpha); - interpolate(makeIndexRange(pop.levelGhostParticlesOld()), pop.particleDensity(), pop.chargeDensity(), pop.flux(), - layout, + interpolate(makeIndexRange(pop.levelGhostParticlesOld()), pop.particleDensity(), + pop.chargeDensity(), pop.flux(), layout, /*coef = */ (1. - alpha)); } } @@ -583,16 +575,16 @@ struct IonUpdaterTest : public ::testing::Test auto& populations = this->ions.getRunTimeResourcesViewList(); auto& protonParticleDensity = populations[0].particleDensity(); - auto& protonChargeDensity = populations[0].chargeDensity(); - auto& protonFx = populations[0].flux().getComponent(Component::X); - auto& protonFy = populations[0].flux().getComponent(Component::Y); - auto& protonFz = populations[0].flux().getComponent(Component::Z); + auto& protonChargeDensity = populations[0].chargeDensity(); + auto& protonFx = populations[0].flux().getComponent(Component::X); + auto& protonFy = populations[0].flux().getComponent(Component::Y); + auto& protonFz = populations[0].flux().getComponent(Component::Z); auto& alphaParticleDensity = populations[1].particleDensity(); - auto& alphaChargeDensity = populations[1].chargeDensity(); - auto& alphaFx = populations[1].flux().getComponent(Component::X); - auto& alphaFy = populations[1].flux().getComponent(Component::Y); - auto& alphaFz = populations[1].flux().getComponent(Component::Z); + auto& alphaChargeDensity = populations[1].chargeDensity(); + auto& alphaFx = populations[1].flux().getComponent(Component::X); + auto& alphaFy = populations[1].flux().getComponent(Component::Y); + auto& alphaFz = populations[1].flux().getComponent(Component::Z); auto ix0 = this->layout.physicalStartIndex(QtyCentering::primal, Direction::X); auto ix1 = this->layout.physicalEndIndex(QtyCentering::primal, Direction::X); @@ -633,7 +625,6 @@ struct IonUpdaterTest : public ::testing::Test check(alphaFy, ionsBufferCpy.alphaF(Component::Y)); check(alphaFz, ionsBufferCpy.alphaF(Component::Z)); - check(ions.particleDensity(), ionsBufferCpy.ionParticleDensity); check(ions.velocity().getComponent(Component::X), ionsBufferCpy.Vi(Component::X)); check(ions.velocity().getComponent(Component::Y), ionsBufferCpy.Vi(Component::Y)); check(ions.velocity().getComponent(Component::Z), ionsBufferCpy.Vi(Component::Z)); @@ -674,7 +665,7 @@ struct IonUpdaterTest : public ::testing::Test } }; - auto& populations = this->ions.getRunTimeResourcesViewList(); + auto& populations = this->ions.getRunTimeResourcesViewList(); auto& protonParticleDensity = populations[0].particleDensity(); auto& alphaParticleDensity = populations[1].particleDensity(); diff --git a/tests/diagnostic/__init__.py b/tests/diagnostic/__init__.py index 0bf944d61..370bfe7a3 100644 --- a/tests/diagnostic/__init__.py +++ b/tests/diagnostic/__init__.py @@ -28,7 +28,12 @@ def dump_all_diags(pops=[], flush_every=100, timestamps=None): # write_timestamps=timestamps # ) - for quantity in ["charge_density", "bulkVelocity", "pressure_tensor"]: + for quantity in [ + "charge_density", + "mass_density", + "bulkVelocity", + "pressure_tensor", + ]: ph.FluidDiagnostics( quantity=quantity, write_timestamps=timestamps, @@ -36,7 +41,7 @@ def dump_all_diags(pops=[], flush_every=100, timestamps=None): ) for pop in pops: - for quantity in ["density", "flux", "pressure_tensor"]: + for quantity in ["density", "charge_density", "flux", "pressure_tensor"]: ph.FluidDiagnostics( quantity=quantity, write_timestamps=timestamps, diff --git a/tests/diagnostic/test_diagnostics.hpp b/tests/diagnostic/test_diagnostics.hpp index 01cf10859..02b24b052 100644 --- a/tests/diagnostic/test_diagnostics.hpp +++ b/tests/diagnostic/test_diagnostics.hpp @@ -226,7 +226,7 @@ void validateAttributes(Simulator& sim, Hi5Diagnostic& hi5) using GridLayout = typename Simulator::PHARETypes::GridLayout_t; constexpr auto dimension = Simulator::dimension; constexpr std::size_t expectedPopNbr = 2; - constexpr std::size_t expectedPopAttrFiles = 5; + constexpr std::size_t expectedPopAttrFiles = 6; std::string const ionsPopPath = "/ions/pop/"; @@ -236,7 +236,8 @@ void validateAttributes(Simulator& sim, Hi5Diagnostic& hi5) auto nbrPop = dict["simulation"]["ions"]["nbrPopulations"].template to(); EXPECT_EQ(nbrPop, expectedPopNbr); - std::vector h5FileTypes{"/EM_B", "/EM_E", "/ions/charge_density", "/ions/bulkVelocity"}; + std::vector h5FileTypes{"/EM_B", "/EM_E", "/ions/charge_density", + "/ions/mass_density", "/ions/bulkVelocity"}; for (std::size_t i = 0; i < nbrPop; ++i) { @@ -247,6 +248,7 @@ void validateAttributes(Simulator& sim, Hi5Diagnostic& hi5) h5FileTypes.emplace_back(ionsPopPath + popName + "/levelGhost"); h5FileTypes.emplace_back(ionsPopPath + popName + "/patchGhost"); h5FileTypes.emplace_back(ionsPopPath + popName + "/density"); + h5FileTypes.emplace_back(ionsPopPath + popName + "/charge_density"); h5FileTypes.emplace_back(ionsPopPath + popName + "/flux"); } diff --git a/tests/simulator/test_advance.py b/tests/simulator/test_advance.py index 40a673123..e7343f243 100644 --- a/tests/simulator/test_advance.py +++ b/tests/simulator/test_advance.py @@ -47,7 +47,7 @@ def getHierarchy( largest_patch_size=20, cells=120, time_step=0.001, - model_init={}, + model_init={"seed": 2}, dl=0.2, extra_diag_options={}, time_step_nbr=1, @@ -289,10 +289,28 @@ def base_test_overlaped_fields_are_equal(self, datahier, coarsest_time): # seems correct considering ghosts are filled with schedules # involving linear/spatial interpolations and so on where # rounding errors may occur.... setting atol to 5.5e-15 - assert_fp_any_all_close(slice1, slice2, atol=5.5e-15, rtol=0) + assert_fp_any_all_close(slice1, slice2, atol=2.5e-14, rtol=0) checks += 1 except AssertionError as e: print("AssertionError", pd1.name, e) + errors = np.where( + ~np.isclose(slice1, slice2, atol=2.5e-14, rtol=0) + ) + errors[0][:] += loc_b1.lower[0] + pd1.ghost_box.lower[0] + errors[1][:] += loc_b1.lower[1] + pd1.ghost_box.lower[1] + + fig = datahier.plot_2d_patches( + ilvl, + collections=[ + {"boxes": [pd1.box], "facecolor": 1}, + {"boxes": [pd2.box], "facecolor": 2}, + {"boxes": [errors], "facecolor": 3}, + ], + xlim=(-5, 50), + ylim=(-5, 50), + ) + fig.savefig(f"pd1.png") + print(f"ilvl {ilvl}") print(pd1.box, pd2.box) print(pd1.x.mean()) print(pd1.y.mean())