Skip to content

Commit

Permalink
add massDensity
Browse files Browse the repository at this point in the history
  • Loading branch information
@nicolasaunai
nicolasaunai committed Nov 22, 2023
1 parent c5cb8db commit e262158
Show file tree
Hide file tree
Showing 2 changed files with 80 additions and 20 deletions.
93 changes: 73 additions & 20 deletions src/core/data/ions/ions.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<particle_array_type, GridLayout>;
Expand All @@ -48,6 +49,8 @@ namespace core
auto const& pop = dict["pop" + std::to_string(ipop)];
populations_.emplace_back(pop);
}
if (allSameMass_())
sameMasses_ = true;
}


Expand Down Expand Up @@ -102,44 +105,70 @@ namespace core

auto& popDensity = pop.density();
std::transform(std::begin(*rho_), std::end(*rho_), std::begin(popDensity),
std::begin(*rho_), std::plus<typename field_type::type>{});
std::begin(*rho_), std::plus<Float>{});
}
}
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);
auto& vz = bulkVelocity_.getComponent(Component::Z);

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<Float(Float, Float)> plus = std::plus<Float>{};
std::function<Float(Float, Float)> plusMass
= [&pop](Float const& v, Float const& f) { return v + f * pop.mass(); };

auto& flux = pop.flux();
auto const& [fx, fy, fz] = flux();

auto fluxSum = [&](auto const& v, auto const& f) { return v + pop.mass() * f; };

std::transform(std::begin(vx), std::end(vx), std::begin(fx), std::begin(vx),
fluxSum);
sameMasses_ ? plus : plusMass);
std::transform(std::begin(vy), std::end(vy), std::begin(fy), std::begin(vy),
fluxSum);
sameMasses_ ? plus : plusMass);
std::transform(std::begin(vz), std::end(vz), std::begin(fz), std::begin(vz),
fluxSum);
sameMasses_ ? plus : plusMass);
}
// auto denominator = [&](auto const& v, auto const& r) { return v / (pop.mass() * r);
// };


std::transform(std::begin(vx), std::end(vx), std::begin(*rho_), std::begin(vx),
std::divides<typename field_type::type>{});
std::transform(std::begin(vy), std::end(vy), std::begin(*rho_), std::begin(vy),
std::divides<typename field_type::type>{});
std::transform(std::begin(vz), std::end(vz), std::begin(*rho_), std::begin(vz),
std::divides<typename field_type::type>{});
std::transform(std::begin(vx), std::end(vx), std::begin(*density), std::begin(vx),
std::divides<Float>{});
std::transform(std::begin(vy), std::end(vy), std::begin(*density), std::begin(vy),
std::divides<Float>{});
std::transform(std::begin(vz), std::end(vz), std::begin(*density), std::begin(vz),
std::divides<Float>{});
}


Expand All @@ -152,6 +181,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();
Expand Down Expand Up @@ -191,7 +222,8 @@ namespace core

NO_DISCARD MomentProperties getFieldNamesAndQuantities() const
{
return {{{densityName(), HybridQuantity::Scalar::rho}}};
return {{{densityName(), HybridQuantity::Scalar::rho},
{"massDensity", HybridQuantity::Scalar::rho}}};
}


Expand All @@ -202,6 +234,10 @@ namespace core
{
rho_ = field;
}
else if (bufferName == "massDensity")
{
massDensity_ = field;
}
else
{
throw std::runtime_error("Error - invalid density buffer name");
Expand Down Expand Up @@ -240,10 +276,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;
});
}




field_type* rho_{nullptr};
field_type* massDensity_{nullptr};
vecfield_type bulkVelocity_;
std::vector<IonPopulation> populations_; // TODO we have to name this so they are unique
// although only 1 Ions should exist.
std::vector<IonPopulation> 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_;
};
} // namespace core
} // namespace PHARE
Expand Down
7 changes: 7 additions & 0 deletions tests/core/numerics/ion_updater/test_updater.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -227,6 +227,7 @@ struct IonsBuffers
using ParticleInitializerFactory = typename PHARETypes::ParticleInitializerFactory;

Field ionDensity;
Field ionMassDensity;
Field protonDensity;
Field alphaDensity;
Field protonFx;
Expand Down Expand Up @@ -258,6 +259,8 @@ struct IonsBuffers
IonsBuffers(GridLayout const& layout)
: ionDensity{"rho", HybridQuantity::Scalar::rho,
layout.allocSize(HybridQuantity::Scalar::rho)}
, ionMassDensity{"rho", 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,
Expand Down Expand Up @@ -304,6 +307,8 @@ struct IonsBuffers
IonsBuffers(IonsBuffers const& source, GridLayout const& layout)
: ionDensity{"rho", HybridQuantity::Scalar::rho,
layout.allocSize(HybridQuantity::Scalar::rho)}
, ionMassDensity{"rho", 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,
Expand Down Expand Up @@ -336,6 +341,7 @@ struct IonsBuffers

{
ionDensity.copyData(source.ionDensity);
ionMassDensity.copyData(source.ionMassDensity);
protonDensity.copyData(source.protonDensity);
alphaDensity.copyData(source.alphaDensity);

Expand Down Expand Up @@ -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);
Expand Down

0 comments on commit e262158

Please sign in to comment.