Skip to content

Commit

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

See merge request gysela-developpers/gyselalibxx!607
  • Loading branch information
EmilyBourne committed Jul 29, 2024
1 parent 9d25977 commit 2674b9f
Show file tree
Hide file tree
Showing 7 changed files with 225 additions and 235 deletions.
14 changes: 8 additions & 6 deletions src/advection/bsl_advection_1d.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -96,13 +96,15 @@ class BslAdvection1D


// Interpolators:
using FunctionPreallocatableInterpolatorType
= interpolator_on_domain_t<IPreallocatableInterpolator, IDimInterest, FunctionDomain>;
using FunctionPreallocatableInterpolatorType = interpolator_on_idx_range_t<
IPreallocatableInterpolator,
IDimInterest,
FunctionDomain>;
using FunctionInterpolatorType
= interpolator_on_domain_t<IInterpolator, IDimInterest, FunctionDomain>;
= interpolator_on_idx_range_t<IInterpolator, IDimInterest, FunctionDomain>;

// Type for the derivatives of the function
using FunctionDerivDomain = typename FunctionInterpolatorType::batched_derivs_domain_type;
using FunctionDerivDomain = typename FunctionInterpolatorType::batched_derivs_idx_range_type;
using FunctionDerivChunk = device_t<ddc::Chunk<double, FunctionDerivDomain>>;

FunctionPreallocatableInterpolatorType const& m_function_interpolator;
Expand Down Expand Up @@ -191,9 +193,9 @@ class BslAdvection1D

// Build derivatives on boundaries and fill with zeros....................................
FunctionDerivChunk function_derivatives_min(
m_function_interpolator.batched_derivs_domain_xmin(function_dom));
m_function_interpolator.batched_derivs_idx_range_xmin(function_dom));
FunctionDerivChunk function_derivatives_max(
m_function_interpolator.batched_derivs_domain_xmax(function_dom));
m_function_interpolator.batched_derivs_idx_range_xmax(function_dom));
ddc::parallel_fill(Kokkos::DefaultExecutionSpace(), function_derivatives_min, 0.);
ddc::parallel_fill(Kokkos::DefaultExecutionSpace(), function_derivatives_max, 0.);

Expand Down
12 changes: 6 additions & 6 deletions src/advection/bsl_advection_vx.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -23,11 +23,11 @@ class BslAdvectionVelocity : public IAdvectionVelocity<Geometry, DDimV>
using CDimV = typename DDimV::continuous_dimension_type;

private:
using PreallocatableInterpolatorType = interpolator_on_domain_t<
using PreallocatableInterpolatorType = interpolator_on_idx_range_t<
IPreallocatableInterpolator,
DDimV,
ddc::cartesian_prod_t<typename Geometry::SpatialDDom, typename Geometry::VelocityDDom>>;
using InterpolatorType = interpolator_on_domain_t<
using InterpolatorType = interpolator_on_idx_range_t<
IInterpolator,
DDimV,
ddc::cartesian_prod_t<typename Geometry::SpatialDDom, typename Geometry::VelocityDDom>>;
Expand Down Expand Up @@ -62,11 +62,11 @@ class BslAdvectionVelocity : public IAdvectionVelocity<Geometry, DDimV>
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_domain_type>>
derivs_min(m_interpolator_v.batched_derivs_domain_xmin(
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_domain_type>>
derivs_max(m_interpolator_v.batched_derivs_domain_xmax(
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)));
ddc::parallel_fill(derivs_min, 0.);
ddc::parallel_fill(derivs_max, 0.);
Expand Down
4 changes: 2 additions & 2 deletions src/advection/bsl_advection_x.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -25,11 +25,11 @@ class BslAdvectionSpatial : public IAdvectionSpatial<Geometry, DDimX>
using CDimV = typename DDimV::continuous_dimension_type;

private:
using PreallocatableInterpolatorType = interpolator_on_domain_t<
using PreallocatableInterpolatorType = interpolator_on_idx_range_t<
IPreallocatableInterpolator,
DDimX,
ddc::cartesian_prod_t<typename Geometry::SpatialDDom, typename Geometry::VelocityDDom>>;
using InterpolatorType = interpolator_on_domain_t<
using InterpolatorType = interpolator_on_idx_range_t<
IInterpolator,
DDimX,
ddc::cartesian_prod_t<typename Geometry::SpatialDDom, typename Geometry::VelocityDDom>>;
Expand Down
133 changes: 66 additions & 67 deletions src/interpolation/Lagrange.hpp
Original file line number Diff line number Diff line change
@@ -1,12 +1,15 @@
// SPDX-License-Identifier: MIT
#pragma once

#include <iostream>
#include <vector>

#include <ddc/ddc.hpp>

#include <ddc_helper.hpp>

#include "ddc_aliases.hpp"

enum class BCond { PERIODIC, DIRICHLET };

/**
Expand All @@ -15,51 +18,47 @@ enum class BCond { PERIODIC, DIRICHLET };
* A simple class which provides the possibility to evaluate
* an interpolation of a function known only over a restricted set of nodes.
*/
template <class Execspace, class DDimI, BCond BcMin, BCond BcMax>
template <class Execspace, class GridInterp, BCond BcMin, BCond BcMax>
class Lagrange
{
static_assert((BcMin == BCond::PERIODIC) == (BcMax == BCond::PERIODIC));
using DElemI = ddc::DiscreteElement<DDimI>;
using DVectI = ddc::DiscreteVector<DDimI>;
using IdxInterp = Idx<GridInterp>;
using IdxStepInterp = IdxStep<GridInterp>;
using IdxRangeInterp = IdxRange<GridInterp>;

using CoordDimI = ddc::Coordinate<typename DDimI::continuous_dimension_type>;
using CoordDimI = Coord<typename GridInterp::continuous_dimension_type>;

private:
ddc::DiscreteDomain<DDimI> m_domain;
ddc::DiscreteDomain<DDimI> m_inner_domain;
ddc::ChunkSpan<
double,
ddc::DiscreteDomain<DDimI>,
std::experimental::layout_right,
typename Execspace::memory_space>
IdxRangeInterp m_idx_range;
IdxRangeInterp m_inner_idx_range;
Field<double, IdxRangeInterp, std::experimental::layout_right, typename Execspace::memory_space>
m_ChunkSpan;
CoordDimI m_left_bound;
CoordDimI m_right_bound;
DVectI m_poly_support;
IdxStepInterp m_poly_support;

public:
/**
* @brief Usual Constructor
*
* @param[in] degree integer which correspond to the degree of interpolation.
* @param[in] x_nodes_fnodes Chunkspan of nodes and associated values of the function.
* @param[in] domain along interest direction, usedful in periodic case
* @param[in] idx_range along interest direction, usedful in periodic case
* @param[in] ghost DiscretVector which gives the number of ghosted points
*/
KOKKOS_FUNCTION Lagrange(
int degree,
ddc::ChunkSpan<
double,
ddc::DiscreteDomain<DDimI>,
std::experimental::layout_right,
typename Execspace::memory_space> x_nodes_fnodes,
ddc::DiscreteDomain<DDimI> domain,
DVectI ghost)
: m_domain(domain)
, m_inner_domain(domain.remove(ghost, ghost))
Field<double,
IdxRangeInterp,
std::experimental::layout_right,
typename Execspace::memory_space> x_nodes_fnodes,
IdxRangeInterp idx_range,
IdxStepInterp ghost)
: m_idx_range(idx_range)
, m_inner_idx_range(idx_range.remove(ghost, ghost))
, m_ChunkSpan(x_nodes_fnodes)
, m_left_bound(ddc::coordinate(m_inner_domain.front()))
, m_right_bound(ddc::coordinate(m_inner_domain.back()))
, m_left_bound(ddc::coordinate(m_inner_idx_range.front()))
, m_right_bound(ddc::coordinate(m_inner_idx_range.back()))
, m_poly_support(degree + 1)
{
}
Expand All @@ -75,40 +74,40 @@ class Lagrange
KOKKOS_FUNCTION double evaluate(CoordDimI x_interp) const;

private:
DElemI getclosest(CoordDimI value) const;
IdxInterp getclosest(CoordDimI value) const;

KOKKOS_FUNCTION DElemI getclosest_binsearch(CoordDimI value) const;
KOKKOS_FUNCTION IdxInterp getclosest_binsearch(CoordDimI value) const;

KOKKOS_FUNCTION double evaluate_lagrange(CoordDimI x_interp) const;

KOKKOS_FUNCTION double apply_bc(CoordDimI x_interp) const;

KOKKOS_FUNCTION double compute_basis(
CoordDimI x_interp,
DElemI j,
ddc::DiscreteDomain<DDimI> polynom_subdomain) const;
IdxInterp j,
IdxRangeInterp polynom_subidx_range) const;
};

template <typename Execspace, class DDimI, BCond BcMin, BCond BcMax>
ddc::DiscreteElement<DDimI> Lagrange<Execspace, DDimI, BcMin, BcMax>::getclosest(
CoordDimI value) const
template <typename Execspace, class GridInterp, BCond BcMin, BCond BcMax>
Idx<GridInterp> Lagrange<Execspace, GridInterp, BcMin, BcMax>::getclosest(CoordDimI value) const
{
assert(value >= m_left_bound && value <= m_right_bound);
auto it = std::min_element(m_domain.begin(), m_domain.end(), [=](DElemI x, DElemI y) {
return std::abs(ddc::coordinate(x) - value) < std::abs(ddc::coordinate(y) - value);
});
auto it = std::
min_element(m_idx_range.begin(), m_idx_range.end(), [=](IdxInterp x, IdxInterp y) {
return std::abs(ddc::coordinate(x) - value) < std::abs(ddc::coordinate(y) - value);
});
return *it;
}
template <typename Execspace, class DDimI, BCond BcMin, BCond BcMax>
KOKKOS_INLINE_FUNCTION ddc::DiscreteElement<DDimI> Lagrange<Execspace, DDimI, BcMin, BcMax>::
template <typename Execspace, class GridInterp, BCond BcMin, BCond BcMax>
KOKKOS_INLINE_FUNCTION Idx<GridInterp> Lagrange<Execspace, GridInterp, BcMin, BcMax>::
getclosest_binsearch(CoordDimI x_interp) const
{
assert(x_interp >= m_left_bound && x_interp <= m_right_bound);
DElemI begin = m_domain.front();
DElemI end = m_domain.back();
DElemI elm_cell = begin + (end - begin) / 2;
IdxInterp begin = m_idx_range.front();
IdxInterp end = m_idx_range.back();
IdxInterp elm_cell = begin + (end - begin) / 2;
while (x_interp < ddc::coordinate(elm_cell)
|| x_interp > ddc::coordinate(elm_cell + DVectI(1))) {
|| x_interp > ddc::coordinate(elm_cell + IdxStepInterp(1))) {
if (x_interp < ddc::coordinate(elm_cell)) {
end = elm_cell;
} else {
Expand All @@ -125,19 +124,19 @@ KOKKOS_INLINE_FUNCTION ddc::DiscreteElement<DDimI> Lagrange<Execspace, DDimI, Bc
*
* @param[in] x_interp a node where we want to evaluate the function.
* @param[in] j index of the basis.
* @param[in] polynom_subdomain a part of the mesh centered around x_interp
* @param[in] polynom_subindex range a part of the mesh centered around x_interp
*
* @return The value of the basis at x_intercept
*/
template <typename Execspace, class DDimI, BCond BcMin, BCond BcMax>
KOKKOS_INLINE_FUNCTION double Lagrange<Execspace, DDimI, BcMin, BcMax>::compute_basis(
template <typename Execspace, class GridInterp, BCond BcMin, BCond BcMax>
KOKKOS_INLINE_FUNCTION double Lagrange<Execspace, GridInterp, BcMin, BcMax>::compute_basis(
CoordDimI x_interp,
DElemI j,
ddc::DiscreteDomain<DDimI> polynom_subdomain) const
IdxInterp j,
IdxRangeInterp polynom_subidx_range) const
{
double w(1);
CoordDimI coord_j = ddc::coordinate(j);
for (DElemI const i : polynom_subdomain) {
for (IdxInterp const i : polynom_subidx_range) {
CoordDimI coord_i = ddc::coordinate(i);
if (coord_i != coord_j) {
w *= (x_interp - coord_i) / (coord_j - coord_i);
Expand All @@ -146,8 +145,8 @@ KOKKOS_INLINE_FUNCTION double Lagrange<Execspace, DDimI, BcMin, BcMax>::compute_
return w;
}

template <typename Execspace, class DDimI, BCond BcMin, BCond BcMax>
KOKKOS_INLINE_FUNCTION double Lagrange<Execspace, DDimI, BcMin, BcMax>::apply_bc(
template <typename Execspace, class GridInterp, BCond BcMin, BCond BcMax>
KOKKOS_INLINE_FUNCTION double Lagrange<Execspace, GridInterp, BcMin, BcMax>::apply_bc(
CoordDimI x_interp) const
{
CoordDimI bc_val = x_interp;
Expand All @@ -164,8 +163,8 @@ KOKKOS_INLINE_FUNCTION double Lagrange<Execspace, DDimI, BcMin, BcMax>::apply_bc
return evaluate_lagrange(bc_val);
}

template <typename Execspace, class DDimI, BCond BcMin, BCond BcMax>
KOKKOS_INLINE_FUNCTION double Lagrange<Execspace, DDimI, BcMin, BcMax>::evaluate(
template <typename Execspace, class GridInterp, BCond BcMin, BCond BcMax>
KOKKOS_INLINE_FUNCTION double Lagrange<Execspace, GridInterp, BcMin, BcMax>::evaluate(
CoordDimI x_interp) const
{
if (x_interp < m_left_bound || m_right_bound < x_interp) {
Expand All @@ -175,37 +174,37 @@ KOKKOS_INLINE_FUNCTION double Lagrange<Execspace, DDimI, BcMin, BcMax>::evaluate
}
}

template <typename Execspace, class DDimI, BCond BcMin, BCond BcMax>
KOKKOS_INLINE_FUNCTION double Lagrange<Execspace, DDimI, BcMin, BcMax>::evaluate_lagrange(
template <typename Execspace, class GridInterp, BCond BcMin, BCond BcMax>
KOKKOS_INLINE_FUNCTION double Lagrange<Execspace, GridInterp, BcMin, BcMax>::evaluate_lagrange(
CoordDimI x_intern) const
{
assert(x_intern >= m_left_bound && x_intern <= m_right_bound);

DElemI begin, end;
DElemI icell = getclosest_binsearch(x_intern);
DElemI mid = icell;
if (mid >= m_inner_domain.back() && BcMax == BCond::PERIODIC) {
IdxInterp begin, end;
IdxInterp icell = getclosest_binsearch(x_intern);
IdxInterp mid = icell;
if (mid >= m_inner_idx_range.back() && BcMax == BCond::PERIODIC) {
begin = mid - m_poly_support / 2;
end = std::min(m_inner_domain.back(), begin + m_poly_support);
} else if (mid <= m_inner_domain.front() && BcMin == BCond::PERIODIC) {
begin = std::max(m_domain.front(), mid - m_poly_support / 2);
end = std::min(m_inner_idx_range.back(), begin + m_poly_support);
} else if (mid <= m_inner_idx_range.front() && BcMin == BCond::PERIODIC) {
begin = std::max(m_idx_range.front(), mid - m_poly_support / 2);
end = begin + m_poly_support;
} else {
if (m_inner_domain.front() + m_poly_support / 2 > mid) {
begin = m_inner_domain.front();
if (m_inner_idx_range.front() + m_poly_support / 2 > mid) {
begin = m_inner_idx_range.front();
} else {
begin = mid - m_poly_support / 2;
}
end = std::min(m_domain.back() + DVectI(1), begin + m_poly_support);
end = std::min(m_idx_range.back() + IdxStepInterp(1), begin + m_poly_support);
}

if (end == m_domain.back())
if (end == m_idx_range.back())
begin = end - m_poly_support;
DVectI npoints_subdomain(end - begin);
ddc::DiscreteDomain<DDimI> subdomain(begin, npoints_subdomain);
IdxStepInterp npoints_subidx_range(end - begin);
IdxRangeInterp subidx_range(begin, npoints_subidx_range);
double p = 0;
for (DElemI const ix : subdomain) {
p += compute_basis(x_intern, ix, subdomain) * m_ChunkSpan(ix);
for (IdxInterp const ix : subidx_range) {
p += compute_basis(x_intern, ix, subidx_range) * m_ChunkSpan(ix);
}
return p;
}
Loading

0 comments on commit 2674b9f

Please sign in to comment.