From e9c927b7240b1c6343ba5badfd97c06e93ed0c80 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Fri, 3 Apr 2026 03:36:02 +0000 Subject: [PATCH 01/30] first attempt at moving window --- pgens/piston_shock/pgen.hpp | 195 +++++++++++++++++++------- src/archetypes/moving_window.h | 247 +++++++++++++++++++++++++++++++++ 2 files changed, 392 insertions(+), 50 deletions(-) create mode 100644 src/archetypes/moving_window.h diff --git a/pgens/piston_shock/pgen.hpp b/pgens/piston_shock/pgen.hpp index 86839f0c9..d5e244023 100644 --- a/pgens/piston_shock/pgen.hpp +++ b/pgens/piston_shock/pgen.hpp @@ -4,27 +4,72 @@ #include "enums.h" #include "global.h" -#include "arch/traits.h" -#include "archetypes/piston.h" +#include "utils/error.h" +#include "utils/numeric.h" -#include "archetypes/particle_injector.h" +#include "archetypes/field_setter.h" #include "archetypes/problem_generator.h" -#include "archetypes/energy_dist.h" -#include "framework/domain/domain.h" -#include "archetypes/spatial_dist.h" -#include "framework/domain/metadomain.h" +#include "archetypes/traits.h" #include "archetypes/utils.h" +#include "archetypes/piston.h" +#include "archetypes/moving_window.h" +#include "framework/domain/metadomain.h" -#include - - +#include +#include namespace user { using namespace ntt; + template + struct InitFields { + /* + Sets up magnetic and electric field components for the simulation. + Must satisfy E = -v x B for Lorentz Force to be zero. + + @param bmag: magnetic field scaling + @param btheta: magnetic field polar angle + @param bphi: magnetic field azimuthal angle + @param drift_ux: drift velocity in the x direction + */ + InitFields(real_t bmag, real_t btheta, real_t bphi, real_t drift_ux) + : Bmag { bmag } + , Btheta { btheta * static_cast(convert::deg2rad) } + , Bphi { bphi * static_cast(convert::deg2rad) } + , Vx { drift_ux } {} + + // magnetic field components + Inline auto bx1(const coord_t&) const -> real_t { + return Bmag * math::cos(Btheta); + } + + Inline auto bx2(const coord_t&) const -> real_t { + return Bmag * math::sin(Btheta) * math::sin(Bphi); + } + + Inline auto bx3(const coord_t&) const -> real_t { + return Bmag * math::sin(Btheta) * math::cos(Bphi); + } + + // electric field components + Inline auto ex1(const coord_t&) const -> real_t { + return ZERO; + } + + Inline auto ex2(const coord_t&) const -> real_t { + return -Vx * Bmag * math::sin(Btheta) * math::cos(Bphi); + } + + Inline auto ex3(const coord_t&) const -> real_t { + return Vx * Bmag * math::sin(Btheta) * math::sin(Bphi); + } + + private: + const real_t Btheta, Bphi, Vx, Bmag; + }; + template struct PGen : public arch::ProblemGenerator { - // compatibility traits for the problem generator static constexpr auto engines { arch::traits::pgen::compatible_with::value @@ -41,42 +86,49 @@ namespace user { using arch::ProblemGenerator::C; using arch::ProblemGenerator::params; - const Metadomain& global_domain; + Metadomain& global_domain; // domain properties const real_t global_xmin, global_xmax; + // gas properties + const real_t temperature, temperature_ratio; + // injector properties + const real_t dt; + // magnetic field properties + real_t Btheta, Bphi, Bmag; + InitFields init_flds; - // plasma - const real_t temperature, temperature_ratio; // piston properties - const real_t piston_velocity, piston_initial_position; + const real_t piston_velocity; + int i_piston; + real_t di_piston, piston_position; - inline PGen(const SimulationParams& p, const Metadomain& global_domain) + inline PGen(const SimulationParams& p, Metadomain& global_domain) : arch::ProblemGenerator { p } - , global_domain { global_domain } + , global_domain { global_domain } , global_xmin { global_domain.mesh().extent(in::x1).first } , global_xmax { global_domain.mesh().extent(in::x1).second } - , piston_velocity {p.template get("setup.piston_velocity", 0.0)} - , piston_initial_position {p.template get("setup.piston_initial_position", 0.0)} - , temperature { p.template get("setup.temperature", 0.0) } - , temperature_ratio { p.template get("setup.temperature_ratio") }{} - - inline void InitPrtls(Domain& local_domain) { - real_t xg_min = global_xmin+piston_initial_position; - real_t xg_max = global_xmax ; - - // define box to inject into - boundaries_t box; - // loop over all dimensions - for (auto d { 0u }; d < (unsigned int)M::Dim; ++d) { - // compute the range for the x-direction - if (d == static_cast(in::x1)) { - box.push_back({ xg_min, xg_max }); - } else { - // inject into full range in other directions - box.push_back(Range::All); - } - } + , drift_ux { p.template get("setup.drift_ux") } + , temperature { p.template get("setup.temperature") } + , temperature_ratio { p.template get("setup.temperature_ratio") } + , Bmag { p.template get("setup.Bmag", ZERO) } + , Btheta { p.template get("setup.Btheta", ZERO) } + , Bphi { p.template get("setup.Bphi", ZERO) } + , init_flds { Bmag, Btheta, Bphi, drift_ux } + , dt { p.template get("algorithms.timestep.dt") } {} + + inline PGen() {} + + auto MatchFields(simtime_t) const -> InitFields { + return init_flds; + } + + inline void InitPrtls(Domain& domain) { + + // set initial position of piston + i_piston = 0; + di_piston = ZERO; + piston_position = ZERO; // define temperatures of species const auto temperatures = std::make_pair(temperature, @@ -88,28 +140,75 @@ namespace user { // inject particles arch::InjectUniformMaxwellians(params, - local_domain, + domain, ONE, temperatures, { 1, 2 }, drifts, - false, - box); + false); } + void CustomPostStep(timestep_t step, simtime_t time, Domain& domain) { + + di_piston += piston_velocity * dt; + if (di_piston >= domain.mesh().CellSize(in::x1)) { + i_piston += static_cast(di_piston / domain.mesh().CellSize(in::x1)); + di_piston = std::fmod(di_piston, domain.mesh().CellSize(in::x1)); + } + + // compute current position of piston + piston_position = static_cast(i_piston) + di_piston; + + // check if the window should be moved + if (i_piston > window_update_frequency) { + + // move the window and all fields and particles in it + arch::MoveWindow(in::x1, domain, window_update_frequency); + + // synch ghost zones after moving the window + global_domain.CommunicateFields(domain, Comm::E | Comm::B); + + /* + Inject slab of fresh plasma + */ + const real_t xmax = global_xmax; + const real_t xmin = xmax - window_update_frequency * domain.mesh().CellSize(in::x1); + // define box to inject into + boundaries_t inj_box; + // loop over all dimension + for (auto d = 0u; d < M::Dim; ++d) { + if (d == 0) { + inj_box.push_back({ xmin, xmax }); + } else { + inj_box.push_back(Range::All); + } + } + + // same maxwell distribution as above + const auto temperatures = std::make_pair(temperature, + temperature_ratio * temperature); + const auto drifts = std::make_pair( + std::vector { ZERO, ZERO, ZERO }, + std::vector { ZERO, ZERO, ZERO }); + arch::InjectUniformMaxwellians(params, + domain, + ONE, + temperatures, + { 1, 2 }, + drifts, + false, + inj_box); + } + } struct CustomPrtlUpdate { real_t piston_position_current; // current position of piston real_t piston_velocity_current; real_t xg_max; - bool is_left; - bool massive; - - template Inline void operator()(index_t p, Coord& xp, PusherKernel& pusher) const { @@ -124,17 +223,13 @@ namespace user { if (arch::CrossesPiston(p, pusher, piston_position_use, piston_velocity_current, is_left)){ arch::PistonUpdate(p, pusher, piston_position_use, piston_velocity_current, massive); } - - } }; template auto CustomParticleUpdate(simtime_t time, spidx_t sp, D& domain) const { - return CustomPrtlUpdate { piston_initial_position+static_cast(time)*piston_velocity, piston_velocity, global_domain.mesh().extent(in::x1).second, true, true}; - + return CustomPrtlUpdate { piston_pos, piston_velocity, global_domain.mesh().extent(in::x1).second, true, true}; }; - }; } // namespace user diff --git a/src/archetypes/moving_window.h b/src/archetypes/moving_window.h new file mode 100644 index 000000000..24174eff2 --- /dev/null +++ b/src/archetypes/moving_window.h @@ -0,0 +1,247 @@ +/** + * @file archetypes/piston.h + * @brief Piston functions for implementing the piston in the CustomPrtlUpdate in pgran + * @implements + * - arch::PistonUpdate<> -> void + * - arch::CrossesPiston<> -> Bool + * @namespaces: + * - arch:: + */ + +#ifndef ARCHETYPES_WINDOW_H +#define ARCHETYPES_WINDOW_H + +#include "enums.h" +#include "global.h" + +#include "archetypes/energy_dist.h" +#include "framework/domain/domain.h" +#include "framework/parameters/parameters.h" +#include "framework/domain/metadomain.h" + +#include + + +namespace arch { + + /** + * @brief Updates particle position and fields in the moving window. + + */ + template + Inline void MoveWindow(dir::direction_t direction, + Domain& domain, + const int window_shift) { + + if constexpr (M::CoordType != Coord::Cart) { + (void)direction; + (void)domain; + (void)tags; + raise::Error( + "Moving window only applicable to cartesian coordinates", + HERE); + } else { + + const auto sign = direction.get_sign(); + const auto dim = direction.get_dim(); + + /* + move particles in the window back by the window size + */ + const auto& mesh = domain.mesh; + + // loop over particle species + for (auto s { 0u }; s < 2; ++s) { + // get particle properties + auto& species = domain.species[s]; + auto i1 = species.i1; + + Kokkos::parallel_for( + "MoveParticles", + species.rangeActiveParticles(), + Lambda(index_t p) { + // shift particle position back by window update frequency + i1(p) -= window_shift; + }); + } + + // shift fields in the window back by the window size + 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(domain.mesh.n_all(dd) - window_shift); + } 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); + + // loop range for shifting fields + 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); + } + + if (dir == in::x1) { + Kokkos::parallel_for( + "ShiftFields", + range, + FieldShift_kernel( + domain.fields.em, + window_shift, + tags)); + } else if (dir == in::x2) { + if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { + Kokkos::parallel_for( + "ShiftFields", + range, + FieldShift_kernel( + domain.fields.em, + window_shift, + tags)); + } else { + raise::Error("Invalid dimension", HERE); + } + } else { + if constexpr (M::Dim == Dim::_3D) { + Kokkos::parallel_for( + "ShiftFields", + range, + FieldShift_kernel( + domain.fields.em, + window_shift, + tags)); + } else { + raise::Error("Invalid dimension", HERE); + } + } + + } + } + + template + struct FieldShift_kernel { + static_assert(static_cast(o) < static_cast(D), + "Invalid component index"); + + ndfield_t Fld; + const BCTags tags; + const int window_shift; + + FieldShift_kernel(ndfield_t Fld, const int window_shift, BCTags tags) + : Fld { Fld } + , window_shift { window_shift } + , tags { tags } {} + + Inline void operator()(index_t i1) const { + if constexpr (D == Dim::_1D) { + if (tags & BC::E) { + Fld(i1, em::ex1) = Fld(i1 + window_shift, em::ex1); + Fld(i1, em::ex2) = Fld(i1 + window_shift, em::ex2); + Fld(i1, em::ex3) = Fld(i1 + window_shift, em::ex3); + } + if (tags & BC::B) { + Fld(i1, em::bx1) = Fld(i1 + window_shift, em::bx1); + Fld(i1, em::bx2) = Fld(i1 + window_shift, em::bx2); + Fld(i1, em::bx3) = Fld(i1 + window_shift, em::bx3); + } + } else { + raise::KernelError( + HERE, + "FieldShift_kernel: 1D implementation called for D != 1"); + } + } + + Inline void operator()(index_t i1, index_t i2) const { + if constexpr (D == Dim::_2D) { + if constexpr (o == in::x1) { + if (tags & BC::E) { + Fld(i1, i2, em::ex1) = Fld(i1 + window_shift, i2, em::ex1); + Fld(i1, i2, em::ex2) = Fld(i1 + window_shift, i2, em::ex2); + Fld(i1, i2, em::ex3) = Fld(i1 + window_shift, i2, em::ex3); + } + if (tags & BC::B) { + Fld(i1, i2, em::bx1) = Fld(i1 + window_shift, i2, em::bx1); + Fld(i1, i2, em::bx2) = Fld(i1 + window_shift, i2, em::bx2); + Fld(i1, i2, em::bx3) = Fld(i1 + window_shift, i2, em::bx3); + } + } else if constexpr (o == in::x2) { + if (tags & BC::E) { + Fld(i1, i2, em::ex1) = Fld(i1, i2 + window_shift, em::ex1); + Fld(i1, i2, em::ex2) = Fld(i1, i2 + window_shift, em::ex2); + Fld(i1, i2, em::ex3) = Fld(i1, i2 + window_shift, em::ex3); + } + if (tags & BC::B) { + Fld(i1, i2, em::bx1) = Fld(i1, i2 + window_shift, em::bx1); + Fld(i1, i2, em::bx2) = Fld(i1, i2 + window_shift, em::bx2); + Fld(i1, i2, em::bx3) = Fld(i1, i2 + window_shift, em::bx3); + } + } + } else { + raise::KernelError( + HERE, + "FieldShift_kernel: 2D implementation called for D != 2"); + } + } + + Inline void operator()(index_t i1, index_t i2, index_t i3) const { + if constexpr (D == Dim::_3D) { + if constexpr (o == in::x1) { + if (tags & BC::E) { + Fld(i1, i2, i3, em::ex1) = Fld(i1 + window_shift, i2, i3, em::ex1); + Fld(i1, i2, i3, em::ex2) = Fld(i1 + window_shift, i2, i3, em::ex2); + Fld(i1, i2, i3, em::ex3) = Fld(i1 + window_shift, i2, i3, em::ex3); + } + if (tags & BC::B) { + Fld(i1, i2, i3, em::bx1) = Fld(i1 + window_shift, i2, i3, em::bx1); + Fld(i1, i2, i3, em::bx2) = Fld(i1 + window_shift, i2, i3, em::bx2); + Fld(i1, i2, i3, em::bx3) = Fld(i1 + window_shift, i2, i3, em::bx3); + } + } else if constexpr (o == in::x2) { + if (tags & BC::E) { + Fld(i1, i2, i3, em::ex1) = Fld(i1, i2 + window_shift, i3, em::ex1); + Fld(i1, i2, i3, em::ex2) = Fld(i1, i2 + window_shift, i3, em::ex2); + Fld(i1, i2, i3, em::ex3) = Fld(i1, i2 + window_shift, i3, em::ex3); + } + if (tags & BC::B) { + Fld(i1, i2, i3, em::bx1) = Fld(i1, i2 + window_shift, i3, em::bx1); + Fld(i1, i2, i3, em::bx2) = Fld(i1, i2 + window_shift, i3, em::bx2); + Fld(i1, i2, i3, em::bx3) = Fld(i1, i2 + window_shift, i3, em::bx3); + } + } else { + if (tags & BC::E) { + Fld(i1, i2, i3, em::ex1) = Fld(i1, i2, i3 + window_shift, em::ex1); + Fld(i1, i2, i3, em::ex2) = Fld(i1, i2, i3 + window_shift, em::ex2); + Fld(i1, i2, i3, em::ex3) = Fld(i1, i2, i3 + window_shift, em::ex3); + } + if (tags & BC::B) { + Fld(i1, i2, i3, em::bx1) = Fld(i1, i2, i3 + window_shift, em::bx1); + Fld(i1, i2, i3, em::bx2) = Fld(i1, i2, i3 + window_shift, em::bx2); + Fld(i1, i2, i3, em::bx3) = Fld(i1, i2, i3 + window_shift, em::bx3); + } + } + } else { + raise::KernelError( + HERE, + "FieldShift_kernel: 3D implementation called for D != 3"); + } + } + }; +} // namespace arch + +#endif // ARCHETYPES_UTILS_H From 76680fe7ac49f4f69d8dfdd73fa35c50e9de348e Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Wed, 8 Apr 2026 02:30:27 +0000 Subject: [PATCH 02/30] moving window almost finished --- pgens/piston_shock/pgen.hpp | 59 +++--- src/archetypes/moving_window.h | 335 ++++++++++++++++++--------------- 2 files changed, 226 insertions(+), 168 deletions(-) diff --git a/pgens/piston_shock/pgen.hpp b/pgens/piston_shock/pgen.hpp index d5e244023..20098b7d7 100644 --- a/pgens/piston_shock/pgen.hpp +++ b/pgens/piston_shock/pgen.hpp @@ -103,19 +103,23 @@ namespace user { int i_piston; real_t di_piston, piston_position; + // window properties + const int window_update_frequency; + inline PGen(const SimulationParams& p, Metadomain& global_domain) : arch::ProblemGenerator { p } , global_domain { global_domain } , global_xmin { global_domain.mesh().extent(in::x1).first } , global_xmax { global_domain.mesh().extent(in::x1).second } - , drift_ux { p.template get("setup.drift_ux") } , temperature { p.template get("setup.temperature") } , temperature_ratio { p.template get("setup.temperature_ratio") } , Bmag { p.template get("setup.Bmag", ZERO) } , Btheta { p.template get("setup.Btheta", ZERO) } , Bphi { p.template get("setup.Bphi", ZERO) } - , init_flds { Bmag, Btheta, Bphi, drift_ux } - , dt { p.template get("algorithms.timestep.dt") } {} + , init_flds { Bmag, Btheta, Bphi, ZERO } + , dt { p.template get("algorithms.timestep.dt") } + , window_update_frequency { p.template get("setup.window_update_frequency", N_GHOSTS) } + , piston_velocity { p.template get("setup.piston_velocity") } {} inline PGen() {} @@ -150,20 +154,23 @@ namespace user { void CustomPostStep(timestep_t step, simtime_t time, Domain& domain) { - di_piston += piston_velocity * dt; - if (di_piston >= domain.mesh().CellSize(in::x1)) { - i_piston += static_cast(di_piston / domain.mesh().CellSize(in::x1)); - di_piston = std::fmod(di_piston, domain.mesh().CellSize(in::x1)); - } - - // compute current position of piston - piston_position = static_cast(i_piston) + di_piston; + /* + update piston position + */ + const real_t gamma_piston = ONE / math::sqrt(ONE - SQR(piston_velocity)); + const real_t v_piston = gamma_piston * piston_velocity; + // piston movement over timestep + di_piston += v_piston * dt; + // check if the piston has moved to the next cell + i_piston += static_cast(di_piston >= ONE); + // keep track of how much the piston has moved into the next cell + di_piston -= (di_piston >= ONE); // check if the window should be moved - if (i_piston > window_update_frequency) { + if (i_piston >= window_update_frequency) { // move the window and all fields and particles in it - arch::MoveWindow(in::x1, domain, window_update_frequency); + arch::MoveWindow(domain, window_update_frequency); // synch ghost zones after moving the window global_domain.CommunicateFields(domain, Comm::E | Comm::B); @@ -171,8 +178,9 @@ namespace user { /* Inject slab of fresh plasma */ + const real_t cell_size = ZERO; // ToDo: get cell size from global domain const real_t xmax = global_xmax; - const real_t xmin = xmax - window_update_frequency * domain.mesh().CellSize(in::x1); + const real_t xmin = xmax - window_update_frequency * cell_size; // define box to inject into boundaries_t inj_box; // loop over all dimension @@ -198,13 +206,19 @@ namespace user { drifts, false, inj_box); + + i_piston -= window_update_frequency; } + + // compute current position of piston + // ToDo: check if this can be used as a global variable to avoid recomputing in piston update + piston_position = static_cast(i_piston) + di_piston; } struct CustomPrtlUpdate { - real_t piston_position_current; // current position of piston - real_t piston_velocity_current; + real_t x_piston; + real_t v_piston; real_t xg_max; bool is_left; bool massive; @@ -214,21 +228,24 @@ namespace user { real_t piston_position_use; - if (piston_position_current(p, pusher, piston_position_use, piston_velocity_current, is_left)){ - arch::PistonUpdate(p, pusher, piston_position_use, piston_velocity_current, massive); + if (arch::CrossesPiston(p, pusher, piston_position_use, v_piston, is_left)){ + arch::PistonUpdate(p, pusher, piston_position_use, v_piston, massive); } } }; template auto CustomParticleUpdate(simtime_t time, spidx_t sp, D& domain) const { - return CustomPrtlUpdate { piston_pos, piston_velocity, global_domain.mesh().extent(in::x1).second, true, true}; + const real_t gamma_piston = ONE / math::sqrt(ONE - SQR(piston_velocity)); + const real_t v_piston = gamma_piston * piston_velocity; + const real_t piston_pos = std::fmod(time * piston_velocity, window_update_frequency); + return CustomPrtlUpdate { piston_pos, v_piston, global_domain.mesh().extent(in::x1).second, true, true}; }; }; diff --git a/src/archetypes/moving_window.h b/src/archetypes/moving_window.h index 24174eff2..eff0c172d 100644 --- a/src/archetypes/moving_window.h +++ b/src/archetypes/moving_window.h @@ -24,126 +24,17 @@ namespace arch { - /** - * @brief Updates particle position and fields in the moving window. - - */ - template - Inline void MoveWindow(dir::direction_t direction, - Domain& domain, - const int window_shift) { - - if constexpr (M::CoordType != Coord::Cart) { - (void)direction; - (void)domain; - (void)tags; - raise::Error( - "Moving window only applicable to cartesian coordinates", - HERE); - } else { - - const auto sign = direction.get_sign(); - const auto dim = direction.get_dim(); - - /* - move particles in the window back by the window size - */ - const auto& mesh = domain.mesh; - - // loop over particle species - for (auto s { 0u }; s < 2; ++s) { - // get particle properties - auto& species = domain.species[s]; - auto i1 = species.i1; - - Kokkos::parallel_for( - "MoveParticles", - species.rangeActiveParticles(), - Lambda(index_t p) { - // shift particle position back by window update frequency - i1(p) -= window_shift; - }); - } - - // shift fields in the window back by the window size - 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(domain.mesh.n_all(dd) - window_shift); - } 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); - - // loop range for shifting fields - 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); - } - - if (dir == in::x1) { - Kokkos::parallel_for( - "ShiftFields", - range, - FieldShift_kernel( - domain.fields.em, - window_shift, - tags)); - } else if (dir == in::x2) { - if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - Kokkos::parallel_for( - "ShiftFields", - range, - FieldShift_kernel( - domain.fields.em, - window_shift, - tags)); - } else { - raise::Error("Invalid dimension", HERE); - } - } else { - if constexpr (M::Dim == Dim::_3D) { - Kokkos::parallel_for( - "ShiftFields", - range, - FieldShift_kernel( - domain.fields.em, - window_shift, - tags)); - } else { - raise::Error("Invalid dimension", HERE); - } - } - - } - } - template struct FieldShift_kernel { static_assert(static_cast(o) < static_cast(D), "Invalid component index"); ndfield_t Fld; + ndfield_t backup_Fld; const BCTags tags; const int window_shift; - FieldShift_kernel(ndfield_t Fld, const int window_shift, BCTags tags) + FieldShift_kernel(ndfield_t Fld, ndfield_t backup_Fld, const int window_shift, BCTags tags) : Fld { Fld } , window_shift { window_shift } , tags { tags } {} @@ -151,14 +42,14 @@ namespace arch { Inline void operator()(index_t i1) const { if constexpr (D == Dim::_1D) { if (tags & BC::E) { - Fld(i1, em::ex1) = Fld(i1 + window_shift, em::ex1); - Fld(i1, em::ex2) = Fld(i1 + window_shift, em::ex2); - Fld(i1, em::ex3) = Fld(i1 + window_shift, em::ex3); + Fld(i1, em::ex1) = backup_Fld(i1 + window_shift, em::ex1); + Fld(i1, em::ex2) = backup_Fld(i1 + window_shift, em::ex2); + Fld(i1, em::ex3) = backup_Fld(i1 + window_shift, em::ex3); } if (tags & BC::B) { - Fld(i1, em::bx1) = Fld(i1 + window_shift, em::bx1); - Fld(i1, em::bx2) = Fld(i1 + window_shift, em::bx2); - Fld(i1, em::bx3) = Fld(i1 + window_shift, em::bx3); + Fld(i1, em::bx1) = backup_Fld(i1 + window_shift, em::bx1); + Fld(i1, em::bx2) = backup_Fld(i1 + window_shift, em::bx2); + Fld(i1, em::bx3) = backup_Fld(i1 + window_shift, em::bx3); } } else { raise::KernelError( @@ -171,25 +62,25 @@ namespace arch { if constexpr (D == Dim::_2D) { if constexpr (o == in::x1) { if (tags & BC::E) { - Fld(i1, i2, em::ex1) = Fld(i1 + window_shift, i2, em::ex1); - Fld(i1, i2, em::ex2) = Fld(i1 + window_shift, i2, em::ex2); - Fld(i1, i2, em::ex3) = Fld(i1 + window_shift, i2, em::ex3); + Fld(i1, i2, em::ex1) = backup_Fld(i1 + window_shift, i2, em::ex1); + Fld(i1, i2, em::ex2) = backup_Fld(i1 + window_shift, i2, em::ex2); + Fld(i1, i2, em::ex3) = backup_Fld(i1 + window_shift, i2, em::ex3); } if (tags & BC::B) { - Fld(i1, i2, em::bx1) = Fld(i1 + window_shift, i2, em::bx1); - Fld(i1, i2, em::bx2) = Fld(i1 + window_shift, i2, em::bx2); - Fld(i1, i2, em::bx3) = Fld(i1 + window_shift, i2, em::bx3); + Fld(i1, i2, em::bx1) = backup_Fld(i1 + window_shift, i2, em::bx1); + Fld(i1, i2, em::bx2) = backup_Fld(i1 + window_shift, i2, em::bx2); + Fld(i1, i2, em::bx3) = backup_Fld(i1 + window_shift, i2, em::bx3); } } else if constexpr (o == in::x2) { if (tags & BC::E) { - Fld(i1, i2, em::ex1) = Fld(i1, i2 + window_shift, em::ex1); - Fld(i1, i2, em::ex2) = Fld(i1, i2 + window_shift, em::ex2); - Fld(i1, i2, em::ex3) = Fld(i1, i2 + window_shift, em::ex3); + Fld(i1, i2, em::ex1) = backup_Fld(i1, i2 + window_shift, em::ex1); + Fld(i1, i2, em::ex2) = backup_Fld(i1, i2 + window_shift, em::ex2); + Fld(i1, i2, em::ex3) = backup_Fld(i1, i2 + window_shift, em::ex3); } if (tags & BC::B) { - Fld(i1, i2, em::bx1) = Fld(i1, i2 + window_shift, em::bx1); - Fld(i1, i2, em::bx2) = Fld(i1, i2 + window_shift, em::bx2); - Fld(i1, i2, em::bx3) = Fld(i1, i2 + window_shift, em::bx3); + Fld(i1, i2, em::bx1) = backup_Fld(i1, i2 + window_shift, em::bx1); + Fld(i1, i2, em::bx2) = backup_Fld(i1, i2 + window_shift, em::bx2); + Fld(i1, i2, em::bx3) = backup_Fld(i1, i2 + window_shift, em::bx3); } } } else { @@ -203,36 +94,36 @@ namespace arch { if constexpr (D == Dim::_3D) { if constexpr (o == in::x1) { if (tags & BC::E) { - Fld(i1, i2, i3, em::ex1) = Fld(i1 + window_shift, i2, i3, em::ex1); - Fld(i1, i2, i3, em::ex2) = Fld(i1 + window_shift, i2, i3, em::ex2); - Fld(i1, i2, i3, em::ex3) = Fld(i1 + window_shift, i2, i3, em::ex3); + Fld(i1, i2, i3, em::ex1) = backup_Fld(i1 + window_shift, i2, i3, em::ex1); + Fld(i1, i2, i3, em::ex2) = backup_Fld(i1 + window_shift, i2, i3, em::ex2); + Fld(i1, i2, i3, em::ex3) = backup_Fld(i1 + window_shift, i2, i3, em::ex3); } if (tags & BC::B) { - Fld(i1, i2, i3, em::bx1) = Fld(i1 + window_shift, i2, i3, em::bx1); - Fld(i1, i2, i3, em::bx2) = Fld(i1 + window_shift, i2, i3, em::bx2); - Fld(i1, i2, i3, em::bx3) = Fld(i1 + window_shift, i2, i3, em::bx3); + Fld(i1, i2, i3, em::bx1) = backup_Fld(i1 + window_shift, i2, i3, em::bx1); + Fld(i1, i2, i3, em::bx2) = backup_Fld(i1 + window_shift, i2, i3, em::bx2); + Fld(i1, i2, i3, em::bx3) = backup_Fld(i1 + window_shift, i2, i3, em::bx3); } } else if constexpr (o == in::x2) { if (tags & BC::E) { - Fld(i1, i2, i3, em::ex1) = Fld(i1, i2 + window_shift, i3, em::ex1); - Fld(i1, i2, i3, em::ex2) = Fld(i1, i2 + window_shift, i3, em::ex2); - Fld(i1, i2, i3, em::ex3) = Fld(i1, i2 + window_shift, i3, em::ex3); + Fld(i1, i2, i3, em::ex1) = backup_Fld(i1, i2 + window_shift, i3, em::ex1); + Fld(i1, i2, i3, em::ex2) = backup_Fld(i1, i2 + window_shift, i3, em::ex2); + Fld(i1, i2, i3, em::ex3) = backup_Fld(i1, i2 + window_shift, i3, em::ex3); } if (tags & BC::B) { - Fld(i1, i2, i3, em::bx1) = Fld(i1, i2 + window_shift, i3, em::bx1); - Fld(i1, i2, i3, em::bx2) = Fld(i1, i2 + window_shift, i3, em::bx2); - Fld(i1, i2, i3, em::bx3) = Fld(i1, i2 + window_shift, i3, em::bx3); + Fld(i1, i2, i3, em::bx1) = backup_Fld(i1, i2 + window_shift, i3, em::bx1); + Fld(i1, i2, i3, em::bx2) = backup_Fld(i1, i2 + window_shift, i3, em::bx2); + Fld(i1, i2, i3, em::bx3) = backup_Fld(i1, i2 + window_shift, i3, em::bx3); } } else { if (tags & BC::E) { - Fld(i1, i2, i3, em::ex1) = Fld(i1, i2, i3 + window_shift, em::ex1); - Fld(i1, i2, i3, em::ex2) = Fld(i1, i2, i3 + window_shift, em::ex2); - Fld(i1, i2, i3, em::ex3) = Fld(i1, i2, i3 + window_shift, em::ex3); + Fld(i1, i2, i3, em::ex1) = backup_Fld(i1, i2, i3 + window_shift, em::ex1); + Fld(i1, i2, i3, em::ex2) = backup_Fld(i1, i2, i3 + window_shift, em::ex2); + Fld(i1, i2, i3, em::ex3) = backup_Fld(i1, i2, i3 + window_shift, em::ex3); } if (tags & BC::B) { - Fld(i1, i2, i3, em::bx1) = Fld(i1, i2, i3 + window_shift, em::bx1); - Fld(i1, i2, i3, em::bx2) = Fld(i1, i2, i3 + window_shift, em::bx2); - Fld(i1, i2, i3, em::bx3) = Fld(i1, i2, i3 + window_shift, em::bx3); + Fld(i1, i2, i3, em::bx1) = backup_Fld(i1, i2, i3 + window_shift, em::bx1); + Fld(i1, i2, i3, em::bx2) = backup_Fld(i1, i2, i3 + window_shift, em::bx2); + Fld(i1, i2, i3, em::bx3) = backup_Fld(i1, i2, i3 + window_shift, em::bx3); } } } else { @@ -242,6 +133,156 @@ namespace arch { } } }; + + /** + * @brief Updates particle position and fields in the moving window. + + */ + template + Inline void MoveWindow(Domain& domain, + const int window_shift) { + + if constexpr (M::CoordType != Coord::Cart) { + raise::Error( + "Moving window only applicable to cartesian coordinates", + HERE); + } else { + + /* + move particles in the window back by the window size + */ + const auto nspec = domain.species.size(); + const auto& mesh = domain.mesh; + if (o == in::x1) { + // loop over particle species + for (auto s { 0u }; s < nspec; ++s) { + // get particle properties + auto& species = domain.species[s]; + auto i1 = species.i1; + Kokkos::parallel_for( + "MoveParticles", + species.rangeActiveParticles(), + Lambda(index_t p) { + // shift particle position back by window update frequency + i1(p) -= window_shift; + } + ); + } + } else if (o == in::x2) { + if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { + // loop over particle species + for (auto s { 0u }; s < nspec; ++s) { + // get particle properties + auto& species = domain.species[s]; + auto i2 = species.i2; + Kokkos::parallel_for( + "MoveParticles", + species.rangeActiveParticles(), + Lambda(index_t p) { + // shift particle position back by window update frequency + i2(p) -= window_shift; + } + ); + } + } else { + raise::Error("Invalid dimension", HERE); + } + } else if (o == in::x3) { + if constexpr (M::Dim == Dim::_3D) { + // loop over particle species + for (auto s { 0u }; s < nspec; ++s) { + // get particle properties + auto& species = domain.species[s]; + auto i3 = species.i3; + Kokkos::parallel_for( + "MoveParticles", + species.rangeActiveParticles(), + Lambda(index_t p) { + // shift particle position back by window update frequency + i3(p) -= window_shift; + } + ); + } + } else { + raise::Error("Invalid direction", HERE); + } + } else { + raise::Error("Invalid direction", HERE); + } + + // shift fields in the window back by the window size + 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 (o == dd) { + xi_min.push_back(0); + xi_max.push_back(domain.mesh.n_all(dd) - window_shift); + } 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); + + // loop range for shifting fields + 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); + } + + // copy fields to backup before shifting + Kokkos::deep_copy(domain.fields.bckp, domain.fields.em); + + if (o == in::x1) { + Kokkos::parallel_for( + "ShiftFields", + range, + FieldShift_kernel( + domain.fields.em, + domain.fields.bckp, + window_shift, + BC::B | BC::E)); + } else if (o == in::x2) { + if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { + Kokkos::parallel_for( + "ShiftFields", + range, + FieldShift_kernel( + domain.fields.em, + domain.fields.bckp, + window_shift, + BC::B | BC::E)); + } else { + raise::Error("Invalid dimension", HERE); + } + } else { + if constexpr (M::Dim == Dim::_3D) { + Kokkos::parallel_for( + "ShiftFields", + range, + FieldShift_kernel( + domain.fields.em, + domain.fields.bckp, + window_shift, + BC::B | BC::E)); + } else { + raise::Error("Invalid dimension", HERE); + } + } + } + } } // namespace arch #endif // ARCHETYPES_UTILS_H From 1413b4c93b3a46059c63dd83da331c1dadcfc881 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Wed, 8 Apr 2026 17:11:33 +0000 Subject: [PATCH 03/30] use globally defined piston position and velocity to avoid duplicate code --- pgens/piston_shock/pgen.hpp | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/pgens/piston_shock/pgen.hpp b/pgens/piston_shock/pgen.hpp index 117ea905b..6f3397c1e 100644 --- a/pgens/piston_shock/pgen.hpp +++ b/pgens/piston_shock/pgen.hpp @@ -119,7 +119,8 @@ namespace user { , init_flds { Bmag, Btheta, Bphi, ZERO } , dt { p.template get("algorithms.timestep.dt") } , window_update_frequency { p.template get("setup.window_update_frequency", N_GHOSTS) } - , piston_velocity { p.template get("setup.piston_velocity") } {} + , piston_velocity { p.template get("setup.piston_velocity", ZERO) + / math::sqrt(ONE - SQR(p.template get("setup.piston_velocity", ZERO)))} {} inline PGen() {} @@ -156,10 +157,8 @@ namespace user { /* update piston position */ - const real_t gamma_piston = ONE / math::sqrt(ONE - SQR(piston_velocity)); - const real_t v_piston = gamma_piston * piston_velocity; // piston movement over timestep - di_piston += v_piston * dt; + di_piston += piston_velocity * dt; // check if the piston has moved to the next cell i_piston += static_cast(di_piston >= ONE); // keep track of how much the piston has moved into the next cell @@ -210,7 +209,6 @@ namespace user { } // compute current position of piston - // ToDo: check if this can be used as a global variable to avoid recomputing in piston update piston_position = static_cast(i_piston) + di_piston; } @@ -241,10 +239,7 @@ namespace user { template auto CustomParticleUpdate(simtime_t time, spidx_t sp, D& domain) const { - const real_t gamma_piston = ONE / math::sqrt(ONE - SQR(piston_velocity)); - const real_t v_piston = gamma_piston * piston_velocity; - const real_t piston_pos = std::fmod(time * piston_velocity, window_update_frequency); - return CustomPrtlUpdate { piston_pos, v_piston, global_domain.mesh().extent(in::x1).second, true, true}; + return CustomPrtlUpdate { piston_position, piston_velocity, global_xmax, true, true}; }; }; From 366366def4da405c9d046613ee955f790e7ce3bf Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Wed, 8 Apr 2026 22:08:26 +0000 Subject: [PATCH 04/30] bugfix in constructor --- src/archetypes/moving_window.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/archetypes/moving_window.h b/src/archetypes/moving_window.h index eff0c172d..50f78ceec 100644 --- a/src/archetypes/moving_window.h +++ b/src/archetypes/moving_window.h @@ -36,6 +36,7 @@ namespace arch { FieldShift_kernel(ndfield_t Fld, ndfield_t backup_Fld, const int window_shift, BCTags tags) : Fld { Fld } + , backup_Fld { backup_Fld } , window_shift { window_shift } , tags { tags } {} From c01953d7155b9818ce6ae59e6d98d639bc1faec5 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Wed, 8 Apr 2026 22:56:49 +0000 Subject: [PATCH 05/30] added unit conversion for piston_velocity and injection of new plasma --- pgens/piston_shock/pgen.hpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/pgens/piston_shock/pgen.hpp b/pgens/piston_shock/pgen.hpp index 6f3397c1e..d066e8328 100644 --- a/pgens/piston_shock/pgen.hpp +++ b/pgens/piston_shock/pgen.hpp @@ -90,6 +90,7 @@ namespace user { // domain properties const real_t global_xmin, global_xmax; + const real_t cell_size; // gas properties const real_t temperature, temperature_ratio; // injector properties @@ -99,9 +100,9 @@ namespace user { InitFields init_flds; // piston properties - const real_t piston_velocity; + const real_t _piston_velocity; int i_piston; - real_t di_piston, piston_position; + real_t di_piston, piston_position, piston_velocity; // window properties const int window_update_frequency; @@ -111,6 +112,7 @@ namespace user { , global_domain { global_domain } , global_xmin { global_domain.mesh().extent(in::x1).first } , global_xmax { global_domain.mesh().extent(in::x1).second } + , cell_size { (global_xmax - global_xmin) / global_domain.mesh().n_all(in::x1) } , temperature { p.template get("setup.temperature") } , temperature_ratio { p.template get("setup.temperature_ratio") } , Bmag { p.template get("setup.Bmag", ZERO) } @@ -119,8 +121,7 @@ namespace user { , init_flds { Bmag, Btheta, Bphi, ZERO } , dt { p.template get("algorithms.timestep.dt") } , window_update_frequency { p.template get("setup.window_update_frequency", N_GHOSTS) } - , piston_velocity { p.template get("setup.piston_velocity", ZERO) - / math::sqrt(ONE - SQR(p.template get("setup.piston_velocity", ZERO)))} {} + , _piston_velocity { p.template get("setup.piston_velocity", ZERO)} {} inline PGen() {} @@ -134,7 +135,10 @@ namespace user { i_piston = 0; di_piston = ZERO; piston_position = ZERO; - + + coord_t xp_Cd { ZERO }; + // piston velocity in code units + piston_velocity = domain.mesh.metric.template transform<1, Idx::XYZ, Idx::U>(xp_Cd, _piston_velocity); // define temperatures of species const auto temperatures = std::make_pair(temperature, temperature_ratio * temperature); @@ -165,18 +169,14 @@ namespace user { di_piston -= (di_piston >= ONE); // check if the window should be moved - if (i_piston >= window_update_frequency) { - + if ((i_piston >= window_update_frequency) && (window_update_frequency > 0)) { // move the window and all fields and particles in it arch::MoveWindow(domain, window_update_frequency); - // synch ghost zones after moving the window global_domain.CommunicateFields(domain, Comm::E | Comm::B); - /* Inject slab of fresh plasma */ - const real_t cell_size = ZERO; // ToDo: get cell size from global domain const real_t xmax = global_xmax; const real_t xmin = xmax - window_update_frequency * cell_size; // define box to inject into From 931ec61852f3b9f33b9a1c3137d7eca2a303aac1 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Thu, 9 Apr 2026 20:32:17 +0000 Subject: [PATCH 06/30] updated description in header --- src/archetypes/moving_window.h | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/archetypes/moving_window.h b/src/archetypes/moving_window.h index 50f78ceec..49865b26e 100644 --- a/src/archetypes/moving_window.h +++ b/src/archetypes/moving_window.h @@ -1,9 +1,8 @@ /** - * @file archetypes/piston.h - * @brief Piston functions for implementing the piston in the CustomPrtlUpdate in pgran + * @file archetypes/moving_window.h + * @brief Moving window functions for implementing the moving window in the CustomPostStep * @implements - * - arch::PistonUpdate<> -> void - * - arch::CrossesPiston<> -> Bool + * - arch::MoveWindow<> -> void * @namespaces: * - arch:: */ From 96651840ef83469723b119208a49de4c940391b3 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Fri, 10 Apr 2026 01:02:09 +0000 Subject: [PATCH 07/30] communicate particles after shift --- pgens/piston_shock/pgen.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pgens/piston_shock/pgen.hpp b/pgens/piston_shock/pgen.hpp index d066e8328..22f373b41 100644 --- a/pgens/piston_shock/pgen.hpp +++ b/pgens/piston_shock/pgen.hpp @@ -174,6 +174,9 @@ namespace user { arch::MoveWindow(domain, window_update_frequency); // synch ghost zones after moving the window global_domain.CommunicateFields(domain, Comm::E | Comm::B); + // communicate particles after moving + global_domain.CommunicateParticles(domain); + /* Inject slab of fresh plasma */ @@ -209,7 +212,7 @@ namespace user { } // compute current position of piston - piston_position = static_cast(i_piston) + di_piston; + piston_position = static_cast(i_piston) + static_cast(di_piston); } @@ -230,7 +233,6 @@ namespace user { } else { piston_position_use = xg_max; } - if (arch::CrossesPiston(p, pusher, piston_position_use, v_piston, is_left)){ arch::PistonUpdate(p, pusher, piston_position_use, v_piston, massive); } From 654436f036f8f377d165efd349aec2576a3c48dd Mon Sep 17 00:00:00 2001 From: A Sullivan Date: Mon, 20 Apr 2026 09:47:11 -0700 Subject: [PATCH 08/30] fixing piston implementation in piston_shock --- pgens/piston_shock/pgen.hpp | 49 ++++++++++++--- pgens/tutorial/pgen.hpp | 112 +++++++++++++++++++++++++++++++++++ pgens/tutorial/tutorial.toml | 49 +++++++++++++++ 3 files changed, 202 insertions(+), 8 deletions(-) create mode 100644 pgens/tutorial/pgen.hpp create mode 100644 pgens/tutorial/tutorial.toml diff --git a/pgens/piston_shock/pgen.hpp b/pgens/piston_shock/pgen.hpp index 22f373b41..dd39622cc 100644 --- a/pgens/piston_shock/pgen.hpp +++ b/pgens/piston_shock/pgen.hpp @@ -18,6 +18,25 @@ #include #include +/* -------------------------------------------------------------------------- */ +/* Local macros (same as in particle_pusher_sr.hpp) */ +/* -------------------------------------------------------------------------- */ +#define from_Xi_to_i(XI, I) \ + { \ + I = static_cast((XI + 1)) - 1; \ + } + +#define from_Xi_to_i_di(XI, I, DI) \ + { \ + from_Xi_to_i((XI), (I)); \ + DI = static_cast((XI)) - static_cast(I); \ + } + +#define i_di_to_Xi(I, DI) static_cast((I)) + static_cast((DI)) + +/* -------------------------------------------------------------------------- */ + + namespace user { using namespace ntt; @@ -100,13 +119,17 @@ namespace user { InitFields init_flds; // piston properties - const real_t _piston_velocity; + const real_t piston_velocity; int i_piston; - real_t di_piston, piston_position, piston_velocity; + real_t di_piston,piston_velocity_cd, piston_position_cd, piston_position; + // piston properties + const real_t piston_initial_position; // window properties const int window_update_frequency; + + inline PGen(const SimulationParams& p, Metadomain& global_domain) : arch::ProblemGenerator { p } , global_domain { global_domain } @@ -121,7 +144,8 @@ namespace user { , init_flds { Bmag, Btheta, Bphi, ZERO } , dt { p.template get("algorithms.timestep.dt") } , window_update_frequency { p.template get("setup.window_update_frequency", N_GHOSTS) } - , _piston_velocity { p.template get("setup.piston_velocity", ZERO)} {} + , piston_velocity {p.template get("setup.piston_velocity", 0.0)} + , piston_initial_position {p.template get("setup.piston_initial_position", 0.0)}{} inline PGen() {} @@ -132,13 +156,19 @@ namespace user { inline void InitPrtls(Domain& domain) { // set initial position of piston - i_piston = 0; - di_piston = ZERO; - piston_position = ZERO; + + + piston_position = piston_initial_position; + piston_position_cd = domain.mesh.metric.template convert<1, Crd::Ph, Crd::Cd>(piston_position); + + from_Xi_to_i_di(piston_position_cd, i_piston, di_piston); + + + coord_t xp_Cd { ZERO }; // piston velocity in code units - piston_velocity = domain.mesh.metric.template transform<1, Idx::XYZ, Idx::U>(xp_Cd, _piston_velocity); + piston_velocity_cd = domain.mesh.metric.template transform<1, Idx::XYZ, Idx::U>(xp_Cd, piston_velocity); // define temperatures of species const auto temperatures = std::make_pair(temperature, temperature_ratio * temperature); @@ -212,7 +242,10 @@ namespace user { } // compute current position of piston - piston_position = static_cast(i_piston) + static_cast(di_piston); + + piston_position_cd = i_di_to_Xi(i_piston, di_piston); + // convert to physical coordinates + piston_position = domain.mesh.metric.template convert<1, Crd::Cd, Crd::Ph>(piston_position_cd); } diff --git a/pgens/tutorial/pgen.hpp b/pgens/tutorial/pgen.hpp new file mode 100644 index 000000000..f90d154de --- /dev/null +++ b/pgens/tutorial/pgen.hpp @@ -0,0 +1,112 @@ +#ifndef PROBLEM_GENERATOR_H +#define PROBLEM_GENERATOR_H + +#include "enums.h" +#include "global.h" + +#include "archetypes/utils.h" + +#include "archetypes/particle_injector.h" +#include "archetypes/traits.h" +#include "archetypes/problem_generator.h" +#include "framework/domain/metadomain.h" +#include "framework/parameters/parameters.h" + +namespace user { + using namespace ntt; + + //let us define a magnetic field + template + // we define a dipole + struct DipoleField{ + // this is a function for the radius from origin + Inline auto radius(const coord_t& x ) const -> real_t{ + if constexpr (D==Dim::_3D){ + return math::sqrt(SQR(x[0])+SQR(x[1])+SQR(x[2])); + } else { + return math::sqrt(SQR(x[0])+SQR(x[1])); + } + + } + Inline auto bx1(const coord_t& x) const -> real_t{ + return THREE*x[0]*x[1]/math::pow(radius(x), 5); + } + Inline auto bx2(const coord_t& x) const -> real_t{ + const auto r = radius(x); + return (THREE*x[1]*x[1]-SQR(r))/math::pow(r, 5); + } + Inline auto bx3(const coord_t& x) const -> real_t{ + if constexpr (D == Dim::_3D) { + return THREE*x[2]*x[1]/math::pow(radius(x), 5); + } else { + return ZERO; + } + + } + }; + + template + struct PGen : public arch::ProblemGenerator { + static constexpr auto engines { + arch::traits::pgen::compatible_with::value + }; + static constexpr auto metrics { + arch::traits::pgen::compatible_with::value + }; + static constexpr auto dimensions { + arch::traits::pgen::compatible_with::value + }; + // let us now define an init_flds + //DipoleField init_flds; + + //Now let us define this dipole as a sort of external field + inline auto ExternalFields(simtime_t, spidx_t, const Domain&) const + -> std::pair> { + return {true, DipoleField {}}; // first part of this is whether external fields need to be applied + // second return is the field structure itself + // in theory this can be a function of time, or species index + } + + const Metadomain& metadomain; + + inline PGen(const SimulationParams& p, const Metadomain& metadomain) + : arch::ProblemGenerator {p} + , metadomain {metadomain}{} + + inline void InitPrtls(Domain& domain) { + arch::InjectUniformMaxwellian(this->params, + domain, + ONE, // number density in ppc + 1e-3, // this is the temperature + {1u, 2u}, // initialize two particles + { {ONE, ZERO, ZERO}, {ONE, ZERO, ZERO}}, // velocity + false, // do not use particle weights + {{-2.0, -1.5}, Range::All, Range::All}); // adding this last ine designates the strip where there will be plasma + // if this last line were not there, I would inject plasma in the whole domain + //std::map> particles { + // {"x1", {} }, + // {"x2", {} }, + // {"x3", {} }, + // {"ux1", {} }, + // {"ux2", {} }, + // {"ux3", {} } + //}; + //const auto npart = 1000u; + //const auto ymin = metadomain.mesh().extent(in::x2).first; // pick min y + //const auto ymax = metadomain.mesh().extent(in::x2).second; // pick max y + //const auto dy = (ymax - ymin) / static_cast(npart - 1u); // define y spacing + //for (auto p { 0u }; p < npart; ++p) { + // particles["x1"].push_back(-1.5); // we stick them at x1 = -1.5 + // particles["x2"].push_back(ymin + p * dy); // here we are evenly spacing the particles in y + // particles["x3"].push_back(0.0); + // particles["ux1"].push_back(1.0); // define them as moving right + // particles["ux2"].push_back(0.0); + // particles["ux3"].push_back(0.0); + //} + //arch::InjectGlobally(metadomain, domain, 1u, particles); + // third entry in this function is the species + } + }; +} + +#endif \ No newline at end of file diff --git a/pgens/tutorial/tutorial.toml b/pgens/tutorial/tutorial.toml new file mode 100644 index 000000000..eb76dc27f --- /dev/null +++ b/pgens/tutorial/tutorial.toml @@ -0,0 +1,49 @@ +[simulation] + name = "my_first_pgen" + engine = "srpic" + runtime = 1.0 + +[grid] + resolution = [1024, 1024] + extent = [[-2.0, 2.0], [-2.0, 2.0]] + + [grid.metric] + metric = "minkowski" + + [grid.boundaries] + fields = [["PERIODIC"], ["PERIODIC"]] + particles = [["PERIODIC"], ["PERIODIC"]] + +[scales] + larmor0 = 0.003 + skindepth0 = 0.01 + +[algorithms.deposit] + enable = false + +[particles] + ppc0 = 8.0 + + [[particles.species]] + label = "e+" + charge = 1.0 + mass = 1.0 + maxnpart = 1e5 + tracking = true + + +[output] + interval_time = 0.1 + + [output.fields] + quantities = ["N_1", "E", "B"] + + [output.particles] + interval_time = 0.01 + stride = 1 + + [output.spectra] + enable = false + +[checkpoint] + keep = 0 \ No newline at end of file From d72d552d4c89eab29d21b6f324efb60b26de44e9 Mon Sep 17 00:00:00 2001 From: A Sullivan Date: Mon, 20 Apr 2026 14:13:15 -0700 Subject: [PATCH 09/30] fixed bug where not all mpi ranks were updating the domain --- pgens/piston_shock/pgen.hpp | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/pgens/piston_shock/pgen.hpp b/pgens/piston_shock/pgen.hpp index dd39622cc..c3367eca2 100644 --- a/pgens/piston_shock/pgen.hpp +++ b/pgens/piston_shock/pgen.hpp @@ -121,6 +121,7 @@ namespace user { // piston properties const real_t piston_velocity; int i_piston; + int initial_i_piston; real_t di_piston,piston_velocity_cd, piston_position_cd, piston_position; // piston properties const real_t piston_initial_position; @@ -160,8 +161,12 @@ namespace user { piston_position = piston_initial_position; piston_position_cd = domain.mesh.metric.template convert<1, Crd::Ph, Crd::Cd>(piston_position); + + from_Xi_to_i_di(piston_position_cd, i_piston, di_piston); + + initial_i_piston = i_piston; @@ -192,14 +197,17 @@ namespace user { update piston position */ // piston movement over timestep - di_piston += piston_velocity * dt; + + + di_piston += piston_velocity_cd * dt; // check if the piston has moved to the next cell i_piston += static_cast(di_piston >= ONE); // keep track of how much the piston has moved into the next cell di_piston -= (di_piston >= ONE); // check if the window should be moved - if ((i_piston >= window_update_frequency) && (window_update_frequency > 0)) { + if ((i_piston-initial_i_piston >= window_update_frequency) && (window_update_frequency > 0)) { + // move the window and all fields and particles in it arch::MoveWindow(domain, window_update_frequency); // synch ghost zones after moving the window @@ -246,6 +254,7 @@ namespace user { piston_position_cd = i_di_to_Xi(i_piston, di_piston); // convert to physical coordinates piston_position = domain.mesh.metric.template convert<1, Crd::Cd, Crd::Ph>(piston_position_cd); + } From 4de520fcfc449d07e6984f9c0221f2a83991b692 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Mon, 20 Apr 2026 21:18:55 +0000 Subject: [PATCH 10/30] moved communication to window function and simplified pgen --- pgens/piston_shock/pgen.hpp | 60 +++++++++------------------------- src/archetypes/moving_window.h | 10 +++++- 2 files changed, 25 insertions(+), 45 deletions(-) diff --git a/pgens/piston_shock/pgen.hpp b/pgens/piston_shock/pgen.hpp index c3367eca2..ebfcde46c 100644 --- a/pgens/piston_shock/pgen.hpp +++ b/pgens/piston_shock/pgen.hpp @@ -120,11 +120,8 @@ namespace user { // piston properties const real_t piston_velocity; - int i_piston; - int initial_i_piston; - real_t di_piston,piston_velocity_cd, piston_position_cd, piston_position; - // piston properties - const real_t piston_initial_position; + int i_piston, initial_i_piston; + real_t di_piston, piston_velocity_cd, piston_position_cd, piston_position; // window properties const int window_update_frequency; @@ -142,11 +139,10 @@ namespace user { , Bmag { p.template get("setup.Bmag", ZERO) } , Btheta { p.template get("setup.Btheta", ZERO) } , Bphi { p.template get("setup.Bphi", ZERO) } - , init_flds { Bmag, Btheta, Bphi, ZERO } , dt { p.template get("algorithms.timestep.dt") } , window_update_frequency { p.template get("setup.window_update_frequency", N_GHOSTS) } - , piston_velocity {p.template get("setup.piston_velocity", 0.0)} - , piston_initial_position {p.template get("setup.piston_initial_position", 0.0)}{} + , piston_velocity {p.template get("setup.piston_velocity", ZERO)} + , init_flds { Bmag, Btheta, Bphi, ZERO } {} inline PGen() {} @@ -156,10 +152,8 @@ namespace user { inline void InitPrtls(Domain& domain) { - // set initial position of piston - - - piston_position = piston_initial_position; + // set initial position of piston + piston_position = ZERO; piston_position_cd = domain.mesh.metric.template convert<1, Crd::Ph, Crd::Cd>(piston_position); @@ -167,28 +161,21 @@ namespace user { from_Xi_to_i_di(piston_position_cd, i_piston, di_piston); initial_i_piston = i_piston; - - - coord_t xp_Cd { ZERO }; // piston velocity in code units piston_velocity_cd = domain.mesh.metric.template transform<1, Idx::XYZ, Idx::U>(xp_Cd, piston_velocity); + // define temperatures of species const auto temperatures = std::make_pair(temperature, temperature_ratio * temperature); - // define drift speed of species - const auto drifts = std::make_pair(std::vector { ZERO, ZERO, ZERO }, - std::vector { ZERO, ZERO, ZERO }); // inject particles arch::InjectUniformMaxwellians(params, domain, ONE, temperatures, - { 1, 2 }, - drifts, - false); + { 1, 2 }); } void CustomPostStep(timestep_t step, simtime_t time, Domain& domain) { @@ -197,8 +184,6 @@ namespace user { update piston position */ // piston movement over timestep - - di_piston += piston_velocity_cd * dt; // check if the piston has moved to the next cell i_piston += static_cast(di_piston >= ONE); @@ -206,19 +191,15 @@ namespace user { di_piston -= (di_piston >= ONE); // check if the window should be moved - if ((i_piston-initial_i_piston >= window_update_frequency) && (window_update_frequency > 0)) { - + if ((window_update_frequency > 0) && ((i_piston-initial_i_piston) >= window_update_frequency)) { + // move the window and all fields and particles in it - arch::MoveWindow(domain, window_update_frequency); - // synch ghost zones after moving the window - global_domain.CommunicateFields(domain, Comm::E | Comm::B); - // communicate particles after moving - global_domain.CommunicateParticles(domain); - + arch::MoveWindow(domain, global_domain, window_update_frequency); + /* Inject slab of fresh plasma */ - const real_t xmax = global_xmax; + const real_t xmax = global_domain.mesh().extent(in::x1).second; const real_t xmin = xmax - window_update_frequency * cell_size; // define box to inject into boundaries_t inj_box; @@ -250,7 +231,6 @@ namespace user { } // compute current position of piston - piston_position_cd = i_di_to_Xi(i_piston, di_piston); // convert to physical coordinates piston_position = domain.mesh.metric.template convert<1, Crd::Cd, Crd::Ph>(piston_position_cd); @@ -261,29 +241,21 @@ namespace user { struct CustomPrtlUpdate { real_t x_piston; real_t v_piston; - real_t xg_max; bool is_left; bool massive; template Inline void operator()(index_t p, Coord& xp, PusherKernel& pusher) const { - real_t piston_position_use; - - if (x_piston < xg_max){ //make sure piston has not reached the right wall - piston_position_use = x_piston; - } else { - piston_position_use = xg_max; - } - if (arch::CrossesPiston(p, pusher, piston_position_use, v_piston, is_left)){ - arch::PistonUpdate(p, pusher, piston_position_use, v_piston, massive); + if (arch::CrossesPiston(p, pusher, x_piston, v_piston, is_left)){ + arch::PistonUpdate(p, pusher, x_piston, v_piston, massive); } } }; template auto CustomParticleUpdate(simtime_t time, spidx_t sp, D& domain) const { - return CustomPrtlUpdate { piston_position, piston_velocity, global_xmax, true, true}; + return CustomPrtlUpdate { piston_position, piston_velocity, true, true}; }; }; diff --git a/src/archetypes/moving_window.h b/src/archetypes/moving_window.h index 49865b26e..92f595cb1 100644 --- a/src/archetypes/moving_window.h +++ b/src/archetypes/moving_window.h @@ -139,7 +139,7 @@ namespace arch { */ template - Inline void MoveWindow(Domain& domain, + Inline void MoveWindow(Domain& domain, Metadomain& global_domain, const int window_shift) { if constexpr (M::CoordType != Coord::Cart) { @@ -281,6 +281,14 @@ namespace arch { raise::Error("Invalid dimension", HERE); } } + + // synch ghost zones after moving the window + global_domain.CommunicateFields(domain, Comm::E | Comm::B); + // communicate particles after moving + global_domain.CommunicateParticles(domain); + + // ToDo: Update metric + } } } // namespace arch From 52b46a9cb400d26988b1cf139d41736f8d02f848 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Tue, 21 Apr 2026 02:33:34 +0000 Subject: [PATCH 11/30] always shift window by N_GHOSTS --- pgens/piston_shock/pgen.hpp | 26 ++++++++++---------------- 1 file changed, 10 insertions(+), 16 deletions(-) diff --git a/pgens/piston_shock/pgen.hpp b/pgens/piston_shock/pgen.hpp index ebfcde46c..cd754c7f6 100644 --- a/pgens/piston_shock/pgen.hpp +++ b/pgens/piston_shock/pgen.hpp @@ -123,11 +123,7 @@ namespace user { int i_piston, initial_i_piston; real_t di_piston, piston_velocity_cd, piston_position_cd, piston_position; - // window properties - const int window_update_frequency; - - - + inline PGen(const SimulationParams& p, Metadomain& global_domain) : arch::ProblemGenerator { p } , global_domain { global_domain } @@ -136,12 +132,11 @@ namespace user { , cell_size { (global_xmax - global_xmin) / global_domain.mesh().n_all(in::x1) } , temperature { p.template get("setup.temperature") } , temperature_ratio { p.template get("setup.temperature_ratio") } + , piston_velocity {p.template get("setup.piston_velocity", ZERO)} , Bmag { p.template get("setup.Bmag", ZERO) } , Btheta { p.template get("setup.Btheta", ZERO) } , Bphi { p.template get("setup.Bphi", ZERO) } , dt { p.template get("algorithms.timestep.dt") } - , window_update_frequency { p.template get("setup.window_update_frequency", N_GHOSTS) } - , piston_velocity {p.template get("setup.piston_velocity", ZERO)} , init_flds { Bmag, Btheta, Bphi, ZERO } {} inline PGen() {} @@ -155,15 +150,13 @@ namespace user { // set initial position of piston piston_position = ZERO; piston_position_cd = domain.mesh.metric.template convert<1, Crd::Ph, Crd::Cd>(piston_position); - - - + // convert to cell units and get cell index and sub-cell position from_Xi_to_i_di(piston_position_cd, i_piston, di_piston); - + // store initial cell index of piston for window updates initial_i_piston = i_piston; - coord_t xp_Cd { ZERO }; // piston velocity in code units + coord_t xp_Cd { ZERO }; piston_velocity_cd = domain.mesh.metric.template transform<1, Idx::XYZ, Idx::U>(xp_Cd, piston_velocity); // define temperatures of species @@ -191,16 +184,16 @@ namespace user { di_piston -= (di_piston >= ONE); // check if the window should be moved - if ((window_update_frequency > 0) && ((i_piston-initial_i_piston) >= window_update_frequency)) { + if ((i_piston-initial_i_piston) >= N_GHOSTS) { // move the window and all fields and particles in it - arch::MoveWindow(domain, global_domain, window_update_frequency); + arch::MoveWindow(domain, global_domain, N_GHOSTS); /* Inject slab of fresh plasma */ const real_t xmax = global_domain.mesh().extent(in::x1).second; - const real_t xmin = xmax - window_update_frequency * cell_size; + const real_t xmin = xmax - N_GHOSTS * cell_size; // define box to inject into boundaries_t inj_box; // loop over all dimension @@ -227,7 +220,8 @@ namespace user { false, inj_box); - i_piston -= window_update_frequency; + // shift the piston indices back by the number of ghost cells to account for window movement + i_piston -= N_GHOSTS; } // compute current position of piston From d3b2641531e3797ef492ed2d514dc85c5422e0b5 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Tue, 21 Apr 2026 02:35:46 +0000 Subject: [PATCH 12/30] applied formatting --- pgens/piston_shock/pgen.hpp | 79 +++++++++++++++--------------- src/archetypes/moving_window.h | 87 +++++++++++++++------------------- 2 files changed, 77 insertions(+), 89 deletions(-) diff --git a/pgens/piston_shock/pgen.hpp b/pgens/piston_shock/pgen.hpp index cd754c7f6..99b430c7a 100644 --- a/pgens/piston_shock/pgen.hpp +++ b/pgens/piston_shock/pgen.hpp @@ -8,11 +8,11 @@ #include "utils/numeric.h" #include "archetypes/field_setter.h" +#include "archetypes/moving_window.h" +#include "archetypes/piston.h" #include "archetypes/problem_generator.h" #include "archetypes/traits.h" #include "archetypes/utils.h" -#include "archetypes/piston.h" -#include "archetypes/moving_window.h" #include "framework/domain/metadomain.h" #include @@ -36,7 +36,6 @@ /* -------------------------------------------------------------------------- */ - namespace user { using namespace ntt; @@ -120,19 +119,19 @@ namespace user { // piston properties const real_t piston_velocity; - int i_piston, initial_i_piston; + int i_piston, initial_i_piston; real_t di_piston, piston_velocity_cd, piston_position_cd, piston_position; - inline PGen(const SimulationParams& p, Metadomain& global_domain) : arch::ProblemGenerator { p } , global_domain { global_domain } , global_xmin { global_domain.mesh().extent(in::x1).first } , global_xmax { global_domain.mesh().extent(in::x1).second } - , cell_size { (global_xmax - global_xmin) / global_domain.mesh().n_all(in::x1) } + , cell_size { (global_xmax - global_xmin) / + global_domain.mesh().n_all(in::x1) } , temperature { p.template get("setup.temperature") } , temperature_ratio { p.template get("setup.temperature_ratio") } - , piston_velocity {p.template get("setup.piston_velocity", ZERO)} + , piston_velocity { p.template get("setup.piston_velocity", ZERO) } , Bmag { p.template get("setup.Bmag", ZERO) } , Btheta { p.template get("setup.Btheta", ZERO) } , Bphi { p.template get("setup.Bphi", ZERO) } @@ -147,9 +146,10 @@ namespace user { inline void InitPrtls(Domain& domain) { - // set initial position of piston + // set initial position of piston piston_position = ZERO; - piston_position_cd = domain.mesh.metric.template convert<1, Crd::Ph, Crd::Cd>(piston_position); + piston_position_cd = domain.mesh.metric.template convert<1, Crd::Ph, Crd::Cd>( + piston_position); // convert to cell units and get cell index and sub-cell position from_Xi_to_i_di(piston_position_cd, i_piston, di_piston); // store initial cell index of piston for window updates @@ -157,43 +157,41 @@ namespace user { // piston velocity in code units coord_t xp_Cd { ZERO }; - piston_velocity_cd = domain.mesh.metric.template transform<1, Idx::XYZ, Idx::U>(xp_Cd, piston_velocity); + piston_velocity_cd = domain.mesh.metric.template transform<1, Idx::XYZ, Idx::U>( + xp_Cd, + piston_velocity); // define temperatures of species const auto temperatures = std::make_pair(temperature, temperature_ratio * temperature); // inject particles - arch::InjectUniformMaxwellians(params, - domain, - ONE, - temperatures, - { 1, 2 }); + arch::InjectUniformMaxwellians(params, domain, ONE, temperatures, { 1, 2 }); } void CustomPostStep(timestep_t step, simtime_t time, Domain& domain) { - /* + /* update piston position */ // piston movement over timestep di_piston += piston_velocity_cd * dt; // check if the piston has moved to the next cell - i_piston += static_cast(di_piston >= ONE); + i_piston += static_cast(di_piston >= ONE); // keep track of how much the piston has moved into the next cell di_piston -= (di_piston >= ONE); // check if the window should be moved - if ((i_piston-initial_i_piston) >= N_GHOSTS) { - + if ((i_piston - initial_i_piston) >= N_GHOSTS) { + // move the window and all fields and particles in it arch::MoveWindow(domain, global_domain, N_GHOSTS); - + /* Inject slab of fresh plasma */ - const real_t xmax = global_domain.mesh().extent(in::x1).second; - const real_t xmin = xmax - N_GHOSTS * cell_size; + const real_t xmax = global_domain.mesh().extent(in::x1).second; + const real_t xmin = xmax - N_GHOSTS * cell_size; // define box to inject into boundaries_t inj_box; // loop over all dimension @@ -207,18 +205,18 @@ namespace user { // same maxwell distribution as above const auto temperatures = std::make_pair(temperature, - temperature_ratio * temperature); - const auto drifts = std::make_pair( + temperature_ratio * temperature); + const auto drifts = std::make_pair( std::vector { ZERO, ZERO, ZERO }, std::vector { ZERO, ZERO, ZERO }); arch::InjectUniformMaxwellians(params, - domain, - ONE, - temperatures, - { 1, 2 }, - drifts, - false, - inj_box); + domain, + ONE, + temperatures, + { 1, 2 }, + drifts, + false, + inj_box); // shift the piston indices back by the number of ghost cells to account for window movement i_piston -= N_GHOSTS; @@ -227,30 +225,29 @@ namespace user { // compute current position of piston piston_position_cd = i_di_to_Xi(i_piston, di_piston); // convert to physical coordinates - piston_position = domain.mesh.metric.template convert<1, Crd::Cd, Crd::Ph>(piston_position_cd); - + piston_position = domain.mesh.metric.template convert<1, Crd::Cd, Crd::Ph>( + piston_position_cd); } - struct CustomPrtlUpdate { real_t x_piston; real_t v_piston; - bool is_left; - bool massive; - + bool is_left; + bool massive; + template Inline void operator()(index_t p, Coord& xp, PusherKernel& pusher) const { - if (arch::CrossesPiston(p, pusher, x_piston, v_piston, is_left)){ - arch::PistonUpdate(p, pusher, x_piston, v_piston, massive); + if (arch::CrossesPiston(p, pusher, x_piston, v_piston, is_left)) { + arch::PistonUpdate(p, pusher, x_piston, v_piston, massive); } } }; template auto CustomParticleUpdate(simtime_t time, spidx_t sp, D& domain) const { - return CustomPrtlUpdate { piston_position, piston_velocity, true, true}; - }; + return CustomPrtlUpdate { piston_position, piston_velocity, true, true }; + }; }; } // namespace user diff --git a/src/archetypes/moving_window.h b/src/archetypes/moving_window.h index 92f595cb1..f8cd51ead 100644 --- a/src/archetypes/moving_window.h +++ b/src/archetypes/moving_window.h @@ -15,12 +15,11 @@ #include "archetypes/energy_dist.h" #include "framework/domain/domain.h" -#include "framework/parameters/parameters.h" #include "framework/domain/metadomain.h" +#include "framework/parameters/parameters.h" #include - namespace arch { template @@ -28,12 +27,15 @@ namespace arch { static_assert(static_cast(o) < static_cast(D), "Invalid component index"); - ndfield_t Fld; - ndfield_t backup_Fld; - const BCTags tags; - const int window_shift; + ndfield_t Fld; + ndfield_t backup_Fld; + const BCTags tags; + const int window_shift; - FieldShift_kernel(ndfield_t Fld, ndfield_t backup_Fld, const int window_shift, BCTags tags) + FieldShift_kernel(ndfield_t Fld, + ndfield_t backup_Fld, + const int window_shift, + BCTags tags) : Fld { Fld } , backup_Fld { backup_Fld } , window_shift { window_shift } @@ -139,20 +141,19 @@ namespace arch { */ template - Inline void MoveWindow(Domain& domain, Metadomain& global_domain, - const int window_shift) { - + Inline void MoveWindow(Domain& domain, + Metadomain& global_domain, + const int window_shift) { + if constexpr (M::CoordType != Coord::Cart) { - raise::Error( - "Moving window only applicable to cartesian coordinates", - HERE); + raise::Error("Moving window only applicable to cartesian coordinates", HERE); } else { /* move particles in the window back by the window size */ - const auto nspec = domain.species.size(); - const auto& mesh = domain.mesh; + const auto nspec = domain.species.size(); + const auto& mesh = domain.mesh; if (o == in::x1) { // loop over particle species for (auto s { 0u }; s < nspec; ++s) { @@ -165,8 +166,7 @@ namespace arch { Lambda(index_t p) { // shift particle position back by window update frequency i1(p) -= window_shift; - } - ); + }); } } else if (o == in::x2) { if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { @@ -181,8 +181,7 @@ namespace arch { Lambda(index_t p) { // shift particle position back by window update frequency i2(p) -= window_shift; - } - ); + }); } } else { raise::Error("Invalid dimension", HERE); @@ -200,8 +199,7 @@ namespace arch { Lambda(index_t p) { // shift particle position back by window update frequency i3(p) -= window_shift; - } - ); + }); } } else { raise::Error("Invalid direction", HERE); @@ -212,7 +210,7 @@ namespace arch { // shift fields in the window back by the window size std::vector xi_min, xi_max; - const std::vector all_dirs { in::x1, in::x2, in::x3 }; + 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 (o == dd) { @@ -243,40 +241,34 @@ namespace arch { } // copy fields to backup before shifting - Kokkos::deep_copy(domain.fields.bckp, domain.fields.em); + Kokkos::deep_copy(domain.fields.bckp, domain.fields.em); if (o == in::x1) { - Kokkos::parallel_for( - "ShiftFields", - range, - FieldShift_kernel( - domain.fields.em, - domain.fields.bckp, - window_shift, - BC::B | BC::E)); + Kokkos::parallel_for("ShiftFields", + range, + FieldShift_kernel(domain.fields.em, + domain.fields.bckp, + window_shift, + BC::B | BC::E)); } else if (o == in::x2) { if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - Kokkos::parallel_for( - "ShiftFields", - range, - FieldShift_kernel( - domain.fields.em, - domain.fields.bckp, - window_shift, - BC::B | BC::E)); + Kokkos::parallel_for("ShiftFields", + range, + FieldShift_kernel(domain.fields.em, + domain.fields.bckp, + window_shift, + BC::B | BC::E)); } else { raise::Error("Invalid dimension", HERE); } } else { if constexpr (M::Dim == Dim::_3D) { - Kokkos::parallel_for( - "ShiftFields", - range, - FieldShift_kernel( - domain.fields.em, - domain.fields.bckp, - window_shift, - BC::B | BC::E)); + Kokkos::parallel_for("ShiftFields", + range, + FieldShift_kernel(domain.fields.em, + domain.fields.bckp, + window_shift, + BC::B | BC::E)); } else { raise::Error("Invalid dimension", HERE); } @@ -288,7 +280,6 @@ namespace arch { global_domain.CommunicateParticles(domain); // ToDo: Update metric - } } } // namespace arch From 45c6650b5fd778a03244cc855d3e690625a3e50b Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Tue, 21 Apr 2026 02:40:31 +0000 Subject: [PATCH 13/30] removed redundant field initialization --- pgens/piston_shock/pgen.hpp | 22 +++++----------------- 1 file changed, 5 insertions(+), 17 deletions(-) diff --git a/pgens/piston_shock/pgen.hpp b/pgens/piston_shock/pgen.hpp index 99b430c7a..14ac92c01 100644 --- a/pgens/piston_shock/pgen.hpp +++ b/pgens/piston_shock/pgen.hpp @@ -50,11 +50,10 @@ namespace user { @param bphi: magnetic field azimuthal angle @param drift_ux: drift velocity in the x direction */ - InitFields(real_t bmag, real_t btheta, real_t bphi, real_t drift_ux) + InitFields(real_t bmag, real_t btheta, real_t bphi) : Bmag { bmag } , Btheta { btheta * static_cast(convert::deg2rad) } - , Bphi { bphi * static_cast(convert::deg2rad) } - , Vx { drift_ux } {} + , Bphi { bphi * static_cast(convert::deg2rad) } {} // magnetic field components Inline auto bx1(const coord_t&) const -> real_t { @@ -69,21 +68,10 @@ namespace user { return Bmag * math::sin(Btheta) * math::cos(Bphi); } - // electric field components - Inline auto ex1(const coord_t&) const -> real_t { - return ZERO; - } - - Inline auto ex2(const coord_t&) const -> real_t { - return -Vx * Bmag * math::sin(Btheta) * math::cos(Bphi); - } - - Inline auto ex3(const coord_t&) const -> real_t { - return Vx * Bmag * math::sin(Btheta) * math::sin(Bphi); - } + // electric field components are all zero private: - const real_t Btheta, Bphi, Vx, Bmag; + const real_t Bmag, Btheta, Bphi; }; template @@ -136,7 +124,7 @@ namespace user { , Btheta { p.template get("setup.Btheta", ZERO) } , Bphi { p.template get("setup.Bphi", ZERO) } , dt { p.template get("algorithms.timestep.dt") } - , init_flds { Bmag, Btheta, Bphi, ZERO } {} + , init_flds { Bmag, Btheta, Bphi } {} inline PGen() {} From dfd0fa82563deb216b4879b009028e291528b716 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Tue, 21 Apr 2026 02:50:31 +0000 Subject: [PATCH 14/30] removed redundant variables --- pgens/piston_shock/pgen.hpp | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/pgens/piston_shock/pgen.hpp b/pgens/piston_shock/pgen.hpp index 14ac92c01..3a0c6ced2 100644 --- a/pgens/piston_shock/pgen.hpp +++ b/pgens/piston_shock/pgen.hpp @@ -43,12 +43,10 @@ namespace user { struct InitFields { /* Sets up magnetic and electric field components for the simulation. - Must satisfy E = -v x B for Lorentz Force to be zero. @param bmag: magnetic field scaling @param btheta: magnetic field polar angle @param bphi: magnetic field azimuthal angle - @param drift_ux: drift velocity in the x direction */ InitFields(real_t bmag, real_t btheta, real_t bphi) : Bmag { bmag } @@ -95,7 +93,6 @@ namespace user { Metadomain& global_domain; // domain properties - const real_t global_xmin, global_xmax; const real_t cell_size; // gas properties const real_t temperature, temperature_ratio; @@ -113,9 +110,8 @@ namespace user { inline PGen(const SimulationParams& p, Metadomain& global_domain) : arch::ProblemGenerator { p } , global_domain { global_domain } - , global_xmin { global_domain.mesh().extent(in::x1).first } - , global_xmax { global_domain.mesh().extent(in::x1).second } - , cell_size { (global_xmax - global_xmin) / + , cell_size { (global_domain.mesh().extent(in::x1).second - + global_domain.mesh().extent(in::x1).first) / global_domain.mesh().n_all(in::x1) } , temperature { p.template get("setup.temperature") } , temperature_ratio { p.template get("setup.temperature_ratio") } From acad86348e3fd64656f304dcfa7443789d4bb6af Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Wed, 22 Apr 2026 17:14:20 +0000 Subject: [PATCH 15/30] updated moving window to match new class definitions --- src/archetypes/moving_window.h | 209 ++++++++++++++++----------------- 1 file changed, 102 insertions(+), 107 deletions(-) diff --git a/src/archetypes/moving_window.h b/src/archetypes/moving_window.h index f8cd51ead..c0967df3d 100644 --- a/src/archetypes/moving_window.h +++ b/src/archetypes/moving_window.h @@ -140,147 +140,142 @@ namespace arch { * @brief Updates particle position and fields in the moving window. */ - template + template Inline void MoveWindow(Domain& domain, Metadomain& global_domain, const int window_shift) { - if constexpr (M::CoordType != Coord::Cart) { - raise::Error("Moving window only applicable to cartesian coordinates", HERE); - } else { - - /* - move particles in the window back by the window size - */ - const auto nspec = domain.species.size(); - const auto& mesh = domain.mesh; - if (o == in::x1) { + /* + move particles in the window back by the window size + */ + const auto nspec = domain.species.size(); + const auto& mesh = domain.mesh; + if (o == in::x1) { + // loop over particle species + for (auto s { 0u }; s < nspec; ++s) { + // get particle properties + auto& species = domain.species[s]; + auto i1 = species.i1; + Kokkos::parallel_for( + "MoveParticles", + species.rangeActiveParticles(), + Lambda(index_t p) { + // shift particle position back by window update frequency + i1(p) -= window_shift; + }); + } + } else if (o == in::x2) { + if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { // loop over particle species for (auto s { 0u }; s < nspec; ++s) { // get particle properties auto& species = domain.species[s]; - auto i1 = species.i1; + auto i2 = species.i2; Kokkos::parallel_for( "MoveParticles", species.rangeActiveParticles(), Lambda(index_t p) { // shift particle position back by window update frequency - i1(p) -= window_shift; + i2(p) -= window_shift; }); } - } else if (o == in::x2) { - if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - // loop over particle species - for (auto s { 0u }; s < nspec; ++s) { - // get particle properties - auto& species = domain.species[s]; - auto i2 = species.i2; - Kokkos::parallel_for( - "MoveParticles", - species.rangeActiveParticles(), - Lambda(index_t p) { - // shift particle position back by window update frequency - i2(p) -= window_shift; - }); - } - } else { - raise::Error("Invalid dimension", HERE); - } - } else if (o == in::x3) { - if constexpr (M::Dim == Dim::_3D) { - // loop over particle species - for (auto s { 0u }; s < nspec; ++s) { - // get particle properties - auto& species = domain.species[s]; - auto i3 = species.i3; - Kokkos::parallel_for( - "MoveParticles", - species.rangeActiveParticles(), - Lambda(index_t p) { - // shift particle position back by window update frequency - i3(p) -= window_shift; - }); - } - } else { - raise::Error("Invalid direction", HERE); + } else { + raise::Error("Invalid dimension", HERE); + } + } else if (o == in::x3) { + if constexpr (M::Dim == Dim::_3D) { + // loop over particle species + for (auto s { 0u }; s < nspec; ++s) { + // get particle properties + auto& species = domain.species[s]; + auto i3 = species.i3; + Kokkos::parallel_for( + "MoveParticles", + species.rangeActiveParticles(), + Lambda(index_t p) { + // shift particle position back by window update frequency + i3(p) -= window_shift; + }); } } else { raise::Error("Invalid direction", HERE); } + } else { + raise::Error("Invalid direction", HERE); + } - // shift fields in the window back by the window size - 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 (o == dd) { - xi_min.push_back(0); - xi_max.push_back(domain.mesh.n_all(dd) - window_shift); - } else { - xi_min.push_back(0); - xi_max.push_back(domain.mesh.n_all(dd)); - } + // shift fields in the window back by the window size + 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 (o == dd) { + xi_min.push_back(0); + xi_max.push_back(domain.mesh.n_all(dd) - window_shift); + } 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); + } + raise::ErrorIf(xi_min.size() != xi_max.size() or + xi_min.size() != static_cast(M::Dim), + "Invalid range size", + HERE); + + // loop range for shifting fields + 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); + } + + // copy fields to backup before shifting + Kokkos::deep_copy(domain.fields.bckp, domain.fields.em); - // loop range for shifting fields - 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] }); + if (o == in::x1) { + Kokkos::parallel_for("ShiftFields", + range, + FieldShift_kernel(domain.fields.em, + domain.fields.bckp, + window_shift, + BC::B | BC::E)); + } else if (o == in::x2) { + if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { + Kokkos::parallel_for("ShiftFields", + range, + FieldShift_kernel(domain.fields.em, + domain.fields.bckp, + window_shift, + BC::B | BC::E)); } else { raise::Error("Invalid dimension", HERE); } - - // copy fields to backup before shifting - Kokkos::deep_copy(domain.fields.bckp, domain.fields.em); - - if (o == in::x1) { + } else { + if constexpr (M::Dim == Dim::_3D) { Kokkos::parallel_for("ShiftFields", range, - FieldShift_kernel(domain.fields.em, + FieldShift_kernel(domain.fields.em, domain.fields.bckp, window_shift, BC::B | BC::E)); - } else if (o == in::x2) { - if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - Kokkos::parallel_for("ShiftFields", - range, - FieldShift_kernel(domain.fields.em, - domain.fields.bckp, - window_shift, - BC::B | BC::E)); - } else { - raise::Error("Invalid dimension", HERE); - } } else { - if constexpr (M::Dim == Dim::_3D) { - Kokkos::parallel_for("ShiftFields", - range, - FieldShift_kernel(domain.fields.em, - domain.fields.bckp, - window_shift, - BC::B | BC::E)); - } else { - raise::Error("Invalid dimension", HERE); - } + raise::Error("Invalid dimension", HERE); } + } - // synch ghost zones after moving the window - global_domain.CommunicateFields(domain, Comm::E | Comm::B); - // communicate particles after moving - global_domain.CommunicateParticles(domain); + // synch ghost zones after moving the window + global_domain.CommunicateFields(domain, Comm::E | Comm::B); + // communicate particles after moving + global_domain.CommunicateParticles(domain); - // ToDo: Update metric - } + // ToDo: Update metric } } // namespace arch From c8fa56336f49bd69ce50655aa8855fe9a089a3a0 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Wed, 22 Apr 2026 18:09:26 +0000 Subject: [PATCH 16/30] moved check for piston crossing inside of piston --- src/archetypes/piston.h | 110 ++++++++++++++++++++++------------------ 1 file changed, 61 insertions(+), 49 deletions(-) diff --git a/src/archetypes/piston.h b/src/archetypes/piston.h index 09bc99d74..20e6f03a4 100644 --- a/src/archetypes/piston.h +++ b/src/archetypes/piston.h @@ -37,6 +37,50 @@ namespace arch { + /** + * @brief Checks whether a particle reflects off a moving piston, called after particle has been moved by regular pusher + * + * @param p Index of particle + * @param dt Timestep + * @param particles Particle data arrays + * @param metric Metric object for coordinate transformations + * @param piston_position Position of the piston at the start of timestep in global coordinates + * @param piston_v Velocity of piston at current timestep + * @param is_left Is piston on the left side of the box or right side of the box + */ + template + Inline bool CrossesPiston(index_t p, + real_t dt, + const kernel::PusherArrays& particles, + const M& metric, + real_t piston_position, + real_t piston_v, + bool is_left) { + const real_t x1_Cd = i_di_to_Xi(particles.i1(p), particles.dx1(p)); + // x1_Cd_wallmove is not the actual particle coordinate + // it is particle position minus how much the wall has moved in this + // timestep This is a computational trick + const real_t x1_Cd_wallmove = + x1_Cd - metric.template transform<1, Idx::XYZ, Idx::U>({}, piston_v) * dt; + const real_t x1_Ph_wallmove = metric.template convert<1, Crd::Cd, Crd::Ph>( + x1_Cd_wallmove); + + if (is_left) { // if piston is moving from left, ask if particle is to the left of piston + if (piston_position > x1_Ph_wallmove) { + return true; + } else { + return false; + } + } else { // if piston is moving from the right, so ask is particle to right of piston + if (piston_position < x1_Ph_wallmove) { + return true; + } else { + return false; + } + } + } + + /** * @brief Updates particle position and velocity if it reflects off a moiving * piston, called to correct patticle position @@ -51,12 +95,23 @@ namespace arch { */ template Inline void PistonUpdate(index_t p, - real_t dt, - const kernel::PusherArrays& particles, - const M& metric, - real_t piston_position, - real_t piston_v, - bool massive) { + real_t dt, + const kernel::PusherArrays& particles, + const M& metric, + real_t piston_position, + real_t piston_v, + bool massive) { + + // check if particle actually crosses the piston, if not return + if (!CrossesPiston(p, + dt, + particles, + metric, + piston_position, + piston_v, + true)) { + return; + } // step 1: calculate the particle 3 velocity const real_t gamma_p { massive ? U2GAMMA(particles.ux1(p), particles.ux2(p), particles.ux3(p)) @@ -125,49 +180,6 @@ namespace arch { particles.dx1(p) += (particles.dx1(p) < ZERO); } - /** - * @brief Checks whether a particle reflects off a moving piston, called after particle has been moved by regular pusher - * - * @param p Index of particle - * @param dt Timestep - * @param particles Particle data arrays - * @param metric Metric object for coordinate transformations - * @param piston_position Position of the piston at the start of timestep in global coordinates - * @param piston_v Velocity of piston at current timestep - * @param is_left Is piston on the left side of the box or right side of the box - */ - template - Inline bool CrossesPiston(index_t p, - real_t dt, - const kernel::PusherArrays& particles, - const M& metric, - real_t piston_position, - real_t piston_v, - bool is_left) { - const real_t x1_Cd = i_di_to_Xi(particles.i1(p), particles.dx1(p)); - // x1_Cd_wallmove is not the actual particle coordinate - // it is particle position minus how much the wall has moved in this - // timestep This is a computational trick - const real_t x1_Cd_wallmove = - x1_Cd - metric.template transform<1, Idx::XYZ, Idx::U>({}, piston_v) * dt; - const real_t x1_Ph_wallmove = metric.template convert<1, Crd::Cd, Crd::Ph>( - x1_Cd_wallmove); - - if (is_left) { // if piston is moving from left, ask if particle is to the left of piston - if (piston_position > x1_Ph_wallmove) { - return true; - } else { - return false; - } - } else { // if piston is moving from the right, so ask is particle to right of piston - if (piston_position < x1_Ph_wallmove) { - return true; - } else { - return false; - } - } - } - } // namespace arch #undef from_Xi_to_i_di From d5f7aeb1803d873f135d7cc449234bc440b00fbf Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Wed, 22 Apr 2026 21:53:42 +0000 Subject: [PATCH 17/30] renamed `PistonUpdate` to `Piston` --- src/archetypes/piston.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/archetypes/piston.h b/src/archetypes/piston.h index 20e6f03a4..e7e2b2f92 100644 --- a/src/archetypes/piston.h +++ b/src/archetypes/piston.h @@ -82,7 +82,7 @@ namespace arch { /** - * @brief Updates particle position and velocity if it reflects off a moiving + * @brief Updates particle position and velocity if it reflects off a moving * piston, called to correct patticle position * * @param p Index of particle @@ -94,7 +94,7 @@ namespace arch { * @param massive Whether the particle is massive or massless (e.g. photon) */ template - Inline void PistonUpdate(index_t p, + Inline void Piston(index_t p, real_t dt, const kernel::PusherArrays& particles, const M& metric, From 1f51356a1074112d9aea2ffdb9b1f3747260d1fd Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Wed, 22 Apr 2026 21:54:03 +0000 Subject: [PATCH 18/30] example for piston --- pgens/examples/piston/pgen.hpp | 159 ++++++++++++++++++++++++++++++ pgens/examples/piston/piston.py | 25 +++++ pgens/examples/piston/piston.toml | 57 +++++++++++ 3 files changed, 241 insertions(+) create mode 100644 pgens/examples/piston/pgen.hpp create mode 100644 pgens/examples/piston/piston.py create mode 100644 pgens/examples/piston/piston.toml diff --git a/pgens/examples/piston/pgen.hpp b/pgens/examples/piston/pgen.hpp new file mode 100644 index 000000000..1908b1514 --- /dev/null +++ b/pgens/examples/piston/pgen.hpp @@ -0,0 +1,159 @@ +#ifndef PROBLEM_GENERATOR_H +#define PROBLEM_GENERATOR_H + +#include "enums.h" +#include "global.h" + +#include "traits/pgen.h" + +#include "archetypes/energy_dist.h" +#include "archetypes/particle_injector.h" +#include "archetypes/piston.h" +#include "archetypes/problem_generator.h" +#include "archetypes/spatial_dist.h" +#include "framework/domain/metadomain.h" +#include "kernels/particle_moments.hpp" + +namespace user { + using namespace ntt; + + template + struct NonUniformTargetDensity { + Inline auto operator()(const coord_t& x_Ph) const -> real_t { + // example of a non-uniform target density that peaks at the center of the domain + real_t r2 { ZERO }; + for (auto d = 0u; d < D; ++d) { + r2 += SQR(x_Ph[d] - 0.5); + } + return std::exp( + -r2 / SQR(0.2)); // <-- characteristic width of the density profile + } + }; + + template + struct PGen : public arch::ProblemGenerator { + + static constexpr auto engines { + ::traits::pgen::compatible_with {} + }; + static constexpr auto metrics { + ::traits::pgen::compatible_with {} + }; + static constexpr auto dimensions { + ::traits::pgen::compatible_with {} + }; + + using arch::ProblemGenerator::D; + using arch::ProblemGenerator::C; + using arch::ProblemGenerator::params; + + Metadomain& global_domain; + const real_t piston_velocity; + + inline PGen(const SimulationParams& p, Metadomain& global_domain) + : arch::ProblemGenerator { p } + , global_domain { global_domain } + , piston_velocity { p.template get("setup.piston_velocity", 0.0) } {} + + inline void InitPrtls(Domain& domain) { + + { // compute density of species #1 and #2 + + // saves the density to domain.fields.buff(:,:,:,0) + const auto ni2 = domain.mesh.n_active(in::x2); + const auto inv_n0 = ONE / params.template get("scales.n0"); + + auto scatter_buff = Kokkos::Experimental::create_scatter_view( + domain.fields.buff); + Kokkos::deep_copy(domain.fields.buff, ZERO); + for (const auto sp : std::vector { 1, 2 }) { + const auto& prtl_spec = domain.species[sp - 1]; + Kokkos::parallel_for("ComputeDensity", + prtl_spec.rangeActiveParticles(), + kernel::ParticleMoments_kernel( + {}, + scatter_buff, + 0u, + 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(), + false, + domain.mesh.metric, + domain.mesh.flds_bc(), + ni2, + inv_n0, + 0u)); + } + Kokkos::Experimental::contribute(domain.fields.buff, scatter_buff); + } + // inject particles + const auto energy_dist = arch::Maxwellian( + domain.mesh.metric, + domain.random_pool(), + 0.2); // <-- target temperature for injection + + const auto target_density_profile = NonUniformTargetDensity {}; + const auto replenish_sdist = + arch::Replenish( + domain.mesh.metric, + domain.fields.buff, + 0u, // <-- index in buff where the density is stored + target_density_profile, + ONE); // <-- target density for replenishment + arch::InjectNonUniform( + params, + domain, + { 1, 2 }, + { energy_dist, energy_dist }, + replenish_sdist, + ONE); + } + + struct CustomPrtlUpdate { + real_t x_piston; // current position of piston + real_t v_piston; + real_t global_xmax; + bool is_left; + bool massive; + + Inline void operator()(index_t p, + const kernel::PusherContext& ctx, + const kernel::PusherBoundaries&, + const kernel::PusherArrays& particles, + const M& metric) const { + + real_t piston_pos; + if (x_piston > global_xmax) { + // piston has moved beyond the domain, set position to a large value so that no particles interact with it + piston_pos = global_xmax; + } else { + piston_pos = x_piston; + } + + arch::Piston(p, ctx.dt, particles, metric, piston_pos, v_piston, massive); + } + }; + + template + auto CustomParticleUpdate(simtime_t time, spidx_t sp, DOM& domain) const { + return CustomPrtlUpdate { piston_velocity * static_cast(time), + piston_velocity, + global_domain.mesh().extent(in::x1).second, + true, + true }; + }; + }; +} // namespace user + +#endif diff --git a/pgens/examples/piston/piston.py b/pgens/examples/piston/piston.py new file mode 100644 index 000000000..2f8c03152 --- /dev/null +++ b/pgens/examples/piston/piston.py @@ -0,0 +1,25 @@ +import nt2 +import matplotlib.pyplot as plt + + +def plot(t, data): + plt.rcParams["text.usetex"] = True + plt.rcParams["font.family"] = "serif" + plt.rcParams["figure.dpi"] = 250 + ax = plt.gca() + data.fields.N_1_2.sel(t=t, method="nearest").plot( + vmin=0, vmax=1.5, ax=ax, cmap="inferno" + ) + plt.gcf().axes[1].set_ylabel(r"$n_-+n_+$") + + ax.set( + xlabel=r"$x$", + ylabel=r"$y$", + title=rf"$t={{{t:.2f}}}$", + xlim=(0, 1), + ylim=(0, 1), + ) + + +piston = nt2.Data("piston") +piston.makeMovie(plot, framerate=10) diff --git a/pgens/examples/piston/piston.toml b/pgens/examples/piston/piston.toml new file mode 100644 index 000000000..0b5273afa --- /dev/null +++ b/pgens/examples/piston/piston.toml @@ -0,0 +1,57 @@ +[simulation] + name = "piston" + engine = "srpic" + runtime = 2.0 + +[grid] + resolution = [128, 128] + extent = [[0.0, 1.0], [0.0, 1.0]] + + [grid.metric] + metric = "minkowski" + + [grid.boundaries] + fields = [["PERIODIC"], ["PERIODIC"]] + particles = [["ABSORB", "ABSORB"], ["PERIODIC"]] + +[scales] + larmor0 = 0.01 + skindepth0 = 0.1 + +[algorithms] + + [algorithms.timestep] + CFL = 0.5 + +[particles] + ppc0 = 8.0 + + [[particles.species]] + label = "e+" + mass = 1.0 + charge = 1.0 + maxnpart = 1e6 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 1e6 + +[setup] + piston_velocity = 0.5 + +[output] + interval = 10 + + [output.fields] + quantities = ["N_1_2"] + + [output.particles] + enable = false + + [output.spectra] + enable = false + +[checkpoint] + keep = 0 From 6502bb613f56bcd8fb52785a1acfaf11984c5827 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Wed, 22 Apr 2026 22:23:13 +0000 Subject: [PATCH 19/30] moving window example --- pgens/examples/moving_window/moving_window.py | 43 ++++ .../examples/moving_window/moving_window.toml | 57 ++++++ pgens/examples/moving_window/pgen.hpp | 185 ++++++++++++++++++ 3 files changed, 285 insertions(+) create mode 100644 pgens/examples/moving_window/moving_window.py create mode 100644 pgens/examples/moving_window/moving_window.toml create mode 100644 pgens/examples/moving_window/pgen.hpp diff --git a/pgens/examples/moving_window/moving_window.py b/pgens/examples/moving_window/moving_window.py new file mode 100644 index 000000000..e2e81f14b --- /dev/null +++ b/pgens/examples/moving_window/moving_window.py @@ -0,0 +1,43 @@ +import nt2 +import matplotlib.pyplot as plt + + +def plot(t, data): + plt.rcParams["text.usetex"] = True + plt.rcParams["font.family"] = "serif" + plt.rcParams["figure.dpi"] = 250 + ax = plt.gca() + data.fields.N_1_2.sel(t=t, method="nearest").plot( + vmin=0, vmax=1.5, ax=ax, cmap="inferno" + ) + plt.gcf().axes[1].set_ylabel(r"$n_-+n_+$") + + ax.annotate( + "", + xy=(0.28, 0.08), + xycoords="axes fraction", + xytext=(0.08, 0.08), + textcoords="axes fraction", + arrowprops=dict(arrowstyle="-|>", color="#e05c2e", lw=1.8), + ) + ax.text( + 0.095, + 0.095, + r"$v_0$", + transform=ax.transAxes, + color="#e05c2e", + fontsize=12, + va="center", + ) + ax.set( + xlabel=r"$x$", + ylabel=r"$y$", + title=rf"$t={{{t:.2f}}}$", + xlim=(0, 1), # Hayk: This needs to be adaptive + ylim=(0, 1), + ) + + +moving_window = nt2.Data("moving_window") +moving_window.makeMovie(plot, framerate=10) + diff --git a/pgens/examples/moving_window/moving_window.toml b/pgens/examples/moving_window/moving_window.toml new file mode 100644 index 000000000..1e07cf5e2 --- /dev/null +++ b/pgens/examples/moving_window/moving_window.toml @@ -0,0 +1,57 @@ +[simulation] + name = "moving_window" + engine = "srpic" + runtime = 2.0 + +[grid] + resolution = [128, 128] + extent = [[0.0, 1.0], [0.0, 1.0]] + + [grid.metric] + metric = "minkowski" + + [grid.boundaries] + fields = [["PERIODIC"], ["PERIODIC"]] + particles = [["ABSORB", "ABSORB"], ["PERIODIC"]] + +[scales] + larmor0 = 0.01 + skindepth0 = 0.1 + +[algorithms] + + [algorithms.timestep] + CFL = 0.5 + +[particles] + ppc0 = 8.0 + + [[particles.species]] + label = "e+" + mass = 1.0 + charge = 1.0 + maxnpart = 1e6 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 1e6 + +[setup] + window_velocity = 0.2 + +[output] + interval = 10 + + [output.fields] + quantities = ["N_1_2"] + + [output.particles] + enable = false + + [output.spectra] + enable = false + +[checkpoint] + keep = 0 diff --git a/pgens/examples/moving_window/pgen.hpp b/pgens/examples/moving_window/pgen.hpp new file mode 100644 index 000000000..aff89ad39 --- /dev/null +++ b/pgens/examples/moving_window/pgen.hpp @@ -0,0 +1,185 @@ +#ifndef PROBLEM_GENERATOR_H +#define PROBLEM_GENERATOR_H + +#include "enums.h" +#include "global.h" + +#include "traits/pgen.h" + +#include "archetypes/moving_window.h" +#include "archetypes/energy_dist.h" +#include "archetypes/particle_injector.h" +#include "archetypes/problem_generator.h" +#include "archetypes/spatial_dist.h" +#include "framework/domain/metadomain.h" +#include "kernels/particle_moments.hpp" + +/* -------------------------------------------------------------------------- */ +/* Local macros (same as in particle_pusher_sr.hpp) */ +/* -------------------------------------------------------------------------- */ +#define from_Xi_to_i(XI, I) \ + { \ + I = static_cast((XI + 1)) - 1; \ + } + +#define from_Xi_to_i_di(XI, I, DI) \ + { \ + from_Xi_to_i((XI), (I)); \ + DI = static_cast((XI)) - static_cast(I); \ + } + +#define i_di_to_Xi(I, DI) static_cast((I)) + static_cast((DI)) + + +namespace user { + using namespace ntt; + + template + struct NonUniformTargetDensity { + Inline auto operator()(const coord_t& x_Ph) const -> real_t { + // example of a non-uniform target density that peaks at the center of the domain + real_t r2 { ZERO }; + for (auto d = 0u; d < D; ++d) { + r2 += SQR(x_Ph[d] - 0.5); + } + return std::exp( + -r2 / SQR(0.2)); // <-- characteristic width of the density profile + } + }; + + template + struct PGen : public arch::ProblemGenerator { + + static constexpr auto engines { + ::traits::pgen::compatible_with {} + }; + static constexpr auto metrics { + ::traits::pgen::compatible_with {} + }; + static constexpr auto dimensions { + ::traits::pgen::compatible_with {} + }; + + using arch::ProblemGenerator::D; + using arch::ProblemGenerator::C; + using arch::ProblemGenerator::params; + + Metadomain& global_domain; + + int i_window, initial_i_window; + real_t di_window, window_velocity_cd, window_position_cd, window_position; + const real_t window_velocity; + + const real_t dt; + + inline PGen(const SimulationParams& p, Metadomain& global_domain) + : arch::ProblemGenerator { p } + , global_domain { global_domain } + , window_velocity { p.template get("setup.window_velocity", 0.0) } + , dt { p.template get("algorithms.timestep.dt") } {} + + inline void InitPrtls(Domain& domain) { + + // set up window properties + window_position = ZERO; + window_position_cd = domain.mesh.metric.template convert<1, Crd::Ph, Crd::Cd>( + window_position); + // convert to cell units and get cell index and sub-cell position + from_Xi_to_i_di(window_position_cd, i_window, di_window); + // store initial cell index of piston for window updates + initial_i_window = i_window; + + // window velocity in code units + window_velocity_cd = domain.mesh.metric.template transform<1, Idx::XYZ, Idx::U>( + {}, + window_velocity); + + { // compute density of species #1 and #2 + + // saves the density to domain.fields.buff(:,:,:,0) + const auto ni2 = domain.mesh.n_active(in::x2); + const auto inv_n0 = ONE / params.template get("scales.n0"); + + auto scatter_buff = Kokkos::Experimental::create_scatter_view( + domain.fields.buff); + Kokkos::deep_copy(domain.fields.buff, ZERO); + for (const auto sp : std::vector { 1, 2 }) { + const auto& prtl_spec = domain.species[sp - 1]; + Kokkos::parallel_for("ComputeDensity", + prtl_spec.rangeActiveParticles(), + kernel::ParticleMoments_kernel( + {}, + scatter_buff, + 0u, + 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(), + false, + domain.mesh.metric, + domain.mesh.flds_bc(), + ni2, + inv_n0, + 0u)); + } + Kokkos::Experimental::contribute(domain.fields.buff, scatter_buff); + } + // inject particles + const auto energy_dist = arch::Maxwellian( + domain.mesh.metric, + domain.random_pool(), + 0.2, + { window_velocity, ZERO, ZERO }); + + const auto target_density_profile = NonUniformTargetDensity {}; + const auto replenish_sdist = + arch::Replenish( + domain.mesh.metric, + domain.fields.buff, + 0u, // <-- index in buff where the density is stored + target_density_profile, + ONE); // <-- target density for replenishment + arch::InjectNonUniform( + params, + domain, + { 1, 2 }, + { energy_dist, energy_dist }, + replenish_sdist, + ONE); + } + + void CustomPostStep(timestep_t step, simtime_t time, Domain& domain) { + // update window position + di_window += window_velocity_cd * dt; + i_window += static_cast(di_window >= ONE); + di_window -= (di_window >= ONE); + + // check if the window should be moved + if ((i_window - initial_i_window) >= N_GHOSTS) { + + // move the window and all fields and particles in it + arch::MoveWindow(domain, global_domain, N_GHOSTS); + + // update window index for next update + i_window -= N_GHOSTS; + } + } + }; + +} // namespace user + +#undef from_Xi_to_i +#undef from_Xi_to_i_di +#undef i_di_to_Xi + +#endif From 546bdf8403f2cd4312204c84b632704d16e96c8b Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Wed, 22 Apr 2026 22:29:56 +0000 Subject: [PATCH 20/30] removed piston_shock pgen --- pgens/piston_shock/pgen.hpp | 252 --------------------------- pgens/piston_shock/piston_shock.toml | 75 -------- pgens/piston_shock/shock.toml | 74 -------- 3 files changed, 401 deletions(-) delete mode 100644 pgens/piston_shock/pgen.hpp delete mode 100644 pgens/piston_shock/piston_shock.toml delete mode 100644 pgens/piston_shock/shock.toml diff --git a/pgens/piston_shock/pgen.hpp b/pgens/piston_shock/pgen.hpp deleted file mode 100644 index a888add83..000000000 --- a/pgens/piston_shock/pgen.hpp +++ /dev/null @@ -1,252 +0,0 @@ -#ifndef PROBLEM_GENERATOR_H -#define PROBLEM_GENERATOR_H - -#include "enums.h" -#include "global.h" - -#include "archetypes/piston.h" -#include "archetypes/moving_window.h" -#include "archetypes/problem_generator.h" -#include "archetypes/utils.h" -#include "framework/domain/metadomain.h" - -#include -#include - -/* -------------------------------------------------------------------------- */ -/* Local macros (same as in particle_pusher_sr.hpp) */ -/* -------------------------------------------------------------------------- */ -#define from_Xi_to_i(XI, I) \ - { \ - I = static_cast((XI + 1)) - 1; \ - } - -#define from_Xi_to_i_di(XI, I, DI) \ - { \ - from_Xi_to_i((XI), (I)); \ - DI = static_cast((XI)) - static_cast(I); \ - } - -#define i_di_to_Xi(I, DI) static_cast((I)) + static_cast((DI)) - -/* -------------------------------------------------------------------------- */ - -namespace user { - using namespace ntt; - - template - struct InitFields { - /* - Sets up magnetic and electric field components for the simulation. - - @param bmag: magnetic field scaling - @param btheta: magnetic field polar angle - @param bphi: magnetic field azimuthal angle - */ - InitFields(real_t bmag, real_t btheta, real_t bphi) - : Bmag { bmag } - , Btheta { btheta * static_cast(convert::deg2rad) } - , Bphi { bphi * static_cast(convert::deg2rad) } {} - - // magnetic field components - Inline auto bx1(const coord_t&) const -> real_t { - return Bmag * math::cos(Btheta); - } - - Inline auto bx2(const coord_t&) const -> real_t { - return Bmag * math::sin(Btheta) * math::sin(Bphi); - } - - Inline auto bx3(const coord_t&) const -> real_t { - return Bmag * math::sin(Btheta) * math::cos(Bphi); - } - - // electric field components are all zero - - private: - const real_t Bmag, Btheta, Bphi; - }; - - template - struct PGen : public arch::ProblemGenerator { - // compatibility traits for the problem generator - static constexpr auto engines { - ::traits::pgen::compatible_with {} - }; - static constexpr auto metrics { - ::traits::pgen::compatible_with {} - }; - static constexpr auto dimensions { - ::traits::pgen::compatible_with {} - }; - - // for easy access to variables in the child class - using arch::ProblemGenerator::D; - using arch::ProblemGenerator::C; - using arch::ProblemGenerator::params; - - Metadomain& global_domain; - - // domain properties - const real_t cell_size; - // gas properties - const real_t temperature, temperature_ratio; - // injector properties - const real_t dt; - // magnetic field properties - real_t Btheta, Bphi, Bmag; - InitFields init_flds; - - // piston properties - const real_t piston_velocity; - int i_piston, initial_i_piston; - real_t di_piston, piston_velocity_cd, piston_position_cd, piston_position; - - inline PGen(const SimulationParams& p, Metadomain& global_domain) - : arch::ProblemGenerator { p } - , global_domain { global_domain } - , cell_size { (global_domain.mesh().extent(in::x1).second - - global_domain.mesh().extent(in::x1).first) / - global_domain.mesh().n_all(in::x1) } - , temperature { p.template get("setup.temperature") } - , temperature_ratio { p.template get("setup.temperature_ratio") } - , piston_velocity { p.template get("setup.piston_velocity", ZERO) } - , Bmag { p.template get("setup.Bmag", ZERO) } - , Btheta { p.template get("setup.Btheta", ZERO) } - , Bphi { p.template get("setup.Bphi", ZERO) } - , dt { p.template get("algorithms.timestep.dt") } - , init_flds { Bmag, Btheta, Bphi } {} - - inline PGen() {} - - auto MatchFields(simtime_t) const -> InitFields { - return init_flds; - } - - inline void InitPrtls(Domain& domain) { - - // set initial position of piston - piston_position = ZERO; - piston_position_cd = domain.mesh.metric.template convert<1, Crd::Ph, Crd::Cd>( - piston_position); - // convert to cell units and get cell index and sub-cell position - from_Xi_to_i_di(piston_position_cd, i_piston, di_piston); - // store initial cell index of piston for window updates - initial_i_piston = i_piston; - - // piston velocity in code units - coord_t xp_Cd { ZERO }; - piston_velocity_cd = domain.mesh.metric.template transform<1, Idx::XYZ, Idx::U>( - xp_Cd, - piston_velocity); - - // define temperatures of species - const auto temperatures = std::make_pair(temperature, - temperature_ratio * temperature); - - // inject particles - arch::InjectUniformMaxwellians(params, domain, ONE, temperatures, { 1, 2 }); - } - - void CustomPostStep(timestep_t step, simtime_t time, Domain& domain) { - - /* - update piston position - */ - // piston movement over timestep - di_piston += piston_velocity_cd * dt; - // check if the piston has moved to the next cell - i_piston += static_cast(di_piston >= ONE); - // keep track of how much the piston has moved into the next cell - di_piston -= (di_piston >= ONE); - - // check if the window should be moved - if ((i_piston - initial_i_piston) >= N_GHOSTS) { - - // move the window and all fields and particles in it - arch::MoveWindow(domain, global_domain, N_GHOSTS); - - /* - Inject slab of fresh plasma - */ - const real_t xmax = global_domain.mesh().extent(in::x1).second; - const real_t xmin = xmax - N_GHOSTS * cell_size; - // define box to inject into - boundaries_t inj_box; - // loop over all dimension - for (auto d = 0u; d < M::Dim; ++d) { - if (d == 0) { - inj_box.push_back({ xmin, xmax }); - } else { - inj_box.push_back(Range::All); - } - } - - // same maxwell distribution as above - const auto temperatures = std::make_pair(temperature, - temperature_ratio * temperature); - const auto drifts = std::make_pair( - std::vector { ZERO, ZERO, ZERO }, - std::vector { ZERO, ZERO, ZERO }); - arch::InjectUniformMaxwellians(params, - domain, - ONE, - temperatures, - { 1, 2 }, - drifts, - false, - inj_box); - - // shift the piston indices back by the number of ghost cells to account for window movement - i_piston -= N_GHOSTS; - } - - // compute current position of piston - piston_position_cd = i_di_to_Xi(i_piston, di_piston); - // convert to physical coordinates - piston_position = domain.mesh.metric.template convert<1, Crd::Cd, Crd::Ph>( - piston_position_cd); - } - struct CustomPrtlUpdate { - real_t x_piston; // current position of piston - real_t v_piston; - bool is_left; - bool massive; - - Inline void operator()(index_t p, - const kernel::PusherContext& ctx, - const kernel::PusherBoundaries&, - const kernel::PusherArrays& particles, - const M& metric) const { - - if (arch::CrossesPiston(p, - ctx.dt, - particles, - metric, - x_piston, - v_piston, - is_left)) { - arch::PistonUpdate(p, - ctx.dt, - particles, - metric, - x_piston, - v_piston, - massive); - } - } - }; - - template - auto CustomParticleUpdate(simtime_t time, spidx_t sp, DOM& domain) const { - return CustomPrtlUpdate { piston_position, - piston_velocity, - true, - true }; - }; - }; - - -} // namespace user - -#endif // PROBLEM_GENERATOR_H \ No newline at end of file diff --git a/pgens/piston_shock/piston_shock.toml b/pgens/piston_shock/piston_shock.toml deleted file mode 100644 index 6a81d1ab5..000000000 --- a/pgens/piston_shock/piston_shock.toml +++ /dev/null @@ -1,75 +0,0 @@ -[simulation] - name = "piston_shock" - engine = "srpic" - runtime = 100.0 - -[grid] - resolution = [300, 24] - extent = [[0.0, 100.0], [0.0, 8.0]] - - [grid.metric] - metric = "minkowski" - - [grid.boundaries] - fields = [["CONDUCTOR", "CONDUCTOR"], ["PERIODIC"]] - particles = [["REFLECT", "REFLECT"], ["PERIODIC"]] - -[scales] - larmor0 = 5.7735 - skindepth0 = 1.0 - -[algorithms] - current_filters = 8 - - [algorithms.timestep] - CFL = 0.7 - -[particles] - ppc0 = 8.0 - - [[particles.species]] - label = "e-" - mass = 1.0 - charge = -1.0 - maxnpart = 2e7 - - [[particles.species]] - label = "e-" - mass = 1.0 - charge = 1.0 - maxnpart = 2e7 - -[setup] - piston_velocity = 0.99 # speed of leftmoving piston - piston_initial_position = 0.0 # speed towards the wall [c] - temperature = 0.168 # temperature of maxwell distribution [kB T / (m_i c^2)] - temperature_ratio = 1.0 # temperature ratio of electrons to protons - #Bmag = 1.0 # magnetic field strength as fraction of magnetisation - #Btheta = 63.0 # magnetic field angle in the plane - #Bphi = 0.0 # magnetic field angle out of plane - #filling_fraction = 0.99 # fraction of the shock piston filled with plasma - #injector_velocity = 0.0 # speed of injector [c] - #injection_start = 0.0 # start time of moving injector - #injection_frequency = 100 - -[output] - interval_time = 1000.0 - format = "BPFile" - - [output.fields] - quantities = ["N", "B", "E"] - - [output.particles] - enable = true - stride = 10 - - [output.spectra] - enable = false - -[diagnostics] - log_level = "WARNING" - blocking_timers = true - -[checkpoint] - interval = 10000 - keep = 1 \ No newline at end of file diff --git a/pgens/piston_shock/shock.toml b/pgens/piston_shock/shock.toml deleted file mode 100644 index 8cde86990..000000000 --- a/pgens/piston_shock/shock.toml +++ /dev/null @@ -1,74 +0,0 @@ -[simulation] - name = "shock" - engine = "srpic" - runtime = 15200.0 - -[grid] - resolution = [32768, 256] - extent = [[0.0, 8192.0], [-32.0, 32.0]] - - [grid.metric] - metric = "minkowski" - - [grid.boundaries] - fields = [["CONDUCTOR", "MATCH"], ["PERIODIC"]] - particles = [["REFLECT", "ABSORB"], ["PERIODIC"]] - -[scales] - larmor0 = 5.7735 - skindepth0 = 1.0 - -[algorithms] - current_filters = 8 - - [algorithms.timestep] - CFL = 0.7 - -[particles] - ppc0 = 8.0 - - [[particles.species]] - label = "e-" - mass = 1.0 - charge = -1.0 - maxnpart = 2e7 - - [[particles.species]] - label = "p+" - mass = 100.0 - charge = 1.0 - maxnpart = 2e7 - -[setup] - drift_ux = 0.15 # speed towards the wall [c] - temperature = 0.168 # temperature of maxwell distribution [kB T / (m_i c^2)] - temperature_ratio = 1.0 # temperature ratio of electrons to protons - Bmag = 1.0 # magnetic field strength as fraction of magnetisation - Btheta = 63.0 # magnetic field angle in the plane - Bphi = 0.0 # magnetic field angle out of plane - filling_fraction = 0.99 # fraction of the shock piston filled with plasma - injector_velocity = 0.0 # speed of injector [c] - injection_start = 0.0 # start time of moving injector - injection_frequency = 100 - -[output] - interval_time = 1000.0 - format = "BPFile" - - [output.fields] - quantities = ["N", "B", "E"] - - [output.particles] - enable = true - stride = 10 - - [output.spectra] - enable = false - -[diagnostics] - log_level = "WARNING" - blocking_timers = true - -[checkpoint] - interval = 10000 - keep = 1 \ No newline at end of file From c788022fd6a90e5e04b21e4cf1780f7fd2e4b918 Mon Sep 17 00:00:00 2001 From: hayk Date: Sat, 25 Apr 2026 18:19:22 -0400 Subject: [PATCH 21/30] fixed & simplified piston & replenish injectors --- pgens/examples/piston/pgen.hpp | 106 ++++++--------------- pgens/examples/replenish_injector/pgen.hpp | 51 ++-------- 2 files changed, 37 insertions(+), 120 deletions(-) diff --git a/pgens/examples/piston/pgen.hpp b/pgens/examples/piston/pgen.hpp index 1908b1514..fe8fb6d4b 100644 --- a/pgens/examples/piston/pgen.hpp +++ b/pgens/examples/piston/pgen.hpp @@ -4,34 +4,35 @@ #include "enums.h" #include "global.h" +#include "arch/kokkos_aliases.h" #include "traits/pgen.h" #include "archetypes/energy_dist.h" #include "archetypes/particle_injector.h" #include "archetypes/piston.h" -#include "archetypes/problem_generator.h" -#include "archetypes/spatial_dist.h" +#include "framework/domain/domain.h" #include "framework/domain/metadomain.h" -#include "kernels/particle_moments.hpp" +#include "framework/parameters/parameters.h" namespace user { using namespace ntt; template - struct NonUniformTargetDensity { + struct NonUniformDensity { Inline auto operator()(const coord_t& x_Ph) const -> real_t { // example of a non-uniform target density that peaks at the center of the domain real_t r2 { ZERO }; for (auto d = 0u; d < D; ++d) { r2 += SQR(x_Ph[d] - 0.5); } - return std::exp( - -r2 / SQR(0.2)); // <-- characteristic width of the density profile + return math::exp(-r2 / SQR(static_cast(0.2))); + // ^ + // characteristic width of the density profile } }; template - struct PGen : public arch::ProblemGenerator { + struct PGen { static constexpr auto engines { ::traits::pgen::compatible_with {} @@ -43,80 +44,29 @@ namespace user { ::traits::pgen::compatible_with {} }; - using arch::ProblemGenerator::D; - using arch::ProblemGenerator::C; - using arch::ProblemGenerator::params; - - Metadomain& global_domain; - const real_t piston_velocity; - - inline PGen(const SimulationParams& p, Metadomain& global_domain) - : arch::ProblemGenerator { p } - , global_domain { global_domain } - , piston_velocity { p.template get("setup.piston_velocity", 0.0) } {} - - inline void InitPrtls(Domain& domain) { - - { // compute density of species #1 and #2 - - // saves the density to domain.fields.buff(:,:,:,0) - const auto ni2 = domain.mesh.n_active(in::x2); - const auto inv_n0 = ONE / params.template get("scales.n0"); - - auto scatter_buff = Kokkos::Experimental::create_scatter_view( - domain.fields.buff); - Kokkos::deep_copy(domain.fields.buff, ZERO); - for (const auto sp : std::vector { 1, 2 }) { - const auto& prtl_spec = domain.species[sp - 1]; - Kokkos::parallel_for("ComputeDensity", - prtl_spec.rangeActiveParticles(), - kernel::ParticleMoments_kernel( - {}, - scatter_buff, - 0u, - 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(), - false, - domain.mesh.metric, - domain.mesh.flds_bc(), - ni2, - inv_n0, - 0u)); - } - Kokkos::Experimental::contribute(domain.fields.buff, scatter_buff); - } + const SimulationParams& params; + Metadomain& metadomain; + const real_t piston_velocity; + + PGen(const SimulationParams& p, Metadomain& m) + : params { p } + , metadomain { m } + , piston_velocity { params.template get("setup.piston_velocity", + 0.0) } {} + + void InitPrtls(Domain& domain) { // inject particles - const auto energy_dist = arch::Maxwellian( - domain.mesh.metric, + const auto energy_dist = arch::energy_dist::Maxwellian( domain.random_pool(), 0.2); // <-- target temperature for injection - const auto target_density_profile = NonUniformTargetDensity {}; - const auto replenish_sdist = - arch::Replenish( - domain.mesh.metric, - domain.fields.buff, - 0u, // <-- index in buff where the density is stored - target_density_profile, - ONE); // <-- target density for replenishment - arch::InjectNonUniform( + const auto density_profile = NonUniformDensity {}; + arch::InjectNonUniform( params, domain, { 1, 2 }, { energy_dist, energy_dist }, - replenish_sdist, + density_profile, ONE); } @@ -127,9 +77,9 @@ namespace user { bool is_left; bool massive; - Inline void operator()(index_t p, - const kernel::PusherContext& ctx, - const kernel::PusherBoundaries&, + Inline void operator()(index_t p, + const kernel::sr::PusherContext& ctx, + const kernel::sr::PusherBoundaries&, const kernel::PusherArrays& particles, const M& metric) const { @@ -146,10 +96,10 @@ namespace user { }; template - auto CustomParticleUpdate(simtime_t time, spidx_t sp, DOM& domain) const { + auto CustomParticleUpdate(simtime_t time, spidx_t /*sp*/, DOM& /*domain*/) const { return CustomPrtlUpdate { piston_velocity * static_cast(time), piston_velocity, - global_domain.mesh().extent(in::x1).second, + metadomain.mesh().extent(in::x1).second, true, true }; }; diff --git a/pgens/examples/replenish_injector/pgen.hpp b/pgens/examples/replenish_injector/pgen.hpp index fb46eaca8..8081f8f03 100644 --- a/pgens/examples/replenish_injector/pgen.hpp +++ b/pgens/examples/replenish_injector/pgen.hpp @@ -4,13 +4,14 @@ #include "enums.h" #include "global.h" +#include "arch/kokkos_aliases.h" #include "traits/pgen.h" #include "archetypes/energy_dist.h" #include "archetypes/particle_injector.h" #include "archetypes/spatial_dist.h" +#include "archetypes/utils.h" #include "framework/domain/metadomain.h" -#include "kernels/particle_moments.hpp" namespace user { using namespace ntt; @@ -34,7 +35,7 @@ namespace user { for (auto d = 0u; d < D; ++d) { r2 += SQR(x_Ph[d] - HALF); } - return std::exp(-r2 / SQR(static_cast(0.2))); + return math::exp(-r2 / SQR(static_cast(0.2))); // ^ // characteristic width of the density profile } @@ -70,46 +71,12 @@ namespace user { } // perform replenishment and injection every 100 timesteps - { // compute density of species #1 and #2 - - // saves the density to domain.fields.buff(:,:,:,0) - const auto ni2 = domain.mesh.n_active(in::x2); - const auto inv_n0 = ONE / params.template get("scales.n0"); - - auto scatter_buff = Kokkos::Experimental::create_scatter_view( - domain.fields.buff); - Kokkos::deep_copy(domain.fields.buff, ZERO); - for (const auto sp : std::vector { 1, 2 }) { - const auto& prtl_spec = domain.species[sp - 1]; - Kokkos::parallel_for("ComputeDensity", - prtl_spec.rangeActiveParticles(), - kernel::ParticleMoments_kernel( - {}, - scatter_buff, - 0u, - 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(), - false, - domain.mesh.metric, - domain.mesh.flds_bc(), - ni2, - inv_n0, - 0u)); - } - Kokkos::Experimental::contribute(domain.fields.buff, scatter_buff); - } + // compute the density of species #1 and #2 + // and save in the field buffer (index 0) + arch::ComputeMomentWithSpecies(params, + domain, + { 1u, 2u }, + domain.fields.buff); const auto energy_dist = arch::energy_dist::Maxwellian( domain.random_pool(), From 4517d933ea30c3f91a14ad2632779b4a68429d8d Mon Sep 17 00:00:00 2001 From: hayk Date: Sat, 25 Apr 2026 18:52:25 -0400 Subject: [PATCH 22/30] fixed moving window pgen --- pgens/examples/moving_window/pgen.hpp | 187 ++++++++------------------ src/archetypes/moving_window.h | 14 +- 2 files changed, 63 insertions(+), 138 deletions(-) diff --git a/pgens/examples/moving_window/pgen.hpp b/pgens/examples/moving_window/pgen.hpp index aff89ad39..7986b9ae5 100644 --- a/pgens/examples/moving_window/pgen.hpp +++ b/pgens/examples/moving_window/pgen.hpp @@ -4,52 +4,35 @@ #include "enums.h" #include "global.h" +#include "arch/kokkos_aliases.h" #include "traits/pgen.h" -#include "archetypes/moving_window.h" #include "archetypes/energy_dist.h" +#include "archetypes/moving_window.h" #include "archetypes/particle_injector.h" -#include "archetypes/problem_generator.h" -#include "archetypes/spatial_dist.h" +#include "framework/domain/domain.h" #include "framework/domain/metadomain.h" -#include "kernels/particle_moments.hpp" - -/* -------------------------------------------------------------------------- */ -/* Local macros (same as in particle_pusher_sr.hpp) */ -/* -------------------------------------------------------------------------- */ -#define from_Xi_to_i(XI, I) \ - { \ - I = static_cast((XI + 1)) - 1; \ - } - -#define from_Xi_to_i_di(XI, I, DI) \ - { \ - from_Xi_to_i((XI), (I)); \ - DI = static_cast((XI)) - static_cast(I); \ - } - -#define i_di_to_Xi(I, DI) static_cast((I)) + static_cast((DI)) - +#include "framework/parameters/parameters.h" namespace user { using namespace ntt; template - struct NonUniformTargetDensity { + struct NonUniformDensity { Inline auto operator()(const coord_t& x_Ph) const -> real_t { // example of a non-uniform target density that peaks at the center of the domain real_t r2 { ZERO }; for (auto d = 0u; d < D; ++d) { r2 += SQR(x_Ph[d] - 0.5); } - return std::exp( - -r2 / SQR(0.2)); // <-- characteristic width of the density profile + return math::exp(-r2 / SQR(static_cast(0.1))); + // ^ + // characteristic width of the density profile } }; template - struct PGen : public arch::ProblemGenerator { - + struct PGen { static constexpr auto engines { ::traits::pgen::compatible_with {} }; @@ -59,127 +42,73 @@ namespace user { static constexpr auto dimensions { ::traits::pgen::compatible_with {} }; + const SimulationParams& params; + const Metadomain& metadomain; + + struct MovingWindow { + int pos_i { -1 }; + int init_pos_i { -1 }; + real_t pos_di { ZERO }; + real_t vel_Cd { ZERO }; + + void init(const M& metric, real_t global_x, real_t velocity) { + const auto pos_Cd = metric.template convert<1, Crd::Ph, Crd::Cd>(global_x); + pos_i = static_cast(pos_Cd + 1) - 1; + init_pos_i = pos_i; + pos_di = pos_Cd - static_cast(pos_i); + vel_Cd = metric.template transform<1, Idx::XYZ, Idx::U>({}, velocity); + } - using arch::ProblemGenerator::D; - using arch::ProblemGenerator::C; - using arch::ProblemGenerator::params; - - Metadomain& global_domain; - - int i_window, initial_i_window; - real_t di_window, window_velocity_cd, window_position_cd, window_position; - const real_t window_velocity; + void update(real_t dt, + int shift, + const Metadomain& metadomain, + Domain& domain) { + pos_di += vel_Cd * dt; + pos_i += static_cast(pos_di >= ONE); + pos_di -= static_cast(pos_di >= ONE); + if ((pos_i - init_pos_i) >= shift) { + // move the window and all fields and particles in it + arch::MoveWindow(domain, metadomain, shift); + // update window index for next update + pos_i -= shift; + } + } + }; + MovingWindow moving_window; const real_t dt; - inline PGen(const SimulationParams& p, Metadomain& global_domain) - : arch::ProblemGenerator { p } - , global_domain { global_domain } - , window_velocity { p.template get("setup.window_velocity", 0.0) } - , dt { p.template get("algorithms.timestep.dt") } {} - - inline void InitPrtls(Domain& domain) { - - // set up window properties - window_position = ZERO; - window_position_cd = domain.mesh.metric.template convert<1, Crd::Ph, Crd::Cd>( - window_position); - // convert to cell units and get cell index and sub-cell position - from_Xi_to_i_di(window_position_cd, i_window, di_window); - // store initial cell index of piston for window updates - initial_i_window = i_window; - - // window velocity in code units - window_velocity_cd = domain.mesh.metric.template transform<1, Idx::XYZ, Idx::U>( - {}, - window_velocity); - - { // compute density of species #1 and #2 - - // saves the density to domain.fields.buff(:,:,:,0) - const auto ni2 = domain.mesh.n_active(in::x2); - const auto inv_n0 = ONE / params.template get("scales.n0"); - - auto scatter_buff = Kokkos::Experimental::create_scatter_view( - domain.fields.buff); - Kokkos::deep_copy(domain.fields.buff, ZERO); - for (const auto sp : std::vector { 1, 2 }) { - const auto& prtl_spec = domain.species[sp - 1]; - Kokkos::parallel_for("ComputeDensity", - prtl_spec.rangeActiveParticles(), - kernel::ParticleMoments_kernel( - {}, - scatter_buff, - 0u, - 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(), - false, - domain.mesh.metric, - domain.mesh.flds_bc(), - ni2, - inv_n0, - 0u)); - } - Kokkos::Experimental::contribute(domain.fields.buff, scatter_buff); - } + PGen(const SimulationParams& p, Metadomain& m) + : params { p } + , metadomain { m } + , dt { params.template get("algorithms.timestep.dt") } {} + + void InitPrtls(Domain& domain) { + const auto window_velocity = params.template get( + "setup.window_velocity", + ZERO); + moving_window.init(domain.mesh.metric, ZERO, window_velocity); // inject particles - const auto energy_dist = arch::Maxwellian( - domain.mesh.metric, + const auto energy_dist = arch::energy_dist::Maxwellian( domain.random_pool(), - 0.2, + 0.001, { window_velocity, ZERO, ZERO }); - const auto target_density_profile = NonUniformTargetDensity {}; - const auto replenish_sdist = - arch::Replenish( - domain.mesh.metric, - domain.fields.buff, - 0u, // <-- index in buff where the density is stored - target_density_profile, - ONE); // <-- target density for replenishment - arch::InjectNonUniform( + const auto density_profile = NonUniformDensity {}; + arch::InjectNonUniform( params, domain, { 1, 2 }, { energy_dist, energy_dist }, - replenish_sdist, + density_profile, ONE); } - void CustomPostStep(timestep_t step, simtime_t time, Domain& domain) { - // update window position - di_window += window_velocity_cd * dt; - i_window += static_cast(di_window >= ONE); - di_window -= (di_window >= ONE); - - // check if the window should be moved - if ((i_window - initial_i_window) >= N_GHOSTS) { - - // move the window and all fields and particles in it - arch::MoveWindow(domain, global_domain, N_GHOSTS); - - // update window index for next update - i_window -= N_GHOSTS; - } + void CustomPostStep(timestep_t /*step*/, simtime_t /*time*/, Domain& domain) { + moving_window.update(dt, N_GHOSTS, metadomain, domain); } }; } // namespace user -#undef from_Xi_to_i -#undef from_Xi_to_i_di -#undef i_di_to_Xi - #endif diff --git a/src/archetypes/moving_window.h b/src/archetypes/moving_window.h index c0967df3d..47a4e8d19 100644 --- a/src/archetypes/moving_window.h +++ b/src/archetypes/moving_window.h @@ -13,12 +13,8 @@ #include "enums.h" #include "global.h" -#include "archetypes/energy_dist.h" #include "framework/domain/domain.h" #include "framework/domain/metadomain.h" -#include "framework/parameters/parameters.h" - -#include namespace arch { @@ -141,9 +137,9 @@ namespace arch { */ template - Inline void MoveWindow(Domain& domain, - Metadomain& global_domain, - const int window_shift) { + Inline void MoveWindow(Domain& domain, + const Metadomain& metadomain, + int window_shift) { /* move particles in the window back by the window size @@ -271,9 +267,9 @@ namespace arch { } // synch ghost zones after moving the window - global_domain.CommunicateFields(domain, Comm::E | Comm::B); + metadomain.CommunicateFields(domain, Comm::E | Comm::B); // communicate particles after moving - global_domain.CommunicateParticles(domain); + metadomain.CommunicateParticles(domain); // ToDo: Update metric } From e3356bb31945e42f80395cec1cda083010bafc80 Mon Sep 17 00:00:00 2001 From: hayk Date: Sat, 25 Apr 2026 20:06:57 -0400 Subject: [PATCH 23/30] moving window now updates metric --- pgens/examples/moving_window/pgen.hpp | 10 ++-- src/archetypes/moving_window.h | 8 ++-- src/framework/CMakeLists.txt | 2 + src/framework/domain/mesh.h | 15 +++++- src/framework/domain/metadomain.h | 12 +++-- src/framework/domain/metadomain_reshape.cpp | 53 +++++++++++++++++++++ src/framework/specialization_registry.h | 5 ++ src/metrics/minkowski.h | 5 ++ 8 files changed, 96 insertions(+), 14 deletions(-) create mode 100644 src/framework/domain/metadomain_reshape.cpp diff --git a/pgens/examples/moving_window/pgen.hpp b/pgens/examples/moving_window/pgen.hpp index 7986b9ae5..58d4ff103 100644 --- a/pgens/examples/moving_window/pgen.hpp +++ b/pgens/examples/moving_window/pgen.hpp @@ -43,7 +43,7 @@ namespace user { ::traits::pgen::compatible_with {} }; const SimulationParams& params; - const Metadomain& metadomain; + Metadomain& metadomain; struct MovingWindow { int pos_i { -1 }; @@ -59,10 +59,10 @@ namespace user { vel_Cd = metric.template transform<1, Idx::XYZ, Idx::U>({}, velocity); } - void update(real_t dt, - int shift, - const Metadomain& metadomain, - Domain& domain) { + void update(real_t dt, + int shift, + Metadomain& metadomain, + Domain& domain) { pos_di += vel_Cd * dt; pos_i += static_cast(pos_di >= ONE); pos_di -= static_cast(pos_di >= ONE); diff --git a/src/archetypes/moving_window.h b/src/archetypes/moving_window.h index 47a4e8d19..b4ef21558 100644 --- a/src/archetypes/moving_window.h +++ b/src/archetypes/moving_window.h @@ -137,9 +137,9 @@ namespace arch { */ template - Inline void MoveWindow(Domain& domain, - const Metadomain& metadomain, - int window_shift) { + Inline void MoveWindow(Domain& domain, + Metadomain& metadomain, + int window_shift) { /* move particles in the window back by the window size @@ -271,7 +271,7 @@ namespace arch { // communicate particles after moving metadomain.CommunicateParticles(domain); - // ToDo: Update metric + metadomain.ShiftByCells(window_shift, o); } } // namespace arch diff --git a/src/framework/CMakeLists.txt b/src/framework/CMakeLists.txt index 032e507fb..48187d470 100644 --- a/src/framework/CMakeLists.txt +++ b/src/framework/CMakeLists.txt @@ -18,6 +18,7 @@ # * domain/metadomain_chckpt.cpp # * domain/metadomain_stats.cpp # * domain/metadomain_io.cpp +# * domain/metadomain_reshape.cpp # * containers/particles.cpp # * containers/particles_comm.cpp # * containers/particles_io.cpp @@ -57,6 +58,7 @@ set(SOURCES ${SRC_DIR}/domain/metadomain_comm.cpp ${SRC_DIR}/domain/metadomain_sort.cpp ${SRC_DIR}/domain/metadomain_stats.cpp + ${SRC_DIR}/domain/metadomain_reshape.cpp ${SRC_DIR}/containers/particles.cpp ${SRC_DIR}/containers/particles_sort.cpp ${SRC_DIR}/containers/fields.cpp) diff --git a/src/framework/domain/mesh.h b/src/framework/domain/mesh.h index 07d1d1500..aa2eea726 100644 --- a/src/framework/domain/mesh.h +++ b/src/framework/domain/mesh.h @@ -45,7 +45,8 @@ namespace ntt { const boundaries_t& ext, const std::map& metric_params) : Grid { res, ext } - , metric { res, ext, metric_params } {} + , metric { res, ext, metric_params } + , m_metric_params_raw { metric_params } {} Mesh(const std::vector& res, const boundaries_t& ext, @@ -53,10 +54,17 @@ namespace ntt { const boundaries_t& flds_bc, const boundaries_t& prtl_bc) : Grid { res, ext, flds_bc, prtl_bc } - , metric { res, ext, metric_params } {} + , metric { res, ext, metric_params } + , m_metric_params_raw { metric_params } {} ~Mesh() = default; + void set_extent(const boundaries_t& new_extent) { + m_extent = new_extent; + metric.~M(); + new (&metric) M { this->m_resolution, new_extent, m_metric_params_raw }; + } + /** * @brief Get the intersection of the mesh with a box * @param box physical extent @@ -193,6 +201,9 @@ namespace ntt { return range; } + + private: + std::map m_metric_params_raw; }; } // namespace ntt diff --git a/src/framework/domain/metadomain.h b/src/framework/domain/metadomain.h index 9c67a25a5..d7971d858 100644 --- a/src/framework/domain/metadomain.h +++ b/src/framework/domain/metadomain.h @@ -10,6 +10,7 @@ * - metadomain_chckpt.cpp * - metadomain_io.cpp * - metadomain_stats.cpp + * - metadomain_reshape.cpp * @namespaces: * - ntt:: * @macros: @@ -95,6 +96,9 @@ namespace ntt { void SynchronizeFields(Domain&, CommTags, const range_tuple_t& = { 0, 0 }) const; +#if defined(MPI_ENABLED) && defined(OUTPUT_ENABLED) + void CommunicateVectorPotential(unsigned short); +#endif void CommunicateParticles(Domain&) const; void SortParticles(simtime_t, timestep_t, @@ -125,6 +129,11 @@ namespace ntt { ~Metadomain() = default; + /* domain update-related ------------------------------------------------ */ + void ShiftByCells(int, in = in::x1) + requires(CartesianMetricClass); + + /* output-related ------------------------------------------------------- */ #if defined(OUTPUT_ENABLED) void InitWriter(adios2::ADIOS*, const SimulationParams&); auto Write(const SimulationParams&, @@ -278,9 +287,6 @@ namespace ntt { #if defined(OUTPUT_ENABLED) out::Writer g_writer; checkpoint::Writer g_checkpoint_writer; - #if defined(MPI_ENABLED) - void CommunicateVectorPotential(unsigned short); - #endif #endif #if defined(MPI_ENABLED) diff --git a/src/framework/domain/metadomain_reshape.cpp b/src/framework/domain/metadomain_reshape.cpp new file mode 100644 index 000000000..5104bcd6f --- /dev/null +++ b/src/framework/domain/metadomain_reshape.cpp @@ -0,0 +1,53 @@ +#include "enums.h" +#include "global.h" + +#include "traits/metric.h" +#include "utils/error.h" + +#include "framework/domain/mesh.h" +#include "framework/domain/metadomain.h" +#include "framework/specialization_registry.h" + +namespace ntt { + + template + void Metadomain::ShiftByCells(int n_cells, in dir) + requires(CartesianMetricClass) + { + raise::ErrorIf( + n_cells == 0, + "ShiftByCells called with n_cells = 0, no shift will be performed", + HERE); + + { + // update metadomain g_mesh + auto new_extent = g_mesh.extent(); + const auto delta_shift = static_cast(n_cells) * + g_mesh.metric.get_dx(); + new_extent[static_cast(dir)].first += delta_shift; + new_extent[static_cast(dir)].second += delta_shift; + + g_mesh.set_extent(new_extent); + } + + // update individual subdomain meshes + for (auto& subdomain : g_subdomains) { + auto new_extent = subdomain.mesh.extent(); + const auto delta_shift = static_cast(n_cells) * + subdomain.mesh.metric.get_dx(); + new_extent[static_cast(dir)].first += delta_shift; + new_extent[static_cast(dir)].second += delta_shift; + + subdomain.mesh.set_extent(new_extent); + } + } + + // NOLINTBEGIN(bugprone-macro-parentheses) +#define METADOMAIN_SHIFT(S, M, D) \ + template void Metadomain>::ShiftByCells(int n_cells, in dir); + + NTT_FOREACH_CARTESIAN_SPECIALIZATION(METADOMAIN_SHIFT) +#undef METADOMAIN_SHIFT + // NOLINTEND(bugprone-macro-parentheses) + +} // namespace ntt \ No newline at end of file diff --git a/src/framework/specialization_registry.h b/src/framework/specialization_registry.h index d9c1628ca..dea5772c6 100644 --- a/src/framework/specialization_registry.h +++ b/src/framework/specialization_registry.h @@ -47,6 +47,11 @@ namespace ntt { MACRO(SimEngine::GRPIC, metric::QKerrSchild, Dim::_2D) \ MACRO(SimEngine::GRPIC, metric::KerrSchild0, Dim::_2D) +#define NTT_FOREACH_CARTESIAN_SPECIALIZATION(MACRO) \ + MACRO(SimEngine::SRPIC, metric::Minkowski, Dim::_1D) \ + MACRO(SimEngine::SRPIC, metric::Minkowski, Dim::_2D) \ + MACRO(SimEngine::SRPIC, metric::Minkowski, Dim::_3D) + #define NTT_BUILD_SPECIALIZATION_ENTRY(S, M, D) SpecializationEntry {}, inline constexpr auto kSpecializations = std::tuple { NTT_FOREACH_SPECIALIZATION( diff --git a/src/metrics/minkowski.h b/src/metrics/minkowski.h index 8d6c87e5f..3a679b389 100644 --- a/src/metrics/minkowski.h +++ b/src/metrics/minkowski.h @@ -70,6 +70,11 @@ namespace metric { } } + [[nodiscard]] + auto get_dx() const -> real_t { + return dx; + } + /** * minimum effective cell size for a given metric (in physical units) */ From cfda16ddb22f1e4283b5600befb99e9ef5d1bc31 Mon Sep 17 00:00:00 2001 From: hayk Date: Sat, 25 Apr 2026 20:24:08 -0400 Subject: [PATCH 24/30] minor --- pgens/examples/moving_window/moving_window.py | 3 - pgens/tutorial/pgen.hpp | 112 ------------------ pgens/tutorial/tutorial.toml | 49 -------- 3 files changed, 164 deletions(-) delete mode 100644 pgens/tutorial/pgen.hpp delete mode 100644 pgens/tutorial/tutorial.toml diff --git a/pgens/examples/moving_window/moving_window.py b/pgens/examples/moving_window/moving_window.py index e2e81f14b..234dc7892 100644 --- a/pgens/examples/moving_window/moving_window.py +++ b/pgens/examples/moving_window/moving_window.py @@ -33,11 +33,8 @@ def plot(t, data): xlabel=r"$x$", ylabel=r"$y$", title=rf"$t={{{t:.2f}}}$", - xlim=(0, 1), # Hayk: This needs to be adaptive - ylim=(0, 1), ) moving_window = nt2.Data("moving_window") moving_window.makeMovie(plot, framerate=10) - diff --git a/pgens/tutorial/pgen.hpp b/pgens/tutorial/pgen.hpp deleted file mode 100644 index f90d154de..000000000 --- a/pgens/tutorial/pgen.hpp +++ /dev/null @@ -1,112 +0,0 @@ -#ifndef PROBLEM_GENERATOR_H -#define PROBLEM_GENERATOR_H - -#include "enums.h" -#include "global.h" - -#include "archetypes/utils.h" - -#include "archetypes/particle_injector.h" -#include "archetypes/traits.h" -#include "archetypes/problem_generator.h" -#include "framework/domain/metadomain.h" -#include "framework/parameters/parameters.h" - -namespace user { - using namespace ntt; - - //let us define a magnetic field - template - // we define a dipole - struct DipoleField{ - // this is a function for the radius from origin - Inline auto radius(const coord_t& x ) const -> real_t{ - if constexpr (D==Dim::_3D){ - return math::sqrt(SQR(x[0])+SQR(x[1])+SQR(x[2])); - } else { - return math::sqrt(SQR(x[0])+SQR(x[1])); - } - - } - Inline auto bx1(const coord_t& x) const -> real_t{ - return THREE*x[0]*x[1]/math::pow(radius(x), 5); - } - Inline auto bx2(const coord_t& x) const -> real_t{ - const auto r = radius(x); - return (THREE*x[1]*x[1]-SQR(r))/math::pow(r, 5); - } - Inline auto bx3(const coord_t& x) const -> real_t{ - if constexpr (D == Dim::_3D) { - return THREE*x[2]*x[1]/math::pow(radius(x), 5); - } else { - return ZERO; - } - - } - }; - - template - struct PGen : public arch::ProblemGenerator { - static constexpr auto engines { - arch::traits::pgen::compatible_with::value - }; - static constexpr auto metrics { - arch::traits::pgen::compatible_with::value - }; - static constexpr auto dimensions { - arch::traits::pgen::compatible_with::value - }; - // let us now define an init_flds - //DipoleField init_flds; - - //Now let us define this dipole as a sort of external field - inline auto ExternalFields(simtime_t, spidx_t, const Domain&) const - -> std::pair> { - return {true, DipoleField {}}; // first part of this is whether external fields need to be applied - // second return is the field structure itself - // in theory this can be a function of time, or species index - } - - const Metadomain& metadomain; - - inline PGen(const SimulationParams& p, const Metadomain& metadomain) - : arch::ProblemGenerator {p} - , metadomain {metadomain}{} - - inline void InitPrtls(Domain& domain) { - arch::InjectUniformMaxwellian(this->params, - domain, - ONE, // number density in ppc - 1e-3, // this is the temperature - {1u, 2u}, // initialize two particles - { {ONE, ZERO, ZERO}, {ONE, ZERO, ZERO}}, // velocity - false, // do not use particle weights - {{-2.0, -1.5}, Range::All, Range::All}); // adding this last ine designates the strip where there will be plasma - // if this last line were not there, I would inject plasma in the whole domain - //std::map> particles { - // {"x1", {} }, - // {"x2", {} }, - // {"x3", {} }, - // {"ux1", {} }, - // {"ux2", {} }, - // {"ux3", {} } - //}; - //const auto npart = 1000u; - //const auto ymin = metadomain.mesh().extent(in::x2).first; // pick min y - //const auto ymax = metadomain.mesh().extent(in::x2).second; // pick max y - //const auto dy = (ymax - ymin) / static_cast(npart - 1u); // define y spacing - //for (auto p { 0u }; p < npart; ++p) { - // particles["x1"].push_back(-1.5); // we stick them at x1 = -1.5 - // particles["x2"].push_back(ymin + p * dy); // here we are evenly spacing the particles in y - // particles["x3"].push_back(0.0); - // particles["ux1"].push_back(1.0); // define them as moving right - // particles["ux2"].push_back(0.0); - // particles["ux3"].push_back(0.0); - //} - //arch::InjectGlobally(metadomain, domain, 1u, particles); - // third entry in this function is the species - } - }; -} - -#endif \ No newline at end of file diff --git a/pgens/tutorial/tutorial.toml b/pgens/tutorial/tutorial.toml deleted file mode 100644 index eb76dc27f..000000000 --- a/pgens/tutorial/tutorial.toml +++ /dev/null @@ -1,49 +0,0 @@ -[simulation] - name = "my_first_pgen" - engine = "srpic" - runtime = 1.0 - -[grid] - resolution = [1024, 1024] - extent = [[-2.0, 2.0], [-2.0, 2.0]] - - [grid.metric] - metric = "minkowski" - - [grid.boundaries] - fields = [["PERIODIC"], ["PERIODIC"]] - particles = [["PERIODIC"], ["PERIODIC"]] - -[scales] - larmor0 = 0.003 - skindepth0 = 0.01 - -[algorithms.deposit] - enable = false - -[particles] - ppc0 = 8.0 - - [[particles.species]] - label = "e+" - charge = 1.0 - mass = 1.0 - maxnpart = 1e5 - tracking = true - - -[output] - interval_time = 0.1 - - [output.fields] - quantities = ["N_1", "E", "B"] - - [output.particles] - interval_time = 0.01 - stride = 1 - - [output.spectra] - enable = false - -[checkpoint] - keep = 0 \ No newline at end of file From bed0625655560bbd83cfdc79de61266f0db2c3f7 Mon Sep 17 00:00:00 2001 From: hayk Date: Sun, 26 Apr 2026 03:26:37 -0400 Subject: [PATCH 25/30] changing size domain compatible with restart --- src/framework/domain/metadomain.h | 2 + src/framework/domain/metadomain_chckpt.cpp | 187 +++++++++++++++- tests/framework/CMakeLists.txt | 10 + tests/framework/checkpoint-data.cpp | 244 +++++++++++++++++++++ tests/framework/checkpoint-extents.cpp | 163 ++++++++++++++ 5 files changed, 603 insertions(+), 3 deletions(-) create mode 100644 tests/framework/checkpoint-data.cpp create mode 100644 tests/framework/checkpoint-extents.cpp diff --git a/src/framework/domain/metadomain.h b/src/framework/domain/metadomain.h index d7971d858..c0a11efa2 100644 --- a/src/framework/domain/metadomain.h +++ b/src/framework/domain/metadomain.h @@ -155,6 +155,8 @@ namespace ntt { simtime_t) -> bool; void ContinueFromCheckpoint(adios2::ADIOS*, const SimulationParams&); + void redecomposeFromCheckpoint(const std::vector>&, + const std::vector>&); #endif void InitStatsWriter(const SimulationParams&, bool); diff --git a/src/framework/domain/metadomain_chckpt.cpp b/src/framework/domain/metadomain_chckpt.cpp index 46aa32b84..e4622f92a 100644 --- a/src/framework/domain/metadomain_chckpt.cpp +++ b/src/framework/domain/metadomain_chckpt.cpp @@ -10,12 +10,15 @@ #include "framework/parameters/parameters.h" #include "framework/specialization_registry.h" #include "output/checkpoint.h" +#include "output/utils/readers.h" +#include "output/utils/writers.h" #if defined(MPI_ENABLED) #include #endif #include +#include #include #include @@ -68,6 +71,23 @@ namespace ntt { for (const auto& species : local_domain->species) { species.CheckpointDeclare(g_checkpoint_writer.io()); } + for (auto d { 0u }; d < M::Dim; ++d) { + g_checkpoint_writer.io().DefineVariable( + fmt::format("subdomain_x%d_min", d + 1), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + g_checkpoint_writer.io().DefineVariable( + fmt::format("subdomain_x%d_max", d + 1), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + g_checkpoint_writer.io().DefineVariable( + fmt::format("subdomain_nx%d", d + 1), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + } } } @@ -111,12 +131,117 @@ namespace ntt { dom_tot, dom_offset); } + for (auto d { 0u }; d < M::Dim; ++d) { + out::WriteVariable(g_checkpoint_writer.io(), + g_checkpoint_writer.writer(), + fmt::format("subdomain_x%d_min", d + 1), + local_domain->mesh.extent()[d].first, + dom_tot, + dom_offset); + out::WriteVariable(g_checkpoint_writer.io(), + g_checkpoint_writer.writer(), + fmt::format("subdomain_x%d_max", d + 1), + local_domain->mesh.extent()[d].second, + dom_tot, + dom_offset); + out::WriteVariable(g_checkpoint_writer.io(), + g_checkpoint_writer.writer(), + fmt::format("subdomain_nx%d", d + 1), + local_domain->mesh.n_active()[d], + dom_tot, + dom_offset); + } } g_checkpoint_writer.endSaving(); logger::Checkpoint("Checkpoint written", HERE); return true; } + template + void Metadomain::redecomposeFromCheckpoint( + const std::vector>& dom_ncells, + const std::vector>& dom_extents) { + + // For each dimension: collect ncells per domain-grid position, + // validate total is unchanged, then compute prefix-sum offsets. + std::vector> offset_ncells_per_dom( + g_ndomains, + std::vector(M::Dim, 0)); + + for (auto d { 0u }; d < M::Dim; ++d) { + std::vector ncells_at_pos(g_ndomains_per_dim[d], 0); + for (unsigned int idx { 0 }; idx < g_ndomains; ++idx) { + ncells_at_pos[g_domain_offsets[idx][d]] = dom_ncells[idx][d]; + } + + ncells_t total { 0 }; + for (const auto& n : ncells_at_pos) { + total += n; + } + raise::ErrorIf( + total != g_mesh.n_active()[d], + fmt::format( + "total cells in dim %d changed between checkpoint (%lu) and current (%lu); " + "changing total domain size is not supported", + d + 1, + total, + g_mesh.n_active()[d]), + HERE); + + ncells_t running { 0 }; + std::vector offset_at_pos(g_ndomains_per_dim[d], 0); + for (unsigned int nd { 0 }; nd < g_ndomains_per_dim[d]; ++nd) { + offset_at_pos[nd] = running; + running += ncells_at_pos[nd]; + } + for (unsigned int idx { 0 }; idx < g_ndomains; ++idx) { + offset_ncells_per_dom[idx][d] = offset_at_pos[g_domain_offsets[idx][d]]; + } + } + + g_subdomains.clear(); + for (unsigned int idx { 0 }; idx < g_ndomains; ++idx) { + const auto& l_offset_ndomains = g_domain_offsets[idx]; + const auto& l_ncells = dom_ncells[idx]; + const auto& l_offset_ncells = offset_ncells_per_dom[idx]; + const auto& l_extent = dom_extents[idx]; + +#if defined(MPI_ENABLED) + const auto local = ((int)idx == g_mpi_rank); + if (not local) { + g_subdomains.emplace_back(false, + idx, + l_offset_ndomains, + l_offset_ncells, + l_ncells, + l_extent, + g_metric_params, + g_species_params); + } else { + g_subdomains.emplace_back(idx, + l_offset_ndomains, + l_offset_ncells, + l_ncells, + l_extent, + g_metric_params, + g_species_params); + } + g_subdomains.back().set_mpi_rank(idx); +#else + g_subdomains.emplace_back(idx, + l_offset_ndomains, + l_offset_ncells, + l_ncells, + l_extent, + g_metric_params, + g_species_params); +#endif + } + + redefineNeighbors(); + redefineBoundaries(); + } + template void Metadomain::ContinueFromCheckpoint(adios2::ADIOS* ptr_adios, const SimulationParams& params) { @@ -140,6 +265,60 @@ namespace ntt { #endif reader.BeginStep(); + + // Phase 1: read all subdomain metadata to detect size changes + std::vector> saved_ncells(g_ndomains, + std::vector(M::Dim)); + std::vector> saved_extents(g_ndomains); + boundaries_t global_extent; + for (auto d { 0u }; d < M::Dim; ++d) { + global_extent.emplace_back(std::numeric_limits::max(), + std::numeric_limits::lowest()); + } + + bool needs_reconstruction = false; + for (unsigned int dom_idx { 0 }; dom_idx < g_ndomains; ++dom_idx) { + for (auto d { 0u }; d < M::Dim; ++d) { + real_t x_min, x_max; + out::ReadVariable(io, + reader, + fmt::format("subdomain_x%d_min", d + 1), + x_min, + dom_idx); + out::ReadVariable(io, + reader, + fmt::format("subdomain_x%d_max", d + 1), + x_max, + dom_idx); + saved_extents[dom_idx].emplace_back(x_min, x_max); + global_extent[d].first = std::min(global_extent[d].first, x_min); + global_extent[d].second = std::max(global_extent[d].second, x_max); + + ncells_t nx; + out::ReadVariable(io, + reader, + fmt::format("subdomain_nx%d", d + 1), + nx, + dom_idx); + saved_ncells[dom_idx][d] = nx; + + if (nx != subdomain_ptr(dom_idx)->mesh.n_active()[d]) { + needs_reconstruction = true; + } + } + } + + // Phase 2: update domain structure + if (needs_reconstruction) { + redecomposeFromCheckpoint(saved_ncells, saved_extents); + } else { + for (unsigned int dom_idx { 0 }; dom_idx < g_ndomains; ++dom_idx) { + subdomain_ptr(dom_idx)->mesh.set_extent(saved_extents[dom_idx]); + } + } + g_mesh.set_extent(global_extent); + + // Phase 3: read field and particle data using the (now-correct) domain layout for (const auto local_domain_idx : l_subdomain_indices()) { auto local_domain = subdomain_ptr(local_domain_idx); @@ -154,8 +333,7 @@ namespace ntt { for (auto& species : local_domain->species) { species.CheckpointRead(io, reader, ndomains(), local_domain_idx); } - - } // local subdomain loop + } reader.EndStep(); reader.Close(); @@ -177,7 +355,10 @@ namespace ntt { simtime_t) -> bool; \ template void Metadomain>::ContinueFromCheckpoint( \ adios2::ADIOS*, \ - const SimulationParams&); + const SimulationParams&); \ + template void Metadomain>::redecomposeFromCheckpoint( \ + const std::vector>&, \ + const std::vector>&); NTT_FOREACH_SPECIALIZATION(METADOMAIN_CHECKPOINTS) #undef METADOMAIN_CHECKPOINTS // NOLINTEND(bugprone-macro-parentheses) diff --git a/tests/framework/CMakeLists.txt b/tests/framework/CMakeLists.txt index ebeadc75c..d6dc295af 100644 --- a/tests/framework/CMakeLists.txt +++ b/tests/framework/CMakeLists.txt @@ -56,6 +56,16 @@ else() gen_test(comm-nompi false) endif() +if(${output}) + if(${mpi}) + gen_test(checkpoint-extents true) + gen_test(checkpoint-data true) + else() + gen_test(checkpoint-extents false) + gen_test(checkpoint-data false) + endif() +endif() + # this test is only run manually to ensure ... # ... command line args are working properly ... # ... and that the logging is done correctly diff --git a/tests/framework/checkpoint-data.cpp b/tests/framework/checkpoint-data.cpp new file mode 100644 index 000000000..3099593de --- /dev/null +++ b/tests/framework/checkpoint-data.cpp @@ -0,0 +1,244 @@ +/** + * Tests that field array sizes and particle coordinates are correctly preserved + * through a checkpoint round-trip. + * + * In MPI mode the checkpoint is written from a metadomain whose subdomain + * sizes have been changed via redecomposeFromCheckpoint (simulating load + * balancing from [32,32] to [40,24] cells), then read back into a uniformly + * decomposed metadomain ([32,32]), exercising the reconstruction path in + * ContinueFromCheckpoint. + */ + +#include "enums.h" +#include "global.h" + +#include "utils/comparators.h" +#include "utils/error.h" +#include "utils/formatting.h" + +#include "framework/domain/metadomain.h" +#include "framework/parameters/parameters.h" + +#include "metrics/minkowski.h" + +#include +#include +#include + +#if defined(MPI_ENABLED) + #include +#endif + +#include +#include +#include +#include + +using namespace ntt; +using namespace metric; + +const std::string CHECKPOINT_DIR = "test_checkpoint_data"; + +void cleanup() { + std::filesystem::remove_all(CHECKPOINT_DIR); +} + +auto main(int argc, char* argv[]) -> int { + GlobalInitialize(argc, argv); + + try { + using M = Minkowski; + + const std::vector res { 64 }; + const boundaries_t init_extent { { static_cast(0.0), + static_cast(10.0) } }; + const boundaries_t fldsbc { { FldsBC::PERIODIC, FldsBC::PERIODIC } }; + const boundaries_t prtlbc { { PrtlBC::PERIODIC, PrtlBC::PERIODIC } }; + const std::vector decomp { -1 }; + + const std::vector species_params { ParticleSpecies { + static_cast(1), + "e-", + 1.0f, + -1.0f, + static_cast(10), + timestep_t { 0 }, + timestep_t { 0 }, + ParticlePusher::BORIS, + false, + RadiativeDrag::NONE, + EmissionType::NONE, + static_cast(0), + static_cast(0) } }; + +#if !defined(MPI_ENABLED) + const unsigned int ndomains { 1 }; + adios2::ADIOS adios; + // non-MPI: single domain, sizes always match (equal-sizes path) + const std::vector saved_ncells_per_dom { 64 }; +#else + int mpi_size, mpi_rank; + MPI_Comm_size(MPI_COMM_WORLD, &mpi_size); + MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); + const unsigned int ndomains { static_cast(mpi_size) }; + adios2::ADIOS adios { MPI_COMM_WORLD }; + raise::ErrorIf(mpi_size != 2, "this test requires exactly 2 MPI ranks", HERE); + // asymmetric split: 40+24=64; exercises redecomposeFromCheckpoint + const std::vector saved_ncells_per_dom { 40, 24 }; + const real_t x_split = init_extent[0].second * + static_cast(saved_ncells_per_dom[0]) / + static_cast(64); +#endif + + const auto checkpoint_step = timestep_t { 2 }; + const int prtl_i1_0 { 3 }; + const real_t prtl_dx1_0 { static_cast(0.25) }; + const int prtl_i1_1 { 7 }; + const real_t prtl_dx1_1 { static_cast(0.75) }; + const npart_t npart_placed { 2 }; + + // ── Write phase ─────────────────────────────────────────────────────────── + { + Metadomain md { + ndomains, decomp, res, init_extent, fldsbc, prtlbc, {}, species_params + }; + +#if defined(MPI_ENABLED) + // Simulate load balancing: [32, 32] → [40, 24] + md.redecomposeFromCheckpoint( + { { saved_ncells_per_dom[0] }, { saved_ncells_per_dom[1] } }, + { { { init_extent[0].first, x_split } }, + { { x_split, init_extent[0].second } } }); +#endif + + auto* local = md.subdomain_ptr(md.l_subdomain_indices()[0]); + const auto n_all = local->mesh.n_all()[0]; + + // set em::ex1(i) = real_t(i) over all cells including ghost cells + { + auto em_h = Kokkos::create_mirror_view(local->fields.em); + for (std::size_t i { 0 }; i < n_all; ++i) { + em_h(i, em::ex1) = static_cast(i); + } + Kokkos::deep_copy(local->fields.em, em_h); + } + + // place 2 particles at known positions (safe within the 24-cell domain) + auto& sp = local->species[0]; + sp.set_npart(npart_placed); + { + auto i1_h = Kokkos::create_mirror_view(sp.i1); + auto dx1_h = Kokkos::create_mirror_view(sp.dx1); + i1_h(0) = prtl_i1_0; + dx1_h(0) = static_cast(prtl_dx1_0); + i1_h(1) = prtl_i1_1; + dx1_h(1) = static_cast(prtl_dx1_1); + Kokkos::deep_copy(sp.i1, i1_h); + Kokkos::deep_copy(sp.dx1, dx1_h); + } + + SimulationParams params; + params.set("checkpoint.write_path", std::string { CHECKPOINT_DIR }); + params.set("checkpoint.interval", timestep_t { 1 }); + params.set("checkpoint.interval_time", simtime_t { -1.0 }); + params.set("checkpoint.keep", 5); + params.set("checkpoint.walltime", std::string { "" }); + params.setRawData(toml::value { toml::table {} }); + + md.InitCheckpointWriter(&adios, params); + + const auto wrote = md.WriteCheckpoint( + params, + checkpoint_step, + checkpoint_step, + simtime_t { 0.0 }, + simtime_t { 0.0 }); + raise::ErrorIf(not wrote, "checkpoint was not written", HERE); + } + + // ── Read phase ──────────────────────────────────────────────────────────── + { + // fresh metadomain with the original (uniform) decomposition + // MPI: [32, 32] — mismatch with saved [40, 24] → reconstruction triggered + Metadomain md2 { + ndomains, decomp, res, init_extent, fldsbc, prtlbc, {}, species_params + }; + + SimulationParams params2; + params2.set("checkpoint.read_path", std::string { CHECKPOINT_DIR }); + params2.set("checkpoint.start_step", checkpoint_step); + + md2.ContinueFromCheckpoint(&adios, params2); + + auto* local = md2.subdomain_ptr(md2.l_subdomain_indices()[0]); + const auto lidx = md2.l_subdomain_indices()[0]; + const auto exp_ncells = saved_ncells_per_dom[lidx]; + const auto exp_n_all = exp_ncells + 2 * N_GHOSTS; + + // 1. field array size reflects saved (possibly reconstructed) ncells + raise::ErrorIf( + local->fields.em.extent(0) != exp_n_all, + fmt::format("field em extent: got %lu, expected %lu", + local->fields.em.extent(0), + exp_n_all), + HERE); + + // 2. field values survive the round-trip + { + auto em_h = Kokkos::create_mirror_view(local->fields.em); + Kokkos::deep_copy(em_h, local->fields.em); + for (std::size_t i { 0 }; i < exp_n_all; ++i) { + raise::ErrorIf( + not cmp::AlmostEqual(em_h(i, em::ex1), static_cast(i)), + fmt::format("em::ex1 mismatch at i=%lu: got %f, expected %f", + i, + static_cast(em_h(i, em::ex1)), + static_cast(i)), + HERE); + } + } + + // 3. particle count and coordinates survive the round-trip + auto& sp = local->species[0]; + raise::ErrorIf( + sp.npart() != npart_placed, + fmt::format("particle count: got %lu, expected %lu", + sp.npart(), + npart_placed), + HERE); + { + auto i1_h = Kokkos::create_mirror_view(sp.i1); + auto dx1_h = Kokkos::create_mirror_view(sp.dx1); + Kokkos::deep_copy(i1_h, sp.i1); + Kokkos::deep_copy(dx1_h, sp.dx1); + + raise::ErrorIf( + i1_h(0) != prtl_i1_0, + fmt::format("particle 0 i1: got %d, expected %d", i1_h(0), prtl_i1_0), + HERE); + raise::ErrorIf( + not cmp::AlmostEqual(static_cast(dx1_h(0)), prtl_dx1_0), + "particle 0 dx1 mismatch", + HERE); + raise::ErrorIf( + i1_h(1) != prtl_i1_1, + fmt::format("particle 1 i1: got %d, expected %d", i1_h(1), prtl_i1_1), + HERE); + raise::ErrorIf( + not cmp::AlmostEqual(static_cast(dx1_h(1)), prtl_dx1_1), + "particle 1 dx1 mismatch", + HERE); + } + } + + } catch (const std::exception& e) { + std::cerr << e.what() << '\n'; + cleanup(); + GlobalFinalize(); + return 1; + } + + cleanup(); + GlobalFinalize(); + return 0; +} diff --git a/tests/framework/checkpoint-extents.cpp b/tests/framework/checkpoint-extents.cpp new file mode 100644 index 000000000..ecf6e5793 --- /dev/null +++ b/tests/framework/checkpoint-extents.cpp @@ -0,0 +1,163 @@ +/** + * Tests that per-subdomain extents and global mesh extents are correctly saved + * to and restored from a checkpoint, as needed by the moving window feature. + */ + +#include "enums.h" +#include "global.h" + +#include "utils/comparators.h" +#include "utils/error.h" + +#include "framework/domain/metadomain.h" +#include "framework/parameters/parameters.h" + +#include "metrics/minkowski.h" + +#include +#include +#include + +#if defined(MPI_ENABLED) + #include +#endif + +#include +#include +#include +#include + +using namespace ntt; +using namespace metric; + +const std::string CHECKPOINT_DIR = "test_checkpoint_extents"; + +void cleanup() { + std::filesystem::remove_all(CHECKPOINT_DIR); +} + +auto extents_match(const boundaries_t& a, + const boundaries_t& b) -> bool { + if (a.size() != b.size()) { + return false; + } + for (std::size_t d { 0 }; d < a.size(); ++d) { + if (not cmp::AlmostEqual(a[d].first, b[d].first) or + not cmp::AlmostEqual(a[d].second, b[d].second)) { + return false; + } + } + return true; +} + +auto main(int argc, char* argv[]) -> int { + GlobalInitialize(argc, argv); + + try { + using M = Minkowski; + + const std::vector res { 64 }; + const boundaries_t init_extent { { static_cast(0.0), + static_cast(10.0) } }; + const boundaries_t fldsbc { + { FldsBC::PERIODIC, FldsBC::PERIODIC } + }; + const boundaries_t prtlbc { + { PrtlBC::PERIODIC, PrtlBC::PERIODIC } + }; + const std::vector decomp { -1 }; + +#if !defined(MPI_ENABLED) + const unsigned int ndomains { 1 }; + adios2::ADIOS adios; +#else + int mpi_size; + MPI_Comm_size(MPI_COMM_WORLD, &mpi_size); + const unsigned int ndomains { static_cast(mpi_size) }; + adios2::ADIOS adios { MPI_COMM_WORLD }; +#endif + + const auto checkpoint_step = timestep_t { 2 }; + const int n_shift = 4; + + // reference extents captured after the window shift + boundaries_t expected_global_extent; + std::vector> expected_subdomain_extents(ndomains); + + { + Metadomain md { + ndomains, decomp, res, init_extent, fldsbc, prtlbc, {}, {} + }; + + SimulationParams params; + params.set("checkpoint.write_path", std::string { CHECKPOINT_DIR }); + params.set("checkpoint.interval", timestep_t { 1 }); + params.set("checkpoint.interval_time", simtime_t { -1.0 }); + params.set("checkpoint.keep", 5); + params.set("checkpoint.walltime", std::string { "" }); + params.setRawData(toml::value { toml::table {} }); + + md.InitCheckpointWriter(&adios, params); + + // shift the window and capture the resulting extents + md.ShiftByCells(n_shift, in::x1); + expected_global_extent = md.mesh().extent(); + for (unsigned int idx { 0 }; idx < ndomains; ++idx) { + expected_subdomain_extents[idx] = md.subdomain(idx).mesh.extent(); + } + + // write checkpoint at step 2 (finished_step must be > 1 to be saved) + const auto wrote = md.WriteCheckpoint( + params, + checkpoint_step, + checkpoint_step, + simtime_t { 0.0 }, + simtime_t { 0.0 }); + raise::ErrorIf(not wrote, "checkpoint was not written", HERE); + } + + { + // construct a fresh metadomain from the original (pre-shift) extent + Metadomain md2 { + ndomains, decomp, res, init_extent, fldsbc, prtlbc, {}, {} + }; + + // sanity: verify global extent is the original one before reading + raise::ErrorIf( + extents_match(md2.mesh().extent(), expected_global_extent), + "global extent should differ from shifted extent before checkpoint read", + HERE); + + SimulationParams params2; + params2.set("checkpoint.read_path", std::string { CHECKPOINT_DIR }); + params2.set("checkpoint.start_step", checkpoint_step); + + md2.ContinueFromCheckpoint(&adios, params2); + + // global mesh extent must match the shifted one + raise::ErrorIf( + not extents_match(md2.mesh().extent(), expected_global_extent), + "global mesh extent mismatch after checkpoint read", + HERE); + + // per-subdomain extents must also match + for (unsigned int idx { 0 }; idx < ndomains; ++idx) { + raise::ErrorIf( + not extents_match(md2.subdomain(idx).mesh.extent(), + expected_subdomain_extents[idx]), + "subdomain extent mismatch after checkpoint read", + HERE); + } + } + + } catch (const std::exception& e) { + std::cerr << e.what() << '\n'; + cleanup(); + GlobalFinalize(); + return 1; + } + + cleanup(); + GlobalFinalize(); + return 0; +} From 1f552ac899be589d4e77cc0de77ed934f84e524c Mon Sep 17 00:00:00 2001 From: hayk Date: Sun, 26 Apr 2026 03:31:05 -0400 Subject: [PATCH 26/30] upd claude md --- CLAUDE.md | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index b96326906..dd36d3fa3 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -49,6 +49,7 @@ entity │ └── entity.cpp # main entry-point ├── tests # unit tests for all submodules ├── .clang-format # code formatting guidelines for clang-format +├── .clang-tidy # configurations for the clang-tidy ├── .gitattributes ├── .gitignore ├── .gitmodules @@ -59,9 +60,7 @@ entity ├── CODE_OF_CONDUCT.md ├── LICENSE ├── README.md -├── conda-entity-nompi.sh ├── dependencies.py # deployment scripts on various machines -├── docker-compose.yml └── input.example.toml # most complete toml file with all possible input options ``` From c27b833c055ea8464434af04d75d9122e5943d09 Mon Sep 17 00:00:00 2001 From: hayk Date: Sun, 26 Apr 2026 18:57:50 -0400 Subject: [PATCH 27/30] media in examples --- pgens/examples/moving_window/moving_window.mp4 | 3 +++ pgens/examples/piston/piston.mp4 | 3 +++ 2 files changed, 6 insertions(+) create mode 100644 pgens/examples/moving_window/moving_window.mp4 create mode 100644 pgens/examples/piston/piston.mp4 diff --git a/pgens/examples/moving_window/moving_window.mp4 b/pgens/examples/moving_window/moving_window.mp4 new file mode 100644 index 000000000..373af2e4a --- /dev/null +++ b/pgens/examples/moving_window/moving_window.mp4 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:280520498f8064d5e4ef406accc9f3c7919f9f564fa76881f481766e70f99928 +size 3041327 diff --git a/pgens/examples/piston/piston.mp4 b/pgens/examples/piston/piston.mp4 new file mode 100644 index 000000000..a4fe967fc --- /dev/null +++ b/pgens/examples/piston/piston.mp4 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e90c9191bfb2c53e0be215ec06fceff673133e1381b0da436f742256619e6022 +size 10899641 From c73ba4199e116b5f53931517daa5630af40bda16 Mon Sep 17 00:00:00 2001 From: hayk Date: Sun, 26 Apr 2026 19:04:57 -0400 Subject: [PATCH 28/30] Fix typos and improve README examples --- pgens/examples/README.md | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/pgens/examples/README.md b/pgens/examples/README.md index 93be01d40..9309c4713 100644 --- a/pgens/examples/README.md +++ b/pgens/examples/README.md @@ -6,11 +6,11 @@ Problem generators in this directory are just examples demonstrating how to use image -- `custom_spatial_distribution`: example of using the non-uniform plasma injector with a custom spatial distirbution +- `custom_spatial_distribution`: example of using the non-uniform plasma injector with a custom spatial distirbution. https://github.com/user-attachments/assets/bea0c290-e7e4-4ec7-b360-68ce76beab5b -- `match_fix_field_boundaries`: example of setting matching and/or fixed (coordinate-independent) field boundaries for the electromagnetic fields +- `match_fix_field_boundaries`: example of setting matching and/or fixed (coordinate-independent) field boundaries for the electromagnetic fields. https://github.com/user-attachments/assets/a1b9ea22-34ce-474f-a9b0-49789a2e52b3 @@ -22,13 +22,23 @@ Problem generators in this directory are just examples demonstrating how to use image -- `atmosphere`: setting up a gravitationally bound "atmosphere" with a constant particle replenisher (the plot also highlights the importance of having a constant replenisher and a gravity force acting on the particles; both are enabled by default when using the `ATMOSPHERE` particle boundary conditions) +- `atmosphere`: setting up a gravitationally bound "atmosphere" with a constant particle replenisher (the plot also highlights the importance of having a constant replenisher and a gravity force acting on the particles; both are enabled by default when using the `ATMOSPHERE` particle boundary conditions). atmosphere -- `replenish_injector`: demonstration of how to use the replenish injector to periodically inject new plasma to a target density (both uniform and non-uniform) in `CustomPostStep` +- `replenish_injector`: demonstration of how to use the replenish injector to periodically inject new plasma to a target density (both uniform and non-uniform) in `CustomPostStep`. https://github.com/user-attachments/assets/4955f05e-4795-439f-a4a5-8196e42e987b https://github.com/user-attachments/assets/ebc303ed-7b21-4f97-9822-d81124fdf962 +- `piston`: implementing a moving piston which continuously pushes the plasma. + + https://github.com/user-attachments/assets/7d2dfc39-5de4-49cb-9d55-be8700cbcaee + +- `moving_window`: demonstrates the domain being updated to follow certain features in your simulation in one of the direction (jerkiness is due to updates taking place every `n`-th timestep, if you limit plotting to a specific physical interval, this goes away). + + https://github.com/user-attachments/assets/00edb648-af80-4d37-b1c2-c49872227d47 + + + From 5e2e660082ac7fbab236c27493e30d7488a44fba Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Mon, 27 Apr 2026 02:46:56 +0000 Subject: [PATCH 29/30] bugfix in piston --- src/archetypes/piston.h | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/src/archetypes/piston.h b/src/archetypes/piston.h index 9f3f031bf..f31dbdffa 100644 --- a/src/archetypes/piston.h +++ b/src/archetypes/piston.h @@ -130,23 +130,20 @@ namespace arch { TWO * piston_v * gamma_p); const real_t remaining_dt_inv_energy { - massive - ? (remaining_dt / - math::sqrt( - ONE + U2GAMMA(particles.ux1(p), particles.ux2(p), particles.ux3(p)))) - : (remaining_dt / math::sqrt(SQR(particles.ux1(p)) + SQR(particles.ux2(p)) + - SQR(particles.ux3(p)))) + massive ? (remaining_dt / + U2GAMMA(particles.ux1(p), particles.ux2(p), particles.ux3(p))) + : (remaining_dt / + NORM(particles.ux1(p), particles.ux2(p), particles.ux3(p))) }; // define piston integer and fractional coordinate - int i_w_coll = i_w; prtldx_t dx_w_coll = dx_w + metric.template transform<1, Idx::XYZ, Idx::U>( {}, piston_v) * dt_to_piston; - i_w_coll += static_cast(dx_w_coll >= ONE) - - static_cast(dx_w_coll < ZERO); + i_w_coll += static_cast(dx_w_coll >= ONE) - + static_cast(dx_w_coll < ZERO); dx_w_coll -= static_cast(dx_w_coll >= ONE); dx_w_coll += static_cast(dx_w_coll < ZERO); @@ -156,8 +153,8 @@ namespace arch { remaining_dt_inv_energy + dx_w_coll; - particles.i1(p) += static_cast(particles.dx1(p) >= ONE) - - static_cast(particles.dx1(p) < ZERO); + particles.i1(p) += static_cast(particles.dx1(p) >= ONE) - + static_cast(particles.dx1(p) < ZERO); particles.dx1(p) -= static_cast(particles.dx1(p) >= ONE); particles.dx1(p) += static_cast(particles.dx1(p) < ZERO); } From 5bdba11462a36cd02111f647e16f2aa3044b1dac Mon Sep 17 00:00:00 2001 From: haykh Date: Mon, 27 Apr 2026 13:20:15 -0400 Subject: [PATCH 30/30] RUNPGENS + RUNTESTS