Skip to content

Commit

Permalink
Split mappings so they represent one equation
Browse files Browse the repository at this point in the history 
Closes #409

See merge request gysela-developpers/gyselalibxx!779

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

Co-authored-by: Emily Bourne <[email protected]>
  • Loading branch information
@Geoflow @EmilyBourne
Geoflow and EmilyBourne committed Dec 3, 2024
1 parent 068e04a commit c3caa14
Show file tree
Hide file tree
Showing 35 changed files with 952 additions and 414 deletions.
32 changes: 19 additions & 13 deletions simulations/geometryRTheta/diocotron/diocotron.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,8 @@

#include <ddc/ddc.hpp>

#include <sll/mapping/cartesian_to_circular.hpp>
#include <sll/mapping/circular_to_cartesian.hpp>
#include <sll/mapping/czarny_to_cartesian.hpp>
#include <sll/mapping/discrete_mapping_builder.hpp>
#include <sll/mapping/discrete_to_cartesian.hpp>

Expand Down Expand Up @@ -51,7 +51,8 @@ using PoissonSolver = PolarSplineFEMPoissonLikeSolver<
SplineRThetaEvaluatorNullBound>;
using DiscreteMappingBuilder
= DiscreteToCartesianBuilder<X, Y, SplineRThetaBuilder, SplineRThetaEvaluatorConstBound>;
using Mapping = CircularToCartesian<X, Y, R, Theta>;
using LogicalToPhysicalMapping = CircularToCartesian<R, Theta, X, Y>;
using PhysicalToLogicalMapping = CartesianToCircular<X, Y, R, Theta>;

namespace fs = std::filesystem;

Expand Down Expand Up @@ -111,10 +112,11 @@ int main(int argc, char** argv)
ddc::PeriodicExtrapolationRule<Theta>(),
ddc::PeriodicExtrapolationRule<Theta>());

const Mapping mapping;
const LogicalToPhysicalMapping to_physical_mapping;
const PhysicalToLogicalMapping to_logical_mapping;
DiscreteMappingBuilder const discrete_mapping_builder(
Kokkos::DefaultHostExecutionSpace(),
mapping,
to_physical_mapping,
builder,
spline_evaluator_extrapol);
DiscreteToCartesian const discrete_mapping = discrete_mapping_builder();
Expand Down Expand Up @@ -161,12 +163,16 @@ int main(int argc, char** argv)

PreallocatableSplineInterpolatorRTheta interpolator(builder, spline_evaluator);

AdvectionPhysicalDomain advection_domain(mapping);
AdvectionPhysicalDomain advection_domain(to_physical_mapping, to_logical_mapping);

SplineFootFinder
find_feet(time_stepper, advection_domain, mapping, builder, spline_evaluator_extrapol);
SplineFootFinder find_feet(
time_stepper,
advection_domain,
to_physical_mapping,
builder,
spline_evaluator_extrapol);

BslAdvectionRTheta advection_operator(interpolator, find_feet, mapping);
BslAdvectionRTheta advection_operator(interpolator, find_feet, to_physical_mapping);



Expand Down Expand Up @@ -195,15 +201,15 @@ int main(int argc, char** argv)
// --- Predictor corrector operator ---------------------------------------------------------------
#if defined(PREDCORR)
BslPredCorrRTheta predcorr_operator(
mapping,
to_physical_mapping,
advection_operator,
builder,
spline_evaluator,
poisson_solver);
#elif defined(EXPLICIT_PREDCORR)
BslExplicitPredCorrRTheta predcorr_operator(
advection_domain,
mapping,
to_physical_mapping,
advection_operator,
mesh_rp,
builder,
Expand All @@ -213,7 +219,7 @@ int main(int argc, char** argv)
#elif defined(IMPLICIT_PREDCORR)
BslImplicitPredCorrRTheta predcorr_operator(
advection_domain,
mapping,
to_physical_mapping,
advection_operator,
mesh_rp,
builder,
Expand Down Expand Up @@ -267,10 +273,10 @@ int main(int argc, char** argv)
host_t<FieldMemRTheta<CoordY>> coords_y(mesh_rp);
host_t<DFieldMemRTheta> jacobian(mesh_rp);
ddc::for_each(mesh_rp, [&](IdxRTheta const irp) {
CoordXY coords_xy = mapping(ddc::coordinate(irp));
CoordXY coords_xy = to_physical_mapping(ddc::coordinate(irp));
coords_x(irp) = ddc::select<X>(coords_xy);
coords_y(irp) = ddc::select<Y>(coords_xy);
jacobian(irp) = mapping.jacobian(ddc::coordinate(irp));
jacobian(irp) = to_physical_mapping.jacobian(ddc::coordinate(irp));
});


Expand Down
31 changes: 19 additions & 12 deletions simulations/geometryRTheta/vortex_merger/vortex_merger.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,8 @@ using PoissonSolver = PolarSplineFEMPoissonLikeSolver<
SplineRThetaEvaluatorNullBound>;
using DiscreteMappingBuilder
= DiscreteToCartesianBuilder<X, Y, SplineRThetaBuilder, SplineRThetaEvaluatorConstBound>;
using CircularMapping = CircularToCartesian<X, Y, R, Theta>;
using LogicalToPhysicalMapping = CircularToCartesian<R, Theta, X, Y>;
using PhysicalToLogicalMapping = CartesianToCircular<X, Y, R, Theta>;

} // end namespace

Expand Down Expand Up @@ -109,10 +110,11 @@ int main(int argc, char** argv)
ddc::PeriodicExtrapolationRule<Theta>(),
ddc::PeriodicExtrapolationRule<Theta>());

const CircularMapping mapping;
const LogicalToPhysicalMapping to_physical_mapping;
const PhysicalToLogicalMapping to_logical_mapping;
DiscreteMappingBuilder const discrete_mapping_builder(
Kokkos::DefaultHostExecutionSpace(),
mapping,
to_physical_mapping,
builder,
spline_evaluator_extrapol);
DiscreteToCartesian const discrete_mapping = discrete_mapping_builder();
Expand All @@ -139,12 +141,16 @@ int main(int argc, char** argv)

PreallocatableSplineInterpolatorRTheta interpolator(builder, spline_evaluator);

AdvectionPhysicalDomain advection_domain(mapping);
AdvectionPhysicalDomain advection_domain(to_physical_mapping, to_logical_mapping);

SplineFootFinder
find_feet(time_stepper, advection_domain, mapping, builder, spline_evaluator_extrapol);
SplineFootFinder find_feet(
time_stepper,
advection_domain,
to_physical_mapping,
builder,
spline_evaluator_extrapol);

BslAdvectionRTheta advection_operator(interpolator, find_feet, mapping);
BslAdvectionRTheta advection_operator(interpolator, find_feet, to_physical_mapping);



Expand Down Expand Up @@ -173,7 +179,7 @@ int main(int argc, char** argv)
// --- Predictor corrector operator ---------------------------------------------------------------
BslImplicitPredCorrRTheta predcorr_operator(
advection_domain,
mapping,
to_physical_mapping,
advection_operator,
grid,
builder,
Expand Down Expand Up @@ -225,23 +231,24 @@ int main(int argc, char** argv)
host_t<FieldMemRTheta<CoordY>> coords_y(grid);
host_t<DFieldMemRTheta> jacobian(grid);
ddc::for_each(grid, [&](IdxRTheta const irp) {
CoordXY coords_xy = mapping(ddc::coordinate(irp));
CoordXY coords_xy = to_physical_mapping(ddc::coordinate(irp));
coords_x(irp) = ddc::select<X>(coords_xy);
coords_y(irp) = ddc::select<Y>(coords_xy);
jacobian(irp) = mapping.jacobian(ddc::coordinate(irp));
jacobian(irp) = to_physical_mapping.jacobian(ddc::coordinate(irp));
});

double const tau(1e-10);
double const phi_max(1.);


VortexMergerEquilibria equilibrium(mapping, grid, builder, spline_evaluator, poisson_solver);
VortexMergerEquilibria
equilibrium(to_physical_mapping, grid, builder, spline_evaluator, poisson_solver);
std::function<double(double const)> const function = [&](double const x) { return x * x; };
host_t<DFieldMemRTheta> rho_eq(grid);
equilibrium.set_equilibrium(rho_eq, function, phi_max, tau);


VortexMergerDensitySolution solution(mapping);
VortexMergerDensitySolution solution(to_physical_mapping);
host_t<DFieldMemRTheta> rho(grid);
solution.set_initialisation(rho, rho_eq, eps, sigma, x_star_1, y_star_1, x_star_2, y_star_2);

Expand Down
48 changes: 30 additions & 18 deletions src/geometryRTheta/advection/advection_domain.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
#include <cassert>
#include <typeinfo>

#include <sll/mapping/cartesian_to_circular.hpp>
#include <sll/mapping/circular_to_cartesian.hpp>
#include <sll/mapping/curvilinear2d_to_cartesian.hpp>

Expand Down Expand Up @@ -35,7 +36,7 @@
* @see BslAdvectionRTheta
* @see IFootFinder
*/
template <class Mapping>
template <class LogicalToPhysicalMapping>
class AdvectionDomain
{
public:
Expand All @@ -62,8 +63,8 @@ class AdvectionDomain
*
*
*/
template <class Mapping>
class AdvectionPhysicalDomain : public AdvectionDomain<Mapping>
template <class LogicalToPhysicalMapping, class PhysicalToLogicalMapping>
class AdvectionPhysicalDomain : public AdvectionDomain<LogicalToPhysicalMapping>
{
public:
/**
Expand All @@ -80,16 +81,23 @@ class AdvectionPhysicalDomain : public AdvectionDomain<Mapping>
using CoordXY_adv = Coord<X_adv, Y_adv>;

private:
Mapping const& m_mapping;
LogicalToPhysicalMapping const& m_to_cartesian_mapping;
PhysicalToLogicalMapping const& m_to_curvilinear_mapping;

public:
/**
* @brief Instantiate a AdvectionPhysicalDomain advection domain.
*
* @param[in] mapping
* @param[in] to_physical_mapping
* The mapping from the logical domain to the physical domain.
* @param[in] to_logical_mapping
* The mapping from the physical domain to the logical domain.
*/
AdvectionPhysicalDomain(Mapping const& mapping) : m_mapping(mapping) {};
AdvectionPhysicalDomain(
LogicalToPhysicalMapping const& to_physical_mapping,
PhysicalToLogicalMapping const& to_logical_mapping)
: m_to_cartesian_mapping(to_physical_mapping)
, m_to_curvilinear_mapping(to_logical_mapping) {};
~AdvectionPhysicalDomain() {};

/**
Expand Down Expand Up @@ -136,18 +144,18 @@ class AdvectionPhysicalDomain : public AdvectionDomain<Mapping>
using namespace ddc;

IdxRangeRTheta const idx_range_rp = get_idx_range<GridR, GridTheta>(feet_coords_rp);
CoordXY coord_center(m_mapping(CoordRTheta(0, 0)));
CoordXY 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_mapping(coord_rp);
CoordXY const coord_xy = m_to_cartesian_mapping(coord_rp);

CoordXY 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);
} else {
feet_coords_rp(irp) = m_mapping(feet_xy);
feet_coords_rp(irp) = m_to_curvilinear_mapping(feet_xy);
ddc::select<Theta>(feet_coords_rp(irp)) = ddcHelper::restrict_to_idx_range(
ddc::select<Theta>(feet_coords_rp(irp)),
IdxRangeTheta(idx_range_rp));
Expand Down Expand Up @@ -203,8 +211,8 @@ class AdvectionPhysicalDomain : public AdvectionDomain<Mapping>
* @see Curvilinear2DToCartesian
* @see DiscreteToCartesian
*/
template <class Mapping>
class AdvectionPseudoCartesianDomain : public AdvectionDomain<Mapping>
template <class LogicalToPhysicalMapping>
class AdvectionPseudoCartesianDomain : public AdvectionDomain<LogicalToPhysicalMapping>
{
public:
/**
Expand All @@ -227,21 +235,23 @@ class AdvectionPseudoCartesianDomain : public AdvectionDomain<Mapping>


private:
Mapping const& m_mapping;
LogicalToPhysicalMapping const& m_to_cartesian_mapping;
double m_epsilon;

public:
/**
* @brief Instantiate an AdvectionPseudoCartesianDomain advection domain.
*
* @param[in] mapping
* @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(Mapping const& mapping, double epsilon = 1e-12)
: m_mapping(mapping)
AdvectionPseudoCartesianDomain(
LogicalToPhysicalMapping const& to_physical_mapping,
double epsilon = 1e-12)
: m_to_cartesian_mapping(to_physical_mapping)
, m_epsilon(epsilon) {};
~AdvectionPseudoCartesianDomain() {};

Expand Down Expand Up @@ -286,10 +296,11 @@ class AdvectionPseudoCartesianDomain : public AdvectionDomain<Mapping>
host_t<DConstVectorFieldRTheta<X_adv, Y_adv>> const& advection_field,
double const dt) const
{
static_assert(!std::is_same_v<Mapping, CircularToCartesian<X, Y, R, Theta>>);
static_assert(
!std::is_same_v<LogicalToPhysicalMapping, CircularToCartesian<R, Theta, X, Y>>);
IdxRangeRTheta const idx_range_rp = get_idx_range(advection_field);

CircularToCartesian<X_adv, Y_adv, R, Theta> const pseudo_Cartesian_mapping;
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.)));

Expand All @@ -302,7 +313,8 @@ class AdvectionPseudoCartesianDomain : public AdvectionDomain<Mapping>
if (norm_inf(feet_xy_pseudo_cart - center_xy_pseudo_cart) < 1e-15) {
feet_coords_rp(irp) = CoordRTheta(0, 0);
} else {
feet_coords_rp(irp) = pseudo_Cartesian_mapping(feet_xy_pseudo_cart);
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));
Expand Down
2 changes: 1 addition & 1 deletion src/geometryRTheta/advection/spline_foot_finder.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ class SplineFootFinder : public IFootFinder
*/
using Y_adv = typename AdvectionDomain::Y_adv;

using CircularToPseudoCartesian = CircularToCartesian<X_adv, Y_adv, R, Theta>;
using CircularToPseudoCartesian = CircularToCartesian<R, Theta, X_adv, Y_adv>;
using AdvectionMapping = CartesianToPseudoCartesian<Mapping, CircularToPseudoCartesian>;

TimeStepper const& m_time_stepper;
Expand Down
Loading

0 comments on commit c3caa14

Please sign in to comment.