diff --git a/docs/CODING_STANDARD.md b/docs/CODING_STANDARD.md index 79309c460..7907d9dbc 100644 --- a/docs/CODING_STANDARD.md +++ b/docs/CODING_STANDARD.md @@ -33,6 +33,7 @@ * we don't rely on case to distinguish between variables * there are two types of DDC objects representing a multidimensional array : `ddc::Chunk` (which possesses the data) and `ddc::ChunkSpan` (which does not own the data but can be captured by `KOKKOS_LAMBDA`). We suffix `Chunk` with `_alloc` if both variables are needed locally. * if a variable is mirrored between host (CPU) and device (GPU) memories, the variable representing data on host is `_host` suffixed +* capturing classes members through `KOKKOS_LAMBDA` or `KOKKOS_CLASS_LAMBDA` may be complicated, we often need to copy-by-reference the member to a local variable, which must be `_proxy` suffixed ## Style * We use the style specified by the `.clang-format` file using clang-format 10 diff --git a/src/geometryXVx/geometry/geometry.hpp b/src/geometryXVx/geometry/geometry.hpp index 177bad152..b08fb2de4 100644 --- a/src/geometryXVx/geometry/geometry.hpp +++ b/src/geometryXVx/geometry/geometry.hpp @@ -144,8 +144,8 @@ using SplineVxEvaluator = ddc::SplineEvaluator< IDimX, IDimVx>; using SplineXBuilder_1d = ddc::SplineBuilder< - Kokkos::DefaultHostExecutionSpace, - Kokkos::DefaultHostExecutionSpace::memory_space, + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, BSplinesX, IDimX, SplineXBoundary, @@ -153,8 +153,8 @@ using SplineXBuilder_1d = ddc::SplineBuilder< ddc::SplineSolver::GINKGO, IDimX>; using SplineXEvaluator_1d = ddc::SplineEvaluator< - Kokkos::DefaultHostExecutionSpace, - Kokkos::DefaultHostExecutionSpace::memory_space, + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, BSplinesX, IDimX, #ifdef PERIODIC_RDIMX @@ -297,6 +297,9 @@ using SpanSpXVx = device_t>; template using SpanSpVx = device_t>; +template +using BSSpanX = device_t>; + template @@ -315,6 +318,7 @@ using DSpanSpXVx = SpanSpXVx; using DSpanSpVx = SpanSpVx; +using DBSSpanX = BSSpanX; template diff --git a/src/geometryXVx/poisson/CMakeLists.txt b/src/geometryXVx/poisson/CMakeLists.txt index 6ac8b4585..c7e6c3135 100644 --- a/src/geometryXVx/poisson/CMakeLists.txt +++ b/src/geometryXVx/poisson/CMakeLists.txt @@ -4,7 +4,6 @@ foreach(GEOMETRY_VARIANT IN LISTS GEOMETRY_XVx_VARIANTS_LIST) add_library("poisson_${GEOMETRY_VARIANT}" STATIC chargedensitycalculator.cpp - electricfield.cpp nullpoissonsolver.cpp ) diff --git a/src/geometryXVx/poisson/electricfield.cpp b/src/geometryXVx/poisson/electricfield.cpp deleted file mode 100644 index b6d77c414..000000000 --- a/src/geometryXVx/poisson/electricfield.cpp +++ /dev/null @@ -1,39 +0,0 @@ -// SPDX-License-Identifier: MIT - -#include - -#include "electricfield.hpp" - -ElectricField::ElectricField( - SplineXBuilder_1d const& spline_x_builder, - SplineXEvaluator_1d const& spline_x_evaluator) - : m_spline_x_builder(spline_x_builder) - , m_spline_x_evaluator(spline_x_evaluator) -{ -} - -//=========================================================================== -// Compute efield = -dPhi/dx where Phi is the electrostatic potential -// input : Phi values -//=========================================================================== -void ElectricField::operator()( - host_t const electric_field, - host_t const electrostatic_potential) const -{ - IDomainX const& x_dom = electric_field.domain(); - ddc::for_each(x_dom, [&](IndexX const ix) { - electric_field(ix) - = -m_spline_x_evaluator.deriv(ddc::coordinate(ix), electrostatic_potential); - }); -} - -void ElectricField::operator()( - host_t const electric_field, - host_t const electrostatic_potential) const -{ - Kokkos::Profiling::pushRegion("ElectricField"); - host_t> elecpot_spline_coef(m_spline_x_builder.bsplines_domain()); - m_spline_x_builder(elecpot_spline_coef.span_view(), electrostatic_potential); - (*this)(electric_field, elecpot_spline_coef); - Kokkos::Profiling::popRegion(); -} diff --git a/src/geometryXVx/poisson/electricfield.hpp b/src/geometryXVx/poisson/electricfield.hpp deleted file mode 100644 index 739a768e9..000000000 --- a/src/geometryXVx/poisson/electricfield.hpp +++ /dev/null @@ -1,50 +0,0 @@ -// SPDX-License-Identifier: MIT - -#include -#include - - -//=========================================================================== -// Compute efield = -dPhi/dx where Phi is the electrostatic potential -// input : Phi values -//=========================================================================== -/** - * @brief An operator which computes the electric field using splines derivation. - * - * An operator which computes the electric field: - * @f$ E = - \frac{d\phi}{dx} @f$ where Phi is the electrostatic potential. - * This operator uses spline derivation. - */ -class ElectricField -{ - SplineXBuilder_1d const& m_spline_x_builder; - - SplineXEvaluator_1d m_spline_x_evaluator; - -public: - /** - * Construct the ElectricField operator. - * - * @param spline_x_builder A spline builder which calculates the coefficients of a spline representation. - * @param spline_x_evaluator A spline evaluator which provides the value of a spline representation from its coefficients. - */ - ElectricField( - SplineXBuilder_1d const& spline_x_builder, - SplineXEvaluator_1d const& spline_x_evaluator); - - /** - * The operator which solves the equation using the method described by the class. - * - * @param[out] electric_field The electric_field, the result of the operation. - * @param[out] electrostatic_potential The electrostatic potential, the input of the operator. - */ - void operator()(host_t electric_field, host_t electrostatic_potential) const; - - /** - * The operator which solves the equation using the method described by the class. - * - * @param[out] electric_field The electric_field, the result of the operation. - * @param[out] electrostatic_potential The electrostatic potential, the input of the operator. - */ - void operator()(host_t electric_field, host_t electrostatic_potential) const; -}; diff --git a/src/geometryXVx/poisson/femnonperiodicpoissonsolver.cpp b/src/geometryXVx/poisson/femnonperiodicpoissonsolver.cpp index 7c815f1fd..16f2c327f 100644 --- a/src/geometryXVx/poisson/femnonperiodicpoissonsolver.cpp +++ b/src/geometryXVx/poisson/femnonperiodicpoissonsolver.cpp @@ -30,7 +30,7 @@ NUBSplineXEvaluator_1d jit_build_nubsplinesx(SplineXEvaluator_1d const& spline_x // Boundary values are never evaluated return NUBSplineXEvaluator_1d( ddc::replace_dim_of( - spline_x_evaluator.spline_domain(), + spline_x_evaluator.bsplines_domain(), ddc::discrete_space().full_domain()), ddc::NullExtrapolationRule(), ddc::NullExtrapolationRule()); @@ -51,7 +51,7 @@ FemNonPeriodicPoissonSolver::FemNonPeriodicPoissonSolver( , m_compute_rho(compute_rho) , m_nbasis(ddc::discrete_space().nbasis()) , m_ncells(ddc::discrete_space().ncells()) - , m_quad_coef(ddc::DiscreteDomain( + , m_quad_coef_alloc(ddc::DiscreteDomain( ddc::DiscreteElement(0), ddc::DiscreteVector(s_npts_gauss * m_ncells))) { @@ -67,10 +67,12 @@ FemNonPeriodicPoissonSolver::FemNonPeriodicPoissonSolver( // Calculate the integration coefficients GaussLegendre const gl(s_npts_gauss); - std::vector> eval_pts_data(m_quad_coef.domain().size()); + std::vector> eval_pts_data(m_quad_coef_alloc.domain().size()); ddc::ChunkSpan, ddc::DiscreteDomain> const - eval_pts(eval_pts_data.data(), m_quad_coef.domain()); - gl.compute_points_and_weights_on_mesh(eval_pts, m_quad_coef.span_view(), knots.span_cview()); + eval_pts(eval_pts_data.data(), m_quad_coef_alloc.domain()); + auto quad_coef_host = ddc::create_mirror_and_copy(m_quad_coef_alloc.span_view()); + gl.compute_points_and_weights_on_mesh(eval_pts, quad_coef_host.span_view(), knots.span_cview()); + ddc::deepcopy(m_quad_coef_alloc, quad_coef_host); ddc::init_discrete_space(eval_pts_data); @@ -94,9 +96,11 @@ void FemNonPeriodicPoissonSolver::build_matrix() m_fem_matrix = Matrix:: make_new_banded(matrix_size, n_lower_diags, n_lower_diags, positive_definite_symmetric); + auto quad_coef_host = ddc::create_mirror_and_copy(m_quad_coef_alloc.span_cview()); + // Fill the banded part of the matrix std::array derivs; - ddc::for_each(m_quad_coef.domain(), [&](ddc::DiscreteElement const ix) { + ddc::for_each(m_quad_coef_alloc.domain(), [&](ddc::DiscreteElement const ix) { ddc::Coordinate const coord = coord_from_quad_point(ddc::coordinate(ix)); ddc::DiscreteElement const jmin = ddc::discrete_space().eval_deriv(derivs, coord); @@ -108,7 +112,7 @@ void FemNonPeriodicPoissonSolver::build_matrix() if (j_idx != -1 && j_idx != matrix_size && k_idx != -1 && k_idx != matrix_size) { double a_jk = m_fem_matrix->get_element(j_idx, k_idx); // Update element - a_jk += derivs[j] * derivs[k] * m_quad_coef(ix); + a_jk += derivs[j] * derivs[k] * quad_coef_host(ix); m_fem_matrix->set_element(j_idx, k_idx, a_jk); } @@ -125,36 +129,46 @@ void FemNonPeriodicPoissonSolver::build_matrix() // Solve the Poisson equation //--------------------------------------------------------------------------- void FemNonPeriodicPoissonSolver::solve_matrix_system( - ddc::ChunkSpan const phi_spline_coef, - ddc::ChunkSpan const rho_spline_coef) const + DNUBSSpanX const phi_spline_coef, + DBSViewX const rho_spline_coef) const { - std::array values; - - for (int i(0); i < m_nbasis; ++i) { - phi_spline_coef(ddc::DiscreteElement(i)) = 0.0; - } + ddc::fill(phi_spline_coef, 0.0); int const rhs_size = m_nbasis - 2; - DSpan1D const phi_rhs(phi_spline_coef.data_handle() + 1, rhs_size); + int const nbasis_proxy = m_nbasis; + SplineXEvaluator_1d spline_x_evaluator_proxy = m_spline_x_evaluator; + ddc::ChunkSpan quad_coef = m_quad_coef_alloc.span_view(); // Fill phi_rhs(i) with \int rho(x) b_i(x) dx // Rk: phi_rhs no longer contains spline coefficients, but is the // RHS of the matrix equation - ddc::for_each(m_quad_coef.domain(), [&](ddc::DiscreteElement const ix) { - ddc::Coordinate const coord = coord_from_quad_point(ddc::coordinate(ix)); - ddc::DiscreteElement const jmin - = ddc::discrete_space().eval_basis(values, coord); - double const rho_val = m_spline_x_evaluator(coord, rho_spline_coef.span_cview()); - for (int j = 0; j < s_degree + 1; ++j) { - int const j_idx = (jmin.uid() + j) % m_nbasis - 1; - if (j_idx != -1 && j_idx != rhs_size) { - phi_rhs(j_idx) += rho_val * values[j] * m_quad_coef(ix); - } - } - }); + ddc::for_each( + ddc::policies::parallel_device, + m_quad_coef_alloc.domain(), + KOKKOS_LAMBDA(ddc::DiscreteElement const ix) { + ddc::Coordinate const coord = coord_from_quad_point(ddc::coordinate(ix)); + std::array values; + ddc::DiscreteElement const jmin + = ddc::discrete_space().eval_basis(values, coord); + double const rho_val + = spline_x_evaluator_proxy(coord, rho_spline_coef.span_cview()); + for (int j = 0; j < s_degree + 1; ++j) { + int const j_idx = (jmin.uid() + j) % nbasis_proxy; + if (j_idx != 0 && j_idx != rhs_size + 1) { + Kokkos::atomic_fetch_add( + &phi_spline_coef(ddc::DiscreteElement(j_idx)), + rho_val * values[j] * quad_coef(ix)); + } + } + }); + + auto phi_spline_coef_host = ddc::create_mirror_and_copy(phi_spline_coef); + DSpan1D const phi_rhs_host(phi_spline_coef_host.data_handle() + 1, rhs_size); // Solve the matrix equation to find the spline coefficients of phi - m_fem_matrix->solve_inplace(phi_rhs); + m_fem_matrix->solve_inplace(phi_rhs_host); + + ddc::deepcopy(phi_spline_coef, phi_spline_coef_host); } @@ -180,32 +194,32 @@ void FemNonPeriodicPoissonSolver::operator()( DViewSpXVx allfdistribu) const { Kokkos::Profiling::pushRegion("PoissonSolver"); - auto electrostatic_potential_host = ddc::create_mirror_and_copy(electrostatic_potential); - auto electric_field_host = ddc::create_mirror_and_copy(electric_field); assert(electrostatic_potential.domain() == ddc::get_domain(allfdistribu)); IDomainX const dom_x = electrostatic_potential.domain(); // Compute the RHS of the Poisson equation - host_t rho_host(dom_x); DFieldX rho(dom_x); m_compute_rho(rho, allfdistribu); - ddc::deepcopy(rho_host, rho); // - ddc::Chunk rho_spline_coef(m_spline_x_builder.spline_domain()); - m_spline_x_builder(rho_spline_coef.span_view(), rho_host.span_cview()); - ddc::Chunk phi_spline_coef( - ddc::discrete_space().full_domain()); - solve_matrix_system(phi_spline_coef, rho_spline_coef); + ddc::Chunk rho_spline_coef(m_spline_x_builder.spline_domain(), ddc::DeviceAllocator()); + m_spline_x_builder(rho_spline_coef.span_view(), rho.span_cview()); + ddc::Chunk phi_spline_coef( + ddc::discrete_space().full_domain(), + ddc::DeviceAllocator()); + solve_matrix_system(phi_spline_coef, rho_spline_coef.span_cview()); // - ddc::for_each(dom_x, [&](IndexX const ix) { - electrostatic_potential_host(ix) - = m_spline_x_nu_evaluator(ddc::coordinate(ix), phi_spline_coef.span_cview()); - electric_field_host(ix) - = -m_spline_x_nu_evaluator.deriv(ddc::coordinate(ix), phi_spline_coef.span_cview()); - }); - ddc::deepcopy(electrostatic_potential, electrostatic_potential_host); - ddc::deepcopy(electric_field, electric_field_host); + NUBSplineXEvaluator_1d spline_x_nu_evaluator_proxy = m_spline_x_nu_evaluator; + ddc::ChunkSpan phi_spline_coef_view = phi_spline_coef.span_cview(); + ddc::for_each( + ddc::policies::parallel_device, + dom_x, + KOKKOS_LAMBDA(IndexX const ix) { + electrostatic_potential(ix) + = spline_x_nu_evaluator_proxy(ddc::coordinate(ix), phi_spline_coef_view); + electric_field(ix) = -spline_x_nu_evaluator_proxy + .deriv(ddc::coordinate(ix), phi_spline_coef_view); + }); Kokkos::Profiling::popRegion(); } diff --git a/src/geometryXVx/poisson/femnonperiodicpoissonsolver.hpp b/src/geometryXVx/poisson/femnonperiodicpoissonsolver.hpp index 9b52f3546..26bc249df 100644 --- a/src/geometryXVx/poisson/femnonperiodicpoissonsolver.hpp +++ b/src/geometryXVx/poisson/femnonperiodicpoissonsolver.hpp @@ -15,8 +15,8 @@ using NUBSDomainX = ddc::DiscreteDomain; using UBSplinesX = ddc::UniformBSplines; using NUBSplineXEvaluator_1d = ddc::SplineEvaluator< - Kokkos::DefaultHostExecutionSpace, - Kokkos::DefaultHostExecutionSpace::memory_space, + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, NUBSplinesX, IDimX, ddc::NullExtrapolationRule, @@ -46,6 +46,8 @@ class FemNonPeriodicPoissonSolver : public IPoissonSolver private: using QMeshX = ddc::NonUniformPointSampling; + using DQFieldX = device_t>>; + using DNUBSSpanX = device_t>; private: // Spline degree in x direction @@ -69,17 +71,19 @@ class FemNonPeriodicPoissonSolver : public IPoissonSolver // Number of cells in x direction int m_ncells; - ddc::Chunk> m_quad_coef; + DQFieldX m_quad_coef_alloc; std::unique_ptr m_fem_matrix; private: - static ddc::Coordinate quad_point_from_coord(ddc::Coordinate const& coord) + static KOKKOS_FUNCTION ddc::Coordinate quad_point_from_coord( + ddc::Coordinate const& coord) { return ddc::Coordinate(ddc::get(coord)); } - static ddc::Coordinate coord_from_quad_point(ddc::Coordinate const& coord) + static KOKKOS_FUNCTION ddc::Coordinate coord_from_quad_point( + ddc::Coordinate const& coord) { return ddc::Coordinate(ddc::get(coord)); } @@ -110,7 +114,12 @@ class FemNonPeriodicPoissonSolver : public IPoissonSolver private: void build_matrix(); - void solve_matrix_system( - ddc::ChunkSpan phi_spline_coef, - ddc::ChunkSpan rho_spline_coef) const; +public: + /** + * [SHOULD BE PRIVATE (Kokkos limitation)] + * + * @param[out] phi_spline_coef + * @param[in] rho_spline_coef + */ + void solve_matrix_system(DNUBSSpanX phi_spline_coef, DBSViewX rho_spline_coef) const; }; diff --git a/src/geometryXVx/poisson/femperiodicpoissonsolver.cpp b/src/geometryXVx/poisson/femperiodicpoissonsolver.cpp index 81aa3d922..34b3503ba 100644 --- a/src/geometryXVx/poisson/femperiodicpoissonsolver.cpp +++ b/src/geometryXVx/poisson/femperiodicpoissonsolver.cpp @@ -18,7 +18,7 @@ FemPeriodicPoissonSolver::FemPeriodicPoissonSolver( , m_compute_rho(compute_rho) , m_nbasis(ddc::discrete_space().nbasis()) , m_ncells(ddc::discrete_space().ncells()) - , m_quad_coef(ddc::DiscreteDomain( + , m_quad_coef_alloc(ddc::DiscreteDomain( ddc::DiscreteElement(0), ddc::DiscreteVector(s_npts_gauss * m_ncells))) { @@ -35,10 +35,12 @@ FemPeriodicPoissonSolver::FemPeriodicPoissonSolver( // Calculate the integration coefficients GaussLegendre const gl(s_npts_gauss); - std::vector> eval_pts_data(m_quad_coef.domain().size()); + std::vector> eval_pts_data(m_quad_coef_alloc.domain().size()); ddc::ChunkSpan, ddc::DiscreteDomain> const - eval_pts(eval_pts_data.data(), m_quad_coef.domain()); - gl.compute_points_and_weights_on_mesh(eval_pts, m_quad_coef.span_view(), knots.span_cview()); + eval_pts(eval_pts_data.data(), m_quad_coef_alloc.domain()); + auto quad_coef_host = ddc::create_mirror_and_copy(m_quad_coef_alloc.span_view()); + gl.compute_points_and_weights_on_mesh(eval_pts, quad_coef_host.span_view(), knots.span_cview()); + ddc::deepcopy(m_quad_coef_alloc, quad_coef_host); ddc::init_discrete_space(eval_pts_data); @@ -65,9 +67,11 @@ void FemPeriodicPoissonSolver::build_matrix() positive_definite_symmetric, n_lower_diags); + auto quad_coef_host = ddc::create_mirror_and_copy(m_quad_coef_alloc.span_cview()); + // Fill the banded part of the matrix std::array derivs; - ddc::for_each(m_quad_coef.domain(), [&](ddc::DiscreteElement const ix) { + ddc::for_each(m_quad_coef_alloc.domain(), [&](ddc::DiscreteElement const ix) { ddc::Coordinate const coord = coord_from_quad_point(ddc::coordinate(ix)); ddc::DiscreteElement const jmin = ddc::discrete_space().eval_deriv(derivs, coord); @@ -77,7 +81,7 @@ void FemPeriodicPoissonSolver::build_matrix() int const k_idx = (k + jmin.uid()) % m_nbasis; double a_jk = m_fem_matrix->get_element(j_idx, k_idx); // Update element - a_jk = a_jk + derivs[j] * derivs[k] * m_quad_coef(ix); + a_jk += derivs[j] * derivs[k] * quad_coef_host(ix); m_fem_matrix->set_element(j_idx, k_idx, a_jk); } @@ -107,39 +111,50 @@ void FemPeriodicPoissonSolver::build_matrix() // Solve the Poisson equation //--------------------------------------------------------------------------- void FemPeriodicPoissonSolver::solve_matrix_system( - ddc::ChunkSpan const phi_spline_coef, - ddc::ChunkSpan const rho_spline_coef) const + DBSSpanX const phi_spline_coef, + DBSViewX const rho_spline_coef) const { - std::array values; + ddc::fill(phi_spline_coef, 0.0); int const rhs_size = m_nbasis + 1; - DSpan1D const phi_rhs(phi_spline_coef.data_handle(), rhs_size); - for (int i(0); i < rhs_size; ++i) { - phi_rhs(i) = 0.0; - } + int const nbasis_proxy = m_nbasis; + SplineXEvaluator_1d spline_x_evaluator_proxy = m_spline_x_evaluator; + ddc::ChunkSpan quad_coef = m_quad_coef_alloc.span_view(); // Fill phi_rhs(i) with \int rho(x) b_i(x) dx // Rk: phi_rhs no longer contains spline coefficients, but is the // RHS of the matrix equation - ddc::for_each(m_quad_coef.domain(), [&](ddc::DiscreteElement const ix) { - ddc::Coordinate const coord = coord_from_quad_point(ddc::coordinate(ix)); - ddc::DiscreteElement const jmin - = ddc::discrete_space().eval_basis(values, coord); - double const rho_val = m_spline_x_evaluator(coord, rho_spline_coef.span_cview()); - for (int j = 0; j < s_degree + 1; ++j) { - int const j_idx = (jmin.uid() + j) % m_nbasis; - phi_rhs(j_idx) = phi_rhs(j_idx) + rho_val * values[j] * m_quad_coef(ix); - } - }); + ddc::for_each( + ddc::policies::parallel_device, + m_quad_coef_alloc.domain(), + KOKKOS_LAMBDA(ddc::DiscreteElement const ix) { + ddc::Coordinate const coord = coord_from_quad_point(ddc::coordinate(ix)); + std::array values; + ddc::DiscreteElement const jmin + = ddc::discrete_space().eval_basis(values, coord); + double const rho_val + = spline_x_evaluator_proxy(coord, rho_spline_coef.span_cview()); + for (int j = 0; j < s_degree + 1; ++j) { + int const j_idx = (jmin.uid() + j) % nbasis_proxy; + Kokkos::atomic_fetch_add( + &phi_spline_coef(ddc::DiscreteElement(j_idx)), + rho_val * values[j] * quad_coef(ix)); + } + }); + + auto phi_spline_coef_host = ddc::create_mirror_and_copy(phi_spline_coef); + DSpan1D const phi_rhs_host(phi_spline_coef_host.data_handle(), rhs_size); // Solve the matrix equation to find the spline coefficients of phi - m_fem_matrix->solve_inplace(phi_rhs); + m_fem_matrix->solve_inplace(phi_rhs_host); // Copy the first d coefficients into the last d coefficients // These coefficients refer to the same BSplines which cross the boundaries for (int i = 0; i < s_degree; i++) { - phi_spline_coef(ddc::DiscreteElement(m_nbasis + i)) = phi_rhs(i); + phi_spline_coef_host(ddc::DiscreteElement(m_nbasis + i)) = phi_rhs_host(i); } + + ddc::deepcopy(phi_spline_coef, phi_spline_coef_host); } @@ -165,8 +180,6 @@ void FemPeriodicPoissonSolver::operator()( DViewSpXVx allfdistribu) const { Kokkos::Profiling::pushRegion("PoissonSolver"); - auto electrostatic_potential_host = ddc::create_mirror_and_copy(electrostatic_potential); - auto electric_field_host = ddc::create_mirror_and_copy(electric_field); assert(electrostatic_potential.domain() == ddc::get_domain(allfdistribu)); IDomainX const dom_x = electrostatic_potential.domain(); @@ -174,22 +187,24 @@ void FemPeriodicPoissonSolver::operator()( host_t rho_host(dom_x); DFieldX rho(dom_x); m_compute_rho(rho, allfdistribu); - ddc::deepcopy(rho_host, rho); // - ddc::Chunk rho_spline_coef(m_spline_x_builder.spline_domain()); - m_spline_x_builder(rho_spline_coef.span_view(), rho_host.span_cview()); - ddc::Chunk phi_spline_coef(m_spline_x_builder.spline_domain()); - solve_matrix_system(phi_spline_coef, rho_spline_coef); + ddc::Chunk rho_spline_coef(m_spline_x_builder.spline_domain(), ddc::DeviceAllocator()); + m_spline_x_builder(rho_spline_coef.span_view(), rho.span_cview()); + ddc::Chunk phi_spline_coef(m_spline_x_builder.spline_domain(), ddc::DeviceAllocator()); + solve_matrix_system(phi_spline_coef, rho_spline_coef.span_cview()); // - ddc::for_each(dom_x, [&](IndexX const ix) { - electrostatic_potential_host(ix) - = m_spline_x_evaluator(ddc::coordinate(ix), phi_spline_coef.span_cview()); - electric_field_host(ix) - = -m_spline_x_evaluator.deriv(ddc::coordinate(ix), phi_spline_coef.span_cview()); - }); - ddc::deepcopy(electrostatic_potential, electrostatic_potential_host); - ddc::deepcopy(electric_field, electric_field_host); + SplineXEvaluator_1d spline_x_evaluator_proxy = m_spline_x_evaluator; + ddc::ChunkSpan phi_spline_coef_view = phi_spline_coef.span_cview(); + ddc::for_each( + ddc::policies::parallel_device, + dom_x, + KOKKOS_LAMBDA(IndexX const ix) { + electrostatic_potential(ix) + = spline_x_evaluator_proxy(ddc::coordinate(ix), phi_spline_coef_view); + electric_field(ix) = -spline_x_evaluator_proxy + .deriv(ddc::coordinate(ix), phi_spline_coef_view); + }); Kokkos::Profiling::popRegion(); } diff --git a/src/geometryXVx/poisson/femperiodicpoissonsolver.hpp b/src/geometryXVx/poisson/femperiodicpoissonsolver.hpp index ab6a9c402..63ea0ff50 100644 --- a/src/geometryXVx/poisson/femperiodicpoissonsolver.hpp +++ b/src/geometryXVx/poisson/femperiodicpoissonsolver.hpp @@ -29,6 +29,7 @@ class FemPeriodicPoissonSolver : public IPoissonSolver private: using QMeshX = ddc::NonUniformPointSampling; + using DQFieldX = device_t>>; private: // Spline degree in x direction @@ -50,17 +51,19 @@ class FemPeriodicPoissonSolver : public IPoissonSolver // Number of cells in x direction int m_ncells; - ddc::Chunk> m_quad_coef; + DQFieldX m_quad_coef_alloc; std::unique_ptr m_fem_matrix; private: - static ddc::Coordinate quad_point_from_coord(ddc::Coordinate const& coord) + static KOKKOS_FUNCTION ddc::Coordinate quad_point_from_coord( + ddc::Coordinate const& coord) { return ddc::Coordinate(ddc::get(coord)); } - static ddc::Coordinate coord_from_quad_point(ddc::Coordinate const& coord) + static KOKKOS_FUNCTION ddc::Coordinate coord_from_quad_point( + ddc::Coordinate const& coord) { return ddc::Coordinate(ddc::get(coord)); } @@ -91,7 +94,12 @@ class FemPeriodicPoissonSolver : public IPoissonSolver private: void build_matrix(); - void solve_matrix_system( - ddc::ChunkSpan phi_spline_coef, - ddc::ChunkSpan rho_spline_coef) const; +public: + /** + * [SHOULD BE PRIVATE (Kokkos limitation)] + * + * @param[out] phi_spline_coef + * @param[in] rho_spline_coef + */ + void solve_matrix_system(DBSSpanX phi_spline_coef, DBSViewX rho_spline_coef) const; }; diff --git a/src/geometryXVx/poisson/fftpoissonsolver.hpp b/src/geometryXVx/poisson/fftpoissonsolver.hpp index ca805687d..1e0dff578 100644 --- a/src/geometryXVx/poisson/fftpoissonsolver.hpp +++ b/src/geometryXVx/poisson/fftpoissonsolver.hpp @@ -2,7 +2,6 @@ #pragma once -#include "electricfield.hpp" #include "ichargedensitycalculator.hpp" #include "ipoissonsolver.hpp" diff --git a/tests/geometryXVx/CMakeLists.txt b/tests/geometryXVx/CMakeLists.txt index 4882e9afb..1bd13389d 100644 --- a/tests/geometryXVx/CMakeLists.txt +++ b/tests/geometryXVx/CMakeLists.txt @@ -26,7 +26,6 @@ target_link_libraries(unit_tests_${GEOMETRY_VARIANT} GTest::gtest GTest::gmock paraconf::paraconf - sll::splines gslx::advection gslx::boltzmann_${GEOMETRY_VARIANT} gslx::initialization_${GEOMETRY_VARIANT}