Skip to content

Commit

Permalink
Merge branch 'ebourne_split_metric_tensor' into 'main'
Browse files Browse the repository at this point in the history
Create metric tensor operator

See merge request gysela-developpers/gyselalibxx!683
  • Loading branch information
EmilyBourne committed Sep 4, 2024
1 parent a357a82 commit 399a8f4
Show file tree
Hide file tree
Showing 7 changed files with 169 additions and 126 deletions.
6 changes: 4 additions & 2 deletions src/geometryRTheta/advection/bsl_advection_rp.hpp
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
// SPDX-License-Identifier: MIT

#pragma once
#include <sll/mapping/analytical_invertible_curvilinear2d_to_cartesian.hpp>
#include <sll/mapping/curvilinear2d_to_cartesian.hpp>
#include <sll/mapping/metric_tensor.hpp>

#include "advection_domain.hpp"
#include "ddc_alias_inline_functions.hpp"
Expand Down Expand Up @@ -165,13 +165,15 @@ class BslAdvectionRTheta : public IAdvectionRTheta
// Convert advection field on RTheta to advection field on XY
DVectorFieldMemRTheta<X, Y> advection_field_xy(grid);

MetricTensor<Mapping, CoordRTheta> metric_tensor(m_mapping);

ddc::for_each(grid_without_Opoint, [&](IdxRTheta const irp) {
CoordRTheta const coord_rp(ddc::coordinate(irp));

std::array<std::array<double, 2>, 2> J; // Jacobian matrix
m_mapping.jacobian_matrix(coord_rp, J);
std::array<std::array<double, 2>, 2> G; // Metric tensor
m_mapping.metric_tensor(coord_rp, G);
metric_tensor(G, coord_rp);

ddcHelper::get<X>(advection_field_xy)(irp)
= ddcHelper::get<R>(advection_field_rp)(irp) * J[1][1] / std::sqrt(G[1][1])
Expand Down
11 changes: 7 additions & 4 deletions src/geometryRTheta/advection_field/advection_field_rp.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
#pragma once
#include <ddc/ddc.hpp>

#include <sll/mapping/metric_tensor.hpp>
#include <sll/polar_spline.hpp>
#include <sll/polar_spline_evaluator.hpp>

Expand Down Expand Up @@ -69,8 +70,8 @@
* 2- In the second case, the advection field along the logical index range axis
* is computed with
* - @f$ \nabla \phi = \sum_{i,j} \partial_{x_i} f g^{ij} \sqrt{g_{jj}} \hat{e}_j@f$,
* - with @f$g^{ij}@f$, the coefficients of the inverse metrix tensor,
* - @f$g_{jj}@f$, the coefficients of the metrix tensor,
* - with @f$g^{ij}@f$, the coefficients of the inverse metric tensor,
* - @f$g_{jj}@f$, the coefficients of the metric tensor,
* - @f$\hat{e}_j@f$, the normalized covariants vectors.
*
* Then, we compute @f$ E = -\nabla \phi @f$ and @f$A = E \wedge e_z@f$.
Expand Down Expand Up @@ -474,16 +475,18 @@ class AdvectionFieldFinder
get_const_field(coords),
get_const_field(electrostatic_potential_coef));

MetricTensor<Mapping, CoordRTheta> metric_tensor(m_mapping);

// > computation of the advection field
ddc::for_each(grid_without_Opoint, [&](IdxRTheta const irp) {
CoordRTheta const coord_rp(ddc::coordinate(irp));

Matrix_2x2 J; // Jacobian matrix
m_mapping.jacobian_matrix(coord_rp, J);
Matrix_2x2 inv_G; // Inverse metric tensor
m_mapping.inverse_metric_tensor(coord_rp, inv_G);
metric_tensor.inverse(inv_G, coord_rp);
Matrix_2x2 G; // Metric tensor
m_mapping.metric_tensor(coord_rp, G);
metric_tensor(G, coord_rp);

// E = -grad phi
double const electric_field_r
Expand Down
6 changes: 5 additions & 1 deletion src/geometryRTheta/poisson/polarpoissonlikesolver.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -5,10 +5,12 @@
#include <ddc/ddc.hpp>

#include <sll/gauss_legendre_integration.hpp>
#include <sll/mapping/metric_tensor.hpp>
#include <sll/math_tools.hpp>
#include <sll/matrix_batch_csr.hpp>
#include <sll/polar_spline.hpp>
#include <sll/polar_spline_evaluator.hpp>
#include <sll/view.hpp>

#include "ddc_alias_inline_functions.hpp"
#include "ddc_aliases.hpp"
Expand Down Expand Up @@ -982,12 +984,14 @@ class PolarSplineFEMPoissonLikeSolver
trial_bspline_val_and_deriv,
trial_bspline_val_and_deriv_theta);

MetricTensor<Mapping, CoordRTheta> metric_tensor(mapping);

// Assemble the weak integral element
return int_volume(ir, ip)
* (alpha
* dot_product(
basis_gradient_test_space,
mapping.to_covariant(basis_gradient_trial_space, coord))
metric_tensor.to_covariant(basis_gradient_trial_space, coord))
+ beta * basis_val_test_space * basis_val_trial_space);
}

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -222,13 +222,15 @@ TEST(AdvectionFieldRThetaComputation, TestAdvectionFieldFinder)

// > Compare the advection field computed on RTheta to the advection field computed on XY
DVectorFieldMemRTheta<X, Y> difference_between_fields_xy_and_rp(grid);

MetricTensor<Mapping, CoordRTheta> metric_tensor(mapping);
ddc::for_each(grid_without_Opoint, [&](IdxRTheta const irp) {
CoordRTheta const coord_rp(ddc::coordinate(irp));

std::array<std::array<double, 2>, 2> J; // Jacobian matrix
mapping.jacobian_matrix(coord_rp, J);
std::array<std::array<double, 2>, 2> G; // Metric tensor
mapping.metric_tensor(coord_rp, G);
metric_tensor(G, coord_rp);

// computation made in BslAdvectionRTheta operator:
ddcHelper::get<X>(advection_field_xy_from_rp)(irp)
Expand Down
71 changes: 0 additions & 71 deletions vendor/sll/include/sll/mapping/curvilinear2d_to_cartesian.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -276,75 +276,4 @@ class Curvilinear2DToCartesian
assert(fabs(jacobian(coord)) >= 1e-15);
return jacobian_11(coord) / jacobian(coord);
}


/**
* @brief Compute the metric tensor assignd to the mapping.
*
* The metric tensor matrix is defined for mapping whose the Jacobian matrix is called
* @f$ J_{\mathcal{F}} @f$ as @f$ G = (J_{\mathcal{F}})^T J_{\mathcal{F}} @f$.
*
* @param[in] coord
* The coordinate where we evaluate the metric tensor.
* @param[out] matrix
* The metric tensor matrix.
*/
virtual void metric_tensor(ddc::Coordinate<R, Theta> const& coord, Matrix_2x2& matrix) const
{
const double J_rr = jacobian_11(coord);
const double J_rp = jacobian_12(coord);
const double J_pr = jacobian_21(coord);
const double J_pp = jacobian_22(coord);
matrix[0][0] = (J_rr * J_rr + J_pr * J_pr);
matrix[0][1] = (J_rr * J_rp + J_pr * J_pp);
matrix[1][0] = (J_rr * J_rp + J_pr * J_pp);
matrix[1][1] = (J_rp * J_rp + J_pp * J_pp);
}

/**
* @brief Compute the inverse metric tensor associated to the mapping.
*
* @param[in] coord
* The coordinate where we evaluate the metric tensor.
* @param[out] matrix
* The metric tensor matrix.
*/
virtual void inverse_metric_tensor(ddc::Coordinate<R, Theta> const& coord, Matrix_2x2& matrix)
const
{
assert(fabs(ddc::get<R>(coord)) >= 1e-15);
const double J_rr = jacobian_11(coord);
const double J_rp = jacobian_12(coord);
const double J_pr = jacobian_21(coord);
const double J_pp = jacobian_22(coord);
const double jacob_2 = jacobian(coord) * jacobian(coord);
matrix[0][0] = (J_rp * J_rp + J_pp * J_pp) / jacob_2;
matrix[0][1] = (-J_rr * J_rp - J_pr * J_pp) / jacob_2;
matrix[1][0] = (-J_rr * J_rp - J_pr * J_pp) / jacob_2;
matrix[1][1] = (J_rr * J_rr + J_pr * J_pr) / jacob_2;
}

/**
* @brief Compute the covariant vector from the contravariant vector
*
* @param[in] contravariant_vector
* The metric tensor matrix.
* @param[in] coord
* The coordinate where we want to compute the convariant vector.
*
* @return A vector of the covariant
*/
std::array<double, 2> to_covariant(
std::array<double, 2> const& contravariant_vector,
ddc::Coordinate<R, Theta> const& coord) const
{
Matrix_2x2 inv_metric_tensor;
inverse_metric_tensor(coord, inv_metric_tensor);
std::array<double, 2> covariant_vector;
covariant_vector[0] = inv_metric_tensor[0][0] * contravariant_vector[0]
+ inv_metric_tensor[0][1] * contravariant_vector[1];
covariant_vector[1] = inv_metric_tensor[1][0] * contravariant_vector[0]
+ inv_metric_tensor[1][1] * contravariant_vector[1];
return covariant_vector;
}
};
97 changes: 97 additions & 0 deletions vendor/sll/include/sll/mapping/metric_tensor.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,97 @@
#pragma once
#include <type_traits>

#include <sll/view.hpp>

/**
* @brief An operator for calculating the metric tensor.
* @tparam Mapping The mapping providing the Jacobian operator.
* @tparam PositionCoordinate The coordinate type where the metric tensor can be evaluated.
*/
template <class Mapping, class PositionCoordinate>
class MetricTensor
{
public:
/// The type of the Jacobian matrix and its inverse
using Matrix_2x2 = std::array<std::array<double, 2>, 2>;

private:
Mapping m_mapping;

public:
/**
* @brief A constructor for the metric tensor operator.
*
* @param[in] mapping The mapping which can be used to calculate the Jacobian.
*/
MetricTensor(Mapping mapping) : m_mapping(mapping) {}

/**
* @brief Compute the metric tensor assignd to the mapping.
*
* The metric tensor matrix is defined as:
* @f$ G = (J_{\mathcal{F}})^T J_{\mathcal{F}} @f$.
* with @f$ J_{\mathcal{F}} @f$ the Jacobian matrix.
*
* @param[in] coord
* The coordinate where we evaluate the metric tensor.
* @param[out] matrix
* The metric tensor matrix.
*/
void operator()(Matrix_2x2& matrix, PositionCoordinate const& coord) const
{
const double J_11 = m_mapping.jacobian_11(coord);
const double J_12 = m_mapping.jacobian_12(coord);
const double J_21 = m_mapping.jacobian_21(coord);
const double J_22 = m_mapping.jacobian_22(coord);
matrix[0][0] = (J_11 * J_11 + J_21 * J_21);
matrix[0][1] = (J_11 * J_12 + J_21 * J_22);
matrix[1][0] = (J_11 * J_12 + J_21 * J_22);
matrix[1][1] = (J_12 * J_12 + J_22 * J_22);
}

/**
* @brief Compute the inverse metric tensor associated to the mapping.
*
* @param[in] coord
* The coordinate where we evaluate the metric tensor.
* @param[out] matrix
* The metric tensor matrix.
*/
void inverse(Matrix_2x2& matrix, PositionCoordinate const& coord) const
{
const double J_11 = m_mapping.jacobian_11(coord);
const double J_12 = m_mapping.jacobian_12(coord);
const double J_21 = m_mapping.jacobian_21(coord);
const double J_22 = m_mapping.jacobian_22(coord);
const double jacob_2 = m_mapping.jacobian(coord) * m_mapping.jacobian(coord);
matrix[0][0] = (J_12 * J_12 + J_22 * J_22) / jacob_2;
matrix[0][1] = (-J_11 * J_12 - J_21 * J_22) / jacob_2;
matrix[1][0] = (-J_11 * J_12 - J_21 * J_22) / jacob_2;
matrix[1][1] = (J_11 * J_11 + J_21 * J_21) / jacob_2;
}

/**
* @brief Compute the covariant vector from the contravariant vector
*
* @param[in] contravariant_vector
* The metric tensor matrix.
* @param[in] coord
* The coordinate where we want to compute the convariant vector.
*
* @return A vector of the covariant
*/
std::array<double, 2> to_covariant(
std::array<double, 2> const& contravariant_vector,
PositionCoordinate const& coord) const
{
Matrix_2x2 inv_metric_tensor;
inverse(inv_metric_tensor, coord);
std::array<double, 2> covariant_vector;
covariant_vector[0] = inv_metric_tensor[0][0] * contravariant_vector[0]
+ inv_metric_tensor[0][1] * contravariant_vector[1];
covariant_vector[1] = inv_metric_tensor[1][0] * contravariant_vector[0]
+ inv_metric_tensor[1][1] * contravariant_vector[1];
return covariant_vector;
}
};
Loading

0 comments on commit 399a8f4

Please sign in to comment.