diff --git a/pgens/examples/custom_particle_update/particle_update.toml b/pgens/examples/custom_particle_update/particle_update.toml new file mode 100644 index 00000000..017ba82b --- /dev/null +++ b/pgens/examples/custom_particle_update/particle_update.toml @@ -0,0 +1,82 @@ +[simulation] + name = "_reflect" + engine = "srpic" + runtime = 100.0 + +[grid] + resolution = [256, 256] + extent = [[0.0, 1.0], [0.0, 1.0]] + + [grid.metric] + metric = "minkowski" + + [grid.boundaries] + fields = [["CONDUCTOR", "CONDUCTOR"], ["PERIODIC"]] + particles = [["PERIODIC"], ["PERIODIC"]] + +[scales] + larmor0 = 1.0 + skindepth0 = 1.0 + +[algorithms] + current_filters = 8 + + [algorithms.timestep] + CFL = 0.7071 + +[particles] + ppc0 = 10.0 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 1e2 + + [[particles.species]] + label = "e+" + mass = 1.0 + charge = 1.0 + maxnpart = 1e2 + +[setup] + x1_e = [0.9] + x2_e = [0.8] + x3_e = [0.0] + ux1_e = [3.162] + ux2_e = [2.0] + ux3_e = [0.0] + + x1_i = [0.1] + x2_i = [0.2] + x3_i = [0.0] + ux1_i = [-3.162] + ux2_i = [-2.0] + ux3_i = [0.0] + + temperature = 0.04 + temperature_gradient = 10.0 + +[output] + format = "BPFile" + interval_time = 0.04 + + [output.fields] + quantities = ["E"] + + [output.particles] + enable = true + stride = 1 + + [output.spectra] + enable = false + + [output.debug] + ghosts = false + +[checkpoint] + keep = 1 + +[diagnostics] + log_level = "WARNING" + blocking_timers = true diff --git a/pgens/examples/custom_particle_update/pgen.hpp b/pgens/examples/custom_particle_update/pgen.hpp new file mode 100644 index 00000000..9d7113d4 --- /dev/null +++ b/pgens/examples/custom_particle_update/pgen.hpp @@ -0,0 +1,243 @@ +#ifndef PROBLEM_GENERATOR_H +#define PROBLEM_GENERATOR_H + +#include "enums.h" +#include "global.h" + +#include "arch/traits.h" + +#include "archetypes/energy_dist.h" +#include "archetypes/particle_injector.h" +#include "archetypes/problem_generator.h" +#include "framework/domain/domain.h" +#include "framework/domain/metadomain.h" + +#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 PGen : public arch::ProblemGenerator { + + // compatibility traits for the problem generator + 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 + }; + + // for easy access to variables in the child class + using arch::ProblemGenerator::D; + using arch::ProblemGenerator::C; + using arch::ProblemGenerator::params; + + const Metadomain& global_domain; + + const real_t temperature, temperature_gradient; + + inline PGen(const SimulationParams& p, const Metadomain& global_domain) + : arch::ProblemGenerator { p } + , global_domain { global_domain } + , temperature { params.template get("setup.temperature", 0.0) } + , temperature_gradient { + params.template get("setup.temperature_gradient", 0.0) + } {} + + inline void InitPrtls(Domain& local_domain) { + const auto empty = std::vector {}; + const auto x1_e = params.template get>("setup.x1_e", + empty); + const auto x2_e = params.template get>("setup.x2_e", + empty); + const auto x3_e = params.template get>("setup.x3_e", + empty); + const auto phi_e = params.template get>("setup.phi_e", + empty); + const auto ux1_e = params.template get>("setup.ux1_e", + empty); + const auto ux2_e = params.template get>("setup.ux2_e", + empty); + const auto ux3_e = params.template get>("setup.ux3_e", + empty); + + const auto x1_i = params.template get>("setup.x1_i", + empty); + const auto x2_i = params.template get>("setup.x2_i", + empty); + const auto x3_i = params.template get>("setup.x3_i", + empty); + const auto phi_i = params.template get>("setup.phi_i", + empty); + const auto ux1_i = params.template get>("setup.ux1_i", + empty); + const auto ux2_i = params.template get>("setup.ux2_i", + empty); + const auto ux3_i = params.template get>("setup.ux3_i", + empty); + std::map> data_e { + { "x1", x1_e }, + { "x2", x2_e }, + { "ux1", ux1_e }, + { "ux2", ux2_e }, + { "ux3", ux3_e } + }; + std::map> data_i { + { "x1", x1_i }, + { "x2", x2_i }, + { "ux1", ux1_i }, + { "ux2", ux2_i }, + { "ux3", ux3_i } + }; + if constexpr (M::CoordType == Coord::Cart or D == Dim::_3D) { + data_e["x3"] = x3_e; + data_i["x3"] = x3_i; + } else if constexpr (D == Dim::_2D) { + data_e["phi"] = phi_e; + data_i["phi"] = phi_i; + } + + arch::InjectGlobally(global_domain, local_domain, (spidx_t)1, data_e); + arch::InjectGlobally(global_domain, local_domain, (spidx_t)2, data_i); + } + + auto FixFieldsConst(const bc_in&, const em&) const -> std::pair { + return { ZERO, false }; + } + + struct CustomPrtlUpdate { + random_number_pool_t pool; + real_t temp_cold, temp_hot; + real_t xmin, xmax; + + template + Inline void operator()(index_t p, Coord& xp, PusherKernel& pusher) const { + + const auto x_Cd = static_cast(pusher.i1(p)) + + static_cast(pusher.dx1(p)); + const auto x_Ph = pusher.metric.template convert<1, Crd::Cd, Crd::XYZ>(x_Cd); + vec_t v { ZERO }; + const coord_t x_dummy { ZERO }; + + // step 1: calculate the particle 3 velocity + const real_t gamma_p = math::sqrt( + ONE + SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p))); + + const real_t beta_x_p = math::abs(pusher.ux1(p)) / gamma_p; + const real_t xp_prev = i_di_to_Xi(pusher.i1_prev(p), pusher.dx1_prev(p)); + + // Reflecting boundary that resamples velocity + if (x_Ph < xmin) { + arch::JuttnerSinge(v, temp_cold, pool); + + // calculate the time for the particle to reach the wall + const int delta_i1_to_wall = pusher.i1_prev(p); + const prtldx_t delta_dx1_to_wall = pusher.dx1_prev(p); + const real_t dx_to_wall = i_di_to_Xi(delta_i1_to_wall, delta_dx1_to_wall); + const real_t dt_to_wall = dx_to_wall / + pusher.metric.template transform<1, Idx::XYZ, Idx::U>(x_dummy, + beta_x_p); + + // update the velocity to the post-collision value (reflection with new speed sampled from Juttner distribution) + pusher.ux1(p) = math::abs(v[0]); + pusher.ux2(p) = v[1]; + pusher.ux3(p) = v[2]; + + // calculate remaining time after the collision + const real_t remaining_dt = pusher.dt - dt_to_wall; + const real_t remaining_dt_inv_energy = remaining_dt / + math::sqrt(ONE + + SQR(pusher.ux1(p)) + + SQR(pusher.ux2(p)) + + SQR(pusher.ux3(p))); + + // update the position to the post-collision value (reflection with new speed sampled from Juttner distribution) + pusher.i1(p) = 0; + pusher.dx1(p) = pusher.metric.template transform<1, Idx::XYZ, Idx::U>( + x_dummy, + pusher.ux1(p)) * + remaining_dt_inv_energy; + + // Re-sync coordinate variable xp for subsequent boundary procedures + xp[0] = static_cast(pusher.i1(p)) + + static_cast(pusher.dx1(p)); + + } else if (x_Ph > xmax) { + arch::JuttnerSinge(v, temp_hot, pool); + + // step 2: calculate the time for the particle to reach the piston + const int delta_i1_to_wall = pusher.ni1 - 1 - pusher.i1_prev(p); + const prtldx_t delta_dx1_to_wall = ONE - pusher.dx1_prev(p); + const real_t dx_to_wall = i_di_to_Xi(delta_i1_to_wall, delta_dx1_to_wall); + const real_t dt_to_wall = dx_to_wall / + pusher.metric.template transform<1, Idx::XYZ, Idx::U>(x_dummy, + beta_x_p); + + // update the velocity to the post-collision value (reflection with new speed sampled from Juttner distribution) + pusher.ux1(p) = -math::abs(v[0]); + pusher.ux2(p) = v[1]; + pusher.ux3(p) = v[2]; + + // step 3: calculate remaining time after the collision + const real_t remaining_dt = pusher.dt - dt_to_wall; + const real_t remaining_dt_inv_energy = remaining_dt / + math::sqrt(ONE + + SQR(pusher.ux1(p)) + + SQR(pusher.ux2(p)) + + SQR(pusher.ux3(p))); + + // update the position to the post-collision value (reflection with new speed sampled from Juttner distribution) + pusher.i1(p) = pusher.ni1 - 2; + pusher.dx1(p) = ONE - + pusher.metric.template transform<1, Idx::XYZ, Idx::U>( + x_dummy, + math::abs(pusher.ux1(p))) * + remaining_dt_inv_energy; + + // Re-sync coordinate variable xp for subsequent boundary procedures + xp[0] = static_cast(pusher.i1(p)) + + static_cast(pusher.dx1(p)); + } + } + }; + + template + auto CustomParticleUpdate(simtime_t time, spidx_t sp, D& domain) const { + + return CustomPrtlUpdate { + domain.random_pool(), + temperature / domain.species[sp - 1].mass(), // sp is 1-indexed + temperature_gradient * temperature / domain.species[sp - 1].mass(), + global_domain.mesh().extent(in::x1).first, // xmin + global_domain.mesh().extent(in::x1).second // xmax + }; + } + }; + +} // namespace user + +#undef from_Xi_to_i_di +#undef from_Xi_to_i +#undef i_di_to_Xi + +#endif // PROBLEM_GENERATOR_H \ No newline at end of file diff --git a/pgens/piston_shock/pgen.hpp b/pgens/piston_shock/pgen.hpp new file mode 100644 index 00000000..745f751d --- /dev/null +++ b/pgens/piston_shock/pgen.hpp @@ -0,0 +1,148 @@ +#ifndef PROBLEM_GENERATOR_H +#define PROBLEM_GENERATOR_H + +#include "enums.h" +#include "global.h" + +#include "arch/traits.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 "archetypes/utils.h" +#include "framework/domain/domain.h" +#include "framework/domain/metadomain.h" + +#include + +namespace user { + using namespace ntt; + + template + struct PGen : public arch::ProblemGenerator { + + // compatibility traits for the problem generator + 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 + }; + + // for easy access to variables in the child class + using arch::ProblemGenerator::D; + using arch::ProblemGenerator::C; + using arch::ProblemGenerator::params; + + const Metadomain& global_domain; + + // domain properties + const real_t global_xmin, global_xmax; + + // plasma + const real_t temperature, temperature_ratio; + // piston properties + const real_t piston_velocity, piston_initial_position; + + inline PGen(const SimulationParams& p, const 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 } + , 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); + } + } + + // 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, + local_domain, + ONE, + temperatures, + { 1, 2 }, + drifts, + false, + 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 { + + real_t piston_position_use; + + // make sure piston has not reached the right wall + if (piston_position_current < xg_max) { + piston_position_use = piston_position_current; + } else { + piston_position_use = xg_max; + } + + 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 }; + }; + }; + +} // 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 new file mode 100644 index 00000000..6a81d1ab --- /dev/null +++ b/pgens/piston_shock/piston_shock.toml @@ -0,0 +1,75 @@ +[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 new file mode 100644 index 00000000..8cde8699 --- /dev/null +++ b/pgens/piston_shock/shock.toml @@ -0,0 +1,74 @@ +[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 diff --git a/src/archetypes/piston.h b/src/archetypes/piston.h new file mode 100644 index 00000000..171d1728 --- /dev/null +++ b/src/archetypes/piston.h @@ -0,0 +1,183 @@ +/** + * @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_PISTON_H +#define ARCHETYPES_PISTON_H + +#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 + +/* -------------------------------------------------------------------------- */ +/* 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 arch { + + /** + * @brief Updates particle position and velocity if it reflects off a moiving + * piston, called to correct patticle position + * + * @param p Index of particle + * @param pusher Particle pusher engine for particle update + * @param piston_position Position of the piston at the start of timestep in + * global coordinates + * @param piston_v Velocity of piston at current timestep + */ + template + Inline void PistonUpdate(const index_t p, + const PusherKernel& pusher, + const real_t piston_position, + const real_t piston_v, + const bool massive) { + + coord_t xp_Cd { ZERO }; + pusher.getParticlePosition(p, xp_Cd); + + // step 1: calculate the particle 3 velocity + const real_t gamma_p { + massive + ? (math::sqrt( + ONE + SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p)))) + : (math::sqrt(SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p)))) + }; + + const real_t beta_x_p = pusher.ux1(p) / gamma_p; + const real_t xp_prev = i_di_to_Xi(pusher.i1_prev(p), pusher.dx1_prev(p)); + + const real_t piston_position_local = + pusher.metric.template convert<1, Crd::Ph, Crd::Cd>(piston_position); + int i_w; + prtldx_t dx_w; + + from_Xi_to_i_di(piston_position_local, i_w, dx_w); + + const real_t piston_gamma = ONE / math::sqrt(ONE - SQR(piston_v)); + + // step 2: calculate the time for the particle to reach the piston + const int delta_i1_to_piston = (i_w - pusher.i1_prev(p)); + const prtldx_t delta_dx1_to_piston = (dx_w - pusher.dx1_prev(p)); + const real_t dx_to_piston = i_di_to_Xi(delta_i1_to_piston, delta_dx1_to_piston); + + const real_t dt_to_piston = dx_to_piston / + pusher.metric.template transform<1, Idx::XYZ, Idx::U>( + xp_Cd, + beta_x_p - piston_v); + // step 3: calculate remaining time after the collision + const real_t remaining_dt = pusher.dt - dt_to_piston; + // step 4: update the particle's velocity and position after the collision + + pusher.ux1(p) = -SQR(piston_gamma) * ((ONE + SQR(piston_v)) * pusher.ux1(p) - + TWO * piston_v * gamma_p); + + const real_t remaining_dt_inv_energy { + massive + ? (remaining_dt / math::sqrt(ONE + SQR(pusher.ux1(p)) + + SQR(pusher.ux2(p)) + SQR(pusher.ux3(p)))) + : (remaining_dt / math::sqrt(SQR(pusher.ux1(p)) + SQR(pusher.ux2(p)) + + SQR(pusher.ux3(p)))) + }; + // define piston integer and fractional coordinate + + int i_w_coll = i_w; + prtldx_t dx_w_coll = dx_w + + pusher.metric.template transform<1, Idx::XYZ, Idx::U>( + xp_Cd, + piston_v) * + dt_to_piston; + + i_w_coll += static_cast(dx_w_coll >= ONE) - + static_cast(dx_w_coll < ZERO); + dx_w_coll -= (dx_w_coll >= ONE); + dx_w_coll += (dx_w_coll < ZERO); + + pusher.i1(p) = i_w_coll; + pusher.dx1(p) = pusher.metric.template transform<1, Idx::XYZ, Idx::U>( + xp_Cd, + pusher.ux1(p)) * + remaining_dt_inv_energy + + dx_w_coll; + + pusher.i1(p) += static_cast(pusher.dx1(p) >= ONE) - + static_cast(pusher.dx1(p) < ZERO); + pusher.dx1(p) -= (pusher.dx1(p) >= ONE); + pusher.dx1(p) += (pusher.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 pusher Particle pusher engine for particle update + * @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(const index_t p, + const PusherKernel& pusher, + const real_t piston_position, + const real_t piston_v, + const bool is_left) { + const real_t x1_Cd = i_di_to_Xi(pusher.i1(p), pusher.dx1(p)); + coord_t xp_Cd { ZERO }; + pusher.getParticlePosition(p, xp_Cd); + // 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 - + pusher.metric.template transform<1, Idx::XYZ, Idx::U>(xp_Cd, piston_v) * + pusher.dt; + const real_t x1_Ph_wallmove = pusher.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 +#undef from_Xi_to_i +#undef i_di_to_Xi + +#endif // ARCHETYPES_UTILS_H diff --git a/src/archetypes/traits.h b/src/archetypes/traits.h index f453c7fc..45ba642b 100644 --- a/src/archetypes/traits.h +++ b/src/archetypes/traits.h @@ -85,6 +85,14 @@ namespace arch { pgen.EmissionPolicy(time, sp, domain); }; + template + concept HasCustomPrtlUpdate = requires(const PG& pgen, + simtime_t time, + spidx_t sp, + D& domain) { + pgen.CustomParticleUpdate(time, sp, domain); + }; + template concept HasInitPrtls = requires(PG& pgen, D& domain) { { pgen.InitPrtls(domain) } -> std::same_as; diff --git a/src/engines/srpic/particle_pusher.h b/src/engines/srpic/particle_pusher.h index 47189e98..a63cd413 100644 --- a/src/engines/srpic/particle_pusher.h +++ b/src/engines/srpic/particle_pusher.h @@ -37,18 +37,31 @@ namespace ntt { const M& metric, const PG& pgen, const F& external_fields) { + auto get_custom_prtl_update = [&]() { + if constexpr ( + arch::traits::pgen::HasCustomPrtlUpdate) { + return pgen.CustomParticleUpdate(pusher_params.time, + pusher_params.species_index, + domain); + } else { + return kernel::sr::NoCustomPrtlUpdate_t {}; + } + }; + const auto custom_prtl_update = get_custom_prtl_update(); + if (emission_policy_flag == EmissionType::NONE) { const auto no_emission = kernel::NoEmissionPolicy_t {}; Kokkos::parallel_for( "ParticlePusher", range, - kernel::sr::Pusher_kernel( + kernel::sr::Pusher_kernel( pusher_params, pusher_arrays, EB, metric, external_fields, - no_emission)); + no_emission, + custom_prtl_update)); } else if (emission_policy_flag == EmissionType::SYNCHROTRON) { const auto photon_species = params.get( "radiation.emission.synchrotron.photon_species"); @@ -78,13 +91,14 @@ namespace ntt { Kokkos::parallel_for( "ParticlePusher", range, - kernel::sr::Pusher_kernel( + kernel::sr::Pusher_kernel( pusher_params, pusher_arrays, EB, metric, external_fields, - emission_policy)); + emission_policy, + custom_prtl_update)); const auto n_inj = emission_policy.numbers_injected(); raise::ErrorIf(n_inj.size() != 1, "Synchrotron emission should only inject one species", @@ -121,13 +135,14 @@ namespace ntt { Kokkos::parallel_for( "ParticlePusher", range, - kernel::sr::Pusher_kernel( + kernel::sr::Pusher_kernel( pusher_params, pusher_arrays, EB, metric, external_fields, - emission_policy)); + emission_policy, + custom_prtl_update)); const auto n_inj = emission_policy.numbers_injected(); raise::ErrorIf(n_inj.size() != 1, "Compton emission should only inject one species", @@ -148,13 +163,14 @@ namespace ntt { Kokkos::parallel_for( "ParticlePusher", range, - kernel::sr::Pusher_kernel( + kernel::sr::Pusher_kernel( pusher_params, pusher_arrays, EB, metric, external_fields, - emission_policy)); + emission_policy, + custom_prtl_update)); const auto emitted_species = emission_policy.emitted_species_indices(); const auto n_inj = emission_policy.numbers_injected(); raise::ErrorIf(emitted_species.size() != n_inj.size(), diff --git a/src/kernels/particle_pusher_sr.hpp b/src/kernels/particle_pusher_sr.hpp index 6706842e..31055eff 100644 --- a/src/kernels/particle_pusher_sr.hpp +++ b/src/kernels/particle_pusher_sr.hpp @@ -101,11 +101,18 @@ namespace kernel::sr { array_t tag; }; + template + struct NoCustomPrtlUpdate_t { + template + Inline void operator()(index_t, coord_t&, const PusherKernel&) const { + } + }; + /** * @tparam M Metric * @tparam F Additional force */ - template + template requires metric::traits::HasD && metric::traits::HasTransformXYZ && metric::traits::HasConvertXYZ && metric::traits::HasTransform_i && metric::traits::HasConvert_i @@ -120,12 +127,15 @@ namespace kernel::sr { std::is_same>::value, "Invalid emission policy E for Pusher_kernel"); - private: const ParticlePusherFlags pusher_flags; const RadiativeDragFlags radiative_drag_flags; const randacc_ndfield_t EB; + public: + const spidx_t sp; + simtime_t time; + const real_t dt; array_t i1, i2, i3; array_t i1_prev, i2_prev, i3_prev; array_t dx1, dx2, dx3; @@ -135,23 +145,26 @@ namespace kernel::sr { array_t weight; array_t tag; const M metric; - const F external_fields; + const int ni1, ni2, ni3; - const E emission_policy; + private: + const F external_fields; + + const E emission_policy; + const PUPD custom_prtl_update; - const real_t coeff, dt; + const real_t coeff; - const int ni1, ni2, ni3; - bool is_absorb_i1min { false }, is_absorb_i1max { false }; - bool is_absorb_i2min { false }, is_absorb_i2max { false }; - bool is_absorb_i3min { false }, is_absorb_i3max { false }; - bool is_periodic_i1min { false }, is_periodic_i1max { false }; - bool is_periodic_i2min { false }, is_periodic_i2max { false }; - bool is_periodic_i3min { false }, is_periodic_i3max { false }; - bool is_reflect_i1min { false }, is_reflect_i1max { false }; - bool is_reflect_i2min { false }, is_reflect_i2max { false }; - bool is_reflect_i3min { false }, is_reflect_i3max { false }; - bool is_axis_i2min { false }, is_axis_i2max { false }; + bool is_absorb_i1min { false }, is_absorb_i1max { false }; + bool is_absorb_i2min { false }, is_absorb_i2max { false }; + bool is_absorb_i3min { false }, is_absorb_i3max { false }; + bool is_periodic_i1min { false }, is_periodic_i1max { false }; + bool is_periodic_i2min { false }, is_periodic_i2max { false }; + bool is_periodic_i3min { false }, is_periodic_i3max { false }; + bool is_reflect_i1min { false }, is_reflect_i1max { false }; + bool is_reflect_i2min { false }, is_reflect_i2max { false }; + bool is_reflect_i3min { false }, is_reflect_i3max { false }; + bool is_axis_i2min { false }, is_axis_i2max { false }; // gca parameters const real_t gca_larmor, gca_EovrB_sqr; @@ -167,10 +180,12 @@ namespace kernel::sr { const randacc_ndfield_t& EB, const M& metric, const F& external_fields, - const E& emission_policy) + const E& emission_policy, + const PUPD& custom_prtl_update) : pusher_flags { pusher_params.pusher_flags } , radiative_drag_flags { pusher_params.radiative_drag_flags } , EB { EB } + , sp { pusher_arrays.sp } , i1 { pusher_arrays.i1 } , i2 { pusher_arrays.i2 } , i3 { pusher_arrays.i3 } @@ -192,6 +207,8 @@ namespace kernel::sr { , metric { metric } , external_fields { external_fields } , emission_policy { emission_policy } + , custom_prtl_update { custom_prtl_update } + , time { pusher_params.time } , coeff { HALF * (pusher_params.charge / pusher_params.mass) * pusher_params.omegaB0 * pusher_params.dt } , dt { pusher_params.dt } @@ -563,8 +580,8 @@ namespace kernel::sr { dx1_prev(p) = dx1(p); dx1(p) += metric.template transform<1, Idx::XYZ, Idx::U>(xp, ux1(p)) * dt_inv_energy; - i1(p) += static_cast(dx1(p) >= ONE) - - static_cast(dx1(p) < ZERO); + i1(p) += static_cast(dx1(p) >= ONE) - + static_cast(dx1(p) < ZERO); dx1(p) -= (dx1(p) >= ONE); dx1(p) += (dx1(p) < ZERO); } @@ -573,8 +590,8 @@ namespace kernel::sr { dx2_prev(p) = dx2(p); dx2(p) += metric.template transform<2, Idx::XYZ, Idx::U>(xp, ux2(p)) * dt_inv_energy; - i2(p) += static_cast(dx2(p) >= ONE) - - static_cast(dx2(p) < ZERO); + i2(p) += static_cast(dx2(p) >= ONE) - + static_cast(dx2(p) < ZERO); dx2(p) -= (dx2(p) >= ONE); dx2(p) += (dx2(p) < ZERO); } @@ -583,8 +600,8 @@ namespace kernel::sr { dx3_prev(p) = dx3(p); dx3(p) += metric.template transform<3, Idx::XYZ, Idx::U>(xp, ux3(p)) * dt_inv_energy; - i3(p) += static_cast(dx3(p) >= ONE) - - static_cast(dx3(p) < ZERO); + i3(p) += static_cast(dx3(p) >= ONE) - + static_cast(dx3(p) < ZERO); dx3(p) -= (dx3(p) >= ONE); dx3(p) += (dx3(p) < ZERO); } @@ -595,8 +612,8 @@ namespace kernel::sr { : ONE / math::sqrt(SQR(ux1(p)) + SQR(ux2(p)) + SQR(ux3(p))) }; vec_t vp_Cart { ux1(p) * inv_energy, - ux2(p) * inv_energy, - ux3(p) * inv_energy }; + ux2(p) * inv_energy, + ux3(p) * inv_energy }; // get cartesian position coord_t xp_Cart { ZERO }; metric.template convert_xyz(xp, xp_Cart); @@ -631,6 +648,11 @@ namespace kernel::sr { from_Xi_to_i_di(xp[2], i3(p), dx3(p)); } } + + // call custom particle position update + custom_prtl_update(p, xp, *this); + + // apply boundary conditions boundaryConditions(p, xp); } @@ -695,11 +717,11 @@ namespace kernel::sr { auto COEFF2 { ONE + NORM_SQR(u1[0], u1[1], u1[2]) - NORM_SQR(b0[0], b0[1], b0[2]) }; - COEFF = ONE / - math::sqrt( - INV_2 * (COEFF2 + math::sqrt(SQR(COEFF2) + - FOUR * (SQR(b0[0]) + SQR(b0[1]) + - SQR(b0[2]) + SQR(COEFF))))); + COEFF = ONE / + math::sqrt( + INV_2 * (COEFF2 + math::sqrt(SQR(COEFF2) + + FOUR * (SQR(b0[0]) + SQR(b0[1]) + + SQR(b0[2]) + SQR(COEFF))))); COEFF2 = ONE / (ONE + SQR(b0[0] * COEFF) + SQR(b0[1] * COEFF) + SQR(b0[2] * COEFF)); @@ -1494,8 +1516,8 @@ namespace kernel::sr { const vec_t& e0, const vec_t& b0) const { real_t gamma_prime_sqr = ONE / math::sqrt(ONE + NORM_SQR(u_prime[0], - u_prime[1], - u_prime[2])); + u_prime[1], + u_prime[2])); u_prime[0] *= gamma_prime_sqr; u_prime[1] *= gamma_prime_sqr; u_prime[2] *= gamma_prime_sqr; @@ -1545,8 +1567,8 @@ namespace kernel::sr { Inline void inverseComptonDrag(index_t p, vec_t& u_prime) const { real_t gamma_prime_sqr = ONE / math::sqrt(ONE + NORM_SQR(u_prime[0], - u_prime[1], - u_prime[2])); + u_prime[1], + u_prime[2])); u_prime[0] *= gamma_prime_sqr; u_prime[1] *= gamma_prime_sqr; u_prime[2] *= gamma_prime_sqr; @@ -1567,7 +1589,7 @@ namespace kernel::sr { { typename E::Payload payload; vec_t delta_u_Ph { ZERO }; - const auto emission_response = emission_policy.shouldEmit(xp_Cd, + const auto emission_response = emission_policy.shouldEmit(xp_Cd, xp_Ph, u_prime, ep_Cart, @@ -1584,8 +1606,8 @@ namespace kernel::sr { if (emission_response.first) { vec_t direction { ZERO }; const auto delta_u_Ph_mag = NORM(delta_u_Ph[0], - delta_u_Ph[1], - delta_u_Ph[2]); + delta_u_Ph[1], + delta_u_Ph[2]); direction[0] = -delta_u_Ph[0] / delta_u_Ph_mag; direction[1] = -delta_u_Ph[1] / delta_u_Ph_mag; direction[2] = -delta_u_Ph[2] / delta_u_Ph_mag; diff --git a/src/kernels/tests/CMakeLists.txt b/src/kernels/tests/CMakeLists.txt index 7b73642a..f1438e10 100644 --- a/src/kernels/tests/CMakeLists.txt +++ b/src/kernels/tests/CMakeLists.txt @@ -37,3 +37,4 @@ gen_test(pusher) gen_test(ext_force) gen_test(reduced_stats) gen_test(custom_emission) +gen_test(custom_prtl_update) diff --git a/src/kernels/tests/custom_emission.cpp b/src/kernels/tests/custom_emission.cpp index 94d63112..dfc88ded 100644 --- a/src/kernels/tests/custom_emission.cpp +++ b/src/kernels/tests/custom_emission.cpp @@ -248,6 +248,7 @@ auto main(int argc, char* argv[]) -> int { }; ndfield_t EB { "EB", 128u + 2u * N_GHOSTS }; + const auto no_custom_update = kernel::sr::NoCustomPrtlUpdate_t> {}; for (auto step = 0u; step < 7u; ++step) { pusher_params.time = static_cast(step) * delta_t; @@ -286,13 +287,14 @@ auto main(int argc, char* argv[]) -> int { Kokkos::parallel_for( "ParticlePusher", 2u, - kernel::sr::Pusher_kernel( + kernel::sr::Pusher_kernel( pusher_params, pusher_arrays, EB, metric, kernel::sr::NoField_t {}, - emission_policy)); + emission_policy, + no_custom_update)); const auto n_injected = emission_policy.numbers_injected(); emitted_species_1.set_counter(emitted_species_1.counter() + n_injected[0]); emitted_species_1.set_npart(emitted_species_1.npart() + n_injected[0]); diff --git a/src/kernels/tests/custom_prtl_update.cpp b/src/kernels/tests/custom_prtl_update.cpp new file mode 100644 index 00000000..c26ad76b --- /dev/null +++ b/src/kernels/tests/custom_prtl_update.cpp @@ -0,0 +1,291 @@ +#include "enums.h" +#include "global.h" + +#include "arch/kokkos_aliases.h" +#include "utils/formatting.h" +#include "utils/numeric.h" + +#include "metrics/minkowski.h" + +#include "kernels/emission/emission.hpp" +#include "kernels/particle_pusher_sr.hpp" + +#include +#include + +#include +#include +#include +#include +#include + +using namespace ntt; +using namespace metric; + +void errorIf(bool condition, const std::string& message = "") { + if (condition) { + throw std::runtime_error(message); + } +} + +Inline auto equal(real_t a, real_t b, const std::string& msg) -> bool { + if (not(math::abs(a - b) < 1e-4)) { + Kokkos::printf("%.12e != %.12e %s\n", a, b, msg.c_str()); + return false; + } + return true; +} + +template +void put_value(array_t& arr, T v, index_t p) { + auto h = Kokkos::create_mirror_view(arr); + Kokkos::deep_copy(h, arr); + h(p) = v; + Kokkos::deep_copy(arr, h); +} + +struct TestCustomPrtlUpdate { + template + Inline void operator()(index_t p, Coord& xp, const PusherKernel& pusher) const { + if constexpr (PusherKernel::D == Dim::_1D || PusherKernel::D == Dim::_2D || PusherKernel::D == Dim::_3D) { + auto invert_vel = false; + if (pusher.i1(p) < 0) { + pusher.i1(p) = 0; + pusher.dx1(p) = ONE - pusher.dx1(p); + invert_vel = true; + } else if (pusher.i1(p) >= pusher.ni1) { + pusher.i1(p) = pusher.ni1 - 1; + pusher.dx1(p) = ONE - pusher.dx1(p); + invert_vel = true; + } + if (invert_vel) { + pusher.ux1(p) = -pusher.ux1(p); + } + } + if constexpr (PusherKernel::D == Dim::_2D || PusherKernel::D == Dim::_3D) { + auto invert_vel = false; + if (pusher.i2(p) < 0) { + pusher.i2(p) = 0; + pusher.dx2(p) = ONE - pusher.dx2(p); + invert_vel = true; + } else if (pusher.i2(p) >= pusher.ni2) { + pusher.i2(p) = pusher.ni2 - 1; + pusher.dx2(p) = ONE - pusher.dx2(p); + invert_vel = true; + } + if (invert_vel) { + pusher.ux2(p) = -pusher.ux2(p); + } + } + if constexpr (PusherKernel::D == Dim::_3D) { + auto invert_vel = false; + if (pusher.i3(p) < 0) { + pusher.i3(p) = 0; + pusher.dx3(p) = ONE - pusher.dx3(p); + invert_vel = true; + } else if (pusher.i3(p) >= pusher.ni3) { + pusher.i3(p) = pusher.ni3 - 1; + pusher.dx3(p) = ONE - pusher.dx3(p); + invert_vel = true; + } + if (invert_vel) { + pusher.ux3(p) = -pusher.ux3(p); + } + } + } +}; + +template +void testCustomPrtlUpdate(const std::vector& res, + const boundaries_t& ext, + const std::map& params = {}) { + + errorIf(res.size() != M::Dim, "res.size() != M::Dim"); + errorIf(M::CoordType != Coord::Cart, "M::CoordType != Coord::Cart"); + + M metric { res, ext, params }; + + int nx1 = 0, nx2 = 0, nx3 = 0; + if (res.size() > 0) nx1 = static_cast(res.at(0)); + if (res.size() > 1) nx2 = static_cast(res.at(1)); + if (res.size() > 2) nx3 = static_cast(res.at(2)); + + const real_t dt = 0.1 * (ext.at(0).second - ext.at(0).first) / static_cast(nx1); + + ndfield_t emfield; + if constexpr (M::Dim == Dim::_1D) { + emfield = ndfield_t { "emfield", res.at(0) + 2 * N_GHOSTS }; + } else if constexpr (M::Dim == Dim::_2D) { + emfield = ndfield_t { "emfield", res.at(0) + 2 * N_GHOSTS, res.at(1) + 2 * N_GHOSTS }; + } else { + emfield = ndfield_t { "emfield", res.at(0) + 2 * N_GHOSTS, res.at(1) + 2 * N_GHOSTS, res.at(2) + 2 * N_GHOSTS }; + } + + // Arrays for REFLECT boundaries + array_t r_i1 { "r_i1", 1 }, r_i2 { "r_i2", 1 }, r_i3 { "r_i3", 1 }; + array_t r_i1_prev { "r_i1_prev", 1 }, r_i2_prev { "r_i2_prev", 1 }, r_i3_prev { "r_i3_prev", 1 }; + array_t r_dx1 { "r_dx1", 1 }, r_dx2 { "r_dx2", 1 }, r_dx3 { "r_dx3", 1 }; + array_t r_dx1_prev { "r_dx1_prev", 1 }, r_dx2_prev { "r_dx2_prev", 1 }, r_dx3_prev { "r_dx3_prev", 1 }; + array_t r_ux1 { "r_ux1", 1 }, r_ux2 { "r_ux2", 1 }, r_ux3 { "r_ux3", 1 }; + array_t r_tag { "r_tag", 1 }; + + // Arrays for CUSTOM boundaries + array_t c_i1 { "c_i1", 1 }, c_i2 { "c_i2", 1 }, c_i3 { "c_i3", 1 }; + array_t c_i1_prev { "c_i1_prev", 1 }, c_i2_prev { "c_i2_prev", 1 }, c_i3_prev { "c_i3_prev", 1 }; + array_t c_dx1 { "c_dx1", 1 }, c_dx2 { "c_dx2", 1 }, c_dx3 { "c_dx3", 1 }; + array_t c_dx1_prev { "c_dx1_prev", 1 }, c_dx2_prev { "c_dx2_prev", 1 }, c_dx3_prev { "c_dx3_prev", 1 }; + array_t c_ux1 { "c_ux1", 1 }, c_ux2 { "c_ux2", 1 }, c_ux3 { "c_ux3", 1 }; + array_t c_tag { "c_tag", 1 }; + + array_t phi; + + real_t ux_1 = 0.569197; + real_t uy_1 = 0.716085; + real_t uz_1 = -0.760101; + put_value(r_ux1, ux_1, 0); put_value(c_ux1, ux_1, 0); + put_value(r_ux2, uy_1, 0); put_value(c_ux2, uy_1, 0); + put_value(r_ux3, uz_1, 0); put_value(c_ux3, uz_1, 0); + + put_value(r_tag, ParticleTag::alive, 0); + put_value(c_tag, ParticleTag::alive, 0); + + if constexpr (M::PrtlDim >= Dim::_1D) { + put_value(r_i1, nx1 - 1, 0); put_value(c_i1, nx1 - 1, 0); + put_value(r_dx1, 0.9, 0); put_value(c_dx1, 0.9, 0); + } + if constexpr (M::PrtlDim >= Dim::_2D) { + put_value(r_i2, nx2 - 1, 0); put_value(c_i2, nx2 - 1, 0); + put_value(r_dx2, 0.9, 0); put_value(c_dx2, 0.9, 0); + } + if constexpr (M::PrtlDim == Dim::_3D) { + put_value(r_i3, 0, 0); put_value(c_i3, 0, 0); + put_value(r_dx3, 0.1, 0); put_value(c_dx3, 0.1, 0); + } + + kernel::sr::PusherParams r_params {}; + r_params.pusher_flags = ParticlePusher::BORIS; + r_params.mass = ONE; + r_params.charge = ONE; + r_params.dt = dt; + r_params.omegaB0 = ONE; + r_params.ni1 = nx1; + r_params.ni2 = nx2; + r_params.ni3 = nx3; + r_params.boundaries = { + { PrtlBC::REFLECT, PrtlBC::REFLECT }, + { PrtlBC::REFLECT, PrtlBC::REFLECT }, + { PrtlBC::REFLECT, PrtlBC::REFLECT } + }; + + kernel::sr::PusherParams c_params = r_params; + // initialize with periodic boundaries so that reflection is only handled by the custom update + c_params.boundaries = { + { PrtlBC::PERIODIC, PrtlBC::PERIODIC }, + { PrtlBC::PERIODIC, PrtlBC::PERIODIC }, + { PrtlBC::PERIODIC, PrtlBC::PERIODIC } + }; + + kernel::sr::PusherArrays r_arrays {}; + r_arrays.sp = 1u; + r_arrays.i1 = r_i1; r_arrays.i2 = r_i2; r_arrays.i3 = r_i3; + r_arrays.i1_prev = r_i1_prev; r_arrays.i2_prev = r_i2_prev; r_arrays.i3_prev = r_i3_prev; + r_arrays.dx1 = r_dx1; r_arrays.dx2 = r_dx2; r_arrays.dx3 = r_dx3; + r_arrays.dx1_prev = r_dx1_prev; r_arrays.dx2_prev = r_dx2_prev; r_arrays.dx3_prev = r_dx3_prev; + r_arrays.ux1 = r_ux1; r_arrays.ux2 = r_ux2; r_arrays.ux3 = r_ux3; + r_arrays.phi = phi; + r_arrays.tag = r_tag; + + kernel::sr::PusherArrays c_arrays {}; + c_arrays.sp = 1u; + c_arrays.i1 = c_i1; c_arrays.i2 = c_i2; c_arrays.i3 = c_i3; + c_arrays.i1_prev = c_i1_prev; c_arrays.i2_prev = c_i2_prev; c_arrays.i3_prev = c_i3_prev; + c_arrays.dx1 = c_dx1; c_arrays.dx2 = c_dx2; c_arrays.dx3 = c_dx3; + c_arrays.dx1_prev = c_dx1_prev; c_arrays.dx2_prev = c_dx2_prev; c_arrays.dx3_prev = c_dx3_prev; + c_arrays.ux1 = c_ux1; c_arrays.ux2 = c_ux2; c_arrays.ux3 = c_ux3; + c_arrays.phi = phi; + c_arrays.tag = c_tag; + + const auto no_emission = kernel::NoEmissionPolicy_t {}; + const auto no_custom_update = kernel::sr::NoCustomPrtlUpdate_t {}; + const auto custom_update = TestCustomPrtlUpdate {}; + + const auto n_iter = 100; + + for (auto n { 0 }; n < n_iter; ++n) { + Kokkos::parallel_for( + "pusher_reflect", + CreateRangePolicy({ 0 }, { 1 }), + kernel::sr::Pusher_kernel( + r_params, r_arrays, emfield, metric, kernel::sr::NoField_t {}, no_emission, no_custom_update) + ); + Kokkos::parallel_for( + "pusher_custom", + CreateRangePolicy({ 0 }, { 1 }), + kernel::sr::Pusher_kernel( + c_params, c_arrays, emfield, metric, kernel::sr::NoField_t {}, no_emission, custom_update) + ); + + auto hr_i1 = Kokkos::create_mirror_view(r_i1); Kokkos::deep_copy(hr_i1, r_i1); + auto hc_i1 = Kokkos::create_mirror_view(c_i1); Kokkos::deep_copy(hc_i1, c_i1); + auto hr_dx1 = Kokkos::create_mirror_view(r_dx1); Kokkos::deep_copy(hr_dx1, r_dx1); + auto hc_dx1 = Kokkos::create_mirror_view(c_dx1); Kokkos::deep_copy(hc_dx1, c_dx1); + auto hr_ux1 = Kokkos::create_mirror_view(r_ux1); Kokkos::deep_copy(hr_ux1, r_ux1); + auto hc_ux1 = Kokkos::create_mirror_view(c_ux1); Kokkos::deep_copy(hc_ux1, c_ux1); + + if constexpr (M::PrtlDim >= Dim::_1D) { + errorIf(hr_i1(0) != hc_i1(0), "i1 mismatch"); + errorIf(not equal(hr_dx1(0), hc_dx1(0), "dx1 mismatch")); + errorIf(not equal(hr_ux1(0), hc_ux1(0), "ux1 mismatch")); + } + + if constexpr (M::PrtlDim >= Dim::_2D) { + auto hr_i2 = Kokkos::create_mirror_view(r_i2); Kokkos::deep_copy(hr_i2, r_i2); + auto hc_i2 = Kokkos::create_mirror_view(c_i2); Kokkos::deep_copy(hc_i2, c_i2); + auto hr_dx2 = Kokkos::create_mirror_view(r_dx2); Kokkos::deep_copy(hr_dx2, r_dx2); + auto hc_dx2 = Kokkos::create_mirror_view(c_dx2); Kokkos::deep_copy(hc_dx2, c_dx2); + auto hr_ux2 = Kokkos::create_mirror_view(r_ux2); Kokkos::deep_copy(hr_ux2, r_ux2); + auto hc_ux2 = Kokkos::create_mirror_view(c_ux2); Kokkos::deep_copy(hc_ux2, c_ux2); + errorIf(hr_i2(0) != hc_i2(0), "i2 mismatch"); + errorIf(not equal(hr_dx2(0), hc_dx2(0), "dx2 mismatch")); + errorIf(not equal(hr_ux2(0), hc_ux2(0), "ux2 mismatch")); + } + + if constexpr (M::PrtlDim == Dim::_3D) { + auto hr_i3 = Kokkos::create_mirror_view(r_i3); Kokkos::deep_copy(hr_i3, r_i3); + auto hc_i3 = Kokkos::create_mirror_view(c_i3); Kokkos::deep_copy(hc_i3, c_i3); + auto hr_dx3 = Kokkos::create_mirror_view(r_dx3); Kokkos::deep_copy(hr_dx3, r_dx3); + auto hc_dx3 = Kokkos::create_mirror_view(c_dx3); Kokkos::deep_copy(hc_dx3, c_dx3); + auto hr_ux3 = Kokkos::create_mirror_view(r_ux3); Kokkos::deep_copy(hr_ux3, r_ux3); + auto hc_ux3 = Kokkos::create_mirror_view(c_ux3); Kokkos::deep_copy(hc_ux3, c_ux3); + errorIf(hr_i3(0) != hc_i3(0), "i3 mismatch"); + errorIf(not equal(hr_dx3(0), hc_dx3(0), "dx3 mismatch")); + errorIf(not equal(hr_ux3(0), hc_ux3(0), "ux3 mismatch")); + } + } +} + +auto main(int argc, char* argv[]) -> int { + Kokkos::initialize(argc, argv); + + try { + using namespace ntt; + + const std::vector res1d { 50 }; + const boundaries_t ext1d { { 0.0, 1000.0 } }; + const std::vector res2d { 30, 20 }; + const boundaries_t ext2d { { -15.0, 15.0 }, { -10.0, 10.0 } }; + const std::vector res3d { 10, 10, 10 }; + const boundaries_t ext3d { { 0.0, 1.0 }, { 0.0, 1.0 }, { 0.0, 1.0 } }; + + testCustomPrtlUpdate>(res1d, ext1d, {}); + testCustomPrtlUpdate>(res2d, ext2d, {}); + testCustomPrtlUpdate>(res3d, ext3d, {}); + + } catch (std::exception& e) { + std::cerr << e.what() << std::endl; + Kokkos::finalize(); + return 1; + } + Kokkos::finalize(); + return 0; +} diff --git a/src/kernels/tests/ext_force.cpp b/src/kernels/tests/ext_force.cpp index 4e1b4ac3..29b22d52 100644 --- a/src/kernels/tests/ext_force.cpp +++ b/src/kernels/tests/ext_force.cpp @@ -198,6 +198,7 @@ void testPusher(const std::vector& res) { const auto no_emission = kernel::NoEmissionPolicy_t> {}; + const auto no_custom_update = kernel::sr::NoCustomPrtlUpdate_t> {}; for (auto t { 0u }; t < 100; ++t) { const real_t time = t * dt; @@ -206,13 +207,14 @@ void testPusher(const std::vector& res) { Kokkos::parallel_for( "pusher", CreateRangePolicy({ 0 }, { 2 }), - kernel::sr::Pusher_kernel, decltype(ext_force), false, decltype(no_emission)>( + kernel::sr::Pusher_kernel, decltype(ext_force), false, decltype(no_emission), decltype(no_custom_update)>( pusher_params, pusher_arrays, emfield, metric, ext_force, - no_emission)); + no_emission, + no_custom_update)); auto i1_prev_ = Kokkos::create_mirror_view(i1_prev); auto i2_prev_ = Kokkos::create_mirror_view(i2_prev); diff --git a/src/kernels/tests/gca_pusher.cpp b/src/kernels/tests/gca_pusher.cpp index a5bd3479..1c198e02 100644 --- a/src/kernels/tests/gca_pusher.cpp +++ b/src/kernels/tests/gca_pusher.cpp @@ -181,6 +181,7 @@ void testPusher(const std::vector& res) { pusher_arrays.tag = tag; const auto no_emission = kernel::NoEmissionPolicy_t> {}; + const auto no_custom_update = kernel::sr::NoCustomPrtlUpdate_t> {}; for (auto t { 0u }; t < 2000; ++t) { pusher_params.time = t * dt; @@ -188,13 +189,14 @@ void testPusher(const std::vector& res) { Kokkos::parallel_for( "pusher", CreateRangePolicy({ 0 }, { 2 }), - kernel::sr::Pusher_kernel, kernel::sr::NoField_t, false, decltype(no_emission)>( + kernel::sr::Pusher_kernel, kernel::sr::NoField_t, false, decltype(no_emission), decltype(no_custom_update)>( pusher_params, pusher_arrays, emfield, metric, kernel::sr::NoField_t {}, - no_emission)); + no_emission, + no_custom_update)); auto ux1_ = Kokkos::create_mirror_view(ux1); auto ux2_ = Kokkos::create_mirror_view(ux2); diff --git a/src/kernels/tests/prtl_bc.cpp b/src/kernels/tests/prtl_bc.cpp index bbdfa57b..c459baa5 100644 --- a/src/kernels/tests/prtl_bc.cpp +++ b/src/kernels/tests/prtl_bc.cpp @@ -266,18 +266,20 @@ void testPeriodicBC(const std::vector& res, pusher_arrays.phi = phi; pusher_arrays.tag = tag; const auto no_emission = kernel::NoEmissionPolicy_t {}; + const auto no_custom_update = kernel::sr::NoCustomPrtlUpdate_t {}; for (auto n { 0 }; n < n_iter; ++n) { Kokkos::parallel_for( "pusher", CreateRangePolicy({ 0 }, { 2 }), - kernel::sr::Pusher_kernel( + kernel::sr::Pusher_kernel( pusher_params, pusher_arrays, emfield, metric, kernel::sr::NoField_t {}, - no_emission)); + no_emission, + no_custom_update)); auto i1_ = Kokkos::create_mirror_view(i1); auto i2_ = Kokkos::create_mirror_view(i2); auto i3_ = Kokkos::create_mirror_view(i3); diff --git a/src/kernels/tests/pusher.cpp b/src/kernels/tests/pusher.cpp index 3bd951de..ef5f69b8 100644 --- a/src/kernels/tests/pusher.cpp +++ b/src/kernels/tests/pusher.cpp @@ -179,6 +179,7 @@ void testPusher(const std::vector& res) { pusher_arrays.tag = tag; const auto no_emission = kernel::NoEmissionPolicy_t> {}; + const auto no_custom_update = kernel::sr::NoCustomPrtlUpdate_t> {}; for (auto t { 0u }; t < 2000; ++t) { const real_t time = t * dt; @@ -188,26 +189,28 @@ void testPusher(const std::vector& res) { Kokkos::parallel_for( "pusher", CreateRangePolicy({ 0 }, { 1 }), - kernel::sr::Pusher_kernel, kernel::sr::NoField_t, false, decltype(no_emission)>( + kernel::sr::Pusher_kernel, kernel::sr::NoField_t, false, decltype(no_emission), decltype(no_custom_update)>( pusher_params, pusher_arrays, emfield, metric, kernel::sr::NoField_t {}, - no_emission)); + no_emission, + no_custom_update)); pusher_params.pusher_flags = ParticlePusher::VAY; Kokkos::parallel_for( "pusher", CreateRangePolicy({ 1 }, { 2 }), - kernel::sr::Pusher_kernel, kernel::sr::NoField_t, false, decltype(no_emission)>( + kernel::sr::Pusher_kernel, kernel::sr::NoField_t, false, decltype(no_emission), decltype(no_custom_update)>( pusher_params, pusher_arrays, emfield, metric, kernel::sr::NoField_t {}, - no_emission)); + no_emission, + no_custom_update)); auto i1_prev_ = Kokkos::create_mirror_view(i1_prev); auto i2_prev_ = Kokkos::create_mirror_view(i2_prev);