Skip to content

Commit

Permalink
Merge branch 'renaming_advection' into 'main'
Browse files Browse the repository at this point in the history
Update naming conventions in advection

See merge request gysela-developpers/gyselalibxx!614
  • Loading branch information
EmilyBourne committed Jul 30, 2024
1 parent 2674b9f commit dd78405
Show file tree
Hide file tree
Showing 12 changed files with 586 additions and 570 deletions.
4 changes: 2 additions & 2 deletions bin/update_naming_conventions.py
Original file line number Diff line number Diff line change
Expand Up @@ -271,11 +271,11 @@ def get_geometry_name_shortcuts(geometry_file):
struct_name = s[2]
if struct_name in species_names:
continue
if struct_name == 'RDim':
if struct_name in ('RDim', 'CDim'):
name_map[struct_name] = 'Dim'
elif struct_name in ('DDim', 'IDim'):
name_map[struct_name] = 'Grid1D'
elif struct_name.startswith('RDim'):
elif struct_name.startswith('RDim') or struct_name.startswith('CDim'):
new_name = struct_name[4:]
if not new_name[0].isalpha():
new_name = 'Dim' + new_name
Expand Down
169 changes: 85 additions & 84 deletions src/advection/bsl_advection_1d.hpp

Large diffs are not rendered by default.

79 changes: 38 additions & 41 deletions src/advection/bsl_advection_vx.hpp
Original file line number Diff line number Diff line change
@@ -1,42 +1,41 @@
// SPDX-License-Identifier: MIT

#pragma once

#include <ddc/ddc.hpp>

#include <ddc_helper.hpp>
#include <iinterpolator.hpp>
#include <species_info.hpp>

#include "ddc_aliases.hpp"
#include "iadvectionvx.hpp"

/**
* @brief A class which computes the velocity advection along the dimension of interest DDimV. Working for every cartesian geometry.
* @brief A class which computes the velocity advection along the dimension of interest GridV. Working for every cartesian geometry.
*/
template <class Geometry, class DDimV>
class BslAdvectionVelocity : public IAdvectionVelocity<Geometry, DDimV>
template <class Geometry, class GridV>
class BslAdvectionVelocity : public IAdvectionVelocity<Geometry, GridV>
{
using FdistribuDDom = typename Geometry::FdistribuDDom;
using SpatialDDom = typename Geometry::SpatialDDom;
using DElemV = ddc::DiscreteElement<DDimV>;
using DElemSp = ddc::DiscreteElement<Species>;
using CDimV = typename DDimV::continuous_dimension_type;
using FdistribuIdxRange = typename Geometry::FdistribuDDom;
using SpatialIdxRange = typename Geometry::SpatialDDom;
using IdxV = Idx<GridV>;
using DimV = typename GridV::continuous_dimension_type;

private:
using PreallocatableInterpolatorType = interpolator_on_idx_range_t<
IPreallocatableInterpolator,
DDimV,
GridV,
ddc::cartesian_prod_t<typename Geometry::SpatialDDom, typename Geometry::VelocityDDom>>;
using InterpolatorType = interpolator_on_idx_range_t<
IInterpolator,
DDimV,
GridV,
ddc::cartesian_prod_t<typename Geometry::SpatialDDom, typename Geometry::VelocityDDom>>;
PreallocatableInterpolatorType const& m_interpolator_v;

public:
/**
* @brief Constructor
* @param[in] interpolator_v interpolator along the DDimV direction which refers to the velocity space.
* @param[in] interpolator_v interpolator along the GridV direction which refers to the velocity space.
*/
explicit BslAdvectionVelocity(PreallocatableInterpolatorType const& interpolator_v)
: m_interpolator_v(interpolator_v)
Expand All @@ -46,68 +45,66 @@ class BslAdvectionVelocity : public IAdvectionVelocity<Geometry, DDimV>
~BslAdvectionVelocity() override = default;

/**
* @brief Advects fdistribu along DDimV for a duration dt.
* @brief Advects fdistribu along GridV for a duration dt.
* @param[in, out] allfdistribu Reference to the whole distribution function for one species, allocated on the device (ie it lets the choice of the location depend on the build configuration).
* @param[in] electric_field Reference to the electric field which derives from electrostatic potential, allocated on the device.
* @param[in] dt Time step
* @return A reference to the allfdistribu array containing the value of the function at the coordinates.
*/
device_t<ddc::ChunkSpan<double, FdistribuDDom>> operator()(
device_t<ddc::ChunkSpan<double, FdistribuDDom>> const allfdistribu,
device_t<ddc::ChunkSpan<const double, SpatialDDom>> const electric_field,
Field<double, FdistribuIdxRange> operator()(
Field<double, FdistribuIdxRange> const allfdistribu,
Field<const double, SpatialIdxRange> const electric_field,
double const dt) const override
{
Kokkos::Profiling::pushRegion("BslAdvectionVelocity");
FdistribuDDom const dom = allfdistribu.domain();
ddc::DiscreteDomain<DDimV> const v_dom = ddc::select<DDimV>(dom);
ddc::DiscreteDomain<Species> const sp_dom = ddc::select<Species>(dom);

device_t<ddc::Chunk<double, typename InterpolatorType::batched_derivs_idx_range_type>>
derivs_min(m_interpolator_v.batched_derivs_idx_range_xmin(
ddc::remove_dims_of<Species>(dom)));
device_t<ddc::Chunk<double, typename InterpolatorType::batched_derivs_idx_range_type>>
derivs_max(m_interpolator_v.batched_derivs_idx_range_xmax(
ddc::remove_dims_of<Species>(dom)));
FdistribuIdxRange const dom = get_idx_range(allfdistribu);
IdxRange<GridV> const v_dom = ddc::select<GridV>(dom);
IdxRange<Species> const sp_dom = ddc::select<Species>(dom);

FieldMem<double, typename InterpolatorType::batched_derivs_idx_range_type> derivs_min(
m_interpolator_v.batched_derivs_idx_range_xmin(ddc::remove_dims_of<Species>(dom)));
FieldMem<double, typename InterpolatorType::batched_derivs_idx_range_type> derivs_max(
m_interpolator_v.batched_derivs_idx_range_xmax(ddc::remove_dims_of<Species>(dom)));
ddc::parallel_fill(derivs_min, 0.);
ddc::parallel_fill(derivs_max, 0.);

// pre-allocate some memory to prevent allocation later in loop
ddc::Chunk feet_coords_alloc(
ddc::remove_dims_of(dom, sp_dom),
ddc::DeviceAllocator<ddc::Coordinate<CDimV>>());
ddc::ChunkSpan feet_coords = feet_coords_alloc.span_view();
ddc::DeviceAllocator<Coord<DimV>>());
ddc::ChunkSpan feet_coords = get_field(feet_coords_alloc);
std::unique_ptr<InterpolatorType> const interpolator_v_ptr = m_interpolator_v.preallocate();
InterpolatorType const& interpolator_v = *interpolator_v_ptr;

SpatialDDom const spatial_dom(allfdistribu.domain());
SpatialIdxRange const spatial_dom(get_idx_range(allfdistribu));

auto c_dom = ddc::remove_dims_of(ddc::remove_dims_of(dom, sp_dom), v_dom);
using DElemC = typename decltype(c_dom)::discrete_element_type;
using DElemSpatial = typename SpatialDDom::discrete_element_type;
auto batch_idx_range = ddc::remove_dims_of(ddc::remove_dims_of(dom, sp_dom), v_dom);
using IdxB = typename decltype(batch_idx_range)::discrete_element_type;
using IdxSpatial = typename SpatialIdxRange::discrete_element_type;

ddc::for_each(sp_dom, [&](DElemSp const isp) {
ddc::for_each(sp_dom, [&](IdxSp const isp) {
double const charge_proxy
= charge(isp); // TODO: consider proper way to access charge from device
double const sqrt_me_on_mspecies = std::sqrt(mass(ielec()) / mass(isp));
ddc::parallel_for_each(
Kokkos::DefaultExecutionSpace(),
c_dom,
KOKKOS_LAMBDA(DElemC const ic) {
DElemSpatial const ix(ic);
batch_idx_range,
KOKKOS_LAMBDA(IdxB const ib) {
IdxSpatial const ix(ib);
// compute the displacement
double const dvx
= charge_proxy * sqrt_me_on_mspecies * dt * electric_field(ix);

// compute the coordinates of the feet
for (DElemV const iv : v_dom) {
feet_coords(iv, ic) = ddc::Coordinate<CDimV>(ddc::coordinate(iv) - dvx);
for (IdxV const iv : v_dom) {
feet_coords(iv, ib) = Coord<DimV>(ddc::coordinate(iv) - dvx);
}
});
interpolator_v(
allfdistribu[isp],
feet_coords.span_cview(),
derivs_min.span_cview(),
derivs_max.span_cview());
get_const_field(feet_coords),
get_const_field(derivs_min),
get_const_field(derivs_max));
});

Kokkos::Profiling::popRegion();
Expand Down
69 changes: 34 additions & 35 deletions src/advection/bsl_advection_x.hpp
Original file line number Diff line number Diff line change
@@ -1,44 +1,42 @@
// SPDX-License-Identifier: MIT

#pragma once

#include <ddc/ddc.hpp>

#include <iinterpolator.hpp>
#include <species_info.hpp>

#include "ddc_aliases.hpp"
#include "iadvectionx.hpp"

/**
* @brief A class which computes the spatial advection along the dimension of interest DDimX. Working for every cartesian geometry.
* @brief A class which computes the spatial advection along the dimension of interest GridX. Working for every cartesian geometry.
*/
template <class Geometry, class DDimX>
class BslAdvectionSpatial : public IAdvectionSpatial<Geometry, DDimX>
template <class Geometry, class GridX>
class BslAdvectionSpatial : public IAdvectionSpatial<Geometry, GridX>
{
using DDimV = typename Geometry::template velocity_dim_for<DDimX>;
using DDom = typename Geometry::FdistribuDDom;
using DElemX = ddc::DiscreteElement<DDimX>;
using DElemV = ddc::DiscreteElement<DDimV>;
using DElemSp = ddc::DiscreteElement<Species>;
using DElemSpV = ddc::DiscreteElement<Species, DDimV>;
using CDimX = typename DDimX::continuous_dimension_type;
using CDimV = typename DDimV::continuous_dimension_type;
using GridV = typename Geometry::template velocity_dim_for<GridX>;
using FdistribIdxRange = typename Geometry::FdistribuDDom;
using IdxX = Idx<GridX>;
using IdxV = Idx<GridV>;
using DimX = typename GridX::continuous_dimension_type;
using DimV = typename GridV::continuous_dimension_type;

private:
using PreallocatableInterpolatorType = interpolator_on_idx_range_t<
IPreallocatableInterpolator,
DDimX,
GridX,
ddc::cartesian_prod_t<typename Geometry::SpatialDDom, typename Geometry::VelocityDDom>>;
using InterpolatorType = interpolator_on_idx_range_t<
IInterpolator,
DDimX,
GridX,
ddc::cartesian_prod_t<typename Geometry::SpatialDDom, typename Geometry::VelocityDDom>>;
PreallocatableInterpolatorType const& m_interpolator_x;

public:
/**
* @brief Constructor
* @param[in] interpolator_x interpolator along the DDimX direction which refers to the spatial space.
* @param[in] interpolator_x interpolator along the GridX direction which refers to the spatial space.
*/
explicit BslAdvectionSpatial(PreallocatableInterpolatorType const& interpolator_x)
: m_interpolator_x(interpolator_x)
Expand All @@ -48,49 +46,50 @@ class BslAdvectionSpatial : public IAdvectionSpatial<Geometry, DDimX>
~BslAdvectionSpatial() override = default;

/**
* @brief Advects fdistribu along DDimX for a duration dt.
* @brief Advects fdistribu along GridX for a duration dt.
* @param[in, out] allfdistribu Reference to the whole distribution function for one species, allocated on the device (ie it lets the choice of the location depend on the build configuration).
* @param[in] dt Time step
* @return A reference to the allfdistribu array containing the value of the function at the coordinates.
*/
device_t<ddc::ChunkSpan<double, DDom>> operator()(
device_t<ddc::ChunkSpan<double, DDom>> const allfdistribu,
Field<double, FdistribIdxRange> operator()(
Field<double, FdistribIdxRange> const allfdistribu,
double const dt) const override
{
Kokkos::Profiling::pushRegion("BslAdvectionSpatial");
DDom const dom = allfdistribu.domain();
ddc::DiscreteDomain<DDimX> const x_dom = ddc::select<DDimX>(dom);
ddc::DiscreteDomain<DDimV> const v_dom = ddc::select<DDimV>(dom);
ddc::DiscreteDomain<Species> const sp_dom = ddc::select<Species>(dom);
FdistribIdxRange const dom = get_idx_range(allfdistribu);
IdxRange<GridX> const x_idx_range = ddc::select<GridX>(dom);
IdxRange<GridV> const v_idx_range = ddc::select<GridV>(dom);
IdxRange<Species> const sp_idx_range = ddc::select<Species>(dom);

// pre-allocate some memory to prevent allocation later in loop
ddc::Chunk feet_coords_alloc(
ddc::remove_dims_of(dom, sp_dom),
ddc::DeviceAllocator<ddc::Coordinate<CDimX>>());
ddc::ChunkSpan feet_coords = feet_coords_alloc.span_view();
ddc::remove_dims_of(dom, sp_idx_range),
ddc::DeviceAllocator<Coord<DimX>>());
ddc::ChunkSpan feet_coords = get_field(feet_coords_alloc);
std::unique_ptr<InterpolatorType> const interpolator_x_ptr = m_interpolator_x.preallocate();
InterpolatorType const& interpolator_x = *interpolator_x_ptr;

auto c_dom = ddc::remove_dims_of(dom, ddc::DiscreteDomain<Species, DDimX>(sp_dom, x_dom));
using DElemC = typename decltype(c_dom)::discrete_element_type;
auto batch_idx_range
= ddc::remove_dims_of(dom, IdxRange<Species, GridX>(sp_idx_range, x_idx_range));
using IdxB = typename decltype(batch_idx_range)::discrete_element_type;

for (DElemSp const isp : sp_dom) {
for (IdxSp const isp : sp_idx_range) {
double const sqrt_me_on_mspecies = std::sqrt(mass(ielec()) / mass(isp));
ddc::parallel_for_each(
Kokkos::DefaultExecutionSpace(),
c_dom,
KOKKOS_LAMBDA(DElemC const ic) {
batch_idx_range,
KOKKOS_LAMBDA(IdxB const ib) {
// compute the displacement
DElemV const iv(ic);
ddc::Coordinate<CDimV> const coord_iv = ddc::coordinate(iv);
IdxV const iv(ib);
Coord<DimV> const coord_iv = ddc::coordinate(iv);
double const dx = sqrt_me_on_mspecies * dt * coord_iv;

// compute the coordinates of the feet
for (DElemX const ix : x_dom) {
feet_coords(ix, ic) = ddc::Coordinate<CDimX>(ddc::coordinate(ix) - dx);
for (IdxX const ix : x_idx_range) {
feet_coords(ix, ib) = Coord<DimX>(ddc::coordinate(ix) - dx);
}
});
interpolator_x(allfdistribu[isp], feet_coords.span_cview());
interpolator_x(allfdistribu[isp], get_const_field(feet_coords));
}

Kokkos::Profiling::popRegion();
Expand Down
14 changes: 7 additions & 7 deletions src/advection/iadvectionvx.hpp
Original file line number Diff line number Diff line change
@@ -1,18 +1,19 @@
// SPDX-License-Identifier: MIT

#pragma once

#include <ddc/ddc.hpp>

#include <ddc_helper.hpp>

#include "ddc_aliases.hpp"

/**
* @brief A class which provides an advection operator.
*
* An abstract class which implements a function that
* applies the transport along a velocity direction of the phase space.
*/
template <class Geometry, class DDimV>
template <class Geometry, class GridV>
class IAdvectionVelocity
{
public:
Expand All @@ -22,14 +23,13 @@ class IAdvectionVelocity
* @brief operates a transport of the distribution function.
*
* @param[in, out] allfdistribu Reference to an array containing the value of the distribution function.
* @param[in] electrostatic_potential The electrostatic potential which is the advection speed.
* @param[in] electric_field The electric field which derives from electrostatic potential and is the advection speed.
* @param[in] dt Time step.
*
* @return A reference to an array containing the value of distribution the function at the updated time t+dt.
*/
virtual device_t<ddc::ChunkSpan<double, typename Geometry::FdistribuDDom>> operator()(
device_t<ddc::ChunkSpan<double, typename Geometry::FdistribuDDom>> allfdistribu,
device_t<ddc::ChunkSpan<const double, typename Geometry::SpatialDDom>>
electrostatic_potential,
virtual Field<double, typename Geometry::FdistribuDDom> operator()(
Field<double, typename Geometry::FdistribuDDom> allfdistribu,
Field<const double, typename Geometry::SpatialDDom> electric_field,
double dt) const = 0;
};
9 changes: 5 additions & 4 deletions src/advection/iadvectionx.hpp
Original file line number Diff line number Diff line change
@@ -1,16 +1,17 @@
// SPDX-License-Identifier: MIT

#pragma once

#include <ddc/ddc.hpp>

#include "ddc_aliases.hpp"

/**
* @brief A class which provides an advection operator.
*
* An abstract class which implements a function that
* applies the transport along a physical space direction of the phase space.
*/
template <class Geometry, class DDimX>
template <class Geometry, class GridX>
class IAdvectionSpatial
{
public:
Expand All @@ -23,7 +24,7 @@ class IAdvectionSpatial
*
* @return A reference to an array containing the value of distribution the function at the updated time T+dt.
*/
virtual device_t<ddc::ChunkSpan<double, typename Geometry::FdistribuDDom>> operator()(
device_t<ddc::ChunkSpan<double, typename Geometry::FdistribuDDom>> allfdistribu,
virtual Field<double, typename Geometry::FdistribuDDom> operator()(
Field<double, typename Geometry::FdistribuDDom> allfdistribu,
double dt) const = 0;
};
Loading

0 comments on commit dd78405

Please sign in to comment.