Skip to content

Commit

Permalink
add momentum tensor to ion and ion populations
Browse files Browse the repository at this point in the history
  • Loading branch information
@nicolasaunai
nicolasaunai committed Nov 7, 2023
1 parent 81ead4d commit 15730ea
Show file tree
Hide file tree
Showing 15 changed files with 420 additions and 83 deletions.
13 changes: 13 additions & 0 deletions src/core/data/grid/gridlayoutdefs.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,19 @@ namespace core
static constexpr std::uint32_t iVy = static_cast<std::uint32_t>(HybridQuantity::Scalar::Vy);
static constexpr std::uint32_t iVz = static_cast<std::uint32_t>(HybridQuantity::Scalar::Vz);

static constexpr std::uint32_t iMxx
= static_cast<std::uint32_t>(HybridQuantity::Scalar::Mxx);
static constexpr std::uint32_t iMxy
= static_cast<std::uint32_t>(HybridQuantity::Scalar::Mxy);
static constexpr std::uint32_t iMxz
= static_cast<std::uint32_t>(HybridQuantity::Scalar::Mxz);
static constexpr std::uint32_t iMyy
= static_cast<std::uint32_t>(HybridQuantity::Scalar::Myy);
static constexpr std::uint32_t iMyz
= static_cast<std::uint32_t>(HybridQuantity::Scalar::Myz);
static constexpr std::uint32_t iMzz
= static_cast<std::uint32_t>(HybridQuantity::Scalar::Mzz);

static constexpr std::uint32_t iP = static_cast<std::uint32_t>(HybridQuantity::Scalar::P);
};
} // namespace core
Expand Down
70 changes: 69 additions & 1 deletion src/core/data/grid/gridlayoutimplyee.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,11 +91,25 @@ namespace core
const std::array<QtyCentering, NBR_COMPO> Vz
= {{data.primal, data.primal, data.primal}};

const std::array<QtyCentering, NBR_COMPO> Mxx
= {{data.primal, data.primal, data.primal}};
const std::array<QtyCentering, NBR_COMPO> Mxy
= {{data.primal, data.primal, data.primal}};
const std::array<QtyCentering, NBR_COMPO> Mxz
= {{data.primal, data.primal, data.primal}};
const std::array<QtyCentering, NBR_COMPO> Myy
= {{data.primal, data.primal, data.primal}};
const std::array<QtyCentering, NBR_COMPO> Myz
= {{data.primal, data.primal, data.primal}};
const std::array<QtyCentering, NBR_COMPO> Mzz
= {{data.primal, data.primal, data.primal}};

const std::array<QtyCentering, NBR_COMPO> P = {{data.primal, data.primal, data.primal}};

const std::array<std::array<QtyCentering, NBR_COMPO>,
static_cast<std::size_t>(HybridQuantity::Scalar::count)>
hybridQtyCentering{Bx, By, Bz, Ex, Ey, Ez, Jx, Jy, Jz, Rho, Vx, Vy, Vz, P};
hybridQtyCentering{Bx, By, Bz, Ex, Ey, Ez, Jx, Jy, Jz, Rho,
Vx, Vy, Vz, P, Mxx, Mxy, Mxz, Myy, Myz, Mzz};


return hybridQtyCentering;
Expand Down Expand Up @@ -150,6 +164,18 @@ namespace core
return {{hybridQtyCentering_[gridData_.iVz][gridData_.idirX]}};
case HybridQuantity::Scalar::P:
return {{hybridQtyCentering_[gridData_.iP][gridData_.idirX]}};
case HybridQuantity::Scalar::Mxx:
return {{hybridQtyCentering_[gridData_.iMxx][gridData_.idirX]}};
case HybridQuantity::Scalar::Mxy:
return {{hybridQtyCentering_[gridData_.iMxy][gridData_.idirX]}};
case HybridQuantity::Scalar::Mxz:
return {{hybridQtyCentering_[gridData_.iMxz][gridData_.idirX]}};
case HybridQuantity::Scalar::Myy:
return {{hybridQtyCentering_[gridData_.iMyy][gridData_.idirX]}};
case HybridQuantity::Scalar::Myz:
return {{hybridQtyCentering_[gridData_.iMyz][gridData_.idirX]}};
case HybridQuantity::Scalar::Mzz:
return {{hybridQtyCentering_[gridData_.iMzz][gridData_.idirX]}};
default: throw std::runtime_error("Wrong hybridQuantity");
}
}
Expand Down Expand Up @@ -200,6 +226,24 @@ namespace core
case HybridQuantity::Scalar::P:
return {{hybridQtyCentering_[gridData_.iP][gridData_.idirX],
hybridQtyCentering_[gridData_.iP][gridData_.idirY]}};
case HybridQuantity::Scalar::Mxx:
return {{hybridQtyCentering_[gridData_.iMxx][gridData_.idirX],
hybridQtyCentering_[gridData_.iMxx][gridData_.idirY]}};
case HybridQuantity::Scalar::Mxy:
return {{hybridQtyCentering_[gridData_.iMxy][gridData_.idirX],
hybridQtyCentering_[gridData_.iMxy][gridData_.idirY]}};
case HybridQuantity::Scalar::Mxz:
return {{hybridQtyCentering_[gridData_.iMxz][gridData_.idirX],
hybridQtyCentering_[gridData_.iMxz][gridData_.idirY]}};
case HybridQuantity::Scalar::Myy:
return {{hybridQtyCentering_[gridData_.iMyy][gridData_.idirX],
hybridQtyCentering_[gridData_.iMyy][gridData_.idirY]}};
case HybridQuantity::Scalar::Myz:
return {{hybridQtyCentering_[gridData_.iMyz][gridData_.idirX],
hybridQtyCentering_[gridData_.iMyz][gridData_.idirY]}};
case HybridQuantity::Scalar::Mzz:
return {{hybridQtyCentering_[gridData_.iMzz][gridData_.idirX],
hybridQtyCentering_[gridData_.iMzz][gridData_.idirY]}};
default: throw std::runtime_error("Wrong hybridQuantity");
}
}
Expand Down Expand Up @@ -264,6 +308,30 @@ namespace core
return {{hybridQtyCentering_[gridData_.iP][gridData_.idirX],
hybridQtyCentering_[gridData_.iP][gridData_.idirY],
hybridQtyCentering_[gridData_.iP][gridData_.idirZ]}};
case HybridQuantity::Scalar::Mxx:
return {{hybridQtyCentering_[gridData_.iMxx][gridData_.idirX],
hybridQtyCentering_[gridData_.iMxx][gridData_.idirY],
hybridQtyCentering_[gridData_.iMxx][gridData_.idirZ]}};
case HybridQuantity::Scalar::Mxy:
return {{hybridQtyCentering_[gridData_.iMxy][gridData_.idirX],
hybridQtyCentering_[gridData_.iMxy][gridData_.idirY],
hybridQtyCentering_[gridData_.iMxy][gridData_.idirZ]}};
case HybridQuantity::Scalar::Mxz:
return {{hybridQtyCentering_[gridData_.iMxz][gridData_.idirX],
hybridQtyCentering_[gridData_.iMxz][gridData_.idirY],
hybridQtyCentering_[gridData_.iMxz][gridData_.idirZ]}};
case HybridQuantity::Scalar::Myy:
return {{hybridQtyCentering_[gridData_.iMyy][gridData_.idirX],
hybridQtyCentering_[gridData_.iMyy][gridData_.idirY],
hybridQtyCentering_[gridData_.iMyy][gridData_.idirZ]}};
case HybridQuantity::Scalar::Myz:
return {{hybridQtyCentering_[gridData_.iMyz][gridData_.idirX],
hybridQtyCentering_[gridData_.iMyz][gridData_.idirY],
hybridQtyCentering_[gridData_.iMyz][gridData_.idirZ]}};
case HybridQuantity::Scalar::Mzz:
return {{hybridQtyCentering_[gridData_.iMzz][gridData_.idirX],
hybridQtyCentering_[gridData_.iMzz][gridData_.idirY],
hybridQtyCentering_[gridData_.iMzz][gridData_.idirZ]}};
default: throw std::runtime_error("Wrong hybridQuantity");
}
}
Expand Down
29 changes: 18 additions & 11 deletions src/core/data/ions/ion_population/ion_population.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,21 +6,19 @@
#include <stdexcept>
#include <string>
#include <tuple>
#include <vector>
#include <array>


#include "core/def.hpp"
#include "core/hybrid/hybrid_quantities.hpp"
#include "core/data/ions/particle_initializers/particle_initializer.hpp"
#include "initializer/data_provider.hpp"
#include "particle_pack.hpp"
#include "core/utilities/algorithm.hpp"

namespace PHARE
{
namespace core
{
template<typename ParticleArray, typename VecField, typename GridLayout>
template<typename ParticleArray, typename VecField, typename TensorField, typename GridLayout>
class IonPopulation
{
public:
Expand All @@ -29,12 +27,14 @@ namespace core
using particle_array_type = ParticleArray;
using particle_resource_type = ParticlesPack<ParticleArray>;
using vecfield_type = VecField;
using tensorfield_type = TensorField;


IonPopulation(initializer::PHAREDict const& initializer)
: name_{initializer["name"].template to<std::string>()}
, mass_{initializer["mass"].template to<double>()}
, flux_{name_ + "_flux", HybridQuantity::Vector::V}
, momentumTensor_{name_ + "_momentumTensor", HybridQuantity::Tensor::M}
, particleInitializerInfo_{initializer["particle_initializer"]}
{
}
Expand All @@ -51,13 +51,15 @@ namespace core

NO_DISCARD bool isUsable() const
{
return particles_ != nullptr && rho_ != nullptr && flux_.isUsable();
return particles_ != nullptr && rho_ != nullptr && flux_.isUsable()
&& momentumTensor_.isUsable();
}


NO_DISCARD bool isSettable() const
{
return particles_ == nullptr && rho_ == nullptr && flux_.isSettable();
return particles_ == nullptr && rho_ == nullptr && flux_.isSettable()
&& momentumTensor_.isSettable();
}


Expand Down Expand Up @@ -187,10 +189,11 @@ namespace core


NO_DISCARD VecField const& flux() const { return flux_; }


NO_DISCARD VecField& flux() { return flux_; }

TensorField const& momentumTensor() const { return momentumTensor_; }
TensorField& momentumTensor() { return momentumTensor_; }



//-------------------------------------------------------------------------
Expand All @@ -205,7 +208,7 @@ namespace core
typename HybridQuantity::Scalar qty;
};

using MomentProperties = std::vector<MomentsProperty>;
using MomentProperties = std::array<MomentsProperty, 1>;



Expand All @@ -223,7 +226,7 @@ namespace core



using ParticleProperties = std::vector<ParticleProperty>;
using ParticleProperties = std::array<ParticleProperty, 1>;

NO_DISCARD ParticleProperties getParticleArrayNames() const { return {{{name_}}}; }

Expand Down Expand Up @@ -255,7 +258,10 @@ namespace core



NO_DISCARD auto getCompileTimeResourcesUserList() { return std::forward_as_tuple(flux_); }
NO_DISCARD auto getCompileTimeResourcesUserList()
{
return std::forward_as_tuple(flux_, momentumTensor_);
}


//-------------------------------------------------------------------------
Expand All @@ -277,6 +283,7 @@ namespace core
std::string name_;
double mass_;
VecField flux_;
TensorField momentumTensor_;
field_type* rho_{nullptr};
ParticlesPack<ParticleArray>* particles_{nullptr};
initializer::PHAREDict const& particleInitializerInfo_;
Expand Down
49 changes: 36 additions & 13 deletions src/core/data/ions/ions.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,8 @@
#include <sstream>
#include <string>
#include <cmath>
#include <vector>
#include <array>


#include "core/def.hpp"
Expand All @@ -24,9 +26,9 @@ namespace core
class Ions
{
public:
using field_type = typename IonPopulation::field_type;
using vecfield_type = typename IonPopulation::vecfield_type;
// using tensorfield_type = typename IonPopulation::tensorfield_type;
using field_type = typename IonPopulation::field_type;
using vecfield_type = typename IonPopulation::vecfield_type;
using tensorfield_type = typename IonPopulation::tensorfield_type;
using particle_array_type = typename IonPopulation::particle_array_type;
using ParticleInitializerFactoryT
= ParticleInitializerFactory<particle_array_type, GridLayout>;
Expand All @@ -38,6 +40,7 @@ namespace core

explicit Ions(PHARE::initializer::PHAREDict const& dict)
: bulkVelocity_{"bulkVel", HybridQuantity::Vector::V}
, momentumTensor_{"momentumTensor", HybridQuantity::Tensor::M}
, populations_{}
{
auto nbrPop = dict["nbrPopulations"].template to<int>();
Expand Down Expand Up @@ -88,6 +91,9 @@ namespace core

NO_DISCARD std::string densityName() const { return "rho"; }

tensorfield_type const& momentumTensor() const { return momentumTensor_; }

tensorfield_type& momentumTensor() { return momentumTensor_; }

void computeDensity()
{
Expand Down Expand Up @@ -137,21 +143,37 @@ namespace core
}


// TODO 3347 compute ion bulk velocity
void computeFullMomentumTensor()
{
momentumTensor_.zero();
auto& mom = momentumTensor_;

for (auto& pop : populations_)
{
auto& p_mom = pop.momentumTensor();
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>{});
}
}
}


NO_DISCARD auto begin() { return std::begin(populations_); }
NO_DISCARD auto end() { return std::end(populations_); }

NO_DISCARD auto begin() const { return std::begin(populations_); }
NO_DISCARD auto end() const { return std::end(populations_); }


NO_DISCARD bool isUsable() const
{
bool usable = rho_ != nullptr && bulkVelocity_.isUsable();
bool usable
= rho_ != nullptr and bulkVelocity_.isUsable() and momentumTensor_.isUsable();
for (auto const& pop : populations_)
{
usable = usable && pop.isUsable();
usable = usable and pop.isUsable();
}
return usable;
}
Expand All @@ -160,10 +182,11 @@ namespace core

NO_DISCARD bool isSettable() const
{
bool settable = rho_ == nullptr && bulkVelocity_.isSettable();
bool settable
= rho_ == nullptr and bulkVelocity_.isSettable() and momentumTensor_.isSettable();
for (auto const& pop : populations_)
{
settable = settable && pop.isSettable();
settable = settable and pop.isSettable();
}
return settable;
}
Expand All @@ -181,7 +204,7 @@ namespace core
typename HybridQuantity::Scalar qty;
};

using MomentProperties = std::vector<MomentsProperty>;
using MomentProperties = std::array<MomentsProperty, 1>;

NO_DISCARD MomentProperties getFieldNamesAndQuantities() const
{
Expand Down Expand Up @@ -211,7 +234,7 @@ namespace core

NO_DISCARD auto getCompileTimeResourcesUserList()
{
return std::forward_as_tuple(bulkVelocity_);
return std::forward_as_tuple(bulkVelocity_, momentumTensor_);
}


Expand All @@ -236,8 +259,8 @@ namespace core
private:
field_type* rho_{nullptr};
vecfield_type bulkVelocity_;
std::vector<IonPopulation> populations_; // TODO we have to name this so they are unique
// although only 1 Ions should exist.
tensorfield_type momentumTensor_;
std::vector<IonPopulation> populations_;
};
} // namespace core
} // namespace PHARE
Expand Down
Loading

0 comments on commit 15730ea

Please sign in to comment.