From e1774afa2ec371ca683f1e147bf49403bbc032b6 Mon Sep 17 00:00:00 2001 From: haykh Date: Tue, 27 Jan 2026 16:39:48 -0500 Subject: [PATCH 1/9] empty commit From 556450870dee7e54bb7455ed1f21aa748dadbdae Mon Sep 17 00:00:00 2001 From: haykh Date: Thu, 29 Jan 2026 17:37:34 -0500 Subject: [PATCH 2/9] srpic engine modularized --- src/engines/engine.hpp | 7 + src/engines/srpic.hpp | 1499 --------------------------- src/engines/srpic/currents.h | 117 +++ src/engines/srpic/fields_bcs.h | 664 ++++++++++++ src/engines/srpic/fieldsolvers.h | 193 ++++ src/engines/srpic/particle_pusher.h | 219 ++++ src/engines/srpic/particles_bcs.h | 307 ++++++ src/engines/srpic/srpic.hpp | 199 ++++ src/engines/srpic/utils.h | 127 +++ src/entity.cpp | 2 +- 10 files changed, 1834 insertions(+), 1500 deletions(-) delete mode 100644 src/engines/srpic.hpp create mode 100644 src/engines/srpic/currents.h create mode 100644 src/engines/srpic/fields_bcs.h create mode 100644 src/engines/srpic/fieldsolvers.h create mode 100644 src/engines/srpic/particle_pusher.h create mode 100644 src/engines/srpic/particles_bcs.h create mode 100644 src/engines/srpic/srpic.hpp create mode 100644 src/engines/srpic/utils.h diff --git a/src/engines/engine.hpp b/src/engines/engine.hpp index 112ff236..31a2269a 100644 --- a/src/engines/engine.hpp +++ b/src/engines/engine.hpp @@ -120,6 +120,13 @@ namespace ntt { virtual void step_forward(timer::Timers&, Domain&) = 0; void run(); + + auto engineParams() const -> prm::Parameters { + auto parameters = prm::Parameters {}; + parameters.set("dt", static_cast(dt)); + parameters.set("time", static_cast(time)); + return parameters; + } }; template diff --git a/src/engines/srpic.hpp b/src/engines/srpic.hpp deleted file mode 100644 index 272faac8..00000000 --- a/src/engines/srpic.hpp +++ /dev/null @@ -1,1499 +0,0 @@ -/** - * @file engines/srpic.hpp - * @brief Simulation engien class which specialized on SRPIC - * @implements - * - ntt::SRPICEngine<> : ntt::Engine<> - * @cpp: - * - srpic.cpp - * @namespaces: - * - ntt:: - * @macros: - */ - -#ifndef ENGINES_SRPIC_SRPIC_H -#define ENGINES_SRPIC_SRPIC_H - -#include "enums.h" -#include "global.h" - -#include "arch/kokkos_aliases.h" -#include "arch/traits.h" -#include "utils/log.h" -#include "utils/numeric.h" -#include "utils/timer.h" -#include "utils/toml.h" - -#include "archetypes/energy_dist.h" -#include "archetypes/particle_injector.h" -#include "archetypes/spatial_dist.h" -#include "engines/traits.h" -#include "framework/domain/domain.h" -#include "framework/parameters/parameters.h" - -#include "engines/engine.hpp" -#include "kernels/ampere_mink.hpp" -#include "kernels/ampere_sr.hpp" -#include "kernels/currents_deposit.hpp" -#include "kernels/digital_filter.hpp" -#include "kernels/faraday_mink.hpp" -#include "kernels/faraday_sr.hpp" -#include "kernels/fields_bcs.hpp" -#include "kernels/particle_moments.hpp" -#include "kernels/particle_pusher_sr.hpp" -#include "pgen.hpp" - -#include -#include - -#include - -namespace ntt { - - template - requires traits::engine::IsCompatibleWithSRPICEngine - class SRPICEngine : public Engine { - - using base_t = Engine; - using pgen_t = user::PGen; - using domain_t = Domain; - // contents - using base_t::m_metadomain; - using base_t::m_params; - using base_t::m_pgen; - // methods - using base_t::init; - // variables - using base_t::dt; - using base_t::max_steps; - using base_t::runtime; - using base_t::step; - using base_t::time; - - public: - static constexpr auto S { SimEngine::SRPIC }; - - SRPICEngine(const SimulationParams& params) : base_t { params } {} - - ~SRPICEngine() = default; - - void step_forward(timer::Timers& timers, domain_t& dom) override { - const auto fieldsolver_enabled = m_params.template get( - "algorithms.fieldsolver.enable"); - const auto deposit_enabled = m_params.template get( - "algorithms.deposit.enable"); - const auto clear_interval = m_params.template get( - "particles.clear_interval"); - - if (step == 0) { - // communicate fields and apply BCs on the first timestep - m_metadomain.CommunicateFields(dom, Comm::B | Comm::E); - FieldBoundaries(dom, BC::B | BC::E); - ParticleInjector(dom); - } - - if (fieldsolver_enabled) { - timers.start("FieldSolver"); - Faraday(dom, HALF); - timers.stop("FieldSolver"); - - timers.start("Communications"); - m_metadomain.CommunicateFields(dom, Comm::B); - timers.stop("Communications"); - - timers.start("FieldBoundaries"); - FieldBoundaries(dom, BC::B); - timers.stop("FieldBoundaries"); - Kokkos::fence(); - } - - { - timers.start("ParticlePusher"); - ParticlePush(dom); - timers.stop("ParticlePusher"); - - if (deposit_enabled) { - timers.start("CurrentDeposit"); - Kokkos::deep_copy(dom.fields.cur, ZERO); - CurrentsDeposit(dom); - timers.stop("CurrentDeposit"); - - timers.start("Communications"); - m_metadomain.SynchronizeFields(dom, Comm::J); - m_metadomain.CommunicateFields(dom, Comm::J); - timers.stop("Communications"); - - timers.start("CurrentFiltering"); - CurrentsFilter(dom); - timers.stop("CurrentFiltering"); - } - - timers.start("Communications"); - m_metadomain.CommunicateParticles(dom); - timers.stop("Communications"); - } - - if (fieldsolver_enabled) { - timers.start("FieldSolver"); - Faraday(dom, HALF); - timers.stop("FieldSolver"); - - timers.start("Communications"); - m_metadomain.CommunicateFields(dom, Comm::B); - timers.stop("Communications"); - - timers.start("FieldBoundaries"); - FieldBoundaries(dom, BC::B); - timers.stop("FieldBoundaries"); - - timers.start("FieldSolver"); - Ampere(dom, ONE); - timers.stop("FieldSolver"); - - if (deposit_enabled) { - timers.start("FieldSolver"); - CurrentsAmpere(dom); - timers.stop("FieldSolver"); - } - - timers.start("Communications"); - m_metadomain.CommunicateFields(dom, Comm::E | Comm::J); - timers.stop("Communications"); - - timers.start("FieldBoundaries"); - FieldBoundaries(dom, BC::E); - timers.stop("FieldBoundaries"); - } - - { - timers.start("Injector"); - ParticleInjector(dom); - timers.stop("Injector"); - } - - if (clear_interval > 0 and step % clear_interval == 0 and step > 0) { - timers.start("PrtlClear"); - m_metadomain.RemoveDeadParticles(dom); - timers.stop("PrtlClear"); - } - } - - /* algorithm substeps --------------------------------------------------- */ - void Faraday(domain_t& domain, real_t fraction = ONE) { - logger::Checkpoint("Launching Faraday kernel", HERE); - const auto dT = fraction * - m_params.template get( - "algorithms.timestep.correction") * - dt; - if constexpr (M::CoordType == Coord::Cart) { - // minkowski case - const auto dx = math::sqrt(domain.mesh.metric.template h_<1, 1>({})); - const auto deltax = m_params.template get( - "algorithms.fieldsolver.delta_x"); - const auto deltay = m_params.template get( - "algorithms.fieldsolver.delta_y"); - const auto betaxy = m_params.template get( - "algorithms.fieldsolver.beta_xy"); - const auto betayx = m_params.template get( - "algorithms.fieldsolver.beta_yx"); - const auto deltaz = m_params.template get( - "algorithms.fieldsolver.delta_z"); - const auto betaxz = m_params.template get( - "algorithms.fieldsolver.beta_xz"); - const auto betazx = m_params.template get( - "algorithms.fieldsolver.beta_zx"); - const auto betayz = m_params.template get( - "algorithms.fieldsolver.beta_yz"); - const auto betazy = m_params.template get( - "algorithms.fieldsolver.beta_zy"); - real_t coeff1, coeff2; - if constexpr (M::Dim == Dim::_2D) { - coeff1 = dT / SQR(dx); - coeff2 = dT; - } else { - coeff1 = dT / dx; - coeff2 = ZERO; - } - Kokkos::parallel_for("Faraday", - domain.mesh.rangeActiveCells(), - kernel::mink::Faraday_kernel(domain.fields.em, - coeff1, - coeff2, - deltax, - deltay, - betaxy, - betayx, - deltaz, - betaxz, - betazx, - betayz, - betazy)); - } else { - Kokkos::parallel_for("Faraday", - domain.mesh.rangeActiveCells(), - kernel::sr::Faraday_kernel(domain.fields.em, - domain.mesh.metric, - dT, - domain.mesh.flds_bc())); - } - } - - void Ampere(domain_t& domain, real_t fraction = ONE) { - logger::Checkpoint("Launching Ampere kernel", HERE); - const auto dT = fraction * - m_params.template get( - "algorithms.timestep.correction") * - dt; - auto range = range_with_axis_BCs(domain); - if constexpr (M::CoordType == Coord::Cart) { - // minkowski case - const auto dx = math::sqrt(domain.mesh.metric.template h_<1, 1>({})); - real_t coeff1, coeff2; - if constexpr (M::Dim == Dim::_2D) { - coeff1 = dT / SQR(dx); - coeff2 = dT; - } else { - coeff1 = dT / dx; - coeff2 = ZERO; - } - - Kokkos::parallel_for( - "Ampere", - range, - kernel::mink::Ampere_kernel(domain.fields.em, coeff1, coeff2)); - } else { - const auto ni2 = domain.mesh.n_active(in::x2); - Kokkos::parallel_for("Ampere", - range, - kernel::sr::Ampere_kernel(domain.fields.em, - domain.mesh.metric, - dT, - ni2, - domain.mesh.flds_bc())); - } - } - - void ParticlePush(domain_t& domain) { - real_t gx1 { ZERO }, gx2 { ZERO }, gx3 { ZERO }, ds { ZERO }; - real_t x_surf { ZERO }; - bool has_atmosphere = false; - for (auto& direction : dir::Directions::orth) { - if (m_metadomain.mesh().prtl_bc_in(direction) == PrtlBC::ATMOSPHERE) { - raise::ErrorIf(has_atmosphere, - "Only one direction is allowed to have atm boundaries", - HERE); - has_atmosphere = true; - const auto g = m_params.template get( - "grid.boundaries.atmosphere.g"); - ds = m_params.template get("grid.boundaries.atmosphere.ds"); - const auto [sign, dim, xg_min, xg_max] = get_atm_extent(direction); - if (dim == in::x1) { - gx1 = sign > 0 ? g : -g; - gx2 = ZERO; - gx3 = ZERO; - } else if (dim == in::x2) { - gx1 = ZERO; - gx2 = sign > 0 ? g : -g; - gx3 = ZERO; - } else if (dim == in::x3) { - gx1 = ZERO; - gx2 = ZERO; - gx3 = sign > 0 ? g : -g; - } else { - raise::Error("Invalid dimension", HERE); - } - if (sign > 0) { - x_surf = xg_min; - } else { - x_surf = xg_max; - } - } - } - for (auto& species : domain.species) { - if ((species.pusher() == ParticlePusher::NONE) or (species.npart() == 0)) { - continue; - } - species.set_unsorted(); - logger::Checkpoint( - fmt::format("Launching particle pusher kernel for %d [%s] : %lu", - species.index(), - species.label().c_str(), - species.npart()), - HERE); - - kernel::sr::PusherParams pusher_params {}; - pusher_params.pusher_flags = species.pusher(); - pusher_params.radiative_drag_flags = species.radiative_drag_flags(); - pusher_params.mass = species.mass(); - pusher_params.charge = species.charge(); - pusher_params.time = time; - pusher_params.dt = dt; - pusher_params.omegaB0 = m_params.template get("scales.omegaB0"); - pusher_params.ni1 = domain.mesh.n_active(in::x1); - pusher_params.ni2 = domain.mesh.n_active(in::x2); - pusher_params.ni3 = domain.mesh.n_active(in::x3); - pusher_params.boundaries = domain.mesh.prtl_bc(); - - if (species.pusher() & ParticlePusher::GCA) { - pusher_params.gca_params.set( - "larmor_max", - m_params.template get("algorithms.gca.larmor_max")); - pusher_params.gca_params.set( - "e_ovr_b_max", - m_params.template get("algorithms.gca.e_ovr_b_max")); - } - - if (species.radiative_drag_flags() & RadiativeDrag::SYNCHROTRON) { - pusher_params.radiative_drag_params.set( - "synchrotron_gamma_rad", - m_params.template get( - "radiation.drag.synchrotron.gamma_rad")); - } - - if (species.radiative_drag_flags() & RadiativeDrag::COMPTON) { - pusher_params.radiative_drag_params.set( - "compton_gamma_rad", - m_params.template get("radiation.drag.compton.gamma_rad")); - } - - kernel::sr::PusherArrays pusher_arrays {}; - pusher_arrays.sp = species.index(); - pusher_arrays.i1 = species.i1; - pusher_arrays.i2 = species.i2; - pusher_arrays.i3 = species.i3; - pusher_arrays.i1_prev = species.i1_prev; - pusher_arrays.i2_prev = species.i2_prev; - pusher_arrays.i3_prev = species.i3_prev; - pusher_arrays.dx1 = species.dx1; - pusher_arrays.dx2 = species.dx2; - pusher_arrays.dx3 = species.dx3; - pusher_arrays.dx1_prev = species.dx1_prev; - pusher_arrays.dx2_prev = species.dx2_prev; - pusher_arrays.dx3_prev = species.dx3_prev; - pusher_arrays.ux1 = species.ux1; - pusher_arrays.ux2 = species.ux2; - pusher_arrays.ux3 = species.ux3; - pusher_arrays.phi = species.phi; - pusher_arrays.tag = species.tag; - - // toggle to indicate whether pgen defines the external force - bool has_extforce = false; - if constexpr (arch::traits::pgen::HasExtForce) { - has_extforce = true; - // toggle to indicate whether the ext force applies to current species - if ( - ::traits::has_member<::traits::species_t, decltype(pgen_t::ext_force)>::value) { - has_extforce &= std::find(m_pgen.ext_force.species.begin(), - m_pgen.ext_force.species.end(), - species.index()) != - m_pgen.ext_force.species.end(); - } - } - - pusher_params.ext_force = has_extforce; - - if (not has_atmosphere and not has_extforce) { - Kokkos::parallel_for("ParticlePusher", - species.rangeActiveParticles(), - kernel::sr::Pusher_kernel(pusher_params, - pusher_arrays, - domain.fields.em, - domain.mesh.metric)); - } else if (has_atmosphere and not has_extforce) { - const auto force = - kernel::sr::Force { - { gx1, gx2, gx3 }, - x_surf, - ds - }; - Kokkos::parallel_for( - "ParticlePusher", - species.rangeActiveParticles(), - kernel::sr::Pusher_kernel(pusher_params, - pusher_arrays, - domain.fields.em, - domain.mesh.metric, - force)); - } else if (not has_atmosphere and has_extforce) { - if constexpr (arch::traits::pgen::HasExtForce) { - const auto force = - kernel::sr::Force { - m_pgen.ext_force - }; - Kokkos::parallel_for( - "ParticlePusher", - species.rangeActiveParticles(), - kernel::sr::Pusher_kernel(pusher_params, - pusher_arrays, - domain.fields.em, - domain.mesh.metric, - force)); - } else { - raise::Error("External force not implemented", HERE); - } - } else { // has_atmosphere and has_extforce - if constexpr (arch::traits::pgen::HasExtForce) { - const auto force = - kernel::sr::Force { - m_pgen.ext_force, - { gx1, gx2, gx3 }, - x_surf, - ds - }; - Kokkos::parallel_for( - "ParticlePusher", - species.rangeActiveParticles(), - kernel::sr::Pusher_kernel(pusher_params, - pusher_arrays, - domain.fields.em, - domain.mesh.metric, - force)); - } else { - raise::Error("External force not implemented", HERE); - } - } - } - } - - void ParticleInjector(domain_t& domain, InjTags tags = Inj::None) { - for (auto& direction : dir::Directions::orth) { - if (m_metadomain.mesh().prtl_bc_in(direction) == PrtlBC::ATMOSPHERE) { - AtmosphereParticlesIn(direction, domain, tags); - } - } - } - - template - void deposit_with(const Particles& species, - const M& metric, - const scatter_ndfield_t& scatter_cur, - real_t dt) { - // 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, - metric, (real_t)(species.charge()), dt)); - // clang-format on - } - - void CurrentsDeposit(domain_t& domain) { - auto scatter_cur = Kokkos::Experimental::create_scatter_view( - domain.fields.cur); - for (auto& species : domain.species) { - if ((species.pusher() == ParticlePusher::NONE) or - (species.npart() == 0) or cmp::AlmostZero_host(species.charge())) { - continue; - } - logger::Checkpoint( - fmt::format("Launching currents deposit kernel for %d [%s] : %lu %f", - species.index(), - species.label().c_str(), - species.npart(), - (double)species.charge()), - HERE); - - deposit_with(species, domain.mesh.metric, scatter_cur, dt); - } - Kokkos::Experimental::contribute(domain.fields.cur, scatter_cur); - } - - void CurrentsAmpere(domain_t& domain) { - logger::Checkpoint("Launching Ampere kernel for adding currents", HERE); - const auto q0 = m_params.template get("scales.q0"); - const auto n0 = m_params.template get("scales.n0"); - const auto B0 = m_params.template get("scales.B0"); - if constexpr (M::CoordType == Coord::Cart) { - // minkowski case - const auto V0 = m_params.template get("scales.V0"); - const auto ppc0 = m_params.template get("particles.ppc0"); - const auto coeff = -dt * q0 / (B0 * V0); - if constexpr (arch::traits::pgen::HasExtCurrent) { - const std::vector xmin { domain.mesh.extent(in::x1).first, - domain.mesh.extent(in::x2).first, - domain.mesh.extent(in::x3).first }; - const auto ext_current = m_pgen.ext_current; - const auto dx = domain.mesh.metric.template sqrt_h_<1, 1>({}); - // clang-format off - Kokkos::parallel_for( - "Ampere", - domain.mesh.rangeActiveCells(), - kernel::mink::CurrentsAmpere_kernel( - domain.fields.em, domain.fields.cur, - coeff, ppc0, ext_current, xmin, dx)); - // clang-format on - } else { - Kokkos::parallel_for( - "Ampere", - domain.mesh.rangeActiveCells(), - kernel::mink::CurrentsAmpere_kernel(domain.fields.em, - domain.fields.cur, - coeff, - ppc0)); - } - } else { - // non-minkowski - const auto coeff = -dt * q0 * n0 / B0; - auto range = range_with_axis_BCs(domain); - const auto ni2 = domain.mesh.n_active(in::x2); - Kokkos::parallel_for( - "Ampere", - range, - kernel::sr::CurrentsAmpere_kernel(domain.fields.em, - domain.fields.cur, - domain.mesh.metric, - coeff, - ONE / n0, - ni2, - domain.mesh.flds_bc())); - } - } - - void CurrentsFilter(domain_t& domain) { - logger::Checkpoint("Launching currents filtering kernels", HERE); - auto range = range_with_axis_BCs(domain); - const auto nfilter = m_params.template get( - "algorithms.current_filters"); - tuple_t size; - if constexpr (M::Dim == Dim::_1D || M::Dim == Dim::_2D || M::Dim == Dim::_3D) { - size[0] = domain.mesh.n_active(in::x1); - } - if constexpr (M::Dim == Dim::_2D || M::Dim == Dim::_3D) { - size[1] = domain.mesh.n_active(in::x2); - } - if constexpr (M::Dim == Dim::_3D) { - size[2] = domain.mesh.n_active(in::x3); - } - // !TODO: this needs to be done more efficiently - for (auto i { 0u }; i < nfilter; ++i) { - Kokkos::deep_copy(domain.fields.buff, domain.fields.cur); - Kokkos::parallel_for("CurrentsFilter", - range, - kernel::DigitalFilter_kernel( - domain.fields.cur, - domain.fields.buff, - size, - domain.mesh.flds_bc())); - m_metadomain.CommunicateFields(domain, Comm::J); - } - } - - void FieldBoundaries(domain_t& domain, BCTags tags) { - for (auto& direction : dir::Directions::orth) { - if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::MATCH) { - MatchFieldsIn(direction, domain, tags); - } else if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::AXIS) { - if (domain.mesh.flds_bc_in(direction) == FldsBC::AXIS) { - AxisFieldsIn(direction, domain, tags); - } - } else if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::ATMOSPHERE) { - AtmosphereFieldsIn(direction, domain, tags); - } else if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::FIXED) { - if (domain.mesh.flds_bc_in(direction) == FldsBC::FIXED) { - FixedFieldsIn(direction, domain, tags); - } - } else if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::CONDUCTOR) { - if (domain.mesh.flds_bc_in(direction) == FldsBC::CONDUCTOR) { - PerfectConductorFieldsIn(direction, domain, tags); - } - } else if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::CUSTOM) { - if (domain.mesh.flds_bc_in(direction) == FldsBC::CUSTOM) { - CustomFieldsIn(direction, domain, tags); - } - } else if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::HORIZON) { - raise::Error("HORIZON BCs only applicable for GR", HERE); - } - } // loop over directions - } - - void MatchFieldsIn(dir::direction_t direction, - domain_t& domain, - BCTags tags) { - /** - * matching boundaries - */ - const auto ds_array = m_params.template get>( - "grid.boundaries.match.ds"); - const auto dim = direction.get_dim(); - real_t xg_min, xg_max, xg_edge; - auto sign = direction.get_sign(); - real_t ds; - if (sign > 0) { // + direction - ds = ds_array[(short)dim].second; - xg_max = m_metadomain.mesh().extent(dim).second; - xg_min = xg_max - ds; - xg_edge = xg_max; - } else { // - direction - ds = ds_array[(short)dim].first; - xg_min = m_metadomain.mesh().extent(dim).first; - xg_max = xg_min + ds; - xg_edge = xg_min; - } - boundaries_t box; - boundaries_t incl_ghosts; - for (dim_t d { 0 }; d < M::Dim; ++d) { - if (d == static_cast(dim)) { - box.push_back({ xg_min, xg_max }); - if (sign > 0) { - incl_ghosts.push_back({ false, true }); - } else { - incl_ghosts.push_back({ true, false }); - } - } 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 (auto d { 0u }; d < M::Dim; ++d) { - range_min[d] = intersect_range[d].first; - range_max[d] = intersect_range[d].second; - } - - if (dim == in::x1) { - if constexpr (arch::traits::pgen::HasMatchFields) { - auto match_fields = m_pgen.MatchFields(time); - call_match_fields(domain.fields.em, - domain.mesh.flds_bc(), - match_fields, - domain.mesh.metric, - xg_edge, - ds, - tags, - range_min, - range_max); - } else if constexpr (arch::traits::pgen::HasMatchFieldsInX1) { - auto match_fields = m_pgen.MatchFieldsInX1(time); - call_match_fields(domain.fields.em, - domain.mesh.flds_bc(), - match_fields, - domain.mesh.metric, - xg_edge, - ds, - tags, - range_min, - range_max); - } - } else if (dim == in::x2) { - if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - if constexpr (arch::traits::pgen::HasMatchFields) { - auto match_fields = m_pgen.MatchFields(time); - call_match_fields(domain.fields.em, - domain.mesh.flds_bc(), - match_fields, - domain.mesh.metric, - xg_edge, - ds, - tags, - range_min, - range_max); - } else if constexpr (arch::traits::pgen::HasMatchFieldsInX2) { - auto match_fields = m_pgen.MatchFieldsInX2(time); - call_match_fields(domain.fields.em, - domain.mesh.flds_bc(), - match_fields, - domain.mesh.metric, - xg_edge, - ds, - tags, - range_min, - range_max); - } - } else { - raise::Error("Invalid dimension", HERE); - } - } else if (dim == in::x3) { - if constexpr (M::Dim == Dim::_3D) { - if constexpr (arch::traits::pgen::HasMatchFields) { - auto match_fields = m_pgen.MatchFields(time); - call_match_fields(domain.fields.em, - domain.mesh.flds_bc(), - match_fields, - domain.mesh.metric, - xg_edge, - ds, - tags, - range_min, - range_max); - } else if constexpr (arch::traits::pgen::HasMatchFieldsInX3) { - auto match_fields = m_pgen.MatchFieldsInX3(time); - call_match_fields(domain.fields.em, - domain.mesh.flds_bc(), - match_fields, - domain.mesh.metric, - xg_edge, - ds, - tags, - range_min, - range_max); - } - } - } else { - raise::Error("Invalid dimension", HERE); - } - } - - void AxisFieldsIn(dir::direction_t direction, - domain_t& domain, - BCTags tags) { - /** - * axis boundaries - */ - if constexpr (M::CoordType != Coord::Cart) { - raise::ErrorIf(direction.get_dim() != in::x2, - "Invalid axis direction, should be x2", - HERE); - const auto i2_min = domain.mesh.i_min(in::x2); - const auto i2_max = domain.mesh.i_max(in::x2); - if (direction.get_sign() < 0) { - Kokkos::parallel_for( - "AxisBCFields", - domain.mesh.n_all(in::x1), - kernel::bc::AxisBoundaries_kernel(domain.fields.em, - i2_min, - tags)); - } else { - Kokkos::parallel_for( - "AxisBCFields", - domain.mesh.n_all(in::x1), - kernel::bc::AxisBoundaries_kernel(domain.fields.em, - i2_max, - tags)); - } - } else { - (void)direction; - (void)domain; - (void)tags; - raise::Error("Invalid coordinate type for axis BCs", HERE); - } - } - - void FixedFieldsIn(dir::direction_t direction, - domain_t& domain, - BCTags tags) { - /** - * fixed field boundaries - */ - const auto sign = direction.get_sign(); - const auto dim = direction.get_dim(); - raise::ErrorIf(dim != in::x1 and M::CoordType != Coord::Cart, - "Fixed BCs only implemented for x1 in " - "non-cartesian coordinates", - HERE); - em normal_b_comp, tang_e_comp1, tang_e_comp2; - if (dim == in::x1) { - normal_b_comp = em::bx1; - tang_e_comp1 = em::ex2; - tang_e_comp2 = em::ex3; - } else if (dim == in::x2) { - normal_b_comp = em::bx2; - tang_e_comp1 = em::ex1; - tang_e_comp2 = em::ex3; - } else if (dim == in::x3) { - normal_b_comp = em::bx3; - tang_e_comp1 = em::ex1; - tang_e_comp2 = em::ex2; - } else { - raise::Error("Invalid dimension", HERE); - } - std::vector xi_min, xi_max; - const std::vector all_dirs { in::x1, in::x2, in::x3 }; - for (dim_t d { 0u }; d < M::Dim; ++d) { - const auto dd = all_dirs[d]; - if (dim == dd) { - if (sign > 0) { // + direction - xi_min.push_back(domain.mesh.n_all(dd) - N_GHOSTS); - xi_max.push_back(domain.mesh.n_all(dd)); - } else { // - direction - xi_min.push_back(0); - xi_max.push_back(N_GHOSTS); - } - } else { - xi_min.push_back(0); - xi_max.push_back(domain.mesh.n_all(dd)); - } - } - raise::ErrorIf(xi_min.size() != xi_max.size() or - xi_min.size() != static_cast(M::Dim), - "Invalid range size", - HERE); - std::vector comps; - if (tags & BC::E) { - comps.push_back(tang_e_comp1); - comps.push_back(tang_e_comp2); - } - if (tags & BC::B) { - comps.push_back(normal_b_comp); - } - if constexpr (arch::traits::pgen::HasFixFieldsConst) { - for (const auto& comp : comps) { - auto value = ZERO; - bool shouldset = false; - // if fix field function present, read from it - const auto newset = m_pgen.FixFieldsConst( - (bc_in)(sign * ((short)dim + 1)), - (em)comp); - value = newset.first; - shouldset = newset.second; - if (shouldset) { - if constexpr (M::Dim == Dim::_1D) { - Kokkos::deep_copy( - Kokkos::subview(domain.fields.em, - std::make_pair(xi_min[0], xi_max[0]), - comp), - value); - } else if constexpr (M::Dim == Dim::_2D) { - Kokkos::deep_copy( - Kokkos::subview(domain.fields.em, - std::make_pair(xi_min[0], xi_max[0]), - std::make_pair(xi_min[1], xi_max[1]), - comp), - value); - } else if constexpr (M::Dim == Dim::_3D) { - Kokkos::deep_copy( - Kokkos::subview(domain.fields.em, - std::make_pair(xi_min[0], xi_max[0]), - std::make_pair(xi_min[1], xi_max[1]), - std::make_pair(xi_min[2], xi_max[2]), - comp), - value); - } else { - raise::Error("Invalid dimension", HERE); - } - } - } - } else { - (void)direction; - (void)domain; - (void)tags; - raise::Error("Fixed fields not present (both const and non-const)", HERE); - } - } - - void PerfectConductorFieldsIn(dir::direction_t direction, - domain_t& domain, - BCTags tags) { - /** - * perfect conductor field boundaries - */ - if constexpr (M::CoordType != Coord::Cart) { - (void)direction; - (void)domain; - (void)tags; - raise::Error( - "Perfect conductor BCs only applicable to cartesian coordinates", - HERE); - } else { - const auto sign = direction.get_sign(); - const auto dim = direction.get_dim(); - - std::vector xi_min, xi_max; - - const std::vector all_dirs { in::x1, in::x2, in::x3 }; - - for (auto d { 0u }; d < M::Dim; ++d) { - const auto dd = all_dirs[d]; - if (dim == dd) { - xi_min.push_back(0); - xi_max.push_back((sign < 0) ? (N_GHOSTS + 1) : N_GHOSTS); - } else { - xi_min.push_back(0); - xi_max.push_back(domain.mesh.n_all(dd)); - } - } - raise::ErrorIf(xi_min.size() != xi_max.size() or - xi_min.size() != static_cast(M::Dim), - "Invalid range size", - HERE); - - range_t range; - if constexpr (M::Dim == Dim::_1D) { - range = CreateRangePolicy({ xi_min[0] }, { xi_max[0] }); - } else if constexpr (M::Dim == Dim::_2D) { - range = CreateRangePolicy({ xi_min[0], xi_min[1] }, - { xi_max[0], xi_max[1] }); - } else if constexpr (M::Dim == Dim::_3D) { - range = CreateRangePolicy({ xi_min[0], xi_min[1], xi_min[2] }, - { xi_max[0], xi_max[1], xi_max[2] }); - } else { - raise::Error("Invalid dimension", HERE); - } - std::size_t i_edge; - if (sign > 0) { - i_edge = domain.mesh.i_max(dim); - } else { - i_edge = domain.mesh.i_min(dim); - } - - if (dim == in::x1) { - if (sign > 0) { - Kokkos::parallel_for( - "ConductorFields", - range, - kernel::bc::ConductorBoundaries_kernel( - domain.fields.em, - i_edge, - tags)); - } else { - Kokkos::parallel_for( - "ConductorFields", - range, - kernel::bc::ConductorBoundaries_kernel( - domain.fields.em, - i_edge, - tags)); - } - } else if (dim == in::x2) { - if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - if (sign > 0) { - Kokkos::parallel_for( - "ConductorFields", - range, - kernel::bc::ConductorBoundaries_kernel( - domain.fields.em, - i_edge, - tags)); - } else { - Kokkos::parallel_for( - "ConductorFields", - range, - kernel::bc::ConductorBoundaries_kernel( - domain.fields.em, - i_edge, - tags)); - } - } else { - raise::Error("Invalid dimension", HERE); - } - } else { - if constexpr (M::Dim == Dim::_3D) { - if (sign > 0) { - Kokkos::parallel_for( - "ConductorFields", - range, - kernel::bc::ConductorBoundaries_kernel( - domain.fields.em, - i_edge, - tags)); - } else { - Kokkos::parallel_for( - "ConductorFields", - range, - kernel::bc::ConductorBoundaries_kernel( - domain.fields.em, - i_edge, - tags)); - } - } else { - raise::Error("Invalid dimension", HERE); - } - } - } - } - - void AtmosphereFieldsIn(dir::direction_t direction, - domain_t& domain, - BCTags tags) { - /** - * atmosphere field boundaries - */ - if constexpr (arch::traits::pgen::HasAtmFields) { - const auto [sign, dim, xg_min, xg_max] = get_atm_extent(direction); - const auto dd = static_cast(dim); - boundaries_t box; - boundaries_t incl_ghosts; - for (auto d { 0u }; d < M::Dim; ++d) { - if (d == dd) { - box.push_back({ xg_min, xg_max }); - if (sign > 0) { - incl_ghosts.push_back({ false, true }); - } else { - incl_ghosts.push_back({ true, false }); - } - } 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 (auto d { 0u }; d < M::Dim; ++d) { - range_min[d] = intersect_range[d].first; - range_max[d] = intersect_range[d].second; - } - auto atm_fields = m_pgen.AtmFields(time); - std::size_t il_edge; - if (sign > 0) { - il_edge = range_min[dd] - N_GHOSTS; - } else { - il_edge = range_max[dd] - 1 - N_GHOSTS; - } - const auto range = CreateRangePolicy(range_min, range_max); - if (dim == in::x1) { - if (sign > 0) { - Kokkos::parallel_for( - "AtmosphereBCFields", - range, - kernel::bc::EnforcedBoundaries_kernel( - domain.fields.em, - atm_fields, - domain.mesh.metric, - il_edge, - tags)); - } else { - Kokkos::parallel_for( - "AtmosphereBCFields", - range, - kernel::bc::EnforcedBoundaries_kernel( - domain.fields.em, - atm_fields, - domain.mesh.metric, - il_edge, - tags)); - } - } else if (dim == in::x2) { - if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - if (sign > 0) { - Kokkos::parallel_for( - "AtmosphereBCFields", - range, - kernel::bc::EnforcedBoundaries_kernel( - domain.fields.em, - atm_fields, - domain.mesh.metric, - il_edge, - tags)); - } else { - Kokkos::parallel_for( - "AtmosphereBCFields", - range, - kernel::bc::EnforcedBoundaries_kernel( - domain.fields.em, - atm_fields, - domain.mesh.metric, - il_edge, - tags)); - } - } else { - raise::Error("Invalid dimension", HERE); - } - } else if (dim == in::x3) { - if constexpr (M::Dim == Dim::_3D) { - if (sign > 0) { - Kokkos::parallel_for( - "AtmosphereBCFields", - range, - kernel::bc::EnforcedBoundaries_kernel( - domain.fields.em, - atm_fields, - domain.mesh.metric, - il_edge, - tags)); - } else { - Kokkos::parallel_for( - "AtmosphereBCFields", - range, - kernel::bc::EnforcedBoundaries_kernel( - domain.fields.em, - atm_fields, - domain.mesh.metric, - il_edge, - tags)); - } - } else { - raise::Error("Invalid dimension", HERE); - } - } else { - raise::Error("Invalid dimension", HERE); - } - } else { - (void)direction; - (void)domain; - (void)tags; - raise::Error("Atm fields not implemented in PGEN for atmosphere BCs", HERE); - } - } - - void CustomFieldsIn(dir::direction_t direction, - domain_t& domain, - BCTags tags) { - (void)direction; - (void)domain; - (void)tags; - raise::Error("Custom boundaries not implemented", HERE); - // if constexpr ( - // traits::has_member::value) { - // const auto [box, custom_fields] = m_pgen.CustomFields(time); - // if (domain.mesh.Intersects(box)) { - // } - // - // } else { - // raise::Error("Custom boundaries not implemented", HERE); - // } - } - - void AtmosphereParticlesIn(const dir::direction_t& direction, - domain_t& domain, - InjTags tags) { - const auto [sign, dim, xg_min, xg_max] = get_atm_extent(direction); - - const auto x_surf = sign > 0 ? xg_min : xg_max; - const auto ds = m_params.template get( - "grid.boundaries.atmosphere.ds"); - const auto temp = m_params.template get( - "grid.boundaries.atmosphere.temperature"); - const auto height = m_params.template get( - "grid.boundaries.atmosphere.height"); - const auto species = m_params.template get>( - "grid.boundaries.atmosphere.species"); - const auto nmax = m_params.template get( - "grid.boundaries.atmosphere.density"); - - Kokkos::deep_copy(domain.fields.bckp, ZERO); - auto scatter_bckp = Kokkos::Experimental::create_scatter_view( - domain.fields.bckp); - const auto use_weights = M::CoordType != Coord::Cart; - const auto ni2 = domain.mesh.n_active(in::x2); - const auto inv_n0 = ONE / m_params.template get("scales.n0"); - - // compute the density of the two species - if (tags & Inj::AssumeEmpty) { - if constexpr (M::Dim == Dim::_1D) { - Kokkos::deep_copy( - Kokkos::subview(domain.fields.bckp, Kokkos::ALL, std::make_pair(0, 1)), - ZERO); - } else if constexpr (M::Dim == Dim::_2D) { - Kokkos::deep_copy(Kokkos::subview(domain.fields.bckp, - Kokkos::ALL, - Kokkos::ALL, - std::make_pair(0, 1)), - ZERO); - } else if constexpr (M::Dim == Dim::_3D) { - Kokkos::deep_copy(Kokkos::subview(domain.fields.bckp, - Kokkos::ALL, - Kokkos::ALL, - Kokkos::ALL, - std::make_pair(0, 1)), - ZERO); - } - } else { - for (const auto& sp : - std::vector { species.first, species.second }) { - auto& prtl_spec = domain.species[sp - 1]; - if (prtl_spec.npart() == 0) { - continue; - } - // clang-format off - Kokkos::parallel_for( - "ComputeMoments", - prtl_spec.rangeActiveParticles(), - kernel::ParticleMoments_kernel( - {}, scatter_bckp, 0, - prtl_spec.i1, prtl_spec.i2, prtl_spec.i3, - prtl_spec.dx1, prtl_spec.dx2, prtl_spec.dx3, - prtl_spec.ux1, prtl_spec.ux2, prtl_spec.ux3, - prtl_spec.phi, prtl_spec.weight, prtl_spec.tag, - prtl_spec.mass(), prtl_spec.charge(), - use_weights, - domain.mesh.metric, domain.mesh.flds_bc(), - ni2, inv_n0, 0)); - // clang-format on - prtl_spec.set_unsorted(); - } - Kokkos::Experimental::contribute(domain.fields.bckp, scatter_bckp); - m_metadomain.SynchronizeFields(domain, Comm::Bckp, { 0, 1 }); - } - - const auto maxwellian = arch::Maxwellian { domain.mesh.metric, - domain.random_pool, - temp }; - - if (dim == in::x1) { - if (sign > 0) { - auto target_density = - arch::AtmosphereDensityProfile { - nmax, - height, - x_surf, - ds - }; - const auto spatial_dist = arch::Replenish { - domain.mesh.metric, - domain.fields.bckp, - 0, - target_density, - nmax - }; - arch::InjectNonUniform( - m_params, - domain, - { species.first, species.second }, - { maxwellian, maxwellian }, - spatial_dist, - nmax, - use_weights); - } else { - auto target_density = - arch::AtmosphereDensityProfile { - nmax, - height, - x_surf, - ds - }; - const auto spatial_dist = arch::Replenish { - domain.mesh.metric, - domain.fields.bckp, - 0, - target_density, - nmax - }; - arch::InjectNonUniform( - m_params, - domain, - { species.first, species.second }, - { maxwellian, maxwellian }, - spatial_dist, - nmax, - use_weights); - } - } else if (dim == in::x2) { - if (sign > 0) { - auto target_density = - arch::AtmosphereDensityProfile { - nmax, - height, - x_surf, - ds - }; - const auto spatial_dist = arch::Replenish { - domain.mesh.metric, - domain.fields.bckp, - 0, - target_density, - nmax - }; - arch::InjectNonUniform( - m_params, - domain, - { species.first, species.second }, - { maxwellian, maxwellian }, - spatial_dist, - nmax, - use_weights); - } else { - auto target_density = - arch::AtmosphereDensityProfile { - nmax, - height, - x_surf, - ds - }; - const auto spatial_dist = arch::Replenish { - domain.mesh.metric, - domain.fields.bckp, - 0, - target_density, - nmax - }; - arch::InjectNonUniform( - m_params, - domain, - { species.first, species.second }, - { maxwellian, maxwellian }, - spatial_dist, - nmax, - use_weights); - } - } else if (dim == in::x3) { - if (sign > 0) { - auto target_density = - arch::AtmosphereDensityProfile { - nmax, - height, - x_surf, - ds - }; - const auto spatial_dist = arch::Replenish { - domain.mesh.metric, - domain.fields.bckp, - 0, - target_density, - nmax - }; - arch::InjectNonUniform( - m_params, - domain, - { species.first, species.second }, - { maxwellian, maxwellian }, - spatial_dist, - nmax, - use_weights); - } else { - auto target_density = - arch::AtmosphereDensityProfile { - nmax, - height, - x_surf, - ds - }; - const auto spatial_dist = arch::Replenish { - domain.mesh.metric, - domain.fields.bckp, - 0, - target_density, - nmax - }; - arch::InjectNonUniform( - m_params, - domain, - { species.first, species.second }, - { maxwellian, maxwellian }, - spatial_dist, - nmax, - use_weights); - } - } else { - raise::Error("Invalid dimension", HERE); - } - return; - } - - private: - /** - * @brief Get the buffer region of the atmosphere and the direction - * @param direction direction in which the atmosphere is applied - * @return tuple: [sign of the direction, the direction (as in::), the min and max extent - * @note xg_min and xg_max are the extents where the fields are set, not the atmosphere itself - * @note i.e. - * - * fields set particles injected - * ghost zone | | - * v v v - * |....|...........|*******************..... -> x1 - * ^ ^ - * xg_min xg_max - * | | | - * |<-- buffer -->|<-- atmosphere -->| - * - * in this case the function returns { -1, in::x1, xg_min, xg_max } - */ - auto get_atm_extent(dir::direction_t direction) const - -> std::tuple { - const auto sign = direction.get_sign(); - const auto dim = direction.get_dim(); - const auto min_buff = m_params.template get( - "algorithms.current_filters") + - 2; - const auto buffer_ncells = min_buff > 5 ? min_buff : 5; - if (M::CoordType != Coord::Cart and (dim != in::x1 or sign > 0)) { - raise::Error("For non-cartesian coordinates atmosphere BCs is " - "possible only in -x1 (@ rmin)", - HERE); - } - real_t xg_min { ZERO }, xg_max { ZERO }; - ncells_t ig_min, ig_max; - if (sign > 0) { // + direction - ig_min = m_metadomain.mesh().n_active(dim) - buffer_ncells; - ig_max = m_metadomain.mesh().n_active(dim); - } else { // - direction - ig_min = 0; - ig_max = buffer_ncells; - } - - if (dim == in::x1) { - xg_min = m_metadomain.mesh().metric.template convert<1, Crd::Cd, Crd::Ph>( - static_cast(ig_min)); - xg_max = m_metadomain.mesh().metric.template convert<1, Crd::Cd, Crd::Ph>( - static_cast(ig_max)); - } else if (dim == in::x2) { - if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - xg_min = m_metadomain.mesh().metric.template convert<2, Crd::Cd, Crd::Ph>( - static_cast(ig_min)); - xg_max = m_metadomain.mesh().metric.template convert<2, Crd::Cd, Crd::Ph>( - static_cast(ig_max)); - } else { - raise::Error("Invalid dimension", HERE); - } - } else if (dim == in::x3) { - if constexpr (M::Dim == Dim::_3D) { - xg_min = m_metadomain.mesh().metric.template convert<3, Crd::Cd, Crd::Ph>( - static_cast(ig_min)); - xg_max = m_metadomain.mesh().metric.template convert<3, Crd::Cd, Crd::Ph>( - static_cast(ig_max)); - } else { - raise::Error("Invalid dimension", HERE); - } - } else { - raise::Error("Invalid dimension", HERE); - } - return { sign, dim, xg_min, xg_max }; - } - - auto range_with_axis_BCs(const domain_t& domain) -> range_t { - auto range = domain.mesh.rangeActiveCells(); - if constexpr (M::CoordType != Coord::Cart) { - /** - * @brief taking one extra cell in the x2 direction if AXIS BCs - */ - if constexpr (M::Dim == Dim::_2D) { - if (domain.mesh.flds_bc_in({ 0, +1 }) == FldsBC::AXIS) { - range = CreateRangePolicy( - { domain.mesh.i_min(in::x1), domain.mesh.i_min(in::x2) }, - { domain.mesh.i_max(in::x1), domain.mesh.i_max(in::x2) + 1 }); - } - } else if constexpr (M::Dim == Dim::_3D) { - if (domain.mesh.flds_bc_in({ 0, +1, 0 }) == FldsBC::AXIS) { - range = CreateRangePolicy({ domain.mesh.i_min(in::x1), - domain.mesh.i_min(in::x2), - domain.mesh.i_min(in::x3) }, - { domain.mesh.i_max(in::x1), - domain.mesh.i_max(in::x2) + 1, - domain.mesh.i_max(in::x3) }); - } - } - } - return range; - } - - template - 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), - kernel::bc::MatchBoundaries_kernel(fields, - match_fields, - metric, - xg_edge, - ds, - tags, - boundaries)); - } - }; - -} // namespace ntt - -#endif // ENGINES_SRPIC_SRPIC_H diff --git a/src/engines/srpic/currents.h b/src/engines/srpic/currents.h new file mode 100644 index 00000000..28b8f9ab --- /dev/null +++ b/src/engines/srpic/currents.h @@ -0,0 +1,117 @@ +#ifndef ENGINES_SRPIC_CURRENTS_H +#define ENGINES_SRPIC_CURRENTS_H + +#include "enums.h" +#include "global.h" + +#include "arch/kokkos_aliases.h" +#include "utils/log.h" +#include "utils/param_container.h" + +#include "metrics/traits.h" + +#include "engines/srpic/utils.h" +#include "framework/domain/domain.h" +#include "framework/domain/metadomain.h" + +#include "kernels/currents_deposit.hpp" +#include "kernels/digital_filter.hpp" + +namespace ntt { + namespace srpic { + + template + requires metric::traits::HasD && metric::traits::HasCoordType + void deposit_with(const Particles& species, + const M& local_metric, + const scatter_ndfield_t& scatter_cur, + real_t dt) { + 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, + local_metric, + (real_t)(species.charge()), + dt)); + } + + template + void CurrentsDeposit(Domain& domain, + const prm::Parameters& engine_params) { + const auto dt = engine_params.get("dt"); + Kokkos::deep_copy(domain.fields.cur, ZERO); + auto scatter_cur = Kokkos::Experimental::create_scatter_view( + domain.fields.cur); + for (auto& species : domain.species) { + if ((species.pusher() == ParticlePusher::NONE) or + (species.npart() == 0) or cmp::AlmostZero_host(species.charge())) { + continue; + } + logger::Checkpoint( + fmt::format("Launching currents deposit kernel for %d [%s] : %lu %f", + species.index(), + species.label().c_str(), + species.npart(), + (double)species.charge()), + HERE); + + deposit_with(species, domain.mesh.metric, scatter_cur, dt); + } + Kokkos::Experimental::contribute(domain.fields.cur, scatter_cur); + } + + template + requires metric::traits::HasD && metric::traits::HasCoordType + void CurrentsFilter(Metadomain& metadomain, + Domain& domain, + const SimulationParams& params) { + logger::Checkpoint("Launching currents filtering kernels", HERE); + auto range = srpic::RangeWithAxisBCs(domain); + const auto nfilter = params.template get( + "algorithms.current_filters"); + tuple_t size; + if constexpr (M::Dim == Dim::_1D || M::Dim == Dim::_2D || M::Dim == Dim::_3D) { + size[0] = domain.mesh.n_active(in::x1); + } + if constexpr (M::Dim == Dim::_2D || M::Dim == Dim::_3D) { + size[1] = domain.mesh.n_active(in::x2); + } + if constexpr (M::Dim == Dim::_3D) { + size[2] = domain.mesh.n_active(in::x3); + } + // !TODO: this needs to be done more efficiently + for (auto i { 0u }; i < nfilter; ++i) { + Kokkos::deep_copy(domain.fields.buff, domain.fields.cur); + Kokkos::parallel_for("CurrentsFilter", + range, + kernel::DigitalFilter_kernel( + domain.fields.cur, + domain.fields.buff, + size, + domain.mesh.flds_bc())); + metadomain.CommunicateFields(domain, Comm::J); + } + } + + } // namespace srpic +} // namespace ntt + +#endif // ENGINES_SRPIC_CURRENTS_H diff --git a/src/engines/srpic/fields_bcs.h b/src/engines/srpic/fields_bcs.h new file mode 100644 index 00000000..a0b159cf --- /dev/null +++ b/src/engines/srpic/fields_bcs.h @@ -0,0 +1,664 @@ +#ifndef ENGINES_SRPIC_FIELDS_BCS_H +#define ENGINES_SRPIC_FIELDS_BCS_H + +#include "enums.h" +#include "global.h" + +#include "arch/directions.h" +#include "utils/numeric.h" + +#include "metrics/traits.h" + +#include "archetypes/traits.h" +#include "engines/srpic/utils.h" +#include "framework/domain/domain.h" +#include "framework/parameters/parameters.h" + +#include "kernels/fields_bcs.hpp" + +namespace ntt { + namespace srpic { + + template + requires metric::traits::HasD + void CallMatchFields(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), + kernel::bc::MatchBoundaries_kernel(fields, + match_fields, + metric, + xg_edge, + ds, + tags, + boundaries)); + } + + template + requires metric::traits::HasD + void MatchFieldsIn(dir::direction_t direction, + Domain& domain, + const Grid& global_grid, + const PG& pgen, + const prm::Parameters& engine_params, + const SimulationParams& params, + BCTags tags) { + const auto time = engine_params.get("time"); + /** + * matching boundaries + */ + const auto ds_array = params.template get>( + "grid.boundaries.match.ds"); + const auto dim = direction.get_dim(); + real_t xg_min, xg_max, xg_edge; + auto sign = direction.get_sign(); + real_t ds; + if (sign > 0) { // + direction + ds = ds_array[(short)dim].second; + xg_max = global_grid.extent(dim).second; + xg_min = xg_max - ds; + xg_edge = xg_max; + } else { // - direction + ds = ds_array[(short)dim].first; + xg_min = global_grid.extent(dim).first; + xg_max = xg_min + ds; + xg_edge = xg_min; + } + boundaries_t box; + boundaries_t incl_ghosts; + for (dim_t d { 0 }; d < M::Dim; ++d) { + if (d == static_cast(dim)) { + box.push_back({ xg_min, xg_max }); + if (sign > 0) { + incl_ghosts.push_back({ false, true }); + } else { + incl_ghosts.push_back({ true, false }); + } + } 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 (auto d { 0u }; d < M::Dim; ++d) { + range_min[d] = intersect_range[d].first; + range_max[d] = intersect_range[d].second; + } + + if (dim == in::x1) { + if constexpr (arch::traits::pgen::HasMatchFields) { + auto match_fields = pgen.MatchFields(time); + CallMatchFields(domain.fields.em, + domain.mesh.flds_bc(), + match_fields, + domain.mesh.metric, + xg_edge, + ds, + tags, + range_min, + range_max); + } else if constexpr (arch::traits::pgen::HasMatchFieldsInX1) { + auto match_fields = pgen.MatchFieldsInX1(time); + CallMatchFields(domain.fields.em, + domain.mesh.flds_bc(), + match_fields, + domain.mesh.metric, + xg_edge, + ds, + tags, + range_min, + range_max); + } + } else if (dim == in::x2) { + if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { + if constexpr (arch::traits::pgen::HasMatchFields) { + auto match_fields = pgen.MatchFields(time); + CallMatchFields( + domain.fields.em, + domain.mesh.flds_bc(), + match_fields, + domain.mesh.metric, + xg_edge, + ds, + tags, + range_min, + range_max); + } else if constexpr (arch::traits::pgen::HasMatchFieldsInX2) { + auto match_fields = pgen.MatchFieldsInX2(time); + CallMatchFields( + domain.fields.em, + domain.mesh.flds_bc(), + match_fields, + domain.mesh.metric, + xg_edge, + ds, + tags, + range_min, + range_max); + } + } else { + raise::Error("Invalid dimension", HERE); + } + } else if (dim == in::x3) { + if constexpr (M::Dim == Dim::_3D) { + if constexpr (arch::traits::pgen::HasMatchFields) { + auto match_fields = pgen.MatchFields(time); + CallMatchFields( + domain.fields.em, + domain.mesh.flds_bc(), + match_fields, + domain.mesh.metric, + xg_edge, + ds, + tags, + range_min, + range_max); + } else if constexpr (arch::traits::pgen::HasMatchFieldsInX3) { + auto match_fields = pgen.MatchFieldsInX3(time); + CallMatchFields( + domain.fields.em, + domain.mesh.flds_bc(), + match_fields, + domain.mesh.metric, + xg_edge, + ds, + tags, + range_min, + range_max); + } + } + } else { + raise::Error("Invalid dimension", HERE); + } + } + + template + requires metric::traits::HasD + void AxisFieldsIn(dir::direction_t direction, + Domain& domain, + BCTags tags) { + /** + * axis boundaries + */ + if constexpr (M::CoordType != Coord::Cart) { + raise::ErrorIf(direction.get_dim() != in::x2, + "Invalid axis direction, should be x2", + HERE); + const auto i2_min = domain.mesh.i_min(in::x2); + const auto i2_max = domain.mesh.i_max(in::x2); + if (direction.get_sign() < 0) { + Kokkos::parallel_for( + "AxisBCFields", + domain.mesh.n_all(in::x1), + kernel::bc::AxisBoundaries_kernel(domain.fields.em, + i2_min, + tags)); + } else { + Kokkos::parallel_for( + "AxisBCFields", + domain.mesh.n_all(in::x1), + kernel::bc::AxisBoundaries_kernel(domain.fields.em, + i2_max, + tags)); + } + } else { + (void)direction; + (void)domain; + (void)tags; + raise::Error("Invalid coordinate type for axis BCs", HERE); + } + } + + template + requires metric::traits::HasD + void FixedFieldsIn(dir::direction_t direction, + Domain& domain, + const PG& pgen, + BCTags tags) { + /** + * fixed field boundaries + */ + const auto sign = direction.get_sign(); + const auto dim = direction.get_dim(); + raise::ErrorIf(dim != in::x1 and M::CoordType != Coord::Cart, + "Fixed BCs only implemented for x1 in " + "non-cartesian coordinates", + HERE); + em normal_b_comp, tang_e_comp1, tang_e_comp2; + if (dim == in::x1) { + normal_b_comp = em::bx1; + tang_e_comp1 = em::ex2; + tang_e_comp2 = em::ex3; + } else if (dim == in::x2) { + normal_b_comp = em::bx2; + tang_e_comp1 = em::ex1; + tang_e_comp2 = em::ex3; + } else if (dim == in::x3) { + normal_b_comp = em::bx3; + tang_e_comp1 = em::ex1; + tang_e_comp2 = em::ex2; + } else { + raise::Error("Invalid dimension", HERE); + } + std::vector xi_min, xi_max; + const std::vector all_dirs { in::x1, in::x2, in::x3 }; + for (dim_t d { 0u }; d < M::Dim; ++d) { + const auto dd = all_dirs[d]; + if (dim == dd) { + if (sign > 0) { // + direction + xi_min.push_back(domain.mesh.n_all(dd) - N_GHOSTS); + xi_max.push_back(domain.mesh.n_all(dd)); + } else { // - direction + xi_min.push_back(0); + xi_max.push_back(N_GHOSTS); + } + } else { + xi_min.push_back(0); + xi_max.push_back(domain.mesh.n_all(dd)); + } + } + raise::ErrorIf(xi_min.size() != xi_max.size() or + xi_min.size() != static_cast(M::Dim), + "Invalid range size", + HERE); + std::vector comps; + if (tags & BC::E) { + comps.push_back(tang_e_comp1); + comps.push_back(tang_e_comp2); + } + if (tags & BC::B) { + comps.push_back(normal_b_comp); + } + if constexpr (arch::traits::pgen::HasFixFieldsConst) { + for (const auto& comp : comps) { + auto value = ZERO; + bool shouldset = false; + // if fix field function present, read from it + const auto newset = pgen.FixFieldsConst((bc_in)(sign * ((short)dim + 1)), + (em)comp); + value = newset.first; + shouldset = newset.second; + if (shouldset) { + if constexpr (M::Dim == Dim::_1D) { + Kokkos::deep_copy( + Kokkos::subview(domain.fields.em, + std::make_pair(xi_min[0], xi_max[0]), + comp), + value); + } else if constexpr (M::Dim == Dim::_2D) { + Kokkos::deep_copy( + Kokkos::subview(domain.fields.em, + std::make_pair(xi_min[0], xi_max[0]), + std::make_pair(xi_min[1], xi_max[1]), + comp), + value); + } else if constexpr (M::Dim == Dim::_3D) { + Kokkos::deep_copy( + Kokkos::subview(domain.fields.em, + std::make_pair(xi_min[0], xi_max[0]), + std::make_pair(xi_min[1], xi_max[1]), + std::make_pair(xi_min[2], xi_max[2]), + comp), + value); + } else { + raise::Error("Invalid dimension", HERE); + } + } + } + } else { + (void)direction; + (void)domain; + (void)tags; + raise::Error("Fixed fields not present (both const and non-const)", HERE); + } + } + + template + requires metric::traits::HasD + void PerfectConductorFieldsIn(dir::direction_t direction, + Domain& domain, + BCTags tags) { + /** + * perfect conductor field boundaries + */ + if constexpr (M::CoordType != Coord::Cart) { + (void)direction; + (void)domain; + (void)tags; + raise::Error( + "Perfect conductor BCs only applicable to cartesian coordinates", + HERE); + } else { + const auto sign = direction.get_sign(); + const auto dim = direction.get_dim(); + + std::vector xi_min, xi_max; + + const std::vector all_dirs { in::x1, in::x2, in::x3 }; + + for (auto d { 0u }; d < M::Dim; ++d) { + const auto dd = all_dirs[d]; + if (dim == dd) { + xi_min.push_back(0); + xi_max.push_back((sign < 0) ? (N_GHOSTS + 1) : N_GHOSTS); + } else { + xi_min.push_back(0); + xi_max.push_back(domain.mesh.n_all(dd)); + } + } + raise::ErrorIf(xi_min.size() != xi_max.size() or + xi_min.size() != static_cast(M::Dim), + "Invalid range size", + HERE); + + range_t range; + if constexpr (M::Dim == Dim::_1D) { + range = CreateRangePolicy({ xi_min[0] }, { xi_max[0] }); + } else if constexpr (M::Dim == Dim::_2D) { + range = CreateRangePolicy({ xi_min[0], xi_min[1] }, + { xi_max[0], xi_max[1] }); + } else if constexpr (M::Dim == Dim::_3D) { + range = CreateRangePolicy({ xi_min[0], xi_min[1], xi_min[2] }, + { xi_max[0], xi_max[1], xi_max[2] }); + } else { + raise::Error("Invalid dimension", HERE); + } + std::size_t i_edge; + if (sign > 0) { + i_edge = domain.mesh.i_max(dim); + } else { + i_edge = domain.mesh.i_min(dim); + } + + if (dim == in::x1) { + if (sign > 0) { + Kokkos::parallel_for( + "ConductorFields", + range, + kernel::bc::ConductorBoundaries_kernel( + domain.fields.em, + i_edge, + tags)); + } else { + Kokkos::parallel_for( + "ConductorFields", + range, + kernel::bc::ConductorBoundaries_kernel( + domain.fields.em, + i_edge, + tags)); + } + } else if (dim == in::x2) { + if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { + if (sign > 0) { + Kokkos::parallel_for( + "ConductorFields", + range, + kernel::bc::ConductorBoundaries_kernel( + domain.fields.em, + i_edge, + tags)); + } else { + Kokkos::parallel_for( + "ConductorFields", + range, + kernel::bc::ConductorBoundaries_kernel( + domain.fields.em, + i_edge, + tags)); + } + } else { + raise::Error("Invalid dimension", HERE); + } + } else { + if constexpr (M::Dim == Dim::_3D) { + if (sign > 0) { + Kokkos::parallel_for( + "ConductorFields", + range, + kernel::bc::ConductorBoundaries_kernel( + domain.fields.em, + i_edge, + tags)); + } else { + Kokkos::parallel_for( + "ConductorFields", + range, + kernel::bc::ConductorBoundaries_kernel( + domain.fields.em, + i_edge, + tags)); + } + } else { + raise::Error("Invalid dimension", HERE); + } + } + } + } + + template + void AtmosphereFieldsIn(dir::direction_t direction, + Domain& domain, + const M& global_metric, + const Grid& global_grid, + const PG& pgen, + const SimulationParams& params, + const prm::Parameters& engine_params, + BCTags tags) { + const auto time = engine_params.get("time"); + /** + * atmosphere field boundaries + */ + if constexpr (arch::traits::pgen::HasAtmFields) { + const auto [sign, dim, xg_min, xg_max] = GetAtmosphereExtent(direction, + global_metric, + global_grid, + params); + // get_atm_extent(direction); + const auto dd = static_cast(dim); + boundaries_t box; + boundaries_t incl_ghosts; + for (auto d { 0u }; d < M::Dim; ++d) { + if (d == dd) { + box.push_back({ xg_min, xg_max }); + if (sign > 0) { + incl_ghosts.push_back({ false, true }); + } else { + incl_ghosts.push_back({ true, false }); + } + } 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 (auto d { 0u }; d < M::Dim; ++d) { + range_min[d] = intersect_range[d].first; + range_max[d] = intersect_range[d].second; + } + auto atm_fields = pgen.AtmFields(time); + std::size_t il_edge; + if (sign > 0) { + il_edge = range_min[dd] - N_GHOSTS; + } else { + il_edge = range_max[dd] - 1 - N_GHOSTS; + } + const auto range = CreateRangePolicy(range_min, range_max); + if (dim == in::x1) { + if (sign > 0) { + Kokkos::parallel_for( + "AtmosphereBCFields", + range, + kernel::bc::EnforcedBoundaries_kernel( + domain.fields.em, + atm_fields, + domain.mesh.metric, + il_edge, + tags)); + } else { + Kokkos::parallel_for( + "AtmosphereBCFields", + range, + kernel::bc::EnforcedBoundaries_kernel( + domain.fields.em, + atm_fields, + domain.mesh.metric, + il_edge, + tags)); + } + } else if (dim == in::x2) { + if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { + if (sign > 0) { + Kokkos::parallel_for( + "AtmosphereBCFields", + range, + kernel::bc::EnforcedBoundaries_kernel( + domain.fields.em, + atm_fields, + domain.mesh.metric, + il_edge, + tags)); + } else { + Kokkos::parallel_for( + "AtmosphereBCFields", + range, + kernel::bc::EnforcedBoundaries_kernel( + domain.fields.em, + atm_fields, + domain.mesh.metric, + il_edge, + tags)); + } + } else { + raise::Error("Invalid dimension", HERE); + } + } else if (dim == in::x3) { + if constexpr (M::Dim == Dim::_3D) { + if (sign > 0) { + Kokkos::parallel_for( + "AtmosphereBCFields", + range, + kernel::bc::EnforcedBoundaries_kernel( + domain.fields.em, + atm_fields, + domain.mesh.metric, + il_edge, + tags)); + } else { + Kokkos::parallel_for( + "AtmosphereBCFields", + range, + kernel::bc::EnforcedBoundaries_kernel( + domain.fields.em, + atm_fields, + domain.mesh.metric, + il_edge, + tags)); + } + } else { + raise::Error("Invalid dimension", HERE); + } + } else { + raise::Error("Invalid dimension", HERE); + } + } else { + (void)direction; + (void)domain; + (void)tags; + raise::Error("Atm fields not implemented in PGEN for atmosphere BCs", HERE); + } + } + + template + requires metric::traits::HasD + void CustomFieldsIn(dir::direction_t direction, + Domain& domain, + BCTags tags) { + (void)direction; + (void)domain; + (void)tags; + raise::Error("Custom boundaries not implemented", HERE); + // if constexpr ( + // traits::has_member::value) { + // const auto [box, custom_fields] = m_pgen.CustomFields(time); + // if (domain.mesh.Intersects(box)) { + // } + // + // } else { + // raise::Error("Custom boundaries not implemented", HERE); + // } + } + + template + requires metric::traits::HasD + void FieldBoundaries(Domain& domain, + const M& global_metric, + const Grid& global_grid, + const PG& pgen, + const prm::Parameters& engine_params, + const SimulationParams& params, + BCTags tags) { + for (auto& direction : dir::Directions::orth) { + if (global_grid.flds_bc_in(direction) == FldsBC::MATCH) { + MatchFieldsIn(direction, + domain, + global_grid, + pgen, + engine_params, + params, + tags); + } else if (global_grid.flds_bc_in(direction) == FldsBC::AXIS) { + if (domain.mesh.flds_bc_in(direction) == FldsBC::AXIS) { + AxisFieldsIn(direction, domain, tags); + } + } else if (global_grid.flds_bc_in(direction) == FldsBC::ATMOSPHERE) { + AtmosphereFieldsIn(direction, + domain, + global_metric, + global_grid, + pgen, + params, + engine_params, + tags); + } else if (global_grid.flds_bc_in(direction) == FldsBC::FIXED) { + if (domain.mesh.flds_bc_in(direction) == FldsBC::FIXED) { + FixedFieldsIn(direction, domain, pgen, tags); + } + } else if (global_grid.flds_bc_in(direction) == FldsBC::CONDUCTOR) { + if (domain.mesh.flds_bc_in(direction) == FldsBC::CONDUCTOR) { + PerfectConductorFieldsIn(direction, domain, tags); + } + } else if (global_grid.flds_bc_in(direction) == FldsBC::CUSTOM) { + if (domain.mesh.flds_bc_in(direction) == FldsBC::CUSTOM) { + CustomFieldsIn(direction, domain, tags); + } + } else if (global_grid.flds_bc_in(direction) == FldsBC::HORIZON) { + raise::Error("HORIZON BCs only applicable for GR", HERE); + } + } // loop over directions + } + + } // namespace srpic +} // namespace ntt + +#endif // ENGINES_SRPIC_FIELDS_BCS_H diff --git a/src/engines/srpic/fieldsolvers.h b/src/engines/srpic/fieldsolvers.h new file mode 100644 index 00000000..04bca274 --- /dev/null +++ b/src/engines/srpic/fieldsolvers.h @@ -0,0 +1,193 @@ +#ifndef ENGINES_SRPIC_FIELDSOLVERS_H +#define ENGINES_SRPIC_FIELDSOLVERS_H + +#include "enums.h" +#include "global.h" + +#include "utils/log.h" +#include "utils/numeric.h" +#include "utils/param_container.h" + +#include "archetypes/traits.h" +#include "engines/srpic/utils.h" +#include "framework/domain/domain.h" +#include "framework/parameters/parameters.h" + +#include "kernels/ampere_mink.hpp" +#include "kernels/ampere_sr.hpp" +#include "kernels/faraday_mink.hpp" +#include "kernels/faraday_sr.hpp" + +namespace ntt { + namespace srpic { + + template + void Faraday(Domain& domain, + const prm::Parameters& engine_params, + const SimulationParams& params, + real_t fraction = ONE) { + logger::Checkpoint("Launching Faraday kernel", HERE); + const auto dt = engine_params.get("dt"); + + const auto dT = fraction * + params.template get( + "algorithms.timestep.correction") * + dt; + if constexpr (M::CoordType == Coord::Cart) { + // minkowski case + const auto dx = math::sqrt(domain.mesh.metric.template h_<1, 1>({})); + const auto deltax = params.template get( + "algorithms.fieldsolver.delta_x"); + const auto deltay = params.template get( + "algorithms.fieldsolver.delta_y"); + const auto betaxy = params.template get( + "algorithms.fieldsolver.beta_xy"); + const auto betayx = params.template get( + "algorithms.fieldsolver.beta_yx"); + const auto deltaz = params.template get( + "algorithms.fieldsolver.delta_z"); + const auto betaxz = params.template get( + "algorithms.fieldsolver.beta_xz"); + const auto betazx = params.template get( + "algorithms.fieldsolver.beta_zx"); + const auto betayz = params.template get( + "algorithms.fieldsolver.beta_yz"); + const auto betazy = params.template get( + "algorithms.fieldsolver.beta_zy"); + real_t coeff1, coeff2; + if constexpr (M::Dim == Dim::_2D) { + coeff1 = dT / SQR(dx); + coeff2 = dT; + } else { + coeff1 = dT / dx; + coeff2 = ZERO; + } + Kokkos::parallel_for("Faraday", + domain.mesh.rangeActiveCells(), + kernel::mink::Faraday_kernel(domain.fields.em, + coeff1, + coeff2, + deltax, + deltay, + betaxy, + betayx, + deltaz, + betaxz, + betazx, + betayz, + betazy)); + } else { + Kokkos::parallel_for("Faraday", + domain.mesh.rangeActiveCells(), + kernel::sr::Faraday_kernel(domain.fields.em, + domain.mesh.metric, + dT, + domain.mesh.flds_bc())); + } + } + + template + void Ampere(Domain& domain, + const prm::Parameters& engine_params, + const SimulationParams& params, + real_t fraction = ONE) { + logger::Checkpoint("Launching Ampere kernel", HERE); + const auto dt = engine_params.get("dt"); + + const auto dT = fraction * + params.template get( + "algorithms.timestep.correction") * + dt; + auto range = RangeWithAxisBCs(domain); + if constexpr (M::CoordType == Coord::Cart) { + // minkowski case + const auto dx = math::sqrt(domain.mesh.metric.template h_<1, 1>({})); + real_t coeff1, coeff2; + if constexpr (M::Dim == Dim::_2D) { + coeff1 = dT / SQR(dx); + coeff2 = dT; + } else { + coeff1 = dT / dx; + coeff2 = ZERO; + } + + Kokkos::parallel_for( + "Ampere", + range, + kernel::mink::Ampere_kernel(domain.fields.em, coeff1, coeff2)); + } else { + const auto ni2 = domain.mesh.n_active(in::x2); + Kokkos::parallel_for("Ampere", + range, + kernel::sr::Ampere_kernel(domain.fields.em, + domain.mesh.metric, + dT, + ni2, + domain.mesh.flds_bc())); + } + } + + template + void CurrentsAmpere(Domain& domain, + const prm::Parameters& engine_params, + const SimulationParams& params, + const PG& pgen) { + logger::Checkpoint("Launching Ampere kernel for adding currents", HERE); + const auto dt = engine_params.get("dt"); + + const auto q0 = params.template get("scales.q0"); + const auto n0 = params.template get("scales.n0"); + const auto B0 = params.template get("scales.B0"); + if constexpr (M::CoordType == Coord::Cart) { + // minkowski case + const auto V0 = params.template get("scales.V0"); + const auto ppc0 = params.template get("particles.ppc0"); + const auto coeff = -dt * q0 / (B0 * V0); + if constexpr (arch::traits::pgen::HasExtCurrent) { + const std::vector xmin { domain.mesh.extent(in::x1).first, + domain.mesh.extent(in::x2).first, + domain.mesh.extent(in::x3).first }; + const auto ext_current = pgen.ext_current; + const auto dx = domain.mesh.metric.template sqrt_h_<1, 1>({}); + Kokkos::parallel_for( + "Ampere", + domain.mesh.rangeActiveCells(), + kernel::mink::CurrentsAmpere_kernel( + domain.fields.em, + domain.fields.cur, + coeff, + ppc0, + ext_current, + xmin, + dx)); + } else { + Kokkos::parallel_for( + "Ampere", + domain.mesh.rangeActiveCells(), + kernel::mink::CurrentsAmpere_kernel(domain.fields.em, + domain.fields.cur, + coeff, + ppc0)); + } + } else { + // non-minkowski + const auto coeff = -dt * q0 * n0 / B0; + auto range = RangeWithAxisBCs(domain); + const auto ni2 = domain.mesh.n_active(in::x2); + Kokkos::parallel_for( + "Ampere", + range, + kernel::sr::CurrentsAmpere_kernel(domain.fields.em, + domain.fields.cur, + domain.mesh.metric, + coeff, + ONE / n0, + ni2, + domain.mesh.flds_bc())); + } + } + + } // namespace srpic +} // namespace ntt + +#endif // ENGINES_SRPIC_FIELDSOLVERS_H diff --git a/src/engines/srpic/particle_pusher.h b/src/engines/srpic/particle_pusher.h new file mode 100644 index 00000000..b263d6ab --- /dev/null +++ b/src/engines/srpic/particle_pusher.h @@ -0,0 +1,219 @@ +#ifndef ENGINES_SRPIC_PARTICLE_PUSHER_H +#define ENGINES_SRPIC_PARTICLE_PUSHER_H + +#include "enums.h" +#include "global.h" + +#include "utils/log.h" +#include "utils/numeric.h" +#include "utils/param_container.h" + +#include "metrics/traits.h" + +#include "archetypes/traits.h" +#include "engines/srpic/utils.h" +#include "framework/domain/domain.h" +#include "framework/domain/grid.h" +#include "framework/parameters/parameters.h" + +#include "kernels/particle_pusher_sr.hpp" + +namespace ntt { + namespace srpic { + + template + requires metric::traits::HasD + void ParticlePush(Domain& domain, + const Grid& global_grid, + const M& global_metric, + const prm::Parameters& engine_params, + const SimulationParams& params, + const PG& pgen) { + const auto dt = engine_params.get("dt"); + const auto time = engine_params.get("time"); + + real_t gx1 { ZERO }, gx2 { ZERO }, gx3 { ZERO }, ds { ZERO }; + real_t x_surf { ZERO }; + bool has_atmosphere = false; + for (auto& direction : dir::Directions::orth) { + if (global_grid.prtl_bc_in(direction) == PrtlBC::ATMOSPHERE) { + raise::ErrorIf(has_atmosphere, + "Only one direction is allowed to have atm boundaries", + HERE); + has_atmosphere = true; + const auto g = params.template get( + "grid.boundaries.atmosphere.g"); + ds = params.template get("grid.boundaries.atmosphere.ds"); + const auto [sign, dim, xg_min, xg_max] = + GetAtmosphereExtent(direction, global_metric, global_grid, params); + if (dim == in::x1) { + gx1 = sign > 0 ? g : -g; + gx2 = ZERO; + gx3 = ZERO; + } else if (dim == in::x2) { + gx1 = ZERO; + gx2 = sign > 0 ? g : -g; + gx3 = ZERO; + } else if (dim == in::x3) { + gx1 = ZERO; + gx2 = ZERO; + gx3 = sign > 0 ? g : -g; + } else { + raise::Error("Invalid dimension", HERE); + } + if (sign > 0) { + x_surf = xg_min; + } else { + x_surf = xg_max; + } + } + } + for (auto& species : domain.species) { + if ((species.pusher() == ParticlePusher::NONE) or (species.npart() == 0)) { + continue; + } + species.set_unsorted(); + logger::Checkpoint( + fmt::format("Launching particle pusher kernel for %d [%s] : %lu", + species.index(), + species.label().c_str(), + species.npart()), + HERE); + + kernel::sr::PusherParams pusher_params {}; + pusher_params.pusher_flags = species.pusher(); + pusher_params.radiative_drag_flags = species.radiative_drag_flags(); + pusher_params.mass = species.mass(); + pusher_params.charge = species.charge(); + pusher_params.time = time; + pusher_params.dt = dt; + pusher_params.omegaB0 = params.template get("scales.omegaB0"); + pusher_params.ni1 = domain.mesh.n_active(in::x1); + pusher_params.ni2 = domain.mesh.n_active(in::x2); + pusher_params.ni3 = domain.mesh.n_active(in::x3); + pusher_params.boundaries = domain.mesh.prtl_bc(); + + if (species.pusher() & ParticlePusher::GCA) { + pusher_params.gca_params.set( + "larmor_max", + params.template get("algorithms.gca.larmor_max")); + pusher_params.gca_params.set( + "e_ovr_b_max", + params.template get("algorithms.gca.e_ovr_b_max")); + } + + if (species.radiative_drag_flags() & RadiativeDrag::SYNCHROTRON) { + pusher_params.radiative_drag_params.set( + "synchrotron_gamma_rad", + params.template get( + "radiation.drag.synchrotron.gamma_rad")); + } + + if (species.radiative_drag_flags() & RadiativeDrag::COMPTON) { + pusher_params.radiative_drag_params.set( + "compton_gamma_rad", + params.template get("radiation.drag.compton.gamma_rad")); + } + + kernel::sr::PusherArrays pusher_arrays {}; + pusher_arrays.sp = species.index(); + pusher_arrays.i1 = species.i1; + pusher_arrays.i2 = species.i2; + pusher_arrays.i3 = species.i3; + pusher_arrays.i1_prev = species.i1_prev; + pusher_arrays.i2_prev = species.i2_prev; + pusher_arrays.i3_prev = species.i3_prev; + pusher_arrays.dx1 = species.dx1; + pusher_arrays.dx2 = species.dx2; + pusher_arrays.dx3 = species.dx3; + pusher_arrays.dx1_prev = species.dx1_prev; + pusher_arrays.dx2_prev = species.dx2_prev; + pusher_arrays.dx3_prev = species.dx3_prev; + pusher_arrays.ux1 = species.ux1; + pusher_arrays.ux2 = species.ux2; + pusher_arrays.ux3 = species.ux3; + pusher_arrays.phi = species.phi; + pusher_arrays.tag = species.tag; + + // toggle to indicate whether pgen defines the external force + bool has_extforce = false; + if constexpr (arch::traits::pgen::HasExtForce) { + has_extforce = true; + // toggle to indicate whether the ext force applies to current species + if (::traits::has_member<::traits::species_t, decltype(PG::ext_force)>::value) { + has_extforce &= std::find(pgen.ext_force.species.begin(), + pgen.ext_force.species.end(), + species.index()) != + pgen.ext_force.species.end(); + } + } + + pusher_params.ext_force = has_extforce; + + if (not has_atmosphere and not has_extforce) { + Kokkos::parallel_for("ParticlePusher", + species.rangeActiveParticles(), + kernel::sr::Pusher_kernel(pusher_params, + pusher_arrays, + domain.fields.em, + domain.mesh.metric)); + } else if (has_atmosphere and not has_extforce) { + const auto force = + kernel::sr::Force { + { gx1, gx2, gx3 }, + x_surf, + ds + }; + Kokkos::parallel_for( + "ParticlePusher", + species.rangeActiveParticles(), + kernel::sr::Pusher_kernel(pusher_params, + pusher_arrays, + domain.fields.em, + domain.mesh.metric, + force)); + } else if (not has_atmosphere and has_extforce) { + if constexpr (arch::traits::pgen::HasExtForce) { + const auto force = + kernel::sr::Force { + pgen.ext_force + }; + Kokkos::parallel_for( + "ParticlePusher", + species.rangeActiveParticles(), + kernel::sr::Pusher_kernel(pusher_params, + pusher_arrays, + domain.fields.em, + domain.mesh.metric, + force)); + } else { + raise::Error("External force not implemented", HERE); + } + } else { // has_atmosphere and has_extforce + if constexpr (arch::traits::pgen::HasExtForce) { + const auto force = + kernel::sr::Force { + pgen.ext_force, + { gx1, gx2, gx3 }, + x_surf, + ds + }; + Kokkos::parallel_for( + "ParticlePusher", + species.rangeActiveParticles(), + kernel::sr::Pusher_kernel(pusher_params, + pusher_arrays, + domain.fields.em, + domain.mesh.metric, + force)); + } else { + raise::Error("External force not implemented", HERE); + } + } + } + } + + } // namespace srpic +} // namespace ntt + +#endif // ENGINES_SRPIC_PARTICLE_PUSHER_H diff --git a/src/engines/srpic/particles_bcs.h b/src/engines/srpic/particles_bcs.h new file mode 100644 index 00000000..72a58cd9 --- /dev/null +++ b/src/engines/srpic/particles_bcs.h @@ -0,0 +1,307 @@ +#ifndef ENGINES_SRPIC_PARTICLES_BCS_H +#define ENGINES_SRPIC_PARTICLES_BCS_H + +#include "enums.h" +#include "global.h" + +#include "arch/directions.h" +#include "utils/numeric.h" + +#include "metrics/traits.h" + +#include "archetypes/energy_dist.h" +#include "archetypes/particle_injector.h" +#include "archetypes/spatial_dist.h" +#include "engines/srpic/utils.h" +#include "framework/domain/domain.h" +#include "framework/domain/metadomain.h" +#include "framework/parameters/parameters.h" + +#include "kernels/particle_moments.hpp" + +namespace ntt { + namespace srpic { + + template + requires metric::traits::HasD && metric::traits::HasCoordType + void AtmosphereParticlesIn(dir::direction_t direction, + Metadomain& metadomain, + Domain& domain, + const SimulationParams& params, + InjTags tags) { + const auto [sign, dim, xg_min, xg_max] = srpic::GetAtmosphereExtent( + direction, + metadomain.mesh().metric, + metadomain.mesh(), + params); + + const auto x_surf = sign > 0 ? xg_min : xg_max; + const auto ds = params.template get( + "grid.boundaries.atmosphere.ds"); + const auto temp = params.template get( + "grid.boundaries.atmosphere.temperature"); + const auto height = params.template get( + "grid.boundaries.atmosphere.height"); + const auto species = params.template get>( + "grid.boundaries.atmosphere.species"); + const auto nmax = params.template get( + "grid.boundaries.atmosphere.density"); + + Kokkos::deep_copy(domain.fields.bckp, ZERO); + auto scatter_bckp = Kokkos::Experimental::create_scatter_view( + domain.fields.bckp); + const auto use_weights = M::CoordType != Coord::Cart; + const auto ni2 = domain.mesh.n_active(in::x2); + const auto inv_n0 = ONE / params.template get("scales.n0"); + + // compute the density of the two species + if (tags & Inj::AssumeEmpty) { + if constexpr (M::Dim == Dim::_1D) { + Kokkos::deep_copy( + Kokkos::subview(domain.fields.bckp, Kokkos::ALL, std::make_pair(0, 1)), + ZERO); + } else if constexpr (M::Dim == Dim::_2D) { + Kokkos::deep_copy(Kokkos::subview(domain.fields.bckp, + Kokkos::ALL, + Kokkos::ALL, + std::make_pair(0, 1)), + ZERO); + } else if constexpr (M::Dim == Dim::_3D) { + Kokkos::deep_copy(Kokkos::subview(domain.fields.bckp, + Kokkos::ALL, + Kokkos::ALL, + Kokkos::ALL, + std::make_pair(0, 1)), + ZERO); + } + } else { + for (const auto& sp : + std::vector { species.first, species.second }) { + auto& prtl_spec = domain.species[sp - 1]; + if (prtl_spec.npart() == 0) { + continue; + } + // clang-format off + Kokkos::parallel_for( + "ComputeMoments", + prtl_spec.rangeActiveParticles(), + kernel::ParticleMoments_kernel( + {}, scatter_bckp, 0, + prtl_spec.i1, prtl_spec.i2, prtl_spec.i3, + prtl_spec.dx1, prtl_spec.dx2, prtl_spec.dx3, + prtl_spec.ux1, prtl_spec.ux2, prtl_spec.ux3, + prtl_spec.phi, prtl_spec.weight, prtl_spec.tag, + prtl_spec.mass(), prtl_spec.charge(), + use_weights, + domain.mesh.metric, domain.mesh.flds_bc(), + ni2, inv_n0, 0)); + // clang-format on + prtl_spec.set_unsorted(); + } + Kokkos::Experimental::contribute(domain.fields.bckp, scatter_bckp); + metadomain.SynchronizeFields(domain, Comm::Bckp, { 0, 1 }); + } + + const auto maxwellian = arch::Maxwellian { + domain.mesh.metric, + domain.random_pool, + temp + }; + + if (dim == in::x1) { + if (sign > 0) { + auto target_density = + arch::AtmosphereDensityProfile { + nmax, + height, + x_surf, + ds + }; + const auto spatial_dist = + arch::Replenish { + domain.mesh.metric, + domain.fields.bckp, + 0, + target_density, + nmax + }; + arch::InjectNonUniform( + params, + domain, + { species.first, species.second }, + { maxwellian, maxwellian }, + spatial_dist, + nmax, + use_weights); + } else { + auto target_density = + arch::AtmosphereDensityProfile { + nmax, + height, + x_surf, + ds + }; + const auto spatial_dist = + arch::Replenish { + domain.mesh.metric, + domain.fields.bckp, + 0, + target_density, + nmax + }; + arch::InjectNonUniform( + params, + domain, + { species.first, species.second }, + { maxwellian, maxwellian }, + spatial_dist, + nmax, + use_weights); + } + } else if (dim == in::x2) { + if (sign > 0) { + auto target_density = + arch::AtmosphereDensityProfile { + nmax, + height, + x_surf, + ds + }; + const auto spatial_dist = + arch::Replenish { + domain.mesh.metric, + domain.fields.bckp, + 0, + target_density, + nmax + }; + arch::InjectNonUniform( + params, + domain, + { species.first, species.second }, + { maxwellian, maxwellian }, + spatial_dist, + nmax, + use_weights); + } else { + auto target_density = + arch::AtmosphereDensityProfile { + nmax, + height, + x_surf, + ds + }; + const auto spatial_dist = + arch::Replenish { + domain.mesh.metric, + domain.fields.bckp, + 0, + target_density, + nmax + }; + arch::InjectNonUniform( + params, + domain, + { species.first, species.second }, + { maxwellian, maxwellian }, + spatial_dist, + nmax, + use_weights); + } + } else if (dim == in::x3) { + if (sign > 0) { + auto target_density = + arch::AtmosphereDensityProfile { + nmax, + height, + x_surf, + ds + }; + const auto spatial_dist = + arch::Replenish { + domain.mesh.metric, + domain.fields.bckp, + 0, + target_density, + nmax + }; + arch::InjectNonUniform( + params, + domain, + { species.first, species.second }, + { maxwellian, maxwellian }, + spatial_dist, + nmax, + use_weights); + } else { + auto target_density = + arch::AtmosphereDensityProfile { + nmax, + height, + x_surf, + ds + }; + const auto spatial_dist = + arch::Replenish { + domain.mesh.metric, + domain.fields.bckp, + 0, + target_density, + nmax + }; + arch::InjectNonUniform( + params, + domain, + { species.first, species.second }, + { maxwellian, maxwellian }, + spatial_dist, + nmax, + use_weights); + } + } else { + raise::Error("Invalid dimension", HERE); + } + return; + } + + template + requires metric::traits::HasD + void ParticleInjector(Metadomain& metadomain, + Domain& domain, + const SimulationParams& params, + InjTags tags = Inj::None) { + for (auto& direction : dir::Directions::orth) { + if (metadomain.mesh().prtl_bc_in(direction) == PrtlBC::ATMOSPHERE) { + AtmosphereParticlesIn(direction, metadomain, domain, params, tags); + } + } + } + + } // namespace srpic +} // namespace ntt + +#endif diff --git a/src/engines/srpic/srpic.hpp b/src/engines/srpic/srpic.hpp new file mode 100644 index 00000000..a531ffe7 --- /dev/null +++ b/src/engines/srpic/srpic.hpp @@ -0,0 +1,199 @@ +/** + * @file engines/srpic.hpp + * @brief Simulation engien class which specialized on SRPIC + * @implements + * - ntt::SRPICEngine<> : ntt::Engine<> + * @cpp: + * - srpic.cpp + * @namespaces: + * - ntt:: + * @macros: + */ + +#ifndef ENGINES_SRPIC_SRPIC_H +#define ENGINES_SRPIC_SRPIC_H + +#include "enums.h" +#include "global.h" + +#include "utils/numeric.h" +#include "utils/timer.h" +#include "utils/toml.h" + +#include "engines/srpic/currents.h" +#include "engines/srpic/fields_bcs.h" +#include "engines/srpic/fieldsolvers.h" +#include "engines/srpic/particle_pusher.h" +#include "engines/srpic/particles_bcs.h" +#include "engines/traits.h" +#include "framework/domain/domain.h" +#include "framework/parameters/parameters.h" + +#include "engines/engine.hpp" +#include "pgen.hpp" + +#include +#include + +namespace ntt { + + template + requires traits::engine::IsCompatibleWithSRPICEngine + class SRPICEngine : public Engine { + + using base_t = Engine; + using pgen_t = user::PGen; + using domain_t = Domain; + // contents + using base_t::m_metadomain; + using base_t::m_params; + using base_t::m_pgen; + // methods + using base_t::init; + // variables + using base_t::dt; + using base_t::max_steps; + using base_t::runtime; + using base_t::step; + using base_t::time; + + public: + static constexpr auto S { SimEngine::SRPIC }; + + SRPICEngine(const SimulationParams& params) : base_t { params } {} + + ~SRPICEngine() = default; + + void step_forward(timer::Timers& timers, domain_t& dom) override { + const auto fieldsolver_enabled = m_params.template get( + "algorithms.fieldsolver.enable"); + const auto deposit_enabled = m_params.template get( + "algorithms.deposit.enable"); + const auto clear_interval = m_params.template get( + "particles.clear_interval"); + + if (step == 0) { + // communicate fields and apply BCs on the first timestep + m_metadomain.CommunicateFields(dom, Comm::B | Comm::E); + srpic::FieldBoundaries(dom, + m_metadomain.mesh().metric, + m_metadomain.mesh(), + m_pgen, + this->engineParams(), + m_params, + BC::B | BC::E); + srpic::ParticleInjector(m_metadomain, dom, m_params); + } + + if (fieldsolver_enabled) { + timers.start("FieldSolver"); + srpic::Faraday(dom, this->engineParams(), m_params, HALF); + timers.stop("FieldSolver"); + + timers.start("Communications"); + m_metadomain.CommunicateFields(dom, Comm::B); + timers.stop("Communications"); + + timers.start("FieldBoundaries"); + srpic::FieldBoundaries(dom, + m_metadomain.mesh().metric, + m_metadomain.mesh(), + m_pgen, + this->engineParams(), + m_params, + BC::B); + timers.stop("FieldBoundaries"); + Kokkos::fence(); + } + + { + timers.start("ParticlePusher"); + srpic::ParticlePush(dom, + m_metadomain.mesh(), + m_metadomain.mesh().metric, + this->engineParams(), + m_params, + m_pgen); + timers.stop("ParticlePusher"); + + if (deposit_enabled) { + timers.start("CurrentDeposit"); + srpic::CurrentsDeposit(dom, this->engineParams()); + timers.stop("CurrentDeposit"); + + timers.start("Communications"); + m_metadomain.SynchronizeFields(dom, Comm::J); + m_metadomain.CommunicateFields(dom, Comm::J); + timers.stop("Communications"); + + timers.start("CurrentFiltering"); + srpic::CurrentsFilter(m_metadomain, dom, m_params); + timers.stop("CurrentFiltering"); + } + + timers.start("Communications"); + m_metadomain.CommunicateParticles(dom); + timers.stop("Communications"); + } + + if (fieldsolver_enabled) { + timers.start("FieldSolver"); + srpic::Faraday(dom, this->engineParams(), m_params, HALF); + timers.stop("FieldSolver"); + + timers.start("Communications"); + m_metadomain.CommunicateFields(dom, Comm::B); + timers.stop("Communications"); + + timers.start("FieldBoundaries"); + srpic::FieldBoundaries(dom, + m_metadomain.mesh().metric, + m_metadomain.mesh(), + m_pgen, + this->engineParams(), + m_params, + BC::B); + timers.stop("FieldBoundaries"); + + timers.start("FieldSolver"); + srpic::Ampere(dom, this->engineParams(), m_params, ONE); + timers.stop("FieldSolver"); + + if (deposit_enabled) { + timers.start("FieldSolver"); + srpic::CurrentsAmpere(dom, this->engineParams(), m_params, m_pgen); + timers.stop("FieldSolver"); + } + + timers.start("Communications"); + m_metadomain.CommunicateFields(dom, Comm::E | Comm::J); + timers.stop("Communications"); + + timers.start("FieldBoundaries"); + srpic::FieldBoundaries(dom, + m_metadomain.mesh().metric, + m_metadomain.mesh(), + m_pgen, + this->engineParams(), + m_params, + BC::E); + timers.stop("FieldBoundaries"); + } + + { + timers.start("Injector"); + srpic::ParticleInjector(m_metadomain, dom, m_params); + timers.stop("Injector"); + } + + if (clear_interval > 0 and step % clear_interval == 0 and step > 0) { + timers.start("PrtlClear"); + m_metadomain.RemoveDeadParticles(dom); + timers.stop("PrtlClear"); + } + } + }; + +} // namespace ntt + +#endif // ENGINES_SRPIC_SRPIC_H diff --git a/src/engines/srpic/utils.h b/src/engines/srpic/utils.h new file mode 100644 index 00000000..9779e55e --- /dev/null +++ b/src/engines/srpic/utils.h @@ -0,0 +1,127 @@ +#ifndef ENGINES_SRPIC_UTILS_H +#define ENGINES_SRPIC_UTILS_H + +#include "enums.h" +#include "global.h" + +#include "arch/directions.h" +#include "utils/numeric.h" + +#include "metrics/traits.h" + +#include "framework/domain/domain.h" +#include "framework/domain/grid.h" +#include "framework/parameters/parameters.h" + +#include + +namespace ntt { + namespace srpic { + + /** + * @brief Get the buffer region of the atmosphere and the direction + * @param direction direction in which the atmosphere is applied + * @return tuple: [sign of the direction, the direction (as in::), the min and max extent + * @note xg_min and xg_max are the extents where the fields are set, not the atmosphere itself + * @note i.e. + * + * fields set particles injected + * ghost zone | | + * v v v + * |....|...........|*******************..... -> x1 + * ^ ^ + * xg_min xg_max + * | | | + * |<-- buffer -->|<-- atmosphere -->| + * + * in this case the function returns { -1, in::x1, xg_min, xg_max } + */ + template + requires metric::traits::HasD && metric::traits::HasConvert + auto GetAtmosphereExtent(dir::direction_t direction, + const M& global_metric, + const Grid& global_grid, + const SimulationParams& params) + -> std::tuple { + const auto sign = direction.get_sign(); + const auto dim = direction.get_dim(); + const auto min_buff = params.template get( + "algorithms.current_filters") + + 2; + const auto buffer_ncells = min_buff > 5 ? min_buff : 5; + if (M::CoordType != Coord::Cart and (dim != in::x1 or sign > 0)) { + raise::Error("For non-cartesian coordinates atmosphere BCs is " + "possible only in -x1 (@ rmin)", + HERE); + } + real_t xg_min { ZERO }, xg_max { ZERO }; + ncells_t ig_min, ig_max; + if (sign > 0) { // + direction + ig_min = global_grid.n_active(dim) - buffer_ncells; + ig_max = global_grid.n_active(dim); + } else { // - direction + ig_min = 0; + ig_max = buffer_ncells; + } + + if (dim == in::x1) { + xg_min = global_metric.template convert<1, Crd::Cd, Crd::Ph>( + static_cast(ig_min)); + xg_max = global_metric.template convert<1, Crd::Cd, Crd::Ph>( + static_cast(ig_max)); + } else if (dim == in::x2) { + if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { + xg_min = global_metric.template convert<2, Crd::Cd, Crd::Ph>( + static_cast(ig_min)); + xg_max = global_metric.template convert<2, Crd::Cd, Crd::Ph>( + static_cast(ig_max)); + } else { + raise::Error("Invalid dimension", HERE); + } + } else if (dim == in::x3) { + if constexpr (M::Dim == Dim::_3D) { + xg_min = global_metric.template convert<3, Crd::Cd, Crd::Ph>( + static_cast(ig_min)); + xg_max = global_metric.template convert<3, Crd::Cd, Crd::Ph>( + static_cast(ig_max)); + } else { + raise::Error("Invalid dimension", HERE); + } + } else { + raise::Error("Invalid dimension", HERE); + } + return { sign, dim, xg_min, xg_max }; + } + + template + auto RangeWithAxisBCs(const Domain& domain) + -> range_t { + auto range = domain.mesh.rangeActiveCells(); + if constexpr (M::CoordType != Coord::Cart) { + /** + * @brief taking one extra cell in the x2 direction if AXIS BCs + */ + if constexpr (M::Dim == Dim::_2D) { + if (domain.mesh.flds_bc_in({ 0, +1 }) == FldsBC::AXIS) { + range = CreateRangePolicy( + { domain.mesh.i_min(in::x1), domain.mesh.i_min(in::x2) }, + { domain.mesh.i_max(in::x1), domain.mesh.i_max(in::x2) + 1 }); + } + } else if constexpr (M::Dim == Dim::_3D) { + if (domain.mesh.flds_bc_in({ 0, +1, 0 }) == FldsBC::AXIS) { + range = CreateRangePolicy({ domain.mesh.i_min(in::x1), + domain.mesh.i_min(in::x2), + domain.mesh.i_min(in::x3) }, + { domain.mesh.i_max(in::x1), + domain.mesh.i_max(in::x2) + 1, + domain.mesh.i_max(in::x3) }); + } + } + } + return range; + } + + } // namespace srpic +} // namespace ntt + +#endif diff --git a/src/entity.cpp b/src/entity.cpp index 969a2637..9afb37cb 100644 --- a/src/entity.cpp +++ b/src/entity.cpp @@ -8,7 +8,7 @@ #include "framework/specialization_registry.h" #include "engines/grpic.hpp" -#include "engines/srpic.hpp" +#include "engines/srpic/srpic.hpp" #include "pgen.hpp" #include From 4e61680e64378b767252a6f349fd0714877c3927 Mon Sep 17 00:00:00 2001 From: haykh Date: Thu, 29 Jan 2026 21:25:42 -0500 Subject: [PATCH 3/9] minor renaming --- dev/nix/kokkos.nix | 8 ++-- src/engines/srpic/currents.h | 10 ++--- src/engines/srpic/particle_pusher.h | 64 +++++++++++++++++------------ 3 files changed, 46 insertions(+), 36 deletions(-) diff --git a/dev/nix/kokkos.nix b/dev/nix/kokkos.nix index fc51647b..ce2ba267 100644 --- a/dev/nix/kokkos.nix +++ b/dev/nix/kokkos.nix @@ -7,7 +7,7 @@ let name = "kokkos"; - pversion = "5.0.0"; + pversion = "5.0.1"; compilerPkgs = { "HIP" = with pkgs.rocmPackages; [ llvm.llvm @@ -17,16 +17,16 @@ let rocprim rocminfo rocm-smi - pkgs.llvmPackages_19.clang-tools + pkgs.clang-tools ]; "CUDA" = with pkgs.cudaPackages; [ - llvmPackages_19.clang-tools + pkgs.clang-tools cudatoolkit cuda_cudart pkgs.gcc13 ]; "NONE" = [ - pkgs.llvmPackages_19.clang-tools + pkgs.clang-tools pkgs.gcc13 ]; }; diff --git a/src/engines/srpic/currents.h b/src/engines/srpic/currents.h index 28b8f9ab..479ba0bc 100644 --- a/src/engines/srpic/currents.h +++ b/src/engines/srpic/currents.h @@ -22,10 +22,10 @@ namespace ntt { template requires metric::traits::HasD && metric::traits::HasCoordType - void deposit_with(const Particles& species, - const M& local_metric, - const scatter_ndfield_t& scatter_cur, - real_t dt) { + void CallDepositKernel(const Particles& species, + const M& local_metric, + const scatter_ndfield_t& scatter_cur, + real_t dt) { Kokkos::parallel_for("CurrentsDeposit", species.rangeActiveParticles(), kernel::DepositCurrents_kernel( @@ -73,7 +73,7 @@ namespace ntt { (double)species.charge()), HERE); - deposit_with(species, domain.mesh.metric, scatter_cur, dt); + CallDepositKernel(species, domain.mesh.metric, scatter_cur, dt); } Kokkos::Experimental::contribute(domain.fields.cur, scatter_cur); } diff --git a/src/engines/srpic/particle_pusher.h b/src/engines/srpic/particle_pusher.h index b263d6ab..e256f2f4 100644 --- a/src/engines/srpic/particle_pusher.h +++ b/src/engines/srpic/particle_pusher.h @@ -21,6 +21,19 @@ namespace ntt { namespace srpic { + template + void CallPusherWithExternalForce(const kernel::sr::PusherParams& pusher_params, + kernel::sr::PusherArrays& pusher_arrays, + const range_t& range, + const ndfield_t& EB, + const M& metric, + const F& force) { + Kokkos::parallel_for( + "ParticlePusher", + range, + kernel::sr::Pusher_kernel(pusher_params, pusher_arrays, EB, metric, force)); + } + template requires metric::traits::HasD void ParticlePush(Domain& domain, @@ -151,12 +164,12 @@ namespace ntt { pusher_params.ext_force = has_extforce; if (not has_atmosphere and not has_extforce) { - Kokkos::parallel_for("ParticlePusher", - species.rangeActiveParticles(), - kernel::sr::Pusher_kernel(pusher_params, - pusher_arrays, - domain.fields.em, - domain.mesh.metric)); + CallPusherWithExternalForce(pusher_params, + pusher_arrays, + species.rangeActiveParticles(), + domain.fields.em, + domain.mesh.metric, + kernel::sr::NoForce_t {}); } else if (has_atmosphere and not has_extforce) { const auto force = kernel::sr::Force { @@ -164,28 +177,26 @@ namespace ntt { x_surf, ds }; - Kokkos::parallel_for( - "ParticlePusher", + CallPusherWithExternalForce( + pusher_params, + pusher_arrays, species.rangeActiveParticles(), - kernel::sr::Pusher_kernel(pusher_params, - pusher_arrays, - domain.fields.em, - domain.mesh.metric, - force)); + domain.fields.em, + domain.mesh.metric, + force); } else if (not has_atmosphere and has_extforce) { if constexpr (arch::traits::pgen::HasExtForce) { const auto force = kernel::sr::Force { pgen.ext_force }; - Kokkos::parallel_for( - "ParticlePusher", + CallPusherWithExternalForce( + pusher_params, + pusher_arrays, species.rangeActiveParticles(), - kernel::sr::Pusher_kernel(pusher_params, - pusher_arrays, - domain.fields.em, - domain.mesh.metric, - force)); + domain.fields.em, + domain.mesh.metric, + force); } else { raise::Error("External force not implemented", HERE); } @@ -198,14 +209,13 @@ namespace ntt { x_surf, ds }; - Kokkos::parallel_for( - "ParticlePusher", + CallPusherWithExternalForce( + pusher_params, + pusher_arrays, species.rangeActiveParticles(), - kernel::sr::Pusher_kernel(pusher_params, - pusher_arrays, - domain.fields.em, - domain.mesh.metric, - force)); + domain.fields.em, + domain.mesh.metric, + force); } else { raise::Error("External force not implemented", HERE); } From f420ca3c9d55063691d6fdd07461591955280bcf Mon Sep 17 00:00:00 2001 From: haykh Date: Thu, 29 Jan 2026 21:52:10 -0500 Subject: [PATCH 4/9] nix upd + O3 flag added + gui option rm --- CMakeLists.txt | 10 +--------- dev/nix/kokkos.nix | 6 +++--- dev/nix/shell.nix | 4 ++-- 3 files changed, 6 insertions(+), 14 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 109a1f2f..9312d0a1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -47,9 +47,6 @@ set(pgen ${default_pgen} CACHE STRING "Problem generator") -set(gui - ${default_gui} - CACHE BOOL "Use GUI [nttiny]") set(output ${default_output} CACHE BOOL "Enable output") @@ -70,7 +67,7 @@ if(${DEBUG} STREQUAL "OFF") set(CMAKE_BUILD_TYPE Release CACHE STRING "CMake build type") - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DNDEBUG") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DNDEBUG -O3") else() set(CMAKE_BUILD_TYPE Debug @@ -177,11 +174,6 @@ elseif(BENCHMARK) # ------------------------------ Benchmark --------------------------------- # include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/benchmark.cmake) else() - # ----------------------------------- GUI ---------------------------------- # - if(${gui}) - find_or_fetch_dependency(nttiny FALSE QUIET) - endif() - # ------------------------------- Main source ------------------------------ # set_problem_generator(${pgen}) add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/src src) diff --git a/dev/nix/kokkos.nix b/dev/nix/kokkos.nix index ce2ba267..f51d2c68 100644 --- a/dev/nix/kokkos.nix +++ b/dev/nix/kokkos.nix @@ -10,7 +10,7 @@ let pversion = "5.0.1"; compilerPkgs = { "HIP" = with pkgs.rocmPackages; [ - llvm.llvm + clang rocm-core clr rocthrust @@ -41,7 +41,7 @@ let "-D Kokkos_ENABLE_HIP=ON" "-D Kokkos_ARCH_${getArch { }}=ON" "-D GPU_TARGETS=${builtins.replaceStrings [ "amd_" ] [ "" ] (pkgs.lib.toLower (getArch { }))}" - "-D CMAKE_CXX_COMPILER=hipcc" + "-D CMAKE_CXX_COMPILER=clang++" ]; "CUDA" = [ "-D Kokkos_ENABLE_CUDA=ON" @@ -57,7 +57,7 @@ pkgs.stdenv.mkDerivation rec { src = pkgs.fetchgit { url = "https://github.com/kokkos/kokkos/"; rev = "${pversion}"; - sha256 = "sha256-C4DarqnEcdF3+19TPqcM0A9bcQSkKTJkB8b7OkzC7T8="; + sha256 = "sha256-ChpwGBwE7sNovjdAM/iCeOqqwGufKxAh5vQ3qK6aFBU="; }; nativeBuildInputs = with pkgs; [ diff --git a/dev/nix/shell.nix b/dev/nix/shell.nix index 5ff7a389..c1a9dd16 100644 --- a/dev/nix/shell.nix +++ b/dev/nix/shell.nix @@ -29,8 +29,8 @@ let CC = "gcc"; }; HIP = { - CXX = "hipcc"; - CC = "hipcc"; + CXX = "clang++"; + CC = "clang"; }; CUDA = { }; }; From 3281a55834c89738a757b1f619c48a92045e4166 Mon Sep 17 00:00:00 2001 From: haykh Date: Thu, 29 Jan 2026 22:18:41 -0500 Subject: [PATCH 5/9] tests now can be compiled with pgen --- CMakeLists.txt | 23 ++++++++++++++++++----- cmake/benchmark.cmake | 19 +++++++++---------- cmake/report.cmake | 35 ++++++++++++++++++++++++++++++----- cmake/tests.cmake | 7 ------- src/CMakeLists.txt | 9 --------- 5 files changed, 57 insertions(+), 36 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 9312d0a1..de79ee90 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -167,16 +167,29 @@ endif() link_libraries(${DEPENDENCIES}) +set(SRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src/) +add_subdirectory(${SRC_DIR}/global ${CMAKE_CURRENT_BINARY_DIR}/global) +add_subdirectory(${SRC_DIR}/metrics ${CMAKE_CURRENT_BINARY_DIR}/metrics) +add_subdirectory(${SRC_DIR}/kernels ${CMAKE_CURRENT_BINARY_DIR}/kernels) +add_subdirectory(${SRC_DIR}/archetypes ${CMAKE_CURRENT_BINARY_DIR}/archetypes) +add_subdirectory(${SRC_DIR}/framework ${CMAKE_CURRENT_BINARY_DIR}/framework) +add_subdirectory(${SRC_DIR}/output ${CMAKE_CURRENT_BINARY_DIR}/output) + +# ------------------------------- Main source ------------------------------ # +if(NOT ${pgen} STREQUAL ${default_pgen}) + set_problem_generator(${pgen}) + add_subdirectory(${SRC_DIR}/engines ${CMAKE_CURRENT_BINARY_DIR}/engines) + add_subdirectory(${SRC_DIR} ${CMAKE_CURRENT_BINARY_DIR}/src) +endif() + if(TESTS) # ---------------------------------- Tests --------------------------------- # include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/tests.cmake) -elseif(BENCHMARK) +endif() + +if(BENCHMARK) # ------------------------------ Benchmark --------------------------------- # include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/benchmark.cmake) -else() - # ------------------------------- Main source ------------------------------ # - set_problem_generator(${pgen}) - add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/src src) endif() include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/report.cmake) diff --git a/cmake/benchmark.cmake b/cmake/benchmark.cmake index fdd8438e..48be3e99 100644 --- a/cmake/benchmark.cmake +++ b/cmake/benchmark.cmake @@ -2,16 +2,15 @@ set(SRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src) -add_subdirectory(${SRC_DIR}/global ${CMAKE_CURRENT_BINARY_DIR}/global) -add_subdirectory(${SRC_DIR}/metrics ${CMAKE_CURRENT_BINARY_DIR}/metrics) -add_subdirectory(${SRC_DIR}/kernels ${CMAKE_CURRENT_BINARY_DIR}/kernels) -add_subdirectory(${SRC_DIR}/archetypes ${CMAKE_CURRENT_BINARY_DIR}/archetypes) -add_subdirectory(${SRC_DIR}/framework ${CMAKE_CURRENT_BINARY_DIR}/framework) - -if(${output}) - add_subdirectory(${SRC_DIR}/output ${CMAKE_CURRENT_BINARY_DIR}/output) - add_subdirectory(${SRC_DIR}/checkpoint ${CMAKE_CURRENT_BINARY_DIR}/checkpoint) -endif() +# add_subdirectory(${SRC_DIR}/global ${CMAKE_CURRENT_BINARY_DIR}/global) +# add_subdirectory(${SRC_DIR}/metrics ${CMAKE_CURRENT_BINARY_DIR}/metrics) +# add_subdirectory(${SRC_DIR}/kernels ${CMAKE_CURRENT_BINARY_DIR}/kernels) +# add_subdirectory(${SRC_DIR}/archetypes ${CMAKE_CURRENT_BINARY_DIR}/archetypes) +# add_subdirectory(${SRC_DIR}/framework ${CMAKE_CURRENT_BINARY_DIR}/framework) +# +# if(${output}) add_subdirectory(${SRC_DIR}/output +# ${CMAKE_CURRENT_BINARY_DIR}/output) add_subdirectory(${SRC_DIR}/checkpoint +# ${CMAKE_CURRENT_BINARY_DIR}/checkpoint) endif() set(exec benchmark.xc) set(src ${CMAKE_CURRENT_SOURCE_DIR}/benchmark/benchmark.cpp) diff --git a/cmake/report.cmake b/cmake/report.cmake index 8f62ac17..0fbd072c 100644 --- a/cmake/report.cmake +++ b/cmake/report.cmake @@ -8,7 +8,9 @@ if(${PGEN_FOUND}) "${Blue}" PGEN_REPORT 0) -elseif(${TESTS}) +endif() + +if(${TESTS}) set(TEST_NAMES "") foreach(test_dir IN LISTS TEST_DIRECTORIES) get_property( @@ -18,14 +20,38 @@ elseif(${TESTS}) list(APPEND TEST_NAMES ${LOCAL_TEST_NAMES}) endforeach() printchoices( - "Test cases" + "Tests" + "TESTS" + "${ON_OFF_VALUES}" + "ON" + "OFF" + "${Green}" + TESTS_REPORT_1 + 46) + printchoices( + "" "" "${TEST_NAMES}" "" "${ColorReset}" "" - TESTS_REPORT + TESTS_REPORT_2 0) + # remove only first line of TESTS_REPORT_2 + string(REPLACE "\n" ";" TESTS_REPORT_2_LIST "${TESTS_REPORT_2}") + list(REMOVE_AT TESTS_REPORT_2_LIST 0) + string(REPLACE ";" "\n" TESTS_REPORT_2 "${TESTS_REPORT_2_LIST}") + set(TESTS_REPORT "${TESTS_REPORT_1}\n${TESTS_REPORT_2}") +else() + printchoices( + "Tests" + "TESTS" + "${ON_OFF_VALUES}" + "OFF" + "OFF" + "${Green}" + TESTS_REPORT + 46) endif() printchoices( @@ -120,9 +146,8 @@ string(APPEND REPORT_TEXT ${DASHED_LINE_SYMBOL} "\n" "Configurations" "\n") if(${PGEN_FOUND}) string(APPEND REPORT_TEXT " " ${PGEN_REPORT} "\n") -elseif(${TESTS}) - string(APPEND REPORT_TEXT " " ${TESTS_REPORT} "\n") endif() +string(APPEND REPORT_TEXT " " ${TESTS_REPORT} "\n") string( APPEND diff --git a/cmake/tests.cmake b/cmake/tests.cmake index 0eb043f7..44bbce44 100644 --- a/cmake/tests.cmake +++ b/cmake/tests.cmake @@ -3,13 +3,6 @@ enable_testing() set(SRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src) -add_subdirectory(${SRC_DIR}/global ${CMAKE_CURRENT_BINARY_DIR}/global) -add_subdirectory(${SRC_DIR}/metrics ${CMAKE_CURRENT_BINARY_DIR}/metrics) -add_subdirectory(${SRC_DIR}/kernels ${CMAKE_CURRENT_BINARY_DIR}/kernels) -add_subdirectory(${SRC_DIR}/archetypes ${CMAKE_CURRENT_BINARY_DIR}/archetypes) -add_subdirectory(${SRC_DIR}/framework ${CMAKE_CURRENT_BINARY_DIR}/framework) -add_subdirectory(${SRC_DIR}/output ${CMAKE_CURRENT_BINARY_DIR}/output) - set(TEST_DIRECTORIES "") if(NOT ${mpi}) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 31114c33..443dc054 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -24,15 +24,6 @@ set(SRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}) -# dependencies -add_subdirectory(${SRC_DIR}/global ${CMAKE_CURRENT_BINARY_DIR}/global) -add_subdirectory(${SRC_DIR}/metrics ${CMAKE_CURRENT_BINARY_DIR}/metrics) -add_subdirectory(${SRC_DIR}/kernels ${CMAKE_CURRENT_BINARY_DIR}/kernels) -add_subdirectory(${SRC_DIR}/archetypes ${CMAKE_CURRENT_BINARY_DIR}/archetypes) -add_subdirectory(${SRC_DIR}/framework ${CMAKE_CURRENT_BINARY_DIR}/framework) -add_subdirectory(${SRC_DIR}/engines ${CMAKE_CURRENT_BINARY_DIR}/engines) -add_subdirectory(${SRC_DIR}/output ${CMAKE_CURRENT_BINARY_DIR}/output) - set(ENTITY ${PROJECT_NAME}.xc) set(SOURCES ${SRC_DIR}/entity.cpp) From 3cf77b637c40266ba4a251f4eec58affeedaa7b5 Mon Sep 17 00:00:00 2001 From: haykh Date: Fri, 30 Jan 2026 22:05:40 -0500 Subject: [PATCH 6/9] prtldim trait for metrics --- src/metrics/traits.h | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/metrics/traits.h b/src/metrics/traits.h index f9f72808..d499e676 100644 --- a/src/metrics/traits.h +++ b/src/metrics/traits.h @@ -37,6 +37,11 @@ namespace metric { { M::Dim } -> std::convertible_to; }; + template + concept HasPrtlDim = requires { + { M::PrtlDim } -> std::convertible_to; + }; + template concept HasCoordType = requires { { M::CoordType } -> std::convertible_to; From c295fb81e5b190d59cee0e5db89f037efbccc049 Mon Sep 17 00:00:00 2001 From: haykh Date: Fri, 30 Jan 2026 22:06:02 -0500 Subject: [PATCH 7/9] support for 2D arrays in setup input --- src/framework/parameters/parameters.cpp | 46 ++++++++++++++++++++++++- 1 file changed, 45 insertions(+), 1 deletion(-) diff --git a/src/framework/parameters/parameters.cpp b/src/framework/parameters/parameters.cpp index 99e24d78..c9d37c65 100644 --- a/src/framework/parameters/parameters.cpp +++ b/src/framework/parameters/parameters.cpp @@ -265,7 +265,51 @@ namespace ntt { } set("setup." + key, vec); } else if (val_arr[0].is_array()) { - raise::Error("only 1D arrays allowed in [setup]", HERE); + if (val_arr[0].as_array().size() == 0) { + raise::Error("empty inner arrays not allowed in [setup]", HERE); + } else if (val_arr[0][0].is_integer()) { + std::vector> vec; + for (const auto& v1 : val_arr) { + std::vector inner_vec; + for (const auto& v2 : v1.as_array()) { + inner_vec.push_back(v2.as_integer()); + } + vec.push_back(inner_vec); + } + set("setup." + key, vec); + } else if (val_arr[0][0].is_floating()) { + std::vector> vec; + for (const auto& v1 : val_arr) { + std::vector inner_vec; + for (const auto& v2 : v1.as_array()) { + inner_vec.push_back(v2.as_floating()); + } + vec.push_back(inner_vec); + } + set("setup." + key, vec); + } else if (val_arr[0][0].is_boolean()) { + std::vector> vec; + for (const auto& v : val_arr) { + std::vector inner_vec; + for (const auto& v2 : v.as_array()) { + inner_vec.push_back(v2.as_boolean()); + } + vec.push_back(inner_vec); + } + set("setup." + key, vec); + } else if (val_arr[0][0].is_string()) { + std::vector> vec; + for (const auto& v : val_arr) { + std::vector inner_vec; + for (const auto& v2 : v.as_array()) { + inner_vec.push_back(v2.as_string()); + } + vec.push_back(inner_vec); + } + set("setup." + key, vec); + } else if (val_arr[0][0].is_array()) { + raise::Error("up to 2D arrays allowed in [setup]", HERE); + } } else { raise::Error("invalid setup variable type", HERE); } From f13a68dbe4d12e7aaa7ce2e92c743098d38366f1 Mon Sep 17 00:00:00 2001 From: haykh Date: Fri, 30 Jan 2026 22:06:25 -0500 Subject: [PATCH 8/9] GR with no init_flds in pgen support --- src/engines/grpic.hpp | 46 ++++++++++++++++++++++--------------------- 1 file changed, 24 insertions(+), 22 deletions(-) diff --git a/src/engines/grpic.hpp b/src/engines/grpic.hpp index 1db35e0a..99cca8ca 100644 --- a/src/engines/grpic.hpp +++ b/src/engines/grpic.hpp @@ -637,28 +637,30 @@ namespace ntt { } if (dim == in::x1) { if (g != gr_bc::curr) { - Kokkos::parallel_for( - "MatchBoundaries", - CreateRangePolicy(range_min, range_max), - kernel::bc::MatchBoundaries_kernel( - domain.fields.em, - m_pgen.init_flds, - domain.mesh.metric, - xg_edge, - ds, - tags, - domain.mesh.flds_bc())); - Kokkos::parallel_for( - "MatchBoundaries", - CreateRangePolicy(range_min, range_max), - kernel::bc::MatchBoundaries_kernel( - domain.fields.em0, - m_pgen.init_flds, - domain.mesh.metric, - xg_edge, - ds, - tags, - domain.mesh.flds_bc())); + if constexpr (arch::traits::pgen::HasInitFlds) { + Kokkos::parallel_for( + "MatchBoundaries", + CreateRangePolicy(range_min, range_max), + kernel::bc::MatchBoundaries_kernel( + domain.fields.em, + m_pgen.init_flds, + domain.mesh.metric, + xg_edge, + ds, + tags, + domain.mesh.flds_bc())); + Kokkos::parallel_for( + "MatchBoundaries", + CreateRangePolicy(range_min, range_max), + kernel::bc::MatchBoundaries_kernel( + domain.fields.em0, + m_pgen.init_flds, + domain.mesh.metric, + xg_edge, + ds, + tags, + domain.mesh.flds_bc())); + } } else { Kokkos::parallel_for( "AbsorbCurrents", From 73e464098ac83b177e4daf9a6099ef9738b1f9ac Mon Sep 17 00:00:00 2001 From: haykh Date: Fri, 30 Jan 2026 22:09:43 -0500 Subject: [PATCH 9/9] rm unnecessary add_subdirectory in benchmark cmake --- cmake/benchmark.cmake | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/cmake/benchmark.cmake b/cmake/benchmark.cmake index 48be3e99..710ae90a 100644 --- a/cmake/benchmark.cmake +++ b/cmake/benchmark.cmake @@ -2,16 +2,6 @@ set(SRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src) -# add_subdirectory(${SRC_DIR}/global ${CMAKE_CURRENT_BINARY_DIR}/global) -# add_subdirectory(${SRC_DIR}/metrics ${CMAKE_CURRENT_BINARY_DIR}/metrics) -# add_subdirectory(${SRC_DIR}/kernels ${CMAKE_CURRENT_BINARY_DIR}/kernels) -# add_subdirectory(${SRC_DIR}/archetypes ${CMAKE_CURRENT_BINARY_DIR}/archetypes) -# add_subdirectory(${SRC_DIR}/framework ${CMAKE_CURRENT_BINARY_DIR}/framework) -# -# if(${output}) add_subdirectory(${SRC_DIR}/output -# ${CMAKE_CURRENT_BINARY_DIR}/output) add_subdirectory(${SRC_DIR}/checkpoint -# ${CMAKE_CURRENT_BINARY_DIR}/checkpoint) endif() - set(exec benchmark.xc) set(src ${CMAKE_CURRENT_SOURCE_DIR}/benchmark/benchmark.cpp)