diff --git a/.gdbinit b/.gdbinit new file mode 100644 index 000000000..9c8c9d7a1 --- /dev/null +++ b/.gdbinit @@ -0,0 +1,3 @@ +break __assert_fail +break abort +catch throw diff --git a/.gitignore b/.gitignore index e31337989..3326515be 100644 --- a/.gitignore +++ b/.gitignore @@ -22,3 +22,4 @@ perf.* .vscode .phare* PHARE_REPORT.zip +.gdbinit diff --git a/pyphare/pyphare/pharein/diagnostics.py b/pyphare/pyphare/pharein/diagnostics.py index 2721c9633..3bd164232 100644 --- a/pyphare/pyphare/pharein/diagnostics.py +++ b/pyphare/pyphare/pharein/diagnostics.py @@ -222,6 +222,7 @@ def population_in_model(population): class FluidDiagnostics_(Diagnostics): fluid_quantities = [ "density", + "charge_density", "mass_density", "flux", "bulkVelocity", diff --git a/pyphare/pyphare/pharesee/hierarchy/__init__.py b/pyphare/pyphare/pharesee/hierarchy/__init__.py index 6bc7e43c8..318b1a24c 100644 --- a/pyphare/pyphare/pharesee/hierarchy/__init__.py +++ b/pyphare/pyphare/pharesee/hierarchy/__init__.py @@ -11,10 +11,11 @@ def hierarchy_from( - simulator=None, qty=None, pop="", h5_filename=None, times=None, hier=None, **kwargs + simulator=None, qty=None, pop="", h5_filename=None, times=None, hier=None, func=None, **kwargs ): from .fromh5 import hierarchy_fromh5 from .fromsim import hierarchy_from_sim + from .fromfunc import hierarchy_from_func """ this function reads an HDF5 PHARE file and returns a PatchHierarchy from @@ -39,4 +40,7 @@ def hierarchy_from( if simulator is not None and qty is not None: return hierarchy_from_sim(simulator, qty, pop=pop) + if func is not None and hier is not None: + return hierarchy_from_func(func, hier, **kwargs) + raise ValueError("can't make hierarchy") diff --git a/pyphare/pyphare/pharesee/hierarchy/fromfunc.py b/pyphare/pyphare/pharesee/hierarchy/fromfunc.py new file mode 100644 index 000000000..616852921 --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/fromfunc.py @@ -0,0 +1,97 @@ +from pyphare.pharesee.hierarchy.hierarchy_utils import compute_hier_from + +import numpy as np + + +def ions_mass_density_func1d(x, **kwargs): + masses = kwargs["masses"] # list of float : the ion pop masses + densities = kwargs["densities"] # list of callable : the ion pop density profiles + + assert len(masses) == len(densities) + funcs = np.zeros((x.size, len(masses))) + + for i, (mass, density) in enumerate(zip(masses, densities)): + funcs[:,i] = mass*density(x) + + return funcs.sum(axis=1) + + +def ions_charge_density_func1d(x, **kwargs): + charges = kwargs["charges"] # list of float : the ion pop charges + densities = kwargs["densities"] # list of callable : the ion pop density profiles + + assert len(charges) == len(densities) + + funcs = np.zeros((x.size, len(charges))) + + for i, (charge, density) in enumerate(zip(charges, densities)): + funcs[:,i] = charge*density(x) + + return funcs.sum(axis=1) + + +def hierarchy_from_func1d(func, hier, **kwargs): + assert hier.ndim == 1 + + def compute_(patch_datas, **kwargs): + ref_name = next(iter(patch_datas.keys())) + x_ = patch_datas[ref_name].x + + return ( + {"name": "value", "data": func(x_, **kwargs), "centering": patch_datas[ref_name].centerings}, + ) + + return compute_hier_from(compute_, hier, **kwargs) + + +def ions_mass_density_func2d(x, y, **kwargs): + masses = kwargs["masses"] # list of float : the ion pop masses + densities = kwargs["densities"] # list of callable : the ion pop density profiles + + yv, xv = np.meshgrid(y, x) + + assert len(masses) == len(densities) + funcs = np.zeros((x.size, y.size, len(masses))) + + for i, (mass, density) in enumerate(zip(masses, densities)): + funcs[:,:,i] = mass*density(xv, yv) + + return funcs.sum(axis=2) + + +def ions_charge_density_func2d(x, y, **kwargs): + charges = kwargs["charges"] # list of float : the ion pop charges + densities = kwargs["densities"] # list of callable : the ion pop density profiles + + yv, xv = np.meshgrid(y, x) + + assert len(charges) == len(densities) + funcs = np.zeros((x.size, y.size, len(charges))) + + for i, (charge, density) in enumerate(zip(charges, densities)): + funcs[:,:,i] = charge*density(xv, yv) + + return funcs.sum(axis=2) + + +def hierarchy_from_func2d(func, hier, **kwargs): + assert hier.ndim == 2 + + def compute_(patch_datas, **kwargs): + ref_name = next(iter(patch_datas.keys())) + x_ = patch_datas[ref_name].x + y_ = patch_datas[ref_name].y + + return ( + {"name": "value", "data": func(x_, y_, **kwargs), "centering": patch_datas[ref_name].centerings}, + ) + + return compute_hier_from(compute_, hier, **kwargs) + + +def hierarchy_from_func(func, hier, **kwargs): + if hier.ndim == 1: + return hierarchy_from_func1d(func, hier, **kwargs) + if hier.ndim == 2: + return hierarchy_from_func2d(func, hier, **kwargs) + diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py index a3ebc30a6..78308760d 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py @@ -360,6 +360,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] ): @@ -370,19 +372,36 @@ 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))) + + 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 - ) - ) + value = collection.get("value", np.nan) + for box in collection["boxes"]: + i0, j0 = box.lower + i1, j1 = box.upper + ij[i0 : i1 + 1, j0 : j1 + 1] = value + if "coords" in collection: + for coords in collection["coords"]: + ij[coords] = collection["value"] + + 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 @@ -390,15 +409,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 @@ -433,7 +447,7 @@ def plot1d(self, **kwargs): qty = pdata_names[0] layout = patch.patch_datas[qty].layout - nbrGhosts = layout.nbrGhostFor(qty) + nbrGhosts = patch.patch_datas[qty].ghosts_nbr val = patch.patch_datas[qty][patch.box] x = patch.patch_datas[qty].x[nbrGhosts[0] : -nbrGhosts[0]] label = "L{level}P{patch}".format(level=lvl_nbr, patch=ip) @@ -544,9 +558,10 @@ def plot2d(self, **kwargs): if "ylim" in kwargs: ax.set_ylim(kwargs["ylim"]) - divider = make_axes_locatable(ax) - cax = divider.append_axes("right", size="5%", pad=0.08) - plt.colorbar(im, ax=ax, cax=cax) + if kwargs.get("cbar", True): + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="5%", pad=0.08) + fig.colorbar(im, ax=ax, cax=cax) if kwargs.get("legend", None) is not None: ax.legend() diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py index 487dc8605..b7fbcc220 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py @@ -36,6 +36,7 @@ "momentum_tensor_zz": "Mzz", "density": "rho", "mass_density": "rho", + "charge_density": "rho", "tags": "tags", } diff --git a/pyphare/pyphare/pharesee/plotting.py b/pyphare/pyphare/pharesee/plotting.py index 180db4375..8f7614919 100644 --- a/pyphare/pyphare/pharesee/plotting.py +++ b/pyphare/pyphare/pharesee/plotting.py @@ -274,7 +274,7 @@ def finest_field_plot(run_path, qty, **kwargs): time = times[0] interpolator, finest_coords = r.GetVi(time, merged=True, interp=interp)[qty] elif qty == "rho": - file = os.path.join(run_path, "ions_density.h5") + file = os.path.join(run_path, "ions_charge_density.h5") if time is None: times = get_times_from_h5(file) time = times[0] diff --git a/pyphare/pyphare/pharesee/run/run.py b/pyphare/pyphare/pharesee/run/run.py index 003943497..30f6a25f8 100644 --- a/pyphare/pyphare/pharesee/run/run.py +++ b/pyphare/pyphare/pharesee/run/run.py @@ -27,7 +27,7 @@ "EM_B": "B", "EM_E": "E", "ions_bulkVelocity": "Vi", - "ions_density": "Ni", + "ions_charge_density": "Ni", "particle_count": "nppc", } @@ -114,7 +114,7 @@ def GetMassDensity(self, time, merged=False, interp="nearest", **kwargs): return ScalarField(self._get(hier, time, merged, interp)) def GetNi(self, time, merged=False, interp="nearest", **kwargs): - hier = self._get_hierarchy(time, "ions_density.h5", **kwargs) + hier = self._get_hierarchy(time, "ions_charge_density.h5", **kwargs) return ScalarField(self._get(hier, time, merged, interp)) def GetN(self, time, pop_name, merged=False, interp="nearest", **kwargs): @@ -151,7 +151,7 @@ def GetPi(self, time, merged=False, interp="nearest", **kwargs): return self._get(Pi, time, merged, interp) # should later be a TensorField def GetPe(self, time, merged=False, interp="nearest", all_primal=True): - hier = self._get_hierarchy(time, "ions_density.h5") + hier = self._get_hierarchy(time, "ions_charge_density.h5") Te = hier.sim.electrons.closure.Te diff --git a/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py b/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py index 29a8249e3..18bfc905d 100644 --- a/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py +++ b/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py @@ -177,11 +177,11 @@ def test_large_patchoverlaps(self, expected): collections=[ { "boxes": overlap_boxes, - "facecolor": "yellow", + "value": 1, }, { "boxes": [p.box for p in hierarchy.level(ilvl).patches], - "facecolor": "grey", + "value": 2, }, ], ) @@ -390,11 +390,11 @@ def test_level_ghost_boxes(self, refinement_boxes, expected): collections=[ { "boxes": ghost_area_box_list, - "facecolor": "yellow", + "value": 1, }, { "boxes": [p.box for p in hierarchy.level(ilvl).patches], - "facecolor": "grey", + "value": 2, }, ], title="".join( @@ -502,7 +502,7 @@ def test_patch_periodicity_copy(self, refinement_boxes, expected): collections=[ { "boxes": [p.box for p in periodic_list], - "facecolor": "grey", + "value": 2, }, ], ) diff --git a/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py b/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py index f65bb39f0..908484607 100644 --- a/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py +++ b/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py @@ -145,7 +145,7 @@ def vthz(x, y): for quantity in ["E", "B"]: ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) - for quantity in ["density", "bulkVelocity"]: + for quantity in ["charge_density", "bulkVelocity"]: ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps) return sim diff --git a/src/amr/level_initializer/hybrid_level_initializer.hpp b/src/amr/level_initializer/hybrid_level_initializer.hpp index e0ef8386a..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.computeDensity(); + ions.computeChargeDensity(); ions.computeBulkVelocity(); } diff --git a/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp b/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp index 6c708b40d..84acb79d3 100644 --- a/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp +++ b/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp @@ -385,20 +385,21 @@ namespace amr { // first thing to do is to project patchGhostParitcles moments auto& patchGhosts = pop.patchGhostParticles(); - auto& density = pop.density(); + auto& particleDensity = pop.particleDensity(); + auto& chargeDensity = pop.chargeDensity(); auto& flux = pop.flux(); - interpolate_(makeRange(patchGhosts), density, flux, layout); + interpolate_(makeRange(patchGhosts), particleDensity, chargeDensity, flux, layout); if (level.getLevelNumber() > 0) // no levelGhost on root level { // then grab levelGhostParticlesOld and levelGhostParticlesNew // and project them with alpha and (1-alpha) coefs, respectively auto& levelGhostOld = pop.levelGhostParticlesOld(); - interpolate_(makeRange(levelGhostOld), density, flux, layout, 1. - alpha); + interpolate_(makeRange(levelGhostOld), particleDensity, chargeDensity, flux, layout, 1. - alpha); auto& levelGhostNew = pop.levelGhostParticlesNew(); - interpolate_(makeRange(levelGhostNew), density, flux, layout, alpha); + interpolate_(makeRange(levelGhostNew), particleDensity, chargeDensity, flux, layout, alpha); } } } @@ -519,7 +520,7 @@ namespace amr auto& J = hybridModel.state.J; auto& Vi = hybridModel.state.ions.velocity(); - auto& Ni = hybridModel.state.ions.density(); + auto& Ni = hybridModel.state.ions.chargeDensity(); Jold_.copyData(J); ViOld_.copyData(Vi); diff --git a/src/amr/physical_models/hybrid_model.hpp b/src/amr/physical_models/hybrid_model.hpp index 449bafc22..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.densityName(); + hybridInfo.modelIonDensity = state.ions.chargeDensityName(); hybridInfo.modelIonBulkVelocity = core::VecFieldNames{state.ions.velocity()}; hybridInfo.modelCurrent = core::VecFieldNames{state.J}; diff --git a/src/amr/tagging/default_hybrid_tagger_strategy.hpp b/src/amr/tagging/default_hybrid_tagger_strategy.hpp index b9a016f7d..741a52b4e 100644 --- a/src/amr/tagging/default_hybrid_tagger_strategy.hpp +++ b/src/amr/tagging/default_hybrid_tagger_strategy.hpp @@ -36,7 +36,7 @@ void DefaultHybridTaggerStrategy::tag(HybridModel& model, auto& By = model.state.electromag.B.getComponent(PHARE::core::Component::Y); auto& Bz = model.state.electromag.B.getComponent(PHARE::core::Component::Z); - auto& N = model.state.ions.density(); + auto& N = model.state.ions.chargeDensity(); // we loop on cell indexes for all qties regardless of their centering auto const& [start_x, _] diff --git a/src/core/data/electrons/electrons.hpp b/src/core/data/electrons/electrons.hpp index a079022da..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 { @@ -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_; }); } diff --git a/src/core/data/ions/ion_population/ion_population.hpp b/src/core/data/ions/ion_population/ion_population.hpp index aa6554e47..33dfea9d8 100644 --- a/src/core/data/ions/ion_population/ion_population.hpp +++ b/src/core/data/ions/ion_population/ion_population.hpp @@ -31,7 +31,8 @@ namespace core , mass_{initializer["mass"].template to()} , flux_{name_ + "_flux", HybridQuantity::Vector::V} , momentumTensor_{name_ + "_momentumTensor", HybridQuantity::Tensor::M} - , rho_{name_ + "_rho", HybridQuantity::Scalar::rho} + , particleDensity_{name_ + "_particleDensity", HybridQuantity::Scalar::rho} + , chargeDensity_{name_ + "_chargeDensity", HybridQuantity::Scalar::rho} , particles_{name_} , particleInitializerInfo_{initializer["particle_initializer"]} { @@ -42,20 +43,16 @@ namespace core NO_DISCARD std::string const& name() const { return name_; } - NO_DISCARD auto const& particleInitializerInfo() const { return particleInitializerInfo_; } - - NO_DISCARD bool isUsable() const { - return core::isUsable(particles_, rho_, flux_, momentumTensor_); + return core::isUsable(particles_, particleDensity_, chargeDensity_, flux_, momentumTensor_); } - NO_DISCARD bool isSettable() const { - return core::isSettable(particles_, rho_, flux_, momentumTensor_); + return core::isSettable(particles_, particleDensity_, chargeDensity_, flux_, momentumTensor_); } NO_DISCARD auto& domainParticles() const { return particles_.domainParticles(); } @@ -79,9 +76,11 @@ namespace core return particles_.levelGhostParticlesNew(); } + NO_DISCARD field_type const& particleDensity() const { return particleDensity_; } + NO_DISCARD field_type& particleDensity() { return particleDensity_; } - NO_DISCARD field_type const& density() const { return rho_; } - NO_DISCARD field_type& density() { return rho_; } + NO_DISCARD field_type const& chargeDensity() const { return chargeDensity_; } + NO_DISCARD field_type& chargeDensity() { return chargeDensity_; } NO_DISCARD VecField const& flux() const { return flux_; } NO_DISCARD VecField& flux() { return flux_; } @@ -100,7 +99,8 @@ namespace core NO_DISCARD auto getCompileTimeResourcesViewList() { - return std::forward_as_tuple(flux_, momentumTensor_, rho_, particles_); + return std::forward_as_tuple(flux_, momentumTensor_, particleDensity_, chargeDensity_, + particles_); } @@ -124,7 +124,8 @@ namespace core double mass_; VecField flux_; TensorField momentumTensor_; - field_type rho_; + field_type particleDensity_; + field_type chargeDensity_; ParticlesPack particles_; initializer::PHAREDict const& particleInitializerInfo_; }; diff --git a/src/core/data/ions/ions.hpp b/src/core/data/ions/ions.hpp index 2701ebdac..a9a536802 100644 --- a/src/core/data/ions/ions.hpp +++ b/src/core/data/ions/ions.hpp @@ -41,15 +41,14 @@ namespace core explicit Ions(PHARE::initializer::PHAREDict const& dict) - : rho_{densityName(), 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( [&dict](auto ipop) { // return IonPopulation{dict["pop" + std::to_string(ipop)]}; }, dict["nbrPopulations"].template to())} - , sameMasses_{allSameMass_()} , momentumTensor_{"momentumTensor", HybridQuantity::Tensor::M} { } @@ -58,29 +57,24 @@ namespace core NO_DISCARD auto nbrPopulations() const { return populations_.size(); } NO_DISCARD auto size() const { return nbrPopulations(); } + NO_DISCARD field_type const& massDensity() const { return massDensity_; } + NO_DISCARD field_type const& massDensity() { return massDensity_; } - NO_DISCARD field_type const& density() const { return rho_; } - NO_DISCARD field_type& density() { return rho_; } - - NO_DISCARD field_type const& massDensity() const - { - return sameMasses_ ? rho_ : massDensity_; - } - NO_DISCARD field_type& massDensity() { return sameMasses_ ? rho_ : massDensity_; } - + NO_DISCARD field_type const& chargeDensity() const { return chargeDensity_; } + NO_DISCARD field_type& chargeDensity() { return chargeDensity_; } NO_DISCARD vecfield_type const& velocity() const { return bulkVelocity_; } NO_DISCARD vecfield_type& velocity() { return bulkVelocity_; } - NO_DISCARD std::string static densityName() { return "rho"; } + 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 computeDensity() + void computeChargeDensity() { - rho_.zero(); + chargeDensity_.zero(); for (auto const& pop : populations_) { @@ -88,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.density(); - std::transform(std::begin(rho_), std::end(rho_), std::begin(popDensity), - std::begin(rho_), 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(); @@ -103,7 +99,7 @@ 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.density(); + auto& popDensity = pop.particleDensity(); std::transform( std::begin(massDensity_), std::end(massDensity_), std::begin(popDensity), std::begin(massDensity_), @@ -114,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 rho_ 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_) ? rho_ : massDensity_; + computeMassDensity(); bulkVelocity_.zero(); auto& vx = bulkVelocity_.getComponent(Component::X); @@ -131,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{}); } @@ -185,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 - = rho_.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_) { @@ -202,11 +186,8 @@ namespace core NO_DISCARD bool isSettable() const { - bool settable - = rho_.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_) { @@ -230,7 +211,8 @@ namespace core NO_DISCARD auto getCompileTimeResourcesViewList() { - return std::forward_as_tuple(bulkVelocity_, momentumTensor_, rho_, massDensity_); + return std::forward_as_tuple(bulkVelocity_, momentumTensor_, chargeDensity_, + massDensity_); } @@ -251,29 +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 rho_; 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 a0a6a5a96..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& density, 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,8 +488,11 @@ class Interpolator : private Weighter currPart->delta); particleToMesh_( - density, *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); particleToMesh_( xFlux, *currPart, [](auto const& part) { return part.v[0]; }, startIndex_, weights_, coef); @@ -502,10 +506,10 @@ class Interpolator : private Weighter PHARE_LOG_STOP(3, "ParticleToMesh::operator()"); } template - inline void operator()(ParticleRange&& range, Field& density, 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, density, flux, layout, coef); + (*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 dbae99b75..6f51a6147 100644 --- a/src/core/numerics/ion_updater/ion_updater.hpp +++ b/src/core/numerics/ion_updater/ion_updater.hpp @@ -107,7 +107,7 @@ void IonUpdater::updatePopulations(Ions& ions, Ele template void IonUpdater::updateIons(Ions& ions) { - ions.computeDensity(); + ions.computeChargeDensity(); ions.computeBulkVelocity(); } @@ -158,7 +158,7 @@ void IonUpdater::updateAndDepositDomain_(Ions& ion inRange, outRange, em, pop.mass(), interpolator_, layout, [](auto& particleRange) { return particleRange; }, inDomainBox); - interpolator_(inDomain, pop.density(), pop.flux(), layout); + interpolator_(inDomain, pop.particleDensity(), pop.chargeDensity(), pop.flux(), layout); // TODO : we can erase here because we know we are working on a state // that has been saved in the solverPPC @@ -181,7 +181,8 @@ void IonUpdater::updateAndDepositDomain_(Ions& ion auto enteredInDomain = pusher_->move(inRange, outRange, em, pop.mass(), interpolator_, layout, inGhostBox, inDomainBox); - interpolator_(enteredInDomain, pop.density(), pop.flux(), layout); + interpolator_(enteredInDomain, pop.particleDensity(), pop.chargeDensity(), pop.flux(), + layout); if (copyInDomain) { @@ -269,7 +270,8 @@ void IonUpdater::updateAndDepositAll_(Ions& ions, pushAndCopyInDomain(makeIndexRange(pop.patchGhostParticles())); pushAndCopyInDomain(makeIndexRange(pop.levelGhostParticles())); - interpolator_(makeIndexRange(domainParticles), pop.density(), 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 d71135614..87b805fa1 100644 --- a/src/core/numerics/moments/moments.hpp +++ b/src/core/numerics/moments/moments.hpp @@ -15,7 +15,8 @@ namespace core { for (auto& pop : ions) { - pop.density().zero(); + pop.particleDensity().zero(); + pop.chargeDensity().zero(); pop.flux().zero(); } } @@ -40,23 +41,24 @@ namespace core { for (auto& pop : ions) { - auto& density = pop.density(); - auto& flux = pop.flux(); + auto& particleDensity = pop.particleDensity(); + auto& chargeDensity = pop.chargeDensity(); + auto& flux = pop.flux(); if constexpr (std::is_same_v) { auto& partArray = pop.domainParticles(); - interpolate(partArray, density, flux, layout); + interpolate(partArray, particleDensity, chargeDensity, flux, layout); } else if constexpr (std::is_same_v) { auto& partArray = pop.patchGhostParticles(); - interpolate(partArray, density, flux, layout); + interpolate(partArray, particleDensity, chargeDensity, flux, layout); } else if constexpr (std::is_same_v) { auto& partArray = pop.levelGhostParticlesOld(); - interpolate(partArray, density, flux, layout); + interpolate(partArray, particleDensity, chargeDensity, flux, layout); } } } diff --git a/src/diagnostic/detail/types/fluid.hpp b/src/diagnostic/detail/types/fluid.hpp index ce8d98c78..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, "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.density(), "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")) @@ -206,8 +209,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, "charge_density")) + infoDS(ions.chargeDensity(), "charge_density", patchAttributes[lvlPatchID]["ion"]); if (isActiveDiag(diagnostic, tree, "mass_density")) infoDS(ions.massDensity(), "mass_density", patchAttributes[lvlPatchID]["ion"]); if (isActiveDiag(diagnostic, tree, "bulkVelocity")) @@ -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")) @@ -275,8 +280,8 @@ void FluidDiagnosticWriter::initDataSets( } std::string tree{"/ions/"}; - if (isActiveDiag(diagnostic, tree, "density")) - initDS(path, attr["ion"], "density", null); + if (isActiveDiag(diagnostic, tree, "charge_density")) + initDS(path, attr["ion"], "charge_density", null); if (isActiveDiag(diagnostic, tree, "mass_density")) initDS(path, attr["ion"], "mass_density", null); if (isActiveDiag(diagnostic, tree, "bulkVelocity")) @@ -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.density()); + 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")) @@ -315,8 +322,8 @@ void FluidDiagnosticWriter::write(DiagnosticProperties& diagnostic) } std::string tree{"/ions/"}; - if (isActiveDiag(diagnostic, tree, "density")) - writeDS(path + "density", ions.density()); + if (isActiveDiag(diagnostic, tree, "charge_density")) + writeDS(path + "charge_density", ions.chargeDensity()); if (isActiveDiag(diagnostic, tree, "mass_density")) writeDS(path + "mass_density", ions.massDensity()); if (isActiveDiag(diagnostic, tree, "bulkVelocity")) @@ -345,6 +352,7 @@ void FluidDiagnosticWriter::writeAttributes( { 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/src/phare/phare_init.py b/src/phare/phare_init.py index 2c90bce3a..d79de42b7 100644 --- a/src/phare/phare_init.py +++ b/src/phare/phare_init.py @@ -101,7 +101,7 @@ def vthz(x): ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) -for quantity in ["density", "bulkVelocity"]: +for quantity in ["charge_density", "bulkVelocity"]: ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps) pops = [ diff --git a/src/phare/phare_init_small.py b/src/phare/phare_init_small.py index b4bccc1cb..f3f8e76e1 100644 --- a/src/phare/phare_init_small.py +++ b/src/phare/phare_init_small.py @@ -104,7 +104,7 @@ def vthz(x): ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) -for quantity in ["density", "bulkVelocity"]: +for quantity in ["charge_density", "bulkVelocity"]: ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps) pops = [ diff --git a/src/python3/patch_level.hpp b/src/python3/patch_level.hpp index 517c03baa..b7600a44a 100644 --- a/src/python3/patch_level.hpp +++ b/src/python3/patch_level.hpp @@ -37,7 +37,7 @@ class PatchLevel auto& ions = model_.state.ions; auto visit = [&](GridLayout& grid, std::string patchID, std::size_t /*iLevel*/) { - setPatchDataFromField(patchDatas.emplace_back(), ions.density(), grid, patchID); + setPatchDataFromField(patchDatas.emplace_back(), ions.chargeDensity(), grid, patchID); }; PHARE::amr::visitLevel(*hierarchy_.getPatchLevel(lvl_), @@ -59,7 +59,7 @@ class PatchLevel if (!pop_data.count(pop.name())) pop_data.emplace(pop.name(), Inner()); - setPatchDataFromField(pop_data.at(pop.name()).emplace_back(), pop.density(), grid, + setPatchDataFromField(pop_data.at(pop.name()).emplace_back(), pop.chargeDensity(), grid, patchID); } }; diff --git a/tests/core/data/electrons/test_electrons.cpp b/tests/core/data/electrons/test_electrons.cpp index d69df65ff..5f5ffc700 100644 --- a/tests/core/data/electrons/test_electrons.cpp +++ b/tests/core/data/electrons/test_electrons.cpp @@ -126,9 +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 densityName = std::string_view{"rho"}; + static constexpr auto dim = typename TypeInfo::first_type{}(); + static constexpr auto interp = typename TypeInfo::second_type{}(); using GridYee = GridLayout>; @@ -149,9 +148,9 @@ struct ElectronsTest : public ::testing::Test Electromag electromag; UsableVecField J, F, Ve, Vi; - UsableTensorField M, protons_M; + UsableTensorField ionTensor, protonTensor; - GridND Nibuffer, NiProtons, Pe; + GridND ionChargeDensity, ionMassDensity, protonParticleDensity, protonChargeDensity, Pe; ParticleArray_t domainParticles{layout.AMRBox()}; ParticleArray_t patchGhostParticles = domainParticles; @@ -164,22 +163,25 @@ struct ElectronsTest : public ::testing::Test template auto static _ions(Args&... args) { - auto const& [Fi, Nibuffer, NiProtons, Vi, M, protons_M, 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, md] = ions.getCompileTimeResourcesViewList(); - d.setBuffer(&Nibuffer); + auto const& [V, m, d_c, d_m] = ions.getCompileTimeResourcesViewList(); + d_c.setBuffer(&ionChargeDensity); + d_m.setBuffer(&ionMassDensity); Vi.set_on(V); - M.set_on(m); + ionTensor.set_on(m); } auto& pops = ions.getRunTimeResourcesViewList(); assert(pops.size() == 1); - auto const& [F, m, d, poppack] = pops[0].getCompileTimeResourcesViewList(); - d.setBuffer(&NiProtons); - Fi.set_on(F); - protons_M.set_on(m); + auto const& [F, m, Np, Nc, poppack] = pops[0].getCompileTimeResourcesViewList(); + Np.setBuffer(&protonParticleDensity); + Nc.setBuffer(&protonChargeDensity); + ionFlux.set_on(F); + protonTensor.set_on(m); poppack.setBuffer(&pack); return ions; } @@ -191,14 +193,19 @@ struct ElectronsTest : public ::testing::Test , F{"protons_flux", layout, HybridQuantity::Vector::V} , Ve{"StandardHybridElectronFluxComputer_Ve", layout, HybridQuantity::Vector::V} , Vi{"bulkVel", layout, HybridQuantity::Vector::V} - , M{"momentumTensor", layout, HybridQuantity::Tensor::M} - , protons_M{"protons_momentumTensor", layout, HybridQuantity::Tensor::M} - , Nibuffer{std::string{densityName}, HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} - , NiProtons{"protons_rho", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} + , ionTensor{"momentumTensor", layout, HybridQuantity::Tensor::M} + , protonTensor{"protons_momentumTensor", layout, HybridQuantity::Tensor::M} + , 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)} + , protonChargeDensity{"protons_chargeDensity", HybridQuantity::Scalar::rho, + layout.allocSize(HybridQuantity::Scalar::rho)} , Pe{"Pe", HybridQuantity::Scalar::P, layout.allocSize(HybridQuantity::Scalar::P)} - , ions{_ions(F, Nibuffer, NiProtons, Vi, M, protons_M, pack)} + , ions{_ions(F, ionChargeDensity, ionMassDensity, protonParticleDensity, + protonChargeDensity, Vi, ionTensor, protonTensor, pack)} , electrons{createDict()["electrons"], ions, J} { auto&& emm = std::get<0>(electrons.getCompileTimeResourcesViewList()); @@ -237,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(Nibuffer, [](double x) { return std::cosh(0.1 * x); }); + fill(ionChargeDensity, [](double x) { return std::cosh(0.1 * x); }); } else if constexpr (dim == 2) { @@ -266,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(Nibuffer, + fill(ionChargeDensity, [](double x, double y) { return std::cosh(0.1 * x) * std::cosh(0.1 * y); }); } else if constexpr (dim == 3) @@ -313,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(Nibuffer, [](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); }); } @@ -363,7 +370,7 @@ TYPED_TEST(ElectronsTest, ThatElectronsDensityEqualIonDensity) electrons.update(layout); auto& Ne = electrons.density(); - auto& Ni = ions.density(); + auto& Ni = ions.chargeDensity(); if constexpr (dim == 1) { diff --git a/tests/core/data/ion_population/test_ion_population_fixtures.hpp b/tests/core/data/ion_population/test_ion_population_fixtures.hpp index e26ee34e3..1fe25487d 100644 --- a/tests/core/data/ion_population/test_ion_population_fixtures.hpp +++ b/tests/core/data/ion_population/test_ion_population_fixtures.hpp @@ -50,7 +50,7 @@ class UsableIonsPopulation : public _defaults::IonPopulation_t public: UsableIonsPopulation(initializer::PHAREDict const& dict, GridLayout_t const& layout) : Super{dict} - , rho{this->name() + "_rho", layout, HybridQuantity::Scalar::rho} + , rho{this->name() + "_particleDensity", layout, HybridQuantity::Scalar::rho} , F{this->name() + "_flux", layout, HybridQuantity::Vector::V} , M{this->name() + "_momentumTensor", layout, HybridQuantity::Tensor::M} , particles{this->name(), layout.AMRBox()} diff --git a/tests/core/data/ions/test_ions.cpp b/tests/core/data/ions/test_ions.cpp index 7fb1b1426..11769fac2 100644 --- a/tests/core/data/ions/test_ions.cpp +++ b/tests/core/data/ions/test_ions.cpp @@ -154,7 +154,7 @@ TEST_F(theIons, areSettableUponConstruction) #ifndef NDEBUG // no throw in release mode! JUST SEGFAULTS! :D TEST_F(theIons, throwIfAccessingDensityWhileNotUsable) { - EXPECT_ANY_THROW(auto& n = ions.density()(0)); + EXPECT_ANY_THROW(auto& n = ions.chargeDensity()(0)); } #endif diff --git a/tests/core/numerics/interpolator/test_main.cpp b/tests/core/numerics/interpolator/test_main.cpp index 28e2fff14..c930ae72f 100644 --- a/tests/core/numerics/interpolator/test_main.cpp +++ b/tests/core/numerics/interpolator/test_main.cpp @@ -516,16 +516,23 @@ class ACollectionOfParticles_1d : public ::testing::Test ParticleArray_t particles; Grid_t rho; + Grid_t rho_c; UsableVecFieldND v; std::array weights; + // For the set of particles herebelow, whatever the interp_order, + // we have a set of particles around idx=20, so that their particle + // density is 1, their charge density is 2 and velocity is (2, -1, 1) + // depending on the interp order, iCell, delta, weight, charge and v + // are then correctly tuned ACollectionOfParticles_1d() : part{} , particles{grow(layout.AMRBox(), safeLayer)} , rho{"field", HybridQuantity::Scalar::rho, nx} + , rho_c{"field", HybridQuantity::Scalar::rho, nx} , v{"v", layout, HybridQuantity::Vector::V} { if constexpr (Interpolator::interp_order == 1) @@ -533,6 +540,7 @@ class ACollectionOfParticles_1d : public ::testing::Test part.iCell[0] = 19; // AMR index part.delta[0] = 0.5; part.weight = 1.0; + part.charge = 2.0; part.v[0] = +2.; part.v[1] = -1.; part.v[2] = +1.; @@ -541,6 +549,7 @@ class ACollectionOfParticles_1d : public ::testing::Test part.iCell[0] = 20; // AMR index part.delta[0] = 0.5; part.weight = 0.4; + part.charge = 1.85; part.v[0] = +2.; part.v[1] = -1.; part.v[2] = +1.; @@ -549,6 +558,7 @@ class ACollectionOfParticles_1d : public ::testing::Test part.iCell[0] = 20; // AMR index part.delta[0] = 0.5; part.weight = 0.6; + part.charge = 2.1; part.v[0] = +2.; part.v[1] = -1.; part.v[2] = +1.; @@ -560,6 +570,7 @@ class ACollectionOfParticles_1d : public ::testing::Test part.iCell[0] = 19; // AMR index part.delta[0] = 0.0; part.weight = 1.0; + part.charge = 2.0; part.v[0] = +2.; part.v[1] = -1.; part.v[2] = +1.; @@ -568,6 +579,7 @@ class ACollectionOfParticles_1d : public ::testing::Test part.iCell[0] = 20; // AMR index part.delta[0] = 0.0; part.weight = 0.2; + part.charge = 3.2; part.v[0] = +2.; part.v[1] = -1.; part.v[2] = +1.; @@ -576,6 +588,7 @@ class ACollectionOfParticles_1d : public ::testing::Test part.iCell[0] = 20; // AMR index part.delta[0] = 0.0; part.weight = 0.8; + part.charge = 1.7; part.v[0] = +2.; part.v[1] = -1.; part.v[2] = +1.; @@ -584,6 +597,7 @@ class ACollectionOfParticles_1d : public ::testing::Test part.iCell[0] = 21; // AMR index part.delta[0] = 0.0; part.weight = 1.0; + part.charge = 2.0; part.v[0] = +2.; part.v[1] = -1.; part.v[2] = +1.; @@ -595,6 +609,7 @@ class ACollectionOfParticles_1d : public ::testing::Test part.iCell[0] = 18; // AMR index part.delta[0] = 0.5; part.weight = 1.0; + part.charge = 2.0; part.v[0] = +2.; part.v[1] = -1.; part.v[2] = +1.; @@ -603,6 +618,7 @@ class ACollectionOfParticles_1d : public ::testing::Test part.iCell[0] = 19; // AMR index part.delta[0] = 0.5; part.weight = 1.0; + part.charge = 2.0; part.v[0] = +2.; part.v[1] = -1.; part.v[2] = +1.; @@ -611,6 +627,7 @@ class ACollectionOfParticles_1d : public ::testing::Test part.iCell[0] = 20; // AMR index part.delta[0] = 0.5; part.weight = 1.0; + part.charge = 2.0; part.v[0] = +2.; part.v[1] = -1.; part.v[2] = +1.; @@ -619,6 +636,7 @@ class ACollectionOfParticles_1d : public ::testing::Test part.iCell[0] = 21; // AMR index part.delta[0] = 0.5; part.weight = 0.1; + part.charge = 3.35; part.v[0] = +2.; part.v[1] = -1.; part.v[2] = +1.; @@ -627,12 +645,13 @@ class ACollectionOfParticles_1d : public ::testing::Test part.iCell[0] = 21; // AMR index part.delta[0] = 0.5; part.weight = 0.9; + part.charge = 1.85; part.v[0] = +2.; part.v[1] = -1.; part.v[2] = +1.; particles.push_back(part); } - interpolator(makeIndexRange(particles), rho, v, layout); + interpolator(makeIndexRange(particles), rho, rho_c, v, layout); } @@ -650,6 +669,7 @@ TYPED_TEST_P(ACollectionOfParticles_1d, DepositCorrectlyTheirWeight_1d) auto const& [vx, vy, vz] = this->v(); EXPECT_DOUBLE_EQ(this->rho(idx), 1.0); + EXPECT_DOUBLE_EQ(this->rho_c(idx), 2.0); EXPECT_DOUBLE_EQ(vx(idx), 2.0); EXPECT_DOUBLE_EQ(vy(idx), -1.0); EXPECT_DOUBLE_EQ(vz(idx), 1.0); @@ -678,12 +698,14 @@ struct ACollectionOfParticles_2d : public ::testing::Test GridLayout_t layout{ConstArray(.1), {nx, ny}, ConstArray(0)}; ParticleArray_t particles; Grid_t rho; + Grid_t rho_c; UsableVecFieldND v; Interpolator interpolator; ACollectionOfParticles_2d() : particles{grow(layout.AMRBox(), safeLayer)} , rho{"field", HybridQuantity::Scalar::rho, nx, ny} + , rho_c{"field", HybridQuantity::Scalar::rho, nx, ny} , v{"v", layout, HybridQuantity::Vector::V} { for (int i = start; i < end; i++) @@ -697,7 +719,7 @@ struct ACollectionOfParticles_2d : public ::testing::Test part.v[1] = -1.; part.v[2] = +1.; } - interpolator(makeIndexRange(particles), rho, v, layout); + interpolator(makeIndexRange(particles), rho, rho_c, v, layout); } }; TYPED_TEST_SUITE_P(ACollectionOfParticles_2d); diff --git a/tests/core/numerics/ion_updater/test_updater.cpp b/tests/core/numerics/ion_updater/test_updater.cpp index 2d939cf7d..09013e517 100644 --- a/tests/core/numerics/ion_updater/test_updater.cpp +++ b/tests/core/numerics/ion_updater/test_updater.cpp @@ -211,10 +211,12 @@ struct IonsBuffers using ParticleArray = typename PHARETypes::ParticleArray_t; using ParticleInitializerFactory = typename PHARETypes::ParticleInitializerFactory; - Grid ionDensity; + Grid ionChargeDensity; Grid ionMassDensity; - Grid protonDensity; - Grid alphaDensity; + Grid protonParticleDensity; + Grid protonChargeDensity; + Grid alphaParticleDensity; + Grid alphaChargeDensity; UsableVecFieldND protonF, alphaF, Vi; UsableTensorField M, alpha_M, protons_M; @@ -237,14 +239,18 @@ struct IonsBuffers ParticlesPack alphaPack; IonsBuffers(GridLayout const& layout) - : ionDensity{"rho", 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)} - , protonDensity{"protons_rho", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} - , alphaDensity{"alpha_rho", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} + , protonParticleDensity{"protons_particleDensity", HybridQuantity::Scalar::rho, + layout.allocSize(HybridQuantity::Scalar::rho)} + , protonChargeDensity{"protons_chargeDensity", HybridQuantity::Scalar::rho, + layout.allocSize(HybridQuantity::Scalar::rho)} + , alphaParticleDensity{"alpha_particleDensity", HybridQuantity::Scalar::rho, + layout.allocSize(HybridQuantity::Scalar::rho)} + , alphaChargeDensity{"alpha_chargeDensity", 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} @@ -270,14 +276,18 @@ struct IonsBuffers IonsBuffers(IonsBuffers const& source, GridLayout const& layout) - : ionDensity{"rho", 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)} - , protonDensity{"protons_rho", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} - , alphaDensity{"alpha_rho", HybridQuantity::Scalar::rho, - layout.allocSize(HybridQuantity::Scalar::rho)} + , protonParticleDensity{"protons_particleDensity", HybridQuantity::Scalar::rho, + layout.allocSize(HybridQuantity::Scalar::rho)} + , protonChargeDensity{"protons_chargeDensity", HybridQuantity::Scalar::rho, + layout.allocSize(HybridQuantity::Scalar::rho)} + , alphaParticleDensity{"alpha_particleDensity", HybridQuantity::Scalar::rho, + layout.allocSize(HybridQuantity::Scalar::rho)} + , alphaChargeDensity{"alpha_chargeDensity", 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} @@ -300,10 +310,12 @@ struct IonsBuffers &alphaLevelGhost, &alphaLevelGhostOld, &alphaLevelGhostNew} { - ionDensity.copyData(source.ionDensity); + ionChargeDensity.copyData(source.ionChargeDensity); ionMassDensity.copyData(source.ionMassDensity); - protonDensity.copyData(source.protonDensity); - alphaDensity.copyData(source.alphaDensity); + protonParticleDensity.copyData(source.protonParticleDensity); + protonChargeDensity.copyData(source.protonChargeDensity); + alphaParticleDensity.copyData(source.alphaParticleDensity); + alphaChargeDensity.copyData(source.alphaChargeDensity); protonF.copyData(source.protonF); alphaF.copyData(source.alphaF); @@ -313,25 +325,27 @@ struct IonsBuffers void setBuffers(Ions& ions) { { - auto const& [V, m, d, md] = ions.getCompileTimeResourcesViewList(); + auto const& [V, m, cd, md] = ions.getCompileTimeResourcesViewList(); Vi.set_on(V); M.set_on(m); - d.setBuffer(&ionDensity); + cd.setBuffer(&ionChargeDensity); md.setBuffer(&ionMassDensity); } auto& pops = ions.getRunTimeResourcesViewList(); { - auto const& [F, M, d, particles] = pops[0].getCompileTimeResourcesViewList(); - d.setBuffer(&protonDensity); + auto const& [F, M, d, c, particles] = pops[0].getCompileTimeResourcesViewList(); + d.setBuffer(&protonParticleDensity); + c.setBuffer(&protonChargeDensity); protons_M.set_on(M); protonF.set_on(F); particles.setBuffer(&protonPack); } { - auto const& [F, M, d, particles] = pops[1].getCompileTimeResourcesViewList(); - d.setBuffer(&alphaDensity); + auto const& [F, M, d, c, particles] = pops[1].getCompileTimeResourcesViewList(); + d.setBuffer(&alphaParticleDensity); + c.setBuffer(&alphaChargeDensity); alpha_M.set_on(M); alphaF.set_on(F); particles.setBuffer(&alphaPack); @@ -515,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{}); @@ -526,7 +540,7 @@ struct IonUpdaterTest : public ::testing::Test PHARE::core::LevelGhostDeposit{}); - ions.computeDensity(); + ions.computeChargeDensity(); ions.computeBulkVelocity(); } // end Ctor @@ -539,17 +553,17 @@ struct IonUpdaterTest : public ::testing::Test for (auto& pop : this->ions) { - interpolate(makeIndexRange(pop.patchGhostParticles()), pop.density(), pop.flux(), - layout); + interpolate(makeIndexRange(pop.patchGhostParticles()), pop.particleDensity(), + pop.chargeDensity(), pop.flux(), layout); double alpha = 0.5; - interpolate(makeIndexRange(pop.levelGhostParticlesNew()), pop.density(), pop.flux(), - layout, + interpolate(makeIndexRange(pop.levelGhostParticlesNew()), pop.particleDensity(), + pop.chargeDensity(), pop.flux(), layout, /*coef = */ alpha); - interpolate(makeIndexRange(pop.levelGhostParticlesOld()), pop.density(), pop.flux(), - layout, + interpolate(makeIndexRange(pop.levelGhostParticlesOld()), pop.particleDensity(), + pop.chargeDensity(), pop.flux(), layout, /*coef = */ (1. - alpha)); } } @@ -560,15 +574,17 @@ struct IonUpdaterTest : public ::testing::Test { auto& populations = this->ions.getRunTimeResourcesViewList(); - auto& protonDensity = populations[0].density(); - 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& 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& alphaDensity = populations[1].density(); - 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& 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 ix0 = this->layout.physicalStartIndex(QtyCentering::primal, Direction::X); auto ix1 = this->layout.physicalEndIndex(QtyCentering::primal, Direction::X); @@ -597,19 +613,18 @@ struct IonUpdaterTest : public ::testing::Test } }; - check(protonDensity, ionsBufferCpy.protonDensity); - - + check(protonParticleDensity, ionsBufferCpy.protonParticleDensity); + check(protonChargeDensity, ionsBufferCpy.protonChargeDensity); check(protonFx, ionsBufferCpy.protonF(Component::X)); check(protonFy, ionsBufferCpy.protonF(Component::Y)); check(protonFz, ionsBufferCpy.protonF(Component::Z)); - check(alphaDensity, ionsBufferCpy.alphaDensity); + check(alphaParticleDensity, ionsBufferCpy.alphaParticleDensity); + check(alphaChargeDensity, ionsBufferCpy.alphaChargeDensity); check(alphaFx, ionsBufferCpy.alphaF(Component::X)); check(alphaFy, ionsBufferCpy.alphaF(Component::Y)); check(alphaFz, ionsBufferCpy.alphaF(Component::Z)); - check(ions.density(), ionsBufferCpy.ionDensity); 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)); @@ -650,12 +665,12 @@ struct IonUpdaterTest : public ::testing::Test } }; - auto& populations = this->ions.getRunTimeResourcesViewList(); - auto& protonDensity = populations[0].density(); - auto& alphaDensity = populations[1].density(); + auto& populations = this->ions.getRunTimeResourcesViewList(); + auto& protonParticleDensity = populations[0].particleDensity(); + auto& alphaParticleDensity = populations[1].particleDensity(); - check(protonDensity, density); - check(alphaDensity, density); + check(protonParticleDensity, density); + check(alphaParticleDensity, density); } }; @@ -909,7 +924,7 @@ TYPED_TEST(IonUpdaterTest, thatNoNaNsExistOnPhysicalNodesMoments) { for (auto ix = ix0; ix <= ix1; ++ix) { - auto& density = pop.density(); + auto& density = pop.particleDensity(); auto& flux = pop.flux(); auto& fx = flux.getComponent(Component::X); diff --git a/tests/diagnostic/__init__.py b/tests/diagnostic/__init__.py index fb7c52d96..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 ["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 efcfbed74..dcd2b1159 100644 --- a/tests/diagnostic/test_diagnostics.hpp +++ b/tests/diagnostic/test_diagnostics.hpp @@ -122,10 +122,10 @@ void validateFluidDump(Simulator& sim, Hi5Diagnostic& hi5) auto& ions = hi5.modelView.getIons(); for (auto& pop : ions) { - checkF(layout, path, "/ions/pop/" + pop.name(), "/density"s, pop.density()); + checkF(layout, path, "/ions/pop/" + pop.name(), "/density"s, pop.chargeDensity()); checkVF(layout, path, "/ions/pop/" + pop.name(), "/flux"s, pop.flux()); } - checkF(layout, path, "/ions"s, "/density"s, ions.density()); + checkF(layout, path, "/ions"s, "/charge_density"s, ions.chargeDensity()); std::string tree{"/ions"}, var{"/bulkVelocity"}; auto hifile = hi5.writer.makeFile(hi5.writer.fileString(tree + var), hi5.flags_); @@ -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/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/diagnostic/test_diagnostics.ipp b/tests/diagnostic/test_diagnostics.ipp index b6813b165..e108c465b 100644 --- a/tests/diagnostic/test_diagnostics.ipp +++ b/tests/diagnostic/test_diagnostics.ipp @@ -14,7 +14,7 @@ void fluid_test(Simulator&& sim, std::string out_dir) { // Scoped to destruct after dump Hi5Diagnostic hi5{hierarchy, hybridModel, out_dir, NEW_HI5_FILE}; - hi5.dMan.addDiagDict(hi5.fluid("/ions/density")) + hi5.dMan.addDiagDict(hi5.fluid("/ions/charge_density")) .addDiagDict(hi5.fluid("/ions/bulkVelocity")) .addDiagDict(hi5.fluid("/ions/momentum_tensor")) .addDiagDict(hi5.fluid("/ions/pop/alpha/density")) diff --git a/tests/functional/alfven_wave/alfven_wave1d.py b/tests/functional/alfven_wave/alfven_wave1d.py index dc830ce13..1fe582bc1 100644 --- a/tests/functional/alfven_wave/alfven_wave1d.py +++ b/tests/functional/alfven_wave/alfven_wave1d.py @@ -96,7 +96,7 @@ def vthz(x): for quantity in ["E", "B"]: ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) - for quantity in ["density", "bulkVelocity"]: + for quantity in ["charge_density", "bulkVelocity"]: ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps) return sim diff --git a/tests/functional/dispersion/dispersion.py b/tests/functional/dispersion/dispersion.py index ad6619984..67ab8d831 100644 --- a/tests/functional/dispersion/dispersion.py +++ b/tests/functional/dispersion/dispersion.py @@ -171,7 +171,7 @@ def vthz(x): quantity=quantity, write_timestamps=timestamps, flush_every=0 ) - for quantity in ["density", "bulkVelocity"]: + for quantity in ["charge_density", "bulkVelocity"]: ph.FluidDiagnostics( quantity=quantity, write_timestamps=timestamps, flush_every=0 ) diff --git a/tests/functional/harris/harris_2d_lb.py b/tests/functional/harris/harris_2d_lb.py index 1153b6e8b..5b553ffb6 100644 --- a/tests/functional/harris/harris_2d_lb.py +++ b/tests/functional/harris/harris_2d_lb.py @@ -134,11 +134,11 @@ def vthxyz(x, y): for quantity in ["E", "B"]: ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) - for quantity in ["density", "bulkVelocity"]: + for quantity in ["charge_density", "bulkVelocity"]: ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps) ph.FluidDiagnostics( - quantity="density", write_timestamps=timestamps, population_name="protons" + quantity="charge_density", write_timestamps=timestamps, population_name="protons" ) ph.InfoDiagnostics(quantity="particle_count") diff --git a/tests/functional/shock/shock.py b/tests/functional/shock/shock.py index 7fd0a6ec2..6cc63f297 100644 --- a/tests/functional/shock/shock.py +++ b/tests/functional/shock/shock.py @@ -107,7 +107,7 @@ def vthz(x): for quantity in ["E", "B"]: ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) - for quantity in ["density", "bulkVelocity"]: + for quantity in ["charge_density", "bulkVelocity"]: ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps) return sim diff --git a/tests/functional/td/td1d.py b/tests/functional/td/td1d.py index 7a715a182..e21d9f28d 100644 --- a/tests/functional/td/td1d.py +++ b/tests/functional/td/td1d.py @@ -100,7 +100,7 @@ def vthz(x): for quantity in ["E", "B"]: ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) - for quantity in ["density", "bulkVelocity"]: + for quantity in ["charge_density", "bulkVelocity"]: ph.FluidDiagnostics( quantity=quantity, write_timestamps=timestamps, diff --git a/tests/functional/tdtagged/td1dtagged.py b/tests/functional/tdtagged/td1dtagged.py index 2d30e59f5..35ab9ec36 100644 --- a/tests/functional/tdtagged/td1dtagged.py +++ b/tests/functional/tdtagged/td1dtagged.py @@ -120,7 +120,7 @@ def config(**options): for quantity in ["E", "B"]: ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) - for quantity in ["density", "bulkVelocity"]: + for quantity in ["charge_density", "bulkVelocity"]: ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps) for pop in sim.model.populations: diff --git a/tests/functional/translation/translat1d.py b/tests/functional/translation/translat1d.py index 53eb7daa0..6a5321698 100644 --- a/tests/functional/translation/translat1d.py +++ b/tests/functional/translation/translat1d.py @@ -87,7 +87,7 @@ def vthz(x): for quantity in ["E", "B"]: ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) - for quantity in ["density", "bulkVelocity"]: + for quantity in ["charge_density", "bulkVelocity"]: ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps) return sim @@ -182,7 +182,7 @@ def vthz(x): for quantity in ["E", "B"]: ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) - for quantity in ["density", "bulkVelocity"]: + for quantity in ["charge_density", "bulkVelocity"]: ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps) return sim diff --git a/tests/simulator/initialize/CMakeLists.txt b/tests/simulator/initialize/CMakeLists.txt index cbaf44d02..abda5a366 100644 --- a/tests/simulator/initialize/CMakeLists.txt +++ b/tests/simulator/initialize/CMakeLists.txt @@ -20,4 +20,6 @@ if(HighFive) endif() + phare_mpi_python3_exec(9 ${PHARE_MPI_PROCS} init-densities density_check.py ${CMAKE_CURRENT_BINARY_DIR}) + endif() diff --git a/tests/simulator/initialize/density_check.py b/tests/simulator/initialize/density_check.py new file mode 100644 index 000000000..bbd54799d --- /dev/null +++ b/tests/simulator/initialize/density_check.py @@ -0,0 +1,281 @@ +import os + +from pyphare.simulator.simulator import Simulator +from pyphare.pharesee.run import Run + +import pyphare.pharein as ph + +from pyphare.pharein import global_vars +from tests.diagnostic import all_timestamps + +from pyphare.pharein.diagnostics import FluidDiagnostics + +import matplotlib.pyplot as plt +import matplotlib as mpl +import numpy as np + +from pyphare.pharesee.hierarchy import hierarchy_from +from pyphare.pharesee.hierarchy.fromfunc import ions_mass_density_func1d +from pyphare.pharesee.hierarchy.fromfunc import ions_charge_density_func1d +from pyphare.pharesee.hierarchy.fromfunc import ions_mass_density_func2d +from pyphare.pharesee.hierarchy.fromfunc import ions_charge_density_func2d + +mpl.use("Agg") + + +ncell = 100 +dl = 0.2 +L = ncell*dl +ts = 0.01 +masses=(2, 3) +charges=(1, 2) + + + +def densityMain_1d(x): + return 1.0 + +def densityBeam_1d(x): + u = x/L-0.5 + return np.exp(-u**2) + +def bx_1d(x): + return 1.0 + +def by_1d(x): + return 0.0 + +def bz_1d(x): + return 0.0 + +def v0_1d(x): + return 0.0 + +def vth_1d(x): + return np.sqrt(1.0) + + +def config_1d(): + + sim = ph.Simulation( + smallest_patch_size=20, + largest_patch_size=60, + time_step=ts, + time_step_nbr=1, + boundary_types="periodic", + cells=ncell, + dl=dl, + diag_options={ + "format": "phareh5", + "options": {"dir": "nCheck_1d", "mode": "overwrite"}, + }, + ) + + v_pop = { + "vbulkx": v0_1d, + "vbulky": v0_1d, + "vbulkz": v0_1d, + "vthx": vth_1d, + "vthy": vth_1d, + "vthz": vth_1d, + } + + ph.MaxwellianFluidModel( + bx=bx_1d, + by=by_1d, + bz=bz_1d, + main={"mass": masses[0], "charge": charges[0], "density": densityMain_1d, "nbr_part_per_cell": 1000, **v_pop}, + beam={"mass": masses[1], "charge": charges[1], "density": densityBeam_1d, "nbr_part_per_cell": 1000, **v_pop}, + ) + + ph.ElectronModel(closure="isothermal", Te=0.0) + + timestamps = all_timestamps(global_vars.sim) + + for quantity in ["charge_density", "mass_density"]: + FluidDiagnostics( + quantity=quantity, + write_timestamps=timestamps + ) + + poplist = ["main", "beam"] + for pop in poplist: + for quantity in ["density", "charge_density"]: + FluidDiagnostics( + quantity=quantity, + write_timestamps=timestamps, + population_name=pop, + ) + + return sim + + + +def densityMain_2d(x, y): + assert len(x) == len(y) + return 1.0*np.ones_like(x) + +def densityBeam_2d(x, y): + assert len(x) == len(y) + u = x/L-0.5 + v = y/L-0.5 + return np.exp(-u**2-v**2) + +def bx_2d(x, y): + return 1.0 + +def by_2d(x, y): + return 0.0 + +def bz_2d(x, y): + return 0.0 + +def v0_2d(x, y): + return 0.0 + +def vth_2d(x, y): + return np.sqrt(1.0) + + +def config_2d(): + + sim = ph.Simulation( + smallest_patch_size=20, + largest_patch_size=60, + time_step=ts, + time_step_nbr=1, + boundary_types=("periodic", "periodic"), + cells=(ncell, ncell), + dl=(dl, dl), + diag_options={ + "format": "phareh5", + "options": {"dir": "nCheck_2d", "mode": "overwrite"}, + }, + ) + + v_pop = { + "vbulkx": v0_2d, + "vbulky": v0_2d, + "vbulkz": v0_2d, + "vthx": vth_2d, + "vthy": vth_2d, + "vthz": vth_2d, + } + + ph.MaxwellianFluidModel( + bx=bx_2d, + by=by_2d, + bz=bz_2d, + main={"mass": masses[0], "charge": charges[0], "density": densityMain_2d, "nbr_part_per_cell": 1000, **v_pop}, + beam={"mass": masses[1], "charge": charges[1], "density": densityBeam_2d, "nbr_part_per_cell": 1000, **v_pop}, + ) + + ph.ElectronModel(closure="isothermal", Te=0.0) + + timestamps = all_timestamps(global_vars.sim) + + for quantity in ["charge_density", "mass_density"]: + FluidDiagnostics( + quantity=quantity, + write_timestamps=timestamps + ) + + poplist = ["main", "beam"] + for pop in poplist: + for quantity in ["density", "charge_density"]: + FluidDiagnostics( + quantity=quantity, + write_timestamps=timestamps, + population_name=pop, + ) + + return sim + + + +def main(): + + Simulator(config_1d()).run().reset() + ph.global_vars.sim = None + Simulator(config_2d()).run().reset() + + def assert_close_enough(h, H): + for lvl_h, lvl_H in zip(h.levels(time).values(), H.levels(time).values()): + for patch_h, patch_H in zip(lvl_h.patches, lvl_H.patches): + pd_h = patch_h.patch_datas["value"] + pd_H = patch_H.patch_datas["value"] + ghosts_num = pd_h.ghosts_nbr[0] + + if pd_H.ndim == 1: + dset_h = pd_h.dataset[ghosts_num:-ghosts_num] + dset_H = pd_H.dataset[ghosts_num:-ghosts_num] + if pd_H.ndim == 2: + dset_h = pd_h.dataset[ghosts_num:-ghosts_num, ghosts_num:-ghosts_num] + dset_H = pd_H.dataset[ghosts_num:-ghosts_num, ghosts_num:-ghosts_num] + + std = np.std(dset_h - dset_H) + print("dim = {}, sigma(user v - actual v) = {}".format(pd_H.ndim, std)) + assert( std < 0.06 ) # empirical value obtained from print just above + + # for h_, H_ in zip(dset_h, dset_H): + # np.testing.assert_almost_equal(h_, H_, decimal=1) + + + fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(3, 2, figsize=(6, 8)) + + # 1d stuffs + run_path = os.path.join(os.curdir, "nCheck_1d") + time = 0.0 + r = Run(run_path) + + h1 = r.GetMassDensity(time) + h2 = r.GetNi(time) + + H1 = hierarchy_from(hier=h1, func=ions_mass_density_func1d, masses=masses, densities=(densityMain_1d, densityBeam_1d)) + H2 = hierarchy_from(hier=h2, func=ions_charge_density_func1d, charges=charges, densities=(densityMain_1d, densityBeam_1d)) + + assert_close_enough(h1, H1) + assert_close_enough(h2, H2) + + cycle = plt.rcParams['axes.prop_cycle'].by_key()['color'] + + h1.plot(ax=ax1, ls="-", lw=2.0, color=cycle[0]) + H1.plot(ax=ax1, ls="-", lw=2.0, color=cycle[1]) + + h2.plot(ax=ax2, ls="-", lw=2.0, color=cycle[0]) + H2.plot(ax=ax2, ls="-", lw=2.0, color=cycle[1]) + + ax1.set_title("mass density : 1d") + ax2.set_title("charge density : 1d") + + + # 2d stuffs + run_path = os.path.join(os.curdir, "nCheck_2d") + time = 0.0 + r = Run(run_path) + + h1 = r.GetMassDensity(time) + h2 = r.GetNi(time) + + H1 = hierarchy_from(hier=h1, func=ions_mass_density_func2d, masses=masses, densities=(densityMain_2d, densityBeam_2d)) + H2 = hierarchy_from(hier=h2, func=ions_charge_density_func2d, charges=charges, densities=(densityMain_2d, densityBeam_2d)) + + assert_close_enough(h1, H1) + assert_close_enough(h2, H2) + + cmap = mpl.colormaps['viridis'] + + h1.plot(ax=ax3, vmin=3.75, vmax=5, cmap=cmap, title="computed mass density : 2d") + H1.plot(ax=ax4, vmin=3.75, vmax=5, cmap=cmap, title="expected mass density : 2d") + h2.plot(ax=ax5, vmin=2.2, vmax=3, cmap=cmap, title="computed charge density : 2d") + H2.plot(ax=ax6, vmin=2.2, vmax=3, cmap=cmap, title="expected charge density : 2d") + + plt.tight_layout() + plt.savefig("nCheck.pdf", dpi=300) + + + + # /home/smets/codes/far/PHARE/tests/simulator/initialize + +if __name__ == "__main__": + main() diff --git a/tests/simulator/refinement/test_2d_10_core.py b/tests/simulator/refinement/test_2d_10_core.py index 031c8cfba..2150c89ec 100644 --- a/tests/simulator/refinement/test_2d_10_core.py +++ b/tests/simulator/refinement/test_2d_10_core.py @@ -101,7 +101,7 @@ def vthz(x, y): for quantity in ["E", "B"]: ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) - for quantity in ["density", "bulkVelocity"]: + for quantity in ["charge_density", "bulkVelocity"]: ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps) return sim @@ -137,7 +137,7 @@ def make_fig(hier, fig_name, ilvl, extra_collections=[]): collections = [ { "boxes": l0_in_l1, - "facecolor": "grey", + "value": 1, } ] if 1 in hier.levels(): @@ -145,7 +145,7 @@ def make_fig(hier, fig_name, ilvl, extra_collections=[]): collections += [ { "boxes": l1_over_l0, - "facecolor": "yellow", + "value": 2, } ] hier.plot_2d_patches( @@ -179,7 +179,7 @@ def post_advance(new_time): extra_collections += [ { "boxes": errors, - "facecolor": "black", + "value": 1, } ] make_fig(L0L1_datahier, L0L1_diags.split("/")[-1], 1, extra_collections) diff --git a/tests/simulator/refinement/test_2d_2_core.py b/tests/simulator/refinement/test_2d_2_core.py index 6a33830ca..e7d72ff6c 100644 --- a/tests/simulator/refinement/test_2d_2_core.py +++ b/tests/simulator/refinement/test_2d_2_core.py @@ -101,7 +101,7 @@ def vthz(x, y): for quantity in ["E", "B"]: ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) - for quantity in ["density", "bulkVelocity"]: + for quantity in ["charge_density", "bulkVelocity"]: ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps) return sim @@ -113,7 +113,7 @@ def make_fig(hier, fig_name, ilvl, extra_collections=[]): collections = [ { "boxes": l0_in_l1, - "facecolor": "grey", + "value": 1, } ] if 1 in hier.levels(): @@ -121,7 +121,7 @@ def make_fig(hier, fig_name, ilvl, extra_collections=[]): collections += [ { "boxes": l1_over_l0, - "facecolor": "yellow", + "value": 2, } ] hier.plot_2d_patches( @@ -171,7 +171,7 @@ def post_advance_1(new_time): extra_collections += [ { "boxes": errors, - "facecolor": "black", + "value": 1, } ] make_fig(L0L1_datahier, L0L1_diags.split("/")[-1], 1, extra_collections) diff --git a/tests/simulator/test_advance.py b/tests/simulator/test_advance.py index 4f385b278..4c0368937 100644 --- a/tests/simulator/test_advance.py +++ b/tests/simulator/test_advance.py @@ -152,7 +152,7 @@ def vthz(*xyz): for quantity in ["E", "B"]: ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) - for quantity in ["density", "bulkVelocity"]: + for quantity in ["charge_density", "bulkVelocity"]: FluidDiagnostics( quantity=quantity, write_timestamps=timestamps, @@ -222,7 +222,7 @@ def vthz(*xyz): h5_filename=diag_outputs + "/EM_B.h5", hier=eb_hier ) mom_hier = hierarchy_from( - h5_filename=diag_outputs + "/ions_density.h5", hier=eb_hier + h5_filename=diag_outputs + "/ions_charge_density.h5", hier=eb_hier ) mom_hier = hierarchy_from( h5_filename=diag_outputs + "/ions_bulkVelocity.h5", hier=mom_hier @@ -231,7 +231,7 @@ def vthz(*xyz): if qty == "moments" or qty == "fields": mom_hier = hierarchy_from( - h5_filename=diag_outputs + "/ions_density.h5", hier=eb_hier + h5_filename=diag_outputs + "/ions_charge_density.h5", hier=eb_hier ) mom_hier = hierarchy_from( h5_filename=diag_outputs + "/ions_bulkVelocity.h5", hier=mom_hier @@ -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], "value": 1}, + {"boxes": [pd2.box], "value": 2}, + {"coords": [errors], "value": 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()) diff --git a/tests/simulator/test_initialization.py b/tests/simulator/test_initialization.py index 2ea135909..d545080f2 100644 --- a/tests/simulator/test_initialization.py +++ b/tests/simulator/test_initialization.py @@ -170,7 +170,7 @@ def vthz(*xyz): quantity=quantity, write_timestamps=np.zeros(time_step_nbr) ) - for quantity in ["density", "bulkVelocity"]: + for quantity in ["charge_density", "bulkVelocity"]: FluidDiagnostics( quantity=quantity, write_timestamps=np.zeros(time_step_nbr) ) @@ -232,7 +232,7 @@ def vthz(*xyz): return particle_hier if qty == "moments": - mom_hier = hierarchy_from(h5_filename=diag_outputs + "/ions_density.h5") + mom_hier = hierarchy_from(h5_filename=diag_outputs + "/ions_charge_density.h5") mom_hier = hierarchy_from( h5_filename=diag_outputs + "/ions_bulkVelocity.h5", hier=mom_hier ) @@ -472,6 +472,7 @@ def _test_density_is_as_provided_by_user(self, dim, interp_order): print("patch {}".format(ip)) ion_density = patch.patch_datas["rho"].dataset[:] + # print(patch.patch_datas.items()) proton_density = patch.patch_datas["protons_rho"].dataset[:] beam_density = patch.patch_datas["beam_rho"].dataset[:] x = patch.patch_datas["rho"].x diff --git a/tests/simulator/test_run.py b/tests/simulator/test_run.py index f935e9768..0b8accc37 100644 --- a/tests/simulator/test_run.py +++ b/tests/simulator/test_run.py @@ -142,7 +142,7 @@ def vthz(x, y): for quantity in ["E", "B"]: ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) - for quantity in ["density", "bulkVelocity"]: + for quantity in ["charge_density", "bulkVelocity"]: ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps) pop = "protons"