Skip to content

Commit

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

See merge request gysela-developpers/gyselalibxx!618
  • Loading branch information
EmilyBourne committed Jul 30, 2024
1 parent 5581835 commit c49e3c3
Show file tree
Hide file tree
Showing 57 changed files with 1,512 additions and 1,354 deletions.
64 changes: 32 additions & 32 deletions simulations/geometryXYVxVy/landau/landau4d_fft.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -53,30 +53,30 @@ int main(int argc, char** argv)

// Reading config
// --> Mesh info
IDomainX const mesh_x = init_spline_dependent_idx_range<
IDimX,
IdxRangeX const mesh_x = init_spline_dependent_idx_range<
GridX,
BSplinesX,
SplineInterpPointsX>(conf_voicexx, "x");
IDomainY const mesh_y = init_spline_dependent_idx_range<
IDimY,
IdxRangeY const mesh_y = init_spline_dependent_idx_range<
GridY,
BSplinesY,
SplineInterpPointsY>(conf_voicexx, "y");
IDomainVx const mesh_vx = init_spline_dependent_idx_range<
IDimVx,
IdxRangeVx const mesh_vx = init_spline_dependent_idx_range<
GridVx,
BSplinesVx,
SplineInterpPointsVx>(conf_voicexx, "vx");
IDomainVy const mesh_vy = init_spline_dependent_idx_range<
IDimVy,
IdxRangeVy const mesh_vy = init_spline_dependent_idx_range<
GridVy,
BSplinesVy,
SplineInterpPointsVy>(conf_voicexx, "vy");
IDomainXY const mesh_xy(mesh_x, mesh_y);
IDomainVxVy mesh_vxvy(mesh_vx, mesh_vy);
IDomainXYVxVy const meshXYVxVy(mesh_x, mesh_y, mesh_vx, mesh_vy);
IdxRangeXY const mesh_xy(mesh_x, mesh_y);
IdxRangeVxVy mesh_vxvy(mesh_vx, mesh_vy);
IdxRangeXYVxVy const meshXYVxVy(mesh_x, mesh_y, mesh_vx, mesh_vy);

IdxRangeSp const dom_kinsp = init_species(conf_voicexx);

IDomainSpVxVy const meshSpVxVy(dom_kinsp, mesh_vx, mesh_vy);
IDomainSpXYVxVy const meshSpXYVxVy(dom_kinsp, meshXYVxVy);
IdxRangeSpVxVy const meshSpVxVy(dom_kinsp, mesh_vx, mesh_vy);
IdxRangeSpXYVxVy const meshSpXYVxVy(dom_kinsp, meshXYVxVy);

SplineXBuilder const builder_x(meshXYVxVy);
SplineYBuilder const builder_y(meshXYVxVy);
Expand All @@ -86,15 +86,15 @@ int main(int argc, char** argv)
SplineVyBuilder_1d const builder_vy_1d(mesh_vy);

// Initialization of the distribution function
DFieldSpVxVy allfequilibrium(meshSpVxVy);
DFieldMemSpVxVy allfequilibrium(meshSpVxVy);
MaxwellianEquilibrium const init_fequilibrium
= MaxwellianEquilibrium::init_from_input(dom_kinsp, conf_voicexx);
init_fequilibrium(allfequilibrium);
DFieldSpXYVxVy allfdistribu(meshSpXYVxVy);
DFieldMemSpXYVxVy allfdistribu(meshSpXYVxVy);
SingleModePerturbInitialization const init = SingleModePerturbInitialization::
init_from_input(allfequilibrium, dom_kinsp, conf_voicexx);
init(allfdistribu);
auto allfequilibrium_host = ddc::create_mirror_view_and_copy(allfequilibrium.span_view());
auto allfequilibrium_host = ddc::create_mirror_view_and_copy(get_field(allfequilibrium));

// --> Algorithm info
double const deltat = PCpp_double(conf_voicexx, ".Algorithm.deltat");
Expand All @@ -105,44 +105,44 @@ int main(int argc, char** argv)
int const nbstep_diag = int(time_diag / deltat);

// Create spline evaluator
ddc::PeriodicExtrapolationRule<RDimX> bv_x_min;
ddc::PeriodicExtrapolationRule<RDimX> bv_x_max;
ddc::PeriodicExtrapolationRule<X> bv_x_min;
ddc::PeriodicExtrapolationRule<X> bv_x_max;
SplineXEvaluator const spline_x_evaluator(bv_x_min, bv_x_max);

PreallocatableSplineInterpolator const spline_x_interpolator(builder_x, spline_x_evaluator);

ddc::PeriodicExtrapolationRule<RDimY> bv_y_min;
ddc::PeriodicExtrapolationRule<RDimY> bv_y_max;
ddc::PeriodicExtrapolationRule<Y> bv_y_min;
ddc::PeriodicExtrapolationRule<Y> bv_y_max;
SplineYEvaluator const spline_y_evaluator(bv_y_min, bv_y_max);

PreallocatableSplineInterpolator const spline_y_interpolator(builder_y, spline_y_evaluator);

ddc::ConstantExtrapolationRule<RDimVx> bv_vx_min(ddc::coordinate(mesh_vx.front()));
ddc::ConstantExtrapolationRule<RDimVx> bv_vx_max(ddc::coordinate(mesh_vx.back()));
ddc::ConstantExtrapolationRule<Vx> bv_vx_min(ddc::coordinate(mesh_vx.front()));
ddc::ConstantExtrapolationRule<Vx> bv_vx_max(ddc::coordinate(mesh_vx.back()));
SplineVxEvaluator const spline_vx_evaluator(bv_vx_min, bv_vx_max);

PreallocatableSplineInterpolator const spline_vx_interpolator(builder_vx, spline_vx_evaluator);

ddc::ConstantExtrapolationRule<RDimVy> bv_vy_min(ddc::coordinate(mesh_vy.front()));
ddc::ConstantExtrapolationRule<RDimVy> bv_vy_max(ddc::coordinate(mesh_vy.back()));
ddc::ConstantExtrapolationRule<Vy> bv_vy_min(ddc::coordinate(mesh_vy.front()));
ddc::ConstantExtrapolationRule<Vy> bv_vy_max(ddc::coordinate(mesh_vy.back()));
SplineVyEvaluator const spline_vy_evaluator(bv_vy_min, bv_vy_max);

PreallocatableSplineInterpolator const spline_vy_interpolator(builder_vy, spline_vy_evaluator);

// Create advection operator
BslAdvectionSpatial<GeometryXYVxVy, IDimX> const advection_x(spline_x_interpolator);
BslAdvectionSpatial<GeometryXYVxVy, IDimY> const advection_y(spline_y_interpolator);
BslAdvectionVelocity<GeometryXYVxVy, IDimVx> const advection_vx(spline_vx_interpolator);
BslAdvectionVelocity<GeometryXYVxVy, IDimVy> const advection_vy(spline_vy_interpolator);
BslAdvectionSpatial<GeometryXYVxVy, GridX> const advection_x(spline_x_interpolator);
BslAdvectionSpatial<GeometryXYVxVy, GridY> const advection_y(spline_y_interpolator);
BslAdvectionVelocity<GeometryXYVxVy, GridVx> const advection_vx(spline_vx_interpolator);
BslAdvectionVelocity<GeometryXYVxVy, GridVy> const advection_vy(spline_vy_interpolator);

SplitVlasovSolver const vlasov(advection_x, advection_y, advection_vx, advection_vy);

host_t<DFieldVxVy> const quadrature_coeffs_host
host_t<DFieldMemVxVy> const quadrature_coeffs_host
= neumann_spline_quadrature_coefficients(mesh_vxvy, builder_vx_1d, builder_vy_1d);
auto quadrature_coeffs = ddc::create_mirror_view_and_copy(
Kokkos::DefaultExecutionSpace(),
quadrature_coeffs_host.span_view());
FFTPoissonSolver<IDomainXY, IDomainXY, Kokkos::DefaultExecutionSpace> fft_poisson_solver(
get_field(quadrature_coeffs_host));
FFTPoissonSolver<IdxRangeXY, IdxRangeXY, Kokkos::DefaultExecutionSpace> fft_poisson_solver(
mesh_xy);
ChargeDensityCalculator const rhs(quadrature_coeffs);
QNSolver const poisson(fft_poisson_solver, rhs);
Expand All @@ -167,7 +167,7 @@ int main(int argc, char** argv)

steady_clock::time_point const start = steady_clock::now();

predcorr(allfdistribu.span_view(), deltat, nbiter);
predcorr(get_field(allfdistribu), deltat, nbiter);

steady_clock::time_point const end = steady_clock::now();

Expand Down
39 changes: 20 additions & 19 deletions src/geometryRTheta/poisson/polarpoissonlikesolver.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -159,7 +159,7 @@ class PolarSplineFEMPoissonLikeSolver
const int nbasis_p;

// Domains
BSDomainPolar fem_non_singular_domain;
BSDomainPolar fem_non_singular_idx_range;
BSDomainR_Polar radial_bsplines;
BSDomainP_Polar polar_bsplines;

Expand Down Expand Up @@ -216,8 +216,8 @@ class PolarSplineFEMPoissonLikeSolver
Mapping const& mapping)
: nbasis_r(ddc::discrete_space<BSplinesR_Polar>().nbasis() - n_overlap_cells - 1)
, nbasis_p(ddc::discrete_space<BSplinesP_Polar>().nbasis())
, fem_non_singular_domain(
ddc::discrete_space<PolarBSplinesRP>().tensor_bspline_domain().remove_last(
, fem_non_singular_idx_range(
ddc::discrete_space<PolarBSplinesRP>().tensor_bspline_idx_range().remove_last(
ddc::DiscreteVector<PolarBSplinesRP> {nbasis_p}))
, radial_bsplines(ddc::discrete_space<BSplinesR_Polar>().full_domain().remove_first(
ddc::DiscreteVector<BSplinesR_Polar> {n_overlap_cells}))
Expand All @@ -240,7 +240,7 @@ class PolarSplineFEMPoissonLikeSolver
, weights_r(quadrature_domain_r)
, weights_p(quadrature_domain_p)
, singular_basis_vals_and_derivs(ddc::DiscreteDomain<PolarBSplinesRP, QDimRMesh, QDimPMesh>(
PolarBSplinesRP::singular_domain<PolarBSplinesRP>(),
PolarBSplinesRP::singular_idx_range<PolarBSplinesRP>(),
ddc::select<QDimRMesh>(quadrature_domain_singular),
ddc::select<QDimPMesh>(quadrature_domain_singular)))
, r_basis_vals_and_derivs(ddc::DiscreteDomain<
Expand Down Expand Up @@ -316,7 +316,7 @@ class PolarSplineFEMPoissonLikeSolver
}
});

auto singular_domain = PolarBSplinesRP::singular_domain<PolarBSplinesRP>();
auto singular_idx_range = PolarBSplinesRP::singular_idx_range<PolarBSplinesRP>();

// Find value and derivative of 2D bsplines covering the singular point
ddc::for_each(quadrature_domain_singular, [&](QuadratureIndexRP const irp) {
Expand All @@ -335,20 +335,20 @@ class PolarSplineFEMPoissonLikeSolver

// Calculate the value
ddc::discrete_space<PolarBSplinesRP>().eval_basis(singular_vals, vals, coord);
for (auto ib : singular_domain) {
for (auto ib : singular_idx_range) {
singular_basis_vals_and_derivs(ib, ir, ip).value = singular_vals[ib.uid()];
}

// Calculate the radial derivative
ddc::discrete_space<PolarBSplinesRP>().eval_deriv_r(singular_vals, vals, coord);
for (auto ib : singular_domain) {
for (auto ib : singular_idx_range) {
singular_basis_vals_and_derivs(ib, ir, ip).radial_derivative
= singular_vals[ib.uid()];
}

// Calculate the poloidal derivative
ddc::discrete_space<PolarBSplinesRP>().eval_deriv_p(singular_vals, vals, coord);
for (auto ib : singular_domain) {
for (auto ib : singular_idx_range) {
singular_basis_vals_and_derivs(ib, ir, ip).poloidal_derivative
= singular_vals[ib.uid()];
}
Expand Down Expand Up @@ -399,8 +399,8 @@ class PolarSplineFEMPoissonLikeSolver

int matrix_idx(0);
// Calculate the matrix elements corresponding to the bsplines which cover the singular point
ddc::for_each(singular_domain, [&](IndexPolarBspl const idx_test) {
ddc::for_each(singular_domain, [&](IndexPolarBspl const idx_trial) {
ddc::for_each(singular_idx_range, [&](IndexPolarBspl const idx_test) {
ddc::for_each(singular_idx_range, [&](IndexPolarBspl const idx_trial) {
// Calculate the weak integral
row_coo_host(matrix_idx) = idx_test.uid();
col_coo_host(matrix_idx) = idx_trial.uid();
Expand Down Expand Up @@ -431,12 +431,13 @@ class PolarSplineFEMPoissonLikeSolver
BSDomainR central_radial_bspline_domain(radial_bsplines.take_first(
ddc::DiscreteVector<BSplinesR_Polar> {BSplinesR_Polar::degree()}));

BSDomainRP non_singular_domain_near_centre(central_radial_bspline_domain, polar_bsplines);
BSDomainRP
non_singular_idx_range_near_centre(central_radial_bspline_domain, polar_bsplines);

// Calculate the matrix elements where bspline products overlap the bsplines which cover the singular point
ddc::for_each(singular_domain, [&](IndexPolarBspl const idx_test) {
ddc::for_each(singular_idx_range, [&](IndexPolarBspl const idx_test) {
ddc::for_each(
non_singular_domain_near_centre,
non_singular_idx_range_near_centre,
[&](IDimBSpline2D_Polar const idx_trial) {
const IndexPolarBspl polar_idx_trial(
PolarBSplinesRP::get_polar_index<PolarBSplinesRP>(idx_trial));
Expand Down Expand Up @@ -520,7 +521,7 @@ class PolarSplineFEMPoissonLikeSolver
assert(matrix_idx == n_elements_singular + n_elements_overlap);

// Calculate the matrix elements following a stencil
ddc::for_each(fem_non_singular_domain, [&](IndexPolarBspl const polar_idx_test) {
ddc::for_each(fem_non_singular_idx_range, [&](IndexPolarBspl const polar_idx_test) {
const IDimBSpline2D_Polar idx_test(PolarBSplinesRP::get_2d_index(polar_idx_test));
const std::size_t r_idx_test(
ddc::select<PolarBSplinesRP::BSplinesR_tag>(idx_test).uid());
Expand Down Expand Up @@ -640,7 +641,7 @@ class PolarSplineFEMPoissonLikeSolver

// Fill b
ddc::for_each(
PolarBSplinesRP::singular_domain<PolarBSplinesRP>(),
PolarBSplinesRP::singular_idx_range<PolarBSplinesRP>(),
[&](IndexPolarBspl const idx) {
b_host(0, idx.uid()) = ddc::transform_reduce(
quadrature_domain_singular,
Expand All @@ -656,7 +657,7 @@ class PolarSplineFEMPoissonLikeSolver
});
});
const std::size_t ncells_r = ddc::discrete_space<BSplinesR_Polar>().ncells();
ddc::for_each(fem_non_singular_domain, [&](IndexPolarBspl const idx) {
ddc::for_each(fem_non_singular_idx_range, [&](IndexPolarBspl const idx) {
const IDimBSpline2D_Polar idx_2d(PolarBSplinesRP::get_2d_index(idx));
const std::size_t r_idx(ddc::select<PolarBSplinesRP::BSplinesR_tag>(idx_2d).uid());
const std::size_t p_idx(ddc::select<PolarBSplinesRP::BSplinesP_tag>(idx_2d).uid());
Expand Down Expand Up @@ -725,11 +726,11 @@ class PolarSplineFEMPoissonLikeSolver

// Fill the spline
ddc::for_each(
PolarBSplinesRP::singular_domain<PolarBSplinesRP>(),
PolarBSplinesRP::singular_idx_range<PolarBSplinesRP>(),
[&](IndexPolarBspl const idx) {
spline.singular_spline_coef(idx) = b_host(0, idx.uid());
});
ddc::for_each(fem_non_singular_domain, [&](IndexPolarBspl const idx) {
ddc::for_each(fem_non_singular_idx_range, [&](IndexPolarBspl const idx) {
const IDimBSpline2D_Polar idx_2d(PolarBSplinesRP::get_2d_index(idx));
spline.spline_coef(idx_2d) = b_host(0, idx.uid());
});
Expand Down Expand Up @@ -774,7 +775,7 @@ class PolarSplineFEMPoissonLikeSolver
{
BSDomainP_Polar polar_domain(ddc::discrete_space<BSplinesP_Polar>().full_domain());
SplinePolar
spline(PolarBSplinesRP::singular_domain<PolarBSplinesRP>(),
spline(PolarBSplinesRP::singular_idx_range<PolarBSplinesRP>(),
BSDomainRP(radial_bsplines, polar_domain));

(*this)(rhs, spline);
Expand Down
2 changes: 1 addition & 1 deletion src/geometryRTheta/time_solver/bsl_predcorr.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -122,7 +122,7 @@ class BslPredCorrRP : public ITimeSolverRP
BSDomainP polar_domain(ddc::discrete_space<BSplinesP>().full_domain());

SplinePolar electrostatic_potential_coef(
PolarBSplinesRP::singular_domain<PolarBSplinesRP>(),
PolarBSplinesRP::singular_idx_range<PolarBSplinesRP>(),
BSDomainRP(radial_bsplines, polar_domain));
ddc::NullExtrapolationRule extrapolation_rule;
PolarSplineEvaluator<PolarBSplinesRP, ddc::NullExtrapolationRule> polar_spline_evaluator(
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -155,7 +155,7 @@ class BslExplicitPredCorrRP : public ITimeSolverRP
DFieldRP electrical_potential(grid);

SplinePolar electrostatic_potential_coef(
PolarBSplinesRP::singular_domain<PolarBSplinesRP>(),
PolarBSplinesRP::singular_idx_range<PolarBSplinesRP>(),
BSDomainRP(radial_bsplines, polar_domain));

ddc::NullExtrapolationRule extrapolation_rule;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -155,7 +155,7 @@ class BslImplicitPredCorrRP : public ITimeSolverRP
DFieldRP electrical_potential(grid);

SplinePolar electrostatic_potential_coef(
PolarBSplinesRP::singular_domain<PolarBSplinesRP>(),
PolarBSplinesRP::singular_idx_range<PolarBSplinesRP>(),
BSDomainRP(radial_bsplines, polar_domain));

ddc::NullExtrapolationRule extrapolation_rule;
Expand Down
2 changes: 1 addition & 1 deletion src/geometryXVx/boltzmann/splitvlasovsolver.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ class IAdvectionVelocity;
* The Vlasov equation is split between two advection equations
* along the X and Vx directions. The splitting involves solving
* the x-direction advection first on a time interval of length dt/2,
* then the vx-direction advection on a tim dt, and then x-direction
* then the vx-direction advection on a time dt, and then x-direction
* again on dt/2.
*/
class SplitVlasovSolver : public IBoltzmannSolver
Expand Down
Loading

0 comments on commit c49e3c3

Please sign in to comment.