diff --git a/CLAUDE.md b/CLAUDE.md index b9632690..dd36d3fa 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -49,6 +49,7 @@ entity │ └── entity.cpp # main entry-point ├── tests # unit tests for all submodules ├── .clang-format # code formatting guidelines for clang-format +├── .clang-tidy # configurations for the clang-tidy ├── .gitattributes ├── .gitignore ├── .gitmodules @@ -59,9 +60,7 @@ entity ├── CODE_OF_CONDUCT.md ├── LICENSE ├── README.md -├── conda-entity-nompi.sh ├── dependencies.py # deployment scripts on various machines -├── docker-compose.yml └── input.example.toml # most complete toml file with all possible input options ``` diff --git a/pgens/examples/README.md b/pgens/examples/README.md index 93be01d4..9309c471 100644 --- a/pgens/examples/README.md +++ b/pgens/examples/README.md @@ -6,11 +6,11 @@ Problem generators in this directory are just examples demonstrating how to use image -- `custom_spatial_distribution`: example of using the non-uniform plasma injector with a custom spatial distirbution +- `custom_spatial_distribution`: example of using the non-uniform plasma injector with a custom spatial distirbution. https://github.com/user-attachments/assets/bea0c290-e7e4-4ec7-b360-68ce76beab5b -- `match_fix_field_boundaries`: example of setting matching and/or fixed (coordinate-independent) field boundaries for the electromagnetic fields +- `match_fix_field_boundaries`: example of setting matching and/or fixed (coordinate-independent) field boundaries for the electromagnetic fields. https://github.com/user-attachments/assets/a1b9ea22-34ce-474f-a9b0-49789a2e52b3 @@ -22,13 +22,23 @@ Problem generators in this directory are just examples demonstrating how to use image -- `atmosphere`: setting up a gravitationally bound "atmosphere" with a constant particle replenisher (the plot also highlights the importance of having a constant replenisher and a gravity force acting on the particles; both are enabled by default when using the `ATMOSPHERE` particle boundary conditions) +- `atmosphere`: setting up a gravitationally bound "atmosphere" with a constant particle replenisher (the plot also highlights the importance of having a constant replenisher and a gravity force acting on the particles; both are enabled by default when using the `ATMOSPHERE` particle boundary conditions). atmosphere -- `replenish_injector`: demonstration of how to use the replenish injector to periodically inject new plasma to a target density (both uniform and non-uniform) in `CustomPostStep` +- `replenish_injector`: demonstration of how to use the replenish injector to periodically inject new plasma to a target density (both uniform and non-uniform) in `CustomPostStep`. https://github.com/user-attachments/assets/4955f05e-4795-439f-a4a5-8196e42e987b https://github.com/user-attachments/assets/ebc303ed-7b21-4f97-9822-d81124fdf962 +- `piston`: implementing a moving piston which continuously pushes the plasma. + + https://github.com/user-attachments/assets/7d2dfc39-5de4-49cb-9d55-be8700cbcaee + +- `moving_window`: demonstrates the domain being updated to follow certain features in your simulation in one of the direction (jerkiness is due to updates taking place every `n`-th timestep, if you limit plotting to a specific physical interval, this goes away). + + https://github.com/user-attachments/assets/00edb648-af80-4d37-b1c2-c49872227d47 + + + diff --git a/pgens/examples/moving_window/moving_window.mp4 b/pgens/examples/moving_window/moving_window.mp4 new file mode 100644 index 00000000..373af2e4 --- /dev/null +++ b/pgens/examples/moving_window/moving_window.mp4 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:280520498f8064d5e4ef406accc9f3c7919f9f564fa76881f481766e70f99928 +size 3041327 diff --git a/pgens/examples/moving_window/moving_window.py b/pgens/examples/moving_window/moving_window.py new file mode 100644 index 00000000..234dc789 --- /dev/null +++ b/pgens/examples/moving_window/moving_window.py @@ -0,0 +1,40 @@ +import nt2 +import matplotlib.pyplot as plt + + +def plot(t, data): + plt.rcParams["text.usetex"] = True + plt.rcParams["font.family"] = "serif" + plt.rcParams["figure.dpi"] = 250 + ax = plt.gca() + data.fields.N_1_2.sel(t=t, method="nearest").plot( + vmin=0, vmax=1.5, ax=ax, cmap="inferno" + ) + plt.gcf().axes[1].set_ylabel(r"$n_-+n_+$") + + ax.annotate( + "", + xy=(0.28, 0.08), + xycoords="axes fraction", + xytext=(0.08, 0.08), + textcoords="axes fraction", + arrowprops=dict(arrowstyle="-|>", color="#e05c2e", lw=1.8), + ) + ax.text( + 0.095, + 0.095, + r"$v_0$", + transform=ax.transAxes, + color="#e05c2e", + fontsize=12, + va="center", + ) + ax.set( + xlabel=r"$x$", + ylabel=r"$y$", + title=rf"$t={{{t:.2f}}}$", + ) + + +moving_window = nt2.Data("moving_window") +moving_window.makeMovie(plot, framerate=10) diff --git a/pgens/examples/moving_window/moving_window.toml b/pgens/examples/moving_window/moving_window.toml new file mode 100644 index 00000000..1e07cf5e --- /dev/null +++ b/pgens/examples/moving_window/moving_window.toml @@ -0,0 +1,57 @@ +[simulation] + name = "moving_window" + engine = "srpic" + runtime = 2.0 + +[grid] + resolution = [128, 128] + extent = [[0.0, 1.0], [0.0, 1.0]] + + [grid.metric] + metric = "minkowski" + + [grid.boundaries] + fields = [["PERIODIC"], ["PERIODIC"]] + particles = [["ABSORB", "ABSORB"], ["PERIODIC"]] + +[scales] + larmor0 = 0.01 + skindepth0 = 0.1 + +[algorithms] + + [algorithms.timestep] + CFL = 0.5 + +[particles] + ppc0 = 8.0 + + [[particles.species]] + label = "e+" + mass = 1.0 + charge = 1.0 + maxnpart = 1e6 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 1e6 + +[setup] + window_velocity = 0.2 + +[output] + interval = 10 + + [output.fields] + quantities = ["N_1_2"] + + [output.particles] + enable = false + + [output.spectra] + enable = false + +[checkpoint] + keep = 0 diff --git a/pgens/examples/moving_window/pgen.hpp b/pgens/examples/moving_window/pgen.hpp new file mode 100644 index 00000000..58d4ff10 --- /dev/null +++ b/pgens/examples/moving_window/pgen.hpp @@ -0,0 +1,114 @@ +#ifndef PROBLEM_GENERATOR_H +#define PROBLEM_GENERATOR_H + +#include "enums.h" +#include "global.h" + +#include "arch/kokkos_aliases.h" +#include "traits/pgen.h" + +#include "archetypes/energy_dist.h" +#include "archetypes/moving_window.h" +#include "archetypes/particle_injector.h" +#include "framework/domain/domain.h" +#include "framework/domain/metadomain.h" +#include "framework/parameters/parameters.h" + +namespace user { + using namespace ntt; + + template + struct NonUniformDensity { + Inline auto operator()(const coord_t& x_Ph) const -> real_t { + // example of a non-uniform target density that peaks at the center of the domain + real_t r2 { ZERO }; + for (auto d = 0u; d < D; ++d) { + r2 += SQR(x_Ph[d] - 0.5); + } + return math::exp(-r2 / SQR(static_cast(0.1))); + // ^ + // characteristic width of the density profile + } + }; + + template + struct PGen { + static constexpr auto engines { + ::traits::pgen::compatible_with {} + }; + static constexpr auto metrics { + ::traits::pgen::compatible_with {} + }; + static constexpr auto dimensions { + ::traits::pgen::compatible_with {} + }; + const SimulationParams& params; + Metadomain& metadomain; + + struct MovingWindow { + int pos_i { -1 }; + int init_pos_i { -1 }; + real_t pos_di { ZERO }; + real_t vel_Cd { ZERO }; + + void init(const M& metric, real_t global_x, real_t velocity) { + const auto pos_Cd = metric.template convert<1, Crd::Ph, Crd::Cd>(global_x); + pos_i = static_cast(pos_Cd + 1) - 1; + init_pos_i = pos_i; + pos_di = pos_Cd - static_cast(pos_i); + vel_Cd = metric.template transform<1, Idx::XYZ, Idx::U>({}, velocity); + } + + void update(real_t dt, + int shift, + Metadomain& metadomain, + Domain& domain) { + pos_di += vel_Cd * dt; + pos_i += static_cast(pos_di >= ONE); + pos_di -= static_cast(pos_di >= ONE); + if ((pos_i - init_pos_i) >= shift) { + // move the window and all fields and particles in it + arch::MoveWindow(domain, metadomain, shift); + // update window index for next update + pos_i -= shift; + } + } + }; + + MovingWindow moving_window; + const real_t dt; + + PGen(const SimulationParams& p, Metadomain& m) + : params { p } + , metadomain { m } + , dt { params.template get("algorithms.timestep.dt") } {} + + void InitPrtls(Domain& domain) { + const auto window_velocity = params.template get( + "setup.window_velocity", + ZERO); + moving_window.init(domain.mesh.metric, ZERO, window_velocity); + // inject particles + const auto energy_dist = arch::energy_dist::Maxwellian( + domain.random_pool(), + 0.001, + { window_velocity, ZERO, ZERO }); + + const auto density_profile = NonUniformDensity {}; + arch::InjectNonUniform( + params, + domain, + { 1, 2 }, + { energy_dist, energy_dist }, + density_profile, + ONE); + } + + void CustomPostStep(timestep_t /*step*/, simtime_t /*time*/, Domain& domain) { + moving_window.update(dt, N_GHOSTS, metadomain, domain); + } + }; + +} // namespace user + +#endif diff --git a/pgens/examples/piston/pgen.hpp b/pgens/examples/piston/pgen.hpp new file mode 100644 index 00000000..fe8fb6d4 --- /dev/null +++ b/pgens/examples/piston/pgen.hpp @@ -0,0 +1,109 @@ +#ifndef PROBLEM_GENERATOR_H +#define PROBLEM_GENERATOR_H + +#include "enums.h" +#include "global.h" + +#include "arch/kokkos_aliases.h" +#include "traits/pgen.h" + +#include "archetypes/energy_dist.h" +#include "archetypes/particle_injector.h" +#include "archetypes/piston.h" +#include "framework/domain/domain.h" +#include "framework/domain/metadomain.h" +#include "framework/parameters/parameters.h" + +namespace user { + using namespace ntt; + + template + struct NonUniformDensity { + Inline auto operator()(const coord_t& x_Ph) const -> real_t { + // example of a non-uniform target density that peaks at the center of the domain + real_t r2 { ZERO }; + for (auto d = 0u; d < D; ++d) { + r2 += SQR(x_Ph[d] - 0.5); + } + return math::exp(-r2 / SQR(static_cast(0.2))); + // ^ + // characteristic width of the density profile + } + }; + + template + struct PGen { + + static constexpr auto engines { + ::traits::pgen::compatible_with {} + }; + static constexpr auto metrics { + ::traits::pgen::compatible_with {} + }; + static constexpr auto dimensions { + ::traits::pgen::compatible_with {} + }; + + const SimulationParams& params; + Metadomain& metadomain; + const real_t piston_velocity; + + PGen(const SimulationParams& p, Metadomain& m) + : params { p } + , metadomain { m } + , piston_velocity { params.template get("setup.piston_velocity", + 0.0) } {} + + void InitPrtls(Domain& domain) { + // inject particles + const auto energy_dist = arch::energy_dist::Maxwellian( + domain.random_pool(), + 0.2); // <-- target temperature for injection + + const auto density_profile = NonUniformDensity {}; + arch::InjectNonUniform( + params, + domain, + { 1, 2 }, + { energy_dist, energy_dist }, + density_profile, + ONE); + } + + struct CustomPrtlUpdate { + real_t x_piston; // current position of piston + real_t v_piston; + real_t global_xmax; + bool is_left; + bool massive; + + Inline void operator()(index_t p, + const kernel::sr::PusherContext& ctx, + const kernel::sr::PusherBoundaries&, + const kernel::PusherArrays& particles, + const M& metric) const { + + real_t piston_pos; + if (x_piston > global_xmax) { + // piston has moved beyond the domain, set position to a large value so that no particles interact with it + piston_pos = global_xmax; + } else { + piston_pos = x_piston; + } + + arch::Piston(p, ctx.dt, particles, metric, piston_pos, v_piston, massive); + } + }; + + template + auto CustomParticleUpdate(simtime_t time, spidx_t /*sp*/, DOM& /*domain*/) const { + return CustomPrtlUpdate { piston_velocity * static_cast(time), + piston_velocity, + metadomain.mesh().extent(in::x1).second, + true, + true }; + }; + }; +} // namespace user + +#endif diff --git a/pgens/examples/piston/piston.mp4 b/pgens/examples/piston/piston.mp4 new file mode 100644 index 00000000..a4fe967f --- /dev/null +++ b/pgens/examples/piston/piston.mp4 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e90c9191bfb2c53e0be215ec06fceff673133e1381b0da436f742256619e6022 +size 10899641 diff --git a/pgens/examples/piston/piston.py b/pgens/examples/piston/piston.py new file mode 100644 index 00000000..2f8c0315 --- /dev/null +++ b/pgens/examples/piston/piston.py @@ -0,0 +1,25 @@ +import nt2 +import matplotlib.pyplot as plt + + +def plot(t, data): + plt.rcParams["text.usetex"] = True + plt.rcParams["font.family"] = "serif" + plt.rcParams["figure.dpi"] = 250 + ax = plt.gca() + data.fields.N_1_2.sel(t=t, method="nearest").plot( + vmin=0, vmax=1.5, ax=ax, cmap="inferno" + ) + plt.gcf().axes[1].set_ylabel(r"$n_-+n_+$") + + ax.set( + xlabel=r"$x$", + ylabel=r"$y$", + title=rf"$t={{{t:.2f}}}$", + xlim=(0, 1), + ylim=(0, 1), + ) + + +piston = nt2.Data("piston") +piston.makeMovie(plot, framerate=10) diff --git a/pgens/examples/piston/piston.toml b/pgens/examples/piston/piston.toml new file mode 100644 index 00000000..0b5273af --- /dev/null +++ b/pgens/examples/piston/piston.toml @@ -0,0 +1,57 @@ +[simulation] + name = "piston" + engine = "srpic" + runtime = 2.0 + +[grid] + resolution = [128, 128] + extent = [[0.0, 1.0], [0.0, 1.0]] + + [grid.metric] + metric = "minkowski" + + [grid.boundaries] + fields = [["PERIODIC"], ["PERIODIC"]] + particles = [["ABSORB", "ABSORB"], ["PERIODIC"]] + +[scales] + larmor0 = 0.01 + skindepth0 = 0.1 + +[algorithms] + + [algorithms.timestep] + CFL = 0.5 + +[particles] + ppc0 = 8.0 + + [[particles.species]] + label = "e+" + mass = 1.0 + charge = 1.0 + maxnpart = 1e6 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 1e6 + +[setup] + piston_velocity = 0.5 + +[output] + interval = 10 + + [output.fields] + quantities = ["N_1_2"] + + [output.particles] + enable = false + + [output.spectra] + enable = false + +[checkpoint] + keep = 0 diff --git a/pgens/examples/replenish_injector/pgen.hpp b/pgens/examples/replenish_injector/pgen.hpp index fb46eaca..8081f8f0 100644 --- a/pgens/examples/replenish_injector/pgen.hpp +++ b/pgens/examples/replenish_injector/pgen.hpp @@ -4,13 +4,14 @@ #include "enums.h" #include "global.h" +#include "arch/kokkos_aliases.h" #include "traits/pgen.h" #include "archetypes/energy_dist.h" #include "archetypes/particle_injector.h" #include "archetypes/spatial_dist.h" +#include "archetypes/utils.h" #include "framework/domain/metadomain.h" -#include "kernels/particle_moments.hpp" namespace user { using namespace ntt; @@ -34,7 +35,7 @@ namespace user { for (auto d = 0u; d < D; ++d) { r2 += SQR(x_Ph[d] - HALF); } - return std::exp(-r2 / SQR(static_cast(0.2))); + return math::exp(-r2 / SQR(static_cast(0.2))); // ^ // characteristic width of the density profile } @@ -70,46 +71,12 @@ namespace user { } // perform replenishment and injection every 100 timesteps - { // compute density of species #1 and #2 - - // saves the density to domain.fields.buff(:,:,:,0) - const auto ni2 = domain.mesh.n_active(in::x2); - const auto inv_n0 = ONE / params.template get("scales.n0"); - - auto scatter_buff = Kokkos::Experimental::create_scatter_view( - domain.fields.buff); - Kokkos::deep_copy(domain.fields.buff, ZERO); - for (const auto sp : std::vector { 1, 2 }) { - const auto& prtl_spec = domain.species[sp - 1]; - Kokkos::parallel_for("ComputeDensity", - prtl_spec.rangeActiveParticles(), - kernel::ParticleMoments_kernel( - {}, - scatter_buff, - 0u, - prtl_spec.i1, - prtl_spec.i2, - prtl_spec.i3, - prtl_spec.dx1, - prtl_spec.dx2, - prtl_spec.dx3, - prtl_spec.ux1, - prtl_spec.ux2, - prtl_spec.ux3, - prtl_spec.phi, - prtl_spec.weight, - prtl_spec.tag, - prtl_spec.mass(), - prtl_spec.charge(), - false, - domain.mesh.metric, - domain.mesh.flds_bc(), - ni2, - inv_n0, - 0u)); - } - Kokkos::Experimental::contribute(domain.fields.buff, scatter_buff); - } + // compute the density of species #1 and #2 + // and save in the field buffer (index 0) + arch::ComputeMomentWithSpecies(params, + domain, + { 1u, 2u }, + domain.fields.buff); const auto energy_dist = arch::energy_dist::Maxwellian( domain.random_pool(), diff --git a/src/archetypes/moving_window.h b/src/archetypes/moving_window.h new file mode 100644 index 00000000..b4ef2155 --- /dev/null +++ b/src/archetypes/moving_window.h @@ -0,0 +1,278 @@ +/** + * @file archetypes/moving_window.h + * @brief Moving window functions for implementing the moving window in the CustomPostStep + * @implements + * - arch::MoveWindow<> -> void + * @namespaces: + * - arch:: + */ + +#ifndef ARCHETYPES_WINDOW_H +#define ARCHETYPES_WINDOW_H + +#include "enums.h" +#include "global.h" + +#include "framework/domain/domain.h" +#include "framework/domain/metadomain.h" + +namespace arch { + + 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, + ndfield_t backup_Fld, + const int window_shift, + BCTags tags) + : Fld { Fld } + , backup_Fld { backup_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) = 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) = 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( + 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) = 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) = 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) = 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) = 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 { + 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) = 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) = 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) = 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) = 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) = 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) = 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 { + raise::KernelError( + HERE, + "FieldShift_kernel: 3D implementation called for D != 3"); + } + } + }; + + /** + * @brief Updates particle position and fields in the moving window. + + */ + template + Inline void MoveWindow(Domain& domain, + Metadomain& metadomain, + int window_shift) { + + /* + 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); + } + } + + // synch ghost zones after moving the window + metadomain.CommunicateFields(domain, Comm::E | Comm::B); + // communicate particles after moving + metadomain.CommunicateParticles(domain); + + metadomain.ShiftByCells(window_shift, o); + } +} // namespace arch + +#endif // ARCHETYPES_UTILS_H diff --git a/src/archetypes/piston.h b/src/archetypes/piston.h index 4cb33bf8..f31dbdff 100644 --- a/src/archetypes/piston.h +++ b/src/archetypes/piston.h @@ -38,7 +38,39 @@ namespace arch { /** - * @brief Updates particle position and velocity if it reflects off a moiving + * @brief Checks whether a particle reflects off a moving piston, called after particle has been moved by regular pusher + * + * @param p Index of particle + * @param dt Timestep + * @param particles Particle data arrays + * @param metric Metric object for coordinate transformations + * @param piston_position Position of the piston at the start of timestep in global coordinates + * @param piston_v Velocity of piston at current timestep + * @param is_left Is piston on the left side of the box or right side of the box + */ + template + Inline bool CrossesPiston(index_t p, + real_t dt, + const kernel::PusherArrays& particles, + const M& metric, + real_t piston_position, + real_t piston_v, + bool is_left) { + const real_t x1_Cd = i_di_to_Xi(particles.i1(p), particles.dx1(p)); + // x1_Cd_wallmove is not the actual particle coordinate + // it is particle position minus how much the wall has moved in this + // timestep This is a computational trick + const real_t x1_Cd_wallmove = + x1_Cd - metric.template transform<1, Idx::XYZ, Idx::U>({}, piston_v) * dt; + const real_t x1_Ph_wallmove = metric.template convert<1, Crd::Cd, Crd::Ph>( + x1_Cd_wallmove); + + return is_left ? piston_position > x1_Ph_wallmove + : piston_position < x1_Ph_wallmove; + } + + /** + * @brief Updates particle position and velocity if it reflects off a moving * piston, called to correct patticle position * * @param p Index of particle @@ -50,13 +82,18 @@ namespace arch { * @param massive Whether the particle is massive or massless (e.g. photon) */ template - Inline void PistonUpdate(index_t p, - real_t dt, - const kernel::PusherArrays& particles, - const M& metric, - real_t piston_position, - real_t piston_v, - bool massive) { + Inline void Piston(index_t p, + real_t dt, + const kernel::PusherArrays& particles, + const M& metric, + real_t piston_position, + real_t piston_v, + bool massive) { + + // check if particle actually crosses the piston, if not return + if (!CrossesPiston(p, dt, particles, metric, piston_position, piston_v, true)) { + return; + } // step 1: calculate the particle 3 velocity const real_t gamma_p { massive ? U2GAMMA(particles.ux1(p), particles.ux2(p), particles.ux3(p)) @@ -93,23 +130,20 @@ namespace arch { TWO * piston_v * gamma_p); const real_t remaining_dt_inv_energy { - massive - ? (remaining_dt / - math::sqrt( - ONE + U2GAMMA(particles.ux1(p), particles.ux2(p), particles.ux3(p)))) - : (remaining_dt / math::sqrt(SQR(particles.ux1(p)) + SQR(particles.ux2(p)) + - SQR(particles.ux3(p)))) + massive ? (remaining_dt / + U2GAMMA(particles.ux1(p), particles.ux2(p), particles.ux3(p))) + : (remaining_dt / + NORM(particles.ux1(p), particles.ux2(p), particles.ux3(p))) }; // define piston integer and fractional coordinate - int i_w_coll = i_w; prtldx_t dx_w_coll = dx_w + metric.template transform<1, Idx::XYZ, Idx::U>( {}, piston_v) * dt_to_piston; - i_w_coll += static_cast(dx_w_coll >= ONE) - - static_cast(dx_w_coll < ZERO); + i_w_coll += static_cast(dx_w_coll >= ONE) - + static_cast(dx_w_coll < ZERO); dx_w_coll -= static_cast(dx_w_coll >= ONE); dx_w_coll += static_cast(dx_w_coll < ZERO); @@ -119,44 +153,12 @@ namespace arch { remaining_dt_inv_energy + dx_w_coll; - particles.i1(p) += static_cast(particles.dx1(p) >= ONE) - - static_cast(particles.dx1(p) < ZERO); + particles.i1(p) += static_cast(particles.dx1(p) >= ONE) - + static_cast(particles.dx1(p) < ZERO); particles.dx1(p) -= static_cast(particles.dx1(p) >= ONE); particles.dx1(p) += static_cast(particles.dx1(p) < ZERO); } - /** - * @brief Checks whether a particle reflects off a moving piston, called after particle has been moved by regular pusher - * - * @param p Index of particle - * @param dt Timestep - * @param particles Particle data arrays - * @param metric Metric object for coordinate transformations - * @param piston_position Position of the piston at the start of timestep in global coordinates - * @param piston_v Velocity of piston at current timestep - * @param is_left Is piston on the left side of the box or right side of the box - */ - template - Inline bool CrossesPiston(index_t p, - real_t dt, - const kernel::PusherArrays& particles, - const M& metric, - real_t piston_position, - real_t piston_v, - bool is_left) { - const real_t x1_Cd = i_di_to_Xi(particles.i1(p), particles.dx1(p)); - // x1_Cd_wallmove is not the actual particle coordinate - // it is particle position minus how much the wall has moved in this - // timestep This is a computational trick - const real_t x1_Cd_wallmove = - x1_Cd - metric.template transform<1, Idx::XYZ, Idx::U>({}, piston_v) * dt; - const real_t x1_Ph_wallmove = metric.template convert<1, Crd::Cd, Crd::Ph>( - x1_Cd_wallmove); - - return is_left ? piston_position > x1_Ph_wallmove - : piston_position < x1_Ph_wallmove; - } - } // namespace arch #undef from_Xi_to_i_di diff --git a/src/framework/CMakeLists.txt b/src/framework/CMakeLists.txt index 371b4be0..5343f3d8 100644 --- a/src/framework/CMakeLists.txt +++ b/src/framework/CMakeLists.txt @@ -18,6 +18,7 @@ # * domain/metadomain_chckpt.cpp # * domain/metadomain_stats.cpp # * domain/metadomain_io.cpp +# * domain/metadomain_reshape.cpp # * containers/particles.cpp # * containers/particles_comm.cpp # * containers/particles_io.cpp @@ -57,6 +58,7 @@ set(SOURCES ${SRC_DIR}/domain/metadomain_comm.cpp ${SRC_DIR}/domain/metadomain_sort.cpp ${SRC_DIR}/domain/metadomain_stats.cpp + ${SRC_DIR}/domain/metadomain_reshape.cpp ${SRC_DIR}/containers/particles.cpp ${SRC_DIR}/containers/particles_sort.cpp ${SRC_DIR}/containers/fields.cpp) diff --git a/src/framework/domain/mesh.h b/src/framework/domain/mesh.h index 07d1d150..aa2eea72 100644 --- a/src/framework/domain/mesh.h +++ b/src/framework/domain/mesh.h @@ -45,7 +45,8 @@ namespace ntt { const boundaries_t& ext, const std::map& metric_params) : Grid { res, ext } - , metric { res, ext, metric_params } {} + , metric { res, ext, metric_params } + , m_metric_params_raw { metric_params } {} Mesh(const std::vector& res, const boundaries_t& ext, @@ -53,10 +54,17 @@ namespace ntt { const boundaries_t& flds_bc, const boundaries_t& prtl_bc) : Grid { res, ext, flds_bc, prtl_bc } - , metric { res, ext, metric_params } {} + , metric { res, ext, metric_params } + , m_metric_params_raw { metric_params } {} ~Mesh() = default; + void set_extent(const boundaries_t& new_extent) { + m_extent = new_extent; + metric.~M(); + new (&metric) M { this->m_resolution, new_extent, m_metric_params_raw }; + } + /** * @brief Get the intersection of the mesh with a box * @param box physical extent @@ -193,6 +201,9 @@ namespace ntt { return range; } + + private: + std::map m_metric_params_raw; }; } // namespace ntt diff --git a/src/framework/domain/metadomain.h b/src/framework/domain/metadomain.h index 9c67a25a..c0a11efa 100644 --- a/src/framework/domain/metadomain.h +++ b/src/framework/domain/metadomain.h @@ -10,6 +10,7 @@ * - metadomain_chckpt.cpp * - metadomain_io.cpp * - metadomain_stats.cpp + * - metadomain_reshape.cpp * @namespaces: * - ntt:: * @macros: @@ -95,6 +96,9 @@ namespace ntt { void SynchronizeFields(Domain&, CommTags, const range_tuple_t& = { 0, 0 }) const; +#if defined(MPI_ENABLED) && defined(OUTPUT_ENABLED) + void CommunicateVectorPotential(unsigned short); +#endif void CommunicateParticles(Domain&) const; void SortParticles(simtime_t, timestep_t, @@ -125,6 +129,11 @@ namespace ntt { ~Metadomain() = default; + /* domain update-related ------------------------------------------------ */ + void ShiftByCells(int, in = in::x1) + requires(CartesianMetricClass); + + /* output-related ------------------------------------------------------- */ #if defined(OUTPUT_ENABLED) void InitWriter(adios2::ADIOS*, const SimulationParams&); auto Write(const SimulationParams&, @@ -146,6 +155,8 @@ namespace ntt { simtime_t) -> bool; void ContinueFromCheckpoint(adios2::ADIOS*, const SimulationParams&); + void redecomposeFromCheckpoint(const std::vector>&, + const std::vector>&); #endif void InitStatsWriter(const SimulationParams&, bool); @@ -278,9 +289,6 @@ namespace ntt { #if defined(OUTPUT_ENABLED) out::Writer g_writer; checkpoint::Writer g_checkpoint_writer; - #if defined(MPI_ENABLED) - void CommunicateVectorPotential(unsigned short); - #endif #endif #if defined(MPI_ENABLED) diff --git a/src/framework/domain/metadomain_chckpt.cpp b/src/framework/domain/metadomain_chckpt.cpp index 46aa32b8..e4622f92 100644 --- a/src/framework/domain/metadomain_chckpt.cpp +++ b/src/framework/domain/metadomain_chckpt.cpp @@ -10,12 +10,15 @@ #include "framework/parameters/parameters.h" #include "framework/specialization_registry.h" #include "output/checkpoint.h" +#include "output/utils/readers.h" +#include "output/utils/writers.h" #if defined(MPI_ENABLED) #include #endif #include +#include #include #include @@ -68,6 +71,23 @@ namespace ntt { for (const auto& species : local_domain->species) { species.CheckpointDeclare(g_checkpoint_writer.io()); } + for (auto d { 0u }; d < M::Dim; ++d) { + g_checkpoint_writer.io().DefineVariable( + fmt::format("subdomain_x%d_min", d + 1), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + g_checkpoint_writer.io().DefineVariable( + fmt::format("subdomain_x%d_max", d + 1), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + g_checkpoint_writer.io().DefineVariable( + fmt::format("subdomain_nx%d", d + 1), + { adios2::UnknownDim }, + { adios2::UnknownDim }, + { adios2::UnknownDim }); + } } } @@ -111,12 +131,117 @@ namespace ntt { dom_tot, dom_offset); } + for (auto d { 0u }; d < M::Dim; ++d) { + out::WriteVariable(g_checkpoint_writer.io(), + g_checkpoint_writer.writer(), + fmt::format("subdomain_x%d_min", d + 1), + local_domain->mesh.extent()[d].first, + dom_tot, + dom_offset); + out::WriteVariable(g_checkpoint_writer.io(), + g_checkpoint_writer.writer(), + fmt::format("subdomain_x%d_max", d + 1), + local_domain->mesh.extent()[d].second, + dom_tot, + dom_offset); + out::WriteVariable(g_checkpoint_writer.io(), + g_checkpoint_writer.writer(), + fmt::format("subdomain_nx%d", d + 1), + local_domain->mesh.n_active()[d], + dom_tot, + dom_offset); + } } g_checkpoint_writer.endSaving(); logger::Checkpoint("Checkpoint written", HERE); return true; } + template + void Metadomain::redecomposeFromCheckpoint( + const std::vector>& dom_ncells, + const std::vector>& dom_extents) { + + // For each dimension: collect ncells per domain-grid position, + // validate total is unchanged, then compute prefix-sum offsets. + std::vector> offset_ncells_per_dom( + g_ndomains, + std::vector(M::Dim, 0)); + + for (auto d { 0u }; d < M::Dim; ++d) { + std::vector ncells_at_pos(g_ndomains_per_dim[d], 0); + for (unsigned int idx { 0 }; idx < g_ndomains; ++idx) { + ncells_at_pos[g_domain_offsets[idx][d]] = dom_ncells[idx][d]; + } + + ncells_t total { 0 }; + for (const auto& n : ncells_at_pos) { + total += n; + } + raise::ErrorIf( + total != g_mesh.n_active()[d], + fmt::format( + "total cells in dim %d changed between checkpoint (%lu) and current (%lu); " + "changing total domain size is not supported", + d + 1, + total, + g_mesh.n_active()[d]), + HERE); + + ncells_t running { 0 }; + std::vector offset_at_pos(g_ndomains_per_dim[d], 0); + for (unsigned int nd { 0 }; nd < g_ndomains_per_dim[d]; ++nd) { + offset_at_pos[nd] = running; + running += ncells_at_pos[nd]; + } + for (unsigned int idx { 0 }; idx < g_ndomains; ++idx) { + offset_ncells_per_dom[idx][d] = offset_at_pos[g_domain_offsets[idx][d]]; + } + } + + g_subdomains.clear(); + for (unsigned int idx { 0 }; idx < g_ndomains; ++idx) { + const auto& l_offset_ndomains = g_domain_offsets[idx]; + const auto& l_ncells = dom_ncells[idx]; + const auto& l_offset_ncells = offset_ncells_per_dom[idx]; + const auto& l_extent = dom_extents[idx]; + +#if defined(MPI_ENABLED) + const auto local = ((int)idx == g_mpi_rank); + if (not local) { + g_subdomains.emplace_back(false, + idx, + l_offset_ndomains, + l_offset_ncells, + l_ncells, + l_extent, + g_metric_params, + g_species_params); + } else { + g_subdomains.emplace_back(idx, + l_offset_ndomains, + l_offset_ncells, + l_ncells, + l_extent, + g_metric_params, + g_species_params); + } + g_subdomains.back().set_mpi_rank(idx); +#else + g_subdomains.emplace_back(idx, + l_offset_ndomains, + l_offset_ncells, + l_ncells, + l_extent, + g_metric_params, + g_species_params); +#endif + } + + redefineNeighbors(); + redefineBoundaries(); + } + template void Metadomain::ContinueFromCheckpoint(adios2::ADIOS* ptr_adios, const SimulationParams& params) { @@ -140,6 +265,60 @@ namespace ntt { #endif reader.BeginStep(); + + // Phase 1: read all subdomain metadata to detect size changes + std::vector> saved_ncells(g_ndomains, + std::vector(M::Dim)); + std::vector> saved_extents(g_ndomains); + boundaries_t global_extent; + for (auto d { 0u }; d < M::Dim; ++d) { + global_extent.emplace_back(std::numeric_limits::max(), + std::numeric_limits::lowest()); + } + + bool needs_reconstruction = false; + for (unsigned int dom_idx { 0 }; dom_idx < g_ndomains; ++dom_idx) { + for (auto d { 0u }; d < M::Dim; ++d) { + real_t x_min, x_max; + out::ReadVariable(io, + reader, + fmt::format("subdomain_x%d_min", d + 1), + x_min, + dom_idx); + out::ReadVariable(io, + reader, + fmt::format("subdomain_x%d_max", d + 1), + x_max, + dom_idx); + saved_extents[dom_idx].emplace_back(x_min, x_max); + global_extent[d].first = std::min(global_extent[d].first, x_min); + global_extent[d].second = std::max(global_extent[d].second, x_max); + + ncells_t nx; + out::ReadVariable(io, + reader, + fmt::format("subdomain_nx%d", d + 1), + nx, + dom_idx); + saved_ncells[dom_idx][d] = nx; + + if (nx != subdomain_ptr(dom_idx)->mesh.n_active()[d]) { + needs_reconstruction = true; + } + } + } + + // Phase 2: update domain structure + if (needs_reconstruction) { + redecomposeFromCheckpoint(saved_ncells, saved_extents); + } else { + for (unsigned int dom_idx { 0 }; dom_idx < g_ndomains; ++dom_idx) { + subdomain_ptr(dom_idx)->mesh.set_extent(saved_extents[dom_idx]); + } + } + g_mesh.set_extent(global_extent); + + // Phase 3: read field and particle data using the (now-correct) domain layout for (const auto local_domain_idx : l_subdomain_indices()) { auto local_domain = subdomain_ptr(local_domain_idx); @@ -154,8 +333,7 @@ namespace ntt { for (auto& species : local_domain->species) { species.CheckpointRead(io, reader, ndomains(), local_domain_idx); } - - } // local subdomain loop + } reader.EndStep(); reader.Close(); @@ -177,7 +355,10 @@ namespace ntt { simtime_t) -> bool; \ template void Metadomain>::ContinueFromCheckpoint( \ adios2::ADIOS*, \ - const SimulationParams&); + const SimulationParams&); \ + template void Metadomain>::redecomposeFromCheckpoint( \ + const std::vector>&, \ + const std::vector>&); NTT_FOREACH_SPECIALIZATION(METADOMAIN_CHECKPOINTS) #undef METADOMAIN_CHECKPOINTS // NOLINTEND(bugprone-macro-parentheses) diff --git a/src/framework/domain/metadomain_reshape.cpp b/src/framework/domain/metadomain_reshape.cpp new file mode 100644 index 00000000..5104bcd6 --- /dev/null +++ b/src/framework/domain/metadomain_reshape.cpp @@ -0,0 +1,53 @@ +#include "enums.h" +#include "global.h" + +#include "traits/metric.h" +#include "utils/error.h" + +#include "framework/domain/mesh.h" +#include "framework/domain/metadomain.h" +#include "framework/specialization_registry.h" + +namespace ntt { + + template + void Metadomain::ShiftByCells(int n_cells, in dir) + requires(CartesianMetricClass) + { + raise::ErrorIf( + n_cells == 0, + "ShiftByCells called with n_cells = 0, no shift will be performed", + HERE); + + { + // update metadomain g_mesh + auto new_extent = g_mesh.extent(); + const auto delta_shift = static_cast(n_cells) * + g_mesh.metric.get_dx(); + new_extent[static_cast(dir)].first += delta_shift; + new_extent[static_cast(dir)].second += delta_shift; + + g_mesh.set_extent(new_extent); + } + + // update individual subdomain meshes + for (auto& subdomain : g_subdomains) { + auto new_extent = subdomain.mesh.extent(); + const auto delta_shift = static_cast(n_cells) * + subdomain.mesh.metric.get_dx(); + new_extent[static_cast(dir)].first += delta_shift; + new_extent[static_cast(dir)].second += delta_shift; + + subdomain.mesh.set_extent(new_extent); + } + } + + // NOLINTBEGIN(bugprone-macro-parentheses) +#define METADOMAIN_SHIFT(S, M, D) \ + template void Metadomain>::ShiftByCells(int n_cells, in dir); + + NTT_FOREACH_CARTESIAN_SPECIALIZATION(METADOMAIN_SHIFT) +#undef METADOMAIN_SHIFT + // NOLINTEND(bugprone-macro-parentheses) + +} // namespace ntt \ No newline at end of file diff --git a/src/framework/specialization_registry.h b/src/framework/specialization_registry.h index d9c1628c..dea5772c 100644 --- a/src/framework/specialization_registry.h +++ b/src/framework/specialization_registry.h @@ -47,6 +47,11 @@ namespace ntt { MACRO(SimEngine::GRPIC, metric::QKerrSchild, Dim::_2D) \ MACRO(SimEngine::GRPIC, metric::KerrSchild0, Dim::_2D) +#define NTT_FOREACH_CARTESIAN_SPECIALIZATION(MACRO) \ + MACRO(SimEngine::SRPIC, metric::Minkowski, Dim::_1D) \ + MACRO(SimEngine::SRPIC, metric::Minkowski, Dim::_2D) \ + MACRO(SimEngine::SRPIC, metric::Minkowski, Dim::_3D) + #define NTT_BUILD_SPECIALIZATION_ENTRY(S, M, D) SpecializationEntry {}, inline constexpr auto kSpecializations = std::tuple { NTT_FOREACH_SPECIALIZATION( diff --git a/src/metrics/minkowski.h b/src/metrics/minkowski.h index 8d6c87e5..3a679b38 100644 --- a/src/metrics/minkowski.h +++ b/src/metrics/minkowski.h @@ -70,6 +70,11 @@ namespace metric { } } + [[nodiscard]] + auto get_dx() const -> real_t { + return dx; + } + /** * minimum effective cell size for a given metric (in physical units) */ diff --git a/tests/framework/CMakeLists.txt b/tests/framework/CMakeLists.txt index ebeadc75..d6dc295a 100644 --- a/tests/framework/CMakeLists.txt +++ b/tests/framework/CMakeLists.txt @@ -56,6 +56,16 @@ else() gen_test(comm-nompi false) endif() +if(${output}) + if(${mpi}) + gen_test(checkpoint-extents true) + gen_test(checkpoint-data true) + else() + gen_test(checkpoint-extents false) + gen_test(checkpoint-data false) + endif() +endif() + # this test is only run manually to ensure ... # ... command line args are working properly ... # ... and that the logging is done correctly diff --git a/tests/framework/checkpoint-data.cpp b/tests/framework/checkpoint-data.cpp new file mode 100644 index 00000000..3099593d --- /dev/null +++ b/tests/framework/checkpoint-data.cpp @@ -0,0 +1,244 @@ +/** + * Tests that field array sizes and particle coordinates are correctly preserved + * through a checkpoint round-trip. + * + * In MPI mode the checkpoint is written from a metadomain whose subdomain + * sizes have been changed via redecomposeFromCheckpoint (simulating load + * balancing from [32,32] to [40,24] cells), then read back into a uniformly + * decomposed metadomain ([32,32]), exercising the reconstruction path in + * ContinueFromCheckpoint. + */ + +#include "enums.h" +#include "global.h" + +#include "utils/comparators.h" +#include "utils/error.h" +#include "utils/formatting.h" + +#include "framework/domain/metadomain.h" +#include "framework/parameters/parameters.h" + +#include "metrics/minkowski.h" + +#include +#include +#include + +#if defined(MPI_ENABLED) + #include +#endif + +#include +#include +#include +#include + +using namespace ntt; +using namespace metric; + +const std::string CHECKPOINT_DIR = "test_checkpoint_data"; + +void cleanup() { + std::filesystem::remove_all(CHECKPOINT_DIR); +} + +auto main(int argc, char* argv[]) -> int { + GlobalInitialize(argc, argv); + + try { + using M = Minkowski; + + const std::vector res { 64 }; + const boundaries_t init_extent { { static_cast(0.0), + static_cast(10.0) } }; + const boundaries_t fldsbc { { FldsBC::PERIODIC, FldsBC::PERIODIC } }; + const boundaries_t prtlbc { { PrtlBC::PERIODIC, PrtlBC::PERIODIC } }; + const std::vector decomp { -1 }; + + const std::vector species_params { ParticleSpecies { + static_cast(1), + "e-", + 1.0f, + -1.0f, + static_cast(10), + timestep_t { 0 }, + timestep_t { 0 }, + ParticlePusher::BORIS, + false, + RadiativeDrag::NONE, + EmissionType::NONE, + static_cast(0), + static_cast(0) } }; + +#if !defined(MPI_ENABLED) + const unsigned int ndomains { 1 }; + adios2::ADIOS adios; + // non-MPI: single domain, sizes always match (equal-sizes path) + const std::vector saved_ncells_per_dom { 64 }; +#else + int mpi_size, mpi_rank; + MPI_Comm_size(MPI_COMM_WORLD, &mpi_size); + MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); + const unsigned int ndomains { static_cast(mpi_size) }; + adios2::ADIOS adios { MPI_COMM_WORLD }; + raise::ErrorIf(mpi_size != 2, "this test requires exactly 2 MPI ranks", HERE); + // asymmetric split: 40+24=64; exercises redecomposeFromCheckpoint + const std::vector saved_ncells_per_dom { 40, 24 }; + const real_t x_split = init_extent[0].second * + static_cast(saved_ncells_per_dom[0]) / + static_cast(64); +#endif + + const auto checkpoint_step = timestep_t { 2 }; + const int prtl_i1_0 { 3 }; + const real_t prtl_dx1_0 { static_cast(0.25) }; + const int prtl_i1_1 { 7 }; + const real_t prtl_dx1_1 { static_cast(0.75) }; + const npart_t npart_placed { 2 }; + + // ── Write phase ─────────────────────────────────────────────────────────── + { + Metadomain md { + ndomains, decomp, res, init_extent, fldsbc, prtlbc, {}, species_params + }; + +#if defined(MPI_ENABLED) + // Simulate load balancing: [32, 32] → [40, 24] + md.redecomposeFromCheckpoint( + { { saved_ncells_per_dom[0] }, { saved_ncells_per_dom[1] } }, + { { { init_extent[0].first, x_split } }, + { { x_split, init_extent[0].second } } }); +#endif + + auto* local = md.subdomain_ptr(md.l_subdomain_indices()[0]); + const auto n_all = local->mesh.n_all()[0]; + + // set em::ex1(i) = real_t(i) over all cells including ghost cells + { + auto em_h = Kokkos::create_mirror_view(local->fields.em); + for (std::size_t i { 0 }; i < n_all; ++i) { + em_h(i, em::ex1) = static_cast(i); + } + Kokkos::deep_copy(local->fields.em, em_h); + } + + // place 2 particles at known positions (safe within the 24-cell domain) + auto& sp = local->species[0]; + sp.set_npart(npart_placed); + { + auto i1_h = Kokkos::create_mirror_view(sp.i1); + auto dx1_h = Kokkos::create_mirror_view(sp.dx1); + i1_h(0) = prtl_i1_0; + dx1_h(0) = static_cast(prtl_dx1_0); + i1_h(1) = prtl_i1_1; + dx1_h(1) = static_cast(prtl_dx1_1); + Kokkos::deep_copy(sp.i1, i1_h); + Kokkos::deep_copy(sp.dx1, dx1_h); + } + + SimulationParams params; + params.set("checkpoint.write_path", std::string { CHECKPOINT_DIR }); + params.set("checkpoint.interval", timestep_t { 1 }); + params.set("checkpoint.interval_time", simtime_t { -1.0 }); + params.set("checkpoint.keep", 5); + params.set("checkpoint.walltime", std::string { "" }); + params.setRawData(toml::value { toml::table {} }); + + md.InitCheckpointWriter(&adios, params); + + const auto wrote = md.WriteCheckpoint( + params, + checkpoint_step, + checkpoint_step, + simtime_t { 0.0 }, + simtime_t { 0.0 }); + raise::ErrorIf(not wrote, "checkpoint was not written", HERE); + } + + // ── Read phase ──────────────────────────────────────────────────────────── + { + // fresh metadomain with the original (uniform) decomposition + // MPI: [32, 32] — mismatch with saved [40, 24] → reconstruction triggered + Metadomain md2 { + ndomains, decomp, res, init_extent, fldsbc, prtlbc, {}, species_params + }; + + SimulationParams params2; + params2.set("checkpoint.read_path", std::string { CHECKPOINT_DIR }); + params2.set("checkpoint.start_step", checkpoint_step); + + md2.ContinueFromCheckpoint(&adios, params2); + + auto* local = md2.subdomain_ptr(md2.l_subdomain_indices()[0]); + const auto lidx = md2.l_subdomain_indices()[0]; + const auto exp_ncells = saved_ncells_per_dom[lidx]; + const auto exp_n_all = exp_ncells + 2 * N_GHOSTS; + + // 1. field array size reflects saved (possibly reconstructed) ncells + raise::ErrorIf( + local->fields.em.extent(0) != exp_n_all, + fmt::format("field em extent: got %lu, expected %lu", + local->fields.em.extent(0), + exp_n_all), + HERE); + + // 2. field values survive the round-trip + { + auto em_h = Kokkos::create_mirror_view(local->fields.em); + Kokkos::deep_copy(em_h, local->fields.em); + for (std::size_t i { 0 }; i < exp_n_all; ++i) { + raise::ErrorIf( + not cmp::AlmostEqual(em_h(i, em::ex1), static_cast(i)), + fmt::format("em::ex1 mismatch at i=%lu: got %f, expected %f", + i, + static_cast(em_h(i, em::ex1)), + static_cast(i)), + HERE); + } + } + + // 3. particle count and coordinates survive the round-trip + auto& sp = local->species[0]; + raise::ErrorIf( + sp.npart() != npart_placed, + fmt::format("particle count: got %lu, expected %lu", + sp.npart(), + npart_placed), + HERE); + { + auto i1_h = Kokkos::create_mirror_view(sp.i1); + auto dx1_h = Kokkos::create_mirror_view(sp.dx1); + Kokkos::deep_copy(i1_h, sp.i1); + Kokkos::deep_copy(dx1_h, sp.dx1); + + raise::ErrorIf( + i1_h(0) != prtl_i1_0, + fmt::format("particle 0 i1: got %d, expected %d", i1_h(0), prtl_i1_0), + HERE); + raise::ErrorIf( + not cmp::AlmostEqual(static_cast(dx1_h(0)), prtl_dx1_0), + "particle 0 dx1 mismatch", + HERE); + raise::ErrorIf( + i1_h(1) != prtl_i1_1, + fmt::format("particle 1 i1: got %d, expected %d", i1_h(1), prtl_i1_1), + HERE); + raise::ErrorIf( + not cmp::AlmostEqual(static_cast(dx1_h(1)), prtl_dx1_1), + "particle 1 dx1 mismatch", + HERE); + } + } + + } catch (const std::exception& e) { + std::cerr << e.what() << '\n'; + cleanup(); + GlobalFinalize(); + return 1; + } + + cleanup(); + GlobalFinalize(); + return 0; +} diff --git a/tests/framework/checkpoint-extents.cpp b/tests/framework/checkpoint-extents.cpp new file mode 100644 index 00000000..ecf6e579 --- /dev/null +++ b/tests/framework/checkpoint-extents.cpp @@ -0,0 +1,163 @@ +/** + * Tests that per-subdomain extents and global mesh extents are correctly saved + * to and restored from a checkpoint, as needed by the moving window feature. + */ + +#include "enums.h" +#include "global.h" + +#include "utils/comparators.h" +#include "utils/error.h" + +#include "framework/domain/metadomain.h" +#include "framework/parameters/parameters.h" + +#include "metrics/minkowski.h" + +#include +#include +#include + +#if defined(MPI_ENABLED) + #include +#endif + +#include +#include +#include +#include + +using namespace ntt; +using namespace metric; + +const std::string CHECKPOINT_DIR = "test_checkpoint_extents"; + +void cleanup() { + std::filesystem::remove_all(CHECKPOINT_DIR); +} + +auto extents_match(const boundaries_t& a, + const boundaries_t& b) -> bool { + if (a.size() != b.size()) { + return false; + } + for (std::size_t d { 0 }; d < a.size(); ++d) { + if (not cmp::AlmostEqual(a[d].first, b[d].first) or + not cmp::AlmostEqual(a[d].second, b[d].second)) { + return false; + } + } + return true; +} + +auto main(int argc, char* argv[]) -> int { + GlobalInitialize(argc, argv); + + try { + using M = Minkowski; + + const std::vector res { 64 }; + const boundaries_t init_extent { { static_cast(0.0), + static_cast(10.0) } }; + const boundaries_t fldsbc { + { FldsBC::PERIODIC, FldsBC::PERIODIC } + }; + const boundaries_t prtlbc { + { PrtlBC::PERIODIC, PrtlBC::PERIODIC } + }; + const std::vector decomp { -1 }; + +#if !defined(MPI_ENABLED) + const unsigned int ndomains { 1 }; + adios2::ADIOS adios; +#else + int mpi_size; + MPI_Comm_size(MPI_COMM_WORLD, &mpi_size); + const unsigned int ndomains { static_cast(mpi_size) }; + adios2::ADIOS adios { MPI_COMM_WORLD }; +#endif + + const auto checkpoint_step = timestep_t { 2 }; + const int n_shift = 4; + + // reference extents captured after the window shift + boundaries_t expected_global_extent; + std::vector> expected_subdomain_extents(ndomains); + + { + Metadomain md { + ndomains, decomp, res, init_extent, fldsbc, prtlbc, {}, {} + }; + + SimulationParams params; + params.set("checkpoint.write_path", std::string { CHECKPOINT_DIR }); + params.set("checkpoint.interval", timestep_t { 1 }); + params.set("checkpoint.interval_time", simtime_t { -1.0 }); + params.set("checkpoint.keep", 5); + params.set("checkpoint.walltime", std::string { "" }); + params.setRawData(toml::value { toml::table {} }); + + md.InitCheckpointWriter(&adios, params); + + // shift the window and capture the resulting extents + md.ShiftByCells(n_shift, in::x1); + expected_global_extent = md.mesh().extent(); + for (unsigned int idx { 0 }; idx < ndomains; ++idx) { + expected_subdomain_extents[idx] = md.subdomain(idx).mesh.extent(); + } + + // write checkpoint at step 2 (finished_step must be > 1 to be saved) + const auto wrote = md.WriteCheckpoint( + params, + checkpoint_step, + checkpoint_step, + simtime_t { 0.0 }, + simtime_t { 0.0 }); + raise::ErrorIf(not wrote, "checkpoint was not written", HERE); + } + + { + // construct a fresh metadomain from the original (pre-shift) extent + Metadomain md2 { + ndomains, decomp, res, init_extent, fldsbc, prtlbc, {}, {} + }; + + // sanity: verify global extent is the original one before reading + raise::ErrorIf( + extents_match(md2.mesh().extent(), expected_global_extent), + "global extent should differ from shifted extent before checkpoint read", + HERE); + + SimulationParams params2; + params2.set("checkpoint.read_path", std::string { CHECKPOINT_DIR }); + params2.set("checkpoint.start_step", checkpoint_step); + + md2.ContinueFromCheckpoint(&adios, params2); + + // global mesh extent must match the shifted one + raise::ErrorIf( + not extents_match(md2.mesh().extent(), expected_global_extent), + "global mesh extent mismatch after checkpoint read", + HERE); + + // per-subdomain extents must also match + for (unsigned int idx { 0 }; idx < ndomains; ++idx) { + raise::ErrorIf( + not extents_match(md2.subdomain(idx).mesh.extent(), + expected_subdomain_extents[idx]), + "subdomain extent mismatch after checkpoint read", + HERE); + } + } + + } catch (const std::exception& e) { + std::cerr << e.what() << '\n'; + cleanup(); + GlobalFinalize(); + return 1; + } + + cleanup(); + GlobalFinalize(); + return 0; +}