diff --git a/ci_tools/gyselalib_static_analysis.py b/ci_tools/gyselalib_static_analysis.py index bcfcf847..fb7de5c3 100644 --- a/ci_tools/gyselalib_static_analysis.py +++ b/ci_tools/gyselalib_static_analysis.py @@ -602,7 +602,8 @@ def check_exec_space_usage(file): scope = scopes[scope_id] if config.data[current_idx]['str'] == 'KOKKOS_LAMBDA': - ast_parent = config.data_xml.find(f".token[@id='{config.data[current_idx-1]['astParent']}']").attrib + ast_parent_id = config.data[current_idx-1].get('astParent', config.data[current_idx-1]['id']) + ast_parent = config.data_xml.find(f".token[@id='{ast_parent_id}']").attrib parallel_func_start_idx = config.data.index(ast_parent) args = [t['str'] for t in config.data[parallel_func_start_idx+1:current_idx]] if 'DefaultHostExecutionSpace' in args: diff --git a/simulations/geometryRTheta/diocotron/diocotron.cpp b/simulations/geometryRTheta/diocotron/diocotron.cpp index b5ee30df..060fe35c 100644 --- a/simulations/geometryRTheta/diocotron/diocotron.cpp +++ b/simulations/geometryRTheta/diocotron/diocotron.cpp @@ -49,8 +49,10 @@ using PoissonSolver = PolarSplineFEMPoissonLikeSolver< GridR, GridTheta, PolarBSplinesRTheta, - SplineRThetaEvaluatorNullBound_host>; -using DiscreteMappingBuilder = DiscreteToCartesianBuilder< + SplineRThetaEvaluatorNullBound>; +using DiscreteMappingBuilder + = DiscreteToCartesianBuilder; +using DiscreteMappingBuilder_host = DiscreteToCartesianBuilder< X, Y, SplineRThetaBuilder_host, @@ -104,7 +106,9 @@ int main(int argc, char** argv) // OPERATORS ====================================================================================== - SplineRThetaBuilder_host const builder(mesh_rp); + SplineRThetaBuilder const builder(mesh_rp); + SplineRThetaBuilder_host const builder_host(mesh_rp); + // --- Define the mapping. ------------------------------------------------------------------------ ddc::ConstantExtrapolationRule boundary_condition_r_left( @@ -112,7 +116,13 @@ int main(int argc, char** argv) ddc::ConstantExtrapolationRule boundary_condition_r_right( ddc::coordinate(mesh_r.back())); - SplineRThetaEvaluatorConstBound_host spline_evaluator_extrapol( + SplineRThetaEvaluatorConstBound spline_evaluator_extrapol( + boundary_condition_r_left, + boundary_condition_r_right, + ddc::PeriodicExtrapolationRule(), + ddc::PeriodicExtrapolationRule()); + + SplineRThetaEvaluatorConstBound_host spline_evaluator_extrapol_host( boundary_condition_r_left, boundary_condition_r_right, ddc::PeriodicExtrapolationRule(), @@ -120,15 +130,22 @@ int main(int argc, char** argv) const LogicalToPhysicalMapping to_physical_mapping; DiscreteMappingBuilder const discrete_mapping_builder( - Kokkos::DefaultHostExecutionSpace(), + Kokkos::DefaultExecutionSpace(), to_physical_mapping, builder, spline_evaluator_extrapol); + DiscreteMappingBuilder_host const discrete_mapping_builder_host( + Kokkos::DefaultHostExecutionSpace(), + to_physical_mapping, + builder_host, + spline_evaluator_extrapol_host); DiscreteToCartesian const discrete_mapping = discrete_mapping_builder(); + DiscreteToCartesian const discrete_mapping_host = discrete_mapping_builder_host(); - ddc::init_discrete_space(discrete_mapping); - IdxRangeBSRTheta const idx_range_bsplinesRTheta = get_spline_idx_range(builder); + ddc::init_discrete_space(discrete_mapping_host); + + IdxRangeBSRTheta const idx_range_bsplinesRTheta = get_spline_idx_range(builder_host); // --- Time integration method -------------------------------------------------------------------- @@ -160,13 +177,18 @@ int main(int argc, char** argv) // --- Advection operator ------------------------------------------------------------------------- ddc::NullExtrapolationRule r_extrapolation_rule; ddc::PeriodicExtrapolationRule p_extrapolation_rule; - SplineRThetaEvaluatorNullBound_host spline_evaluator( + SplineRThetaEvaluatorNullBound spline_evaluator( + r_extrapolation_rule, + r_extrapolation_rule, + p_extrapolation_rule, + p_extrapolation_rule); + SplineRThetaEvaluatorNullBound_host spline_evaluator_host( r_extrapolation_rule, r_extrapolation_rule, p_extrapolation_rule, p_extrapolation_rule); - PreallocatableSplineInterpolatorRTheta interpolator(builder, spline_evaluator); + PreallocatableSplineInterpolatorRTheta interpolator(builder_host, spline_evaluator_host); AdvectionDomain advection_domain(to_physical_mapping); @@ -174,8 +196,8 @@ int main(int argc, char** argv) time_stepper, advection_domain, to_physical_mapping, - builder, - spline_evaluator_extrapol); + builder_host, + spline_evaluator_extrapol_host); BslAdvectionRTheta advection_operator(interpolator, find_feet, to_physical_mapping); @@ -183,16 +205,14 @@ int main(int argc, char** argv) // --- Poisson solver ----------------------------------------------------------------------------- // Coefficients alpha and beta of the Poisson equation: - host_t coeff_alpha(mesh_rp); - host_t coeff_beta(mesh_rp); + DFieldMemRTheta coeff_alpha(mesh_rp); + DFieldMemRTheta coeff_beta(mesh_rp); + ddc::parallel_fill(coeff_alpha, -1); + ddc::parallel_fill(coeff_beta, 0); - ddc::for_each(mesh_rp, [&](IdxRTheta const irp) { - coeff_alpha(irp) = -1.0; - coeff_beta(irp) = 0.0; - }); - host_t coeff_alpha_spline(idx_range_bsplinesRTheta); - host_t coeff_beta_spline(idx_range_bsplinesRTheta); + Spline2DMem coeff_alpha_spline(idx_range_bsplinesRTheta); + Spline2DMem coeff_beta_spline(idx_range_bsplinesRTheta); builder(get_field(coeff_alpha_spline), get_const_field(coeff_alpha)); builder(get_field(coeff_beta_spline), get_const_field(coeff_beta)); @@ -208,8 +228,8 @@ int main(int argc, char** argv) BslPredCorrRTheta predcorr_operator( to_physical_mapping, advection_operator, - builder, - spline_evaluator, + builder_host, + spline_evaluator_host, poisson_solver); #elif defined(EXPLICIT_PREDCORR) BslExplicitPredCorrRTheta predcorr_operator( @@ -217,20 +237,20 @@ int main(int argc, char** argv) to_physical_mapping, advection_operator, mesh_rp, - builder, - spline_evaluator, + builder_host, + spline_evaluator_host, poisson_solver, - spline_evaluator_extrapol); + spline_evaluator_extrapol_host); #elif defined(IMPLICIT_PREDCORR) BslImplicitPredCorrRTheta predcorr_operator( advection_domain, to_physical_mapping, advection_operator, mesh_rp, - builder, - spline_evaluator, + builder_host, + spline_evaluator_host, poisson_solver, - spline_evaluator_extrapol); + spline_evaluator_extrapol_host); #endif // ================================================================================================ @@ -296,12 +316,13 @@ int main(int argc, char** argv) }); // Compute phi equilibrium phi_eq from Poisson solver. *********** - host_t phi_eq(mesh_rp); + DFieldMemRTheta phi_eq(mesh_rp); + host_t phi_eq_host(mesh_rp); host_t rho_coef_eq(idx_range_bsplinesRTheta); - builder(get_field(rho_coef_eq), get_const_field(rho_eq)); - PoissonLikeRHSFunction poisson_rhs_eq(get_const_field(rho_coef_eq), spline_evaluator); + builder_host(get_field(rho_coef_eq), get_const_field(rho_eq)); + PoissonLikeRHSFunction poisson_rhs_eq(get_const_field(rho_coef_eq), spline_evaluator_host); poisson_solver(poisson_rhs_eq, get_field(phi_eq)); - + ddc::parallel_deepcopy(phi_eq, phi_eq_host); // --- Save initial data -------------------------------------------------------------------------- ddc::PdiEvent("initialization") @@ -309,7 +330,7 @@ int main(int argc, char** argv) .and_with("y_coords", coords_y) .and_with("jacobian", jacobian) .and_with("density_eq", rho_eq) - .and_with("electrical_potential_eq", phi_eq); + .and_with("electrical_potential_eq", phi_eq_host); // ================================================================================================ diff --git a/simulations/geometryRTheta/vortex_merger/vortex_merger.cpp b/simulations/geometryRTheta/vortex_merger/vortex_merger.cpp index a01cdd02..3e206793 100644 --- a/simulations/geometryRTheta/vortex_merger/vortex_merger.cpp +++ b/simulations/geometryRTheta/vortex_merger/vortex_merger.cpp @@ -49,8 +49,10 @@ using PoissonSolver = PolarSplineFEMPoissonLikeSolver< GridR, GridTheta, PolarBSplinesRTheta, - SplineRThetaEvaluatorNullBound_host>; -using DiscreteMappingBuilder = DiscreteToCartesianBuilder< + SplineRThetaEvaluatorNullBound>; +using DiscreteMappingBuilder + = DiscreteToCartesianBuilder; +using DiscreteMappingBuilder_host = DiscreteToCartesianBuilder< X, Y, SplineRThetaBuilder_host, @@ -99,7 +101,8 @@ int main(int argc, char** argv) // OPERATORS ====================================================================================== - SplineRThetaBuilder_host const builder(grid); + SplineRThetaBuilder const builder(grid); + SplineRThetaBuilder_host const builder_host(grid); // --- Define the mapping. ------------------------------------------------------------------------ ddc::ConstantExtrapolationRule boundary_condition_r_left( @@ -107,7 +110,12 @@ int main(int argc, char** argv) ddc::ConstantExtrapolationRule boundary_condition_r_right( ddc::coordinate(mesh_r.back())); - SplineRThetaEvaluatorConstBound_host spline_evaluator_extrapol( + SplineRThetaEvaluatorConstBound spline_evaluator_extrapol( + boundary_condition_r_left, + boundary_condition_r_right, + ddc::PeriodicExtrapolationRule(), + ddc::PeriodicExtrapolationRule()); + SplineRThetaEvaluatorConstBound_host spline_evaluator_extrapol_host( boundary_condition_r_left, boundary_condition_r_right, ddc::PeriodicExtrapolationRule(), @@ -115,15 +123,21 @@ int main(int argc, char** argv) const LogicalToPhysicalMapping to_physical_mapping; DiscreteMappingBuilder const discrete_mapping_builder( - Kokkos::DefaultHostExecutionSpace(), + Kokkos::DefaultExecutionSpace(), to_physical_mapping, builder, spline_evaluator_extrapol); + DiscreteMappingBuilder_host const discrete_mapping_builder_host( + Kokkos::DefaultHostExecutionSpace(), + to_physical_mapping, + builder_host, + spline_evaluator_extrapol_host); DiscreteToCartesian const discrete_mapping = discrete_mapping_builder(); + DiscreteToCartesian const discrete_mapping_host = discrete_mapping_builder_host(); - ddc::init_discrete_space(discrete_mapping); + ddc::init_discrete_space(discrete_mapping_host); - IdxRangeBSRTheta const idx_range_bsplinesRTheta = get_spline_idx_range(builder); + IdxRangeBSRTheta const idx_range_bsplinesRTheta = get_spline_idx_range(builder_host); // --- Time integration method -------------------------------------------------------------------- @@ -135,13 +149,18 @@ int main(int argc, char** argv) // --- Advection operator ------------------------------------------------------------------------- ddc::NullExtrapolationRule r_extrapolation_rule; ddc::PeriodicExtrapolationRule p_extrapolation_rule; - SplineRThetaEvaluatorNullBound_host spline_evaluator( + SplineRThetaEvaluatorNullBound spline_evaluator( + r_extrapolation_rule, + r_extrapolation_rule, + p_extrapolation_rule, + p_extrapolation_rule); + SplineRThetaEvaluatorNullBound_host spline_evaluator_host( r_extrapolation_rule, r_extrapolation_rule, p_extrapolation_rule, p_extrapolation_rule); - PreallocatableSplineInterpolatorRTheta interpolator(builder, spline_evaluator); + PreallocatableSplineInterpolatorRTheta interpolator(builder_host, spline_evaluator_host); AdvectionDomain advection_domain(to_physical_mapping); @@ -149,8 +168,8 @@ int main(int argc, char** argv) time_stepper, advection_domain, to_physical_mapping, - builder, - spline_evaluator_extrapol); + builder_host, + spline_evaluator_extrapol_host); BslAdvectionRTheta advection_operator(interpolator, find_feet, to_physical_mapping); @@ -158,16 +177,14 @@ int main(int argc, char** argv) // --- Poisson solver ----------------------------------------------------------------------------- // Coefficients alpha and beta of the Poisson equation: - host_t coeff_alpha(grid); - host_t coeff_beta(grid); + DFieldMemRTheta coeff_alpha(grid); + DFieldMemRTheta coeff_beta(grid); - ddc::for_each(grid, [&](IdxRTheta const irp) { - coeff_alpha(irp) = -1.0; - coeff_beta(irp) = 0.0; - }); + ddc::parallel_fill(coeff_alpha, -1); + ddc::parallel_fill(coeff_beta, 0); - host_t coeff_alpha_spline(idx_range_bsplinesRTheta); - host_t coeff_beta_spline(idx_range_bsplinesRTheta); + Spline2DMem coeff_alpha_spline(idx_range_bsplinesRTheta); + Spline2DMem coeff_beta_spline(idx_range_bsplinesRTheta); builder(get_field(coeff_alpha_spline), get_const_field(coeff_alpha)); builder(get_field(coeff_beta_spline), get_const_field(coeff_beta)); @@ -184,10 +201,10 @@ int main(int argc, char** argv) to_physical_mapping, advection_operator, grid, - builder, - spline_evaluator, + builder_host, + spline_evaluator_host, poisson_solver, - spline_evaluator_extrapol); + spline_evaluator_extrapol_host); @@ -243,8 +260,12 @@ int main(int argc, char** argv) double const phi_max(1.); - VortexMergerEquilibria - equilibrium(to_physical_mapping, grid, builder, spline_evaluator, poisson_solver); + VortexMergerEquilibria equilibrium( + to_physical_mapping, + grid, + builder_host, + spline_evaluator_host, + poisson_solver); std::function const function = [&](double const x) { return x * x; }; host_t rho_eq(grid); equilibrium.set_equilibrium(get_field(rho_eq), function, phi_max, tau); @@ -264,11 +285,13 @@ int main(int argc, char** argv) // Compute phi equilibrium phi_eq from Poisson solver. *********** - host_t phi_eq(grid); + DFieldMemRTheta phi_eq(grid); + host_t phi_eq_host(grid); host_t rho_coef_eq(idx_range_bsplinesRTheta); - builder(get_field(rho_coef_eq), get_const_field(rho_eq)); - PoissonLikeRHSFunction poisson_rhs_eq(get_const_field(rho_coef_eq), spline_evaluator); + builder_host(get_field(rho_coef_eq), get_const_field(rho_eq)); + PoissonLikeRHSFunction poisson_rhs_eq(get_const_field(rho_coef_eq), spline_evaluator_host); poisson_solver(poisson_rhs_eq, get_field(phi_eq)); + ddc::parallel_deepcopy(phi_eq, phi_eq_host); // --- Save initial data -------------------------------------------------------------------------- diff --git a/src/geometryRTheta/geometry/geometry.hpp b/src/geometryRTheta/geometry/geometry.hpp index 1aa9b812..589b1123 100644 --- a/src/geometryRTheta/geometry/geometry.hpp +++ b/src/geometryRTheta/geometry/geometry.hpp @@ -163,7 +163,48 @@ using SplineRThetaEvaluatorNullBound_host = ddc::SplineEvaluator2D< GridR, GridTheta>; +using SplineRThetaBuilder = ddc::SplineBuilder2D< + Kokkos::DefaultExecutionSpace, + typename Kokkos::DefaultExecutionSpace::memory_space, + BSplinesR, + BSplinesTheta, + GridR, + GridTheta, + SplineRBoundary, // boundary at r=0 + SplineRBoundary, // boundary at rmax + SplinePBoundary, + SplinePBoundary, + ddc::SplineSolver::LAPACK, + GridR, + GridTheta>; + +using SplineRThetaEvaluatorConstBound = ddc::SplineEvaluator2D< + Kokkos::DefaultExecutionSpace, + typename Kokkos::DefaultExecutionSpace::memory_space, + BSplinesR, + BSplinesTheta, + GridR, + GridTheta, + ddc::ConstantExtrapolationRule, // boundary at r=0 + ddc::ConstantExtrapolationRule, // boundary at rmax + ddc::PeriodicExtrapolationRule, + ddc::PeriodicExtrapolationRule, + GridR, + GridTheta>; +using SplineRThetaEvaluatorNullBound = ddc::SplineEvaluator2D< + Kokkos::DefaultExecutionSpace, + typename Kokkos::DefaultExecutionSpace::memory_space, + BSplinesR, + BSplinesTheta, + GridR, + GridTheta, + ddc::NullExtrapolationRule, // boundary at r=0 + ddc::NullExtrapolationRule, // boundary at rmax + ddc::PeriodicExtrapolationRule, + ddc::PeriodicExtrapolationRule, + GridR, + GridTheta>; // --- Index definitions using IdxR = Idx; using IdxTheta = Idx; diff --git a/src/geometryRTheta/initialization/vortex_merger_equilibrium.hpp b/src/geometryRTheta/initialization/vortex_merger_equilibrium.hpp index 55a0300d..42d81c06 100644 --- a/src/geometryRTheta/initialization/vortex_merger_equilibrium.hpp +++ b/src/geometryRTheta/initialization/vortex_merger_equilibrium.hpp @@ -31,7 +31,7 @@ class VortexMergerEquilibria GridR, GridTheta, PolarBSplinesRTheta, - SplineRThetaEvaluatorNullBound_host> const& m_poisson_solver; + SplineRThetaEvaluatorNullBound> const& m_poisson_solver; public: /** @@ -60,7 +60,7 @@ class VortexMergerEquilibria GridR, GridTheta, PolarBSplinesRTheta, - SplineRThetaEvaluatorNullBound_host> const& poisson_solver) + SplineRThetaEvaluatorNullBound> const& poisson_solver) : m_mapping(mapping) , m_grid(grid) , m_builder(builder) @@ -110,7 +110,7 @@ class VortexMergerEquilibria double const tau, int count_max = 25) const { - host_t phi_star(m_grid); + DFieldMemRTheta phi_star(m_grid); host_t ci(m_grid); IdxRangeBSRTheta idx_range_bsplinesRTheta = get_spline_idx_range(m_builder); @@ -131,12 +131,13 @@ class VortexMergerEquilibria m_builder(get_field(rho_coef), get_const_field(rho_eq)); PoissonLikeRHSFunction poisson_rhs(get_const_field(rho_coef), m_evaluator); m_poisson_solver(poisson_rhs, get_field(phi_star)); + auto phi_star_host = ddc::create_mirror_view_and_copy(get_field(phi_star)); // STEP 3: compute c^i // If phi_max is given: double norm_Linf_phi_star(0.); ddc::for_each(m_grid, [&](IdxRTheta const irp) { - double const abs_phi_star = fabs(phi_star(irp)); + double const abs_phi_star = fabs(phi_star_host(irp)); norm_Linf_phi_star = norm_Linf_phi_star > abs_phi_star ? norm_Linf_phi_star : abs_phi_star; }); @@ -154,7 +155,7 @@ class VortexMergerEquilibria = difference_sigma > abs_diff_sigma ? difference_sigma : abs_diff_sigma; sigma(irp) = ci(irp) * sigma(irp); - phi_eq(irp) = ci(irp) * phi_star(irp); + phi_eq(irp) = ci(irp) * phi_star_host(irp); }); } while ((difference_sigma > tau) and (count < count_max)); diff --git a/src/geometryRTheta/time_solver/bsl_predcorr.hpp b/src/geometryRTheta/time_solver/bsl_predcorr.hpp index 04585148..7681639c 100644 --- a/src/geometryRTheta/time_solver/bsl_predcorr.hpp +++ b/src/geometryRTheta/time_solver/bsl_predcorr.hpp @@ -66,7 +66,7 @@ class BslPredCorrRTheta : public ITimeSolverRTheta GridR, GridTheta, PolarBSplinesRTheta, - SplineRThetaEvaluatorNullBound_host> const& m_poisson_solver; + SplineRThetaEvaluatorNullBound> const& m_poisson_solver; SplineRThetaBuilder_host const& m_builder; SplineRThetaEvaluatorNullBound_host const& m_spline_evaluator; @@ -99,7 +99,7 @@ class BslPredCorrRTheta : public ITimeSolverRTheta GridR, GridTheta, PolarBSplinesRTheta, - SplineRThetaEvaluatorNullBound_host> const& poisson_solver) + SplineRThetaEvaluatorNullBound> const& poisson_solver) : m_mapping(mapping) , m_advection_solver(advection_solver) , m_poisson_solver(poisson_solver) @@ -134,19 +134,19 @@ class BslPredCorrRTheta : public ITimeSolverRTheta polar_spline_evaluator(extrapolation_rule); - host_t electrical_potential0(grid); - + host_t electrical_potential0_host(grid); + DFieldMemRTheta electrical_potential0(grid); host_t allfdistribu_coef(get_spline_idx_range(m_builder)); m_builder(get_field(allfdistribu_coef), get_const_field(allfdistribu)); PoissonLikeRHSFunction const charge_density_coord(get_const_field(allfdistribu_coef), m_spline_evaluator); m_poisson_solver(charge_density_coord, get_field(electrical_potential0)); - + ddc::parallel_deepcopy(electrical_potential0, electrical_potential0_host); ddc::PdiEvent("iteration") .with("iter", 0) .and_with("time", 0) .and_with("density", allfdistribu) - .and_with("electrical_potential", electrical_potential0); + .and_with("electrical_potential", electrical_potential0_host); std::function>, host_t)> @@ -174,7 +174,8 @@ class BslPredCorrRTheta : public ITimeSolverRTheta host_t>, Kokkos::DefaultHostExecutionSpace> time_stepper(grid); - host_t electrical_potential(grid); + DFieldMemRTheta electrical_potential(grid); + host_t electrical_potential_host(grid); start_time = std::chrono::system_clock::now(); for (int iter(0); iter < steps; ++iter) { time_stepper @@ -189,12 +190,12 @@ class BslPredCorrRTheta : public ITimeSolverRTheta PoissonLikeRHSFunction const charge_density_coord(get_const_field(allfdistribu_coef), m_spline_evaluator); m_poisson_solver(charge_density_coord, get_field(electrical_potential)); - + ddc::parallel_deepcopy(electrical_potential_host, electrical_potential); ddc::PdiEvent("iteration") .with("iter", iter + 1) .and_with("time", iter * dt) .and_with("density", allfdistribu) - .and_with("electrical_potential", electrical_potential); + .and_with("electrical_potential", electrical_potential_host); } end_time = std::chrono::system_clock::now(); diff --git a/src/geometryRTheta/time_solver/bsl_predcorr_second_order_explicit.hpp b/src/geometryRTheta/time_solver/bsl_predcorr_second_order_explicit.hpp index 3fdfac03..99814566 100644 --- a/src/geometryRTheta/time_solver/bsl_predcorr_second_order_explicit.hpp +++ b/src/geometryRTheta/time_solver/bsl_predcorr_second_order_explicit.hpp @@ -85,7 +85,7 @@ class BslExplicitPredCorrRTheta : public ITimeSolverRTheta GridR, GridTheta, PolarBSplinesRTheta, - SplineRThetaEvaluatorNullBound_host> const& m_poisson_solver; + SplineRThetaEvaluatorNullBound> const& m_poisson_solver; SplineRThetaBuilder_host const& m_builder; SplineRThetaEvaluatorConstBound_host const& m_evaluator; @@ -128,7 +128,7 @@ class BslExplicitPredCorrRTheta : public ITimeSolverRTheta GridR, GridTheta, PolarBSplinesRTheta, - SplineRThetaEvaluatorNullBound_host> const& poisson_solver, + SplineRThetaEvaluatorNullBound> const& poisson_solver, SplineRThetaEvaluatorConstBound_host const& advection_evaluator) : m_mapping(mapping) , m_advection_solver(advection_solver) @@ -162,7 +162,8 @@ class BslExplicitPredCorrRTheta : public ITimeSolverRTheta IdxRangeBSTheta polar_idx_range(ddc::discrete_space().full_domain()); // --- Electrostatic potential (phi). ------------------------------------------------------------- - host_t electrical_potential(grid); + DFieldMemRTheta electrical_potential(grid); + host_t electrical_potential_host(grid); host_t electrostatic_potential_coef( PolarBSplinesRTheta::singular_idx_range(), @@ -199,7 +200,7 @@ class BslExplicitPredCorrRTheta : public ITimeSolverRTheta m_poisson_solver(charge_density_coord_1, electrostatic_potential_coef); polar_spline_evaluator( - get_field(electrical_potential), + get_field(electrical_potential_host), get_const_field(coords), electrostatic_potential_coef); @@ -207,7 +208,7 @@ class BslExplicitPredCorrRTheta : public ITimeSolverRTheta .with("iter", iter) .and_with("time", time) .and_with("density", allfdistribu) - .and_with("electrical_potential", electrical_potential); + .and_with("electrical_potential", electrical_potential_host); // STEP 2: From phi^n, we compute A^n: advection_field_computer(electrostatic_potential_coef, advection_field); @@ -281,12 +282,12 @@ class BslExplicitPredCorrRTheta : public ITimeSolverRTheta PoissonLikeRHSFunction const charge_density_coord(get_const_field(allfdistribu_coef), m_evaluator); m_poisson_solver(charge_density_coord, get_field(electrical_potential)); - + ddc::parallel_deepcopy(electrical_potential_host, electrical_potential); ddc::PdiEvent("last_iteration") .with("iter", steps) .and_with("time", steps * dt) .and_with("density", allfdistribu) - .and_with("electrical_potential", electrical_potential); + .and_with("electrical_potential", electrical_potential_host); end_time = std::chrono::system_clock::now(); diff --git a/src/geometryRTheta/time_solver/bsl_predcorr_second_order_implicit.hpp b/src/geometryRTheta/time_solver/bsl_predcorr_second_order_implicit.hpp index 23ca5d0b..d9458be5 100644 --- a/src/geometryRTheta/time_solver/bsl_predcorr_second_order_implicit.hpp +++ b/src/geometryRTheta/time_solver/bsl_predcorr_second_order_implicit.hpp @@ -84,7 +84,7 @@ class BslImplicitPredCorrRTheta : public ITimeSolverRTheta GridR, GridTheta, PolarBSplinesRTheta, - SplineRThetaEvaluatorNullBound_host> const& m_poisson_solver; + SplineRThetaEvaluatorNullBound> const& m_poisson_solver; SplineRThetaBuilder_host const& m_builder; SplineRThetaEvaluatorConstBound_host const& m_evaluator; @@ -127,7 +127,7 @@ class BslImplicitPredCorrRTheta : public ITimeSolverRTheta GridR, GridTheta, PolarBSplinesRTheta, - SplineRThetaEvaluatorNullBound_host> const& poisson_solver, + SplineRThetaEvaluatorNullBound> const& poisson_solver, SplineRThetaEvaluatorConstBound_host const& advection_evaluator) : m_mapping(mapping) , m_advection_solver(advection_solver) @@ -161,7 +161,8 @@ class BslImplicitPredCorrRTheta : public ITimeSolverRTheta IdxRangeBSTheta polar_idx_range(ddc::discrete_space().full_domain()); // --- Electrostatic potential (phi). ------------------------------------------------------------- - host_t electrical_potential(grid); + DFieldMemRTheta electrical_potential(grid); + host_t electrical_potential_host(grid); host_t electrostatic_potential_coef( PolarBSplinesRTheta::singular_idx_range(), @@ -189,7 +190,7 @@ class BslImplicitPredCorrRTheta : public ITimeSolverRTheta m_poisson_solver(charge_density_coord_1, electrostatic_potential_coef); polar_spline_evaluator( - get_field(electrical_potential), + get_field(electrical_potential_host), get_const_field(coords), electrostatic_potential_coef); @@ -197,7 +198,7 @@ class BslImplicitPredCorrRTheta : public ITimeSolverRTheta .with("iter", iter) .and_with("time", iter * dt) .and_with("density", allfdistribu) - .and_with("electrical_potential", electrical_potential); + .and_with("electrical_potential", electrical_potential_host); // STEP 2: From phi^n, we compute A^n: @@ -335,12 +336,13 @@ class BslImplicitPredCorrRTheta : public ITimeSolverRTheta PoissonLikeRHSFunction const charge_density_coord(get_const_field(allfdistribu_coef), m_evaluator); m_poisson_solver(charge_density_coord, get_field(electrical_potential)); + ddc::parallel_deepcopy(electrical_potential_host, electrical_potential); ddc::PdiEvent("last_iteration") .with("iter", steps) .and_with("time", steps * dt) .and_with("density", allfdistribu) - .and_with("electrical_potential", electrical_potential); + .and_with("electrical_potential", electrical_potential_host); end_time = std::chrono::system_clock::now(); display_time_difference("Iterations time: ", start_time, end_time); diff --git a/src/pde_solvers/polarpoissonlikesolver.hpp b/src/pde_solvers/polarpoissonlikesolver.hpp index c42a7e49..27789b0a 100644 --- a/src/pde_solvers/polarpoissonlikesolver.hpp +++ b/src/pde_solvers/polarpoissonlikesolver.hpp @@ -15,6 +15,7 @@ #include "ddc_aliases.hpp" #include "matrix_batch_csr.hpp" + /** * @brief Define a polar PDE solver for a Poisson-like equation. * @@ -36,15 +37,15 @@ * @tparam GridR The radial grid type. * @tparam GridR The poloidal grid type. * @tparam PolarBSplinesRTheta The type of the 2D polar bsplines (on the coordinate - * system (r,theta) including bsplines which traverse the O point). - * @tparam SplineRThetaEvaluatorNullBound_host The type of the 2D (cross-product) spline evaluator. + * system @f$(r,\theta)@f$ including bsplines which traverse the O point). + * @tparam SplineRThetaEvaluatorNullBound The type of the 2D (cross-product) spline evaluator. * @tparam IdxRangeFull The full index range of @f$ \phi @f$ including any batch dimensions. */ template < class GridR, class GridTheta, class PolarBSplinesRTheta, - class SplineRThetaEvaluatorNullBound_host, + class SplineRThetaEvaluatorNullBound, class IdxRangeFull = IdxRange> class PolarSplineFEMPoissonLikeSolver { @@ -90,7 +91,6 @@ class PolarSplineFEMPoissonLikeSolver private: using CoordRTheta = Coord; - /// The 1D bsplines in the radial direction using BSplinesR = typename PolarBSplinesRTheta::BSplinesR_tag; /// The 1D bsplines in the poloidal direction @@ -159,19 +159,21 @@ class PolarSplineFEMPoissonLikeSolver using CoordFieldRTheta = Field; using DFieldRTheta = DField; +public: /** - * @brief Object storing a value and a value of the derivative - * of a 1D function. - */ + * @brief Object storing a value and a value of the derivative + * of a 1D function. + */ struct EvalDeriv1DType { double value; double derivative; }; + /** - * @brief Object storing a value and a value of the derivatives - * in each direction of a 2D function. - */ + * @brief Object storing a value and a value of the derivatives + * in each direction of a 2D function. + */ struct EvalDeriv2DType { double value; @@ -207,6 +209,9 @@ class PolarSplineFEMPoissonLikeSolver const int m_nbasis_r; const int m_nbasis_theta; + // Matrix size is equal to the number of Polar bspline + const int m_matrix_size; + // Domains IdxRangeBSPolar m_idxrange_fem_non_singular; IdxRangeBSR m_idxrange_bsplines_r; @@ -225,11 +230,11 @@ class PolarSplineFEMPoissonLikeSolver // Basis Spline values and derivatives at Gauss-Legendre points host_t>> m_singular_basis_vals_and_derivs; - host_t>> r_basis_vals_and_derivs; + host_t>> m_r_basis_vals_and_derivs; host_t>> m_theta_basis_vals_and_derivs; - host_t> m_int_volume; + FieldMem m_int_volume; PolarSplineEvaluator m_polar_spline_evaluator; @@ -239,7 +244,6 @@ class PolarSplineFEMPoissonLikeSolver Kokkos::View m_x_init; const int m_batch_idx {0}; // TODO: Remove when batching is supported - public: /** * @brief Instantiate a polar Poisson-like solver using FEM with B-splines. @@ -260,18 +264,19 @@ class PolarSplineFEMPoissonLikeSolver * The mapping from the logical index range to the physical index range where * the equation is defined. * @param[in] spline_evaluator - * An evaluator for evaluating 2D splines on (r, theta) + * An evaluator for evaluating 2D splines on @f$(r,\theta)@f$. * * @tparam Mapping A class describing a mapping from curvilinear coordinates to cartesian coordinates. */ template PolarSplineFEMPoissonLikeSolver( - host_t coeff_alpha, - host_t coeff_beta, + ConstSpline2D coeff_alpha, + ConstSpline2D coeff_beta, Mapping const& mapping, - SplineRThetaEvaluatorNullBound_host const& spline_evaluator) + SplineRThetaEvaluatorNullBound const& spline_evaluator) : m_nbasis_r(ddc::discrete_space().nbasis() - m_n_overlap_cells - 1) , m_nbasis_theta(ddc::discrete_space().nbasis()) + , m_matrix_size(ddc::discrete_space().nbasis() - m_nbasis_theta) , m_idxrange_fem_non_singular( ddc::discrete_space().tensor_bspline_idx_range().remove_last( IdxStep {m_nbasis_theta})) @@ -299,7 +304,7 @@ class PolarSplineFEMPoissonLikeSolver PolarBSplinesRTheta::template singular_idx_range(), ddc::select(m_idxrange_quadrature_singular), ddc::select(m_idxrange_quadrature_singular))) - , r_basis_vals_and_derivs( + , m_r_basis_vals_and_derivs( IdxRange(m_non_zero_bases_r, m_idxrange_quadrature_r)) , m_theta_basis_vals_and_derivs( IdxRange< @@ -367,8 +372,8 @@ class PolarSplineFEMPoissonLikeSolver .eval_basis_and_n_derivs(vals, ddc::coordinate(idx_r), 1); for (auto ib : m_non_zero_bases_r) { const int ib_idx = ib - m_non_zero_bases_r.front(); - r_basis_vals_and_derivs(ib, idx_r).value = vals(ib_idx, 0); - r_basis_vals_and_derivs(ib, idx_r).derivative = vals(ib_idx, 1); + m_r_basis_vals_and_derivs(ib, idx_r).value = vals(ib_idx, 0); + m_r_basis_vals_and_derivs(ib, idx_r).derivative = vals(ib_idx, 1); } }); @@ -398,8 +403,8 @@ class PolarSplineFEMPoissonLikeSolver // Values of the polar basis splines, that do not cover the singular point, // at a given coordinate DSpan2D vals(data.data(), m_n_non_zero_bases_r, m_n_non_zero_bases_theta); - IdxQuadratureR idx_r = ddc::select(irp); - IdxQuadratureTheta idx_theta = ddc::select(irp); + IdxQuadratureR const idx_r(irp); + IdxQuadratureTheta const idx_theta(irp); const CoordRTheta coord(ddc::coordinate(irp)); @@ -425,17 +430,6 @@ class PolarSplineFEMPoissonLikeSolver } }); - // Find the integral volume associated with each point used in the quadrature scheme - IdxRangeQuadratureRTheta - all_quad_points(m_idxrange_quadrature_r, m_idxrange_quadrature_theta); - ddc::for_each(all_quad_points, [&](IdxQuadratureRTheta const irp) { - IdxQuadratureR const idx_r = ddc::select(irp); - IdxQuadratureTheta const idx_theta = ddc::select(irp); - CoordRTheta coord(ddc::coordinate(irp)); - m_int_volume(idx_r, idx_theta) = abs(mapping.jacobian(coord)) * m_weights_r(idx_r) - * m_weights_theta(idx_theta); - }); - // Number of elements in the matrix that correspond to the splines // that cover the singular point constexpr int n_elements_singular @@ -454,9 +448,7 @@ class PolarSplineFEMPoissonLikeSolver const int n_elements_stencil = n_stencil_r * n_stencil_theta; const int batch_size = 1; - // Matrix size is equal to the number Polar bspline - const int matrix_size - = ddc::discrete_space().nbasis() - m_nbasis_theta; + const int n_matrix_elements = n_elements_singular + n_elements_overlap + n_elements_stencil; //CSR data storage @@ -465,35 +457,159 @@ class PolarSplineFEMPoissonLikeSolver Kokkos::View col_idx_csr_host("idx_csr", n_matrix_elements); Kokkos::View - nnz_per_row_csr("nnz_per_row_csr", matrix_size + 1); + nnz_per_row_csr("nnz_per_row_csr", m_matrix_size + 1); Kokkos::View - nnz_per_row_csr_host("nnz_per_row_csr", matrix_size + 1); + nnz_per_row_csr_host("nnz_per_row_csr", m_matrix_size + 1); + + fill_int_volume(mapping); + + m_gko_matrix = std::make_unique>(1, m_matrix_size, n_matrix_elements); + auto [values, col_idx, nnz_per_row] = m_gko_matrix->get_batch_csr(); + init_nnz_per_line(nnz_per_row); + Kokkos::deep_copy(nnz_per_row_csr_host, nnz_per_row); + + compute_singular_elements( + coeff_alpha, + coeff_beta, + mapping, + spline_evaluator, + values_csr_host, + col_idx_csr_host, + nnz_per_row_csr_host); + compute_overlapping_singular_elements( + coeff_alpha, + coeff_beta, + mapping, + spline_evaluator, + values_csr_host, + col_idx_csr_host, + nnz_per_row_csr_host); + compute_stencil_elements( + coeff_alpha, + coeff_beta, + mapping, + spline_evaluator, + values_csr_host, + col_idx_csr_host, + nnz_per_row_csr_host); + + assert(nnz_per_row_csr_host(m_matrix_size) == n_matrix_elements); + Kokkos::deep_copy(values, values_csr_host); + Kokkos::deep_copy(col_idx, col_idx_csr_host); + Kokkos::deep_copy(nnz_per_row, nnz_per_row_csr_host); + m_gko_matrix->setup_solver(); + } + + /** + * @brief Compute the volume integrals and stores the values in a member variable. + * + * @param[in] mapping + * The mapping from the logical index range to the physical index range where + * the equation is defined. + * + * @tparam Mapping A class describing a mapping from curvilinear coordinates to cartesian coordinates. + */ + template + void fill_int_volume(Mapping const& mapping) + { + auto weights_r_alloc = ddc::create_mirror_view_and_copy( + Kokkos::DefaultExecutionSpace(), + get_field(m_weights_r)); + auto weights_theta_alloc = ddc::create_mirror_view_and_copy( + Kokkos::DefaultExecutionSpace(), + get_field(m_weights_theta)); + DField weights_r = get_field(weights_r_alloc); + DField weights_theta = get_field(weights_theta_alloc); + // Find the integral volume associated with each point used in the quadrature scheme + IdxRangeQuadratureRTheta const + all_quad_points(m_idxrange_quadrature_r, m_idxrange_quadrature_theta); + auto int_volume_alloc = ddc::create_mirror_view_and_copy( + Kokkos::DefaultExecutionSpace(), + get_field(m_int_volume)); + auto int_volume = get_field(int_volume_alloc); + ddc::parallel_for_each( + Kokkos::DefaultExecutionSpace(), + all_quad_points, + KOKKOS_LAMBDA(IdxQuadratureRTheta const irp) { + IdxQuadratureR const idx_r(irp); + IdxQuadratureTheta const idx_theta(irp); + CoordRTheta coord(ddc::coordinate(irp)); + int_volume(idx_r, idx_theta) = Kokkos::abs(mapping.jacobian(coord)) + * weights_r(idx_r) * weights_theta(idx_theta); + }); + } - init_nnz_per_line(nnz_per_row_csr); - Kokkos::deep_copy(nnz_per_row_csr_host, nnz_per_row_csr); + /** + * @brief Computes the matrix element corresponding to the singular area. + * ie: the region enclosing the O-point. + * + * @param[in] coeff_alpha + * The spline representation of the @f$ \alpha @f$ function in the + * definition of the Poisson-like equation. + * @param[in] coeff_beta + * The spline representation of the @f$ \beta @f$ function in the + * definition of the Poisson-like equation. + * @param[in] mapping + * The mapping from the logical index range to the physical index range where + * the equation is defined. + * @param[in] spline_evaluator + * An evaluator for evaluating 2D splines on @f$(r,\theta)@f$. + * @param[out] values_csr_host + * A 2D Kokkos view which stores the values of non-zero elements for the whole batch. + * @param[out] col_idx_csr_host + * A 1D Kokkos view which stores the column indices for each non-zero component.(only for one matrix). + * @param[inout] nnz_per_row_csr_host + * A 1D Kokkos view of length matrix_size+1 which stores the count of the non-zeros along the lines of the matrix. + */ + template + void compute_singular_elements( + ConstSpline2D coeff_alpha, + ConstSpline2D coeff_beta, + Mapping const& mapping, + SplineRThetaEvaluatorNullBound const& spline_evaluator, + Kokkos::View const values_csr_host, + Kokkos::View const col_idx_csr_host, + Kokkos::View const nnz_per_row_csr_host) + { + IdxRangeBSPolar idxrange_singular + = PolarBSplinesRTheta::template singular_idx_range(); + IdxRangeQuadratureRTheta idxrange_quadrature_singular = m_idxrange_quadrature_singular; + + auto singular_basis_vals_and_derivs_alloc = ddc::create_mirror_view_and_copy( + Kokkos::DefaultExecutionSpace(), + get_field(m_singular_basis_vals_and_derivs)); + auto r_basis_vals_and_derivs_alloc = ddc::create_mirror_view_and_copy( + Kokkos::DefaultExecutionSpace(), + get_field(m_r_basis_vals_and_derivs)); + Field> + singular_basis_vals_and_derivs = get_field(singular_basis_vals_and_derivs_alloc); + DField int_volume_proxy = get_field(m_int_volume); Kokkos::Profiling::pushRegion("PolarPoissonFillFemMatrix"); // Calculate the matrix elements corresponding to the bsplines which cover the singular point ddc::for_each(idxrange_singular, [&](IdxBSPolar const idx_test) { ddc::for_each(idxrange_singular, [&](IdxBSPolar const idx_trial) { // Calculate the weak integral - double const element = ddc::transform_reduce( - m_idxrange_quadrature_singular, + double const element = ddc::parallel_transform_reduce( + Kokkos::DefaultExecutionSpace(), + idxrange_quadrature_singular, 0.0, ddc::reducer::sum(), - [&](IdxQuadratureRTheta const idx_quad) { - IdxQuadratureR const idx_r = ddc::select(idx_quad); - IdxQuadratureTheta const idx_theta - = ddc::select(idx_quad); + KOKKOS_LAMBDA(Idx const& idx_quad) { + Idx const idx_r(idx_quad); + Idx const idx_theta(idx_quad); return weak_integral_element( idx_r, idx_theta, - m_singular_basis_vals_and_derivs(idx_test, idx_r, idx_theta), - m_singular_basis_vals_and_derivs(idx_trial, idx_r, idx_theta), + singular_basis_vals_and_derivs(idx_test, idx_r, idx_theta), + singular_basis_vals_and_derivs(idx_trial, idx_r, idx_theta), coeff_alpha, coeff_beta, spline_evaluator, - mapping); + mapping, + int_volume_proxy); }); const int row_idx = idx_test - idxrange_singular.front(); const int col_idx = idx_trial - idxrange_singular.front(); @@ -504,8 +620,43 @@ class PolarSplineFEMPoissonLikeSolver nnz_per_row_csr_host(row_idx + 1)++; }); }); + } + /** + * @brief Computes the matrix element corresponding to singular elements overlapping + * with regular grid. + * + * @param[in] coeff_alpha + * The spline representation of the @f$ \alpha @f$ function in the + * definition of the Poisson-like equation. + * @param[in] coeff_beta + * The spline representation of the @f$ \beta @f$ function in the + * definition of the Poisson-like equation. + * @param[in] mapping + * The mapping from the logical index range to the physical index range where + * the equation is defined. + * @param[in] spline_evaluator + * An evaluator for evaluating 2D splines on @f$(r,\theta)@f$. + * @param[out] values_csr_host + * A 2D Kokkos view which stores the values of non-zero elements for the whole batch. + * @param[out] col_idx_csr_host + * A 1D Kokkos view which stores the column indices for each non-zero component.(only for one matrix) + * @param[inout] nnz_per_row_csr_host + * A 1D Kokkos view of length matrix_size+1 which stores the count of the non-zeros along the lines of the matrix. + */ + template + void compute_overlapping_singular_elements( + ConstSpline2D coeff_alpha, + ConstSpline2D coeff_beta, + Mapping const& mapping, + SplineRThetaEvaluatorNullBound const& spline_evaluator, + Kokkos::View const values_csr_host, + Kokkos::View const col_idx_csr_host, + Kokkos::View const nnz_per_row_csr_host) + { // Create index ranges associated with the 2D splines + IdxRangeBSPolar idxrange_singular + = PolarBSplinesRTheta::template singular_idx_range(); IdxRangeBSR central_radial_bspline_idx_range( m_idxrange_bsplines_r.take_first(IdxStep {BSplinesR::degree()})); @@ -513,6 +664,22 @@ class PolarSplineFEMPoissonLikeSolver central_radial_bspline_idx_range, m_idxrange_bsplines_theta); + auto singular_basis_vals_and_derivs_alloc = ddc::create_mirror_view_and_copy( + Kokkos::DefaultExecutionSpace(), + get_field(m_singular_basis_vals_and_derivs)); + auto r_basis_vals_and_derivs_alloc = ddc::create_mirror_view_and_copy( + Kokkos::DefaultExecutionSpace(), + get_field(m_r_basis_vals_and_derivs)); + auto theta_basis_vals_and_derivs_alloc = ddc::create_mirror_view_and_copy( + Kokkos::DefaultExecutionSpace(), + get_field(m_theta_basis_vals_and_derivs)); + DField int_volume_proxy = get_field(m_int_volume); + Field> + singular_basis_vals_and_derivs = get_field(singular_basis_vals_and_derivs_alloc); + Field> r_basis_vals_and_derivs + = get_field(r_basis_vals_and_derivs_alloc); + Field> + theta_basis_vals_and_derivs = get_field(theta_basis_vals_and_derivs_alloc); // Calculate the matrix elements where bspline products overlap the bsplines which cover the singular point ddc::for_each(idxrange_singular, [&](IdxBSPolar const idx_test) { ddc::for_each(idxrange_non_singular_near_centre, [&](IdxBSRTheta const idx_trial) { @@ -554,29 +721,28 @@ class PolarSplineFEMPoissonLikeSolver Idx ib_trial_theta( theta_mod(idx_trial_theta.uid() - cell_idx_theta)); // Calculate the weak integral - element += ddc::transform_reduce( + element += ddc::parallel_transform_reduce( + Kokkos::DefaultExecutionSpace(), cell_quad_points, 0.0, ddc::reducer::sum(), - [&](IdxQuadratureRTheta const idx_quad) { - IdxQuadratureR const idx_r = ddc::select(idx_quad); - IdxQuadratureTheta const idx_theta - = ddc::select(idx_quad); + KOKKOS_LAMBDA(IdxQuadratureRTheta const idx_quad) { + IdxQuadratureR const idx_r(idx_quad); + IdxQuadratureTheta const idx_theta(idx_quad); return weak_integral_element( idx_r, idx_theta, - m_singular_basis_vals_and_derivs( + singular_basis_vals_and_derivs( idx_test, idx_r, idx_theta), r_basis_vals_and_derivs(ib_trial_r, idx_r), - m_theta_basis_vals_and_derivs( - ib_trial_theta, - idx_theta), + theta_basis_vals_and_derivs(ib_trial_theta, idx_theta), coeff_alpha, coeff_beta, spline_evaluator, - mapping); + mapping, + int_volume_proxy); }); }); @@ -593,6 +759,42 @@ class PolarSplineFEMPoissonLikeSolver } }); }); + } + + /** + * @brief Computes the matrix element corresponding to the regular stencil + * ie: out to singular or overlapping areas. + * + * @param[in] coeff_alpha + * The spline representation of the @f$ \alpha @f$ function in the + * definition of the Poisson-like equation. + * @param[in] coeff_beta + * The spline representation of the @f$ \beta @f$ function in the + * definition of the Poisson-like equation. + * @param[in] mapping + * The mapping from the logical index range to the physical index range where + * the equation is defined. + * @param[in] spline_evaluator + * An evaluator for evaluating 2D splines on @f$(r,\theta)@f$. + * @param[out] values_csr_host + * A 2D Kokkos view which stores the values of non-zero elements for the whole batch. + * @param[out] col_idx_csr_host + * A 1D Kokkos view which stores the column indices for each non-zero component.(only for one matrix) + * @param[inout] nnz_per_row_csr_host + * A 1D Kokkos view of length matrix_size+1 which stores the count of the non-zeros along the lines of the matrix. + */ + template + void compute_stencil_elements( + ConstSpline2D coeff_alpha, + ConstSpline2D coeff_beta, + Mapping const& mapping, + SplineRThetaEvaluatorNullBound const& spline_evaluator, + Kokkos::View const values_csr_host, + Kokkos::View const col_idx_csr_host, + Kokkos::View const nnz_per_row_csr_host) + { + IdxRangeBSPolar idxrange_singular + = PolarBSplinesRTheta::template singular_idx_range(); // Calculate the matrix elements following a stencil ddc::for_each(m_idxrange_fem_non_singular, [&](IdxBSPolar const idx_test_polar) { @@ -682,19 +884,9 @@ class PolarSplineFEMPoissonLikeSolver } }); }); - assert(nnz_per_row_csr_host(matrix_size) == n_matrix_elements); - m_gko_matrix = std::make_unique>(1, matrix_size, n_matrix_elements); - auto [values, col_idx, nnz_per_row] = m_gko_matrix->get_batch_csr(); - Kokkos::deep_copy(values, values_csr_host); - Kokkos::deep_copy(col_idx, col_idx_csr_host); - Kokkos::deep_copy(nnz_per_row, nnz_per_row_csr_host); - m_gko_matrix->setup_solver(); Kokkos::Profiling::popRegion(); } - /** * @brief Solve the Poisson-like equation. * @@ -727,6 +919,7 @@ class PolarSplineFEMPoissonLikeSolver Kokkos::View x_init_host("x_init_host", batch_size, b_size); // Fill b + auto int_volume_host = ddc::create_mirror_view_and_copy(get_field(m_int_volume)); ddc::for_each( PolarBSplinesRTheta::template singular_idx_range(), [&](IdxBSPolar const idx) { @@ -739,17 +932,17 @@ class PolarSplineFEMPoissonLikeSolver 0.0, ddc::reducer::sum(), [&](IdxQuadratureRTheta const idx_quad) { - IdxQuadratureR const idx_r = ddc::select(idx_quad); - IdxQuadratureTheta const idx_theta - = ddc::select(idx_quad); + IdxQuadratureR const idx_r(idx_quad); + IdxQuadratureTheta const idx_theta(idx_quad); CoordRTheta coord(ddc::coordinate(idx_quad)); return rhs(coord) * m_singular_basis_vals_and_derivs(idx, idx_r, idx_theta) .value - * m_int_volume(idx_r, idx_theta); + * int_volume_host(idx_r, idx_theta); }); }); const std::size_t ncells_r = ddc::discrete_space().ncells(); + ddc::for_each(m_idxrange_fem_non_singular, [&](IdxBSPolar const idx) { const IdxBSRTheta idx_2d(PolarBSplinesRTheta::get_2d_index(idx)); const std::size_t idx_r(ddc::select(idx_2d).uid()); @@ -791,13 +984,12 @@ class PolarSplineFEMPoissonLikeSolver 0.0, ddc::reducer::sum(), [&](IdxQuadratureRTheta const idx_quad) { - IdxQuadratureR const idx_r = ddc::select(idx_quad); - IdxQuadratureTheta const idx_theta - = ddc::select(idx_quad); + IdxQuadratureR const idx_r(idx_quad); + IdxQuadratureTheta const idx_theta(idx_quad); CoordRTheta coord(ddc::coordinate(idx_quad)); - double rb = r_basis_vals_and_derivs(ib_r, idx_r).value; + double rb = m_r_basis_vals_and_derivs(ib_r, idx_r).value; double pb = m_theta_basis_vals_and_derivs(ib_theta, idx_theta).value; - return rhs(coord) * rb * pb * m_int_volume(idx_r, idx_theta); + return rhs(coord) * rb * pb * int_volume_host(idx_r, idx_theta); }); }); const std::size_t singular_index @@ -820,7 +1012,6 @@ class PolarSplineFEMPoissonLikeSolver m_idxrange_bsplines_theta); IdxRangeBSTheta idxrange_polar(ddc::discrete_space().full_domain()); - // Fill the spline ddc::for_each( PolarBSplinesRTheta::template singular_idx_range(), @@ -854,7 +1045,6 @@ class PolarSplineFEMPoissonLikeSolver Kokkos::Profiling::popRegion(); } - /** * @brief Solve the Poisson-like equation. * @@ -869,24 +1059,39 @@ class PolarSplineFEMPoissonLikeSolver * The values of the solution @f$\phi@f$ on the given coords_eval, also used as initial data for the iterative solver. */ template - void operator()(RHSFunction const& rhs, host_t phi) const + void operator()(RHSFunction const& rhs, DFieldRTheta phi) const { static_assert( std::is_invocable_r_v, "RHSFunction must have an operator() which takes a coordinate and returns a " "double"); - (*this)(rhs, m_phi_spline_coef); - host_t coords_eval_alloc(get_idx_range(phi)); - host_t coords_eval(get_field(coords_eval_alloc)); - ddc::for_each(get_idx_range(phi), [&](IdxRTheta idx) { - coords_eval(idx) = ddc::coordinate(idx); - }); - m_polar_spline_evaluator(phi, get_const_field(coords_eval), m_phi_spline_coef); + CoordFieldMemRTheta coords_eval_alloc(get_idx_range(phi)); + CoordFieldRTheta coords_eval(get_field(coords_eval_alloc)); + ddc::parallel_for_each( + Kokkos::DefaultExecutionSpace(), + get_idx_range(phi), + KOKKOS_LAMBDA(IdxRTheta idx) { coords_eval(idx) = ddc::coordinate(idx); }); + auto coords_eval_host = ddc::create_mirror_and_copy(coords_eval); + auto phi_host = ddc::create_mirror_and_copy(phi); + m_polar_spline_evaluator( + get_field(phi_host), + get_const_field(coords_eval_host), + get_const_field(m_phi_spline_coef)); + ddc::parallel_deepcopy(phi, phi_host); } -private: + /** + * @brief compute the quadrature range for a given pair of indices + * + * @param[in] cell_idx_r + * The index for radial direction + * @param[in] cell_idx_theta + * The index for poloidal direction + * @return + * The quadrature range corresponding to the @f$(r,\theta)@f$ indices. + */ static KOKKOS_FUNCTION IdxRangeQuadratureRTheta get_quadrature_points_in_cell(int cell_idx_r, int cell_idx_theta) { @@ -899,16 +1104,44 @@ class PolarSplineFEMPoissonLikeSolver return IdxRangeQuadratureRTheta(quad_points_r, quad_points_theta); } + /** + * @brief compute the weak integral value. + * + * @param[in] idx_r + * The index for radial direction. + * @param[in] idx_theta + * The index for poloidal direction + * @param[in] test_bspline_val_and_deriv + * The data structure containing the derivatives over radial and poloidal directions for test space. + * @param[in] trial_bspline_val_and_deriv + * The data structure containing the derivatives over radial and poloidal directions for trial space. + * @param[in] coeff_alpha + * The spline representation of the @f$ \alpha @f$ function in the + * definition of the Poisson-like equation. + * @param[in] coeff_beta + * The spline representation of the @f$ \beta @f$ function in the + * definition of the Poisson-like equation. + * @param[in] mapping + * The mapping from the logical index range to the physical index range where + * the equation is defined. + * @param[in] evaluator + * An evaluator for evaluating 2D splines on @f$(r,\theta)@f$. + * @param[in] int_volume + * The integral volume associated with each point used in the quadrature scheme. + * @return + * The value of the weak integral. + */ template - double weak_integral_element( + static KOKKOS_FUNCTION double weak_integral_element( IdxQuadratureR idx_r, IdxQuadratureTheta idx_theta, EvalDeriv2DType const& test_bspline_val_and_deriv, EvalDeriv2DType const& trial_bspline_val_and_deriv, - host_t coeff_alpha, - host_t coeff_beta, - SplineRThetaEvaluatorNullBound_host const& evaluator, - Mapping const& mapping) + ConstSpline2D coeff_alpha, + ConstSpline2D coeff_beta, + SplineRThetaEvaluatorNullBound const& evaluator, + Mapping const& mapping, + DField int_volume) { return templated_weak_integral_element( idx_r, @@ -920,20 +1153,23 @@ class PolarSplineFEMPoissonLikeSolver coeff_alpha, coeff_beta, evaluator, - mapping); + mapping, + int_volume); } + ///@cond template - double weak_integral_element( + static KOKKOS_FUNCTION double weak_integral_element( IdxQuadratureR idx_r, IdxQuadratureTheta idx_theta, EvalDeriv2DType const& test_bspline_val_and_deriv, EvalDeriv1DType const& trial_bspline_val_and_deriv_r, EvalDeriv1DType const& trial_bspline_val_and_deriv_theta, - host_t coeff_alpha, - host_t coeff_beta, - SplineRThetaEvaluatorNullBound_host const& evaluator, - Mapping const& mapping) + ConstSpline2D coeff_alpha, + ConstSpline2D coeff_beta, + SplineRThetaEvaluatorNullBound const& evaluator, + Mapping const& mapping, + DField int_volume) { return templated_weak_integral_element( idx_r, @@ -945,20 +1181,22 @@ class PolarSplineFEMPoissonLikeSolver coeff_alpha, coeff_beta, evaluator, - mapping); + mapping, + int_volume); } template - double weak_integral_element( + static KOKKOS_FUNCTION double weak_integral_element( IdxQuadratureR idx_r, IdxQuadratureTheta idx_theta, EvalDeriv1DType const& test_bspline_val_and_deriv_r, EvalDeriv2DType const& trial_bspline_val_and_deriv, EvalDeriv1DType const& test_bspline_val_and_deriv_theta, - host_t coeff_alpha, - host_t coeff_beta, - SplineRThetaEvaluatorNullBound_host const& evaluator, - Mapping const& mapping) + ConstSpline2D coeff_alpha, + ConstSpline2D coeff_beta, + SplineRThetaEvaluatorNullBound const& evaluator, + Mapping const& mapping, + DField int_volume) { return templated_weak_integral_element( idx_r, @@ -970,21 +1208,23 @@ class PolarSplineFEMPoissonLikeSolver coeff_alpha, coeff_beta, evaluator, - mapping); + mapping, + int_volume); } template - double weak_integral_element( + static KOKKOS_FUNCTION double weak_integral_element( IdxQuadratureR idx_r, IdxQuadratureTheta idx_theta, EvalDeriv1DType const& test_bspline_val_and_deriv_r, EvalDeriv1DType const& trial_bspline_val_and_deriv_r, EvalDeriv1DType const& test_bspline_val_and_deriv_theta, EvalDeriv1DType const& trial_bspline_val_and_deriv_theta, - host_t coeff_alpha, - host_t coeff_beta, - SplineRThetaEvaluatorNullBound_host const& evaluator, - Mapping const& mapping) + ConstSpline2D coeff_alpha, + ConstSpline2D coeff_beta, + SplineRThetaEvaluatorNullBound const& evaluator, + Mapping const& mapping, + DField int_volume) { return templated_weak_integral_element( idx_r, @@ -996,24 +1236,47 @@ class PolarSplineFEMPoissonLikeSolver coeff_alpha, coeff_beta, evaluator, - mapping); + mapping, + int_volume); } + ///@endcond - inline void get_value_and_gradient( + /** + * @brief Computes the value and gradient from r_basis and theta_basis inputs. + * + * @param[out] value The product of radial and poloidal values. + * + * @param[out] gradient derivatives over @f$ (r, \theta) @f$ directions. + * + * @param[in] r_basis A data structure containing values and derivative over radial direction. + * + * @param[in] theta_basis A data structure containing values and derivative over poloidal direction. + */ + static KOKKOS_INLINE_FUNCTION void get_value_and_gradient( double& value, std::array& gradient, EvalDeriv1DType const& r_basis, - EvalDeriv1DType const& theta_basis) const + EvalDeriv1DType const& theta_basis) { value = r_basis.value * theta_basis.value; gradient = {r_basis.derivative * theta_basis.value, r_basis.value * theta_basis.derivative}; } - inline void get_value_and_gradient( + /** + * @brief Computes the value and gradient from r_basis and theta_basis inputs. + * + * @param[out] value The product of radial and poloidal values. + * + * @param[out] gradient derivatives over @f$ (r, \theta) @f$ directions. + * + * @param[in] basis A data structure containing values and derivative over radial and poloidal directions. + * + */ + static KOKKOS_INLINE_FUNCTION void get_value_and_gradient( double& value, std::array& gradient, EvalDeriv2DType const& basis, - EvalDeriv2DType const&) const // Last argument is duplicate + EvalDeriv2DType const&) // Last argument is duplicate { value = basis.value; gradient = {basis.radial_derivative, basis.poloidal_derivative}; @@ -1021,25 +1284,51 @@ class PolarSplineFEMPoissonLikeSolver /** * @brief Computes a quadrature summand corresponding to the - * inner product - * - * The inner product of the test and trial spline is computed using a + * inner product. + * + * @param[in] idx_r + * The index for radial direction. + * @param[in] idx_theta + * The index for poloidal direction + * @param[in] test_bspline_val_and_deriv + * The data structure containing the derivatives over radial and poloidal directions for test space. + * @param[in] trial_bspline_val_and_deriv + * The data structure containing the derivatives over radial and poloidal directions for trial space. + * @param[in] test_bspline_val_and_deriv_theta + * The data structure containing the value and derivative along poloidal direction for test space. + * @param[in] trial_bspline_val_and_deriv_theta + * The data structure containing the value and derivative along poloidal direction for trial space. + * @param[in] coeff_alpha + * The spline representation of the @f$ \alpha @f$ function in the + * definition of the Poisson-like equation. + * @param[in] coeff_beta + * The spline representation of the @f$ \beta @f$ function in the + * definition of the Poisson-like equation. + * @param[in] spline_evaluator + * An evaluator for evaluating 2D splines on @f$(r,\theta)@f$. + * @param[in] mapping + * The mapping from the logical index range to the physical index range where + * the equation is defined. + * @param[in] int_volume + * The integral volume associated with each point used in the quadrature scheme. + * @return + inner product of the test and trial spline is computed using a * quadrature. This function returns one summand of the quadrature for * the quadrature point given by the indices. */ - template - double templated_weak_integral_element( + static KOKKOS_FUNCTION double templated_weak_integral_element( IdxQuadratureR idx_r, IdxQuadratureTheta idx_theta, TestValDerivType const& test_bspline_val_and_deriv, TrialValDerivType const& trial_bspline_val_and_deriv, TestValDerivType const& test_bspline_val_and_deriv_theta, TrialValDerivType const& trial_bspline_val_and_deriv_theta, - host_t coeff_alpha, - host_t coeff_beta, - SplineRThetaEvaluatorNullBound_host const& spline_evaluator, - Mapping const& mapping) + ConstSpline2D coeff_alpha, + ConstSpline2D coeff_beta, + SplineRThetaEvaluatorNullBound const& spline_evaluator, + Mapping const& mapping, + DField int_volume) { static_assert( std::is_same_v< @@ -1074,7 +1363,7 @@ class PolarSplineFEMPoissonLikeSolver MetricTensor metric_tensor(mapping); // Assemble the weak integral element - return m_int_volume(idx_r, idx_theta) + return int_volume(idx_r, idx_theta) * (alpha * dot_product( basis_gradient_test_space, @@ -1085,14 +1374,32 @@ class PolarSplineFEMPoissonLikeSolver /** * @brief Computes the matrix element corresponding to two tensor product splines * with index idx_test and idx_trial + * + * @param[in] idx_test + * The index for polar B-spline in the test space. + * @param[in] idx_trial + * The index for polar B-spline in the trial space. + * @param[in] coeff_alpha + * The spline representation of the @f$ \alpha @f$ function in the + * definition of the Poisson-like equation. + * @param[in] coeff_beta + * The spline representation of the @f$ \beta @f$ function in the + * definition of the Poisson-like equation. + * @param[in] evaluator + * An evaluator for evaluating 2D splines on @f$ (r, \theta) @f$. + * @param[in] mapping + * The mapping from the logical index range to the physical index range where + * the equation is defined. + * @return + * The value of the matrix element. */ template double get_matrix_stencil_element( IdxBSRTheta idx_test, IdxBSRTheta idx_trial, - host_t coeff_alpha, - host_t coeff_beta, - SplineRThetaEvaluatorNullBound_host const& evaluator, + ConstSpline2D coeff_alpha, + ConstSpline2D coeff_beta, + SplineRThetaEvaluatorNullBound const& evaluator, Mapping const& mapping) { // 0 <= idx_test_r < 8 @@ -1151,6 +1458,19 @@ class PolarSplineFEMPoissonLikeSolver const IdxRange theta_cells(first_overlap_element_theta, n_overlap_theta); const IdxRange non_zero_cells(r_cells, theta_cells); + auto r_basis_vals_and_derivs_alloc = ddc::create_mirror_view_and_copy( + Kokkos::DefaultExecutionSpace(), + get_field(m_r_basis_vals_and_derivs)); + auto theta_basis_vals_and_derivs_alloc = ddc::create_mirror_view_and_copy( + Kokkos::DefaultExecutionSpace(), + get_field(m_theta_basis_vals_and_derivs)); + + Field> r_basis_vals_and_derivs + = get_field(r_basis_vals_and_derivs_alloc); + Field> + theta_basis_vals_and_derivs = get_field(theta_basis_vals_and_derivs_alloc); + DField int_volume_proxy = get_field(m_int_volume); + assert(n_overlap_r * n_overlap_theta > 0); return ddc::transform_reduce( non_zero_cells, @@ -1178,29 +1498,37 @@ class PolarSplineFEMPoissonLikeSolver assert(ib_trial_theta.uid() < BSplinesTheta::degree() + 1); // Calculate the weak integral - return ddc::transform_reduce( + return ddc::parallel_transform_reduce( + Kokkos::DefaultExecutionSpace(), cell_quad_points, 0.0, ddc::reducer::sum(), - [&](IdxQuadratureRTheta const idx_quad) { - IdxQuadratureR const idx_r = ddc::select(idx_quad); - IdxQuadratureTheta const idx_theta - = ddc::select(idx_quad); + KOKKOS_LAMBDA(IdxQuadratureRTheta const idx_quad) { + IdxQuadratureR const idx_r(idx_quad); + IdxQuadratureTheta const idx_theta(idx_quad); return weak_integral_element( idx_r, idx_theta, r_basis_vals_and_derivs(ib_test_r, idx_r), r_basis_vals_and_derivs(ib_trial_r, idx_r), - m_theta_basis_vals_and_derivs(ib_test_theta, idx_theta), - m_theta_basis_vals_and_derivs(ib_trial_theta, idx_theta), + theta_basis_vals_and_derivs(ib_test_theta, idx_theta), + theta_basis_vals_and_derivs(ib_trial_theta, idx_theta), coeff_alpha, coeff_beta, evaluator, - mapping); + mapping, + int_volume_proxy); }); }); } + /** + * @brief Calculates the modulo idx_theta in relation to cells number along @f$ \theta @f$ direction . + * + * @param[in] idx_theta @f$ \theta @f$ index. + * + * @return The corresponding indice modulo @f$ \theta @f$ direction cells number + */ static KOKKOS_FUNCTION int theta_mod(int idx_theta) { int ncells_theta = ddc::discrete_space().ncells(); @@ -1211,20 +1539,19 @@ class PolarSplineFEMPoissonLikeSolver return idx_theta; } -public: /** - * @brief Fills the nnz data structure by computing the number of non-zero per line. - * This number is linked to the weak formulation and depends on @f$(r,\theta)@f$ splines. - * After this function the array will contain: - * nnz[0] = 0. - * nnz[1] = 0. - * nnz[2] = number of non-zero elements in line 0. - * nnz[3] = number of non-zero elements in lines 0-1. - * ... - * nnz[matrix_size] = number of non-zero elements in lines 0-(matrix_size-1). - * - * @param[out] nnz A 1D Kokkos view of length matrix_size+1 which stores the count of the non-zeros along the lines of the matrix. - */ + * @brief Fills the nnz data structure by computing the number of non-zero per line. + * This number is linked to the weak formulation and depends on @f$(r,\theta)@f$ splines. + * After this function the array will contain: + * nnz[0] = 0. + * nnz[1] = 0. + * nnz[2] = number of non-zero elements in line 0. + * nnz[3] = number of non-zero elements in lines 0-1. + * ... + * nnz[matrix_size] = number of non-zero elements in lines 0-(matrix_size-1). + * + * @param[out] nnz A 1D Kokkos view of length matrix_size+1 which stores the count of the non-zeros along the lines of the matrix. + */ void init_nnz_per_line(Kokkos::View nnz) const { Kokkos::Profiling::pushRegion("PolarPoissonInitNnz"); diff --git a/tests/geometryRTheta/polar_poisson/CMakeLists.txt b/tests/geometryRTheta/polar_poisson/CMakeLists.txt index 08965f25..ecd5a474 100644 --- a/tests/geometryRTheta/polar_poisson/CMakeLists.txt +++ b/tests/geometryRTheta/polar_poisson/CMakeLists.txt @@ -26,7 +26,7 @@ foreach(MAPPING_TYPE "CIRCULAR_MAPPING" "CZARNY_MAPPING") add_test(NAME TestPoissonConvergence_${MAPPING_TYPE}_${SOLUTION} COMMAND "$" "${CMAKE_CURRENT_SOURCE_DIR}/test_poisson.py" "$") - set_property(TEST TestPoissonConvergence_${MAPPING_TYPE}_${SOLUTION} PROPERTY TIMEOUT 200) + set_property(TEST TestPoissonConvergence_${MAPPING_TYPE}_${SOLUTION} PROPERTY TIMEOUT 300) set_property(TEST TestPoissonConvergence_${MAPPING_TYPE}_${SOLUTION} PROPERTY COST 100) endforeach() endforeach() diff --git a/tests/geometryRTheta/polar_poisson/polarpoissonfemsolver.cpp b/tests/geometryRTheta/polar_poisson/polarpoissonfemsolver.cpp index e5d16321..919908f8 100644 --- a/tests/geometryRTheta/polar_poisson/polarpoissonfemsolver.cpp +++ b/tests/geometryRTheta/polar_poisson/polarpoissonfemsolver.cpp @@ -22,18 +22,22 @@ #include "polarpoissonlikesolver.hpp" #include "test_cases.hpp" + using PoissonSolver = PolarSplineFEMPoissonLikeSolver< GridR, GridTheta, PolarBSplinesRTheta, - SplineRThetaEvaluatorNullBound_host>; + SplineRThetaEvaluatorNullBound>; #if defined(CIRCULAR_MAPPING) using Mapping = CircularToCartesian; #elif defined(CZARNY_MAPPING) using Mapping = CzarnyToCartesian; #endif -using DiscreteMappingBuilder = DiscreteToCartesianBuilder< +using DiscreteMappingBuilder + = DiscreteToCartesianBuilder; + +using DiscreteMappingBuilder_host = DiscreteToCartesianBuilder< X, Y, SplineRThetaBuilder_host, @@ -97,56 +101,74 @@ int main(int argc, char** argv) IdxRangeTheta interpolation_idx_range_P(SplineInterpPointsTheta::get_domain()); IdxRangeRTheta grid(interpolation_idx_range_R, interpolation_idx_range_P); - SplineRThetaBuilder_host const builder(grid); + SplineRThetaBuilder const builder(grid); + SplineRThetaBuilder_host const builder_host(grid); + #if defined(CIRCULAR_MAPPING) const Mapping mapping; #elif defined(CZARNY_MAPPING) const Mapping mapping(0.3, 1.4); #endif + ddc::NullExtrapolationRule bv_r_min; ddc::NullExtrapolationRule bv_r_max; ddc::PeriodicExtrapolationRule bv_p_min; ddc::PeriodicExtrapolationRule bv_p_max; - SplineRThetaEvaluatorNullBound_host evaluator(bv_r_min, bv_r_max, bv_p_min, bv_p_max); - DiscreteMappingBuilder const discrete_mapping_builder( + SplineRThetaEvaluatorNullBound evaluator(bv_r_min, bv_r_max, bv_p_min, bv_p_max); + SplineRThetaEvaluatorNullBound_host evaluator_host(bv_r_min, bv_r_max, bv_p_min, bv_p_max); + + DiscreteMappingBuilder_host const discrete_mapping_builder_host( Kokkos::DefaultHostExecutionSpace(), mapping, - builder, - evaluator); - DiscreteToCartesian const discrete_mapping = discrete_mapping_builder(); + builder_host, + evaluator_host); - ddc::init_discrete_space(discrete_mapping); - IdxRangeBSRTheta idx_range_bsplinesRTheta = get_spline_idx_range(builder); + DiscreteMappingBuilder const + discrete_mapping_builder(Kokkos::DefaultExecutionSpace(), mapping, builder, evaluator); + DiscreteToCartesian const discrete_mapping = discrete_mapping_builder(); + DiscreteToCartesian const discrete_mapping_host = discrete_mapping_builder_host(); - host_t coeff_alpha(grid); // values of the coefficent alpha - host_t coeff_beta(grid); - host_t x(grid); - host_t y(grid); + ddc::init_discrete_space(discrete_mapping_host); - ddc::for_each(grid, [&](IdxRTheta const irp) { - coeff_alpha(irp) - = std::exp(-std::tanh((ddc::coordinate(ddc::select(irp)) - 0.7) / 0.05)); - coeff_beta(irp) = 1.0 / coeff_alpha(irp); - Coord - coord(ddc::coordinate(ddc::select(irp)), - ddc::coordinate(ddc::select(irp))); - Coord cartesian_coord = mapping(coord); - x(irp) = ddc::get(cartesian_coord); - y(irp) = ddc::get(cartesian_coord); - }); + IdxRangeBSRTheta idx_range_bsplinesRTheta = get_spline_idx_range(builder); - host_t coeff_alpha_spline(idx_range_bsplinesRTheta); - host_t coeff_beta_spline(idx_range_bsplinesRTheta); + DFieldMemRTheta coeff_alpha_alloc(grid); // values of the coefficent alpha + DFieldMemRTheta coeff_beta_alloc(grid); + DFieldMemRTheta x_alloc(grid); + DFieldMemRTheta y_alloc(grid); + + DFieldRTheta coeff_alpha = get_field(coeff_alpha_alloc); // values of the coefficent alpha + DFieldRTheta coeff_beta = get_field(coeff_beta_alloc); + DFieldRTheta x = get_field(x_alloc); + DFieldRTheta y = get_field(y_alloc); + + ddc::parallel_for_each( + Kokkos::DefaultExecutionSpace(), + grid, + KOKKOS_LAMBDA(IdxRTheta const irp) { + coeff_alpha(irp) = Kokkos::exp( + -Kokkos::tanh((ddc::coordinate(ddc::select(irp)) - 0.7) / 0.05)); + coeff_beta(irp) = 1.0 / coeff_alpha(irp); + Coord + coord(ddc::coordinate(ddc::select(irp)), + ddc::coordinate(ddc::select(irp))); + Coord cartesian_coord = mapping(coord); + x(irp) = ddc::get(cartesian_coord); + y(irp) = ddc::get(cartesian_coord); + }); + + Spline2DMem coeff_alpha_spline(idx_range_bsplinesRTheta); + Spline2DMem coeff_beta_spline(idx_range_bsplinesRTheta); builder(get_field(coeff_alpha_spline), get_const_field(coeff_alpha)); // coeff_alpha_spline are the coefficients // of the spline representation of the values given by coeff_alpha. builder(get_field(coeff_beta_spline), get_const_field(coeff_beta)); - host_t x_spline_representation(idx_range_bsplinesRTheta); - host_t y_spline_representation(idx_range_bsplinesRTheta); + Spline2DMem x_spline_representation(idx_range_bsplinesRTheta); + Spline2DMem y_spline_representation(idx_range_bsplinesRTheta); builder(get_field(x_spline_representation), get_const_field(x)); builder(get_field(y_spline_representation), get_const_field(y)); @@ -173,7 +195,9 @@ int main(int argc, char** argv) LHSFunction lhs(mapping); RHSFunction rhs(mapping); host_t> coords(grid); - host_t result(grid); + DFieldMemRTheta result_alloc(grid); + DField result = get_field(result_alloc); + ddc::for_each(grid, [&](IdxRTheta const irp) { coords(irp) = CoordRTheta( ddc::coordinate(ddc::select(irp)), @@ -181,30 +205,30 @@ int main(int argc, char** argv) }); if (discrete_rhs) { // Build the spline approximation of the rhs + Spline2DMem rhs_spline(idx_range_bsplinesRTheta); + host_t rhs_vals_host(grid); - host_t rhs_spline(idx_range_bsplinesRTheta); - host_t rhs_vals(grid); - ddc::for_each(grid, [&](IdxRTheta const irp) { rhs_vals(irp) = rhs(coords(irp)); }); - builder(get_field(rhs_spline), get_const_field(rhs_vals)); - - + ddc::for_each(grid, [&](IdxRTheta const irp) { rhs_vals_host(irp) = rhs(coords(irp)); }); + auto rhs_vals = ddc::create_mirror_view_and_copy( + Kokkos::DefaultExecutionSpace(), + get_field(rhs_vals_host)); + builder(get_field(rhs_spline), get_const_field(rhs_vals)); + ConstSpline2D rhs_spline_field = get_const_field(rhs_spline); start_time = std::chrono::system_clock::now(); - ddc::NullExtrapolationRule r_extrapolation_rule; - ddc::PeriodicExtrapolationRule p_extrapolation_rule; - SplineRThetaEvaluatorNullBound_host - eval(r_extrapolation_rule, - r_extrapolation_rule, - p_extrapolation_rule, - p_extrapolation_rule); - solver([&](CoordRTheta const& coord) { return eval(coord, get_const_field(rhs_spline)); }, - get_field(result)); + solver( + KOKKOS_LAMBDA(CoordRTheta const& coord) { + return evaluator(coord, rhs_spline_field); + }, + get_field(result)); end_time = std::chrono::system_clock::now(); } else { start_time = std::chrono::system_clock::now(); solver(rhs, get_field(result)); end_time = std::chrono::system_clock::now(); } + auto result_alloc_host = ddc::create_mirror_view_and_copy(result); + host_t result_host = get_field(result_alloc_host); std::cout << "Solver time : " << std::chrono::duration_cast(end_time - start_time) .count() @@ -212,7 +236,7 @@ int main(int argc, char** argv) double max_err = 0.0; ddc::for_each(grid, [&](IdxRTheta const irp) { - const double err = result(irp) - lhs(coords(irp)); + const double err = result_host(irp) - lhs(coords(irp)); if (err > 0) { max_err = max_err > err ? max_err : err; } else { diff --git a/vendor/sll/include/sll/math_tools.hpp b/vendor/sll/include/sll/math_tools.hpp index bf6e1e69..c6ce6fb4 100644 --- a/vendor/sll/include/sll/math_tools.hpp +++ b/vendor/sll/include/sll/math_tools.hpp @@ -89,7 +89,7 @@ inline std::size_t factorial(std::size_t f) } template -inline T dot_product(std::array const& a, std::array const& b) +KOKKOS_INLINE_FUNCTION T dot_product(std::array const& a, std::array const& b) { T result = 0; for (std::size_t i(0); i < D; ++i) { diff --git a/vendor/sll/include/sll/polar_bsplines.hpp b/vendor/sll/include/sll/polar_bsplines.hpp index f69de5a8..a62c22d6 100644 --- a/vendor/sll/include/sll/polar_bsplines.hpp +++ b/vendor/sll/include/sll/polar_bsplines.hpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include @@ -252,6 +253,7 @@ class PolarBSplines template Impl(const DiscreteMapping& curvilinear_to_cartesian) { + static_assert(is_accessible_v); static_assert(std::is_same_v); using DimX = typename DiscreteMapping::cartesian_tag_x; using DimY = typename DiscreteMapping::cartesian_tag_y;