From e9c927b7240b1c6343ba5badfd97c06e93ed0c80 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Fri, 3 Apr 2026 03:36:02 +0000 Subject: [PATCH 1/7] 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 86839f0c..d5e24402 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 00000000..24174eff --- /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 2/7] 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 d5e24402..20098b7d 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 24174eff..eff0c172 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 3/7] 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 117ea905..6f3397c1 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 4/7] 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 eff0c172..50f78cee 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 5/7] 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 6f3397c1..d066e832 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 6/7] 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 50f78cee..49865b26 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 7/7] 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 d066e832..22f373b4 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); }