diff --git a/CITATION b/CITATION index f685dd7b6..8b1d4947d 100644 --- a/CITATION +++ b/CITATION @@ -41,3 +41,23 @@ archivePrefix = {arXiv}, adsnote = {Provided by the SAO/NASA Astrophysics Data System} } ``` + +For the higher-order shape functions and generalized field stencils, please cite the following paper: + +```latex +@ARTICLE{EntityHighOrder_2026, + author = {{B{\"o}ss}, Ludwig M. and {Vanthieghem}, Arno and {Hakobyan}, Hayk and {Gorbunov}, Evgeny A. and {Caprioli}, Damiano}, + title = "{Entity -- Hardware-agnostic Particle-in-Cell Code for Plasma Astrophysics. III: Higher-order shape functions \& generalized field stencils}", + journal = {arXiv e-prints}, + keywords = {High Energy Astrophysical Phenomena}, + year = 2026, + month = may, + eid = {arXiv:2605.15260}, + pages = {arXiv:2605.15260}, +archivePrefix = {arXiv}, + eprint = {2605.15260}, + primaryClass = {astro-ph.HE}, + adsurl = {https://ui.adsabs.harvard.edu/abs/2026arXiv260515260B}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} +``` diff --git a/README.md b/README.md index 10d1634dd..41ee4c7bc 100644 --- a/README.md +++ b/README.md @@ -12,12 +12,16 @@ Our [detailed documentation](https://entity-toolkit.github.io/) includes everyth ## News +- [May 2026]: our **method paper** for higher-order shape functions and generalized field stencils is [online](https://ui.adsabs.harvard.edu/abs/2026arXiv260515260B/abstract) +- [Apr 2026]: user-defined **custom particle update** functionality [PR #198](https://github.com/entity-toolkit/entity/pull/198) +- [Apr 2026]: **moving window** [PR #196](https://github.com/entity-toolkit/entity/pull/196) +- [Apr 2026]: **built-in piston** [PR #192](https://github.com/entity-toolkit/entity/pull/192) - [Mar 2026]: **single-particle emission** [PR #174](https://github.com/entity-toolkit/entity/pull/174) and [PR #188](https://github.com/entity-toolkit/entity/pull/188) - [Mar 2026]: **spatial sorting of particles** [PR #181](https://github.com/entity-toolkit/entity/pull/181) - [Mar 2026]: **external EM & force fields** [PR #183](https://github.com/entity-toolkit/entity/pull/183) - [Mar 2026]: **examples** of most commonly used custom patterns (see `pgens/examples`) -- [Dec 2025]: **high-order** shape functions [PR #109](https://github.com/entity-toolkit/entity/pull/109) and advanced field stencils [PR #103](https://github.com/entity-toolkit/entity/pull/103) are now supported. -- [Dec 2025]: **particle tracking** is now fully supported via [PR #144](https://github.com/entity-toolkit/entity/pull/144). +- [Dec 2025]: **high-order** shape functions [PR #109](https://github.com/entity-toolkit/entity/pull/109) and advanced field stencils [PR #103](https://github.com/entity-toolkit/entity/pull/103) are now supported +- [Dec 2025]: **particle tracking** is now fully supported via [PR #144](https://github.com/entity-toolkit/entity/pull/144) - [Nov 2025]: our **method papers** are online: [Special relativistic module](https://ui.adsabs.harvard.edu/abs/2025arXiv251117710H/abstract), [GR module](https://ui.adsabs.harvard.edu/abs/2025arXiv251117701G/abstract)! ## Citation diff --git a/input.example.toml b/input.example.toml index dbd52c119..741d5c2d2 100644 --- a/input.example.toml +++ b/input.example.toml @@ -514,10 +514,6 @@ # @type: array # @default: [] custom = "" - # Smoothing window for the output of moments ("Rho", "Charge", "T", ...) - # @type: ushort - # @default: 0 - mom_smooth = "" # Number of timesteps between field outputs # @type: uint # @default: 0 @@ -537,6 +533,19 @@ # @note: If a scalar is given, it is applied to all directions downsampling = "" + [output.fields.smoothing] + # Smoothing order for the output of moments ("Rho", "Charge", "T", ...) + # @type: ushort + # @default: 0 + order = "" + # Smoothing algorithm + # @type: string + # @enum: "const", "spline" + # @default: "spline" + # @note: When using "spline", `order` corresponds to the order of the polynomial used + # @note: When using "const", the smoothing window is `ceil(order / 2)` in both directions + method = "" + [output.particles] # Toggle for the particles output # @type: bool diff --git a/pgens/accretion/pgen.hpp b/pgens/accretion/pgen.hpp index a9dbd65bc..bc792a502 100644 --- a/pgens/accretion/pgen.hpp +++ b/pgens/accretion/pgen.hpp @@ -11,9 +11,9 @@ #include "archetypes/energy_dist.h" #include "archetypes/particle_injector.h" +#include "archetypes/utils.h" #include "framework/domain/metadomain.h" #include "framework/parameters/parameters.h" -#include "kernels/particle_moments.hpp" namespace user { using namespace ntt; @@ -117,39 +117,18 @@ namespace user { std::copy(xi_min.begin(), xi_min.end(), x_min); std::copy(xi_max.begin(), xi_max.end(), x_max); - std::vector specs {}; + std::vector specs {}; for (auto& sp : domain_ptr->species) { if (sp.mass() > 0) { specs.push_back(sp.index()); } } - Kokkos::deep_copy(density, ZERO); - auto scatter_buff = Kokkos::Experimental::create_scatter_view(density); - // some parameters - auto& mesh = domain_ptr->mesh; - const auto use_weights = params.template get( - "particles.use_weights"); - const auto ni2 = mesh.n_active(in::x2); - - for (const auto& sp : specs) { - auto& prtl_spec = domain_ptr->species[sp - 1]; - // clang-format off - Kokkos::parallel_for( - "ComputeMoments", - 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(), - use_weights, - metric, mesh.flds_bc(), - ni2, inv_n0, ZERO)); - // clang-format on - } - Kokkos::Experimental::contribute(density, scatter_buff); + arch::ComputeMomentWithSpecies( + params, + *domain_ptr, + specs, + density); } Inline auto sigma_crit(const coord_t& x_Ph) const -> bool { diff --git a/pgens/reconnection/pgen.hpp b/pgens/reconnection/pgen.hpp index b0450b600..16579ce07 100644 --- a/pgens/reconnection/pgen.hpp +++ b/pgens/reconnection/pgen.hpp @@ -13,7 +13,6 @@ #include "archetypes/spatial_dist.h" #include "archetypes/utils.h" #include "framework/domain/metadomain.h" -#include "kernels/particle_moments.hpp" namespace user { using namespace ntt; @@ -247,35 +246,11 @@ namespace user { inj_box_down.push_back(Range::All); } - { - // compute density of species #1 and #2 - const auto use_weights = params.template get( - "particles.use_weights"); - 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]; - // clang-format off - Kokkos::parallel_for( - "ComputeMoments", - 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(), - use_weights, - domain.mesh.metric, domain.mesh.flds_bc(), - ni2, inv_n0, 0u)); - // clang-format on - } - Kokkos::Experimental::contribute(domain.fields.buff, scatter_buff); - } + // compute density of species #1 and #2 + arch::ComputeMomentWithSpecies(params, + domain, + { 1, 2 }, + domain.fields.buff); const auto replenish_sdist = arch::spatial_dist::ReplenishUniform( domain.mesh.metric, diff --git a/src/archetypes/moving_window.h b/src/archetypes/moving_window.h index 42fd8f9ca..770d9c8c6 100644 --- a/src/archetypes/moving_window.h +++ b/src/archetypes/moving_window.h @@ -28,10 +28,10 @@ namespace arch { const BCTags tags; const int window_shift; - FieldShift_kernel(ndfield_t Fld, - ndfield_t backup_Fld, - const int window_shift, - BCTags tags) + FieldShift_kernel(ndfield_t& Fld, + const ndfield_t& backup_Fld, + int window_shift, + BCTags tags) : Fld { Fld } , backup_Fld { backup_Fld } , window_shift { window_shift } diff --git a/src/archetypes/utils.h b/src/archetypes/utils.h index 51f083aab..1ae775309 100644 --- a/src/archetypes/utils.h +++ b/src/archetypes/utils.h @@ -124,7 +124,8 @@ namespace arch { * @param components Vector of field components to compute (e.g. {} for N, {0} * {1} {2} for V, {0, 1} for T, etc., default: empty, i.e. scalar) * @param buffer_idx Index of the field component in the buffer to save the result to - * @param window Window size for smoothing (in number of cells, default: 0, i.e. no smoothing) + * @param smoothing_order Smoothing order + * @param smoothing_method Smoothing algorithm (e.g. SPLINE, CONST, etc.) * * @tparam S Simulation engine type * @tparam M Metric type @@ -132,13 +133,15 @@ namespace arch { * @tparam N Last dimension of the buffer (e.g. 3 or 6) */ template - inline void ComputeMomentWithSpecies(const SimulationParams& params, - Domain& domain, - const std::vector& species, - ndfield_t& buffer, - const std::vector& components = {}, - idx_t buffer_idx = 0u, - unsigned short window = 0u) { + inline void ComputeMomentWithSpecies( + const SimulationParams& params, + Domain& domain, + const std::vector& species, + ndfield_t& buffer, + const std::vector& components = {}, + idx_t buffer_idx = 0u, + uint8_t smoothing_order = 0u, + OutputSmoothingTypeFlag smoothing_method = OutputSmoothingType::SPLINE) { const auto ni2 = domain.mesh.n_active(in::x2); const auto inv_n0 = ONE / params.template get("scales.n0"); const auto use_weights = params.template get("particles.use_weights"); @@ -153,26 +156,14 @@ namespace arch { kernel::ParticleMoments_kernel(components, scatter_buff, buffer_idx, - 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(), + prtl_spec, use_weights, domain.mesh.metric, domain.mesh.flds_bc(), ni2, inv_n0, - window)); + smoothing_order, + smoothing_method)); } Kokkos::Experimental::contribute(buffer, scatter_buff); } diff --git a/src/engines/grpic/fields_bcs.h b/src/engines/grpic/fields_bcs.h index 3297a17df..06a708dc7 100644 --- a/src/engines/grpic/fields_bcs.h +++ b/src/engines/grpic/fields_bcs.h @@ -39,13 +39,13 @@ namespace ntt { }; template PG> - void MatchFieldsIn(dir::direction_t direction, - Domain& domain, - const Grid& global_grid, - const PG& pgen, - const SimulationParams& params, - BCTags tags, - const gr_bc& g) { + void MatchFieldsIn(const dir::direction_t& direction, + Domain& domain, + const Grid& global_grid, + const PG& pgen, + const SimulationParams& params, + BCTags tags, + const gr_bc& g) { /** * match boundaries */ @@ -138,11 +138,11 @@ namespace ntt { } template - void HorizonFieldsIn(dir::direction_t direction, - Domain& domain, - const SimulationParams& params, - BCTags tags, - const gr_bc& g) { + void HorizonFieldsIn(const dir::direction_t& direction, + Domain& domain, + const SimulationParams& params, + BCTags tags, + const gr_bc& g) { /** * open boundaries */ @@ -176,9 +176,9 @@ namespace ntt { } template - void AxisFieldsIn(dir::direction_t direction, - Domain& domain, - BCTags tags) { + void AxisFieldsIn(const dir::direction_t& direction, + Domain& domain, + BCTags tags) { /** * axis boundaries */ @@ -220,10 +220,10 @@ namespace ntt { } template - void CustomFieldsIn(dir::direction_t direction, - Domain& domain, - BCTags tags, - const gr_bc& g) { + void CustomFieldsIn(const dir::direction_t& direction, + Domain& domain, + BCTags tags, + const gr_bc& g) { (void)direction; (void)domain; (void)tags; diff --git a/src/engines/srpic/particles_bcs.h b/src/engines/srpic/particles_bcs.h index d7b3b9a47..5e79f3b78 100644 --- a/src/engines/srpic/particles_bcs.h +++ b/src/engines/srpic/particles_bcs.h @@ -88,21 +88,20 @@ namespace ntt { if (prtl_spec.npart() == 0) { continue; } - // clang-format off Kokkos::parallel_for( "ComputeMoments", prtl_spec.rangeActiveParticles(), kernel::ParticleMoments_kernel( - {}, scatter_bckp, 0, - 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(), + {}, + scatter_bckp, + 0, + prtl_spec, use_weights, - domain.mesh.metric, domain.mesh.flds_bc(), - ni2, inv_n0, 0)); - // clang-format on + domain.mesh.metric, + domain.mesh.flds_bc(), + ni2, + inv_n0, + 0u)); prtl_spec.set_unsorted(); } Kokkos::Experimental::contribute(domain.fields.bckp, scatter_bckp); diff --git a/src/framework/domain/metadomain_io.cpp b/src/framework/domain/metadomain_io.cpp index 04cdcf1f6..3fc0f5f57 100644 --- a/src/framework/domain/metadomain_io.cpp +++ b/src/framework/domain/metadomain_io.cpp @@ -1,3 +1,4 @@ +#include "defaults.h" #include "enums.h" #include "global.h" @@ -146,25 +147,27 @@ namespace ntt { const auto use_weights = params.template get("particles.use_weights"); const auto ni2 = mesh.n_active(in::x2); const auto inv_n0 = ONE / params.template get("scales.n0"); - const auto window = params.template get( - "output.fields.mom_smooth"); + const auto smooth_order = params.template get( + "output.fields.smoothing.order"); + const auto smooth_method = OutputSmoothingType::from_string( + params.template get("output.fields.smoothing.method")); for (const auto& sp : specs) { auto& prtl_spec = prtl_species[sp - 1]; - // clang-format off Kokkos::parallel_for( "ComputeMoments", prtl_spec.rangeActiveParticles(), - kernel::ParticleMoments_kernel(components, scatter_buff, buff_idx, - 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(), + kernel::ParticleMoments_kernel(components, + scatter_buff, + buff_idx, + prtl_spec, use_weights, - mesh.metric, mesh.flds_bc(), - ni2, inv_n0, window)); - // clang-format on + mesh.metric, + mesh.flds_bc(), + ni2, + inv_n0, + smooth_order, + smooth_method)); } Kokkos::Experimental::contribute(buffer, scatter_buff); } diff --git a/src/framework/domain/metadomain_stats.cpp b/src/framework/domain/metadomain_stats.cpp index 1df1b3510..35ff26e75 100644 --- a/src/framework/domain/metadomain_stats.cpp +++ b/src/framework/domain/metadomain_stats.cpp @@ -19,6 +19,7 @@ #include #include +#include #include #include #include @@ -103,15 +104,10 @@ namespace ntt { Kokkos::parallel_reduce( "ComputeMoments", prtl_spec.rangeActiveParticles(), - // clang-format off kernel::ReducedParticleMoments_kernel(components, - 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(), - use_weights, mesh.metric), - // clang-format on + prtl_spec, + use_weights, + mesh.metric), temp_buff); buffer += temp_buff; } diff --git a/src/framework/parameters/output.cpp b/src/framework/parameters/output.cpp index a926ed585..6c1858e9d 100644 --- a/src/framework/parameters/output.cpp +++ b/src/framework/parameters/output.cpp @@ -20,9 +20,9 @@ namespace ntt { void Output::read(Dimension dim, std::size_t nspec, const toml::value& toml_data) { format = toml::find_or(toml_data, "output", "format", defaults::output::format); global_interval = toml::find_or(toml_data, - "output", - "interval", - defaults::output::interval); + "output", + "interval", + defaults::output::interval); global_interval_time = toml::find_or(toml_data, "output", "interval_time", @@ -34,12 +34,12 @@ namespace ntt { categories.emplace(); for (const auto& category : { "fields", "particles", "spectra", "stats" }) { - const auto q_int = toml::find_or(toml_data, - "output", - category, - "interval", - 0); - const auto q_int_time = toml::find_or(toml_data, + const auto q_int = toml::find_or(toml_data, + "output", + category, + "interval", + 0); + const auto q_int_time = toml::find_or(toml_data, "output", category, "interval_time", @@ -60,10 +60,10 @@ namespace ntt { /* Fields --------------------------------------------------------------- */ const auto flds_out = toml::find_or(toml_data, - "output", - "fields", - "quantities", - std::vector {}); + "output", + "fields", + "quantities", + std::vector {}); const auto custom_flds_out = toml::find_or(toml_data, "output", "fields", @@ -74,11 +74,18 @@ namespace ntt { } fields_quantities = flds_out; fields_custom_quantities = custom_flds_out; - fields_mom_smooth = toml::find_or(toml_data, - "output", - "fields", - "mom_smooth", - defaults::output::mom_smooth); + fields_smooth_order = toml::find_or(toml_data, + "output", + "fields", + "smoothing", + "order", + defaults::output::smooth_order); + fields_smooth_method = toml::find_or(toml_data, + "output", + "fields", + "smoothing", + "method", + defaults::output::smooth_method); fields_downsampling.emplace(); try { auto field_dwn_ = toml::find>(toml_data, @@ -125,35 +132,35 @@ namespace ntt { "species", all_specs); particles_stride = toml::find_or(toml_data, - "output", - "particles", - "stride", - defaults::output::prtl_stride); + "output", + "particles", + "stride", + defaults::output::prtl_stride); /* Spectra -------------------------------------------------------------- */ spectra_e_min = toml::find_or(toml_data, - "output", - "spectra", - "e_min", - defaults::output::spec_emin); + "output", + "spectra", + "e_min", + defaults::output::spec_emin); spectra_e_max = toml::find_or(toml_data, - "output", - "spectra", - "e_max", - defaults::output::spec_emax); + "output", + "spectra", + "e_max", + defaults::output::spec_emax); spectra_log_bins = toml::find_or(toml_data, "output", "spectra", "log_bins", defaults::output::spec_log); spectra_n_bins = toml::find_or(toml_data, - "output", - "spectra", - "n_bins", - defaults::output::spec_nbins); + "output", + "spectra", + "n_bins", + defaults::output::spec_nbins); /* Stats ---------------------------------------------------------------- */ - stats_quantities = toml::find_or(toml_data, + stats_quantities = toml::find_or(toml_data, "output", "stats", "quantities", @@ -190,7 +197,8 @@ namespace ntt { params->set("output.fields.quantities", fields_quantities.value()); params->set("output.fields.custom", fields_custom_quantities.value()); - params->set("output.fields.mom_smooth", fields_mom_smooth.value()); + params->set("output.fields.smoothing.order", fields_smooth_order.value()); + params->set("output.fields.smoothing.method", fields_smooth_method.value()); params->set("output.fields.downsampling", fields_downsampling.value()); params->set("output.particles.species", particles_species.value()); diff --git a/src/framework/parameters/output.h b/src/framework/parameters/output.h index 94afff324..27b077bc5 100644 --- a/src/framework/parameters/output.h +++ b/src/framework/parameters/output.h @@ -43,16 +43,17 @@ namespace ntt { std::optional> fields_quantities; std::optional> fields_custom_quantities; - std::optional fields_mom_smooth; std::optional> fields_downsampling; + std::optional fields_smooth_order; + std::optional fields_smooth_method; std::optional> particles_species; std::optional particles_stride; - std::optional spectra_e_min; - std::optional spectra_e_max; - std::optional spectra_log_bins; - std::optional spectra_n_bins; + std::optional spectra_e_min; + std::optional spectra_e_max; + std::optional spectra_log_bins; + std::optional spectra_n_bins; std::optional> stats_quantities; std::optional> stats_custom_quantities; diff --git a/src/global/defaults.h b/src/global/defaults.h index 5558a6c9f..e1387677e 100644 --- a/src/global/defaults.h +++ b/src/global/defaults.h @@ -67,7 +67,8 @@ namespace ntt::defaults { namespace output { const std::string format = "BPFile"; const timestep_t interval = 100; - const unsigned short mom_smooth = 0; + const unsigned short smooth_order = 0; + const std::string smooth_method = "spline"; const npart_t prtl_stride = 100; const real_t spec_emin = 1e-3; const real_t spec_emax = 1e3; diff --git a/src/global/enums.h b/src/global/enums.h index 68eefb75b..5700a3a28 100644 --- a/src/global/enums.h +++ b/src/global/enums.h @@ -18,6 +18,7 @@ * - enum ntt::ParticlePusher // photon, boris, vay, gca, none * - enum ntt::RadiativeDrag // compton, synchrotron, none * - enum ntt::EmissionType // none, synchrotron, inversecompton, custom + * - enum ntt::OutputSmoothingTypeFlag // none, const, spline * * @namespaces: * - ntt:: @@ -401,6 +402,40 @@ namespace ntt { using EmissionTypeFlag = uint8_t; + namespace OutputSmoothingType { + enum OutputSmoothingTypeFlag_ : uint8_t { + NONE = 0, + CONST = 1, + SPLINE = 2, + }; + + inline auto to_string(uint8_t flags) -> std::string { + switch (flags) { + case NONE: + return "none"; + case CONST: + return "const"; + case SPLINE: + return "spline"; + default: + return "unknown"; + } + } + + inline auto from_string(const std::string& str) -> uint8_t { + const auto str_lo = fmt::toLower(str); + if (str_lo == "const") { + return CONST; + } else if (str_lo == "spline") { + return SPLINE; + } else { + return NONE; + } + } + } // namespace OutputSmoothingType + + using OutputSmoothingTypeFlag = uint8_t; + } // namespace ntt #endif // GLOBAL_ENUMS_H diff --git a/src/kernels/fields_bcs.hpp b/src/kernels/fields_bcs.hpp index 73c023484..729a09d10 100644 --- a/src/kernels/fields_bcs.hpp +++ b/src/kernels/fields_bcs.hpp @@ -827,7 +827,7 @@ namespace kernel::bc { const ncells_t i_edge; const bool setE { false }, setB { false }; - AxisBoundaries_kernel(ndfield_t Fld, ncells_t i_edge, BCTags tags) + AxisBoundaries_kernel(ndfield_t& Fld, ncells_t i_edge, BCTags tags) : Fld { Fld } , i_edge { i_edge } , setE { tags & BC::Ex1 or tags & BC::Ex2 or tags & BC::Ex3 } @@ -1191,10 +1191,10 @@ namespace kernel::bc { const bool setE { false }, setB { false }; const ncells_t nfilter; - HorizonBoundaries_kernel(ndfield_t Fld, - ncells_t i1_min, - BCTags tags, - ncells_t nfilter) + HorizonBoundaries_kernel(ndfield_t& Fld, + ncells_t i1_min, + BCTags tags, + ncells_t nfilter) : Fld { Fld } , i1_min { i1_min } , setE { (tags & BC::Ex1 or tags & BC::Ex2 or tags & BC::Ex3) or @@ -1249,10 +1249,10 @@ namespace kernel::bc { const real_t xg_edge; const real_t dx_abs; - AbsorbCurrents_kernel(ndfield_t J, - const M& metric, - real_t xg_edge, - real_t dx_abs) + AbsorbCurrents_kernel(ndfield_t& J, + const M& metric, + real_t xg_edge, + real_t dx_abs) : J { J } , metric { metric } , xg_edge { xg_edge } diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp index b2919c9f1..1244ba21d 100644 --- a/src/kernels/particle_moments.hpp +++ b/src/kernels/particle_moments.hpp @@ -23,6 +23,7 @@ #include "utils/numeric.h" #include "framework/containers/particles.h" +#include "kernels/particle_shapes.hpp" #include #include @@ -49,75 +50,52 @@ namespace kernel { (F == FldsID::Nppc) || (F == FldsID::T) || (F == FldsID::V), "Invalid field ID"); - const uint8_t c1, c2; - scatter_ndfield_t Buff; - const idx_t buff_idx; - const array_t i1, i2, i3; - const array_t dx1, dx2, dx3; - const array_t ux1, ux2, ux3; - const array_t phi; - const array_t weight; - const array_t tag; - const float mass; - const float charge; - const bool use_weights; - const M metric; - const int ni2; - const uint8_t window; + const uint8_t c1, c2; + scatter_ndfield_t Buff; + const idx_t buff_idx; + const ParticleArrays particles; + const float mass; + const float charge; + const bool use_weights; + const M metric; + const int ni2; + const real_t inv_n0; + + const uint8_t order; + const uint8_t window; + const OutputSmoothingTypeFlag smoothing; const real_t contrib; - const real_t smooth; bool is_axis_i2min { false }, is_axis_i2max { false }; public: - ParticleMoments_kernel(const std::vector& components, - const scatter_ndfield_t& scatter_buff, - idx_t buff_idx, - const array_t& i1, - const array_t& i2, - const array_t& i3, - const array_t& dx1, - const array_t& dx2, - const array_t& dx3, - const array_t& ux1, - const array_t& ux2, - const array_t& ux3, - const array_t& phi, - const array_t& weight, - const array_t& tag, - float mass, - float charge, - bool use_weights, - const M& metric, - const boundaries_t& boundaries, - ncells_t ni2, - real_t inv_n0, - uint8_t window) + ParticleMoments_kernel( + const std::vector& components, + const scatter_ndfield_t& scatter_buff, + idx_t buff_idx, + const Particles& particles, + bool use_weights, + const M& metric, + const boundaries_t& boundaries, + ncells_t ni2, + real_t inv_n0, + uint8_t order = 0u, + OutputSmoothingTypeFlag smoothing = OutputSmoothingType::SPLINE) : c1 { not components.empty() ? components[0] : static_cast(0) } , c2 { (components.size() == 2) ? components[1] : static_cast(0) } , Buff { scatter_buff } , buff_idx { buff_idx } - , i1 { i1 } - , i2 { i2 } - , i3 { i3 } - , dx1 { dx1 } - , dx2 { dx2 } - , dx3 { dx3 } - , ux1 { ux1 } - , ux2 { ux2 } - , ux3 { ux3 } - , phi { phi } - , weight { weight } - , tag { tag } - , mass { mass } - , charge { charge } + , particles { static_cast(particles) } + , mass { particles.mass() } + , charge { particles.charge() } , use_weights { use_weights } , metric { metric } , ni2 { static_cast(ni2) } - , window { window } - , contrib { get_contrib(mass, charge) } - , smooth { inv_n0 / (real_t)(math::pow(TWO * (real_t)window + ONE, - static_cast(D))) } { + , inv_n0 { inv_n0 } + , order { order } + , smoothing { smoothing } + , window { static_cast(math::ceil(static_cast(order) / 2.0f)) } + , contrib { get_contrib(mass, charge) } { raise::ErrorIf(buff_idx >= N, "Invalid buffer index", HERE); raise::ErrorIf(window > N_GHOSTS, "Window size too large", HERE); raise::ErrorIf(((F == FldsID::Rho) || (F == FldsID::Charge)) && (mass == ZERO), @@ -131,28 +109,70 @@ namespace kernel { } } + Inline auto shapeFunction(real_t delta_x) const -> real_t { + if (smoothing == OutputSmoothingType::SPLINE) { + if (order == 0) { + return ONE; + } else if (order == 1) { + return prtl_shape::S1(delta_x); + } else if (order == 2) { + return prtl_shape::S2(delta_x); + } else if (order == 3) { + return prtl_shape::S3(delta_x); + } else if (order == 4) { + return prtl_shape::S4(delta_x); + } else if (order == 5) { + return prtl_shape::S5(delta_x); + } else if (order == 6) { + return prtl_shape::S6(delta_x); + } else if (order == 7) { + return prtl_shape::S7(delta_x); + } else if (order == 8) { + return prtl_shape::S8(delta_x); + } else if (order == 9) { + return prtl_shape::S9(delta_x); + } else if (order == 10) { + return prtl_shape::S10(delta_x); + } else if (order == 11) { + return prtl_shape::S11(delta_x); + } else { + raise::KernelError(HERE, "Unsupported shape function order"); + return ZERO; + } + } else if (smoothing == OutputSmoothingType::CONST) { + return ONE / (TWO * static_cast(window) + ONE); + } else { + raise::KernelError(HERE, "Unsupported smoothing method"); + return ZERO; + } + } + Inline auto computeStressEnergyComponent(prtlidx_t p) const -> real_t { real_t u0 { ZERO }; vec_t u_Phys { ZERO }; if constexpr (S == SimEngine::SRPIC) { // stress-energy tensor for SR is computed in the tetrad (hatted) basis if constexpr (M::CoordType == Coord::Cartesian) { - u_Phys[0] = ux1(p); - u_Phys[1] = ux2(p); - u_Phys[2] = ux3(p); + u_Phys[0] = particles.ux1(p); + u_Phys[1] = particles.ux2(p); + u_Phys[2] = particles.ux3(p); } else { static_assert(D != Dim::_1D, "non-Cartesian SRPIC 1D"); coord_t x_Code { ZERO }; - x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); - x_Code[1] = static_cast(i2(p)) + static_cast(dx2(p)); + x_Code[0] = static_cast(particles.i1(p)) + + static_cast(particles.dx1(p)); + x_Code[1] = static_cast(particles.i2(p)) + + static_cast(particles.dx2(p)); if constexpr (D == Dim::_3D) { - x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); + x_Code[2] = static_cast(particles.i3(p)) + + static_cast(particles.dx3(p)); } else { - x_Code[2] = phi(p); + x_Code[2] = particles.phi(p); } - metric.template transform_xyz(x_Code, - { ux1(p), ux2(p), ux3(p) }, - u_Phys); + metric.template transform_xyz( + x_Code, + { particles.ux1(p), particles.ux2(p), particles.ux3(p) }, + u_Phys); } u0 = (mass == ZERO) ? (NORM(u_Phys[0], u_Phys[1], u_Phys[2])) @@ -161,20 +181,24 @@ namespace kernel { // stress-energy tensor for GR is computed in contravariant basis static_assert(D != Dim::_1D, "GRPIC 1D"); coord_t x_Code { ZERO }; - x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); - x_Code[1] = static_cast(i2(p)) + static_cast(dx2(p)); + x_Code[0] = static_cast(particles.i1(p)) + + static_cast(particles.dx1(p)); + x_Code[1] = static_cast(particles.i2(p)) + + static_cast(particles.dx2(p)); if constexpr (D == Dim::_3D) { - x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); + x_Code[2] = static_cast(particles.i3(p)) + + static_cast(particles.dx3(p)); } // raise full covariant 4-vector to get correct contravariant u^i // u^i != h^{ij} u_j - const real_t u_0_cov { metric.u_0(x_Code, - { ux1(p), ux2(p), ux3(p) }, - (mass == ZERO) ? ZERO : ONE) }; + const real_t u_0_cov { metric.u_0( + x_Code, + { particles.ux1(p), particles.ux2(p), particles.ux3(p) }, + (mass == ZERO) ? ZERO : ONE) }; vec_t u_cntrv_4d { ZERO }; metric.template transform_4d( x_Code, - { u_0_cov, ux1(p), ux2(p), ux3(p) }, + { u_0_cov, particles.ux1(p), particles.ux2(p), particles.ux3(p) }, u_cntrv_4d); // in GR: u^0 = Gamma/alpha u0 = u_cntrv_4d[0]; @@ -203,21 +227,25 @@ namespace kernel { // for bulk 3vel (tetrad basis) vec_t u_Phys { ZERO }; if constexpr (M::CoordType == Coord::Cartesian) { - u_Phys[0] = ux1(p); - u_Phys[1] = ux2(p); - u_Phys[2] = ux3(p); + u_Phys[0] = particles.ux1(p); + u_Phys[1] = particles.ux2(p); + u_Phys[2] = particles.ux3(p); } else { coord_t x_Code { ZERO }; - x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); - x_Code[1] = static_cast(i2(p)) + static_cast(dx2(p)); + x_Code[0] = static_cast(particles.i1(p)) + + static_cast(particles.dx1(p)); + x_Code[1] = static_cast(particles.i2(p)) + + static_cast(particles.dx2(p)); if constexpr (D == Dim::_3D) { - x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); + x_Code[2] = static_cast(particles.i3(p)) + + static_cast(particles.dx3(p)); } else { - x_Code[2] = phi(p); + x_Code[2] = particles.phi(p); } - metric.template transform_xyz(x_Code, - { ux1(p), ux2(p), ux3(p) }, - u_Phys); + metric.template transform_xyz( + x_Code, + { particles.ux1(p), particles.ux2(p), particles.ux3(p) }, + u_Phys); } if (mass == ZERO) { u0 = NORM(u_Phys[0], u_Phys[1], u_Phys[2]); @@ -231,20 +259,24 @@ namespace kernel { // GR: Eckart frame flux N^μ = m * u^μ / u^0 static_assert(D != Dim::_1D, "GRPIC 1D"); coord_t x_Code { ZERO }; - x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); - x_Code[1] = static_cast(i2(p)) + static_cast(dx2(p)); + x_Code[0] = static_cast(particles.i1(p)) + + static_cast(particles.dx1(p)); + x_Code[1] = static_cast(particles.i2(p)) + + static_cast(particles.dx2(p)); if constexpr (D == Dim::_3D) { - x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); + x_Code[2] = static_cast(particles.i3(p)) + + static_cast(particles.dx3(p)); } // raise full covariant 4-vector to get correct contravariant u^μ // u^i != h^{ij} u_j - const real_t u_0_cov { metric.u_0(x_Code, - { ux1(p), ux2(p), ux3(p) }, - (mass == ZERO) ? ZERO : ONE) }; + const real_t u_0_cov { metric.u_0( + x_Code, + { particles.ux1(p), particles.ux2(p), particles.ux3(p) }, + (mass == ZERO) ? ZERO : ONE) }; vec_t u_cntrv_4d { ZERO }; metric.template transform_4d( x_Code, - { u_0_cov, ux1(p), ux2(p), ux3(p) }, + { u_0_cov, particles.ux1(p), particles.ux2(p), particles.ux3(p) }, u_cntrv_4d); const real_t u0 { u_cntrv_4d[0] }; // Deposit flux N^μ = mass * u^μ / u^0 @@ -256,7 +288,7 @@ namespace kernel { } Inline void operator()(prtlidx_t p) const { - if (tag(p) == ParticleTag::dead) { + if (particles.tag(p) == ParticleTag::dead) { return; } real_t coeff { ZERO }; @@ -276,48 +308,57 @@ namespace kernel { // for nppc calculation ... // ... do not take volume, weights or smoothing into account if constexpr (D == Dim::_1D) { - coeff *= smooth / - metric.sqrt_det_h({ static_cast(i1(p)) + HALF }); + coeff *= inv_n0 / metric.sqrt_det_h( + { static_cast(particles.i1(p)) + HALF }); } else if constexpr (D == Dim::_2D) { - coeff *= smooth / - metric.sqrt_det_h({ static_cast(i1(p)) + HALF, - static_cast(i2(p)) + HALF }); + coeff *= inv_n0 / metric.sqrt_det_h( + { static_cast(particles.i1(p)) + HALF, + static_cast(particles.i2(p)) + HALF }); } else if constexpr (D == Dim::_3D) { - coeff *= smooth / - metric.sqrt_det_h({ static_cast(i1(p)) + HALF, - static_cast(i2(p)) + HALF, - static_cast(i3(p)) + HALF }); + coeff *= inv_n0 / metric.sqrt_det_h( + { static_cast(particles.i1(p)) + HALF, + static_cast(particles.i2(p)) + HALF, + static_cast(particles.i3(p)) + HALF }); } if (use_weights) { - coeff *= weight(p); + coeff *= particles.weight(p); } } auto buff_access = Buff.access(); if constexpr (D == Dim::_1D) { for (auto di1 { -window }; di1 <= window; ++di1) { - buff_access(i1(p) + di1 + N_GHOSTS, buff_idx) += coeff; + const real_t delta_x1 = math::abs(static_cast(particles.dx1(p)) - + (static_cast(di1) + HALF)); + buff_access(particles.i1(p) + di1 + N_GHOSTS, + buff_idx) += coeff * shapeFunction(delta_x1); } } else if constexpr (D == Dim::_2D) { for (auto di2 { -window }; di2 <= window; ++di2) { for (auto di1 { -window }; di1 <= window; ++di1) { + const real_t delta_x1 = math::abs(static_cast(particles.dx1(p)) - + (static_cast(di1) + HALF)); + const real_t delta_x2 = math::abs(static_cast(particles.dx2(p)) - + (static_cast(di2) + HALF)); + const auto shape_coeff = shapeFunction(delta_x1) * + shapeFunction(delta_x2); if constexpr (M::CoordType == Coord::Cartesian) { - buff_access(i1(p) + di1 + N_GHOSTS, - i2(p) + di2 + N_GHOSTS, - buff_idx) += coeff; + buff_access(particles.i1(p) + di1 + N_GHOSTS, + particles.i2(p) + di2 + N_GHOSTS, + buff_idx) += coeff * shape_coeff; } else { // reflect contribution at axes - if (is_axis_i2min && (i2(p) + di2 < 0)) { - buff_access(i1(p) + di1 + N_GHOSTS, - N_GHOSTS - (i2(p) + di2), - buff_idx) += coeff; - } else if (is_axis_i2max && (i2(p) + di2 >= ni2)) { - buff_access(i1(p) + di1 + N_GHOSTS, - 2 * ni2 - (i2(p) + di2) + N_GHOSTS, - buff_idx) += coeff; + if (is_axis_i2min && (particles.i2(p) + di2 < 0)) { + buff_access(particles.i1(p) + di1 + N_GHOSTS, + N_GHOSTS - (particles.i2(p) + di2), + buff_idx) += coeff * shape_coeff; + } else if (is_axis_i2max && (particles.i2(p) + di2 >= ni2)) { + buff_access(particles.i1(p) + di1 + N_GHOSTS, + 2 * ni2 - (particles.i2(p) + di2) + N_GHOSTS, + buff_idx) += coeff * shape_coeff; } else { - buff_access(i1(p) + di1 + N_GHOSTS, - i2(p) + di2 + N_GHOSTS, - buff_idx) += coeff; + buff_access(particles.i1(p) + di1 + N_GHOSTS, + particles.i2(p) + di2 + N_GHOSTS, + buff_idx) += coeff * shape_coeff; } } } @@ -326,28 +367,37 @@ namespace kernel { for (auto di3 { -window }; di3 <= window; ++di3) { for (auto di2 { -window }; di2 <= window; ++di2) { for (auto di1 { -window }; di1 <= window; ++di1) { + const auto delta_x1 = math::abs(static_cast(particles.dx1(p)) - + (static_cast(di1) + HALF)); + const auto delta_x2 = math::abs(static_cast(particles.dx2(p)) - + (static_cast(di2) + HALF)); + const auto delta_x3 = math::abs(static_cast(particles.dx3(p)) - + (static_cast(di3) + HALF)); + const auto shape_coeff = shapeFunction(delta_x1) * + shapeFunction(delta_x2) * + shapeFunction(delta_x3); if constexpr (M::CoordType == Coord::Cartesian) { - buff_access(i1(p) + di1 + N_GHOSTS, - i2(p) + di2 + N_GHOSTS, - i3(p) + di3 + N_GHOSTS, - buff_idx) += coeff; + buff_access(particles.i1(p) + di1 + N_GHOSTS, + particles.i2(p) + di2 + N_GHOSTS, + particles.i3(p) + di3 + N_GHOSTS, + buff_idx) += coeff * shape_coeff; } else { // reflect contribution at axes - if (is_axis_i2min && (i2(p) + di2 < 0)) { - buff_access(i1(p) + di1 + N_GHOSTS, - N_GHOSTS - (i2(p) + di2), - i3(p) + di3 + N_GHOSTS, - buff_idx) += coeff; - } else if (is_axis_i2max && (i2(p) + di2 >= ni2)) { - buff_access(i1(p) + di1 + N_GHOSTS, - 2 * ni2 - (i2(p) + di2) + N_GHOSTS, - i3(p) + di3 + N_GHOSTS, - buff_idx) += coeff; + if (is_axis_i2min && (particles.i2(p) + di2 < 0)) { + buff_access(particles.i1(p) + di1 + N_GHOSTS, + N_GHOSTS - (particles.i2(p) + di2), + particles.i3(p) + di3 + N_GHOSTS, + buff_idx) += coeff * shape_coeff; + } else if (is_axis_i2max && (particles.i2(p) + di2 >= ni2)) { + buff_access(particles.i1(p) + di1 + N_GHOSTS, + 2 * ni2 - (particles.i2(p) + di2) + N_GHOSTS, + particles.i3(p) + di3 + N_GHOSTS, + buff_idx) += coeff * shape_coeff; } else { - buff_access(i1(p) + di1 + N_GHOSTS, - i2(p) + di2 + N_GHOSTS, - i3(p) + di3 + N_GHOSTS, - buff_idx) += coeff; + buff_access(particles.i1(p) + di1 + N_GHOSTS, + particles.i2(p) + di2 + N_GHOSTS, + particles.i3(p) + di3 + N_GHOSTS, + buff_idx) += coeff * shape_coeff; } } } @@ -680,7 +730,7 @@ namespace kernel { bool log_bins, size_t n_bins, const M& metric) - : particles { static_cast(particles) } + : particles { static_cast(particles) } , is_massive { (particles.mass() != 0.0f) } , dn_scatter { dn_scatter } , e_min { e_min } diff --git a/src/kernels/particle_shapes.hpp b/src/kernels/particle_shapes.hpp index 13e345625..8e155b795 100644 --- a/src/kernels/particle_shapes.hpp +++ b/src/kernels/particle_shapes.hpp @@ -17,6 +17,34 @@ namespace prtl_shape { + // clang-format off + // S(x) = 1 - |x| |x| < 1 + // 0.0 |x| ≥ 1 + // + // clang-format on + Inline real_t S1(const real_t x) { + if (x < ONE) { + return ONE - x; + } else { + return ZERO; + } + } + + // clang-format off + // 3/4 - |x|^2 |x| < 1/2 + // S(x) = 1/2 * (3/2 - |x|)^2 1/2 ≤ |x| < 3/2 + // 0.0 |x| ≥ 3/2 + // clang-format on + Inline real_t S2(const real_t x) { + if (x < HALF) { + return static_cast(0.75) - SQR(x); + } else if (x < THREE_HALFS) { + return HALF * SQR(THREE_HALFS - x); + } else { + return ZERO; + } + } + // clang-format off // 115/192 - (5/8) * |x|^2 + (1/4) * |x|^4 |x| < 1/2 // S(x) = 55/96 + (5/24) * |x| - (5/4) * |x|^2 + (5/6) * |x|^3 - (1/6) * |x|^4 1/2 ≤ |x| < 3/2 diff --git a/src/kernels/reduced_stats.hpp b/src/kernels/reduced_stats.hpp index e15924b93..a6411f681 100644 --- a/src/kernels/reduced_stats.hpp +++ b/src/kernels/reduced_stats.hpp @@ -18,6 +18,8 @@ #include "traits/metric.h" #include "utils/numeric.h" +#include "framework/containers/particles.h" + namespace kernel { using namespace ntt; @@ -401,54 +403,25 @@ namespace kernel { class ReducedParticleMoments_kernel { static constexpr auto D = M::Dim; - const uint8_t c1, c2; - const array_t i1, i2, i3; - const array_t dx1, dx2, dx3; - const array_t ux1, ux2, ux3; - const array_t phi; - const array_t weight; - const array_t tag; - const float mass; - const float charge; - const bool use_weights; - const M metric; + const uint8_t c1, c2; + const ParticleArrays particles; + const float mass; + const float charge; + const bool use_weights; + const M metric; const real_t contrib; public: ReducedParticleMoments_kernel(const std::vector& components, - const array_t& i1, - const array_t& i2, - const array_t& i3, - const array_t& dx1, - const array_t& dx2, - const array_t& dx3, - const array_t& ux1, - const array_t& ux2, - const array_t& ux3, - const array_t& phi, - const array_t& weight, - const array_t& tag, - float mass, - float charge, - bool use_weights, - const M& metric) + const Particles& particles, + bool use_weights, + const M& metric) : c1 { not components.empty() ? components[0] : static_cast(0) } , c2 { (components.size() == 2) ? components[1] : static_cast(0) } - , i1 { i1 } - , i2 { i2 } - , i3 { i3 } - , dx1 { dx1 } - , dx2 { dx2 } - , dx3 { dx3 } - , ux1 { ux1 } - , ux2 { ux2 } - , ux3 { ux3 } - , phi { phi } - , weight { weight } - , tag { tag } - , mass { mass } - , charge { charge } + , particles { static_cast(particles) } + , mass { particles.mass() } + , charge { particles.charge() } , use_weights { use_weights } , metric { metric } , contrib { get_contrib

(mass, charge) } { @@ -459,7 +432,7 @@ namespace kernel { } Inline void operator()(prtlidx_t p, real_t& buff) const { - if (tag(p) != ParticleTag::alive) { + if (particles.tag(p) != ParticleTag::alive) { return; } auto dV = ONE; @@ -470,17 +443,20 @@ namespace kernel { } else { coord_t x_Code { ZERO }; if constexpr ((D == Dim::_1D) or (D == Dim::_2D) or (D == Dim::_3D)) { - x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); + x_Code[0] = static_cast(particles.i1(p)) + + static_cast(particles.dx1(p)); } if constexpr ((D == Dim::_2D) or (D == Dim::_3D)) { - x_Code[1] = static_cast(i2(p)) + static_cast(dx2(p)); + x_Code[1] = static_cast(particles.i2(p)) + + static_cast(particles.dx2(p)); } if constexpr (D == Dim::_3D) { - x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); + x_Code[2] = static_cast(particles.i3(p)) + + static_cast(particles.dx3(p)); } dV = metric.sqrt_det_h(x_Code); if constexpr (P == StatsID::N or P == StatsID::Rho or P == StatsID::Charge) { - buff += dV * (use_weights ? weight(p) : contrib); + buff += dV * (use_weights ? particles.weight(p) : contrib); } else { // for stress-energy tensor real_t energy { ZERO }; @@ -489,23 +465,25 @@ namespace kernel { // SR // stress-energy tensor for SR is computed in the tetrad (hatted) basis if constexpr (M::CoordType == Coord::Cartesian) { - u_Phys[0] = ux1(p); - u_Phys[1] = ux2(p); - u_Phys[2] = ux3(p); + u_Phys[0] = particles.ux1(p); + u_Phys[1] = particles.ux2(p); + u_Phys[2] = particles.ux3(p); } else { static_assert(D != Dim::_1D, "non-Cartesian SRPIC 1D"); coord_t x_Code { ZERO }; - x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); - x_Code[1] = static_cast(i2(p)) + static_cast(dx2(p)); + x_Code[0] = static_cast(particles.i1(p)) + + static_cast(particles.dx1(p)); + x_Code[1] = static_cast(particles.i2(p)) + + static_cast(particles.dx2(p)); if constexpr (D == Dim::_3D) { - x_Code[2] = static_cast(i3(p)) + - static_cast(dx3(p)); + x_Code[2] = static_cast(particles.i3(p)) + + static_cast(particles.dx3(p)); } else { - x_Code[2] = phi(p); + x_Code[2] = particles.phi(p); } metric.template transform_xyz( x_Code, - { ux1(p), ux2(p), ux3(p) }, + { particles.ux1(p), particles.ux2(p), particles.ux3(p) }, u_Phys); } if (mass == ZERO) { @@ -519,18 +497,22 @@ namespace kernel { // stress-energy tensor for GR is computed in contravariant basis static_assert(D != Dim::_1D, "GRPIC 1D"); coord_t x_Code { ZERO }; - x_Code[0] = static_cast(i1(p)) + static_cast(dx1(p)); - x_Code[1] = static_cast(i2(p)) + static_cast(dx2(p)); + x_Code[0] = static_cast(particles.i1(p)) + + static_cast(particles.dx1(p)); + x_Code[1] = static_cast(particles.i2(p)) + + static_cast(particles.dx2(p)); if constexpr (D == Dim::_3D) { - x_Code[2] = static_cast(i3(p)) + static_cast(dx3(p)); + x_Code[2] = static_cast(particles.i3(p)) + + static_cast(particles.dx3(p)); } vec_t u_Cntrv { ZERO }; // compute u_i u^i for energy - metric.template transform(x_Code, - { ux1(p), ux2(p), ux3(p) }, - u_Cntrv); - energy = u_Cntrv[0] * ux1(p) + u_Cntrv[1] * ux2(p) + - u_Cntrv[2] * ux3(p); + metric.template transform( + x_Code, + { particles.ux1(p), particles.ux2(p), particles.ux3(p) }, + u_Cntrv); + energy = u_Cntrv[0] * particles.ux1(p) + + u_Cntrv[1] * particles.ux2(p) + u_Cntrv[2] * particles.ux3(p); if (mass == ZERO) { energy = math::sqrt(energy); } else { diff --git a/src/output/utils/interpret_prompt.cpp b/src/output/utils/interpret_prompt.cpp index ee85299a5..27c027315 100644 --- a/src/output/utils/interpret_prompt.cpp +++ b/src/output/utils/interpret_prompt.cpp @@ -5,6 +5,7 @@ #include "utils/error.h" #include "utils/formatting.h" +#include #include #include diff --git a/tests/framework/parameters.cpp b/tests/framework/parameters.cpp index 82821b899..edb054709 100644 --- a/tests/framework/parameters.cpp +++ b/tests/framework/parameters.cpp @@ -92,10 +92,12 @@ const auto mink_1d = u8R"( [output.fields] quantities = ["Rho", "J", "B"] - mom_smooth = 2 downsampling = [4, 5] interval = 100 + [output.fields.smoothing] + order = 2 + [output.particles] species = [1, 2] stride = 100 diff --git a/tests/kernels/particle_moments.cpp b/tests/kernels/particle_moments.cpp index 23d417748..cbda7546f 100644 --- a/tests/kernels/particle_moments.cpp +++ b/tests/kernels/particle_moments.cpp @@ -9,9 +9,8 @@ #include "utils/numeric.h" #include "metrics/minkowski.h" -#include "metrics/qspherical.h" -#include "metrics/spherical.h" +#include "framework/containers/particles.h" #include "kernels/reduced_stats.hpp" #include @@ -67,42 +66,44 @@ void testParticleMoments(const std::vector& res, const auto nx2 = res[1]; ndfield_t buff { "buff", nx1 + 2 * N_GHOSTS, nx2 + 2 * N_GHOSTS }; - array_t i1 { "i1", 10 }; - array_t i2 { "i2", 10 }; - array_t i3 { "i3", 0 }; - array_t dx1 { "dx1", 10 }; - array_t dx2 { "dx2", 10 }; - array_t dx3 { "dx3", 0 }; - array_t ux1 { "ux1", 10 }; - array_t ux2 { "ux2", 10 }; - array_t ux3 { "ux3", 10 }; - array_t phi { "phi", 10 }; - array_t weight { "weight", 10 }; - array_t tag { "tag", 10 }; - const float mass = 1.0; - const float charge = 1.0; - const bool use_weights = false; - const real_t inv_n0 = 1.0; - - put_value(i1, 5, 0); - put_value(i2, 4, 0); - put_value(dx1, (prtldx_t)(0.15), 0); - put_value(dx2, (prtldx_t)(0.85), 0); - put_value(ux1, (real_t)(1.0), 0); - put_value(ux2, (real_t)(-2.0), 0); - put_value(ux3, (real_t)(3.0), 0); - put_value(tag, ParticleTag::alive, 0); - put_value(weight, 1.0, 0); - - put_value(i1, 2, 4); - put_value(i2, 2, 4); - put_value(dx1, (prtldx_t)(0.22), 4); - put_value(dx2, (prtldx_t)(0.55), 4); - put_value(ux1, (real_t)(-3.0), 4); - put_value(ux2, (real_t)(2.0), 4); - put_value(ux3, (real_t)(-1.0), 4); - put_value(tag, ParticleTag::alive, 4); - put_value(weight, 1.0, 4); + Particles particles { 1u, + "test", + 1.0f, + 1.0f, + 10, + 0, + 0, + ParticlePusher::NONE, + false, + RadiativeDrag::NONE, + EmissionType::NONE, + 0, + 0 }; + + const float mass = 1.0; + const float charge = 1.0; + const bool use_weights = false; + const real_t inv_n0 = 1.0; + + put_value(particles.i1, 5, 0); + put_value(particles.i2, 4, 0); + put_value(particles.dx1, (prtldx_t)(0.15), 0); + put_value(particles.dx2, (prtldx_t)(0.85), 0); + put_value(particles.ux1, (real_t)(1.0), 0); + put_value(particles.ux2, (real_t)(-2.0), 0); + put_value(particles.ux3, (real_t)(3.0), 0); + put_value(particles.tag, ParticleTag::alive, 0); + put_value(particles.weight, 1.0, 0); + + put_value(particles.i1, 2, 4); + put_value(particles.i2, 2, 4); + put_value(particles.dx1, (prtldx_t)(0.22), 4); + put_value(particles.dx2, (prtldx_t)(0.55), 4); + put_value(particles.ux1, (real_t)(-3.0), 4); + put_value(particles.ux2, (real_t)(2.0), 4); + put_value(particles.ux3, (real_t)(-1.0), 4); + put_value(particles.tag, ParticleTag::alive, 4); + put_value(particles.weight, 1.0, 4); auto boundaries = boundaries_t {}; if constexpr (M::CoordType != Coord::Cartesian) { @@ -118,91 +119,87 @@ void testParticleMoments(const std::vector& res, const uint8_t window = 1; auto scatter_buff = Kokkos::Experimental::create_scatter_view(buff); - // clang-format off Kokkos::parallel_for( - "ParticleMoments", 10, - kernel::ParticleMoments_kernel(comp1, scatter_buff, 0, - i1, i2, i3, - dx1, dx2, dx3, - ux1, ux2, ux3, - phi, weight, tag, - mass, charge, + "ParticleMoments", + 10, + kernel::ParticleMoments_kernel(comp1, + scatter_buff, + 0, + particles, use_weights, metric, - boundaries, nx2, inv_n0, window)); + boundaries, + nx2, + inv_n0, + window)); Kokkos::parallel_for( - "ParticleMoments", 10, - kernel::ParticleMoments_kernel(comp2, scatter_buff, 1, - i1, i2, i3, - dx1, dx2, dx3, - ux1, ux2, ux3, - phi, weight, tag, - mass, charge, + "ParticleMoments", + 10, + kernel::ParticleMoments_kernel(comp2, + scatter_buff, + 1, + particles, use_weights, metric, - boundaries, nx2, inv_n0, window)); + boundaries, + nx2, + inv_n0, + window)); Kokkos::parallel_for( - "ParticleMoments", 10, - kernel::ParticleMoments_kernel(comp3, scatter_buff, 2, - i1, i2, i3, - dx1, dx2, dx3, - ux1, ux2, ux3, - phi, weight, tag, - mass, charge, + "ParticleMoments", + 10, + kernel::ParticleMoments_kernel(comp3, + scatter_buff, + 2, + particles, use_weights, metric, - boundaries, nx2, inv_n0, window)); + boundaries, + nx2, + inv_n0, + window)); real_t n = ZERO, npart = ZERO, rho = ZERO, t00 = ZERO; Kokkos::parallel_reduce( - "ReducedParticleMoments", 10, - kernel::ReducedParticleMoments_kernel({}, - i1, i2, i3, - dx1, dx2, dx3, - ux1, ux2, ux3, - phi, weight, tag, - mass, charge, + "ReducedParticleMoments", + 10, + kernel::ReducedParticleMoments_kernel({}, + particles, use_weights, - metric), n); + metric), + n); Kokkos::parallel_reduce( - "ReducedParticleMoments", 10, - kernel::ReducedParticleMoments_kernel({}, - i1, i2, i3, - dx1, dx2, dx3, - ux1, ux2, ux3, - phi, weight, tag, - mass, charge, + "ReducedParticleMoments", + 10, + kernel::ReducedParticleMoments_kernel({}, + particles, use_weights, - metric), npart); + metric), + npart); Kokkos::parallel_reduce( - "ReducedParticleMoments", 10, - kernel::ReducedParticleMoments_kernel({}, - i1, i2, i3, - dx1, dx2, dx3, - ux1, ux2, ux3, - phi, weight, tag, - mass, charge, + "ReducedParticleMoments", + 10, + kernel::ReducedParticleMoments_kernel({}, + particles, use_weights, - metric), rho); + metric), + rho); Kokkos::parallel_reduce( - "ReducedParticleMoments", 10, - kernel::ReducedParticleMoments_kernel({0u, 0u}, - i1, i2, i3, - dx1, dx2, dx3, - ux1, ux2, ux3, - phi, weight, tag, - mass, charge, + "ReducedParticleMoments", + 10, + kernel::ReducedParticleMoments_kernel({ 0u, 0u }, + particles, use_weights, - metric), t00); - // clang-format on + metric), + t00); Kokkos::Experimental::contribute(buff, scatter_buff); - auto i1_h = Kokkos::create_mirror_view(i1); - auto i2_h = Kokkos::create_mirror_view(i2); - Kokkos::deep_copy(i1_h, i1); - Kokkos::deep_copy(i2_h, i2); + auto i1_h = Kokkos::create_mirror_view(particles.i1); + auto i2_h = Kokkos::create_mirror_view(particles.i2); + Kokkos::deep_copy(i1_h, particles.i1); + Kokkos::deep_copy(i2_h, particles.i2); const auto h1 = metric.sqrt_det_h({ static_cast(i1_h(0)) + HALF, static_cast(i2_h(0)) + HALF }); @@ -234,7 +231,7 @@ void testParticleMoments(const std::vector& res, const real_t gammaSQR_2_expect = 15.0; const real_t n_expect = 2.0; - const real_t t00_expect = 2.0 * math::sqrt(15.0); + const real_t t00_expect = static_cast(2.0 * math::sqrt(15.0)); errorIf(not cmp::AlmostEqual_host(gammaSQR_1, gammaSQR_1_expect, epsilon * acc), fmt::format("wrong gamma_1 %.12e %.12e for %dD %s",