From b46a0562f4d6a5da9e4b6e7096f6f8994d1ab68b Mon Sep 17 00:00:00 2001 From: Andrew Sullivan Date: Tue, 17 Feb 2026 10:24:06 -0800 Subject: [PATCH] spect3D included --- pgens/stripedwind_new/pgen.hpp | 580 +++++++++++++ pgens/stripedwind_new/stripedwind.toml | 73 ++ src/archetypes/particle_injector.h | 4 +- src/archetypes/utils.h | 4 +- src/engines/engine_init.cpp | 21 +- src/engines/engine_printer.cpp | 27 +- src/engines/engine_run.cpp | 21 +- src/engines/engine_traits.h | 30 - src/engines/grpic.hpp | 4 +- src/engines/srpic.hpp | 64 +- src/entity.cpp | 134 ++- src/framework/containers/fields_io.cpp | 6 +- src/framework/containers/particles_io.cpp | 24 +- src/framework/domain/checkpoint.cpp | 39 +- src/framework/domain/communications.cpp | 34 +- src/framework/domain/domain.h | 19 +- src/framework/domain/grid.cpp | 1 + src/framework/domain/metadomain.cpp | 26 +- src/framework/domain/output.cpp | 290 ++++++- src/framework/domain/stats.cpp | 38 +- src/framework/domain/test.cpp | 43 + src/framework/parameters.cpp | 59 +- src/framework/specialization_registry.h | 60 -- src/global/arch/traits.h | 3 - src/global/defaults.h | 13 +- src/global/enums.h | 4 +- src/global/global.h | 9 +- src/global/utils/numeric.h | 64 +- src/kernels/currents_deposit.hpp | 806 +++++------------ src/kernels/digital_filter.hpp | 49 +- src/kernels/faraday_mink.hpp | 1 + src/kernels/injectors.hpp | 23 +- src/kernels/particle_moments.hpp | 2 +- src/kernels/particle_pusher_gr.hpp | 4 +- src/kernels/particle_pusher_sr.hpp | 783 ++++++----------- src/kernels/particle_shapes.hpp | 999 ---------------------- src/kernels/tests/deposit.cpp | 185 ++-- src/kernels/tests/faraday_mink.cpp | 1 - src/output/utils/interpret_prompt.h | 4 +- src/output/utils/readers.cpp | 21 +- src/output/utils/writers.cpp | 18 +- src/output/utils/writers.h | 9 + src/output/writer.cpp | 154 +++- src/output/writer.h | 2 + 44 files changed, 2077 insertions(+), 2678 deletions(-) create mode 100644 pgens/stripedwind_new/pgen.hpp create mode 100644 pgens/stripedwind_new/stripedwind.toml delete mode 100644 src/engines/engine_traits.h create mode 100644 src/framework/domain/test.cpp delete mode 100644 src/framework/specialization_registry.h delete mode 100644 src/kernels/particle_shapes.hpp diff --git a/pgens/stripedwind_new/pgen.hpp b/pgens/stripedwind_new/pgen.hpp new file mode 100644 index 000000000..d5e34958a --- /dev/null +++ b/pgens/stripedwind_new/pgen.hpp @@ -0,0 +1,580 @@ +#ifndef PROBLEM_GENERATOR_H +#define PROBLEM_GENERATOR_H + +#include "enums.h" +#include "global.h" + +#include "arch/traits.h" +#include "utils/error.h" +#include "utils/numeric.h" + +#include "archetypes/field_setter.h" +#include "archetypes/problem_generator.h" +#include "archetypes/utils.h" +#include "archetypes/spatial_dist.h" +#include "framework/domain/metadomain.h" + +#include +#include + +namespace user { + using namespace ntt; + + template + struct SWCurrentSheet : public arch::SpatialDistribution { + SWCurrentSheet(const M& metric, + real_t alpha_s, + real_t lambda_s, + real_t delta_s, + real_t drift_beta_x, + real_t time) + : arch::SpatialDistribution { metric } + , alpha_s { alpha_s } + , lambda_s { lambda_s } + , delta_s { delta_s } + , drift_beta_x { drift_beta_x } + , time { time } {} + + Inline auto operator()(const coord_t& x_Ph) const -> real_t { + return ONE / SQR(math::cosh( + lambda_s / (constant::TWO_PI * delta_s) * + (alpha_s + + math::cos(constant::TWO_PI * (x_Ph[0] + drift_beta_x * time) / lambda_s)))); + } + + private: + real_t alpha_s, lambda_s, delta_s, drift_beta_x, time; + }; + + template + struct InitBfields { + InitBfields(real_t bmag, + real_t btheta, + real_t bphi, + real_t drift_beta_x, + real_t lambda_s, + real_t delta_s, + real_t alpha_s, + real_t time = ZERO) + : Bmag { bmag } + , Btheta { btheta * static_cast(convert::deg2rad) } + , Bphi { bphi * static_cast(convert::deg2rad) } + , drift_beta_x { drift_beta_x } + , Lambda_s { lambda_s } + , Delta_s { delta_s } + , Alpha_s { alpha_s } + , time { time } {} + + // magnetic field components + + Inline auto bx1(const coord_t&) const -> real_t { + return ZERO; + } + + Inline auto bx2(const coord_t& x_Ph) const -> real_t { + + return Bmag * + math::tanh(Lambda_s / (constant::TWO_PI * Delta_s) * + (Alpha_s + math::cos(constant::TWO_PI * + (x_Ph[0] + time * drift_beta_x) / Lambda_s))); + } + + Inline auto bx3(const coord_t&) const -> real_t { + return ZERO; + } + + const real_t Btheta, Bphi, drift_beta_x, Bmag, Delta_s, Alpha_s, Lambda_s; + const real_t time; + }; + + template + struct InitFields : public InitBfields { + /* + 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_beta_x, + real_t lambda_s, + real_t delta_s, + real_t alpha_s, + real_t time = ZERO) + : InitBfields { bmag, btheta, bphi, drift_beta_x, + lambda_s, delta_s, alpha_s, time } {} + + // magnetic field components + + // electric field components + Inline auto ex1(const coord_t&) const -> real_t { + return ZERO; + } + + Inline auto ex2(const coord_t&) const -> real_t { + return ZERO; + } + + Inline auto ex3(const coord_t& x_Ph) const -> real_t { + return -this->drift_beta_x * this->Bmag * + math::tanh( + this->Lambda_s / (constant::TWO_PI * this->Delta_s) * + (this->Alpha_s + math::cos(constant::TWO_PI * (x_Ph[0]) / this->Lambda_s))); + } + }; + + template + struct PGen : public arch::ProblemGenerator { + // compatibility traits for the problem generator + static constexpr auto engines { traits::compatible_with::value }; + static constexpr auto metrics { traits::compatible_with::value }; + static constexpr auto dimensions { + traits::compatible_with::value + }; + + // 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 global_xmin, global_xmax; + // gas properties + const real_t drift_ux, drift_gamma, drift_beta_x, temperature, temperature_ratio, filling_fraction; + // injector properties + const real_t injector_velocity, injection_start, dt; + const int injector_interval, delta_x_resetb; + // magnetic field properties + real_t Btheta, Bphi, Bmag, Alpha_s, Lambda_s, Delta_s, eta_s; + InitFields init_flds; + + real_t injector_x; + + 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) } + , Delta_s { p.template get("setup.delta_s", 5.0) } + , Lambda_s { p.template get("setup.lambda_s", 100.0) } + , Alpha_s { p.template get("setup.alpha_s", ZERO) } + , eta_s {p.template get("setup.eta_s", 3.0) } // this is the ratio of the density in current sheet to cold wind + , drift_gamma {static_cast(math::sqrt(1.+SQR(drift_ux)))} + , drift_beta_x {drift_ux/drift_gamma} + , init_flds { Bmag, Btheta, Bphi, drift_beta_x, Lambda_s, Delta_s, Alpha_s } + , filling_fraction { p.template get("setup.filling_fraction", 1.0) } + , injector_velocity { p.template get("setup.injector_velocity", 1.0) } + , injection_start { p.template get("setup.injection_start", 0.0) } + , injector_interval { p.template get("setup.injector_interval", 100) } + , delta_x_resetb {p.template get("setup.delta_x_resetb", 50) } + , dt { p.template get("algorithms.timestep.dt") } {} + + + inline PGen() {} + + //real_t drift_gamma = math::sqrt(1.+SQR(drift_ux)); + + //real_t drift_beta_x = drift_ux / drift_gamma; + + auto MatchFields(real_t time) const -> InitBfields { + const auto init_bflds = + InitBfields(Bmag, Btheta, Bphi, drift_beta_x, Lambda_s, Delta_s, Alpha_s, time); + return init_bflds; + } + + // auto FixFieldsConst(const bc_in&, + // const em& comp) const -> std::pair { + // if (comp == em::ex1) { + // return { init_flds.ex1({ ZERO }), true }; + // } else if (comp == em::ex2) { + // return { ZERO, true }; + // } else if (comp == em::ex3) { + // return { ZERO, true }; + // } else if (comp == em::bx1) { + // return { init_flds.bx1({ ZERO }), true }; + // } else if (comp == em::bx2) { + // return { init_flds.bx2({ ZERO }), true }; + // } else if (comp == em::bx3) { + // return { init_flds.bx3({ ZERO }), true }; + // } else { + // raise::Error("Invalid component", HERE); + // return { ZERO, false }; + // } + // } + + inline void InitPrtls(Domain& domain) { + /* + * Plasma setup as partially filled box + * + * Plasma setup: + * + * global_xmin global_xmax + * | | + * V V + * |:::::::::::|..........................| + * ^ + * | + * filling_fraction + */ + + // minimum and maximum position of particles + + const auto sigma = params.template get("scales.sigma0"); + const auto c_omp = params.template get("scales.skindepth0"); + + const auto cs_temperature = HALF * sigma / eta_s; + + real_t xg_min = global_xmin; + real_t xg_max = global_xmin + filling_fraction * (global_xmax - global_xmin); + + // 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 { -drift_ux, ZERO, ZERO }, + std::vector { -drift_ux, ZERO, ZERO }); + + // inject particles + arch::InjectUniformMaxwellians(params, + domain, + ONE, + temperatures, + { 1, 2 }, + drifts, + false, + box); + + + + + const auto edist_cs = arch::Maxwellian(domain.mesh.metric, domain.random_pool, cs_temperature, { -drift_ux, ZERO, ZERO } ); + + const auto sdist_cs = SWCurrentSheet(domain.mesh.metric, Alpha_s, Lambda_s, Delta_s, drift_beta_x, ZERO); + + // inject particles in current sheet, need to include classes of energy distribution and spatial distribution + arch::InjectNonUniform( + params, + domain, + {1, 2}, + {edist_cs, edist_cs}, + sdist_cs, + eta_s, + false, + box); + injector_x = xg_max; + } + + void CustomPostStep(timestep_t step, simtime_t time, Domain& domain) { + + if (step % injector_interval == 0) { + const auto dt = params.template get("algorithms.timestep.dt"); + const auto new_injector_x = injector_x + + injector_velocity * injector_interval * dt; + const auto xmin = injector_x - injector_interval * dt; + /** + * tag particles inside the injection zone as dead + **/ + 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; + auto dx1 = species.dx1; + auto tag = species.tag; + + Kokkos::parallel_for( + "RemoveParticles", + species.rangeActiveParticles(), + Lambda(index_t p) { + // check if the particle is already dead + if (tag(p) == ParticleTag::dead) { + return; + } + // convert particle position to grid coordinates + const auto x_Cd = static_cast(i1(p)) + + static_cast(dx1(p)); + // convert grid coordinates to physical coordinates + const auto x_Ph = mesh.metric.template convert<1, Crd::Cd, Crd::XYZ>( + x_Cd); + + // if the particle position is to the right of xmin, tag it as dead + if (x_Ph > xmin) { + tag(p) = ParticleTag::dead; + } + }); + } + + /* + Reset the fields inside the purged region + */ + // define indices range to reset fields + // (not including ghost zones in either direction) + boundaries_t incl_ghosts; + for (auto d = 0; d < M::Dim; ++d) { + incl_ghosts.push_back({ false, false }); + } + + // define the rectangular box region where fields are reset + boundaries_t purge_box; + // loop over all dimension + for (auto d = 0u; d < M::Dim; ++d) { + if (d == 0) { + purge_box.push_back({ xmin, new_injector_x }); + } else { + purge_box.push_back(Range::All); + } + } + + // convert physical extent to a range of cells + const auto extent = domain.mesh.ExtentToRange(purge_box, incl_ghosts); + + // record the range min/max boundaries in each dimension + tuple_t x_min { 0 }, x_max { 0 }; + + // define a tuple to set the volume over which the magnetic field is reset to the right of the injector + tuple_t x_min_bzone { 0 }, x_max_bzone { 0 }; + + for (auto d = 0; d < M::Dim; ++d) { + x_min[d] = extent[d].first ; + x_max[d] = extent[d].second ; + if (d == 0){ + x_min_bzone[d] = x_max[d]+1 ; + x_max_bzone[d] = x_max[d]+delta_x_resetb ; + } else { + x_min_bzone[d] = x_min[d]; + x_max_bzone[d] = x_max[d]; + } + + } + // I am re-setting the fields in the region where the plasma is being injected + const auto init_bflds = InitBfields(Bmag, Btheta, Bphi, drift_beta_x, Lambda_s, Delta_s, Alpha_s, time); + + Kokkos::parallel_for("ResetFields", + CreateRangePolicy(x_min_bzone, x_max_bzone), + arch::SetEMFields_kernel { + domain.fields.em, // ^ + init_bflds, // <-- but proper injector fields + domain.mesh.metric }); + + // ADD FRESH PLASMA INJECTOR (COPY FROM BELOW) + // This is code taken from below + // inject particles now + boundaries_t inj_box; + //loop over all dimensions + for (auto d = 0u; d {-drift_ux, ZERO, ZERO}, + std::vector {-drift_ux, ZERO, ZERO}); + + arch::InjectUniformMaxwellians(params, + domain, + ONE, + temperatures, + {1, 2}, + drifts, + false, + inj_box); + + // Now inject the current sheet + const auto sigma = params.template get("scales.sigma0"); + const auto c_omp = params.template get("scales.skindepth0"); + + const auto cs_temperature = HALF * sigma / eta_s; + + const auto edist_cs = arch::Maxwellian(domain.mesh.metric, domain.random_pool, cs_temperature, { -drift_ux, ZERO, ZERO } ); + + const auto sdist_cs = SWCurrentSheet(domain.mesh.metric, Alpha_s, Lambda_s, Delta_s, drift_beta_x, ZERO); + + // inject particles in current sheet, need to include classes of energy distribution and spatial distribution + arch::InjectNonUniform( + params, + domain, + {1, 2}, + {edist_cs, edist_cs}, + sdist_cs, + eta_s, + false, + inj_box); + + injector_x = new_injector_x; + // + } + } + + // + // /* + // * Replenish plasma in a moving injector + // * + // * Injector setup: + // * + // * global_xmin purge/replenish global_xmax + // * | x_init | | + // * V v V V + // * |:::::::::::;::::::::::|\\\\\\\\|......| + // * xmin xmax + // * ^ + // * | + // * moving injector + // */ + // + // // check if the injector should be active + // if (step % injection_frequency != 0) { + // return; + // } + // + // // initial position of injector + // const auto x_init = global_xmin + + // filling_fraction * (global_xmax - global_xmin); + // + // // compute the position of the injector after the current timestep + // auto xmax = x_init + injector_velocity * + // (std::max(time - injection_start, ZERO) + dt); + // if (xmax >= global_xmax) { + // xmax = global_xmax; + // } + // + // // compute the beginning of the injected region + // auto xmin = xmax - injection_frequency * dt; + // if (xmin <= global_xmin) { + // xmin = global_xmin; + // } + // + // // define indice range to reset fields + // boundaries_t incl_ghosts; + // for (auto d = 0; d < M::Dim; ++d) { + // incl_ghosts.push_back({ false, false }); + // } + // + // // define box to reset fields + // boundaries_t purge_box; + // // loop over all dimension + // for (auto d = 0u; d < M::Dim; ++d) { + // if (d == 0) { + // purge_box.push_back({ xmin, global_xmax }); + // } else { + // purge_box.push_back(Range::All); + // } + // } + // + // const auto extent = domain.mesh.ExtentToRange(purge_box, incl_ghosts); + // tuple_t x_min { 0 }, x_max { 0 }; + // for (auto d = 0; d < M::Dim; ++d) { + // x_min[d] = extent[d].first; + // x_max[d] = extent[d].second; + // } + // + // Kokkos::parallel_for("ResetFields", + // CreateRangePolicy(x_min, x_max), + // arch::SetEMFields_kernel { + // domain.fields.em, + // init_flds, + // domain.mesh.metric }); + // global_domain.CommunicateFields(domain, Comm::E | Comm::B); + // + // /* + // tag particles inside the injection zone as dead + // */ + // 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; + // auto dx1 = species.dx1; + // auto tag = species.tag; + // + // Kokkos::parallel_for( + // "RemoveParticles", + // species.rangeActiveParticles(), + // Lambda(index_t p) { + // // check if the particle is already dead + // if (tag(p) == ParticleTag::dead) { + // return; + // } + // const auto x_Cd = static_cast(i1(p)) + + // static_cast(dx1(p)); + // const auto x_Ph = mesh.metric.template convert<1, Crd::Cd, Crd::XYZ>( + // x_Cd); + // + // if (x_Ph > xmin) { + // tag(p) = ParticleTag::dead; + // } + // }); + // } + // + // /* + // Inject slab of fresh plasma + // */ + // + // // 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 { -drift_ux, ZERO, ZERO }, + // std::vector { -drift_ux, ZERO, ZERO }); + // arch::InjectUniformMaxwellians(params, + // domain, + // ONE, + // temperatures, + // { 1, 2 }, + // drifts, + // false, + // inj_box); + }; + +} // namespace user +#endif \ No newline at end of file diff --git a/pgens/stripedwind_new/stripedwind.toml b/pgens/stripedwind_new/stripedwind.toml new file mode 100644 index 000000000..88721b342 --- /dev/null +++ b/pgens/stripedwind_new/stripedwind.toml @@ -0,0 +1,73 @@ +[simulation] + name = "stripedwind" + engine = "srpic" + runtime = 50.0 + + [simulation.domain] + decomposition = [1, -1] + +[grid] + resolution = [4096, 128] + extent = [[0.0, 409.6], [0.0, 12.8]] + + [grid.metric] + metric = "minkowski" + + [grid.boundaries] + fields = [["CONDUCTOR", "MATCH"], ["PERIODIC"]] + particles = [["REFLECT", "ABSORB"], ["PERIODIC"]] + + +[scales] + larmor0 = 0.316227766017 + skindepth0 = 1.0 + +[algorithms] + current_filters = 8 + + [algorithms.timestep] + CFL = 0.5 + +[particles] + ppc0 = 8.0 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 8e7 + + [[particles.species]] + label = "e+" + mass = 1.0 + charge = 1.0 + maxnpart = 8e7 + +[setup] + drift_ux = 0.0 # speed towards the wall [c] + temperature = 0.0001 # 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 = 0.0 # magnetic field angle in the plane + Bphi = 0.0 # magnetic field angle out of plane + alpha_s = 0.1 # stripe averaged field (see Sironi & Spitkovsky 2011) + lambda_s = 20.0 # striped wind wavelength + delta_s = 1.0 # width of region where magnetic field flips + filling_fraction = 0.1 # fraction of the shock piston filled with plasma + injector_velocity = 0.99 # speed of injector [c] + injection_start = 0.0 # start time of moving injector + injector_interval = 1 # inject particles every 1 timestep + +[output] + interval_time = 0.1 + format = "BPFile" + + [output.fields] + quantities = ["N_1", "N_2", "B", "E"] + + [output.particles] + enable = true + stride = 10 + + [output.spectra] + enable = false \ No newline at end of file diff --git a/src/archetypes/particle_injector.h b/src/archetypes/particle_injector.h index 24af59965..79743ce65 100644 --- a/src/archetypes/particle_injector.h +++ b/src/archetypes/particle_injector.h @@ -254,7 +254,7 @@ namespace arch { energy_dists.first, energy_dists.second, ONE / params.template get("scales.V0"), - domain.random_pool())); + domain.random_pool)); domain.species[species.first - 1].set_npart( domain.species[species.first - 1].npart() + nparticles); domain.species[species.second - 1].set_npart( @@ -372,7 +372,7 @@ namespace arch { energy_dists.second, spatial_dist, ONE / params.template get("scales.V0"), - domain.random_pool()); + domain.random_pool); Kokkos::parallel_for("InjectNonUniformNumberDensity", cell_range, injector_kernel); diff --git a/src/archetypes/utils.h b/src/archetypes/utils.h index b1a29cb71..3447558fb 100644 --- a/src/archetypes/utils.h +++ b/src/archetypes/utils.h @@ -53,11 +53,11 @@ namespace arch { const auto temperature_2 = temperatures.second / mass_2; const auto maxwellian_1 = arch::Maxwellian(domain.mesh.metric, - domain.random_pool(), + domain.random_pool, temperature_1, drift_four_vels.first); const auto maxwellian_2 = arch::Maxwellian(domain.mesh.metric, - domain.random_pool(), + domain.random_pool, temperature_2, drift_four_vels.second); diff --git a/src/engines/engine_init.cpp b/src/engines/engine_init.cpp index f98b117dd..7a963b615 100644 --- a/src/engines/engine_init.cpp +++ b/src/engines/engine_init.cpp @@ -3,8 +3,14 @@ #include "arch/traits.h" +#include "metrics/kerr_schild.h" +#include "metrics/kerr_schild_0.h" +#include "metrics/minkowski.h" +#include "metrics/qkerr_schild.h" +#include "metrics/qspherical.h" +#include "metrics/spherical.h" + #include "archetypes/field_setter.h" -#include "framework/specialization_registry.h" #include "engines/engine.hpp" @@ -65,10 +71,13 @@ namespace ntt { print_report(); } -#define ENGINE_INIT(S, M, D) template class Engine>; - - NTT_FOREACH_SPECIALIZATION(ENGINE_INIT) - -#undef ENGINE_INIT + template class Engine>; + template class Engine>; + template class Engine>; + template class Engine>; + template class Engine>; + template class Engine>; + template class Engine>; + template class Engine>; } // namespace ntt diff --git a/src/engines/engine_printer.cpp b/src/engines/engine_printer.cpp index f3d6f3b38..cbfde9304 100644 --- a/src/engines/engine_printer.cpp +++ b/src/engines/engine_printer.cpp @@ -6,7 +6,12 @@ #include "utils/colors.h" #include "utils/formatting.h" -#include "framework/specialization_registry.h" +#include "metrics/kerr_schild.h" +#include "metrics/kerr_schild_0.h" +#include "metrics/minkowski.h" +#include "metrics/qkerr_schild.h" +#include "metrics/qspherical.h" +#include "metrics/spherical.h" #include "engines/engine.hpp" @@ -301,12 +306,6 @@ namespace ntt { add_param(report, 4, "Problem generator", "%s", pgen.c_str()); add_param(report, 4, "Engine", "%s", SimEngine(S).to_string()); add_param(report, 4, "Metric", "%s", Metric(M::MetricType).to_string()); -#if SHAPE_ORDER == 0 - add_param(report, 4, "Deposit", "%s", "zigzag"); -#else - add_param(report, 4, "Deposit", "%s", "esirkepov"); - add_param(report, 4, "Interpolation order", "%i", SHAPE_ORDER); -#endif add_param(report, 4, "Timestep [dt]", "%.3e", dt); add_param(report, 4, "Runtime", "%.3e [%d steps]", runtime, max_steps); report += "\n"; @@ -490,11 +489,13 @@ namespace ntt { } } -#define ENGINE_PRINTER(S, M, D) \ - template void Engine>::print_report() const; - - NTT_FOREACH_SPECIALIZATION(ENGINE_PRINTER) - -#undef ENGINE_PRINTER + template void Engine>::print_report() const; + template void Engine>::print_report() const; + template void Engine>::print_report() const; + template void Engine>::print_report() const; + template void Engine>::print_report() const; + template void Engine>::print_report() const; + template void Engine>::print_report() const; + template void Engine>::print_report() const; } // namespace ntt diff --git a/src/engines/engine_run.cpp b/src/engines/engine_run.cpp index 961e39206..65c87e939 100644 --- a/src/engines/engine_run.cpp +++ b/src/engines/engine_run.cpp @@ -3,8 +3,14 @@ #include "arch/traits.h" #include "utils/diag.h" +#include "metrics/kerr_schild.h" +#include "metrics/kerr_schild_0.h" +#include "metrics/minkowski.h" +#include "metrics/qkerr_schild.h" +#include "metrics/qspherical.h" +#include "metrics/spherical.h" + #include "framework/domain/domain.h" -#include "framework/specialization_registry.h" #include "engines/engine.hpp" @@ -137,10 +143,13 @@ namespace ntt { } } -#define ENGINE_RUN(S, M, D) template void Engine>::run(); - - NTT_FOREACH_SPECIALIZATION(ENGINE_RUN) - -#undef ENGINE_RUN + template void Engine>::run(); + template void Engine>::run(); + template void Engine>::run(); + template void Engine>::run(); + template void Engine>::run(); + template void Engine>::run(); + template void Engine>::run(); + template void Engine>::run(); } // namespace ntt diff --git a/src/engines/engine_traits.h b/src/engines/engine_traits.h deleted file mode 100644 index 28a74bea3..000000000 --- a/src/engines/engine_traits.h +++ /dev/null @@ -1,30 +0,0 @@ -// SPDX-License-Identifier: BSD-3-Clause - -#ifndef ENGINES_ENGINE_TRAITS_H -#define ENGINES_ENGINE_TRAITS_H - -#include "enums.h" - -#include "engines/grpic.hpp" -#include "engines/srpic.hpp" - -namespace ntt { - - template - struct EngineSelector; - - template <> - struct EngineSelector { - template - using type = SRPICEngine; - }; - - template <> - struct EngineSelector { - template - using type = GRPICEngine; - }; - -} // namespace ntt - -#endif // ENGINES_ENGINE_TRAITS_H diff --git a/src/engines/grpic.hpp b/src/engines/grpic.hpp index 4d9d89d92..35e544a99 100644 --- a/src/engines/grpic.hpp +++ b/src/engines/grpic.hpp @@ -94,9 +94,9 @@ namespace ntt { void step_forward(timer::Timers& timers, domain_t& dom) override { const auto fieldsolver_enabled = m_params.template get( - "algorithms.fieldsolver.enable"); + "algorithms.toggles.fieldsolver"); const auto deposit_enabled = m_params.template get( - "algorithms.deposit.enable"); + "algorithms.toggles.deposit"); const auto clear_interval = m_params.template get( "particles.clear_interval"); diff --git a/src/engines/srpic.hpp b/src/engines/srpic.hpp index 0dad75e91..c2b9894d6 100644 --- a/src/engines/srpic.hpp +++ b/src/engines/srpic.hpp @@ -78,9 +78,9 @@ namespace ntt { void step_forward(timer::Timers& timers, domain_t& dom) override { const auto fieldsolver_enabled = m_params.template get( - "algorithms.fieldsolver.enable"); + "algorithms.toggles.fieldsolver"); const auto deposit_enabled = m_params.template get( - "algorithms.deposit.enable"); + "algorithms.toggles.deposit"); const auto clear_interval = m_params.template get( "particles.clear_interval"); @@ -348,7 +348,7 @@ namespace ntt { : ZERO; // cooling const auto has_synchrotron = (cooling == Cooling::SYNCHROTRON); - const auto has_compton = (cooling == Cooling::COMPTON); + const auto has_compton = (cooling == Cooling::COMPTON); const auto sync_grad = has_synchrotron ? m_params.template get( "algorithms.synchrotron.gamma_rad") @@ -359,16 +359,17 @@ namespace ntt { "scales.omegaB0") / (SQR(sync_grad) * species.mass()) : ZERO; - const auto comp_grad = has_compton ? m_params.template get( - "algorithms.compton.gamma_rad") - : ZERO; - const auto comp_coeff = has_compton ? (real_t)(0.1) * dt * - m_params.template get( - "scales.omegaB0") / - (SQR(comp_grad) * species.mass()) - : ZERO; + const auto comp_grad = has_compton + ? m_params.template get( + "algorithms.compton.gamma_rad") + : ZERO; + const auto comp_coeff = has_compton + ? (real_t)(0.1) * dt * + m_params.template get( + "scales.omegaB0") / (SQR(comp_grad) * species.mass()) + : ZERO; // toggle to indicate whether pgen defines the external force - bool has_extforce = false; + bool has_extforce = false; if constexpr (traits::has_member::value) { has_extforce = true; // toggle to indicate whether the ext force applies to current species @@ -518,30 +519,9 @@ namespace ntt { } } - template - void deposit_with(const Particles& species, - const M& metric, - const scatter_ndfield_t& scatter_cur, - real_t dt) { - // clang-format off - Kokkos::parallel_for("CurrentsDeposit", - species.rangeActiveParticles(), - kernel::DepositCurrents_kernel( - scatter_cur, - species.i1, species.i2, species.i3, - species.i1_prev, species.i2_prev, species.i3_prev, - species.dx1, species.dx2, species.dx3, - species.dx1_prev, species.dx2_prev, species.dx3_prev, - species.ux1, species.ux2, species.ux3, - species.phi, species.weight, species.tag, - metric, (real_t)(species.charge()), dt)); - // clang-format on - } - void CurrentsDeposit(domain_t& domain) { auto scatter_cur = Kokkos::Experimental::create_scatter_view( domain.fields.cur); - auto shape_order = m_params.template get("algorithms.deposit.order"); for (auto& species : domain.species) { if ((species.pusher() == PrtlPusher::NONE) or (species.npart() == 0) or cmp::AlmostZero_host(species.charge())) { @@ -554,8 +534,20 @@ namespace ntt { species.npart(), (double)species.charge()), HERE); - - deposit_with(species, domain.mesh.metric, scatter_cur, dt); + // clang-format off + Kokkos::parallel_for("CurrentsDeposit", + species.rangeActiveParticles(), + kernel::DepositCurrents_kernel( + scatter_cur, + species.i1, species.i2, species.i3, + species.i1_prev, species.i2_prev, species.i3_prev, + species.dx1, species.dx2, species.dx3, + species.dx1_prev, species.dx2_prev, species.dx3_prev, + species.ux1, species.ux2, species.ux3, + species.phi, species.weight, species.tag, + domain.mesh.metric, + (real_t)(species.charge()), dt)); + // clang-format on } Kokkos::Experimental::contribute(domain.fields.cur, scatter_cur); } @@ -1291,7 +1283,7 @@ namespace ntt { } const auto maxwellian = arch::Maxwellian { domain.mesh.metric, - domain.random_pool(), + domain.random_pool, temp }; if (dim == in::x1) { diff --git a/src/entity.cpp b/src/entity.cpp index 8d88e2654..79b2f1335 100644 --- a/src/entity.cpp +++ b/src/entity.cpp @@ -3,13 +3,18 @@ #include "arch/traits.h" #include "utils/error.h" -#include "engines/engine_traits.h" +#include "metrics/kerr_schild.h" +#include "metrics/kerr_schild_0.h" +#include "metrics/minkowski.h" +#include "metrics/qkerr_schild.h" +#include "metrics/qspherical.h" +#include "metrics/spherical.h" + #include "framework/simulation.h" -#include "framework/specialization_registry.h" +#include "engines/grpic.hpp" +#include "engines/srpic.hpp" #include "pgen.hpp" - -#include template class M, Dimension D> static constexpr bool should_compile { @@ -18,54 +23,95 @@ static constexpr bool should_compile { traits::check_compatibility::value(user::PGen>::dimensions) }; -template class M, Dimension D> -void dispatch_engine(ntt::Simulation& sim) { - if constexpr (S == SimEngine::SRPIC) { - sim.run::template type, M, D>(); - } else if constexpr (S == SimEngine::GRPIC) { - sim.run::template type, M, D>(); - } else { - static_assert( - traits::always_false>::value, - "Unsupported engine"); +template