Skip to content

Commit

Permalink
add compute the total ion momentum tensor diag
Browse files Browse the repository at this point in the history
  • Loading branch information
@nicolasaunai
nicolasaunai committed Nov 8, 2023
1 parent 0c885ea commit fcfb6da
Show file tree
Hide file tree
Showing 4 changed files with 122 additions and 30 deletions.
17 changes: 9 additions & 8 deletions src/amr/physical_models/hybrid_model.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,14 +29,15 @@ class HybridModel : public IPhysicalModel<AMR_Types>
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<gridlayout_type>;
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<gridlayout_type>;
static constexpr auto dimension = GridLayoutT::dimension;
static constexpr auto interp_order = GridLayoutT::interp_order;
using ParticleInitializerFactory
= core::ParticleInitializerFactory<particle_array_type, gridlayout_type>;

Expand Down
6 changes: 3 additions & 3 deletions src/core/data/ions/ions.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<typename field_type::type>{});
std::transform(std::begin(**mij), std::end(**mij), std::begin(**p_mij),
std::begin(**mij), std::plus<typename field_type::type>{});
}
}
}
Expand Down
56 changes: 48 additions & 8 deletions src/core/numerics/interpolator/interpolator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
{
Expand Down Expand Up @@ -369,6 +370,7 @@ class ParticleToMesh<3>
template<std::size_t dim, std::size_t interpOrder>
class Interpolator : private Weighter<interpOrder>
{
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
Expand Down Expand Up @@ -478,18 +480,11 @@ class Interpolator : private Weighter<interpOrder>
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_<QtyCentering, QtyCentering::primal>(layout, currPart->iCell,
currPart->delta);

Expand All @@ -510,6 +505,8 @@ class Interpolator : private Weighter<interpOrder>
}




/**
* @brief Given a delta and an interpolation order, deduce which lower index to start
* traversing from
Expand Down Expand Up @@ -547,7 +544,7 @@ class Interpolator : private Weighter<interpOrder>
}


private:
protected:
static_assert(dimension <= 3 && dimension > 0 && interpOrder >= 1 && interpOrder <= 3, "error");

using Starts = std::array<std::uint32_t, dimension>;
Expand All @@ -564,6 +561,49 @@ class Interpolator : private Weighter<interpOrder>
Weights primal_weights_;
};

template<std::size_t dim, std::size_t interpOrder>
class MomentumTensorInterpolator : public Interpolator<dim, interpOrder>
{
public:
template<typename ParticleRange, typename TensorField, typename GridLayout>
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_<QtyCentering, QtyCentering::primal>(
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

Expand Down
73 changes: 62 additions & 11 deletions src/diagnostic/detail/types/fluid.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand Down Expand Up @@ -33,6 +34,10 @@ class FluidDiagnosticWriter : public H5TypeWriter<H5Writer>
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}
{
Expand Down Expand Up @@ -68,25 +73,71 @@ class FluidDiagnosticWriter : public H5TypeWriter<H5Writer>
template<typename H5Writer>
void FluidDiagnosticWriter<H5Writer>::compute(DiagnosticProperties& diagnostic)
{
auto& h5Writer = this->h5Writer_;
auto& ions = h5Writer.modelView().getIons();
core::MomentumTensorInterpolator<dimension, interp_order> 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<typename H5Writer>
void FluidDiagnosticWriter<H5Writer>::createFiles(DiagnosticProperties& diagnostic)
{
Expand Down Expand Up @@ -147,17 +198,17 @@ void FluidDiagnosticWriter<H5Writer>::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/"};
if (isActiveDiag(diagnostic, tree, "density"))
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"]);
}


Expand Down

0 comments on commit fcfb6da

Please sign in to comment.