Skip to content

Commit

Permalink
Add is_2d_mapping_v method
Browse files Browse the repository at this point in the history
Add is_2d_mapping_v method which checks that a mapping is fully defined (all jacobian methods and the operator are present). Fixes #408

See merge request gysela-developpers/gyselalibxx!780

--------------------------------------------
  • Loading branch information
EmilyBourne committed Nov 22, 2024
1 parent d4245b0 commit dffd882
Show file tree
Hide file tree
Showing 15 changed files with 719 additions and 229 deletions.
2 changes: 1 addition & 1 deletion docs/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ set(DOXYGEN_DISABLE_INDEX YES)
set(DOXYGEN_ENUM_VALUES_PER_LINE 1)
set(DOXYGEN_EXPAND_ONLY_PREDEF NO)
set(DOXYGEN_EXTRACT_ALL NO)
set(DOXYGEN_EXCLUDE_SYMBOLS "detail" "details")
set(DOXYGEN_EXCLUDE_SYMBOLS "detail" "details" "mapping_detail")
set(DOXYGEN_EXTRACT_LOCAL_CLASSES YES)
set(DOXYGEN_FULL_PATH_NAMES YES)
set(DOXYGEN_GENERATE_TREEVIEW YES)
Expand Down
5 changes: 2 additions & 3 deletions src/mapping/vector_mapper.hpp
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@
// SPDX-License-Identifier: MIT
#pragma once

#include <sll/view.hpp>

#include "directional_tag.hpp"
#include "vector_field.hpp"
#include "vector_field_mem.hpp"
Expand Down Expand Up @@ -35,9 +37,6 @@ class VectorMapper<NDTag<XIn, YIn>, NDTag<XOut, YOut>, Mapping, ExecSpace>
/// The vector type in the coordinate system returned as output.
using vector_element_type_out = typename ddc::detail::TaggedVector<double, XOut, YOut>;

private:
using Matrix_2x2 = typename Mapping::Matrix_2x2;

private:
Mapping m_mapping;

Expand Down
2 changes: 2 additions & 0 deletions src/pde_solvers/polarpoissonlikesolver.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
#include <ddc/ddc.hpp>

#include <sll/gauss_legendre_integration.hpp>
#include <sll/mapping/mapping_tools.hpp>
#include <sll/mapping/metric_tensor.hpp>
#include <sll/math_tools.hpp>
#include <sll/matrix_batch_csr.hpp>
Expand Down Expand Up @@ -305,6 +306,7 @@ class PolarSplineFEMPoissonLikeSolver
IdxRangeQuadratureRTheta(m_idxrange_quadrature_r, m_idxrange_quadrature_theta))
, m_polar_spline_evaluator(ddc::NullExtrapolationRule())
{
static_assert(has_2d_jacobian_v<Mapping, CoordRTheta>);
// Get break points
IdxRange<KnotsR> idxrange_r_edges = ddc::discrete_space<BSplinesR>().break_point_domain();
IdxRange<KnotsTheta> idxrange_theta_edges
Expand Down
3 changes: 3 additions & 0 deletions src/quadrature/volume_quadrature_nd.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@
#pragma once
#include <ddc/ddc.hpp>

#include <sll/mapping/mapping_tools.hpp>

#include "ddc_alias_inline_functions.hpp"
#include "ddc_aliases.hpp"
#include "quadrature.hpp"
Expand Down Expand Up @@ -46,6 +48,7 @@ DFieldMem<IdxRangeCoeffs, typename ExecSpace::memory_space> compute_coeffs_on_ma
Mapping& mapping,
DFieldMem<IdxRangeCoeffs, typename ExecSpace::memory_space>&& coefficients_alloc)
{
static_assert(has_2d_jacobian_v<Mapping, Coord<R, Theta>>);
using R = typename Mapping::curvilinear_tag_r;
using Theta = typename Mapping::curvilinear_tag_theta;
using IdxCoeffs = typename IdxRangeCoeffs::discrete_element_type;
Expand Down
111 changes: 36 additions & 75 deletions vendor/sll/include/sll/mapping/cartesian_to_pseudo_cartesian.hpp
Original file line number Diff line number Diff line change
@@ -1,16 +1,15 @@
// SPDX-License-Identifier: MIT
#pragma once

#include <sll/view.hpp>

/**
* @brief A class describing a mapping from a Cartesian domain to a pseudo-Cartesian.
* This operation is carried out using 2 curvilinear to Cartesian mappings.
* The Cartesian to pseudo-Cartesian mapping is obtained by combining these 2 mappings.
*/
template <class MappingToCartesian, class MappingToPseudoCartesian>
class CartesianToPseudoCartesian
: public Jacobian<ddc::Coordinate<
typename MappingToCartesian::cartesian_tag_x,
typename MappingToCartesian::cartesian_tag_y>>
{
static_assert(
std::is_same_v<
Expand All @@ -29,9 +28,6 @@ class CartesianToPseudoCartesian
typename MappingToCartesian::cartesian_tag_x,
typename MappingToCartesian::cartesian_tag_y>;

/// The type of the Jacobian matrix and its inverse.
using Matrix_2x2 = typename Jacobian<CartesianCoordinate>::Matrix_2x2;

/// The X dimension in the Cartesian geometry.
using cartesian_tag_x = typename MappingToCartesian::cartesian_tag_x;
/// The Y dimension in the Cartesian geometry.
Expand All @@ -42,6 +38,11 @@ class CartesianToPseudoCartesian
/// The Y dimension in the pseudo-Cartesian geometry.
using pseudo_cartesian_tag_y = typename MappingToPseudoCartesian::cartesian_tag_y;

/// The type of the argument of the function described by this mapping
using CoordArg = ddc::Coordinate<cartesian_tag_x, cartesian_tag_y>;
/// The type of the result of the function described by this mapping
using CoordResult = ddc::Coordinate<pseudo_cartesian_tag_x, pseudo_cartesian_tag_y>;

private:
using R = typename MappingToCartesian::curvilinear_tag_r;
using Theta = typename MappingToCartesian::curvilinear_tag_theta;
Expand Down Expand Up @@ -113,88 +114,68 @@ class CartesianToPseudoCartesian
}
}

/**
* @brief Compute the full Jacobian matrix from a Cartesian coordinate.
* @see jacobian_matrix
*
* @param[in] coord_cart The coordinate where we evaluate the Jacobian matrix.
* @param[out] J The Jacobian matrix returned.
*/
KOKKOS_INLINE_FUNCTION void jacobian_matrix(
CartesianCoordinate const& coord_cart,
Matrix_2x2& J) const final
{
if constexpr (std::is_invocable_r_v<CoordRTheta, MappingToCartesian, CartesianCoordinate>) {
CoordRTheta coord_rtheta = m_mapping_to_cartesian(coord_cart);
jacobian_matrix(coord_rtheta, J);
} else {
Kokkos::abort("The provided MappingToCartesian class does not allow the Jacobian "
"matrix to be calculated from a Cartesian coordinate.");
}
}

/**
* @brief Compute the determinant of the Jacobian matrix.
* @see jacobian_matrix
* @param[in] coord_cart The coordinate where we evaluate the Jacobian matrix.
* @param[in] coord_rtheta The coordinate where we evaluate the Jacobian matrix.
* @returns The determinant of the Jacobian matrix.
*/
KOKKOS_FUNCTION double jacobian(CartesianCoordinate const& coord_cart) const final
KOKKOS_FUNCTION double jacobian(CoordRTheta const& coord_rtheta) const
{
Matrix_2x2 J;
jacobian_matrix(coord_cart, J);
jacobian_matrix(coord_rtheta, J);
return determinant(J);
}

/**
* @brief Compute the (1,1) coefficient of the Jacobian matrix.
* @see jacobian_matrix
* @param[in] coord_cart The coordinate where we evaluate the Jacobian matrix.
* @param[in] coord_rtheta The coordinate where we evaluate the Jacobian matrix.
* @return The (1,1) coefficient of the Jacobian matrix.
*/
KOKKOS_FUNCTION double jacobian_11(CartesianCoordinate const& coord_cart) const final
KOKKOS_FUNCTION double jacobian_11(CoordRTheta const& coord_rtheta) const
{
Matrix_2x2 J;
jacobian_matrix(coord_cart, J);
jacobian_matrix(coord_rtheta, J);
return J[0][0];
}

/**
* @brief Compute the (1,2) coefficient of the Jacobian matrix.
* @see jacobian_matrix
* @param[in] coord_cart The coordinate where we evaluate the Jacobian matrix.
* @param[in] coord_rtheta The coordinate where we evaluate the Jacobian matrix.
* @return The (1,2) coefficient of the Jacobian matrix.
*/
KOKKOS_FUNCTION double jacobian_12(CartesianCoordinate const& coord_cart) const final
KOKKOS_FUNCTION double jacobian_12(CoordRTheta const& coord_rtheta) const
{
Matrix_2x2 J;
jacobian_matrix(coord_cart, J);
jacobian_matrix(coord_rtheta, J);
return J[0][1];
}

/**
* @brief Compute the (2,1) coefficient of the Jacobian matrix.
* @see jacobian_matrix
* @param[in] coord_cart The coordinate where we evaluate the Jacobian matrix.
* @param[in] coord_rtheta The coordinate where we evaluate the Jacobian matrix.
* @return The (2,1) coefficient of the Jacobian matrix.
*/
KOKKOS_FUNCTION double jacobian_21(CartesianCoordinate const& coord_cart) const final
KOKKOS_FUNCTION double jacobian_21(CoordRTheta const& coord_rtheta) const
{
Matrix_2x2 J;
jacobian_matrix(coord_cart, J);
jacobian_matrix(coord_rtheta, J);
return J[1][0];
}

/**
* @brief Compute the (2,2) coefficient of the Jacobian matrix.
* @see jacobian_matrix
* @param[in] coord_cart The coordinate where we evaluate the Jacobian matrix.
* @param[in] coord_rtheta The coordinate where we evaluate the Jacobian matrix.
* @return The (2,2) coefficient of the Jacobian matrix.
*/
KOKKOS_FUNCTION double jacobian_22(CartesianCoordinate const& coord_cart) const final
KOKKOS_FUNCTION double jacobian_22(CoordRTheta const& coord_rtheta) const
{
Matrix_2x2 J;
jacobian_matrix(coord_cart, J);
jacobian_matrix(coord_rtheta, J);
return J[1][1];
}

Expand Down Expand Up @@ -241,80 +222,60 @@ class CartesianToPseudoCartesian
}
}

/**
* @brief Compute the full inverse of the Jacobian matrix from a Cartesian coordinate.
* @see jacobian_matrix
*
* @param[in] coord_cart The coordinate where we evaluate the Jacobian matrix.
* @param[out] J The Jacobian matrix returned.
*/
KOKKOS_INLINE_FUNCTION void inv_jacobian_matrix(
CartesianCoordinate const& coord_cart,
Matrix_2x2& J) const final
{
if constexpr (std::is_invocable_r_v<CoordRTheta, MappingToCartesian, CartesianCoordinate>) {
CoordRTheta coord_rtheta = m_mapping_to_cartesian(coord_cart);
jacobian_matrix(coord_rtheta, J);
} else {
Kokkos::abort("The provided MappingToCartesian class does not allow the Jacobian "
"matrix to be calculated from a Cartesian coordinate.");
}
}

/**
* @brief Compute the (1,1) coefficient of the inverse of the Jacobian matrix.
* @see inv_jacobian_matrix
* @param[in] coord_cart The coordinate where we evaluate the inverse of the Jacobian matrix.
* @param[in] coord_rtheta The coordinate where we evaluate the inverse of the Jacobian matrix.
* @return The (1,1) coefficient of the inverse of the Jacobian matrix.
*/
KOKKOS_FUNCTION double inv_jacobian_11(CartesianCoordinate const& coord_cart) const final
KOKKOS_FUNCTION double inv_jacobian_11(CoordRTheta const& coord_rtheta) const
{
Matrix_2x2 J;
inv_jacobian_matrix(coord_cart, J);
inv_jacobian_matrix(coord_rtheta, J);
return J[0][0];
}

/**
* @brief Compute the (1,2) coefficient of the inverse of the Jacobian matrix.
* @see inv_jacobian_matrix
* @param[in] coord_cart The coordinate where we evaluate the inverse of the Jacobian matrix.
* @param[in] coord_rtheta The coordinate where we evaluate the inverse of the Jacobian matrix.
* @return The (1,2) coefficient of the inverse of the Jacobian matrix.
*/
KOKKOS_FUNCTION double inv_jacobian_12(CartesianCoordinate const& coord_cart) const final
KOKKOS_FUNCTION double inv_jacobian_12(CoordRTheta const& coord_rtheta) const
{
Matrix_2x2 J;
inv_jacobian_matrix(coord_cart, J);
inv_jacobian_matrix(coord_rtheta, J);
return J[0][1];
}

/**
* @brief Compute the (2,1) coefficient of the inverse of the Jacobian matrix.
* @see inv_jacobian_matrix
* @param[in] coord_cart The coordinate where we evaluate the inverse of the Jacobian matrix.
* @param[in] coord_rtheta The coordinate where we evaluate the inverse of the Jacobian matrix.
* @return The (2,1) coefficient of the inverse of the Jacobian matrix.
*/
KOKKOS_FUNCTION double inv_jacobian_21(CartesianCoordinate const& coord_cart) const final
KOKKOS_FUNCTION double inv_jacobian_21(CoordRTheta const& coord_rtheta) const
{
Matrix_2x2 J;
inv_jacobian_matrix(coord_cart, J);
inv_jacobian_matrix(coord_rtheta, J);
return J[1][0];
}

/**
* @brief Compute the (2,2) coefficient of the inverse of the Jacobian matrix.
* @see inv_jacobian_matrix
* @param[in] coord_cart The coordinate where we evaluate the inverse of the Jacobian matrix.
* @param[in] coord_rtheta The coordinate where we evaluate the inverse of the Jacobian matrix.
* @return The (2,2) coefficient of the inverse of the Jacobian matrix.
*/
KOKKOS_FUNCTION double inv_jacobian_22(CartesianCoordinate const& coord_cart) const final
KOKKOS_FUNCTION double inv_jacobian_22(CoordRTheta const& coord_rtheta) const
{
Matrix_2x2 J;
inv_jacobian_matrix(coord_cart, J);
inv_jacobian_matrix(coord_rtheta, J);
return J[1][1];
}
};

namespace detail {
namespace mapping_detail {
template <class ExecSpace, class MappingToPseudoCartesian, class MappingToCartesian>
struct MappingAccessibility<
ExecSpace,
Expand All @@ -324,4 +285,4 @@ struct MappingAccessibility<
&& MappingAccessibility<ExecSpace, MappingToCartesian>::value;
};

} // namespace detail
} // namespace mapping_detail
Loading

0 comments on commit dffd882

Please sign in to comment.