From fcfb6da6df31a5fef2ba095d626e665579473a10 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Tue, 7 Nov 2023 15:00:34 +0100 Subject: [PATCH] add compute the total ion momentum tensor diag --- src/amr/physical_models/hybrid_model.hpp | 17 +++-- src/core/data/ions/ions.hpp | 6 +- .../numerics/interpolator/interpolator.hpp | 56 ++++++++++++-- src/diagnostic/detail/types/fluid.hpp | 73 ++++++++++++++++--- 4 files changed, 122 insertions(+), 30 deletions(-) diff --git a/src/amr/physical_models/hybrid_model.hpp b/src/amr/physical_models/hybrid_model.hpp index 9f6308666..7c78e09ea 100644 --- a/src/amr/physical_models/hybrid_model.hpp +++ b/src/amr/physical_models/hybrid_model.hpp @@ -29,14 +29,15 @@ class HybridModel : public IPhysicalModel using patch_t = typename AMR_Types::patch_t; using level_t = typename AMR_Types::level_t; static const std::string model_name; - using gridlayout_type = GridLayoutT; - using electromag_type = Electromag; - using vecfield_type = typename Electromag::vecfield_type; - using field_type = typename vecfield_type::field_type; - using ions_type = Ions; - using particle_array_type = typename Ions::particle_array_type; - using resources_manager_type = amr::ResourcesManager; - static constexpr auto dimension = GridLayoutT::dimension; + using gridlayout_type = GridLayoutT; + using electromag_type = Electromag; + using vecfield_type = typename Electromag::vecfield_type; + using field_type = typename vecfield_type::field_type; + using ions_type = Ions; + using particle_array_type = typename Ions::particle_array_type; + using resources_manager_type = amr::ResourcesManager; + static constexpr auto dimension = GridLayoutT::dimension; + static constexpr auto interp_order = GridLayoutT::interp_order; using ParticleInitializerFactory = core::ParticleInitializerFactory; diff --git a/src/core/data/ions/ions.hpp b/src/core/data/ions/ions.hpp index 8a388ab9b..b948e0ab5 100644 --- a/src/core/data/ions/ions.hpp +++ b/src/core/data/ions/ions.hpp @@ -151,11 +151,11 @@ namespace core for (auto& pop : populations_) { auto& p_mom = pop.momentumTensor(); - for (auto &p_mij = p_mom.begin(), mij = mom.begin(); p_mij != p_mom.end(); + for (auto p_mij = p_mom.begin(), mij = mom.begin(); p_mij != p_mom.end(); ++p_mij, ++mij) { - std::transform(std::begin(mij), std::end(mij), std::begin(p_mij), - std::begin(mij), std::plus{}); + std::transform(std::begin(**mij), std::end(**mij), std::begin(**p_mij), + std::begin(**mij), std::plus{}); } } } diff --git a/src/core/numerics/interpolator/interpolator.hpp b/src/core/numerics/interpolator/interpolator.hpp index c11a8281a..66c3148f3 100644 --- a/src/core/numerics/interpolator/interpolator.hpp +++ b/src/core/numerics/interpolator/interpolator.hpp @@ -11,6 +11,7 @@ #include "core/logger.hpp" #include "core/hybrid/hybrid_quantities.hpp" #include "core/data/grid/gridlayoutdefs.hpp" +#include "core/utilities/range/range.hpp" namespace PHARE::core { @@ -369,6 +370,7 @@ class ParticleToMesh<3> template class Interpolator : private Weighter { +protected: // this calculates the startIndex and the nbrPointsSupport() weights for // dual field interpolation and puts this at the corresponding location // in 'startIndex' and 'weights'. For dual fields, the normalizedPosition @@ -478,18 +480,11 @@ class Interpolator : private Weighter auto& weights_ = primal_weights_; auto const& [xFlux, yFlux, zFlux] = flux(); - // for each particle, first calculate the startIndex and weights - // for dual and primal quantities. - // then, knowing the centering (primal or dual) of each electromagnetic - // component, we use Interpol to actually perform the interpolation. - // the trick here is that the StartIndex and weights have only been calculated - // twice, and not for each E,B component. PHARE_LOG_START("ParticleToMesh::operator()"); for (auto currPart = begin; currPart != end; ++currPart) { - // TODO #3375 indexAndWeights_(layout, currPart->iCell, currPart->delta); @@ -510,6 +505,8 @@ class Interpolator : private Weighter } + + /** * @brief Given a delta and an interpolation order, deduce which lower index to start * traversing from @@ -547,7 +544,7 @@ class Interpolator : private Weighter } -private: +protected: static_assert(dimension <= 3 && dimension > 0 && interpOrder >= 1 && interpOrder <= 3, "error"); using Starts = std::array; @@ -564,6 +561,49 @@ class Interpolator : private Weighter Weights primal_weights_; }; +template +class MomentumTensorInterpolator : public Interpolator +{ +public: + template + inline void operator()(ParticleRange&& particleRange, TensorField& momentumTensor, + GridLayout const& layout) + { + auto begin = particleRange.begin(); + auto end = particleRange.end(); + auto& startIndex_ = this->primal_startIndex_; + auto& weights_ = this->primal_weights_; + auto const& [xx, xy, xz, yy, yz, zz] = momentumTensor(); + + PHARE_LOG_START("ParticleToMesh::operator()"); + + for (auto currPart = begin; currPart != end; ++currPart) + { + this->template indexAndWeights_( + layout, currPart->iCell, currPart->delta); + + this->particleToMesh_( + xx, *currPart, [](auto const& part) { return part.v[0] * part.v[0]; }, startIndex_, + weights_); + this->particleToMesh_( + xy, *currPart, [](auto const& part) { return part.v[0] * part.v[1]; }, startIndex_, + weights_); + this->particleToMesh_( + xz, *currPart, [](auto const& part) { return part.v[0] * part.v[2]; }, startIndex_, + weights_); + this->particleToMesh_( + yy, *currPart, [](auto const& part) { return part.v[1] * part.v[1]; }, startIndex_, + weights_); + this->particleToMesh_( + yz, *currPart, [](auto const& part) { return part.v[1] * part.v[2]; }, startIndex_, + weights_); + this->particleToMesh_( + zz, *currPart, [](auto const& part) { return part.v[2] * part.v[2]; }, startIndex_, + weights_); + } + PHARE_LOG_STOP("ParticleToMesh::operator()"); + } +}; } // namespace PHARE::core diff --git a/src/diagnostic/detail/types/fluid.hpp b/src/diagnostic/detail/types/fluid.hpp index 7f964a65f..7033bfa53 100644 --- a/src/diagnostic/detail/types/fluid.hpp +++ b/src/diagnostic/detail/types/fluid.hpp @@ -2,6 +2,7 @@ #define PHARE_DIAGNOSTIC_DETAIL_TYPES_FLUID_HPP #include "diagnostic/detail/h5typewriter.hpp" +#include "core/numerics/interpolator/interpolator.hpp" #include "core/data/vecfield/vecfield_component.hpp" @@ -33,6 +34,10 @@ class FluidDiagnosticWriter : public H5TypeWriter using GridLayout = typename H5Writer::GridLayout; using FloatType = typename H5Writer::FloatType; + static constexpr auto dimension = GridLayout::dimension; + static constexpr auto interp_order = GridLayout::interp_order; + + FluidDiagnosticWriter(H5Writer& h5Writer) : Super{h5Writer} { @@ -68,25 +73,71 @@ class FluidDiagnosticWriter : public H5TypeWriter template void FluidDiagnosticWriter::compute(DiagnosticProperties& diagnostic) { - auto& h5Writer = this->h5Writer_; - auto& ions = h5Writer.modelView().getIons(); + core::MomentumTensorInterpolator interpolator; + + auto& h5Writer = this->h5Writer_; + auto& modelView = h5Writer.modelView(); + auto& ions = modelView.getIons(); + auto minLvl = this->h5Writer_.minLevel; + auto maxLvl = this->h5Writer_.maxLevel; // compute the momentum tensor for each population that requires it // compute for all ions but that requires the computation of all pop std::string tree{"/ions/"}; if (isActiveDiag(diagnostic, tree, "momentum_tensor")) + { + auto computeMomentumTensor + = [&](GridLayout& layout, std::string patchID, std::size_t iLvel) { + for (auto& pop : ions) + { + std::string tree{"/ions/pop/" + pop.name() + "/"}; + auto& pop_momentum_tensor = pop.momentumTensor(); + auto domainParts = core::makeIndexRange(pop.domainParticles()); + auto patchGhostParts = core::makeIndexRange(pop.patchGhostParticles()); + + // use levelGhostPartsNew because levelGhostParts may be empty + // and diags are typically done at sync time where levelGhostPartsNew is + // defined + auto levelGhostParts = core::makeIndexRange(pop.levelGhostParticlesNew()); + + interpolator(domainParts, pop_momentum_tensor, layout); + interpolator(patchGhostParts, pop_momentum_tensor, layout); + interpolator(levelGhostParts, pop_momentum_tensor, layout); + } + ions.computeFullMomentumTensor(); + }; + modelView.visitHierarchy(computeMomentumTensor, minLvl, maxLvl); + } + else // if not computing total momentum tensor, user may want to compute it for some pop + { for (auto& pop : ions) { std::string tree{"/ions/pop/" + pop.name() + "/"}; - auto& pop_momentum_tensor = pop.momentumTensor(); - // use some interpolator to calculate the momentum tensor - // from the domain particles and ghost particles + + auto computePopMomentumTensor + = [&](GridLayout& layout, std::string patchID, std::size_t iLvel) { + auto& pop_momentum_tensor = pop.momentumTensor(); + auto domainParts = core::makeIndexRange(pop.domainParticles()); + auto patchGhostParts = core::makeIndexRange(pop.patchGhostParticles()); + + // use levelGhostPartsNew because levelGhostParts may be empty + // and diags are typically done at sync time where levelGhostPartsNew is + // defined + auto levelGhostParts = core::makeIndexRange(pop.levelGhostParticlesNew()); + + interpolator(domainParts, pop_momentum_tensor, layout); + interpolator(patchGhostParts, pop_momentum_tensor, layout); + interpolator(levelGhostParts, pop_momentum_tensor, layout); + }; + if (isActiveDiag(diagnostic, tree, "momentum_tensor")) + { + modelView.visitHierarchy(computePopMomentumTensor, minLvl, maxLvl); + } } - ions.computeFulMomentumTensor(); + } } - template void FluidDiagnosticWriter::createFiles(DiagnosticProperties& diagnostic) { @@ -147,8 +198,8 @@ void FluidDiagnosticWriter::getDataSetInfo(DiagnosticProperties& diagn infoDS(pop.density(), "density", popAttr); if (isActiveDiag(diagnostic, tree, "flux")) infoVF(pop.flux(), "flux", popAttr); - // if (isActiveDiag(diagnostic, tree, "momentum_tensor")) - // infoTF(momentum_tensor_[pop.name()], "momentum_tensor", popAttr); + if (isActiveDiag(diagnostic, tree, "momentum_tensor")) + infoTF(pop.momentumTensor(), "momentum_tensor", popAttr); } std::string tree{"/ions/"}; @@ -156,8 +207,8 @@ void FluidDiagnosticWriter::getDataSetInfo(DiagnosticProperties& diagn infoDS(ions.density(), "density", patchAttributes[lvlPatchID]["ion"]); if (isActiveDiag(diagnostic, tree, "bulkVelocity")) infoVF(ions.velocity(), "bulkVelocity", patchAttributes[lvlPatchID]["ion"]); - // if (isActiveDiag(diagnostic, tree, "momentum_tensor")) - // infoTF(momentum_tensor_["ion"], "momentum_tensor", patchAttributes[lvlPatchID]["ion"]); + if (isActiveDiag(diagnostic, tree, "momentum_tensor")) + infoTF(ions.momentumTensor(), "momentum_tensor", patchAttributes[lvlPatchID]["ion"]); }