From cfd3b7747a6506d43def1c87c03fd6149ddc27d0 Mon Sep 17 00:00:00 2001 From: haykh Date: Tue, 17 Jun 2025 08:46:39 -0400 Subject: [PATCH 1/3] new branch declaration From 6a36364e2837562eeae613489d65fe73d3373c93 Mon Sep 17 00:00:00 2001 From: haykh Date: Tue, 17 Jun 2025 16:11:37 -0400 Subject: [PATCH 2/3] passing flds_bc to kernel --- src/engines/grpic.hpp | 54 +++------------------------------- src/engines/srpic.hpp | 55 +++++++++++++++++------------------ src/kernels/tests/flds_bc.cpp | 13 ++++++--- 3 files changed, 39 insertions(+), 83 deletions(-) diff --git a/src/engines/grpic.hpp b/src/engines/grpic.hpp index 78966a25..168f2c73 100644 --- a/src/engines/grpic.hpp +++ b/src/engines/grpic.hpp @@ -647,7 +647,8 @@ namespace ntt { domain.mesh.metric, xg_edge, ds, - tags)); + tags, + domain.mesh.flds_bc())); Kokkos::parallel_for( "MatchBoundaries", CreateRangePolicy(range_min, range_max), @@ -657,7 +658,8 @@ namespace ntt { domain.mesh.metric, xg_edge, ds, - tags)); + tags, + domain.mesh.flds_bc())); } else { Kokkos::parallel_for( "AbsorbCurrents", @@ -672,54 +674,6 @@ namespace ntt { } } - // void AbsorbCurrentsIn(dir::direction_t direction, - // domain_t& domain, - // BCTags tags) { - // /** - // * absorbing currents at the boundaries - // */ - // // @TODO: also in MatchFieldsIn (need a better way) - // const auto ds_array = m_params.template get>( - // "grid.boundaries.match.ds"); - // const auto dim = in::x1; - // real_t xg_min, xg_max, xg_edge; - // real_t ds; - // ds = ds_array[(short)dim].second; - // xg_max = m_metadomain.mesh().extent(dim).second; - // xg_min = xg_max - ds; - // xg_edge = xg_max; - // - // boundaries_t box; - // boundaries_t incl_ghosts; - // for (unsigned short d { 0 }; d < M::Dim; ++d) { - // if (d == static_cast(dim)) { - // box.push_back({ xg_min, xg_max }); - // incl_ghosts.push_back({ false, true }); - // } else { - // box.push_back(Range::All); - // incl_ghosts.push_back({ true, true }); - // } - // } - // if (not domain.mesh.Intersects(box)) { - // return; - // } - // const auto intersect_range = domain.mesh.ExtentToRange(box, incl_ghosts); - // tuple_t range_min { 0 }; - // tuple_t range_max { 0 }; - // - // for (unsigned short d { 0 }; d < M::Dim; ++d) { - // range_min[d] = intersect_range[d].first; - // range_max[d] = intersect_range[d].second; - // } - // Kokkos::parallel_for( - // "AbsorbCurrentsGR", - // CreateRangePolicy(range_min, range_max), - // kernel::bc::AbsorbCurrentsGR_kernel(domain.fields.cur0, - // domain.mesh.metric, - // xg_edge, - // ds)); - // } - void HorizonFieldsIn(dir::direction_t direction, domain_t& domain, BCTags tags, diff --git a/src/engines/srpic.hpp b/src/engines/srpic.hpp index 0cdedb9b..6b6a5203 100644 --- a/src/engines/srpic.hpp +++ b/src/engines/srpic.hpp @@ -492,31 +492,20 @@ namespace ntt { species.npart(), (double)species.charge()), HERE); + // clang-format off Kokkos::parallel_for("CurrentsDeposit", species.rangeActiveParticles(), kernel::DepositCurrents_kernel( scatter_cur, - species.i1, - species.i2, - species.i3, - species.i1_prev, - species.i2_prev, - species.i3_prev, - species.dx1, - species.dx2, - species.dx3, - species.dx1_prev, - species.dx2_prev, - species.dx3_prev, - species.ux1, - species.ux2, - species.ux3, - species.phi, - species.weight, - species.tag, + species.i1, species.i2, species.i3, + species.i1_prev, species.i2_prev, species.i3_prev, + species.dx1, species.dx2, species.dx3, + species.dx1_prev, species.dx2_prev, species.dx3_prev, + species.ux1, species.ux2, species.ux3, + species.phi, species.weight, species.tag, domain.mesh.metric, - (real_t)(species.charge()), - dt)); + (real_t)(species.charge()), dt)); + // clang-format on } Kokkos::Experimental::contribute(domain.fields.cur, scatter_cur); } @@ -685,6 +674,7 @@ namespace ntt { traits::has_member::value) { auto match_fields = m_pgen.MatchFields(time); call_match_fields(domain.fields.em, + domain.mesh.flds_bc(), match_fields, domain.mesh.metric, xg_edge, @@ -696,6 +686,7 @@ namespace ntt { traits::has_member::value) { auto match_fields = m_pgen.MatchFieldsInX1(time); call_match_fields(domain.fields.em, + domain.mesh.flds_bc(), match_fields, domain.mesh.metric, xg_edge, @@ -710,6 +701,7 @@ namespace ntt { traits::has_member::value) { auto match_fields = m_pgen.MatchFields(time); call_match_fields(domain.fields.em, + domain.mesh.flds_bc(), match_fields, domain.mesh.metric, xg_edge, @@ -721,6 +713,7 @@ namespace ntt { traits::has_member::value) { auto match_fields = m_pgen.MatchFieldsInX2(time); call_match_fields(domain.fields.em, + domain.mesh.flds_bc(), match_fields, domain.mesh.metric, xg_edge, @@ -738,6 +731,7 @@ namespace ntt { traits::has_member::value) { auto match_fields = m_pgen.MatchFields(time); call_match_fields(domain.fields.em, + domain.mesh.flds_bc(), match_fields, domain.mesh.metric, xg_edge, @@ -749,6 +743,7 @@ namespace ntt { traits::has_member::value) { auto match_fields = m_pgen.MatchFieldsInX3(time); call_match_fields(domain.fields.em, + domain.mesh.flds_bc(), match_fields, domain.mesh.metric, xg_edge, @@ -1463,14 +1458,15 @@ namespace ntt { } template - void call_match_fields(ndfield_t& fields, - const T& match_fields, - const M& metric, - real_t xg_edge, - real_t ds, - BCTags tags, - tuple_t& range_min, - tuple_t& range_max) { + void call_match_fields(ndfield_t& fields, + const boundaries_t& boundaries, + const T& match_fields, + const M& metric, + real_t xg_edge, + real_t ds, + BCTags tags, + tuple_t& range_min, + tuple_t& range_max) { Kokkos::parallel_for( "MatchFields", CreateRangePolicy(range_min, range_max), @@ -1479,7 +1475,8 @@ namespace ntt { metric, xg_edge, ds, - tags)); + tags, + boundaries)); } }; diff --git a/src/kernels/tests/flds_bc.cpp b/src/kernels/tests/flds_bc.cpp index c5675cde..5de006d1 100644 --- a/src/kernels/tests/flds_bc.cpp +++ b/src/kernels/tests/flds_bc.cpp @@ -89,9 +89,13 @@ void testFldsBCs(const std::vector& res) { { res[0] + 2 * N_GHOSTS, res[1] + N_GHOSTS, res[2] + N_GHOSTS }); } - const auto xg_edge = (real_t)(sx[0].second); - const auto dx_abs = (real_t)(res[0] / 10.0); - + const auto xg_edge = (real_t)(sx[0].second); + const auto dx_abs = (real_t)(res[0] / 10.0); + boundaries_t flds_bc { + { FldsBC::PERIODIC, FldsBC::PERIODIC }, + { FldsBC::PERIODIC, FldsBC::PERIODIC }, + { FldsBC::PERIODIC, FldsBC::PERIODIC } + }; Kokkos::parallel_for( "MatchBoundaries_kernel", range, @@ -101,7 +105,8 @@ void testFldsBCs(const std::vector& res) { metric, xg_edge, dx_abs, - BC::E | BC::B)); + BC::E | BC::B, + flds_bc)); if constexpr (D == Dim::_1D) { Kokkos::parallel_for( From b3e72bbe41eae0a2313583a58fbd9ffac1652134 Mon Sep 17 00:00:00 2001 From: haykh Date: Tue, 17 Jun 2025 16:12:23 -0400 Subject: [PATCH 3/3] special treatment of axis --- src/kernels/fields_bcs.hpp | 60 ++++++++++++++++++++++++++------------ 1 file changed, 41 insertions(+), 19 deletions(-) diff --git a/src/kernels/fields_bcs.hpp b/src/kernels/fields_bcs.hpp index 40dd3d13..04e63052 100644 --- a/src/kernels/fields_bcs.hpp +++ b/src/kernels/fields_bcs.hpp @@ -42,6 +42,7 @@ namespace kernel::bc { static_assert(M::is_metric, "M must be a metric class"); static_assert(static_cast(o) < static_cast(M::Dim), "Invalid component index"); + static constexpr auto D = M::Dim; static constexpr idx_t i = static_cast(o) + 1u; static constexpr bool defines_dx1 = traits::has_method::value; static constexpr bool defines_dx2 = traits::has_method::value; @@ -66,18 +67,30 @@ namespace kernel::bc { const real_t dx_abs; const BCTags tags; - MatchBoundaries_kernel(ndfield_t Fld, - const I& fset, - const M& metric, - real_t xg_edge, - real_t dx_abs, - BCTags tags) + ncells_t extent_2 { 0u }; + bool is_axis_i2min { false }, is_axis_i2max { false }; + + MatchBoundaries_kernel(ndfield_t Fld, + const I& fset, + const M& metric, + real_t xg_edge, + real_t dx_abs, + BCTags tags, + const boundaries_t& boundaries) : Fld { Fld } , fset { fset } , metric { metric } , xg_edge { xg_edge } , dx_abs { dx_abs } - , tags { tags } {} + , tags { tags } { + if constexpr ((M::CoordType != Coord::Cart) && + ((D == Dim::_2D) || (D == Dim::_3D))) { + raise::ErrorIf(boundaries.size() < 2, "boundaries defined incorrectly", HERE); + is_axis_i2min = (boundaries[1].first == FldsBC::AXIS); + is_axis_i2max = (boundaries[1].second == FldsBC::AXIS); + extent_2 = static_cast(Fld.extent(1)); + } + } Inline auto shape(const real_t& dx) const -> real_t { return math::tanh(dx * FOUR / dx_abs); @@ -281,14 +294,20 @@ namespace kernel::bc { metric.template convert({ i1_, i2_ }, x_Ph_00); if constexpr (defines_ex3 and S == SimEngine::SRPIC) { - Fld(i1, i2, em::ex3) = s * Fld(i1, i2, em::ex3) + - (ONE - s) * - metric.template transform<3, Idx::T, Idx::U>( - { i1_, i2_ }, - fset.ex3(x_Ph_00)); + Fld(i1, i2, em::ex3) = s * Fld(i1, i2, em::ex3); + if ((!is_axis_i2min or (i2 > N_GHOSTS)) and + (!is_axis_i2max or (i2 < extent_2 - N_GHOSTS))) { + Fld(i1, i2, em::ex3) += (ONE - s) * + metric.template transform<3, Idx::T, Idx::U>( + { i1_, i2_ }, + fset.ex3(x_Ph_00)); + } } else if constexpr (defines_dx3 and S == SimEngine::GRPIC) { - Fld(i1, i2, em::dx3) = s * Fld(i1, i2, em::dx3) + - (ONE - s) * fset.dx3(x_Ph_00); + Fld(i1, i2, em::dx3) = s * Fld(i1, i2, em::dx3); + if ((!is_axis_i2min or (i2 > N_GHOSTS)) and + (!is_axis_i2max or (i2 < extent_2 - N_GHOSTS))) { + Fld(i1, i2, em::dx3) += (ONE - s) * fset.dx3(x_Ph_00); + } } } } @@ -403,11 +422,14 @@ namespace kernel::bc { coord_t x_Ph_00H { ZERO }; metric.template convert({ i1_, i2_, i3_ + HALF }, x_Ph_00H); - Fld(i1, i2, i3, em::ex3) = - s * Fld(i1, i2, i3, em::ex3) + - (ONE - s) * metric.template transform<3, Idx::T, Idx::U>( - { i1_, i2_, i3_ + HALF }, - fset.ex3(x_Ph_00H)); + Fld(i1, i2, i3, em::ex3) = s * Fld(i1, i2, i3, em::ex3); + if ((!is_axis_i2min or (i2 > N_GHOSTS)) and + (!is_axis_i2max or (i2 < extent_2 - N_GHOSTS))) { + Fld(i1, i2, i3, em::ex3) += + (ONE - s) * metric.template transform<3, Idx::T, Idx::U>( + { i1_, i2_, i3_ + HALF }, + fset.ex3(x_Ph_00H)); + } } } }