Skip to content

Commit

Permalink
Merge branch 'ebourne_fix_github_docker' into 'main'
Browse files Browse the repository at this point in the history 
Correct SHA tag for docker deployment

See merge request gysela-developpers/gyselalibxx!803

--------------------------------------------

Co-authored-by: Emily Bourne <[email protected]>
  • Loading branch information
@tpadioleau @EmilyBourne
tpadioleau and EmilyBourne committed Dec 11, 2024
1 parent 8a72d29 commit 66e5206
Show file tree
Hide file tree
Showing 7 changed files with 45 additions and 228 deletions.
2 changes: 1 addition & 1 deletion simulations/geometryRTheta/diocotron/diocotron.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,7 +164,7 @@ int main(int argc, char** argv)

PreallocatableSplineInterpolatorRTheta interpolator(builder, spline_evaluator);

AdvectionPhysicalDomain advection_domain(to_physical_mapping);
AdvectionDomain advection_domain(to_physical_mapping);

SplineFootFinder find_feet(
time_stepper,
Expand Down
2 changes: 1 addition & 1 deletion simulations/geometryRTheta/vortex_merger/vortex_merger.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,7 @@ int main(int argc, char** argv)

PreallocatableSplineInterpolatorRTheta interpolator(builder, spline_evaluator);

AdvectionPhysicalDomain advection_domain(to_physical_mapping);
AdvectionDomain advection_domain(to_physical_mapping);

SplineFootFinder find_feet(
time_stepper,
Expand Down
210 changes: 12 additions & 198 deletions src/geometryRTheta/advection/advection_domain.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@
/**
* @brief Define a domain for the advection.
*
* The natural advection domain is the physical domain (AdvectionPhysicalDomain),
* The natural advection domain is the physical domain (AdvectionDomain),
* where the studied equation is given.
* However, not all the mappings used are analytically invertible and inverting
* the Jacobian matrix of the mapping could be costly. That is why, we also
Expand All @@ -37,65 +37,38 @@
*/
template <class LogicalToPhysicalMapping>
class AdvectionDomain
{
public:
virtual ~AdvectionDomain() {};
};


/**
* @brief Define the physical domain for the advection.
*
* The physical domain can only be used for analytically invertible mappings.
*
* Here the advection field in the physical domain and the advection domain is the same.
* So the AdvectionPhysicalDomain::compute_advection_field return the same advection field.
*
* The method to advect is here given by:
* - compute @f$ (x,y)_{ij} = \mathcal{F} (r,\theta)_{ij} @f$,
* - advect
* - @f$ \text{feet}_{x, ij} = x_{ij} - dt A_x(x_{ij}, y_{ij}) @f$,
* - @f$ \text{feet}_{y, ij} = y_{ij} - dt A_y(x_{ij}, y_{ij}) @f$,
* - compute @f$ (\text{feet}_{r, ij}, \text{feet}_{\theta, ij})
* = \mathcal{F}^{-1} (\text{feet}_{x, ij}, \text{feet}_{y, ij}) @f$
* as @f$\mathcal{F} @f$ easily invertible.
*
*
*/
template <class LogicalToPhysicalMapping>
class AdvectionPhysicalDomain : public AdvectionDomain<LogicalToPhysicalMapping>
{
static_assert(is_analytical_mapping_v<LogicalToPhysicalMapping>);

private:
using PhysicalToLogicalMapping = inverse_mapping_t<LogicalToPhysicalMapping>;

public:
/**
* @brief The first dimension in the advection domain.
*/
using X_adv = X;
using X_adv = typename PhysicalToLogicalMapping::cartesian_tag_x;
/**
* @brief The second dimension in the advection domain.
*/
using Y_adv = Y;
using Y_adv = typename PhysicalToLogicalMapping::cartesian_tag_y;
/**
* @brief The coordinate type associated to the dimensions in the advection domain.
*/
using CoordXY_adv = Coord<X_adv, Y_adv>;

private:
using PhysicalToLogicalMapping = inverse_mapping_t<LogicalToPhysicalMapping>;

private:
LogicalToPhysicalMapping m_to_cartesian_mapping;
PhysicalToLogicalMapping m_to_curvilinear_mapping;

public:
/**
* @brief Instantiate a AdvectionPhysicalDomain advection domain.
* @brief Instantiate a AdvectionDomain advection domain.
*
* @param[in] to_physical_mapping
* The mapping from the logical domain to the physical domain.
*/
AdvectionPhysicalDomain(LogicalToPhysicalMapping const& to_physical_mapping)
AdvectionDomain(LogicalToPhysicalMapping const& to_physical_mapping)
: m_to_cartesian_mapping(to_physical_mapping)
, m_to_curvilinear_mapping(to_physical_mapping.get_inverse_mapping())
{
Expand Down Expand Up @@ -142,16 +115,15 @@ class AdvectionPhysicalDomain : public AdvectionDomain<LogicalToPhysicalMapping>
host_t<DConstVectorFieldRTheta<X_adv, Y_adv>> advection_field,
double dt) const
{
using namespace ddc;

IdxRangeRTheta const idx_range_rp = get_idx_range<GridR, GridTheta>(feet_coords_rp);
CoordXY coord_center(m_to_cartesian_mapping(CoordRTheta(0, 0)));

CoordXY_adv coord_center(m_to_cartesian_mapping(CoordRTheta(0, 0)));

ddc::for_each(idx_range_rp, [&](IdxRTheta const irp) {
CoordRTheta const coord_rp(feet_coords_rp(irp));
CoordXY const coord_xy = m_to_cartesian_mapping(coord_rp);
CoordXY_adv const coord_xy = m_to_cartesian_mapping(coord_rp);

CoordXY const feet_xy = coord_xy - dt * advection_field(irp);
CoordXY_adv const feet_xy = coord_xy - dt * advection_field(irp);

if (norm_inf(feet_xy - coord_center) < 1e-15) {
feet_coords_rp(irp) = CoordRTheta(0, 0);
Expand All @@ -164,161 +136,3 @@ class AdvectionPhysicalDomain : public AdvectionDomain<LogicalToPhysicalMapping>
});
}
};


/**
* @brief Define the pseudo-Cartesian domain for the advection.
*
* The pseudo-Cartesian domain is recommended for not analytically
* invertible mappings.
*
* We introduce a circular mapping @f$ \mathcal{G}@f$
* from the logical domain to the pseudo-Cartesian domain.
* The circular mapping is analytically invertible.
*
* The advection field is defined in the physical domain, so we need to
* define it in the pseudo-Cartesian domain before advecting
* (AdvectionPseudoCartesianDomain::compute_advection_field).
* To do so, we use the Jacobian matrix @f$ J_{\mathcal{F}}J_{\mathcal{G}}^{-1} @f$.
* This matrix is invertible :
*
* @f$ A_{\text{cart}} = (J_{\mathcal{F}}J_{\mathcal{G}}^{-1})^{-1} A @f$.
*
* The main difficulty is the center point. So for @f$ r \in [0, \varepsilon] @f$, we
* linearize the Jacobian matrix:
*
* @f$ (J_{\mathcal{F}}J_{\mathcal{G}}^{-1})^{-1} (r, \theta)
* = \left(1 - \frac{r}{\varepsilon} \right) (J_{\mathcal{F}}J_{\mathcal{G}}^{-1})^{-1} (0, \theta)
* + \frac{r}{\varepsilon} (J_{\mathcal{F}}J_{\mathcal{G}}^{-1})^{-1} (\varepsilon, \theta)
* @f$.
*
* The @f$ (J_{\mathcal{F}}J_{\mathcal{G}}^{-1})^{-1} (0, \theta) @f$ matrices are implemented
* in the Curvilinear2DToCartesian::to_pseudo_cartesian_jacobian_center_matrix child classes.
*
*
* The method to advect is here given by:
* - compute @f$ (x_{\text{cart}},y_{\text{cart}})_{ij} = \mathcal{G} (r,\theta)_{ij} @f$,
* - advect
* - @f$ \text{feet}_{x_{\text{cart}}, ij} = x_{\text{cart}, ij} - dt A_{\text{cart}, x}(x_{\text{cart}, ij}, y_{\text{cart}, ij}) @f$,
* - @f$ \text{feet}_{y_{\text{cart}}, ij} = y_{\text{cart}, ij} - dt A_{\text{cart}, y}(x_{\text{cart}, ij}, y_{\text{cart}, ij}) @f$,
* - compute @f$ (\text{feet}_{r, ij}, \text{feet}_{\theta, ij})
* = \mathcal{G}^{-1} (\text{feet}_{x_{\text{cart}}, ij}, \text{feet}_{y_{\text{cart}}, ij}) @f$.
*
*
*
* More details can be found in Edoardo Zoni's article
* (https://doi.org/10.1016/j.jcp.2019.108889).
*
* @see DiscreteToCartesian
*/
template <class LogicalToPhysicalMapping>
class AdvectionPseudoCartesianDomain : public AdvectionDomain<LogicalToPhysicalMapping>
{
public:
/**
* @brief Define a 2x2 matrix with an 2D array of an 2D array.
*/
using Matrix2x2 = std::array<std::array<double, 2>, 2>;

/**
* @brief The first dimension in the advection domain.
*/
using X_adv = X_pC;
/**
* @brief The second dimension in the advection domain.
*/
using Y_adv = Y_pC;
/**
* @brief The coordinate type associated to the dimensions in the advection domain.
*/
using CoordXY_adv = Coord<X_adv, Y_adv>;


private:
LogicalToPhysicalMapping const& m_to_cartesian_mapping;
double m_epsilon;

public:
/**
* @brief Instantiate an AdvectionPseudoCartesianDomain advection domain.
*
* @param[in] to_physical_mapping
* The mapping from the logical domain to the physical domain.
* @param[in] epsilon
* @f$ \varepsilon @f$ parameter used for the linearization of the
* advection field around the central point.
*/
AdvectionPseudoCartesianDomain(
LogicalToPhysicalMapping const& to_physical_mapping,
double epsilon = 1e-12)
: m_to_cartesian_mapping(to_physical_mapping)
, m_epsilon(epsilon) {};
~AdvectionPseudoCartesianDomain() {};

/**
* @brief Advect the characteristic feet.
*
* In the Backward Semi-Lagrangian method, the advection of a function
* uses the conservation along the characteristic property. So, we firstly
* compute the characteristic feet and then interpolate the function at this
* characteristic feet.
*
* The function implemented here deals with the computation of the characteristic feet.
* The IFootFinder class uses a time integration method to solve the characteristic
* equation.
* The BslAdvectionRTheta class calls advect_feet to compute the characteristic feet
* and interpolate the function we want to advect.
*
* The advect_feet implemented here computes only
*
* - @f$ (\text{feet}_r, \text{feet}_\theta) = \mathcal{G}^{-1} (\text{feet}_x, \text{feet}_y) @f$
*
* with
* - @f$ \text{feet}_x = x_{ij} - dt A_x(x_{ij}, y_{ij}) @f$,
* - @f$ \text{feet}_y = y_{ij} - dt A_y(x_{ij}, y_{ij}) @f$,
*
* and
* - @f$ (x,y)_{ij} = \mathcal{G}(r_{ij}, \theta_{ij})@f$, with @f$\{(r, \theta)_{ij}\}_{ij} @f$ the logical mesh points,
* - @f$ A @f$ the advection field in the advection domain,
* - @f$ \mathcal{G} @f$ the mapping from the logical domain to the advection domain.
*
* @param[in, out] feet_coords_rp
* The computed characteristic feet in the logical domain.
* On input: the points we want to advect.
* On output: the characteristic feet.
* @param[in] advection_field
* The advection field defined in the advection domain.
* @param[in] dt
* The time step.
*/
void advect_feet(
host_t<FieldRTheta<CoordRTheta>> feet_coords_rp,
host_t<DConstVectorFieldRTheta<X_adv, Y_adv>> const& advection_field,
double const dt) const
{
static_assert(
!std::is_same_v<LogicalToPhysicalMapping, CircularToCartesian<R, Theta, X, Y>>);
IdxRangeRTheta const idx_range_rp = get_idx_range(advection_field);

CircularToCartesian<R, Theta, X_adv, Y_adv> const pseudo_Cartesian_mapping;
CoordXY_adv const center_xy_pseudo_cart
= CoordXY_adv(pseudo_Cartesian_mapping(CoordRTheta(0., 0.)));

ddc::for_each(idx_range_rp, [&](IdxRTheta const irp) {
CoordRTheta const coord_rp(feet_coords_rp(irp));
CoordXY_adv const coord_xy_pseudo_cart = pseudo_Cartesian_mapping(coord_rp);
CoordXY_adv const feet_xy_pseudo_cart
= coord_xy_pseudo_cart - dt * advection_field(irp);

if (norm_inf(feet_xy_pseudo_cart - center_xy_pseudo_cart) < 1e-15) {
feet_coords_rp(irp) = CoordRTheta(0, 0);
} else {
CartesianToCircular<X_adv, Y_adv, R, Theta> const inv_pseudo_Cartesian_mapping;
feet_coords_rp(irp) = inv_pseudo_Cartesian_mapping(feet_xy_pseudo_cart);
ddc::select<Theta>(feet_coords_rp(irp)) = ddcHelper::restrict_to_idx_range(
ddc::select<Theta>(feet_coords_rp(irp)),
IdxRangeTheta(idx_range_rp));
}
});
}
};
10 changes: 6 additions & 4 deletions tests/geometryRTheta/advection_2d_rp/advection_all_tests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,7 @@ using CircularToCartMapping = CircularToCartesian<R, Theta, X, Y>;
using CzarnyToCartMapping = CzarnyToCartesian<R, Theta, X, Y>;
using CartToCircularMapping = CartesianToCircular<X, Y, R, Theta>;
using CartToCzarnyMapping = CartesianToCzarny<X, Y, R, Theta>;
using CircularToPseudoCartMapping = CircularToCartesian<R, Theta, X_pC, Y_pC>;
using DiscreteMappingBuilder
= DiscreteToCartesianBuilder<X, Y, SplineRThetaBuilder, SplineRThetaEvaluatorConstBound>;

Expand Down Expand Up @@ -321,6 +322,7 @@ int main(int argc, char** argv)

// SET THE DIFFERENT PARAMETERS OF THE TESTS ------------------------------------------------
CircularToCartMapping const from_circ_map;
CircularToPseudoCartMapping const to_pseudo_circ_map;
CartesianToCircular<X, Y, R, Theta> to_circ_map;
CzarnyToCartMapping const from_czarny_map(0.3, 1.4);
CartesianToCzarny<X, Y, R, Theta> const to_czarny_map(0.3, 1.4);
Expand All @@ -331,10 +333,10 @@ int main(int argc, char** argv)
spline_evaluator_extrapol);
DiscreteToCartesian const discrete_czarny_map = discrete_czarny_map_builder();

AdvectionPhysicalDomain<CircularToCartMapping> const physical_circular_mapping(from_circ_map);
AdvectionPhysicalDomain<CzarnyToCartMapping> const physical_czarny_mapping(from_czarny_map);
AdvectionPseudoCartesianDomain<CzarnyToCartMapping> const pseudo_cartesian_czarny_mapping(
from_czarny_map);
AdvectionDomain<CircularToCartMapping> const physical_circular_mapping(from_circ_map);
AdvectionDomain<CzarnyToCartMapping> const physical_czarny_mapping(from_czarny_map);
AdvectionDomain<CircularToPseudoCartMapping> const pseudo_cartesian_czarny_mapping(
to_pseudo_circ_map);

std::tuple simulations = std::make_tuple(
SimulationParameters(
Expand Down
29 changes: 15 additions & 14 deletions tests/geometryRTheta/advection_2d_rp/advection_selected_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
#include "ddc_helper.hpp"
#include "euler.hpp"
#include "geometry.hpp"
#include "geometry_pseudo_cartesian.hpp"
#include "input.hpp"
#include "itimestepper.hpp"
#include "mesh_builder.hpp"
Expand All @@ -42,21 +43,19 @@ namespace fs = std::filesystem;

namespace {
#if defined(CIRCULAR_MAPPING_PHYSICAL)
using X_adv = typename AdvectionPhysicalDomain<CircularToCartesian<R, Theta, X, Y>>::X_adv;
using Y_adv = typename AdvectionPhysicalDomain<CircularToCartesian<R, Theta, X, Y>>::Y_adv;
using X_adv = X;
using Y_adv = Y;
#elif defined(CZARNY_MAPPING_PHYSICAL)
using X_adv = typename AdvectionPhysicalDomain<CzarnyToCartesian<R, Theta, X, Y>>::X_adv;
using Y_adv = typename AdvectionPhysicalDomain<CzarnyToCartesian<R, Theta, X, Y>>::Y_adv;
using X_adv = X;
using Y_adv = Y;

#elif defined(CZARNY_MAPPING_PSEUDO_CARTESIAN)
using X_adv = typename AdvectionPseudoCartesianDomain<CzarnyToCartesian<R, Theta, X, Y>>::X_adv;
using Y_adv = typename AdvectionPseudoCartesianDomain<CzarnyToCartesian<R, Theta, X, Y>>::Y_adv;
using X_adv = X_pC;
using Y_adv = Y_pC;

#elif defined(DISCRETE_MAPPING_PSEUDO_CARTESIAN)
using X_adv = typename AdvectionPseudoCartesianDomain<
DiscreteToCartesian<X, Y, SplineRThetaEvaluatorConstBound>>::X_adv;
using Y_adv = typename AdvectionPseudoCartesianDomain<
DiscreteToCartesian<X, Y, SplineRThetaEvaluatorConstBound>>::Y_adv;
using X_adv = X_pC;
using Y_adv = Y_pC;
#endif

} //end namespace
Expand Down Expand Up @@ -167,7 +166,7 @@ int main(int argc, char** argv)
CircularToCartesian<R, Theta, X, Y> analytical_mapping;
CircularToCartesian<R, Theta, X, Y> to_physical_mapping;
CartesianToCircular<X, Y, R, Theta> to_logical_mapping;
AdvectionPhysicalDomain advection_domain(analytical_mapping);
AdvectionDomain advection_domain(analytical_mapping);
std::string const mapping_name = "CIRCULAR";
std::string const adv_domain_name = "PHYSICAL";
key += "circular_physical";
Expand All @@ -180,7 +179,7 @@ int main(int argc, char** argv)
CzarnyToCartesian<R, Theta, X, Y> analytical_mapping(czarny_e, czarny_epsilon);
CzarnyToCartesian<R, Theta, X, Y> to_physical_mapping(czarny_e, czarny_epsilon);
CartesianToCzarny<X, Y, R, Theta> to_logical_mapping(czarny_e, czarny_epsilon);
AdvectionPhysicalDomain advection_domain(analytical_mapping);
AdvectionDomain advection_domain(analytical_mapping);
std::string const mapping_name = "CZARNY";
std::string const adv_domain_name = "PHYSICAL";
key += "czarny_physical";
Expand All @@ -189,7 +188,8 @@ int main(int argc, char** argv)
CzarnyToCartesian<R, Theta, X, Y> analytical_mapping(czarny_e, czarny_epsilon);
CzarnyToCartesian<R, Theta, X, Y> to_physical_mapping(czarny_e, czarny_epsilon);
CartesianToCzarny<X, Y, R, Theta> to_logical_mapping(czarny_e, czarny_epsilon);
AdvectionPseudoCartesianDomain advection_domain(to_physical_mapping);
CircularToCartesian<R, Theta, X_pC, Y_pC> logical_to_pseudo_cart_mapping;
AdvectionDomain advection_domain(logical_to_pseudo_cart_mapping);
std::string const mapping_name = "CZARNY";
std::string const adv_domain_name = "PSEUDO CARTESIAN";
key += "czarny_pseudo_cartesian";
Expand All @@ -204,7 +204,8 @@ int main(int argc, char** argv)
builder,
spline_evaluator_extrapol);
DiscreteToCartesian to_physical_mapping = mapping_builder();
AdvectionPseudoCartesianDomain advection_domain(to_physical_mapping);
CircularToCartesian<R, Theta, X_pC, Y_pC> logical_to_pseudo_cart_mapping;
AdvectionDomain advection_domain(logical_to_pseudo_cart_mapping);
std::string const mapping_name = "DISCRETE";
std::string const adv_domain_name = "PSEUDO CARTESIAN";
key += "discrete_pseudo_cartesian";
Expand Down
Loading

0 comments on commit 66e5206

Please sign in to comment.