Skip to content

Commit

Permalink
Save a span in Quadrature to allow copies to be made
Browse files Browse the repository at this point in the history
  • Loading branch information
EmilyBourne committed Nov 16, 2023
1 parent 1c81b31 commit 20668da
Show file tree
Hide file tree
Showing 21 changed files with 166 additions and 75 deletions.
2 changes: 1 addition & 1 deletion src/geometryXVx/README.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
The `geoemtryXVx` folder contains all the code describing methods which are specific to a geometry with 1 spatial dimension and 1 velocity dimension. It is broken up into the following sub-folders:

<!-- - [boltzmann](./botzmann/README.md) - -->
- [geometry](./geometry/README.md) --> - All the dimension tags used for a simulation in the geoemtry.
- [geometry](./geometry/README.md) --> - All the dimension tags used for a simulation in the geometry.
<!-- - [initialization](./initialization/README.md) - -->
<!-- - [poisson](./poisson/README.md) - Code describing the polar Poisson solver. -->
<!-- - [rhs](./rhs/README.md) - Code describing the operators on the right hand side of the Boltzmann equation; namely sources, sinks and collisions.-->
Expand Down
4 changes: 3 additions & 1 deletion src/geometryXVx/rhs/collisions_inter.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,9 @@ void CollisionsInter::compute_rhs(DSpanSpXVx const rhs, DViewSpXVx const allfdis
IDomainSp const dom_sp(ddc::get_domain<IDimSp>(allfdistribu));
IDomainVx const gridvx(ddc::get_domain<IDimVx>(allfdistribu));

FluidMoments moments(Quadrature<IDimVx>(trapezoid_quadrature_coefficients(gridvx)));
DFieldVx const quadrature_coeffs = trapezoid_quadrature_coefficients(gridvx);
Quadrature<IDimVx> integrate(quadrature_coeffs);
FluidMoments moments(integrate);

ddc::for_each(ddc::get_domain<IDimX>(allfdistribu), [&](IndexX const ix) {
//Moments computation
Expand Down
7 changes: 5 additions & 2 deletions src/geometryXVx/rhs/collisions_intra.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -254,8 +254,11 @@ DSpanSpXVx CollisionsIntra::operator()(DSpanSpXVx allfdistribu, double dt) const
DFieldSpX density(ddc::get_domain<IDimSp, IDimX>(allfdistribu));
DFieldSpX mean_velocity(ddc::get_domain<IDimSp, IDimX>(allfdistribu));
DFieldSpX temperature(ddc::get_domain<IDimSp, IDimX>(allfdistribu));
FluidMoments moments(Quadrature<IDimVx>(
trapezoid_quadrature_coefficients(ddc::get_domain<IDimVx>(allfdistribu))));

DFieldVx const quadrature_coeffs
= trapezoid_quadrature_coefficients(ddc::get_domain<IDimVx>(allfdistribu));
Quadrature<IDimVx> integrate(quadrature_coeffs);
FluidMoments moments(integrate);

moments(density.span_view(), allfdistribu.span_cview(), FluidMoments::s_density);
moments(mean_velocity.span_view(),
Expand Down
5 changes: 3 additions & 2 deletions src/geometryXVx/rhs/collisions_utils.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -106,8 +106,9 @@ void compute_Vcoll_Tcoll(
ddc::ChunkSpan<double const, ddc::DiscreteDomain<IDimSp, IDimX, IDimension>> Dcoll,
ddc::ChunkSpan<double const, ddc::DiscreteDomain<IDimSp, IDimX, IDimension>> dvDcoll)
{
Quadrature<IDimVx> const integrate_v(
trapezoid_quadrature_coefficients(ddc::get_domain<IDimVx>(allfdistribu)));
DFieldVx const quadrature_coeffs
= trapezoid_quadrature_coefficients(ddc::get_domain<IDimVx>(allfdistribu));
Quadrature<IDimVx> const integrate_v(quadrature_coeffs);

// computation of the integrands
DFieldSpXVx I0mean_integrand(allfdistribu.domain());
Expand Down
3 changes: 2 additions & 1 deletion src/geometryXVx/rhs/krook_source_adaptive.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,8 @@ void KrookSourceAdaptive::get_amplitudes(DSpanSp amplitudes, DViewSpVx const all
amplitudes(iion) = m_amplitude;

IDomainVx const gridvx = allfdistribu.domain<IDimVx>();
Quadrature<IDimVx> const integrate_v(trapezoid_quadrature_coefficients(gridvx));
DFieldVx const quadrature_coeffs = trapezoid_quadrature_coefficients(gridvx);
Quadrature<IDimVx> const integrate_v(quadrature_coeffs);
double const density_ion = integrate_v(allfdistribu[iion]);
double const density_electron = integrate_v(allfdistribu[ielec()]);

Expand Down
3 changes: 2 additions & 1 deletion src/geometryXVx/rhs/mask_tanh.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,8 @@ DFieldX mask_tanh(
}

if (normalized) {
Quadrature<IDimX> const integrate_x(trapezoid_quadrature_coefficients(gridx));
DFieldX const quadrature_coeffs = trapezoid_quadrature_coefficients(gridx);
Quadrature<IDimX> const integrate_x(quadrature_coeffs);
double const coeff_norm = integrate_x(mask);
ddc::for_each(gridx, [&](IndexX const ix) { mask(ix) = mask(ix) / coeff_norm; });
}
Expand Down
10 changes: 4 additions & 6 deletions src/geometryXVx/utils/fluid_moments.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,7 @@
#include <quadrature.hpp>
#include <trapezoid_quadrature.hpp>

FluidMoments::FluidMoments(Quadrature<IDimVx> integrate_v) : m_integrate_v(std::move(integrate_v))
{
}
FluidMoments::FluidMoments(Quadrature<IDimVx> integrate_v) : m_integrate_v(integrate_v) {}

/*
* Computes the density of fdistribu
Expand Down Expand Up @@ -37,7 +35,7 @@ void FluidMoments::operator()(
void FluidMoments::operator()(
double& mean_velocity,
DViewVx const fdistribu,
double const& density,
double density,
FluidMoments::MomentVelocity)
{
DFieldVx integrand(fdistribu.domain());
Expand Down Expand Up @@ -74,8 +72,8 @@ void FluidMoments::operator()(
void FluidMoments::operator()(
double& temperature,
DViewVx const fdistribu,
double const& density,
double const& mean_velocity,
double density,
double mean_velocity,
FluidMoments::MomentTemperature)
{
DFieldVx integrand(fdistribu.domain());
Expand Down
95 changes: 86 additions & 9 deletions src/geometryXVx/utils/fluid_moments.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -15,53 +15,130 @@ class FluidMoments
Quadrature<IDimVx> m_integrate_v;

public:
/**
* A tag type to indicate that the density should be calculated.
*/
struct MomentDensity
{
};

/**
* A tag type to indicate that the velocity should be calculated.
*/
struct MomentVelocity
{
};

/**
* A tag type to indicate that the temperature should be calculated.
*/
struct MomentTemperature
{
};

/**
* A static instance of MomentDensity that can be used to indicated to the operator()
* that the density should be calculated.
*/
static constexpr MomentDensity s_density = MomentDensity();
/**
* A static instance of MomentVelocity that can be used to indicated to the operator()
* that the velocity should be calculated.
*/
static constexpr MomentVelocity s_velocity = MomentVelocity();
/**
* A static instance of MomentTemperature that can be used to indicated to the operator()
* that the temperature should be calculated.
*/
static constexpr MomentTemperature s_temperature = MomentTemperature();

/**
* The constructor for the operator.
*
* @param[in] integrate_v A quadrature method which integrates over the velocity space.
*/
FluidMoments(Quadrature<IDimVx> integrate_v);

~FluidMoments() = default;

void operator()(double& density, DViewVx allfdistribu, MomentDensity);
/**
* Calculate the density at a specific point of the distribution function.
*
* @param[out] density The density at the point.
* @param[in] fdistribu A slice in velocity space of the distribution function
* at the given point.
* @param[in] moment_density A tag to ensure that the correct operator is called.
*/
void operator()(double& density, DViewVx fdistribu, MomentDensity moment_density);

void operator()(DSpanSpX density, DViewSpXVx allfdistribu, MomentDensity);
/**
* Calculate the density of the distribution function.
*
* @param[out] density The density at various points for different species.
* @param[in] allfdistribu The distribution function.
* @param[in] moment_density A tag to ensure that the correct operator is called.
*/
void operator()(DSpanSpX density, DViewSpXVx allfdistribu, MomentDensity moment_density);

/**
* Calculate the mean velocity at a specific point of the distribution function.
*
* @param[out] mean_velocity The mean velocity at the point.
* @param[in] fdistribu A slice in velocity space of the distribution function
* at the given point.
* @param[in] density The density at the point.
* @param[in] moment_velocity A tag to ensure that the correct operator is called.
*/
void operator()(
double& mean_velocity,
DViewVx fdistribu,
double const& density,
MomentVelocity);
double density,
MomentVelocity moment_velocity);

/**
* Calculate the mean velocity of the distribution function.
*
* @param[out] mean_velocity The mean velocity at various points for different species.
* @param[in] allfdistribu The distribution function.
* @param[in] density The density at various points for different species.
* @param[in] moment_velocity A tag to ensure that the correct operator is called.
*/
void operator()(
DSpanSpX mean_velocity,
DViewSpXVx allfdistribu,
DViewSpX density,
MomentVelocity);
MomentVelocity moment_velocity);

/**
* Calculate the temperature at a specific point of the distribution function.
*
* @param[out] temperature The mean temperature at the point.
* @param[in] fdistribu A slice in velocity space of the distribution function
* at the given point.
* @param[in] density The density at the point.
* @param[in] mean_velocity The mean velocity at the point.
* @param[in] moment_temperature A tag to ensure that the correct operator is called.
*/
void operator()(
double& temperature,
DViewVx fdistribu,
double const& density,
double const& mean_velocity,
MomentTemperature);
double density,
double mean_velocity,
MomentTemperature moment_temperature);

/**
* Calculate the mean temperature of the distribution function.
*
* @param[out] temperature The mean temperature at various points for different species.
* @param[in] allfdistribu The distribution function.
* @param[in] density The density at various points for different species.
* @param[in] mean_velocity The mean velocity at various points for different species.
* @param[in] moment_temperature A tag to ensure that the correct operator is called.
*/
void operator()(
DSpanSpX temperature,
DViewSpXVx allfdistribu,
DViewSpX density,
DViewSpX mean_velocity,
MomentTemperature);
MomentTemperature moment_temperature);
};
4 changes: 2 additions & 2 deletions src/quadrature/compute_norms.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@
*/
template <class... IDim>
double compute_L1_norm(
Quadrature<IDim...>& quadrature,
Quadrature<IDim...> quadrature,
ddc::ChunkSpan<double, ddc::DiscreteDomain<IDim...>> function)
{
auto grid = ddc::get_domain<IDim...>(function);
Expand Down Expand Up @@ -53,7 +53,7 @@ double compute_L1_norm(
*/
template <class... IDim>
double compute_L2_norm(
Quadrature<IDim...>& quadrature,
Quadrature<IDim...> quadrature,
ddc::ChunkSpan<double, ddc::DiscreteDomain<IDim...>> function)
{
auto grid = ddc::get_domain<IDim...>(function);
Expand Down
14 changes: 3 additions & 11 deletions src/quadrature/quadrature.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -13,27 +13,19 @@ template <class... IDim>
class Quadrature
{
private:
ddc::Chunk<double, ddc::DiscreteDomain<IDim...>> m_coefficients;
ddc::ChunkView<double, ddc::DiscreteDomain<IDim...>> m_coefficients;

public:
/**
* @brief Create a Quadrature object.
* @param coeffs
* The coefficients of the quadrature.
*/
explicit Quadrature(ddc::Chunk<double, ddc::DiscreteDomain<IDim...>>&& coeffs)
: m_coefficients(std::move(coeffs))
explicit Quadrature(ddc::ChunkView<double, ddc::DiscreteDomain<IDim...>> const& coeffs)
: m_coefficients(coeffs)
{
}

/**
* @brief Create a Quadrature object by copy.
* @param rhs The object being copied.
*/
Quadrature(Quadrature&& rhs) = default;

~Quadrature() = default;

/**
* @brief An operator for calculating the integral of a function defined on a discrete domain.
* @param[in] values
Expand Down
14 changes: 6 additions & 8 deletions tests/geometryRTheta/quadrature/tests_L1_and_L2_norms.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ namespace {
* The error tolerance @f$ \varepsilon @f$.
*/
void check_norm_L1(
Quadrature<IDimR, IDimP>& quadrature,
Quadrature<IDimR, IDimP> quadrature,
DSpanRP function,
double const expected_norm,
double const TOL)
Expand Down Expand Up @@ -66,7 +66,7 @@ void check_norm_L1(
* The error tolerance @f$ \varepsilon @f$.
*/
void check_norm_L2(
Quadrature<IDimR, IDimP>& quadrature,
Quadrature<IDimR, IDimP> quadrature,
DSpanRP function,
double const expected_norm,
double const TOL)
Expand All @@ -92,7 +92,7 @@ void check_norm_L2(
* A 2x1 array with the error tolerance @f$ \varepsilon @f$.
*/
void check_norms(
Quadrature<IDimR, IDimP>& quadrature,
Quadrature<IDimR, IDimP> quadrature,
DSpanRP function,
std::array<double, 2> const& expected_norms,
std::array<double, 2> const& TOLs)
Expand All @@ -112,12 +112,10 @@ void launch_tests(
{
// Test spline quadrature: ------------------------------------------------------------------------
// Instantiate a quadrature with coefficients where we added the Jacobian determinant.
Quadrature<IDimR, IDimP> quadrature(compute_coeffs_on_mapping(
DFieldRP const quadrature_coeffs = compute_coeffs_on_mapping(
mapping,
spline_quadrature_coefficients(
grid,
builder.get_builder_1(),
builder.get_builder_2())));
spline_quadrature_coefficients(grid, builder.get_builder_1(), builder.get_builder_2()));
Quadrature<IDimR, IDimP> quadrature(quadrature_coeffs);


DFieldRP test(grid);
Expand Down
4 changes: 3 additions & 1 deletion tests/geometryXVx/collisions_inter.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,9 @@ TEST(CollisionsInter, CollisionsInter)
double const nustar0(0.1);
CollisionsInter collisions(mesh, nustar0);

FluidMoments moments(Quadrature<IDimVx>(trapezoid_quadrature_coefficients(gridvx)));
DFieldVx const quadrature_coeffs = trapezoid_quadrature_coefficients(gridvx);
Quadrature<IDimVx> integrate(quadrature_coeffs);
FluidMoments moments(integrate);

DFieldSpX nustar_profile(ddc::get_domain<IDimSp, IDimX>(allfdistribu));
compute_nustar_profile(nustar_profile.span_view(), nustar0);
Expand Down
6 changes: 4 additions & 2 deletions tests/geometryXVx/collisions_intra_maxwellian.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -197,8 +197,10 @@ TEST(CollisionsIntraMaxwellian, CollisionsIntraMaxwellian)
DFieldSpX density_res(ddc::get_domain<IDimSp, IDimX>(allfdistribu));
DFieldSpX mean_velocity_res(ddc::get_domain<IDimSp, IDimX>(allfdistribu));
DFieldSpX temperature_res(ddc::get_domain<IDimSp, IDimX>(allfdistribu));
FluidMoments moments(Quadrature<IDimVx>(
trapezoid_quadrature_coefficients(ddc::get_domain<IDimVx>(allfdistribu))));
DFieldVx const quadrature_coeffs
= trapezoid_quadrature_coefficients(ddc::get_domain<IDimVx>(allfdistribu));
Quadrature<IDimVx> integrate(quadrature_coeffs);
FluidMoments moments(integrate);

moments(density_res.span_view(), allfdistribu.span_cview(), FluidMoments::s_density);
moments(mean_velocity_res.span_view(),
Expand Down
6 changes: 4 additions & 2 deletions tests/geometryXVx/fluid_moments.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -106,8 +106,10 @@ TEST(Physics, FluidMoments)
DFieldSpX density_computed(ddc::get_domain<IDimSp, IDimX>(allfdistribu));
DFieldSpX mean_velocity_computed(ddc::get_domain<IDimSp, IDimX>(allfdistribu));
DFieldSpX temperature_computed(ddc::get_domain<IDimSp, IDimX>(allfdistribu));
FluidMoments moments(Quadrature<IDimVx>(
trapezoid_quadrature_coefficients(ddc::get_domain<IDimVx>(allfdistribu))));
DFieldVx const quadrature_coeffs
= trapezoid_quadrature_coefficients(ddc::get_domain<IDimVx>(allfdistribu));
Quadrature<IDimVx> integrate(quadrature_coeffs);
FluidMoments moments(integrate);

moments(density_computed.span_view(), allfdistribu.span_cview(), FluidMoments::s_density);
moments(mean_velocity_computed.span_view(),
Expand Down
6 changes: 4 additions & 2 deletions tests/geometryXVx/kineticsource.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -50,8 +50,10 @@ TEST(KineticSource, Moments)
IDomainVx const gridvx = builder_vx.interpolation_domain();
IDomainSpXVx const mesh(IDomainSp(my_iion, IVectSp(1)), gridx, gridvx);

Quadrature<IDimX> const integrate_x(trapezoid_quadrature_coefficients(gridx));
Quadrature<IDimVx> const integrate_v(trapezoid_quadrature_coefficients(gridvx));
DFieldX quadrature_coeffs_x = trapezoid_quadrature_coefficients(gridx);
DFieldVx quadrature_coeffs_vx = trapezoid_quadrature_coefficients(gridvx);
Quadrature<IDimX> const integrate_x(quadrature_coeffs_x);
Quadrature<IDimVx> const integrate_v(quadrature_coeffs_vx);

FieldSp<int> charges(dom_sp);
charges(my_ielec) = -1;
Expand Down
6 changes: 4 additions & 2 deletions tests/geometryXVx/krooksource.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -60,8 +60,10 @@ TEST(KrookSource, Adaptive)
IDomainSp const gridsp = dom_sp;
IDomainSpXVx const mesh(gridsp, gridx, gridvx);

Quadrature<IDimX> const integrate_x(trapezoid_quadrature_coefficients(gridx));
Quadrature<IDimVx> const integrate_v(trapezoid_quadrature_coefficients(gridvx));
DFieldX const quadrature_coeffs_x = trapezoid_quadrature_coefficients(gridx);
DFieldVx const quadrature_coeffs_vx = trapezoid_quadrature_coefficients(gridvx);
Quadrature<IDimX> const integrate_x(quadrature_coeffs_x);
Quadrature<IDimVx> const integrate_v(quadrature_coeffs_vx);

FieldSp<int> charges(dom_sp);
DFieldSp masses(dom_sp);
Expand Down
Loading

0 comments on commit 20668da

Please sign in to comment.