From 3abf8c5c1953e96a31f155c98e4bb21a2d7a608e Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Tue, 21 Nov 2023 16:37:54 +0100 Subject: [PATCH] account for population masses to bulk velocity --- .../hybrid_hybrid_messenger_strategy.hpp | 1 + src/core/data/ions/ions.hpp | 106 +++++++++++++----- .../numerics/ion_updater/test_updater.cpp | 7 ++ tests/simulator/test_restarts.py | 37 +++++- 4 files changed, 117 insertions(+), 34 deletions(-) diff --git a/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp b/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp index 5fe33d77f..5e790bd25 100644 --- a/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp +++ b/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp @@ -673,6 +673,7 @@ namespace amr NiOldUser_.name, fieldRefineOp_, fieldTimeOp_, info->modelIonDensity); + velGhostsRefiners_.addTimeRefiners(info->ghostBulkVelocity, info->modelIonBulkVelocity, core::VecFieldNames{ViOld_}, fieldRefineOp_, fieldTimeOp_); diff --git a/src/core/data/ions/ions.hpp b/src/core/data/ions/ions.hpp index 3bd96cc3b..98abf73fd 100644 --- a/src/core/data/ions/ions.hpp +++ b/src/core/data/ions/ions.hpp @@ -28,6 +28,7 @@ namespace core public: using field_type = typename IonPopulation::field_type; using vecfield_type = typename IonPopulation::vecfield_type; + using Float = typename field_type::type; using particle_array_type = typename IonPopulation::particle_array_type; using ParticleInitializerFactoryT = ParticleInitializerFactory; @@ -48,6 +49,7 @@ namespace core auto const& pop = dict["pop" + std::to_string(ipop)]; populations_.emplace_back(pop); } + sameMasses_ = allSameMass_(); } @@ -57,13 +59,9 @@ namespace core NO_DISCARD field_type const& density() const { if (isUsable()) - { return *rho_; - } else - { throw std::runtime_error("Error - cannot access density data"); - } } @@ -71,23 +69,19 @@ namespace core NO_DISCARD field_type& density() { if (isUsable()) - { return *rho_; - } else - { throw std::runtime_error("Error - cannot access density data"); - } } - NO_DISCARD vecfield_type const& velocity() const { return bulkVelocity_; } NO_DISCARD vecfield_type& velocity() { return bulkVelocity_; } NO_DISCARD std::string densityName() const { return "rho"; } + NO_DISCARD std::string massDensityName() const { return "massDensity"; } void computeDensity() @@ -102,13 +96,39 @@ namespace core auto& popDensity = pop.density(); std::transform(std::begin(*rho_), std::end(*rho_), std::begin(popDensity), - std::begin(*rho_), std::plus{}); + std::begin(*rho_), std::plus{}); + } + } + void computeMassDensity() + { + massDensity_->zero(); + + for (auto const& pop : populations_) + { + // we sum over all nodes contiguously, including ghosts + // 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(*massDensity_), std::end(*massDensity_), std::begin(popDensity), + std::begin(*massDensity_), + [&pop](auto const& n, auto const& pop_n) { return n + pop_n * pop.mass(); }); } } 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_; + bulkVelocity_.zero(); auto& vx = bulkVelocity_.getComponent(Component::X); auto& vy = bulkVelocity_.getComponent(Component::Y); @@ -116,29 +136,33 @@ namespace core for (auto& pop : populations_) { - auto& flux = pop.flux(); - auto& fx = flux.getComponent(Component::X); - auto& fy = flux.getComponent(Component::Y); - auto& fz = flux.getComponent(Component::Z); + // 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), - std::plus{}); + sameMasses_ ? plus : plusMass); std::transform(std::begin(vy), std::end(vy), std::begin(fy), std::begin(vy), - std::plus{}); + sameMasses_ ? plus : plusMass); std::transform(std::begin(vz), std::end(vz), std::begin(fz), std::begin(vz), - std::plus{}); + sameMasses_ ? plus : plusMass); } - std::transform(std::begin(vx), std::end(vx), std::begin(*rho_), std::begin(vx), - std::divides{}); - std::transform(std::begin(vy), std::end(vy), std::begin(*rho_), std::begin(vy), - std::divides{}); - std::transform(std::begin(vz), std::end(vz), std::begin(*rho_), std::begin(vz), - std::divides{}); + + std::transform(std::begin(vx), std::end(vx), std::begin(*density), std::begin(vx), + std::divides{}); + std::transform(std::begin(vy), std::end(vy), std::begin(*density), std::begin(vy), + std::divides{}); + std::transform(std::begin(vz), std::end(vz), std::begin(*density), std::begin(vz), + std::divides{}); } - // TODO 3347 compute ion bulk velocity NO_DISCARD auto begin() { return std::begin(populations_); } NO_DISCARD auto end() { return std::end(populations_); } @@ -147,6 +171,8 @@ namespace core NO_DISCARD auto end() const { return std::end(populations_); } + // in the following isUsable and isSettable the massDensity_ is not checked + // because it is for internal use only so no object will ever need to access it. NO_DISCARD bool isUsable() const { bool usable = rho_ != nullptr && bulkVelocity_.isUsable(); @@ -186,20 +212,25 @@ namespace core NO_DISCARD MomentProperties getFieldNamesAndQuantities() const { - return {{{densityName(), HybridQuantity::Scalar::rho}}}; + return {{{densityName(), HybridQuantity::Scalar::rho}, + {massDensityName(), HybridQuantity::Scalar::rho}}}; } void setBuffer(std::string const& bufferName, field_type* field) { - if (bufferName == "rho") + if (bufferName == densityName()) { rho_ = field; } + else if (bufferName == massDensityName()) + { + massDensity_ = field; + } else { - throw std::runtime_error("Error - invalid density buffer name"); + throw std::runtime_error("Error - invalid density buffer name : " + bufferName); } } @@ -235,10 +266,27 @@ namespace core private: + bool allSameMass_() const + { + auto const& firstPop = populations_.front(); + return std::all_of( + std::begin(populations_), std::end(populations_), [&firstPop](auto const& pop) { + return std::abs(pop.mass() - firstPop.mass()) < 1e-10; // arbitrary small diff + }); + } + + + + field_type* rho_{nullptr}; + field_type* massDensity_{nullptr}; vecfield_type bulkVelocity_; - std::vector populations_; // TODO we have to name this so they are unique - // although only 1 Ions should exist. + 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}; }; } // namespace core } // namespace PHARE diff --git a/tests/core/numerics/ion_updater/test_updater.cpp b/tests/core/numerics/ion_updater/test_updater.cpp index f7d177830..47d507e6b 100644 --- a/tests/core/numerics/ion_updater/test_updater.cpp +++ b/tests/core/numerics/ion_updater/test_updater.cpp @@ -227,6 +227,7 @@ struct IonsBuffers using ParticleInitializerFactory = typename PHARETypes::ParticleInitializerFactory; Field ionDensity; + Field ionMassDensity; Field protonDensity; Field alphaDensity; Field protonFx; @@ -258,6 +259,8 @@ struct IonsBuffers IonsBuffers(GridLayout const& layout) : ionDensity{"rho", 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, @@ -304,6 +307,8 @@ struct IonsBuffers IonsBuffers(IonsBuffers const& source, GridLayout const& layout) : ionDensity{"rho", 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, @@ -336,6 +341,7 @@ struct IonsBuffers { ionDensity.copyData(source.ionDensity); + ionMassDensity.copyData(source.ionMassDensity); protonDensity.copyData(source.protonDensity); alphaDensity.copyData(source.alphaDensity); @@ -366,6 +372,7 @@ struct IonsBuffers void setBuffers(Ions& ions) { ions.setBuffer("rho", &ionDensity); + ions.setBuffer("massDensity", &ionMassDensity); auto& v = ions.velocity(); v.setBuffer("bulkVel_x", &Vx); v.setBuffer("bulkVel_y", &Vy); diff --git a/tests/simulator/test_restarts.py b/tests/simulator/test_restarts.py index 140476de4..adc2e769d 100644 --- a/tests/simulator/test_restarts.py +++ b/tests/simulator/test_restarts.py @@ -20,7 +20,7 @@ def permute(dic, expected_num_levels): # from pyphare.pharein.simulation import supported_dimensions # eventually - dims = [1] # supported_dimensions() + dims = [1] # supported_dimensions() return [ [dim, interp, dic, expected_num_levels] for dim in dims for interp in [1, 2, 3] ] @@ -60,6 +60,9 @@ def vy(x): def vz(x): return 0.0 + def vxalpha(x): + return 3.0 + def vthxyz(x): return T(x) @@ -71,6 +74,14 @@ def vthxyz(x): "vthy": vthxyz, "vthz": vthxyz, } + vvvalpha = { + "vbulkx": vxalpha, + "vbulky": vy, + "vbulkz": vz, + "vthx": vthxyz, + "vthy": vthxyz, + "vthz": vthxyz, + } model = ph.MaxwellianFluidModel( bx=bx, by=by, @@ -84,10 +95,10 @@ def vthxyz(x): "init": {"seed": 1337}, }, alpha={ - "mass": 4, + "mass": 4.0, "charge": 1, "density": density, - **vvv, + **vvvalpha, "nbr_part_per_cell": ppc, "init": {"seed": 2334}, }, @@ -160,6 +171,17 @@ def count_levels_and_patches(qty): datahier0.level(0).patches[0].patch_datas.keys() ) + self.assertEqual(len(datahier0.levels()), len(datahier1.levels())) + for ilvl in range(len(datahier0.levels())): + self.assertEqual( + len(datahier0.level(ilvl).patches), + len(datahier1.level(ilvl).patches), + ) + for patch0, patch1 in zip( + datahier0.level(ilvl).patches, datahier1.level(ilvl).patches + ): + self.assertEqual(patch0.box, patch1.box) + self.assertGreater(len(datahier0.levels()), 0) for ilvl, lvl0 in datahier0.levels().items(): @@ -170,7 +192,12 @@ def count_levels_and_patches(qty): pd1 = patch1.patch_datas[pd_key] self.assertNotEqual(id(pd0), id(pd1)) if "domain" in pd_key: - self.assertEqual(pd0.dataset, pd1.dataset) + try: + self.assertEqual(pd0.dataset, pd1.dataset) + except AssertionError: + print( + f"FAILED domain particles at time {time} {ilvl} {patch1.box} {patch0.box}" + ) else: np.testing.assert_equal(pd0.dataset[:], pd1.dataset[:]) checks += 1 @@ -221,7 +248,7 @@ def test_restarts(self, ndim, interp, simInput, expected_num_levels): restart_idx = 4 restart_time = time_step * restart_idx - timestamps = [time_step * restart_idx, time_step * time_step_nbr] + timestamps = [restart_time, time_step * time_step_nbr] # first simulation local_out = self.unique_diag_dir_for_test_case(f"{out}/test", ndim, interp)