From 58ca99f9c57cd85fcb8a97bb094b481dd35522db Mon Sep 17 00:00:00 2001 From: haykh Date: Sat, 15 Mar 2025 14:42:26 -0400 Subject: [PATCH 01/25] deposit test pgens --- setups/tests/deposit/deposit-gr.toml | 53 +++++++++++ setups/tests/deposit/deposit-mink.toml | 69 +++++++++++++++ setups/tests/deposit/deposit-sr.toml | 53 +++++++++++ setups/tests/deposit/pgen.hpp | 116 +++++++++---------------- 4 files changed, 218 insertions(+), 73 deletions(-) create mode 100644 setups/tests/deposit/deposit-gr.toml create mode 100644 setups/tests/deposit/deposit-mink.toml create mode 100644 setups/tests/deposit/deposit-sr.toml diff --git a/setups/tests/deposit/deposit-gr.toml b/setups/tests/deposit/deposit-gr.toml new file mode 100644 index 000000000..077e3e59a --- /dev/null +++ b/setups/tests/deposit/deposit-gr.toml @@ -0,0 +1,53 @@ +[simulation] + name = "deposit-test" + engine = "srpic" + runtime = 10.0 + +[grid] + resolution = [256, 256] + extent = [[0.0, 1.0], [0.0, 1.0]] + + [grid.metric] + metric = "minkowski" + + [grid.boundaries] + fields = [["PERIODIC"], ["PERIODIC"]] + particles = [["PERIODIC"], ["PERIODIC"]] + +[scales] + larmor0 = 0.1 + skindepth0 = 0.1 + +[algorithms] + current_filters = 4 + + [algorithms.timestep] + CFL = 0.5 + +[particles] + ppc0 = 10.0 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 1e2 + + [[particles.species]] + label = "e+" + mass = 1.0 + charge = 1.0 + maxnpart = 1e2 + +[setup] + +[output] + format = "hdf5" + interval_time = 0.01 + + [output.quantities] + quantities = ["N_1", "N_2", "E", "B", "J"] + +[diagnostics] + colored_stdout = true + blocking_timers = true diff --git a/setups/tests/deposit/deposit-mink.toml b/setups/tests/deposit/deposit-mink.toml new file mode 100644 index 000000000..df70f7552 --- /dev/null +++ b/setups/tests/deposit/deposit-mink.toml @@ -0,0 +1,69 @@ +[simulation] + name = "deposit-test-mink" + engine = "srpic" + runtime = 10.0 + +[grid] + resolution = [512, 512] + extent = [[0.0, 1.0], [0.0, 1.0]] + + [grid.metric] + metric = "minkowski" + + [grid.boundaries] + fields = [["PERIODIC"], ["PERIODIC"]] + particles = [["PERIODIC"], ["PERIODIC"]] + +[scales] + larmor0 = 0.1 + skindepth0 = 0.1 + +[algorithms] + current_filters = 4 + + [algorithms.timestep] + CFL = 0.5 + +[particles] + ppc0 = 10.0 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 1e2 + + [[particles.species]] + label = "e+" + mass = 1.0 + charge = 1.0 + maxnpart = 1e2 + +[setup] + x1s = [0.25] + y1s = [0.85] + z1s = [0.33] + ux1s = [0.6] + uy1s = [-0.3] + uz1s = [-0.2] + + x2s = [0.25] + y2s = [0.85] + z2s = [0.33] + ux2s = [-0.2] + uy2s = [-0.2] + uz2s = [0.1] + +[output] + format = "hdf5" + interval_time = 0.01 + + [output.fields] + quantities = ["N_1", "N_2", "E", "B", "J"] + +[checkpoint] + keep = 0 + +[diagnostics] + colored_stdout = true + blocking_timers = true diff --git a/setups/tests/deposit/deposit-sr.toml b/setups/tests/deposit/deposit-sr.toml new file mode 100644 index 000000000..077e3e59a --- /dev/null +++ b/setups/tests/deposit/deposit-sr.toml @@ -0,0 +1,53 @@ +[simulation] + name = "deposit-test" + engine = "srpic" + runtime = 10.0 + +[grid] + resolution = [256, 256] + extent = [[0.0, 1.0], [0.0, 1.0]] + + [grid.metric] + metric = "minkowski" + + [grid.boundaries] + fields = [["PERIODIC"], ["PERIODIC"]] + particles = [["PERIODIC"], ["PERIODIC"]] + +[scales] + larmor0 = 0.1 + skindepth0 = 0.1 + +[algorithms] + current_filters = 4 + + [algorithms.timestep] + CFL = 0.5 + +[particles] + ppc0 = 10.0 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 1e2 + + [[particles.species]] + label = "e+" + mass = 1.0 + charge = 1.0 + maxnpart = 1e2 + +[setup] + +[output] + format = "hdf5" + interval_time = 0.01 + + [output.quantities] + quantities = ["N_1", "N_2", "E", "B", "J"] + +[diagnostics] + colored_stdout = true + blocking_timers = true diff --git a/setups/tests/deposit/pgen.hpp b/setups/tests/deposit/pgen.hpp index fd9a41c2e..1be4deb69 100644 --- a/setups/tests/deposit/pgen.hpp +++ b/setups/tests/deposit/pgen.hpp @@ -4,14 +4,8 @@ #include "enums.h" #include "global.h" -#include "arch/kokkos_aliases.h" #include "arch/traits.h" -#include "utils/comparators.h" -#include "utils/formatting.h" -#include "utils/log.h" -#include "utils/numeric.h" -#include "archetypes/energy_dist.h" #include "archetypes/particle_injector.h" #include "archetypes/problem_generator.h" #include "framework/domain/domain.h" @@ -26,10 +20,20 @@ namespace user { struct PGen : public arch::ProblemGenerator { // compatibility traits for the problem generator - static constexpr auto engines = traits::compatible_with::value; - static constexpr auto metrics = traits::compatible_with::value; - static constexpr auto dimensions = - traits::compatible_with::value; + static constexpr auto engines { + traits::compatible_with::value + }; + static constexpr auto metrics { + traits::compatible_with::value + }; + static constexpr auto dimensions { + traits::compatible_with::value + }; // for easy access to variables in the child class using arch::ProblemGenerator::D; @@ -47,85 +51,51 @@ namespace user { const auto x1s = params.template get>("setup.x1s", empty); const auto y1s = params.template get>("setup.y1s", empty); const auto z1s = params.template get>("setup.z1s", empty); - const auto ux1s = params.template get>("setup.ux1s", + const auto phi1s = params.template get>("setup.phi1s", + empty); + const auto ux1s = params.template get>("setup.ux1s", empty); - const auto uy1s = params.template get>("setup.uy1s", + const auto uy1s = params.template get>("setup.uy1s", empty); - const auto uz1s = params.template get>("setup.uz1s", + const auto uz1s = params.template get>("setup.uz1s", empty); const auto x2s = params.template get>("setup.x2s", empty); const auto y2s = params.template get>("setup.y2s", empty); const auto z2s = params.template get>("setup.z2s", empty); - const auto ux2s = params.template get>("setup.ux2s", + const auto ux2s = params.template get>("setup.ux2s", empty); - const auto uy2s = params.template get>("setup.uy2s", + const auto uy2s = params.template get>("setup.uy2s", empty); - const auto uz2s = params.template get>("setup.uz2s", + const auto uz2s = params.template get>("setup.uz2s", empty); - // std::vector x, y, z, ux_1, uy_1, uz_1, ux_2, uy_2, uz_2; - // x.push_back(0.85); - // x.push_back(0.123); - // if constexpr (D == Dim::_2D || D == Dim::_3D) { - // y.push_back(0.32); - // y.push_back(0.321); - // } - // if constexpr (D == Dim::_3D) { - // z.push_back(0.231); - // z.push_back(0.687); - // } - // ux_1.push_back(1.0); - // uy_1.push_back(-1.0); - // uz_1.push_back(0.0); - // ux_1.push_back(1.0); - // uy_1.push_back(-2.0); - // uz_1.push_back(1.0); - // - // ux_2.push_back(1.0); - // uy_2.push_back(1.0); - // uz_2.push_back(0.0); - // ux_2.push_back(-2.0); - // uy_2.push_back(3.0); - // uz_2.push_back(-1.0); - // - const std::map> data_1 { - { "x1", x1s}, - { "x2", y1s}, - { "x3", z1s}, - {"ux1", ux1s}, - {"ux2", uy1s}, - {"ux3", uz1s} + const auto phi2s = params.template get>("setup.phi2s", + empty); + std::map> data_1 { + { "x1", x1s }, + { "x2", y1s }, + { "ux1", ux1s }, + { "ux2", uy1s }, + { "ux3", uz1s } }; - const std::map> data_2 { - { "x1", x2s}, - { "x2", y2s}, - { "x3", z2s}, - {"ux1", ux2s}, - {"ux2", uy2s}, - {"ux3", uz2s} + std::map> data_2 { + { "x1", x2s }, + { "x2", y2s }, + { "ux1", ux2s }, + { "ux2", uy2s }, + { "ux3", uz2s } }; + if constexpr (M::CoordType == Coord::Cart or D == Dim::_3D) { + data_1["x3"] = z1s; + data_2["x3"] = z2s; + } else if constexpr (D == Dim::_2D) { + data_1["phi"] = phi1s; + data_2["phi"] = phi2s; + } arch::InjectGlobally(global_domain, local_domain, (arch::spidx_t)1, data_1); arch::InjectGlobally(global_domain, local_domain, (arch::spidx_t)2, data_2); } - - // void CustomPostStep(std::size_t, long double time, Domain& domain) { - // if (time >= 0.1) { - // for (auto& species : domain.species) { - // auto ux1 = species.ux1; - // auto ux2 = species.ux2; - // auto ux3 = species.ux3; - // Kokkos::parallel_for( - // "Stop", - // species.rangeActiveParticles(), - // Lambda(index_t p) { - // ux1(p) = ZERO; - // ux2(p) = ZERO; - // ux3(p) = ZERO; - // }); - // } - // } - // } }; } // namespace user From 9fcca0b05f4b6378509cc268131a7bdd20c92bbb Mon Sep 17 00:00:00 2001 From: haykh Date: Sat, 15 Mar 2025 14:43:12 -0400 Subject: [PATCH 02/25] output domain shapes --- src/framework/domain/communications.cpp | 18 +++++----- src/framework/domain/output.cpp | 18 +++++++++- src/output/tests/writer-mpi.cpp | 1 + src/output/tests/writer-nompi.cpp | 1 + src/output/writer.cpp | 48 ++++++++++++++----------- src/output/writer.h | 6 +++- 6 files changed, 60 insertions(+), 32 deletions(-) diff --git a/src/framework/domain/communications.cpp b/src/framework/domain/communications.cpp index 7dc5d285a..0bd06b82e 100644 --- a/src/framework/domain/communications.cpp +++ b/src/framework/domain/communications.cpp @@ -36,10 +36,10 @@ namespace ntt { using comm_params_t = std::pair>; template - auto GetSendRecvRanks(Metadomain* metadomain, - Domain& domain, - dir::direction_t direction) - -> std::pair { + auto GetSendRecvRanks( + Metadomain* metadomain, + Domain& domain, + dir::direction_t direction) -> std::pair { Domain* send_to_nghbr_ptr = nullptr; Domain* recv_from_nghbr_ptr = nullptr; // set pointers to the correct send/recv domains @@ -119,11 +119,11 @@ namespace ntt { } template - auto GetSendRecvParams(Metadomain* metadomain, - Domain& domain, - dir::direction_t direction, - bool synchronize) - -> std::pair { + auto GetSendRecvParams( + Metadomain* metadomain, + Domain& domain, + dir::direction_t direction, + bool synchronize) -> std::pair { const auto [send_indrank, recv_indrank] = GetSendRecvRanks(metadomain, domain, direction); const auto [send_ind, send_rank] = send_indrank; diff --git a/src/framework/domain/output.cpp b/src/framework/domain/output.cpp index 3dc95c18a..c51894e55 100644 --- a/src/framework/domain/output.cpp +++ b/src/framework/domain/output.cpp @@ -70,6 +70,7 @@ namespace ntt { g_writer.defineMeshLayout(glob_shape_with_ghosts, off_ncells_with_ghosts, loc_shape_with_ghosts, + { local_domain->index(), ndomains() }, params.template get>( "output.fields.downsampling"), incl_ghosts, @@ -227,6 +228,17 @@ namespace ntt { const auto dwn = params.template get>( "output.fields.downsampling"); + auto off_ncells_with_ghosts = local_domain->offset_ncells(); + auto loc_shape_with_ghosts = local_domain->mesh.n_active(); + { // compute positions/sizes of meshblocks in cells in all dimensions + const auto off_ndomains = local_domain->offset_ndomains(); + if (incl_ghosts) { + for (auto d { 0 }; d <= M::Dim; ++d) { + off_ncells_with_ghosts[d] += 2 * N_GHOSTS * off_ndomains[d]; + loc_shape_with_ghosts[d] += 2 * N_GHOSTS; + } + } + } for (unsigned short dim = 0; dim < M::Dim; ++dim) { const auto l_size = local_domain->mesh.n_active()[dim]; const auto l_offset = local_domain->offset_ncells()[dim]; @@ -275,7 +287,11 @@ namespace ntt { xe(offset + i_dwn + 1) = x_Ph[dim]; } }); - g_writer.writeMesh(dim, xc, xe); + g_writer.writeMesh( + dim, + xc, + xe, + { off_ncells_with_ghosts[dim], loc_shape_with_ghosts[dim] }); } const auto output_asis = params.template get("output.debug.as_is"); // !TODO: this can probably be optimized to dump things at once diff --git a/src/output/tests/writer-mpi.cpp b/src/output/tests/writer-mpi.cpp index f6d3ee88a..4ec6941ac 100644 --- a/src/output/tests/writer-mpi.cpp +++ b/src/output/tests/writer-mpi.cpp @@ -64,6 +64,7 @@ auto main(int argc, char* argv[]) -> int { writer.defineMeshLayout({ static_cast(mpi_size) * nx1 }, { static_cast(mpi_rank) * nx1 }, { nx1 }, + { mpi_rank, mpi_size }, { dwn1 }, false, Coord::Cart); diff --git a/src/output/tests/writer-nompi.cpp b/src/output/tests/writer-nompi.cpp index 8fb2ac026..7e795c2a9 100644 --- a/src/output/tests/writer-nompi.cpp +++ b/src/output/tests/writer-nompi.cpp @@ -74,6 +74,7 @@ auto main(int argc, char* argv[]) -> int { writer.defineMeshLayout({ nx1, nx2, nx3 }, { 0, 0, 0 }, { nx1, nx2, nx3 }, + { 0, 1 }, { dwn1, dwn2, dwn3 }, false, Coord::Cart); diff --git a/src/output/writer.cpp b/src/output/writer.cpp index 95965c864..e2e094535 100644 --- a/src/output/writer.cpp +++ b/src/output/writer.cpp @@ -63,12 +63,14 @@ namespace out { m_mode = mode; } - void Writer::defineMeshLayout(const std::vector& glob_shape, - const std::vector& loc_corner, - const std::vector& loc_shape, - const std::vector& dwn, - bool incl_ghosts, - Coord coords) { + void Writer::defineMeshLayout( + const std::vector& glob_shape, + const std::vector& loc_corner, + const std::vector& loc_shape, + const std::pair& domain_idx, + const std::vector& dwn, + bool incl_ghosts, + Coord coords) { m_flds_ghosts = incl_ghosts; m_dwn = dwn; @@ -99,24 +101,24 @@ namespace out { for (std::size_t i { 0 }; i < m_flds_g_shape.size(); ++i) { // cell-centers - adios2::Dims g_shape = { m_flds_g_shape_dwn[i] }; - adios2::Dims l_corner = { m_flds_l_corner_dwn[i] }; - adios2::Dims l_shape = { m_flds_l_shape_dwn[i] }; m_io.DefineVariable("X" + std::to_string(i + 1), - g_shape, - l_corner, - l_shape, + { m_flds_g_shape_dwn[i] }, + { m_flds_l_corner_dwn[i] }, + { m_flds_l_shape_dwn[i] }, adios2::ConstantDims); // cell-edges - const auto is_last = (m_flds_l_corner[i] + m_flds_l_shape[i] == + const auto is_last = (m_flds_l_corner[i] + m_flds_l_shape[i] == m_flds_g_shape[i]); - adios2::Dims g_shape1 = { m_flds_g_shape_dwn[i] + 1 }; - adios2::Dims l_shape1 = { m_flds_l_shape_dwn[i] + (is_last ? 1 : 0) }; m_io.DefineVariable("X" + std::to_string(i + 1) + "e", - g_shape1, - l_corner, - l_shape1, + { m_flds_g_shape_dwn[i] + 1 }, + { m_flds_l_corner_dwn[i] }, + { m_flds_l_shape_dwn[i] + (is_last ? 1 : 0) }, adios2::ConstantDims); + m_io.DefineVariable("N" + std::to_string(i + 1) + "l", + { 2 * domain_idx.second }, + { 2 * domain_idx.first }, + { 2 }, + adios2::ConstantDims); } if constexpr (std::is_same::array_layout, @@ -402,9 +404,10 @@ namespace out { m_writer.Put(var, e_bins_h); } - void Writer::writeMesh(unsigned short dim, - const array_t& xc, - const array_t& xe) { + void Writer::writeMesh(unsigned short dim, + const array_t& xc, + const array_t& xe, + const std::vector& loc_off_sz) { auto varc = m_io.InquireVariable("X" + std::to_string(dim + 1)); auto vare = m_io.InquireVariable("X" + std::to_string(dim + 1) + "e"); auto xc_h = Kokkos::create_mirror_view(xc); @@ -413,6 +416,9 @@ namespace out { Kokkos::deep_copy(xe_h, xe); m_writer.Put(varc, xc_h); m_writer.Put(vare, xe_h); + auto vard = m_io.InquireVariable( + "N" + std::to_string(dim + 1) + "l"); + m_writer.Put(vard, loc_off_sz.data()); } void Writer::beginWriting(WriteModeTags write_mode, diff --git a/src/output/writer.h b/src/output/writer.h index a8abf4b12..0b95dfb0e 100644 --- a/src/output/writer.h +++ b/src/output/writer.h @@ -86,6 +86,7 @@ namespace out { void defineMeshLayout(const std::vector&, const std::vector&, const std::vector&, + const std::pair&, const std::vector&, bool, Coord); @@ -94,7 +95,10 @@ namespace out { void defineParticleOutputs(Dimension, const std::vector&); void defineSpectraOutputs(const std::vector&); - void writeMesh(unsigned short, const array_t&, const array_t&); + void writeMesh(unsigned short, + const array_t&, + const array_t&, + const std::vector&); template void writeField(const std::vector&, From 1e18874d3c870198ba1c74ead9ab59c1f0eca51b Mon Sep 17 00:00:00 2001 From: haykh Date: Thu, 27 Mar 2025 17:45:57 -0400 Subject: [PATCH 03/25] i+di coord push for cart --- src/kernels/particle_pusher_sr.hpp | 112 +++++++++++++++++++---------- 1 file changed, 76 insertions(+), 36 deletions(-) diff --git a/src/kernels/particle_pusher_sr.hpp b/src/kernels/particle_pusher_sr.hpp index 2e8a5f652..50fbec5b0 100644 --- a/src/kernels/particle_pusher_sr.hpp +++ b/src/kernels/particle_pusher_sr.hpp @@ -562,45 +562,85 @@ namespace kernel::sr { Inline void posUpd(bool massive, index_t& p, coord_t& xp) const { // get cartesian velocity - const real_t inv_energy { - massive ? ONE / math::sqrt(ONE + SQR(ux1(p)) + SQR(ux2(p)) + SQR(ux3(p))) - : ONE / math::sqrt(SQR(ux1(p)) + SQR(ux2(p)) + SQR(ux3(p))) - }; - vec_t vp_Cart { ux1(p) * inv_energy, - ux2(p) * inv_energy, - ux3(p) * inv_energy }; - // get cartesian position - coord_t xp_Cart { ZERO }; - metric.template convert_xyz(xp, xp_Cart); - // update cartesian position - for (auto d = 0u; d < M::PrtlDim; ++d) { - xp_Cart[d] += vp_Cart[d] * dt; - } - // transform back to code - metric.template convert_xyz(xp_Cart, xp); - - // update x1 - if constexpr (D == Dim::_1D || D == Dim::_2D || D == Dim::_3D) { - i1_prev(p) = i1(p); - dx1_prev(p) = dx1(p); - from_Xi_to_i_di(xp[0], i1(p), dx1(p)); - } + if constexpr (M::CoordType == Coord::Cart) { + // i+di push for Cartesian basis + const real_t dt_inv_energy { + massive + ? (dt / math::sqrt(ONE + SQR(ux1(p)) + SQR(ux2(p)) + SQR(ux3(p)))) + : (dt / math::sqrt(SQR(ux1(p)) + SQR(ux2(p)) + SQR(ux3(p)))) + }; + if constexpr (D == Dim::_1D || D == Dim::_2D || D == Dim::_3D) { + i1_prev(p) = i1(p); + dx1_prev(p) = dx1(p); + dx1(p) += metric.template transform<1, Idx::XYZ, Idx::U>(xp, ux1(p)) * + dt_inv_energy; + i1(p) += static_cast(dx1(p) >= ONE) - + static_cast(dx1(p) < ZERO); + dx1(p) -= (dx1(p) >= ONE); + dx1(p) += (dx1(p) < ZERO); + } + if constexpr (D == Dim::_2D || D == Dim::_3D) { + i2_prev(p) = i2(p); + dx2_prev(p) = dx2(p); + dx2(p) += metric.template transform<2, Idx::XYZ, Idx::U>(xp, ux2(p)) * + dt_inv_energy; + i2(p) += static_cast(dx2(p) >= ONE) - + static_cast(dx2(p) < ZERO); + dx2(p) -= (dx2(p) >= ONE); + dx2(p) += (dx2(p) < ZERO); + } + if constexpr (D == Dim::_3D) { + i3_prev(p) = i3(p); + dx3_prev(p) = dx3(p); + dx3(p) += metric.template transform<3, Idx::XYZ, Idx::U>(xp, ux3(p)) * + dt_inv_energy; + i3(p) += static_cast(dx3(p) >= ONE) - + static_cast(dx3(p) < ZERO); + dx3(p) -= (dx3(p) >= ONE); + dx3(p) += (dx3(p) < ZERO); + } + } else { + // full Cartesian coordinate push in non-Cartesian basis + const real_t inv_energy { + massive ? ONE / math::sqrt(ONE + SQR(ux1(p)) + SQR(ux2(p)) + SQR(ux3(p))) + : ONE / math::sqrt(SQR(ux1(p)) + SQR(ux2(p)) + SQR(ux3(p))) + }; + vec_t vp_Cart { ux1(p) * inv_energy, + ux2(p) * inv_energy, + ux3(p) * inv_energy }; + // get cartesian position + coord_t xp_Cart { ZERO }; + metric.template convert_xyz(xp, xp_Cart); + // update cartesian position + for (auto d = 0u; d < M::PrtlDim; ++d) { + xp_Cart[d] += vp_Cart[d] * dt; + } + // transform back to code + metric.template convert_xyz(xp_Cart, xp); + + // update x1 + if constexpr (D == Dim::_1D || D == Dim::_2D || D == Dim::_3D) { + i1_prev(p) = i1(p); + dx1_prev(p) = dx1(p); + from_Xi_to_i_di(xp[0], i1(p), dx1(p)); + } - // update x2 & phi - if constexpr (D == Dim::_2D || D == Dim::_3D) { - i2_prev(p) = i2(p); - dx2_prev(p) = dx2(p); - from_Xi_to_i_di(xp[1], i2(p), dx2(p)); - if constexpr (D == Dim::_2D && M::PrtlDim == Dim::_3D) { - phi(p) = xp[2]; + // update x2 & phi + if constexpr (D == Dim::_2D || D == Dim::_3D) { + i2_prev(p) = i2(p); + dx2_prev(p) = dx2(p); + from_Xi_to_i_di(xp[1], i2(p), dx2(p)); + if constexpr (D == Dim::_2D && M::PrtlDim == Dim::_3D) { + phi(p) = xp[2]; + } } - } - // update x3 - if constexpr (D == Dim::_3D) { - i3_prev(p) = i3(p); - dx3_prev(p) = dx3(p); - from_Xi_to_i_di(xp[2], i3(p), dx3(p)); + // update x3 + if constexpr (D == Dim::_3D) { + i3_prev(p) = i3(p); + dx3_prev(p) = dx3(p); + from_Xi_to_i_di(xp[2], i3(p), dx3(p)); + } } boundaryConditions(p, xp); } From 0ee06ce516691783e22527ec296231cf56e617f9 Mon Sep 17 00:00:00 2001 From: haykh Date: Wed, 9 Apr 2025 07:43:26 -0400 Subject: [PATCH 04/25] new deposit scheme --- src/framework/parameters.cpp | 16 +- src/kernels/currents_deposit.hpp | 460 +++++++++++++++---------------- 2 files changed, 220 insertions(+), 256 deletions(-) diff --git a/src/framework/parameters.cpp b/src/framework/parameters.cpp index 79f8e333e..e4a40522b 100644 --- a/src/framework/parameters.cpp +++ b/src/framework/parameters.cpp @@ -31,10 +31,10 @@ namespace ntt { template - auto get_dx0_V0(const std::vector& resolution, - const boundaries_t& extent, - const std::map& params) - -> std::pair { + auto get_dx0_V0( + const std::vector& resolution, + const boundaries_t& extent, + const std::map& params) -> std::pair { const auto metric = M(resolution, extent, params); const auto dx0 = metric.dxMin(); coord_t x_corner { ZERO }; @@ -738,13 +738,7 @@ namespace ntt { for (const auto& e : extent_pairwise) { min_extent = std::min(min_extent, e.second - e.first); } - const auto default_ds = min_extent * defaults::bc::match::ds_frac; - std::size_t n_match_bcs { 0u }; - for (const auto& bcs : flds_bc_pairwise) { - if (bcs.first == FldsBC::MATCH or bcs.second == FldsBC::MATCH) { - n_match_bcs += 1; - } - } + const auto default_ds = min_extent * defaults::bc::match::ds_frac; boundaries_t ds_array; try { auto ds = toml::find(toml_data, "grid", "boundaries", "match", "ds"); diff --git a/src/kernels/currents_deposit.hpp b/src/kernels/currents_deposit.hpp index ca9a94878..98d00a9b0 100644 --- a/src/kernels/currents_deposit.hpp +++ b/src/kernels/currents_deposit.hpp @@ -67,8 +67,8 @@ namespace kernel { const array_t& weight, const array_t& tag, const M& metric, - const real_t& charge, - const real_t& dt) + real_t charge, + real_t dt) : J { scatter_cur } , i1 { i1 } , i2 { i2 } @@ -100,45 +100,70 @@ namespace kernel { if (tag(p) == ParticleTag::dead) { return; } - // _f = final, _i = initial - tuple_t Ip_f, Ip_i; - coord_t xp_f, xp_i, xp_r; + // recover particle velocity to deposit in unsimulated direction vec_t vp { ZERO }; + { + coord_t xp { ZERO }; + if constexpr (D == Dim::_1D) { + xp[0] = i_di_to_Xi(i1(p), dx1(p)); + } else if constexpr (D == Dim::_2D) { + if constexpr (M::PrtlDim == Dim::_3D) { + xp[0] = i_di_to_Xi(i1(p), dx1(p)); + xp[1] = i_di_to_Xi(i2(p), dx2(p)); + xp[2] = phi(p); + } else { + xp[0] = i_di_to_Xi(i1(p), dx1(p)); + xp[1] = i_di_to_Xi(i2(p), dx2(p)); + } + } else { + xp[0] = i_di_to_Xi(i1(p), dx1(p)); + xp[1] = i_di_to_Xi(i2(p), dx2(p)); + xp[2] = i_di_to_Xi(i3(p), dx3(p)); + } + auto inv_energy { ZERO }; + if constexpr (S == SimEngine::SRPIC) { + metric.template transform_xyz(xp, + { ux1(p), ux2(p), ux3(p) }, + vp); + inv_energy = ONE / math::sqrt(ONE + NORM_SQR(ux1(p), ux2(p), ux3(p))); + } else { + metric.template transform(xp, + { ux1(p), ux2(p), ux3(p) }, + vp); + inv_energy = ONE / math::sqrt(ONE + ux1(p) * vp[0] + ux2(p) * vp[1] + + ux3(p) * vp[2]); + } + if (Kokkos::isnan(vp[2]) || Kokkos::isinf(vp[2])) { + vp[2] = ZERO; + } + vp[0] *= inv_energy; + vp[1] *= inv_energy; + vp[2] *= inv_energy; + } - // get [i, di]_init and [i, di]_final (per dimension) - getDepositInterval(p, Ip_f, Ip_i, xp_f, xp_i, xp_r); - // recover particle velocity to deposit in unsimulated direction - getPrtl3Vel(p, vp); const real_t coeff { weight(p) * charge }; - depositCurrentsFromParticle(coeff, vp, Ip_f, Ip_i, xp_f, xp_i, xp_r); - } - /** - * @brief Deposit currents from a single particle. - * @param[in] coeff Particle weight x charge. - * @param[in] vp Particle 3-velocity. - * @param[in] Ip_f Final position of the particle (cell index). - * @param[in] Ip_i Initial position of the particle (cell index). - * @param[in] xp_f Final position. - * @param[in] xp_i Previous step position. - * @param[in] xp_r Intermediate point used in zig-zag deposit. - */ - Inline auto depositCurrentsFromParticle(const real_t& coeff, - const vec_t& vp, - const tuple_t& Ip_f, - const tuple_t& Ip_i, - const coord_t& xp_f, - const coord_t& xp_i, - const coord_t& xp_r) const -> void { - const real_t Wx1_1 { HALF * (xp_i[0] + xp_r[0]) - - static_cast(Ip_i[0]) }; - const real_t Wx1_2 { HALF * (xp_f[0] + xp_r[0]) - - static_cast(Ip_f[0]) }; - const real_t Fx1_1 { (xp_r[0] - xp_i[0]) * coeff * inv_dt }; - const real_t Fx1_2 { (xp_f[0] - xp_r[0]) * coeff * inv_dt }; + const auto dxp_r_1 { static_cast(i1(p) == i1_prev(p)) * + (dx1(p) + dx1_prev(p)) * static_cast(INV_2) }; + + const real_t Wx1_1 { INV_2 * (dxp_r_1 + dx1_prev(p) + + static_cast(i1(p) > i1_prev(p))) }; + const real_t Wx1_2 { INV_2 * (dx1(p) + dxp_r_1 + + static_cast( + static_cast(i1(p) > i1_prev(p)) + + i1_prev(p) - i1(p))) }; + const real_t Fx1_1 { (static_cast(i1(p) > i1_prev(p)) + dxp_r_1 - + dx1_prev(p)) * + coeff * inv_dt }; + const real_t Fx1_2 { (static_cast( + i1(p) - i1_prev(p) - + static_cast(i1(p) > i1_prev(p))) + + dx1(p) - dxp_r_1) * + coeff * inv_dt }; auto J_acc = J.access(); + // tuple_t dxp_r; if constexpr (D == Dim::_1D) { const real_t Fx2_1 { HALF * vp[1] * coeff }; const real_t Fx2_2 { HALF * vp[1] * coeff }; @@ -146,265 +171,210 @@ namespace kernel { const real_t Fx3_1 { HALF * vp[2] * coeff }; const real_t Fx3_2 { HALF * vp[2] * coeff }; - J_acc(Ip_i[0] + N_GHOSTS, cur::jx1) += Fx1_1; - J_acc(Ip_f[0] + N_GHOSTS, cur::jx1) += Fx1_2; + J_acc(i1_prev(p) + N_GHOSTS, cur::jx1) += Fx1_1; + J_acc(i1(p) + N_GHOSTS, cur::jx1) += Fx1_2; - J_acc(Ip_i[0] + N_GHOSTS, cur::jx2) += Fx2_1 * (ONE - Wx1_1); - J_acc(Ip_i[0] + N_GHOSTS + 1, cur::jx2) += Fx2_1 * Wx1_1; - J_acc(Ip_f[0] + N_GHOSTS, cur::jx2) += Fx2_2 * (ONE - Wx1_2); - J_acc(Ip_f[0] + N_GHOSTS + 1, cur::jx2) += Fx2_2 * Wx1_2; + J_acc(i1_prev(p) + N_GHOSTS, cur::jx2) += Fx2_1 * (ONE - Wx1_1); + J_acc(i1_prev(p) + N_GHOSTS + 1, cur::jx2) += Fx2_1 * Wx1_1; + J_acc(i1(p) + N_GHOSTS, cur::jx2) += Fx2_2 * (ONE - Wx1_2); + J_acc(i1(p) + N_GHOSTS + 1, cur::jx2) += Fx2_2 * Wx1_2; - J_acc(Ip_i[0] + N_GHOSTS, cur::jx3) += Fx3_1 * (ONE - Wx1_1); - J_acc(Ip_i[0] + N_GHOSTS + 1, cur::jx3) += Fx3_1 * Wx1_1; - J_acc(Ip_f[0] + N_GHOSTS, cur::jx3) += Fx3_2 * (ONE - Wx1_2); - J_acc(Ip_f[0] + N_GHOSTS + 1, cur::jx3) += Fx3_2 * Wx1_2; + J_acc(i1_prev(p) + N_GHOSTS, cur::jx3) += Fx3_1 * (ONE - Wx1_1); + J_acc(i1_prev(p) + N_GHOSTS + 1, cur::jx3) += Fx3_1 * Wx1_1; + J_acc(i1(p) + N_GHOSTS, cur::jx3) += Fx3_2 * (ONE - Wx1_2); + J_acc(i1(p) + N_GHOSTS + 1, cur::jx3) += Fx3_2 * Wx1_2; } else if constexpr (D == Dim::_2D || D == Dim::_3D) { - const real_t Wx2_1 { HALF * (xp_i[1] + xp_r[1]) - - static_cast(Ip_i[1]) }; - const real_t Wx2_2 { HALF * (xp_f[1] + xp_r[1]) - - static_cast(Ip_f[1]) }; - const real_t Fx2_1 { (xp_r[1] - xp_i[1]) * coeff * inv_dt }; - const real_t Fx2_2 { (xp_f[1] - xp_r[1]) * coeff * inv_dt }; + const auto dxp_r_2 { static_cast(i2(p) == i2_prev(p)) * + (dx2(p) + dx2_prev(p)) * + static_cast(INV_2) }; + + const real_t Wx2_1 { INV_2 * (dxp_r_2 + dx2_prev(p) + + static_cast(i2(p) > i2_prev(p))) }; + const real_t Wx2_2 { INV_2 * (dx2(p) + dxp_r_2 + + static_cast( + static_cast(i2(p) > i2_prev(p)) + + i2_prev(p) - i2(p))) }; + const real_t Fx2_1 { (static_cast(i2(p) > i2_prev(p)) + + dxp_r_2 - dx2_prev(p)) * + coeff * inv_dt }; + const real_t Fx2_2 { (static_cast( + i2(p) - i2_prev(p) - + static_cast(i2(p) > i2_prev(p))) + + dx2(p) - dxp_r_2) * + coeff * inv_dt }; if constexpr (D == Dim::_2D) { const real_t Fx3_1 { HALF * vp[2] * coeff }; const real_t Fx3_2 { HALF * vp[2] * coeff }; - J_acc(Ip_i[0] + N_GHOSTS, Ip_i[1] + N_GHOSTS, cur::jx1) += Fx1_1 * - (ONE - Wx2_1); - J_acc(Ip_i[0] + N_GHOSTS, Ip_i[1] + N_GHOSTS + 1, cur::jx1) += Fx1_1 * - Wx2_1; - J_acc(Ip_f[0] + N_GHOSTS, Ip_f[1] + N_GHOSTS, cur::jx1) += Fx1_2 * - (ONE - Wx2_2); - J_acc(Ip_f[0] + N_GHOSTS, Ip_f[1] + N_GHOSTS + 1, cur::jx1) += Fx1_2 * - Wx2_2; - - J_acc(Ip_i[0] + N_GHOSTS, Ip_i[1] + N_GHOSTS, cur::jx2) += Fx2_1 * - (ONE - Wx1_1); - J_acc(Ip_i[0] + N_GHOSTS + 1, Ip_i[1] + N_GHOSTS, cur::jx2) += Fx2_1 * - Wx1_1; - J_acc(Ip_f[0] + N_GHOSTS, Ip_f[1] + N_GHOSTS, cur::jx2) += Fx2_2 * - (ONE - Wx1_2); - J_acc(Ip_f[0] + N_GHOSTS + 1, Ip_f[1] + N_GHOSTS, cur::jx2) += Fx2_2 * - Wx1_2; - - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS, + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS, + cur::jx1) += Fx1_1 * (ONE - Wx2_1); + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS + 1, + cur::jx1) += Fx1_1 * Wx2_1; + J_acc(i1(p) + N_GHOSTS, i2(p) + N_GHOSTS, cur::jx1) += Fx1_2 * + (ONE - Wx2_2); + J_acc(i1(p) + N_GHOSTS, i2(p) + N_GHOSTS + 1, cur::jx1) += Fx1_2 * Wx2_2; + + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS, + cur::jx2) += Fx2_1 * (ONE - Wx1_1); + J_acc(i1_prev(p) + N_GHOSTS + 1, + i2_prev(p) + N_GHOSTS, + cur::jx2) += Fx2_1 * Wx1_1; + J_acc(i1(p) + N_GHOSTS, i2(p) + N_GHOSTS, cur::jx2) += Fx2_2 * + (ONE - Wx1_2); + J_acc(i1(p) + N_GHOSTS + 1, i2(p) + N_GHOSTS, cur::jx2) += Fx2_2 * Wx1_2; + + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS, cur::jx3) += Fx3_1 * (ONE - Wx1_1) * (ONE - Wx2_1); - J_acc(Ip_i[0] + N_GHOSTS + 1, - Ip_i[1] + N_GHOSTS, + J_acc(i1_prev(p) + N_GHOSTS + 1, + i2_prev(p) + N_GHOSTS, cur::jx3) += Fx3_1 * Wx1_2 * (ONE - Wx2_1); - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS + 1, + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS + 1, cur::jx3) += Fx3_1 * (ONE - Wx1_1) * Wx2_1; - J_acc(Ip_i[0] + N_GHOSTS + 1, - Ip_i[1] + N_GHOSTS + 1, + J_acc(i1_prev(p) + N_GHOSTS + 1, + i2_prev(p) + N_GHOSTS + 1, cur::jx3) += Fx3_1 * Wx1_1 * Wx2_1; - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS, - cur::jx3) += Fx3_2 * (ONE - Wx1_2) * (ONE - Wx2_2); - J_acc(Ip_f[0] + N_GHOSTS + 1, - Ip_f[1] + N_GHOSTS, + J_acc(i1(p) + N_GHOSTS, i2(p) + N_GHOSTS, cur::jx3) += Fx3_2 * + (ONE - Wx1_2) * + (ONE - Wx2_2); + J_acc(i1(p) + N_GHOSTS + 1, + i2(p) + N_GHOSTS, cur::jx3) += Fx3_2 * Wx1_2 * (ONE - Wx2_2); - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS + 1, + J_acc(i1(p) + N_GHOSTS, + i2(p) + N_GHOSTS + 1, cur::jx3) += Fx3_2 * (ONE - Wx1_2) * Wx2_2; - J_acc(Ip_f[0] + N_GHOSTS + 1, - Ip_f[1] + N_GHOSTS + 1, - cur::jx3) += Fx3_2 * Wx1_2 * Wx2_2; + J_acc(i1(p) + N_GHOSTS + 1, i2(p) + N_GHOSTS + 1, cur::jx3) += Fx3_2 * + Wx1_2 * + Wx2_2; } else { - const real_t Wx3_1 { HALF * (xp_i[2] + xp_r[2]) - - static_cast(Ip_i[2]) }; - const real_t Wx3_2 { HALF * (xp_f[2] + xp_r[2]) - - static_cast(Ip_f[2]) }; - const real_t Fx3_1 { (xp_r[2] - xp_i[2]) * coeff * inv_dt }; - const real_t Fx3_2 { (xp_f[2] - xp_r[2]) * coeff * inv_dt }; - - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS, - Ip_i[2] + N_GHOSTS, + const auto dxp_r_3 { static_cast(i3(p) == i3_prev(p)) * + (dx3(p) + dx3_prev(p)) * + static_cast(INV_2) }; + const real_t Wx3_1 { INV_2 * (dxp_r_3 + dx3_prev(p) + + static_cast(i3(p) > i3_prev(p))) }; + const real_t Wx3_2 { INV_2 * (dx3(p) + dxp_r_3 + + static_cast( + static_cast(i3(p) > i3_prev(p)) + + i3_prev(p) - i3(p))) }; + const real_t Fx3_1 { (static_cast(i3(p) > i3_prev(p)) + + dxp_r_3 - dx3_prev(p)) * + coeff * inv_dt }; + const real_t Fx3_2 { (static_cast( + i3(p) - i3_prev(p) - + static_cast(i3(p) > i3_prev(p))) + + dx3(p) - dxp_r_3) * + coeff * inv_dt }; + + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS, + i3_prev(p) + N_GHOSTS, cur::jx1) += Fx1_1 * (ONE - Wx2_1) * (ONE - Wx3_1); - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS + 1, - Ip_i[2] + N_GHOSTS, + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS + 1, + i3_prev(p) + N_GHOSTS, cur::jx1) += Fx1_1 * Wx2_1 * (ONE - Wx3_1); - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS, - Ip_i[2] + N_GHOSTS + 1, + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS, + i3_prev(p) + N_GHOSTS + 1, cur::jx1) += Fx1_1 * (ONE - Wx2_1) * Wx3_1; - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS + 1, - Ip_i[2] + N_GHOSTS + 1, + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS + 1, + i3_prev(p) + N_GHOSTS + 1, cur::jx1) += Fx1_1 * Wx2_1 * Wx3_1; - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS, - Ip_f[2] + N_GHOSTS, + J_acc(i1(p) + N_GHOSTS, + i2(p) + N_GHOSTS, + i3(p) + N_GHOSTS, cur::jx1) += Fx1_2 * (ONE - Wx2_2) * (ONE - Wx3_2); - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS + 1, - Ip_f[2] + N_GHOSTS, + J_acc(i1(p) + N_GHOSTS, + i2(p) + N_GHOSTS + 1, + i3(p) + N_GHOSTS, cur::jx1) += Fx1_2 * Wx2_2 * (ONE - Wx3_2); - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS, - Ip_f[2] + N_GHOSTS + 1, + J_acc(i1(p) + N_GHOSTS, + i2(p) + N_GHOSTS, + i3(p) + N_GHOSTS + 1, cur::jx1) += Fx1_2 * (ONE - Wx2_2) * Wx3_2; - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS + 1, - Ip_f[2] + N_GHOSTS + 1, + J_acc(i1(p) + N_GHOSTS, + i2(p) + N_GHOSTS + 1, + i3(p) + N_GHOSTS + 1, cur::jx1) += Fx1_2 * Wx2_2 * Wx3_2; - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS, - Ip_i[2] + N_GHOSTS, + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS, + i3_prev(p) + N_GHOSTS, cur::jx2) += Fx2_1 * (ONE - Wx1_1) * (ONE - Wx3_1); - J_acc(Ip_i[0] + N_GHOSTS + 1, - Ip_i[1] + N_GHOSTS, - Ip_i[2] + N_GHOSTS, + J_acc(i1_prev(p) + N_GHOSTS + 1, + i2_prev(p) + N_GHOSTS, + i3_prev(p) + N_GHOSTS, cur::jx2) += Fx2_1 * Wx1_1 * (ONE - Wx3_1); - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS, - Ip_i[2] + N_GHOSTS + 1, + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS, + i3_prev(p) + N_GHOSTS + 1, cur::jx2) += Fx2_1 * (ONE - Wx1_1) * Wx3_1; - J_acc(Ip_i[0] + N_GHOSTS + 1, - Ip_i[1] + N_GHOSTS, - Ip_i[2] + N_GHOSTS + 1, + J_acc(i1_prev(p) + N_GHOSTS + 1, + i2_prev(p) + N_GHOSTS, + i3_prev(p) + N_GHOSTS + 1, cur::jx2) += Fx2_1 * Wx1_1 * Wx3_1; - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS, - Ip_f[2] + N_GHOSTS, + J_acc(i1(p) + N_GHOSTS, + i2(p) + N_GHOSTS, + i3(p) + N_GHOSTS, cur::jx2) += Fx2_2 * (ONE - Wx1_2) * (ONE - Wx3_2); - J_acc(Ip_f[0] + N_GHOSTS + 1, - Ip_f[1] + N_GHOSTS, - Ip_f[2] + N_GHOSTS, + J_acc(i1(p) + N_GHOSTS + 1, + i2(p) + N_GHOSTS, + i3(p) + N_GHOSTS, cur::jx2) += Fx2_2 * Wx1_2 * (ONE - Wx3_2); - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS, - Ip_f[2] + N_GHOSTS + 1, + J_acc(i1(p) + N_GHOSTS, + i2(p) + N_GHOSTS, + i3(p) + N_GHOSTS + 1, cur::jx2) += Fx2_2 * (ONE - Wx1_2) * Wx3_2; - J_acc(Ip_f[0] + N_GHOSTS + 1, - Ip_f[1] + N_GHOSTS, - Ip_f[2] + N_GHOSTS + 1, + J_acc(i1(p) + N_GHOSTS + 1, + i2(p) + N_GHOSTS, + i3(p) + N_GHOSTS + 1, cur::jx2) += Fx2_2 * Wx1_2 * Wx3_2; - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS, - Ip_i[2] + N_GHOSTS, + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS, + i3_prev(p) + N_GHOSTS, cur::jx3) += Fx3_1 * (ONE - Wx1_1) * (ONE - Wx2_1); - J_acc(Ip_i[0] + N_GHOSTS + 1, - Ip_i[1] + N_GHOSTS, - Ip_i[2] + N_GHOSTS, + J_acc(i1_prev(p) + N_GHOSTS + 1, + i2_prev(p) + N_GHOSTS, + i3_prev(p) + N_GHOSTS, cur::jx3) += Fx3_1 * Wx1_1 * (ONE - Wx2_1); - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS + 1, - Ip_i[2] + N_GHOSTS, + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS + 1, + i3_prev(p) + N_GHOSTS, cur::jx3) += Fx3_1 * (ONE - Wx1_1) * Wx2_1; - J_acc(Ip_i[0] + N_GHOSTS + 1, - Ip_i[1] + N_GHOSTS + 1, - Ip_i[2] + N_GHOSTS, + J_acc(i1_prev(p) + N_GHOSTS + 1, + i2_prev(p) + N_GHOSTS + 1, + i3_prev(p) + N_GHOSTS, cur::jx3) += Fx3_1 * Wx1_1 * Wx2_1; - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS, - Ip_f[2] + N_GHOSTS, + J_acc(i1(p) + N_GHOSTS, + i2(p) + N_GHOSTS, + i3(p) + N_GHOSTS, cur::jx3) += Fx3_2 * (ONE - Wx1_2) * (ONE - Wx2_2); - J_acc(Ip_f[0] + N_GHOSTS + 1, - Ip_f[1] + N_GHOSTS, - Ip_f[2] + N_GHOSTS, + J_acc(i1(p) + N_GHOSTS + 1, + i2(p) + N_GHOSTS, + i3(p) + N_GHOSTS, cur::jx3) += Fx3_2 * Wx1_2 * (ONE - Wx2_2); - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS + 1, - Ip_f[2] + N_GHOSTS, + J_acc(i1(p) + N_GHOSTS, + i2(p) + N_GHOSTS + 1, + i3(p) + N_GHOSTS, cur::jx3) += Fx3_2 * (ONE - Wx1_2) * Wx2_2; - J_acc(Ip_f[0] + N_GHOSTS + 1, - Ip_f[1] + N_GHOSTS + 1, - Ip_f[2] + N_GHOSTS, + J_acc(i1(p) + N_GHOSTS + 1, + i2(p) + N_GHOSTS + 1, + i3(p) + N_GHOSTS, cur::jx3) += Fx3_2 * Wx1_2 * Wx2_2; } } } - - /** - * @brief Get particle position in `coord_t` form. - * @param[in] p Index of particle. - * @param[out] Ip_f Final position of the particle (cell index). - * @param[out] Ip_i Initial position of the particle (cell index). - * @param[out] xp_f Final position. - * @param[out] xp_i Previous step position. - * @param[out] xp_r Intermediate point used in zig-zag deposit. - */ - Inline auto getDepositInterval(index_t& p, - tuple_t& Ip_f, - tuple_t& Ip_i, - coord_t& xp_f, - coord_t& xp_i, - coord_t& xp_r) const -> void { - Ip_f[0] = i1(p); - Ip_i[0] = i1_prev(p); - xp_f[0] = i_di_to_Xi(Ip_f[0], dx1(p)); - xp_i[0] = i_di_to_Xi(Ip_i[0], dx1_prev(p)); - if constexpr (D == Dim::_2D || D == Dim::_3D) { - Ip_f[1] = i2(p); - Ip_i[1] = i2_prev(p); - xp_f[1] = i_di_to_Xi(Ip_f[1], dx2(p)); - xp_i[1] = i_di_to_Xi(Ip_i[1], dx2_prev(p)); - } - if constexpr (D == Dim::_3D) { - Ip_f[2] = i3(p); - Ip_i[2] = i3_prev(p); - xp_f[2] = i_di_to_Xi(Ip_f[2], dx3(p)); - xp_i[2] = i_di_to_Xi(Ip_i[2], dx3_prev(p)); - } - for (auto i = 0u; i < D; ++i) { - xp_r[i] = math::fmin(static_cast(IMIN(Ip_i[i], Ip_f[i]) + 1), - math::fmax(static_cast(IMAX(Ip_i[i], Ip_f[i])), - HALF * (xp_i[i] + xp_f[i]))); - } - } - - // Getters - Inline void getPrtlPos(index_t& p, coord_t& xp) const { - if constexpr (D == Dim::_1D) { - xp[0] = i_di_to_Xi(i1(p), dx1(p)); - } else if constexpr (D == Dim::_2D) { - if constexpr (M::PrtlDim == Dim::_3D) { - xp[0] = i_di_to_Xi(i1(p), dx1(p)); - xp[1] = i_di_to_Xi(i2(p), dx2(p)); - xp[2] = phi(p); - } else { - xp[0] = i_di_to_Xi(i1(p), dx1(p)); - xp[1] = i_di_to_Xi(i2(p), dx2(p)); - } - } else { - xp[0] = i_di_to_Xi(i1(p), dx1(p)); - xp[1] = i_di_to_Xi(i2(p), dx2(p)); - xp[2] = i_di_to_Xi(i3(p), dx3(p)); - } - } - - Inline void getPrtl3Vel(index_t& p, vec_t& vp) const { - coord_t xp { ZERO }; - getPrtlPos(p, xp); - auto inv_energy { ZERO }; - if constexpr (S == SimEngine::SRPIC) { - metric.template transform_xyz(xp, - { ux1(p), ux2(p), ux3(p) }, - vp); - inv_energy = ONE / math::sqrt(ONE + NORM_SQR(ux1(p), ux2(p), ux3(p))); - } else { - metric.template transform(xp, { ux1(p), ux2(p), ux3(p) }, vp); - inv_energy = ONE / math::sqrt(ONE + ux1(p) * vp[0] + ux2(p) * vp[1] + - ux3(p) * vp[2]); - } - if (Kokkos::isnan(vp[2]) || Kokkos::isinf(vp[2])) { - vp[2] = ZERO; - } - vp[0] *= inv_energy; - vp[1] *= inv_energy; - vp[2] *= inv_energy; - } }; } // namespace kernel From 9a17030497bddee6ce0b9bcce93ecf00f9f17036 Mon Sep 17 00:00:00 2001 From: hayk Date: Wed, 9 Apr 2025 15:58:31 -0400 Subject: [PATCH 05/25] keepconstinjector works --- setups/wip/reconnection/pgen.hpp | 159 +++++++++++----------- src/archetypes/particle_injector.h | 209 +++++++++++++++-------------- src/kernels/injectors.hpp | 16 +-- 3 files changed, 197 insertions(+), 187 deletions(-) diff --git a/setups/wip/reconnection/pgen.hpp b/setups/wip/reconnection/pgen.hpp index 34314ca68..7ae550756 100644 --- a/setups/wip/reconnection/pgen.hpp +++ b/setups/wip/reconnection/pgen.hpp @@ -159,8 +159,9 @@ namespace user { using arch::ProblemGenerator::C; using arch::ProblemGenerator::params; - const real_t bg_B, bg_Bguide, bg_temperature; + const real_t bg_B, bg_Bguide, bg_temperature, inj_ypad; const real_t cs_width, cs_overdensity, cs_x, cs_y; + const real_t ymin, ymax; InitFields init_flds; @@ -169,6 +170,7 @@ namespace user { , bg_B { p.template get("setup.bg_B", 1.0) } , bg_Bguide { p.template get("setup.bg_Bguide", 0.0) } , bg_temperature { p.template get("setup.bg_temperature", 0.001) } + , inj_ypad { p.template get("setup.inj_ypad", (real_t)0.05) } , cs_width { p.template get("setup.cs_width") } , cs_overdensity { p.template get("setup.cs_overdensity") } , cs_x { m.mesh().extent(in::x1).first + @@ -177,6 +179,8 @@ namespace user { , cs_y { m.mesh().extent(in::x2).first + INV_2 * (m.mesh().extent(in::x2).second - m.mesh().extent(in::x2).first) } + , ymin { m.mesh().extent(in::x2).first } + , ymax { m.mesh().extent(in::x2).second } , init_flds { bg_B, bg_Bguide, cs_width, cs_y } {} inline PGen() {} @@ -232,81 +236,84 @@ namespace user { cs_overdensity); } - void CustomPostStep(std::size_t step, long double time, Domain& domain) { - // // 0. define target density profile and box where it has to be - // reached - // // 1. compute density - // // 2. define spatial distribution - // // 3. define energy distribution - // // 4. define particle injector - // // 5. inject particles - // - // // step 0. - // // defining the regions of interest (lower and upper - // y-boundaries) boundaries_t box_upper, box_lower; const - // auto [xmin, xmax] = domain.mesh.extent(in::x1); const auto [ymin, - // ymax] = domain.mesh.extent(in::x2); box_upper.push_back({ xmin, - // xmax }); box_lower.push_back({ xmin, xmax }); - // - // box_upper.push_back({ ymax - dy, ymax }); - // box_lower.push_back({ ymin, ymin + dy }); - // - // if constexpr (M::Dim == Dim::_3D) { - // const auto [zmin, zmax] = domain.mesh.extent(in::x3); - // box_upper.push_back({ zmin, zmax }); - // box_lower.push_back({ zmin, zmax }); - // } - // - // const auto const_dens = ConstDens(); - // const auto inv_n0 = ONE / n0; - // - // // step 1 compute density - // auto scatter_bckp = Kokkos::Experimental::create_scatter_view( - // domain.fields.bckp); - // for (auto& prtl_spec : domain.species) { - // // 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(), - // false, - // domain.mesh.metric, - // domain.mesh.flds_bc(), - // 0, inv_n0, 0)); - // } - // Kokkos::Experimental::contribute(domain.fields.bckp, scatter_bckp); - // - // - // //step 2. define spatial distribution - // // const auto spatial_dist = arch::Replenish>(domain.mesh.metric, domain.fields.bckp, 0, - // const_dens, ONE); - // const auto spatial_dist = spatial_dist_t(domain.mesh.metric, - // domain.fields.bckp,0,const_dens,ONE); - // //step 3. define energy distribution - // const auto energy_dist = arch::Maxwellian(domain.mesh.metric, - // domain.random_pool, up_temperature); - // //step 4. define particle injector - // const auto injector = arch::NonUniformInjector(energy_dist, spatial_dist, {1,2}); - // //step 5. inject particles - // arch::InjectNonUniform(params, domain, - // injector, ONE, false, box_upper); //upper boudary arch::InjectNonUniform(params, domain, injector, ONE, false, - // box_lower); //lower boundary - // - // + void CustomPostStep(timestep_t, simtime_t, Domain& domain) { + const auto energy_dist = arch::Maxwellian(domain.mesh.metric, + domain.random_pool, + bg_temperature); + + const auto dx = domain.mesh.metric.template sqrt_h_<1, 1>({}); + + boundaries_t inj_box_up, inj_box_down; + boundaries_t probe_box_up, probe_box_down; + inj_box_up.push_back(Range::All); + inj_box_down.push_back(Range::All); + probe_box_up.push_back(Range::All); + probe_box_down.push_back(Range::All); + inj_box_up.push_back({ ymax - inj_ypad - 10 * dx, ymax - inj_ypad }); + inj_box_down.push_back({ ymin + inj_ypad, ymin + inj_ypad + 10 * dx }); + probe_box_up.push_back({ ymax - inj_ypad - 10 * dx, ymax - inj_ypad }); + probe_box_down.push_back({ ymin + inj_ypad, ymin + inj_ypad + 10 * dx }); + + if constexpr (M::Dim == Dim::_3D) { + inj_box_up.push_back(Range::All); + 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); + } + + const auto injector_up = arch::KeepConstantInjector( + energy_dist, + { 1, 2 }, + 0u, + probe_box_up); + const auto injector_down = arch::KeepConstantInjector( + energy_dist, + { 1, 2 }, + 0u, + probe_box_down); + + arch::InjectUniform( + params, + domain, + injector_up, + ONE, + params.template get("particles.use_weights"), + inj_box_up); + arch::InjectUniform( + params, + domain, + injector_down, + ONE, + params.template get("particles.use_weights"), + inj_box_down); } }; diff --git a/src/archetypes/particle_injector.h b/src/archetypes/particle_injector.h index 39c588519..88b0db58b 100644 --- a/src/archetypes/particle_injector.h +++ b/src/archetypes/particle_injector.h @@ -66,45 +66,43 @@ namespace arch { ~UniformInjector() = default; - auto DeduceInjectRegion(const Domain& domain, - const boundaries_t& box) const - -> std::tuple, array_t> { - auto i_min = array_t { "i_min", M::Dim }; - auto i_max = array_t { "i_max", M::Dim }; - + auto DeduceRegion(const Domain& domain, const boundaries_t& box) const + -> std::tuple, array_t> { if (not domain.mesh.Intersects(box)) { - return { false, (ncells_t)0, i_min, i_max }; + return { false, array_t {}, array_t {} }; } + coord_t xCorner_min_Ph { ZERO }; + coord_t xCorner_max_Ph { ZERO }; + coord_t xCorner_min_Cd { ZERO }; + coord_t xCorner_max_Cd { ZERO }; - tuple_t range_min { 0 }; - tuple_t range_max { 0 }; - - boundaries_t incl_ghosts; for (auto d { 0u }; d < M::Dim; ++d) { - incl_ghosts.push_back({ false, false }); + const auto local_xi_min = domain.mesh.extent(static_cast(d)).first; + const auto local_xi_max = domain.mesh.extent(static_cast(d)).second; + const auto extent_min = std::min(std::max(local_xi_min, box[d].first), + local_xi_max); + const auto extent_max = std::max(std::min(local_xi_max, box[d].second), + local_xi_min); + xCorner_min_Ph[d] = extent_min; + xCorner_max_Ph[d] = extent_max; } - const auto intersect_range = domain.mesh.ExtentToRange(box, incl_ghosts); + domain.mesh.metric.template convert(xCorner_min_Ph, + xCorner_min_Cd); + domain.mesh.metric.template convert(xCorner_max_Ph, + xCorner_max_Cd); - for (auto d { 0u }; d < M::Dim; ++d) { - range_min[d] = intersect_range[d].first; - range_max[d] = intersect_range[d].second; - } + array_t xi_min { "xi_min", M::Dim }, xi_max { "xi_max", M::Dim }; - ncells_t ncells = 1; - for (auto d = 0u; d < M::Dim; ++d) { - ncells *= (range_max[d] - range_min[d]); - } - - auto i_min_h = Kokkos::create_mirror_view(i_min); - auto i_max_h = Kokkos::create_mirror_view(i_max); - for (auto d = 0u; d < M::Dim; ++d) { - i_min_h(d) = (real_t)(range_min[d]); - i_max_h(d) = (real_t)(range_max[d]); + auto xi_min_h = Kokkos::create_mirror_view(xi_min); + auto xi_max_h = Kokkos::create_mirror_view(xi_max); + for (auto d { 0u }; d < M::Dim; ++d) { + xi_min_h(d) = xCorner_min_Cd[d]; + xi_max_h(d) = xCorner_max_Cd[d]; } + Kokkos::deep_copy(xi_min, xi_min_h); + Kokkos::deep_copy(xi_max, xi_max_h); - Kokkos::deep_copy(i_min, i_min_h); - Kokkos::deep_copy(i_max, i_max_h); - return { true, ncells, i_min, i_max }; + return { true, xi_min, xi_max }; } auto ComputeNumInject(const SimulationParams& params, @@ -112,20 +110,28 @@ namespace arch { real_t number_density, const boundaries_t& box) const -> std::tuple, array_t> { - const auto result = DeduceInjectRegion(domain, box); - const auto i_min = std::get<2>(result); - const auto i_max = std::get<3>(result); - + const auto result = DeduceRegion(domain, box); if (not std::get<0>(result)) { - return { false, (npart_t)0, i_min, i_max }; + return { false, (npart_t)0, array_t {}, array_t {} }; + } + const auto xi_min = std::get<1>(result); + const auto xi_max = std::get<2>(result); + auto xi_min_h = Kokkos::create_mirror_view(xi_min); + auto xi_max_h = Kokkos::create_mirror_view(xi_max); + Kokkos::deep_copy(xi_min_h, xi_min); + Kokkos::deep_copy(xi_max_h, xi_max); + + long double num_cells { 1.0 }; + for (auto d { 0u }; d < M::Dim; ++d) { + num_cells *= static_cast(xi_max_h(d)) - + static_cast(xi_min_h(d)); } - const auto ncells = std::get<1>(result); const auto ppc0 = params.template get("particles.ppc0"); const auto nparticles = static_cast( - (long double)(ppc0 * number_density * 0.5) * (long double)(ncells)); + (long double)(ppc0 * number_density * 0.5) * num_cells); - return { true, nparticles, i_min, i_max }; + return { true, nparticles, xi_min, xi_max }; } }; @@ -135,12 +141,15 @@ namespace arch { using UniformInjector::D; using UniformInjector::C; + const unsigned short density_buff_idx; boundaries_t probe_box; KeepConstantInjector(const energy_dist_t& energy_dist, const std::pair& species, + unsigned short density_buff_idx, boundaries_t box = {}) - : UniformInjector { energy_dist, species } { + : UniformInjector { energy_dist, species } + , density_buff_idx { density_buff_idx } { for (auto d { 0u }; d < M::Dim; ++d) { if (d < box.size()) { probe_box.push_back({ box[d].first, box[d].second }); @@ -153,98 +162,92 @@ namespace arch { ~KeepConstantInjector() = default; auto ComputeAvgDensity(const SimulationParams& params, - const Domain& domain, - boundaries_t box) const -> real_t { - 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); - - for (const auto& sp : std::vector { this->species.first, - this->species.second }) { - raise::ErrorIf(sp >= domain.species.size(), - "Species index out of bounds", - HERE); - const auto& prtl_spec = domain.species[sp - 1]; - Kokkos::deep_copy(domain.fields.buff, ZERO); - // clang-format off - Kokkos::parallel_for( - "ComputeMoments", - prtl_spec.rangeActiveParticles(), - kernel::ParticleMoments_kernel({}, scatter_buff, 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(), - use_weights, - domain.mesh.metric, domain.mesh.flds_bc(), - ni2, inv_n0, 0)); - // clang-format on + Domain& domain) const -> real_t { + const auto result = this->DeduceRegion(domain, probe_box); + const auto should_probe = std::get<0>(result); + if (not should_probe) { + return ZERO; } - Kokkos::Experimental::contribute(domain.fields.buff, scatter_buff); + const auto xi_min_arr = std::get<1>(result); + const auto xi_max_arr = std::get<2>(result); - real_t dens { ZERO }; + tuple_t i_min { 0 }; + tuple_t i_max { 0 }; - const auto result = DeduceInjectRegion(domain, box); - const auto should_probe = std::get<0>(result); - const auto ncells = std::get<1>(result); - const auto i_min = std::get<2>(result); - const auto i_max = std::get<3>(result); + auto xi_min_h = Kokkos::create_mirror_view(xi_min_arr); + auto xi_max_h = Kokkos::create_mirror_view(xi_max_arr); + Kokkos::deep_copy(xi_min_h, xi_min_arr); + Kokkos::deep_copy(xi_max_h, xi_max_arr); + + ncells_t num_cells = 1u; + for (auto d { 0u }; d < M::Dim; ++d) { + i_min[d] = std::floor(xi_min_h(d)) + N_GHOSTS; + i_max[d] = std::ceil(xi_max_h(d)) + N_GHOSTS; + num_cells *= (i_max[d] - i_min[d]); + } + real_t dens { ZERO }; if (should_probe) { - range_t probe_range = CreateRangePolicy(i_min, i_max); Kokkos::parallel_reduce( "AvgDensity", - probe_range, - kernel::ComputeSum_kernel(domain.fields.buff, 0), + CreateRangePolicy(i_min, i_max), + kernel::ComputeSum_kernel(domain.fields.buff, density_buff_idx), dens); } #if defined(MPI_ENABLED) real_t tot_dens { ZERO }; - ncells_t tot_ncells { 0 }; + ncells_t tot_num_cells { 0 }; MPI_Allreduce(&dens, &tot_dens, 1, mpi::get_type(), MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&ncells, - &tot_ncells, + MPI_Allreduce(&num_cells, + &tot_num_cells, 1, mpi::get_type(), MPI_SUM, MPI_COMM_WORLD); - dens = tot_dens; - ncells = tot_ncells; + dens = tot_dens; + num_cells = tot_num_cells; #endif - if (ncells > 0) { - return dens / (real_t)(ncells); + if (num_cells > 0) { + return dens / (real_t)(num_cells); } else { return ZERO; } } auto ComputeNumInject(const SimulationParams& params, - const Domain& domain, + Domain& domain, real_t number_density, const boundaries_t& box) const -> std::tuple, array_t> { const auto computed_avg_density = ComputeAvgDensity(params, domain); - const auto result = DeduceInjectRegion(domain, box); - const auto i_min = std::get<2>(result); - const auto i_max = std::get<3>(result); - + const auto result = this->DeduceRegion(domain, box); if (not std::get<0>(result)) { - return { false, (npart_t)0, i_min, i_max }; + return { false, (npart_t)0, array_t {}, array_t {} }; } - const auto ncells = std::get<1>(result); - const auto ppc0 = params.template get("particles.ppc0"); - const auto nparticles = static_cast( - (long double)(ppc0 * (number_density - computed_avg_density) * 0.5) * - (long double)(ncells)); + const auto xi_min = std::get<1>(result); + const auto xi_max = std::get<2>(result); + auto xi_min_h = Kokkos::create_mirror_view(xi_min); + auto xi_max_h = Kokkos::create_mirror_view(xi_max); + Kokkos::deep_copy(xi_min_h, xi_min); + Kokkos::deep_copy(xi_max_h, xi_max); + + long double num_cells { 1.0 }; + for (auto d { 0u }; d < M::Dim; ++d) { + num_cells *= static_cast(xi_max_h(d)) - + static_cast(xi_min_h(d)); + } + + const auto ppc0 = params.template get("particles.ppc0"); + npart_t nparticles { 0u }; + if (number_density > computed_avg_density) { + nparticles = static_cast( + (long double)(ppc0 * (number_density - computed_avg_density) * 0.5) * + num_cells); + } - return { true, nparticles, i_min, i_max }; + return { nparticles != 0u, nparticles, xi_min, xi_max }; } }; @@ -471,8 +474,8 @@ namespace arch { return; } const auto nparticles = std::get<1>(result); - const auto i_min = std::get<2>(result); - const auto i_max = std::get<3>(result); + const auto xi_min = std::get<2>(result); + const auto xi_max = std::get<3>(result); Kokkos::parallel_for( "InjectUniform", @@ -485,8 +488,8 @@ namespace arch { domain.species[injector.species.first - 1].npart(), domain.species[injector.species.second - 1].npart(), domain.mesh.metric, - i_min, - i_max, + xi_min, + xi_max, injector.energy_dist, ONE / params.template get("scales.V0"), domain.random_pool)); diff --git a/src/kernels/injectors.hpp b/src/kernels/injectors.hpp index 7dcd017a3..d50b5c24c 100644 --- a/src/kernels/injectors.hpp +++ b/src/kernels/injectors.hpp @@ -49,7 +49,7 @@ namespace kernel { npart_t offset1, offset2; const M metric; - const array_t i_min, i_max; + const array_t xi_min, xi_max; const ED energy_dist; const real_t inv_V0; random_number_pool_t random_pool; @@ -61,8 +61,8 @@ namespace kernel { npart_t offset1, npart_t offset2, const M& metric, - const array_t& i_min, - const array_t& i_max, + const array_t& xi_min, + const array_t& xi_max, const ED& energy_dist, real_t inv_V0, random_number_pool_t& random_pool) @@ -95,8 +95,8 @@ namespace kernel { , offset1 { offset1 } , offset2 { offset2 } , metric { metric } - , i_min { i_min } - , i_max { i_max } + , xi_min { xi_min } + , xi_max { xi_max } , energy_dist { energy_dist } , inv_V0 { inv_V0 } , random_pool { random_pool } {} @@ -106,12 +106,12 @@ namespace kernel { vec_t v1 { ZERO }, v2 { ZERO }; { // generate a random coordinate auto rand_gen = random_pool.get_state(); - x_Cd[0] = i_min(0) + Random(rand_gen) * (i_max(0) - i_min(0)); + x_Cd[0] = xi_min(0) + Random(rand_gen) * (xi_max(0) - xi_min(0)); if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { - x_Cd[1] = i_min(1) + Random(rand_gen) * (i_max(1) - i_min(1)); + x_Cd[1] = xi_min(1) + Random(rand_gen) * (xi_max(1) - xi_min(1)); } if constexpr (M::Dim == Dim::_3D) { - x_Cd[2] = i_min(2) + Random(rand_gen) * (i_max(2) - i_min(2)); + x_Cd[2] = xi_min(2) + Random(rand_gen) * (xi_max(2) - xi_min(2)); } random_pool.free_state(rand_gen); } From c2d44a22bf93b2a28bda391798100a02a79ac026 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludwig=20B=C3=B6ss?= Date: Sat, 19 Apr 2025 10:03:37 -0500 Subject: [PATCH 06/25] bugfixes --- setups/srpic/shock/pgen.hpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/setups/srpic/shock/pgen.hpp b/setups/srpic/shock/pgen.hpp index 738072f27..e354262b7 100644 --- a/setups/srpic/shock/pgen.hpp +++ b/setups/srpic/shock/pgen.hpp @@ -92,7 +92,6 @@ namespace user { // magnetic field properties real_t Btheta, Bphi, Bmag; InitFields init_flds; - BCFields bc_flds; inline PGen(const SimulationParams& p, const Metadomain& global_domain) : arch::ProblemGenerator { p } @@ -105,7 +104,6 @@ namespace user { , Btheta { p.template get("setup.Btheta", ZERO) } , Bphi { p.template get("setup.Bphi", ZERO) } , init_flds { Bmag, Btheta, Bphi, drift_ux } - , bc_flds { Bmag, Btheta, Bphi, drift_ux } , filling_fraction { p.template get("setup.filling_fraction", 1.0) } , injector_velocity { p.template get("setup.injector_velocity", 1.0) } , injection_start { p.template get("setup.injection_start", 0.0) } @@ -114,7 +112,7 @@ namespace user { inline PGen() {} - auto MatchFields(real_t time) const -> BCFields { + auto MatchFields(real_t time) const -> InitFields { return init_flds; } From a8a32061cde710cde30a39d54f4400db71d53d5c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludwig=20B=C3=B6ss?= Date: Sat, 19 Apr 2025 10:44:57 -0500 Subject: [PATCH 07/25] bugfix in temp calculation --- setups/srpic/shock/pgen.hpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/setups/srpic/shock/pgen.hpp b/setups/srpic/shock/pgen.hpp index e354262b7..6bd6f21a9 100644 --- a/setups/srpic/shock/pgen.hpp +++ b/setups/srpic/shock/pgen.hpp @@ -176,8 +176,7 @@ namespace user { const auto energy_dist = arch::TwoTemperatureMaxwellian( local_domain.mesh.metric, local_domain.random_pool, - { temperature_ratio * temperature * - (local_domain.species[2].mass() / local_domain.species[1].mass()), + { temperature_ratio * temperature * local_domain.species[1].mass() , temperature }, { 1, 2 }, -drift_ux, @@ -319,8 +318,7 @@ namespace user { const auto energy_dist = arch::TwoTemperatureMaxwellian( domain.mesh.metric, domain.random_pool, - { temperature_ratio * temperature * - (domain.species[2].mass() / domain.species[1].mass()), + { temperature_ratio * temperature * domain.species[1].mass(), temperature }, { 1, 2 }, -drift_ux, From 62d1ce432c8344507ddfd84c2339034e24c7a6e5 Mon Sep 17 00:00:00 2001 From: hayk Date: Sat, 19 Apr 2025 20:55:13 -0400 Subject: [PATCH 08/25] change BCs routine --- src/framework/domain/domain.h | 24 ++--- src/framework/domain/metadomain.cpp | 136 ++++++++++++++++++++++++++++ src/framework/domain/metadomain.h | 2 + 3 files changed, 147 insertions(+), 15 deletions(-) diff --git a/src/framework/domain/domain.h b/src/framework/domain/domain.h index 7966cdb54..c26f5995f 100644 --- a/src/framework/domain/domain.h +++ b/src/framework/domain/domain.h @@ -146,7 +146,8 @@ namespace ntt { } /* setters -------------------------------------------------------------- */ - auto set_neighbor_idx(const dir::direction_t& dir, unsigned int idx) -> void { + auto set_neighbor_idx(const dir::direction_t& dir, unsigned int idx) + -> void { m_neighbor_idx[dir] = idx; } @@ -164,8 +165,8 @@ namespace ntt { }; template - inline auto operator<<(std::ostream& os, - const Domain& domain) -> std::ostream& { + inline auto operator<<(std::ostream& os, const Domain& domain) + -> std::ostream& { os << "Domain #" << domain.index(); #if defined(MPI_ENABLED) os << " [MPI rank: " << domain.mpi_rank() << "]"; @@ -184,23 +185,16 @@ namespace ntt { } os << "\n"; os << std::setw(19) << std::left << " physical extent: "; - for (auto dim = 0; dim < M::Dim; ++dim) { + for (auto dim { 0u }; dim < M::Dim; ++dim) { os << std::setw(15) << std::left << fmt::format("{%.2f; %.2f}", - domain.mesh.extent(dim).first, - domain.mesh.extent(dim).second); + domain.mesh.extent(static_cast(dim)).first, + domain.mesh.extent(static_cast(dim)).second); } os << "\n neighbors:\n"; for (auto& direction : dir::Directions::all) { - auto neighbor = domain.neighbor_in(direction); - os << " " << direction; - if (neighbor != nullptr) { - os << " -> #" << neighbor->index() << "\n"; - } else { - os << " -> " - << "N/A" - << "\n"; - } + auto neighbor_idx = domain.neighbor_idx_in(direction); + os << " " << direction << " -> #" << neighbor_idx << "\n"; } os << " field boundaries:\n"; for (auto& direction : dir::Directions::orth) { diff --git a/src/framework/domain/metadomain.cpp b/src/framework/domain/metadomain.cpp index fc3c59d4c..4b9057e23 100644 --- a/src/framework/domain/metadomain.cpp +++ b/src/framework/domain/metadomain.cpp @@ -403,6 +403,142 @@ namespace ntt { #endif } + template + void Metadomain::setFldsBC(const bc_in& dir, const FldsBC& new_bcs) { + if (dir == bc_in::Mx1) { + if constexpr (M::Dim == Dim::_1D) { + g_mesh.set_flds_bc({ -1 }, new_bcs); + } else if constexpr (M::Dim == Dim::_2D) { + g_mesh.set_flds_bc({ -1, 0 }, new_bcs); + } else if constexpr (M::Dim == Dim::_3D) { + g_mesh.set_flds_bc({ -1, 0, 0 }, new_bcs); + } else { + raise::Error("Invalid dimension", HERE); + } + } else if (dir == bc_in::Px1) { + if constexpr (M::Dim == Dim::_1D) { + g_mesh.set_flds_bc({ +1 }, new_bcs); + } else if constexpr (M::Dim == Dim::_2D) { + g_mesh.set_flds_bc({ +1, 0 }, new_bcs); + } else if constexpr (M::Dim == Dim::_3D) { + g_mesh.set_flds_bc({ +1, 0, 0 }, new_bcs); + } else { + raise::Error("Invalid dimension", HERE); + } + } else if (dir == bc_in::Mx2) { + if constexpr (M::Dim == Dim::_1D) { + raise::Error("Cannot set -x2 BCs for 1D", HERE); + } else if constexpr (M::Dim == Dim::_2D) { + g_mesh.set_flds_bc({ -1, 0 }, new_bcs); + } else if constexpr (M::Dim == Dim::_3D) { + g_mesh.set_flds_bc({ -1, 0, 0 }, new_bcs); + } else { + raise::Error("Invalid dimension", HERE); + } + } else if (dir == bc_in::Px2) { + if constexpr (M::Dim == Dim::_1D) { + raise::Error("Cannot set +x2 BCs for 1D", HERE); + } else if constexpr (M::Dim == Dim::_2D) { + g_mesh.set_flds_bc({ +1, 0 }, new_bcs); + } else if constexpr (M::Dim == Dim::_3D) { + g_mesh.set_flds_bc({ +1, 0, 0 }, new_bcs); + } else { + raise::Error("Invalid dimension", HERE); + } + } else if (dir == bc_in::Mx3) { + if constexpr (M::Dim == Dim::_1D) { + raise::Error("Cannot set -x3 BCs for 1D", HERE); + } else if constexpr (M::Dim == Dim::_2D) { + raise::Error("Cannot set -x3 BCs for 2D", HERE); + } else if constexpr (M::Dim == Dim::_3D) { + g_mesh.set_flds_bc({ 0, 0, -1 }, new_bcs); + } else { + raise::Error("Invalid dimension", HERE); + } + } else if (dir == bc_in::Px3) { + if constexpr (M::Dim == Dim::_1D) { + raise::Error("Cannot set +x3 BCs for 1D", HERE); + } else if constexpr (M::Dim == Dim::_2D) { + raise::Error("Cannot set +x3 BCs for 2D", HERE); + } else if constexpr (M::Dim == Dim::_3D) { + g_mesh.set_flds_bc({ 0, 0, +1 }, new_bcs); + } else { + raise::Error("Invalid dimension", HERE); + } + } else { + raise::Error("Invalid direction", HERE); + } + redefineBoundaries(); + } + + template + void Metadomain::setPrtlBC(const bc_in& dir, const PrtlBC& new_bcs) { + if (dir == bc_in::Mx1) { + if constexpr (M::Dim == Dim::_1D) { + g_mesh.set_prtl_bc({ -1 }, new_bcs); + } else if constexpr (M::Dim == Dim::_2D) { + g_mesh.set_prtl_bc({ -1, 0 }, new_bcs); + } else if constexpr (M::Dim == Dim::_3D) { + g_mesh.set_prtl_bc({ -1, 0, 0 }, new_bcs); + } else { + raise::Error("Invalid dimension", HERE); + } + } else if (dir == bc_in::Px1) { + if constexpr (M::Dim == Dim::_1D) { + g_mesh.set_prtl_bc({ +1 }, new_bcs); + } else if constexpr (M::Dim == Dim::_2D) { + g_mesh.set_prtl_bc({ +1, 0 }, new_bcs); + } else if constexpr (M::Dim == Dim::_3D) { + g_mesh.set_prtl_bc({ +1, 0, 0 }, new_bcs); + } else { + raise::Error("Invalid dimension", HERE); + } + } else if (dir == bc_in::Mx2) { + if constexpr (M::Dim == Dim::_1D) { + raise::Error("Cannot set -x2 BCs for 1D", HERE); + } else if constexpr (M::Dim == Dim::_2D) { + g_mesh.set_prtl_bc({ -1, 0 }, new_bcs); + } else if constexpr (M::Dim == Dim::_3D) { + g_mesh.set_prtl_bc({ -1, 0, 0 }, new_bcs); + } else { + raise::Error("Invalid dimension", HERE); + } + } else if (dir == bc_in::Px2) { + if constexpr (M::Dim == Dim::_1D) { + raise::Error("Cannot set +x2 BCs for 1D", HERE); + } else if constexpr (M::Dim == Dim::_2D) { + g_mesh.set_prtl_bc({ +1, 0 }, new_bcs); + } else if constexpr (M::Dim == Dim::_3D) { + g_mesh.set_prtl_bc({ +1, 0, 0 }, new_bcs); + } else { + raise::Error("Invalid dimension", HERE); + } + } else if (dir == bc_in::Mx3) { + if constexpr (M::Dim == Dim::_1D) { + raise::Error("Cannot set -x3 BCs for 1D", HERE); + } else if constexpr (M::Dim == Dim::_2D) { + raise::Error("Cannot set -x3 BCs for 2D", HERE); + } else if constexpr (M::Dim == Dim::_3D) { + g_mesh.set_prtl_bc({ 0, 0, -1 }, new_bcs); + } else { + raise::Error("Invalid dimension", HERE); + } + } else if (dir == bc_in::Px3) { + if constexpr (M::Dim == Dim::_1D) { + raise::Error("Cannot set +x3 BCs for 1D", HERE); + } else if constexpr (M::Dim == Dim::_2D) { + raise::Error("Cannot set +x3 BCs for 2D", HERE); + } else if constexpr (M::Dim == Dim::_3D) { + g_mesh.set_prtl_bc({ 0, 0, +1 }, new_bcs); + } else { + raise::Error("Invalid dimension", HERE); + } + } else { + raise::Error("Invalid direction", HERE); + } + redefineBoundaries(); + } + template struct Metadomain>; template struct Metadomain>; template struct Metadomain>; diff --git a/src/framework/domain/metadomain.h b/src/framework/domain/metadomain.h index 1728fbd65..1fb6ce007 100644 --- a/src/framework/domain/metadomain.h +++ b/src/framework/domain/metadomain.h @@ -142,6 +142,8 @@ namespace ntt { -> bool; /* setters -------------------------------------------------------------- */ + void setFldsBC(const bc_in&, const FldsBC&); + void setPrtlBC(const bc_in&, const PrtlBC&); /* getters -------------------------------------------------------------- */ [[nodiscard]] From e0a7b18ee9a66ce56394fc01a07189bd58b1b8b8 Mon Sep 17 00:00:00 2001 From: hayk Date: Mon, 21 Apr 2025 06:35:19 -0400 Subject: [PATCH 09/25] rec setup BC + tested with MPI --- setups/wip/reconnection/pgen.hpp | 71 ++++++++++++++++++++------------ 1 file changed, 45 insertions(+), 26 deletions(-) diff --git a/setups/wip/reconnection/pgen.hpp b/setups/wip/reconnection/pgen.hpp index 70bab7b32..e92de5847 100644 --- a/setups/wip/reconnection/pgen.hpp +++ b/setups/wip/reconnection/pgen.hpp @@ -22,17 +22,19 @@ namespace user { template struct CurrentLayer : public arch::SpatialDistribution { - CurrentLayer(const M& metric, real_t cs_width, real_t cs_y) + CurrentLayer(const M& metric, real_t cs_width, real_t center_x, real_t cs_y) : arch::SpatialDistribution { metric } , cs_width { cs_width } + , center_x { center_x } , cs_y { cs_y } {} Inline auto operator()(const coord_t& x_Ph) const -> real_t override { - return ONE / SQR(math::cosh((x_Ph[1] - cs_y) / cs_width)); + return ONE / SQR(math::cosh((x_Ph[1] - cs_y) / cs_width)) * + (ONE - math::exp(-SQR((x_Ph[0] - center_x) / cs_width))); } private: - const real_t cs_y, cs_width; + const real_t cs_width, center_x, cs_y; }; // field initializer @@ -159,13 +161,17 @@ namespace user { using arch::ProblemGenerator::C; using arch::ProblemGenerator::params; - const real_t bg_B, bg_Bguide, bg_temperature, inj_ypad; - const real_t cs_width, cs_overdensity, cs_x, cs_y; - const real_t ymin, ymax; + const real_t bg_B, bg_Bguide, bg_temperature, inj_ypad; + const real_t cs_width, cs_overdensity, cs_x, cs_y; + const real_t ymin, ymax; + const simtime_t t_open; + bool bc_opened { false }; + + Metadomain& metadomain; InitFields init_flds; - inline PGen(const SimulationParams& p, const Metadomain& m) + inline PGen(const SimulationParams& p, Metadomain& m) : arch::ProblemGenerator(p) , bg_B { p.template get("setup.bg_B", 1.0) } , bg_Bguide { p.template get("setup.bg_Bguide", 0.0) } @@ -173,14 +179,17 @@ namespace user { , inj_ypad { p.template get("setup.inj_ypad", (real_t)0.05) } , cs_width { p.template get("setup.cs_width") } , cs_overdensity { p.template get("setup.cs_overdensity") } - , cs_x { m.mesh().extent(in::x1).first + - INV_2 * (m.mesh().extent(in::x1).second - - m.mesh().extent(in::x1).first) } - , cs_y { m.mesh().extent(in::x2).first + - INV_2 * (m.mesh().extent(in::x2).second - - m.mesh().extent(in::x2).first) } + , cs_x { INV_2 * + (m.mesh().extent(in::x1).second + m.mesh().extent(in::x1).first) } + , cs_y { INV_2 * + (m.mesh().extent(in::x2).second + m.mesh().extent(in::x2).first) } , ymin { m.mesh().extent(in::x2).first } , ymax { m.mesh().extent(in::x2).second } + , t_open { p.template get( + "setup.t_open", + 1.5 * HALF * + (m.mesh().extent(in::x1).second - m.mesh().extent(in::x1).first)) } + , metadomain { m } , init_flds { bg_B, bg_Bguide, cs_width, cs_y } {} inline PGen() {} @@ -225,6 +234,7 @@ namespace user { false); const auto sdist_cs = CurrentLayer(local_domain.mesh.metric, cs_width, + cs_x, cs_y); const auto inj_cs = arch::NonUniformInjector( edist_cs, @@ -236,7 +246,16 @@ namespace user { cs_overdensity); } - void CustomPostStep(timestep_t, simtime_t, Domain& domain) { + void CustomPostStep(timestep_t, simtime_t time, Domain& domain) { + // open boundaries if not yet opened at time = t_open + if ((t_open > 0.0) and (not bc_opened) and (time > t_open)) { + bc_opened = true; + metadomain.setFldsBC(bc_in::Mx1, FldsBC::MATCH); + metadomain.setPrtlBC(bc_in::Mx1, PrtlBC::ABSORB); + metadomain.setFldsBC(bc_in::Px1, FldsBC::MATCH); + metadomain.setPrtlBC(bc_in::Px1, PrtlBC::ABSORB); + } + const auto energy_dist = arch::Maxwellian(domain.mesh.metric, domain.random_pool, bg_temperature); @@ -272,18 +291,18 @@ namespace user { 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)); + 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); From ec5c65ce51a799689ebe79e55024d6cdfac50fb1 Mon Sep 17 00:00:00 2001 From: hayk Date: Mon, 21 Apr 2025 06:35:37 -0400 Subject: [PATCH 10/25] inlined writer --- src/output/stats.cpp | 31 ------------------------------- src/output/stats.h | 16 +++++++++++++++- 2 files changed, 15 insertions(+), 32 deletions(-) diff --git a/src/output/stats.cpp b/src/output/stats.cpp index 2dfac5c9b..ba3e85c40 100644 --- a/src/output/stats.cpp +++ b/src/output/stats.cpp @@ -11,8 +11,6 @@ #include -#include -#include #include #include @@ -90,21 +88,6 @@ namespace stats { m_stat_writers); } - template - void Writer::write(const T& value) { -#if defined(MPI_ENABLED) - // @TODO: reduce -#endif - CallOnce( - [](auto& fname, auto& value) { - std::fstream StatsOut(fname, std::fstream::out | std::fstream::app); - StatsOut << value << ","; - StatsOut.close(); - }, - m_fname, - value); - } - void Writer::endWriting() { CallOnce( [](auto& fname) { @@ -115,18 +98,4 @@ namespace stats { m_fname); } - template void Writer::write(const char&); - template void Writer::write(const std::string&); - template void Writer::write(const float&); - template void Writer::write(const double&); - template void Writer::write(const int&); - template void Writer::write(const unsigned int&); - template void Writer::write(const short&); - template void Writer::write(const unsigned short&); - template void Writer::write(const long int&); - template void Writer::write(const unsigned long int&); - template void Writer::write(const long long int&); - template void Writer::write( - const unsigned long long int&); - } // namespace stats diff --git a/src/output/stats.h b/src/output/stats.h index dbc874bda..c81b9b7b3 100644 --- a/src/output/stats.h +++ b/src/output/stats.h @@ -18,6 +18,8 @@ #include "utils/error.h" #include "utils/tools.h" +#include +#include #include #include @@ -141,7 +143,19 @@ namespace stats { auto shouldWrite(timestep_t, simtime_t) -> bool; template - void write(const T&); + inline void write(const T& value) const { +#if defined(MPI_ENABLED) + // @TODO: reduce +#endif + CallOnce( + [](auto& fname, auto& value) { + std::fstream StatsOut(fname, std::fstream::out | std::fstream::app); + StatsOut << value << ","; + StatsOut.close(); + }, + m_fname, + value); + } void endWriting(); From 1f8db6ca0c3d453b5ff60909cccec701b0465fbe Mon Sep 17 00:00:00 2001 From: hayk Date: Mon, 21 Apr 2025 06:35:48 -0400 Subject: [PATCH 11/25] minor print rm --- cmake/report.cmake | 1 - 1 file changed, 1 deletion(-) diff --git a/cmake/report.cmake b/cmake/report.cmake index 5a38b0dd5..626de15a2 100644 --- a/cmake/report.cmake +++ b/cmake/report.cmake @@ -155,7 +155,6 @@ if(${Kokkos_DEVICES} MATCHES "CUDA") COMMAND bash -c ${cmd} OUTPUT_VARIABLE CUDACOMP_VERSION OUTPUT_STRIP_TRAILING_WHITESPACE) - message(STATUS "CUDACOMP: ${CUDACOMP_VERSION}") string( APPEND REPORT_TEXT From f82b681c571c33c41229b30839f0101b099824d5 Mon Sep 17 00:00:00 2001 From: haykh Date: Wed, 9 Apr 2025 07:43:26 -0400 Subject: [PATCH 12/25] new deposit scheme --- src/framework/parameters.cpp | 8 +- src/kernels/currents_deposit.hpp | 460 +++++++++++++++---------------- 2 files changed, 219 insertions(+), 249 deletions(-) diff --git a/src/framework/parameters.cpp b/src/framework/parameters.cpp index 26f542956..8a30f52d9 100644 --- a/src/framework/parameters.cpp +++ b/src/framework/parameters.cpp @@ -31,10 +31,10 @@ namespace ntt { template - auto get_dx0_V0(const std::vector& resolution, - const boundaries_t& extent, - const std::map& params) - -> std::pair { + auto get_dx0_V0( + const std::vector& resolution, + const boundaries_t& extent, + const std::map& params) -> std::pair { const auto metric = M(resolution, extent, params); const auto dx0 = metric.dxMin(); coord_t x_corner { ZERO }; diff --git a/src/kernels/currents_deposit.hpp b/src/kernels/currents_deposit.hpp index ca9a94878..98d00a9b0 100644 --- a/src/kernels/currents_deposit.hpp +++ b/src/kernels/currents_deposit.hpp @@ -67,8 +67,8 @@ namespace kernel { const array_t& weight, const array_t& tag, const M& metric, - const real_t& charge, - const real_t& dt) + real_t charge, + real_t dt) : J { scatter_cur } , i1 { i1 } , i2 { i2 } @@ -100,45 +100,70 @@ namespace kernel { if (tag(p) == ParticleTag::dead) { return; } - // _f = final, _i = initial - tuple_t Ip_f, Ip_i; - coord_t xp_f, xp_i, xp_r; + // recover particle velocity to deposit in unsimulated direction vec_t vp { ZERO }; + { + coord_t xp { ZERO }; + if constexpr (D == Dim::_1D) { + xp[0] = i_di_to_Xi(i1(p), dx1(p)); + } else if constexpr (D == Dim::_2D) { + if constexpr (M::PrtlDim == Dim::_3D) { + xp[0] = i_di_to_Xi(i1(p), dx1(p)); + xp[1] = i_di_to_Xi(i2(p), dx2(p)); + xp[2] = phi(p); + } else { + xp[0] = i_di_to_Xi(i1(p), dx1(p)); + xp[1] = i_di_to_Xi(i2(p), dx2(p)); + } + } else { + xp[0] = i_di_to_Xi(i1(p), dx1(p)); + xp[1] = i_di_to_Xi(i2(p), dx2(p)); + xp[2] = i_di_to_Xi(i3(p), dx3(p)); + } + auto inv_energy { ZERO }; + if constexpr (S == SimEngine::SRPIC) { + metric.template transform_xyz(xp, + { ux1(p), ux2(p), ux3(p) }, + vp); + inv_energy = ONE / math::sqrt(ONE + NORM_SQR(ux1(p), ux2(p), ux3(p))); + } else { + metric.template transform(xp, + { ux1(p), ux2(p), ux3(p) }, + vp); + inv_energy = ONE / math::sqrt(ONE + ux1(p) * vp[0] + ux2(p) * vp[1] + + ux3(p) * vp[2]); + } + if (Kokkos::isnan(vp[2]) || Kokkos::isinf(vp[2])) { + vp[2] = ZERO; + } + vp[0] *= inv_energy; + vp[1] *= inv_energy; + vp[2] *= inv_energy; + } - // get [i, di]_init and [i, di]_final (per dimension) - getDepositInterval(p, Ip_f, Ip_i, xp_f, xp_i, xp_r); - // recover particle velocity to deposit in unsimulated direction - getPrtl3Vel(p, vp); const real_t coeff { weight(p) * charge }; - depositCurrentsFromParticle(coeff, vp, Ip_f, Ip_i, xp_f, xp_i, xp_r); - } - /** - * @brief Deposit currents from a single particle. - * @param[in] coeff Particle weight x charge. - * @param[in] vp Particle 3-velocity. - * @param[in] Ip_f Final position of the particle (cell index). - * @param[in] Ip_i Initial position of the particle (cell index). - * @param[in] xp_f Final position. - * @param[in] xp_i Previous step position. - * @param[in] xp_r Intermediate point used in zig-zag deposit. - */ - Inline auto depositCurrentsFromParticle(const real_t& coeff, - const vec_t& vp, - const tuple_t& Ip_f, - const tuple_t& Ip_i, - const coord_t& xp_f, - const coord_t& xp_i, - const coord_t& xp_r) const -> void { - const real_t Wx1_1 { HALF * (xp_i[0] + xp_r[0]) - - static_cast(Ip_i[0]) }; - const real_t Wx1_2 { HALF * (xp_f[0] + xp_r[0]) - - static_cast(Ip_f[0]) }; - const real_t Fx1_1 { (xp_r[0] - xp_i[0]) * coeff * inv_dt }; - const real_t Fx1_2 { (xp_f[0] - xp_r[0]) * coeff * inv_dt }; + const auto dxp_r_1 { static_cast(i1(p) == i1_prev(p)) * + (dx1(p) + dx1_prev(p)) * static_cast(INV_2) }; + + const real_t Wx1_1 { INV_2 * (dxp_r_1 + dx1_prev(p) + + static_cast(i1(p) > i1_prev(p))) }; + const real_t Wx1_2 { INV_2 * (dx1(p) + dxp_r_1 + + static_cast( + static_cast(i1(p) > i1_prev(p)) + + i1_prev(p) - i1(p))) }; + const real_t Fx1_1 { (static_cast(i1(p) > i1_prev(p)) + dxp_r_1 - + dx1_prev(p)) * + coeff * inv_dt }; + const real_t Fx1_2 { (static_cast( + i1(p) - i1_prev(p) - + static_cast(i1(p) > i1_prev(p))) + + dx1(p) - dxp_r_1) * + coeff * inv_dt }; auto J_acc = J.access(); + // tuple_t dxp_r; if constexpr (D == Dim::_1D) { const real_t Fx2_1 { HALF * vp[1] * coeff }; const real_t Fx2_2 { HALF * vp[1] * coeff }; @@ -146,265 +171,210 @@ namespace kernel { const real_t Fx3_1 { HALF * vp[2] * coeff }; const real_t Fx3_2 { HALF * vp[2] * coeff }; - J_acc(Ip_i[0] + N_GHOSTS, cur::jx1) += Fx1_1; - J_acc(Ip_f[0] + N_GHOSTS, cur::jx1) += Fx1_2; + J_acc(i1_prev(p) + N_GHOSTS, cur::jx1) += Fx1_1; + J_acc(i1(p) + N_GHOSTS, cur::jx1) += Fx1_2; - J_acc(Ip_i[0] + N_GHOSTS, cur::jx2) += Fx2_1 * (ONE - Wx1_1); - J_acc(Ip_i[0] + N_GHOSTS + 1, cur::jx2) += Fx2_1 * Wx1_1; - J_acc(Ip_f[0] + N_GHOSTS, cur::jx2) += Fx2_2 * (ONE - Wx1_2); - J_acc(Ip_f[0] + N_GHOSTS + 1, cur::jx2) += Fx2_2 * Wx1_2; + J_acc(i1_prev(p) + N_GHOSTS, cur::jx2) += Fx2_1 * (ONE - Wx1_1); + J_acc(i1_prev(p) + N_GHOSTS + 1, cur::jx2) += Fx2_1 * Wx1_1; + J_acc(i1(p) + N_GHOSTS, cur::jx2) += Fx2_2 * (ONE - Wx1_2); + J_acc(i1(p) + N_GHOSTS + 1, cur::jx2) += Fx2_2 * Wx1_2; - J_acc(Ip_i[0] + N_GHOSTS, cur::jx3) += Fx3_1 * (ONE - Wx1_1); - J_acc(Ip_i[0] + N_GHOSTS + 1, cur::jx3) += Fx3_1 * Wx1_1; - J_acc(Ip_f[0] + N_GHOSTS, cur::jx3) += Fx3_2 * (ONE - Wx1_2); - J_acc(Ip_f[0] + N_GHOSTS + 1, cur::jx3) += Fx3_2 * Wx1_2; + J_acc(i1_prev(p) + N_GHOSTS, cur::jx3) += Fx3_1 * (ONE - Wx1_1); + J_acc(i1_prev(p) + N_GHOSTS + 1, cur::jx3) += Fx3_1 * Wx1_1; + J_acc(i1(p) + N_GHOSTS, cur::jx3) += Fx3_2 * (ONE - Wx1_2); + J_acc(i1(p) + N_GHOSTS + 1, cur::jx3) += Fx3_2 * Wx1_2; } else if constexpr (D == Dim::_2D || D == Dim::_3D) { - const real_t Wx2_1 { HALF * (xp_i[1] + xp_r[1]) - - static_cast(Ip_i[1]) }; - const real_t Wx2_2 { HALF * (xp_f[1] + xp_r[1]) - - static_cast(Ip_f[1]) }; - const real_t Fx2_1 { (xp_r[1] - xp_i[1]) * coeff * inv_dt }; - const real_t Fx2_2 { (xp_f[1] - xp_r[1]) * coeff * inv_dt }; + const auto dxp_r_2 { static_cast(i2(p) == i2_prev(p)) * + (dx2(p) + dx2_prev(p)) * + static_cast(INV_2) }; + + const real_t Wx2_1 { INV_2 * (dxp_r_2 + dx2_prev(p) + + static_cast(i2(p) > i2_prev(p))) }; + const real_t Wx2_2 { INV_2 * (dx2(p) + dxp_r_2 + + static_cast( + static_cast(i2(p) > i2_prev(p)) + + i2_prev(p) - i2(p))) }; + const real_t Fx2_1 { (static_cast(i2(p) > i2_prev(p)) + + dxp_r_2 - dx2_prev(p)) * + coeff * inv_dt }; + const real_t Fx2_2 { (static_cast( + i2(p) - i2_prev(p) - + static_cast(i2(p) > i2_prev(p))) + + dx2(p) - dxp_r_2) * + coeff * inv_dt }; if constexpr (D == Dim::_2D) { const real_t Fx3_1 { HALF * vp[2] * coeff }; const real_t Fx3_2 { HALF * vp[2] * coeff }; - J_acc(Ip_i[0] + N_GHOSTS, Ip_i[1] + N_GHOSTS, cur::jx1) += Fx1_1 * - (ONE - Wx2_1); - J_acc(Ip_i[0] + N_GHOSTS, Ip_i[1] + N_GHOSTS + 1, cur::jx1) += Fx1_1 * - Wx2_1; - J_acc(Ip_f[0] + N_GHOSTS, Ip_f[1] + N_GHOSTS, cur::jx1) += Fx1_2 * - (ONE - Wx2_2); - J_acc(Ip_f[0] + N_GHOSTS, Ip_f[1] + N_GHOSTS + 1, cur::jx1) += Fx1_2 * - Wx2_2; - - J_acc(Ip_i[0] + N_GHOSTS, Ip_i[1] + N_GHOSTS, cur::jx2) += Fx2_1 * - (ONE - Wx1_1); - J_acc(Ip_i[0] + N_GHOSTS + 1, Ip_i[1] + N_GHOSTS, cur::jx2) += Fx2_1 * - Wx1_1; - J_acc(Ip_f[0] + N_GHOSTS, Ip_f[1] + N_GHOSTS, cur::jx2) += Fx2_2 * - (ONE - Wx1_2); - J_acc(Ip_f[0] + N_GHOSTS + 1, Ip_f[1] + N_GHOSTS, cur::jx2) += Fx2_2 * - Wx1_2; - - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS, + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS, + cur::jx1) += Fx1_1 * (ONE - Wx2_1); + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS + 1, + cur::jx1) += Fx1_1 * Wx2_1; + J_acc(i1(p) + N_GHOSTS, i2(p) + N_GHOSTS, cur::jx1) += Fx1_2 * + (ONE - Wx2_2); + J_acc(i1(p) + N_GHOSTS, i2(p) + N_GHOSTS + 1, cur::jx1) += Fx1_2 * Wx2_2; + + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS, + cur::jx2) += Fx2_1 * (ONE - Wx1_1); + J_acc(i1_prev(p) + N_GHOSTS + 1, + i2_prev(p) + N_GHOSTS, + cur::jx2) += Fx2_1 * Wx1_1; + J_acc(i1(p) + N_GHOSTS, i2(p) + N_GHOSTS, cur::jx2) += Fx2_2 * + (ONE - Wx1_2); + J_acc(i1(p) + N_GHOSTS + 1, i2(p) + N_GHOSTS, cur::jx2) += Fx2_2 * Wx1_2; + + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS, cur::jx3) += Fx3_1 * (ONE - Wx1_1) * (ONE - Wx2_1); - J_acc(Ip_i[0] + N_GHOSTS + 1, - Ip_i[1] + N_GHOSTS, + J_acc(i1_prev(p) + N_GHOSTS + 1, + i2_prev(p) + N_GHOSTS, cur::jx3) += Fx3_1 * Wx1_2 * (ONE - Wx2_1); - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS + 1, + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS + 1, cur::jx3) += Fx3_1 * (ONE - Wx1_1) * Wx2_1; - J_acc(Ip_i[0] + N_GHOSTS + 1, - Ip_i[1] + N_GHOSTS + 1, + J_acc(i1_prev(p) + N_GHOSTS + 1, + i2_prev(p) + N_GHOSTS + 1, cur::jx3) += Fx3_1 * Wx1_1 * Wx2_1; - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS, - cur::jx3) += Fx3_2 * (ONE - Wx1_2) * (ONE - Wx2_2); - J_acc(Ip_f[0] + N_GHOSTS + 1, - Ip_f[1] + N_GHOSTS, + J_acc(i1(p) + N_GHOSTS, i2(p) + N_GHOSTS, cur::jx3) += Fx3_2 * + (ONE - Wx1_2) * + (ONE - Wx2_2); + J_acc(i1(p) + N_GHOSTS + 1, + i2(p) + N_GHOSTS, cur::jx3) += Fx3_2 * Wx1_2 * (ONE - Wx2_2); - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS + 1, + J_acc(i1(p) + N_GHOSTS, + i2(p) + N_GHOSTS + 1, cur::jx3) += Fx3_2 * (ONE - Wx1_2) * Wx2_2; - J_acc(Ip_f[0] + N_GHOSTS + 1, - Ip_f[1] + N_GHOSTS + 1, - cur::jx3) += Fx3_2 * Wx1_2 * Wx2_2; + J_acc(i1(p) + N_GHOSTS + 1, i2(p) + N_GHOSTS + 1, cur::jx3) += Fx3_2 * + Wx1_2 * + Wx2_2; } else { - const real_t Wx3_1 { HALF * (xp_i[2] + xp_r[2]) - - static_cast(Ip_i[2]) }; - const real_t Wx3_2 { HALF * (xp_f[2] + xp_r[2]) - - static_cast(Ip_f[2]) }; - const real_t Fx3_1 { (xp_r[2] - xp_i[2]) * coeff * inv_dt }; - const real_t Fx3_2 { (xp_f[2] - xp_r[2]) * coeff * inv_dt }; - - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS, - Ip_i[2] + N_GHOSTS, + const auto dxp_r_3 { static_cast(i3(p) == i3_prev(p)) * + (dx3(p) + dx3_prev(p)) * + static_cast(INV_2) }; + const real_t Wx3_1 { INV_2 * (dxp_r_3 + dx3_prev(p) + + static_cast(i3(p) > i3_prev(p))) }; + const real_t Wx3_2 { INV_2 * (dx3(p) + dxp_r_3 + + static_cast( + static_cast(i3(p) > i3_prev(p)) + + i3_prev(p) - i3(p))) }; + const real_t Fx3_1 { (static_cast(i3(p) > i3_prev(p)) + + dxp_r_3 - dx3_prev(p)) * + coeff * inv_dt }; + const real_t Fx3_2 { (static_cast( + i3(p) - i3_prev(p) - + static_cast(i3(p) > i3_prev(p))) + + dx3(p) - dxp_r_3) * + coeff * inv_dt }; + + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS, + i3_prev(p) + N_GHOSTS, cur::jx1) += Fx1_1 * (ONE - Wx2_1) * (ONE - Wx3_1); - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS + 1, - Ip_i[2] + N_GHOSTS, + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS + 1, + i3_prev(p) + N_GHOSTS, cur::jx1) += Fx1_1 * Wx2_1 * (ONE - Wx3_1); - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS, - Ip_i[2] + N_GHOSTS + 1, + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS, + i3_prev(p) + N_GHOSTS + 1, cur::jx1) += Fx1_1 * (ONE - Wx2_1) * Wx3_1; - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS + 1, - Ip_i[2] + N_GHOSTS + 1, + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS + 1, + i3_prev(p) + N_GHOSTS + 1, cur::jx1) += Fx1_1 * Wx2_1 * Wx3_1; - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS, - Ip_f[2] + N_GHOSTS, + J_acc(i1(p) + N_GHOSTS, + i2(p) + N_GHOSTS, + i3(p) + N_GHOSTS, cur::jx1) += Fx1_2 * (ONE - Wx2_2) * (ONE - Wx3_2); - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS + 1, - Ip_f[2] + N_GHOSTS, + J_acc(i1(p) + N_GHOSTS, + i2(p) + N_GHOSTS + 1, + i3(p) + N_GHOSTS, cur::jx1) += Fx1_2 * Wx2_2 * (ONE - Wx3_2); - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS, - Ip_f[2] + N_GHOSTS + 1, + J_acc(i1(p) + N_GHOSTS, + i2(p) + N_GHOSTS, + i3(p) + N_GHOSTS + 1, cur::jx1) += Fx1_2 * (ONE - Wx2_2) * Wx3_2; - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS + 1, - Ip_f[2] + N_GHOSTS + 1, + J_acc(i1(p) + N_GHOSTS, + i2(p) + N_GHOSTS + 1, + i3(p) + N_GHOSTS + 1, cur::jx1) += Fx1_2 * Wx2_2 * Wx3_2; - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS, - Ip_i[2] + N_GHOSTS, + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS, + i3_prev(p) + N_GHOSTS, cur::jx2) += Fx2_1 * (ONE - Wx1_1) * (ONE - Wx3_1); - J_acc(Ip_i[0] + N_GHOSTS + 1, - Ip_i[1] + N_GHOSTS, - Ip_i[2] + N_GHOSTS, + J_acc(i1_prev(p) + N_GHOSTS + 1, + i2_prev(p) + N_GHOSTS, + i3_prev(p) + N_GHOSTS, cur::jx2) += Fx2_1 * Wx1_1 * (ONE - Wx3_1); - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS, - Ip_i[2] + N_GHOSTS + 1, + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS, + i3_prev(p) + N_GHOSTS + 1, cur::jx2) += Fx2_1 * (ONE - Wx1_1) * Wx3_1; - J_acc(Ip_i[0] + N_GHOSTS + 1, - Ip_i[1] + N_GHOSTS, - Ip_i[2] + N_GHOSTS + 1, + J_acc(i1_prev(p) + N_GHOSTS + 1, + i2_prev(p) + N_GHOSTS, + i3_prev(p) + N_GHOSTS + 1, cur::jx2) += Fx2_1 * Wx1_1 * Wx3_1; - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS, - Ip_f[2] + N_GHOSTS, + J_acc(i1(p) + N_GHOSTS, + i2(p) + N_GHOSTS, + i3(p) + N_GHOSTS, cur::jx2) += Fx2_2 * (ONE - Wx1_2) * (ONE - Wx3_2); - J_acc(Ip_f[0] + N_GHOSTS + 1, - Ip_f[1] + N_GHOSTS, - Ip_f[2] + N_GHOSTS, + J_acc(i1(p) + N_GHOSTS + 1, + i2(p) + N_GHOSTS, + i3(p) + N_GHOSTS, cur::jx2) += Fx2_2 * Wx1_2 * (ONE - Wx3_2); - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS, - Ip_f[2] + N_GHOSTS + 1, + J_acc(i1(p) + N_GHOSTS, + i2(p) + N_GHOSTS, + i3(p) + N_GHOSTS + 1, cur::jx2) += Fx2_2 * (ONE - Wx1_2) * Wx3_2; - J_acc(Ip_f[0] + N_GHOSTS + 1, - Ip_f[1] + N_GHOSTS, - Ip_f[2] + N_GHOSTS + 1, + J_acc(i1(p) + N_GHOSTS + 1, + i2(p) + N_GHOSTS, + i3(p) + N_GHOSTS + 1, cur::jx2) += Fx2_2 * Wx1_2 * Wx3_2; - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS, - Ip_i[2] + N_GHOSTS, + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS, + i3_prev(p) + N_GHOSTS, cur::jx3) += Fx3_1 * (ONE - Wx1_1) * (ONE - Wx2_1); - J_acc(Ip_i[0] + N_GHOSTS + 1, - Ip_i[1] + N_GHOSTS, - Ip_i[2] + N_GHOSTS, + J_acc(i1_prev(p) + N_GHOSTS + 1, + i2_prev(p) + N_GHOSTS, + i3_prev(p) + N_GHOSTS, cur::jx3) += Fx3_1 * Wx1_1 * (ONE - Wx2_1); - J_acc(Ip_i[0] + N_GHOSTS, - Ip_i[1] + N_GHOSTS + 1, - Ip_i[2] + N_GHOSTS, + J_acc(i1_prev(p) + N_GHOSTS, + i2_prev(p) + N_GHOSTS + 1, + i3_prev(p) + N_GHOSTS, cur::jx3) += Fx3_1 * (ONE - Wx1_1) * Wx2_1; - J_acc(Ip_i[0] + N_GHOSTS + 1, - Ip_i[1] + N_GHOSTS + 1, - Ip_i[2] + N_GHOSTS, + J_acc(i1_prev(p) + N_GHOSTS + 1, + i2_prev(p) + N_GHOSTS + 1, + i3_prev(p) + N_GHOSTS, cur::jx3) += Fx3_1 * Wx1_1 * Wx2_1; - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS, - Ip_f[2] + N_GHOSTS, + J_acc(i1(p) + N_GHOSTS, + i2(p) + N_GHOSTS, + i3(p) + N_GHOSTS, cur::jx3) += Fx3_2 * (ONE - Wx1_2) * (ONE - Wx2_2); - J_acc(Ip_f[0] + N_GHOSTS + 1, - Ip_f[1] + N_GHOSTS, - Ip_f[2] + N_GHOSTS, + J_acc(i1(p) + N_GHOSTS + 1, + i2(p) + N_GHOSTS, + i3(p) + N_GHOSTS, cur::jx3) += Fx3_2 * Wx1_2 * (ONE - Wx2_2); - J_acc(Ip_f[0] + N_GHOSTS, - Ip_f[1] + N_GHOSTS + 1, - Ip_f[2] + N_GHOSTS, + J_acc(i1(p) + N_GHOSTS, + i2(p) + N_GHOSTS + 1, + i3(p) + N_GHOSTS, cur::jx3) += Fx3_2 * (ONE - Wx1_2) * Wx2_2; - J_acc(Ip_f[0] + N_GHOSTS + 1, - Ip_f[1] + N_GHOSTS + 1, - Ip_f[2] + N_GHOSTS, + J_acc(i1(p) + N_GHOSTS + 1, + i2(p) + N_GHOSTS + 1, + i3(p) + N_GHOSTS, cur::jx3) += Fx3_2 * Wx1_2 * Wx2_2; } } } - - /** - * @brief Get particle position in `coord_t` form. - * @param[in] p Index of particle. - * @param[out] Ip_f Final position of the particle (cell index). - * @param[out] Ip_i Initial position of the particle (cell index). - * @param[out] xp_f Final position. - * @param[out] xp_i Previous step position. - * @param[out] xp_r Intermediate point used in zig-zag deposit. - */ - Inline auto getDepositInterval(index_t& p, - tuple_t& Ip_f, - tuple_t& Ip_i, - coord_t& xp_f, - coord_t& xp_i, - coord_t& xp_r) const -> void { - Ip_f[0] = i1(p); - Ip_i[0] = i1_prev(p); - xp_f[0] = i_di_to_Xi(Ip_f[0], dx1(p)); - xp_i[0] = i_di_to_Xi(Ip_i[0], dx1_prev(p)); - if constexpr (D == Dim::_2D || D == Dim::_3D) { - Ip_f[1] = i2(p); - Ip_i[1] = i2_prev(p); - xp_f[1] = i_di_to_Xi(Ip_f[1], dx2(p)); - xp_i[1] = i_di_to_Xi(Ip_i[1], dx2_prev(p)); - } - if constexpr (D == Dim::_3D) { - Ip_f[2] = i3(p); - Ip_i[2] = i3_prev(p); - xp_f[2] = i_di_to_Xi(Ip_f[2], dx3(p)); - xp_i[2] = i_di_to_Xi(Ip_i[2], dx3_prev(p)); - } - for (auto i = 0u; i < D; ++i) { - xp_r[i] = math::fmin(static_cast(IMIN(Ip_i[i], Ip_f[i]) + 1), - math::fmax(static_cast(IMAX(Ip_i[i], Ip_f[i])), - HALF * (xp_i[i] + xp_f[i]))); - } - } - - // Getters - Inline void getPrtlPos(index_t& p, coord_t& xp) const { - if constexpr (D == Dim::_1D) { - xp[0] = i_di_to_Xi(i1(p), dx1(p)); - } else if constexpr (D == Dim::_2D) { - if constexpr (M::PrtlDim == Dim::_3D) { - xp[0] = i_di_to_Xi(i1(p), dx1(p)); - xp[1] = i_di_to_Xi(i2(p), dx2(p)); - xp[2] = phi(p); - } else { - xp[0] = i_di_to_Xi(i1(p), dx1(p)); - xp[1] = i_di_to_Xi(i2(p), dx2(p)); - } - } else { - xp[0] = i_di_to_Xi(i1(p), dx1(p)); - xp[1] = i_di_to_Xi(i2(p), dx2(p)); - xp[2] = i_di_to_Xi(i3(p), dx3(p)); - } - } - - Inline void getPrtl3Vel(index_t& p, vec_t& vp) const { - coord_t xp { ZERO }; - getPrtlPos(p, xp); - auto inv_energy { ZERO }; - if constexpr (S == SimEngine::SRPIC) { - metric.template transform_xyz(xp, - { ux1(p), ux2(p), ux3(p) }, - vp); - inv_energy = ONE / math::sqrt(ONE + NORM_SQR(ux1(p), ux2(p), ux3(p))); - } else { - metric.template transform(xp, { ux1(p), ux2(p), ux3(p) }, vp); - inv_energy = ONE / math::sqrt(ONE + ux1(p) * vp[0] + ux2(p) * vp[1] + - ux3(p) * vp[2]); - } - if (Kokkos::isnan(vp[2]) || Kokkos::isinf(vp[2])) { - vp[2] = ZERO; - } - vp[0] *= inv_energy; - vp[1] *= inv_energy; - vp[2] *= inv_energy; - } }; } // namespace kernel From d66c2b5ef730fc5b78b1909308f1b49dbff8e541 Mon Sep 17 00:00:00 2001 From: hayk Date: Mon, 21 Apr 2025 11:48:21 -0400 Subject: [PATCH 13/25] (RUNTEST) From 9bc0f2c66067badccd73625e310fc844b764280f Mon Sep 17 00:00:00 2001 From: hayk Date: Mon, 21 Apr 2025 11:54:17 -0400 Subject: [PATCH 14/25] minor bug in enum test (RUNTEST) --- src/global/tests/enums.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/global/tests/enums.cpp b/src/global/tests/enums.cpp index 9208ed8d3..1190417ef 100644 --- a/src/global/tests/enums.cpp +++ b/src/global/tests/enums.cpp @@ -71,8 +71,8 @@ auto main() -> int { "h", "j", "a", "t", "rho", "charge", "n", "nppc", "v", "custom" }; - enum_str_t all_out_stats = { "b^2", "e^2", "exb", "j.e", "t", - "rho", "charge", "n", "npart", "v" }; + enum_str_t all_out_stats = { "b^2", "e^2", "exb", "j.e", "t", + "rho", "charge", "n", "npart" }; checkEnum(all_coords); checkEnum(all_metrics); From 9c58aa4eafdfaf039b7cb0a79a3cb2957d79bfe0 Mon Sep 17 00:00:00 2001 From: hayk Date: Mon, 21 Apr 2025 11:58:45 -0400 Subject: [PATCH 15/25] accuracy in deposit test (RUNTEST) --- src/kernels/tests/deposit.cpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/kernels/tests/deposit.cpp b/src/kernels/tests/deposit.cpp index ec364a313..63d65a619 100644 --- a/src/kernels/tests/deposit.cpp +++ b/src/kernels/tests/deposit.cpp @@ -29,7 +29,8 @@ void errorIf(bool condition, const std::string& message) { inline static constexpr auto epsilon = std::numeric_limits::epsilon(); -Inline auto equal(real_t a, real_t b, const char* msg = "", real_t acc = ONE) -> bool { +Inline auto equal(real_t a, real_t b, const char* msg = "", real_t acc = ONE) + -> bool { const auto eps = epsilon * acc; if (not cmp::AlmostEqual(a, b, eps)) { printf("%.12e != %.12e %s\n", a, b, msg); @@ -176,7 +177,7 @@ auto main(int argc, char* argv[]) -> int { }, { { 0.0, 55.0 }, { 0.0, 55.0 } }, {}, - 30); + 100); testDeposit, SimEngine::SRPIC>( { @@ -185,7 +186,7 @@ auto main(int argc, char* argv[]) -> int { }, { { 1.0, 100.0 } }, {}, - 30); + 100); testDeposit, SimEngine::SRPIC>( { @@ -194,7 +195,7 @@ auto main(int argc, char* argv[]) -> int { }, { { 1.0, 100.0 } }, { { "r0", 0.0 }, { "h", 0.25 } }, - 30); + 100); testDeposit, SimEngine::GRPIC>( { @@ -203,7 +204,7 @@ auto main(int argc, char* argv[]) -> int { }, { { 1.0, 100.0 } }, { { "a", 0.9 } }, - 30); + 100); testDeposit, SimEngine::GRPIC>( { @@ -212,7 +213,7 @@ auto main(int argc, char* argv[]) -> int { }, { { 1.0, 100.0 } }, { { "r0", 0.0 }, { "h", 0.25 }, { "a", 0.9 } }, - 30); + 100); testDeposit, SimEngine::GRPIC>( { @@ -221,7 +222,7 @@ auto main(int argc, char* argv[]) -> int { }, { { 1.0, 100.0 } }, { { "a", 0.9 } }, - 30); + 100); } catch (std::exception& e) { std::cerr << e.what() << std::endl; From 479425489a01054d13aadff6927028d990a3dc47 Mon Sep 17 00:00:00 2001 From: hayk Date: Mon, 21 Apr 2025 13:27:48 -0400 Subject: [PATCH 16/25] long double -> simtime_t --- src/checkpoint/writer.cpp | 2 +- src/checkpoint/writer.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/checkpoint/writer.cpp b/src/checkpoint/writer.cpp index fe77d3b56..b8571c246 100644 --- a/src/checkpoint/writer.cpp +++ b/src/checkpoint/writer.cpp @@ -200,7 +200,7 @@ namespace checkpoint { m_writer.Put(var, &data); } - void Writer::saveAttrs(const ntt::SimulationParams& params, long double time) { + void Writer::saveAttrs(const ntt::SimulationParams& params, simtime_t time) { CallOnce([&]() { std::ofstream metadata; if (m_written.empty()) { diff --git a/src/checkpoint/writer.h b/src/checkpoint/writer.h index 992c54c96..91c9d7a41 100644 --- a/src/checkpoint/writer.h +++ b/src/checkpoint/writer.h @@ -54,7 +54,7 @@ namespace checkpoint { void beginSaving(timestep_t, simtime_t); void endSaving(); - void saveAttrs(const ntt::SimulationParams&, long double); + void saveAttrs(const ntt::SimulationParams&, simtime_t); template void savePerDomainVariable(const std::string&, std::size_t, std::size_t, T); From e1ae5989e1b1015d4e307b0cf38c3d0c61358486 Mon Sep 17 00:00:00 2001 From: hayk Date: Mon, 21 Apr 2025 13:28:23 -0400 Subject: [PATCH 17/25] acc in dep test --- src/kernels/tests/deposit.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/kernels/tests/deposit.cpp b/src/kernels/tests/deposit.cpp index 63d65a619..e6967eb14 100644 --- a/src/kernels/tests/deposit.cpp +++ b/src/kernels/tests/deposit.cpp @@ -177,7 +177,7 @@ auto main(int argc, char* argv[]) -> int { }, { { 0.0, 55.0 }, { 0.0, 55.0 } }, {}, - 100); + 500); testDeposit, SimEngine::SRPIC>( { @@ -186,7 +186,7 @@ auto main(int argc, char* argv[]) -> int { }, { { 1.0, 100.0 } }, {}, - 100); + 500); testDeposit, SimEngine::SRPIC>( { @@ -195,7 +195,7 @@ auto main(int argc, char* argv[]) -> int { }, { { 1.0, 100.0 } }, { { "r0", 0.0 }, { "h", 0.25 } }, - 100); + 500); testDeposit, SimEngine::GRPIC>( { @@ -204,7 +204,7 @@ auto main(int argc, char* argv[]) -> int { }, { { 1.0, 100.0 } }, { { "a", 0.9 } }, - 100); + 500); testDeposit, SimEngine::GRPIC>( { @@ -213,7 +213,7 @@ auto main(int argc, char* argv[]) -> int { }, { { 1.0, 100.0 } }, { { "r0", 0.0 }, { "h", 0.25 }, { "a", 0.9 } }, - 100); + 500); testDeposit, SimEngine::GRPIC>( { @@ -222,7 +222,7 @@ auto main(int argc, char* argv[]) -> int { }, { { 1.0, 100.0 } }, { { "a", 0.9 } }, - 100); + 500); } catch (std::exception& e) { std::cerr << e.what() << std::endl; From ab854542327fc0f6670256307c1e336852fda843 Mon Sep 17 00:00:00 2001 From: hayk Date: Mon, 21 Apr 2025 13:28:32 -0400 Subject: [PATCH 18/25] (RUNTEST) From d5cdbd9656b3d114b4715a4e5f89d1a95fbeaa5c Mon Sep 17 00:00:00 2001 From: hayk Date: Tue, 22 Apr 2025 17:16:49 -0400 Subject: [PATCH 19/25] proper comm test --- src/framework/tests/comm_mpi.cpp | 237 +++++++++++++++++++++++++---- src/framework/tests/comm_nompi.cpp | 10 +- 2 files changed, 216 insertions(+), 31 deletions(-) diff --git a/src/framework/tests/comm_mpi.cpp b/src/framework/tests/comm_mpi.cpp index 5d2c8d4f0..487976f73 100644 --- a/src/framework/tests/comm_mpi.cpp +++ b/src/framework/tests/comm_mpi.cpp @@ -5,6 +5,11 @@ #include "arch/directions.h" #include "arch/kokkos_aliases.h" +#include "utils/error.h" +#include "utils/numeric.h" + +#include +#include #include #include @@ -13,49 +18,227 @@ using namespace ntt; auto main(int argc, char* argv[]) -> int { Kokkos::initialize(argc, argv); + MPI_Init(&argc, &argv); try { - const ncells_t nx1 = 15, nx2 = 15; - ndfield_t fld_b1 { "fld", nx1 + 2 * N_GHOSTS, nx2 + 2 * N_GHOSTS }; - ndfield_t fld_b2 { "fld", nx1 + 2 * N_GHOSTS, nx2 + 2 * N_GHOSTS }; + int rank, size; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &size); + + const ncells_t nx1 = 11, nx2 = 15; + ndfield_t fld { "fld", nx1 + 2 * N_GHOSTS, nx2 + 2 * N_GHOSTS }; Kokkos::parallel_for( "Fill", CreateRangePolicy({ 0, 0 }, { nx1 + 2 * N_GHOSTS, nx2 + 2 * N_GHOSTS }), Lambda(index_t i1, index_t i2) { - if ((i1 >= 2 * N_GHOSTS) and (i1 < nx1) and (i2 >= 2 * N_GHOSTS) and - (i2 < nx2)) { - fld_b1(i1, i2, 0) = 4.0; - fld_b1(i1, i2, 1) = 12.0; - fld_b1(i1, i2, 2) = 20.0; - fld_b2(i1, i2, 0) = 4.0; - fld_b2(i1, i2, 1) = 12.0; - fld_b2(i1, i2, 2) = 20.0; - } else if ( - ((i1 < 2 * N_GHOSTS or i1 >= nx1) and (i2 >= 2 * N_GHOSTS and i2 < nx2)) or - ((i2 < 2 * N_GHOSTS or i2 >= nx2) and (i1 >= 2 * N_GHOSTS and i1 < nx1))) { - fld_b1(i1, i2, 0) = 2.0; - fld_b1(i1, i2, 1) = 6.0; - fld_b1(i1, i2, 2) = 10.0; - fld_b2(i1, i2, 0) = 2.0; - fld_b2(i1, i2, 1) = 6.0; - fld_b2(i1, i2, 2) = 10.0; - } else { - fld_b1(i1, i2, 0) = 1.0; - fld_b1(i1, i2, 1) = 3.0; - fld_b1(i1, i2, 2) = 5.0; - fld_b2(i1, i2, 0) = 1.0; - fld_b2(i1, i2, 1) = 3.0; - fld_b2(i1, i2, 2) = 5.0; + if ((i1 >= N_GHOSTS) and (i1 < N_GHOSTS + nx1) and (i2 >= N_GHOSTS) and + (i2 < N_GHOSTS + nx2)) { + fld(i1, i2, 0) = static_cast(rank + 1) + 4.0; + fld(i1, i2, 1) = static_cast(rank + 1) + 12.0; + fld(i1, i2, 2) = static_cast(rank + 1) + 20.0; + } + }); + + { + // send right, recv left + const int send_idx = (rank + 1) % size; + const int recv_idx = (rank - 1 + size) % size; + const unsigned int send_rank = (unsigned int)send_idx; + const unsigned int recv_rank = (unsigned int)recv_idx; + + const std::vector send_slice { + { nx1, nx1 + N_GHOSTS }, + { N_GHOSTS, nx2 + N_GHOSTS } + }; + const std::vector recv_slice { + { 0, N_GHOSTS }, + { N_GHOSTS, nx2 + N_GHOSTS } + }; + const range_tuple_t comp_slice { 0, 3 }; + comm::CommunicateField((unsigned int)(rank), + fld, + fld, + send_idx, + recv_idx, + send_rank, + recv_rank, + send_slice, + recv_slice, + comp_slice, + false); + } + { + // recv right, send left + const int send_idx = (rank - 1 + size) % size; + const int recv_idx = (rank + 1) % size; + const unsigned int send_rank = (unsigned int)send_idx; + const unsigned int recv_rank = (unsigned int)recv_idx; + + const std::vector send_slice { + { N_GHOSTS, N_GHOSTS + 2 }, + { N_GHOSTS, nx2 + N_GHOSTS } + }; + const std::vector recv_slice { + { nx1 + N_GHOSTS, nx1 + 2 * N_GHOSTS }, + { N_GHOSTS, nx2 + N_GHOSTS } + }; + const range_tuple_t comp_slice { 0, 3 }; + comm::CommunicateField((unsigned int)(rank), + fld, + fld, + send_idx, + recv_idx, + send_rank, + recv_rank, + send_slice, + recv_slice, + comp_slice, + false); + } + + { + const auto left_expect = static_cast((rank - 1 + size) % size + 1); + const auto right_expect = static_cast((rank + 1) % size + 1); + + Kokkos::parallel_for( + "Check", + CreateRangePolicy({ N_GHOSTS }, { nx2 + N_GHOSTS }), + Lambda(index_t i2) { + for (auto i1 { 0u }; i1 < N_GHOSTS; ++i1) { + if (fld(i1, i2, 0) != left_expect + 4.0) { + raise::KernelError(HERE, "Left boundary not correct for #0"); + } + if (fld(i1, i2, 1) != left_expect + 12.0) { + raise::KernelError(HERE, "Left boundary not correct for #1"); + } + if (fld(i1, i2, 2) != left_expect + 20.0) { + raise::KernelError(HERE, "Left boundary not correct for #2"); + } + } + for (auto i1 { nx1 + N_GHOSTS }; i1 < nx1 + 2 * N_GHOSTS; ++i1) { + if (fld(i1, i2, 0) != right_expect + 4.0) { + raise::KernelError(HERE, "Right boundary not correct for #0"); + } + if (fld(i1, i2, 1) != right_expect + 12.0) { + raise::KernelError(HERE, "Right boundary not correct for #1"); + } + if (fld(i1, i2, 2) != right_expect + 20.0) { + raise::KernelError(HERE, "Right boundary not correct for #2"); + } + } + }); + } + + Kokkos::parallel_for( + "Carve", + CreateRangePolicy({ 0, 0 }, + { nx1 + 2 * N_GHOSTS, nx2 + 2 * N_GHOSTS }), + Lambda(index_t i1, index_t i2) { + if (((i1 >= N_GHOSTS) and (i1 < 2 * N_GHOSTS)) or + ((i1 >= nx1) and (i1 < nx1 + N_GHOSTS))) { + fld(i1, i2, 0) = ZERO; + fld(i1, i2, 1) = ZERO; + fld(i1, i2, 2) = ZERO; } }); + + { + // send right, recv left + const int send_idx = (rank + 1) % size; + const int recv_idx = (rank - 1 + size) % size; + const unsigned int send_rank = (unsigned int)send_idx; + const unsigned int recv_rank = (unsigned int)recv_idx; + + const std::vector send_slice { + { nx1 + N_GHOSTS, nx1 + 2 * N_GHOSTS }, + { N_GHOSTS, nx2 + N_GHOSTS } + }; + const std::vector recv_slice { + { N_GHOSTS, 2 * N_GHOSTS }, + { N_GHOSTS, nx2 + N_GHOSTS } + }; + const range_tuple_t comp_slice { 0, 3 }; + comm::CommunicateField((unsigned int)(rank), + fld, + fld, + send_idx, + recv_idx, + send_rank, + recv_rank, + send_slice, + recv_slice, + comp_slice, + true); + } + { + // recv right, send left + const int send_idx = (rank - 1 + size) % size; + const int recv_idx = (rank + 1) % size; + const unsigned int send_rank = (unsigned int)send_idx; + const unsigned int recv_rank = (unsigned int)recv_idx; + + const std::vector send_slice { + { 0, N_GHOSTS }, + { N_GHOSTS, nx2 + N_GHOSTS } + }; + const std::vector recv_slice { + { nx1, nx1 + N_GHOSTS }, + { N_GHOSTS, nx2 + N_GHOSTS } + }; + const range_tuple_t comp_slice { 0, 3 }; + comm::CommunicateField((unsigned int)(rank), + fld, + fld, + send_idx, + recv_idx, + send_rank, + recv_rank, + send_slice, + recv_slice, + comp_slice, + true); + } + + { + const auto expect = static_cast(rank + 1); + Kokkos::parallel_for( + "Check", + CreateRangePolicy({ N_GHOSTS }, { nx2 + N_GHOSTS }), + Lambda(index_t i2) { + for (auto i1 { N_GHOSTS }; i1 < 2 * N_GHOSTS; ++i1) { + if (fld(i1, i2, 0) != expect + 4.0) { + raise::KernelError(HERE, "Left boundary not correct for #0"); + } + if (fld(i1, i2, 1) != expect + 12.0) { + raise::KernelError(HERE, "Left boundary not correct for #1"); + } + if (fld(i1, i2, 2) != expect + 20.0) { + raise::KernelError(HERE, "Left boundary not correct for #2"); + } + } + for (auto i1 { nx1 }; i1 < nx1 + N_GHOSTS; ++i1) { + if (fld(i1, i2, 0) != expect + 4.0) { + raise::KernelError(HERE, "Right boundary not correct for #0"); + } + if (fld(i1, i2, 1) != expect + 12.0) { + raise::KernelError(HERE, "Right boundary not correct for #1"); + } + if (fld(i1, i2, 2) != expect + 20.0) { + raise::KernelError(HERE, "Right boundary not correct for #2"); + } + } + }); + } } catch (std::exception& e) { std::cerr << "Exception: " << e.what() << std::endl; + MPI_Finalize(); Kokkos::finalize(); return 1; } + MPI_Finalize(); Kokkos::finalize(); return 0; } diff --git a/src/framework/tests/comm_nompi.cpp b/src/framework/tests/comm_nompi.cpp index 05d54d589..f9581c1e1 100644 --- a/src/framework/tests/comm_nompi.cpp +++ b/src/framework/tests/comm_nompi.cpp @@ -7,6 +7,8 @@ #include "arch/kokkos_aliases.h" #include "utils/numeric.h" +#include "framework/domain/comm_mpi.hpp" + #include #include @@ -45,12 +47,12 @@ auto main(int argc, char* argv[]) -> int { Kokkos::deep_copy(buff, ZERO); const auto send_slice = std::vector { - {nx1 + N_GHOSTS, nx1 + 2 * N_GHOSTS}, - {nx2 + N_GHOSTS, nx2 + 2 * N_GHOSTS} + { nx1 + N_GHOSTS, nx1 + 2 * N_GHOSTS }, + { nx2 + N_GHOSTS, nx2 + 2 * N_GHOSTS } }; const auto recv_slice = std::vector { - {N_GHOSTS, 2 * N_GHOSTS}, - {N_GHOSTS, 2 * N_GHOSTS} + { N_GHOSTS, 2 * N_GHOSTS }, + { N_GHOSTS, 2 * N_GHOSTS } }; const auto comp_slice = range_tuple_t(cur::jx1, cur::jx3 + 1); From 9b262e3c4b0c5a4a41a2fe58cebf193bfdebb1d4 Mon Sep 17 00:00:00 2001 From: hayk Date: Wed, 23 Apr 2025 09:47:28 -0400 Subject: [PATCH 20/25] deptest fixed --- setups/tests/deposit/pgen.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/setups/tests/deposit/pgen.hpp b/setups/tests/deposit/pgen.hpp index 1be4deb69..dc1bd6380 100644 --- a/setups/tests/deposit/pgen.hpp +++ b/setups/tests/deposit/pgen.hpp @@ -93,8 +93,8 @@ namespace user { data_2["phi"] = phi2s; } - arch::InjectGlobally(global_domain, local_domain, (arch::spidx_t)1, data_1); - arch::InjectGlobally(global_domain, local_domain, (arch::spidx_t)2, data_2); + arch::InjectGlobally(global_domain, local_domain, (spidx_t)1, data_1); + arch::InjectGlobally(global_domain, local_domain, (spidx_t)2, data_2); } }; From b1ebdf17180d1aefcf92b257645d4e1b1ebf35b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludwig=20B=C3=B6ss?= Date: Wed, 23 Apr 2025 16:15:22 -0500 Subject: [PATCH 21/25] higher order Faraday kernel implemented by @vanthieg --- src/engines/srpic.hpp | 36 +++++++++-- src/framework/parameters.cpp | 63 +++++++++++++++++-- src/global/defaults.h | 14 +++++ src/kernels/faraday_mink.hpp | 119 ++++++++++++++++++++++++++++------- 4 files changed, 200 insertions(+), 32 deletions(-) diff --git a/src/engines/srpic.hpp b/src/engines/srpic.hpp index a10070c34..8f8215f1e 100644 --- a/src/engines/srpic.hpp +++ b/src/engines/srpic.hpp @@ -185,6 +185,24 @@ namespace ntt { if constexpr (M::CoordType == Coord::Cart) { // minkowski case const auto dx = math::sqrt(domain.mesh.metric.template h_<1, 1>({})); + const auto deltax = m_params.template get( + "algorithms.fieldsolver.deltax"); + const auto deltay = m_params.template get( + "algorithms.fieldsolver.deltay"); + const auto betaxy = m_params.template get( + "algorithms.fieldsolver.betaxy"); + const auto betayx = m_params.template get( + "algorithms.fieldsolver.betayx"); + const auto deltaz = m_params.template get( + "algorithms.fieldsolver.deltaz"); + const auto betaxz = m_params.template get( + "algorithms.fieldsolver.betaxz"); + const auto betazx = m_params.template get( + "algorithms.fieldsolver.betazx"); + const auto betayz = m_params.template get( + "algorithms.fieldsolver.betayz"); + const auto betazy = m_params.template get( + "algorithms.fieldsolver.betazy"); real_t coeff1, coeff2; if constexpr (M::Dim == Dim::_2D) { coeff1 = dT / SQR(dx); @@ -193,10 +211,20 @@ namespace ntt { coeff1 = dT / dx; coeff2 = ZERO; } - Kokkos::parallel_for( - "Faraday", - domain.mesh.rangeActiveCells(), - kernel::mink::Faraday_kernel(domain.fields.em, coeff1, coeff2)); + Kokkos::parallel_for("Faraday", + domain.mesh.rangeActiveCells(), + kernel::mink::Faraday_kernel(domain.fields.em, + coeff1, + coeff2, + deltax, + deltay, + betaxy, + betayx, + deltaz, + betaxz, + betazx, + betayz, + betazy)); } else { Kokkos::parallel_for("Faraday", domain.mesh.rangeActiveCells(), diff --git a/src/framework/parameters.cpp b/src/framework/parameters.cpp index 8a30f52d9..42655b207 100644 --- a/src/framework/parameters.cpp +++ b/src/framework/parameters.cpp @@ -31,10 +31,10 @@ namespace ntt { template - auto get_dx0_V0( - const std::vector& resolution, - const boundaries_t& extent, - const std::map& params) -> std::pair { + auto get_dx0_V0(const std::vector& resolution, + const boundaries_t& extent, + const std::map& params) + -> std::pair { const auto metric = M(resolution, extent, params); const auto dx0 = metric.dxMin(); coord_t x_corner { ZERO }; @@ -416,6 +416,61 @@ namespace ntt { set("algorithms.toggles.deposit", toml::find_or(toml_data, "algorithms", "toggles", "deposit", true)); + /* [algorithms.fieldsolver] --------------------------------------------- */ + set("algorithms.fieldsolver.deltax", + toml::find_or(toml_data, + "algorithms", + "fieldsolver", + "deltax", + defaults::fieldsolver::deltax)); + set("algorithms.fieldsolver.deltay", + toml::find_or(toml_data, + "algorithms", + "fieldsolver", + "deltay", + defaults::fieldsolver::deltay)); + set("algorithms.fieldsolver.deltaz", + toml::find_or(toml_data, + "algorithms", + "fieldsolver", + "deltaz", + defaults::fieldsolver::deltaz)); + set("algorithms.fieldsolver.betaxy", + toml::find_or(toml_data, + "algorithms", + "fieldsolver", + "betaxy", + defaults::fieldsolver::betaxy)); + set("algorithms.fieldsolver.betayx", + toml::find_or(toml_data, + "algorithms", + "fieldsolver", + "betayx", + defaults::fieldsolver::betayx)); + set("algorithms.fieldsolver.betaxz", + toml::find_or(toml_data, + "algorithms", + "fieldsolver", + "betaxz", + defaults::fieldsolver::betaxz)); + set("algorithms.fieldsolver.betazx", + toml::find_or(toml_data, + "algorithms", + "fieldsolver", + "betazx", + defaults::fieldsolver::betazx)); + set("algorithms.fieldsolver.betayz", + toml::find_or(toml_data, + "algorithms", + "fieldsolver", + "betayz", + defaults::fieldsolver::betayz)); + set("algorithms.fieldsolver.betazy", + toml::find_or(toml_data, + "algorithms", + "fieldsolver", + "betazy", + defaults::fieldsolver::betazy)); /* [algorithms.timestep] ------------------------------------------------ */ set("algorithms.timestep.CFL", toml::find_or(toml_data, "algorithms", "timestep", "CFL", defaults::cfl)); diff --git a/src/global/defaults.h b/src/global/defaults.h index e44103ed0..b5120d6e4 100644 --- a/src/global/defaults.h +++ b/src/global/defaults.h @@ -26,6 +26,20 @@ namespace ntt::defaults { const std::string ph_pusher = "Photon"; const timestep_t clear_interval = 100; + namespace fieldsolver { + const real_t deltax = 0.0; + + const real_t deltay = 0.0; + const real_t betaxy = 0.0; + const real_t betayx = 0.0; + + const real_t deltaz = 0.0; + const real_t betaxz = 0.0; + const real_t betazx = 0.0; + const real_t betayz = 0.0; + const real_t betazy = 0.0; + } // namespace fieldsolver + namespace qsph { const real_t r0 = 0.0; const real_t h = 0.0; diff --git a/src/kernels/faraday_mink.hpp b/src/kernels/faraday_mink.hpp index 6d249b999..1112e56e7 100644 --- a/src/kernels/faraday_mink.hpp +++ b/src/kernels/faraday_mink.hpp @@ -26,6 +26,15 @@ namespace kernel::mink { ndfield_t EB; const real_t coeff1; const real_t coeff2; + const real_t deltax; + const real_t deltay; + const real_t betaxy; + const real_t betayx; + const real_t deltaz; + const real_t betaxz; + const real_t betazx; + const real_t betayz; + const real_t betazy; public: /** @@ -33,15 +42,34 @@ namespace kernel::mink { * ! 2D: coeff1 = dt / dx^2, coeff2 = dt * ! 3D: coeff1 = dt / dx */ - Faraday_kernel(const ndfield_t& EB, real_t coeff1, real_t coeff2) + Faraday_kernel(const ndfield_t& EB, real_t coeff1, real_t coeff2 + , real_t deltax, real_t deltay, real_t betaxy, real_t betayx + , real_t deltaz, real_t betaxz, real_t betazx, real_t betayz + , real_t betazy) : EB { EB } , coeff1 { coeff1 } - , coeff2 { coeff2 } {} + , coeff2 { coeff2 } + , deltax { deltax } + , deltay { deltay } + , betaxy { betaxy } + , betayx { betayx } + , deltaz { deltaz } + , betaxz { betaxz } + , betazx { betazx } + , betayz { betayz } + , betazy { betazy } {} + + Inline void operator()(index_t i1) const { if constexpr (D == Dim::_1D) { - EB(i1, em::bx2) += coeff1 * (EB(i1 + 1, em::ex3) - EB(i1, em::ex3)); - EB(i1, em::bx3) += coeff1 * (EB(i1, em::ex2) - EB(i1 + 1, em::ex2)); + const auto alphax = ONE - THREE * deltax; + EB(i1, em::bx2) += coeff1 * ( + + alphax * (EB(i1 + 1, em::ex3) - EB(i1 , em::ex3)) + + deltax * (EB(i1 + 2, em::ex3) - EB(i1 - 1, em::ex3))); + EB(i1, em::bx3) += coeff1 * ( + - alphax * (EB(i1 + 1, em::ex2) - EB(i1 , em::ex2)) + - deltax * (EB(i1 + 2, em::ex2) - EB(i1 - 1, em::ex2))); } else { raise::KernelError(HERE, "Faraday_kernel: 1D implementation called for D != 1"); } @@ -49,13 +77,27 @@ namespace kernel::mink { Inline void operator()(index_t i1, index_t i2) const { if constexpr (D == Dim::_2D) { - EB(i1, i2, em::bx1) += coeff1 * - (EB(i1, i2, em::ex3) - EB(i1, i2 + 1, em::ex3)); - EB(i1, i2, em::bx2) += coeff1 * - (EB(i1 + 1, i2, em::ex3) - EB(i1, i2, em::ex3)); - EB(i1, i2, em::bx3) += coeff2 * - (EB(i1, i2 + 1, em::ex1) - EB(i1, i2, em::ex1) + - EB(i1, i2, em::ex2) - EB(i1 + 1, i2, em::ex2)); + const auto alphax = ONE - TWO * betaxy - THREE * deltax; + const auto alphay = ONE - TWO * betayx - THREE * deltay; + EB(i1, i2, em::bx1) += coeff1 * ( + - alphay * (EB(i1 , i2 + 1, em::ex3) - EB(i1 , i2 , em::ex3)) + - deltay * (EB(i1 , i2 + 2, em::ex3) - EB(i1 , i2 - 1, em::ex3)) + - betayx * (EB(i1 + 1, i2 + 1, em::ex3) - EB(i1 + 1, i2 , em::ex3)) + - betayx * (EB(i1 - 1, i2 + 1, em::ex3) - EB(i1 - 1, i2 , em::ex3))); + EB(i1, i2, em::bx2) += coeff1 * ( + + alphax * (EB(i1 + 1, i2 , em::ex3) - EB(i1 , i2 , em::ex3)) + + deltax * (EB(i1 + 2, i2 , em::ex3) - EB(i1 - 1, i2 , em::ex3)) + + betaxy * (EB(i1 + 1, i2 + 1, em::ex3) - EB(i1 , i2 + 1, em::ex3)) + + betaxy * (EB(i1 + 1, i2 - 1, em::ex3) - EB(i1 , i2 - 1, em::ex3))); + EB(i1, i2, em::bx3) += coeff2 * ( + + alphay * (EB(i1 , i2 + 1, em::ex1) - EB(i1 , i2 , em::ex1)) + + deltay * (EB(i1 , i2 + 2, em::ex1) - EB(i1 , i2 - 1, em::ex1)) + + betayx * (EB(i1 + 1, i2 + 1, em::ex1) - EB(i1 + 1, i2 , em::ex1)) + + betayx * (EB(i1 - 1, i2 + 1, em::ex1) - EB(i1 - 1, i2 , em::ex1)) + - alphax * (EB(i1 + 1, i2 , em::ex2) - EB(i1 , i2 , em::ex2)) + - deltax * (EB(i1 + 2, i2 , em::ex2) - EB(i1 - 1, i2 , em::ex2)) + - betaxy * (EB(i1 + 1, i2 + 1, em::ex2) - EB(i1 , i2 + 1, em::ex2)) + - betaxy * (EB(i1 + 1, i2 - 1, em::ex2) - EB(i1 , i2 - 1, em::ex2))); } else { raise::KernelError(HERE, "Faraday_kernel: 2D implementation called for D != 2"); @@ -64,19 +106,48 @@ namespace kernel::mink { Inline void operator()(index_t i1, index_t i2, index_t i3) const { if constexpr (D == Dim::_3D) { - EB(i1, i2, i3, em::bx1) += coeff1 * (EB(i1, i2, i3 + 1, em::ex2) - - EB(i1, i2, i3, em::ex2) + - EB(i1, i2, i3, em::ex3) - - EB(i1, i2 + 1, i3, em::ex3)); - EB(i1, i2, i3, em::bx2) += coeff1 * (EB(i1 + 1, i2, i3, em::ex3) - - EB(i1, i2, i3, em::ex3) + - EB(i1, i2, i3, em::ex1) - - EB(i1, i2, i3 + 1, em::ex1)); - EB(i1, i2, i3, em::bx3) += coeff1 * (EB(i1, i2 + 1, i3, em::ex1) - - EB(i1, i2, i3, em::ex1) + - EB(i1, i2, i3, em::ex2) - - EB(i1 + 1, i2, i3, em::ex2)); - + const auto alphax = ONE - TWO * betaxy - TWO * betaxz - THREE * deltax; + const auto alphay = ONE - TWO * betayx - TWO * betayz - THREE * deltay; + const auto alphaz = ONE - TWO * betazx - TWO * betazy - THREE * deltaz; + EB(i1, i2, i3, em::bx1) += coeff1 * ( + + alphaz * (EB(i1 , i2 , i3 + 1, em::ex2) - EB(i1 , i2 , i3 , em::ex2)) + + deltaz * (EB(i1 , i2 , i3 + 2, em::ex2) - EB(i1 , i2 , i3 - 1, em::ex2)) + + betazx * (EB(i1 + 1, i2 , i3 + 1, em::ex2) - EB(i1 + 1, i2 , i3 , em::ex2)) + + betazx * (EB(i1 - 1, i2 , i3 + 1, em::ex2) - EB(i1 - 1, i2 , i3 , em::ex2)) + + betazy * (EB(i1 , i2 + 1, i3 + 1, em::ex2) - EB(i1 , i2 + 1, i3 , em::ex2)) + + betazy * (EB(i1 , i2 - 1, i3 + 1, em::ex2) - EB(i1 , i2 - 1, i3 , em::ex2)) + - alphay * (EB(i1 , i2 + 1, i3 , em::ex3) - EB(i1 , i2 , i3 , em::ex3)) + - deltay * (EB(i1 , i2 + 2, i3 , em::ex3) - EB(i1 , i2 - 1, i3 , em::ex3)) + - betayx * (EB(i1 + 1, i2 + 1, i3 , em::ex3) - EB(i1 + 1, i2 , i3 , em::ex3)) + - betayx * (EB(i1 - 1, i2 + 1, i3 , em::ex3) - EB(i1 - 1, i2 , i3 , em::ex3)) + - betayz * (EB(i1 , i2 + 1, i3 + 1, em::ex3) - EB(i1 , i2 , i3 + 1, em::ex3)) + - betayz * (EB(i1 , i2 + 1, i3 - 1, em::ex3) - EB(i1 , i2 , i3 - 1, em::ex3))); + EB(i1, i2, i3, em::bx2) += coeff1 * ( + + alphax * (EB(i1 + 1, i2 , i3 , em::ex3) - EB(i1 , i2 , i3 , em::ex3)) + + deltax * (EB(i1 + 2, i2 , i3 , em::ex3) - EB(i1 - 1, i2 , i3 , em::ex3)) + + betaxy * (EB(i1 + 1, i2 + 1, i3 , em::ex3) - EB(i1 , i2 + 1, i3 , em::ex3)) + + betaxy * (EB(i1 + 1, i2 - 1, i3 , em::ex3) - EB(i1 , i2 - 1, i3 , em::ex3)) + + betaxz * (EB(i1 + 1, i2 , i3 + 1, em::ex3) - EB(i1 , i2 , i3 + 1, em::ex3)) + + betaxz * (EB(i1 + 1, i2 , i3 - 1, em::ex3) - EB(i1 , i2 , i3 - 1, em::ex3)) + - alphaz * (EB(i1 , i2 , i3 + 1, em::ex1) - EB(i1 , i2 , i3 , em::ex1)) + - deltaz * (EB(i1 , i2 , i3 + 2, em::ex1) - EB(i1 , i2 , i3 - 1, em::ex1)) + - betazx * (EB(i1 + 1, i2 , i3 + 1, em::ex1) - EB(i1 + 1, i2 , i3 , em::ex1)) + - betazx * (EB(i1 - 1, i2 , i3 + 1, em::ex1) - EB(i1 - 1, i2 , i3 , em::ex1)) + - betazy * (EB(i1 , i2 + 1, i3 + 1, em::ex1) - EB(i1 , i2 + 1, i3 , em::ex1)) + - betazy * (EB(i1 , i2 - 1, i3 + 1, em::ex1) - EB(i1 , i2 - 1, i3 , em::ex1))); + EB(i1, i2, i3, em::bx3) += coeff1 * ( + + alphay * (EB(i1 , i2 + 1, i3 , em::ex1) - EB(i1 , i2 , i3 , em::ex1)) + + deltay * (EB(i1 , i2 + 2, i3 , em::ex1) - EB(i1 , i2 - 1, i3 , em::ex1)) + + betayx * (EB(i1 + 1, i2 + 1, i3 , em::ex1) - EB(i1 + 1, i2 , i3 , em::ex1)) + + betayx * (EB(i1 - 1, i2 + 1, i3 , em::ex1) - EB(i1 - 1, i2 , i3 , em::ex1)) + + betayz * (EB(i1 , i2 + 1, i3 + 1, em::ex1) - EB(i1 , i2 , i3 + 1, em::ex1)) + + betayz * (EB(i1 , i2 + 1, i3 - 1, em::ex1) - EB(i1 , i2 , i3 - 1, em::ex1)) + - alphax * (EB(i1 + 1, i2 , i3 , em::ex2) - EB(i1 , i2 , i3 , em::ex2)) + - deltax * (EB(i1 + 2, i2 , i3 , em::ex2) - EB(i1 - 1, i2 , i3 , em::ex2)) + - betaxy * (EB(i1 + 1, i2 + 1, i3 , em::ex2) - EB(i1 , i2 + 1, i3 , em::ex2)) + - betaxy * (EB(i1 + 1, i2 - 1, i3 , em::ex2) - EB(i1 , i2 - 1, i3 , em::ex2)) + - betaxz * (EB(i1 + 1, i2 , i3 + 1, em::ex2) - EB(i1 , i2 , i3 + 1, em::ex2)) + - betaxz * (EB(i1 + 1, i2 , i3 - 1, em::ex2) - EB(i1 , i2 , i3 - 1, em::ex2))); } else { raise::KernelError(HERE, "Faraday_kernel: 3D implementation called for D != 3"); } From 0515a9d9d7abdecc8ab0ccf8a8b0ba46f0d5bec8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludwig=20B=C3=B6ss?= Date: Wed, 23 Apr 2025 16:17:48 -0500 Subject: [PATCH 22/25] example settings for higher order field solver --- input.example.toml | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/input.example.toml b/input.example.toml index bd5d52b06..960d4b7b2 100644 --- a/input.example.toml +++ b/input.example.toml @@ -253,6 +253,21 @@ # @type: float: > 0 gamma_rad = "" + [algorithms.fieldsolver] + # Yee - all 0.0 - default + # 1D + deltax = -0.065 + # 2D + deltay = -0.065 + betaxy = -0.065 + betayx = -0.065 + # 3D - not yet tested + deltaz = 0.0 + betaxz = 0.0 + betazx = 0.0 + betayz = 0.0 + betazy = 0.0 + [particles] # Fiducial number of particles per cell: # @required From 88050f2dd2ffcf6674d08f2755df19295332899b Mon Sep 17 00:00:00 2001 From: hayk Date: Thu, 24 Apr 2025 01:08:43 -0400 Subject: [PATCH 23/25] deposit test pgens --- setups/tests/deposit/deposit-gr.toml | 53 ------------------ setups/tests/deposit/deposit-mink.toml | 44 ++++++++------- setups/tests/deposit/deposit-sr.toml | 42 ++++++++++---- setups/tests/deposit/deposit.toml | 53 ------------------ setups/tests/deposit/pgen.hpp | 76 +++++++++++++++----------- setups/tests/deposit/plot-mink.py | 62 +++++++++++++++++++++ setups/tests/deposit/plot-sr.py | 29 ++++++++++ setups/tests/deposit/run-mink.sh | 14 +++++ 8 files changed, 201 insertions(+), 172 deletions(-) delete mode 100644 setups/tests/deposit/deposit-gr.toml delete mode 100644 setups/tests/deposit/deposit.toml create mode 100644 setups/tests/deposit/plot-mink.py create mode 100644 setups/tests/deposit/plot-sr.py create mode 100755 setups/tests/deposit/run-mink.sh diff --git a/setups/tests/deposit/deposit-gr.toml b/setups/tests/deposit/deposit-gr.toml deleted file mode 100644 index 077e3e59a..000000000 --- a/setups/tests/deposit/deposit-gr.toml +++ /dev/null @@ -1,53 +0,0 @@ -[simulation] - name = "deposit-test" - engine = "srpic" - runtime = 10.0 - -[grid] - resolution = [256, 256] - extent = [[0.0, 1.0], [0.0, 1.0]] - - [grid.metric] - metric = "minkowski" - - [grid.boundaries] - fields = [["PERIODIC"], ["PERIODIC"]] - particles = [["PERIODIC"], ["PERIODIC"]] - -[scales] - larmor0 = 0.1 - skindepth0 = 0.1 - -[algorithms] - current_filters = 4 - - [algorithms.timestep] - CFL = 0.5 - -[particles] - ppc0 = 10.0 - - [[particles.species]] - label = "e-" - mass = 1.0 - charge = -1.0 - maxnpart = 1e2 - - [[particles.species]] - label = "e+" - mass = 1.0 - charge = 1.0 - maxnpart = 1e2 - -[setup] - -[output] - format = "hdf5" - interval_time = 0.01 - - [output.quantities] - quantities = ["N_1", "N_2", "E", "B", "J"] - -[diagnostics] - colored_stdout = true - blocking_timers = true diff --git a/setups/tests/deposit/deposit-mink.toml b/setups/tests/deposit/deposit-mink.toml index df70f7552..2dc953896 100644 --- a/setups/tests/deposit/deposit-mink.toml +++ b/setups/tests/deposit/deposit-mink.toml @@ -1,10 +1,10 @@ [simulation] name = "deposit-test-mink" engine = "srpic" - runtime = 10.0 + runtime = 5.0 [grid] - resolution = [512, 512] + resolution = [32, 32] extent = [[0.0, 1.0], [0.0, 1.0]] [grid.metric] @@ -40,30 +40,32 @@ maxnpart = 1e2 [setup] - x1s = [0.25] - y1s = [0.85] - z1s = [0.33] - ux1s = [0.6] - uy1s = [-0.3] - uz1s = [-0.2] - - x2s = [0.25] - y2s = [0.85] - z2s = [0.33] - ux2s = [-0.2] - uy2s = [-0.2] - uz2s = [0.1] + x1_e = [0.25] + x2_e = [0.85] + x3_e = [0.33] + ux1_e = [0.6] + ux2_e = [-0.3] + ux3_e = [-0.2] + + x1_i = [0.25] + x2_i = [0.85] + x3_i = [0.33] + ux1_i = [-0.6] + ux2_i = [0.3] + ux3_i = [0.2] [output] - format = "hdf5" - interval_time = 0.01 + format = "hdf5" + interval = 5 [output.fields] quantities = ["N_1", "N_2", "E", "B", "J"] + [output.particles] + enable = false + + [output.spectra] + enable = false + [checkpoint] keep = 0 - -[diagnostics] - colored_stdout = true - blocking_timers = true diff --git a/setups/tests/deposit/deposit-sr.toml b/setups/tests/deposit/deposit-sr.toml index 077e3e59a..0e1648d12 100644 --- a/setups/tests/deposit/deposit-sr.toml +++ b/setups/tests/deposit/deposit-sr.toml @@ -1,18 +1,18 @@ [simulation] - name = "deposit-test" + name = "deposit-sr" engine = "srpic" runtime = 10.0 [grid] - resolution = [256, 256] - extent = [[0.0, 1.0], [0.0, 1.0]] + resolution = [64, 64] + extent = [[1.0, 5.0]] [grid.metric] - metric = "minkowski" + metric = "qspherical" [grid.boundaries] - fields = [["PERIODIC"], ["PERIODIC"]] - particles = [["PERIODIC"], ["PERIODIC"]] + fields = [["FIXED", "FIXED"]] + particles = [["REFLECT", "REFLECT"]] [scales] larmor0 = 0.1 @@ -40,14 +40,32 @@ maxnpart = 1e2 [setup] + x1_e = [2.25] + x2_e = [1.25] + phi_e = [0.0] + ux1_e = [0.6] + ux2_e = [-0.3] + ux3_e = [-0.2] + + x1_i = [2.25] + x2_i = [1.25] + phi_i = [0.0] + ux1_i = [-0.6] + ux2_i = [0.3] + ux3_i = [0.2] [output] - format = "hdf5" - interval_time = 0.01 + format = "hdf5" + interval = 5 - [output.quantities] + [output.fields] quantities = ["N_1", "N_2", "E", "B", "J"] -[diagnostics] - colored_stdout = true - blocking_timers = true + [output.particles] + enable = false + + [output.spectra] + enable = false + +[checkpoint] + keep = 0 diff --git a/setups/tests/deposit/deposit.toml b/setups/tests/deposit/deposit.toml deleted file mode 100644 index 04c23ce7d..000000000 --- a/setups/tests/deposit/deposit.toml +++ /dev/null @@ -1,53 +0,0 @@ -[simulation] - name = "deposit-test" - engine = "srpic" - runtime = 1.0 - -[grid] - resolution = [256, 256] - extent = [[0.0, 1.0], [0.0, 1.0]] - - [grid.metric] - metric = "minkowski" - - [grid.boundaries] - fields = [["PERIODIC"], ["PERIODIC"]] - particles = [["PERIODIC"], ["PERIODIC"]] - -[scales] - larmor0 = 0.1 - skindepth0 = 0.1 - -[algorithms] - current_filters = 4 - - [algorithms.timestep] - CFL = 0.5 - -[particles] - ppc0 = 10.0 - - [[particles.species]] - label = "e-" - mass = 1.0 - charge = -1.0 - maxnpart = 1e2 - - [[particles.species]] - label = "e+" - mass = 1.0 - charge = 1.0 - maxnpart = 1e2 - -[setup] - -[output] - format = "hdf5" - interval_time = 0.01 - - [output.quantities] - quantities = ["N_1", "N_2", "E", "B", "J"] - -[diagnostics] - colored_stdout = true - blocking_timers = true diff --git a/setups/tests/deposit/pgen.hpp b/setups/tests/deposit/pgen.hpp index dc1bd6380..0080af8fd 100644 --- a/setups/tests/deposit/pgen.hpp +++ b/setups/tests/deposit/pgen.hpp @@ -48,53 +48,63 @@ namespace user { inline void InitPrtls(Domain& local_domain) { const auto empty = std::vector {}; - const auto x1s = params.template get>("setup.x1s", empty); - const auto y1s = params.template get>("setup.y1s", empty); - const auto z1s = params.template get>("setup.z1s", empty); - const auto phi1s = params.template get>("setup.phi1s", - empty); - const auto ux1s = params.template get>("setup.ux1s", + const auto x1_e = params.template get>("setup.x1_e", empty); - const auto uy1s = params.template get>("setup.uy1s", + const auto x2_e = params.template get>("setup.x2_e", empty); - const auto uz1s = params.template get>("setup.uz1s", + const auto x3_e = params.template get>("setup.x3_e", empty); + const auto phi_e = params.template get>("setup.phi_e", + empty); + const auto ux1_e = params.template get>("setup.ux1_e", + empty); + const auto ux2_e = params.template get>("setup.ux2_e", + empty); + const auto ux3_e = params.template get>("setup.ux3_e", + empty); - const auto x2s = params.template get>("setup.x2s", empty); - const auto y2s = params.template get>("setup.y2s", empty); - const auto z2s = params.template get>("setup.z2s", empty); - const auto ux2s = params.template get>("setup.ux2s", + const auto x1_i = params.template get>("setup.x1_i", empty); - const auto uy2s = params.template get>("setup.uy2s", + const auto x2_i = params.template get>("setup.x2_i", empty); - const auto uz2s = params.template get>("setup.uz2s", + const auto x3_i = params.template get>("setup.x3_i", empty); - const auto phi2s = params.template get>("setup.phi2s", + const auto phi_i = params.template get>("setup.phi_i", empty); - std::map> data_1 { - { "x1", x1s }, - { "x2", y1s }, - { "ux1", ux1s }, - { "ux2", uy1s }, - { "ux3", uz1s } + const auto ux1_i = params.template get>("setup.ux1_i", + empty); + const auto ux2_i = params.template get>("setup.ux2_i", + empty); + const auto ux3_i = params.template get>("setup.ux3_i", + empty); + std::map> data_e { + { "x1", x1_e }, + { "x2", x2_e }, + { "ux1", ux1_e }, + { "ux2", ux2_e }, + { "ux3", ux3_e } }; - std::map> data_2 { - { "x1", x2s }, - { "x2", y2s }, - { "ux1", ux2s }, - { "ux2", uy2s }, - { "ux3", uz2s } + std::map> data_i { + { "x1", x1_i }, + { "x2", x2_i }, + { "ux1", ux1_i }, + { "ux2", ux2_i }, + { "ux3", ux3_i } }; if constexpr (M::CoordType == Coord::Cart or D == Dim::_3D) { - data_1["x3"] = z1s; - data_2["x3"] = z2s; + data_e["x3"] = x3_e; + data_i["x3"] = x3_i; } else if constexpr (D == Dim::_2D) { - data_1["phi"] = phi1s; - data_2["phi"] = phi2s; + data_e["phi"] = phi_e; + data_i["phi"] = phi_i; } - arch::InjectGlobally(global_domain, local_domain, (spidx_t)1, data_1); - arch::InjectGlobally(global_domain, local_domain, (spidx_t)2, data_2); + arch::InjectGlobally(global_domain, local_domain, (spidx_t)1, data_e); + arch::InjectGlobally(global_domain, local_domain, (spidx_t)2, data_i); + } + + auto FixFieldsConst(const bc_in&, const em&) const -> std::pair { + return { ZERO, false }; } }; diff --git a/setups/tests/deposit/plot-mink.py b/setups/tests/deposit/plot-mink.py new file mode 100644 index 000000000..9d6760613 --- /dev/null +++ b/setups/tests/deposit/plot-mink.py @@ -0,0 +1,62 @@ +import nt2 +import matplotlib.pyplot as plt +import matplotlib as mpl + +datas = [] +cpus = [1, 2, 3, 4, 5, 6, 8] +for i in cpus: + datas.append(nt2.Data(path=f"mink-np{i}")) + + +def plot(ti): + fig = plt.figure(figsize=(16, 7), dpi=300) + gs = mpl.gridspec.GridSpec(3, 7, figure=fig) + + for p, quant in enumerate(["Jx", "Jy", "Jz"]): + axs = [fig.add_subplot(gs[p, i]) for i in range(7)] + (datas[0].fields[quant]).isel(t=ti).plot( + ax=axs[0], + cmap="seismic", + add_colorbar=False, + norm=mpl.colors.SymLogNorm( + linthresh=1e-5, + linscale=1, + vmin=-1e-2, + vmax=1e-2, + ), + ) + for i, (d, ax) in enumerate(zip(datas[1:], axs[1:])): + (d.fields[quant] - datas[0].fields[quant]).isel(t=ti).plot( + ax=ax, + cmap="seismic", + add_colorbar=False, + norm=mpl.colors.SymLogNorm( + linthresh=1e-10, + linscale=1, + vmin=-1e-7, + vmax=1e-7, + ), + ) + + for i, ax in enumerate(axs): + ax.set_aspect(1) + if i > 0: + if p == 0: + ax.set_title(f"np{cpus[i]} - np1") + else: + ax.set_title(None) + ax.set_yticklabels([]) + ax.set_ylabel(None) + else: + if p == 0: + ax.set_title(f"np1") + else: + ax.set_title(None) + + if p != 2: + ax.set_xticklabels([]) + ax.set_xlabel(None) + + +nt2.export.makeFrames(plot, datas[0].fields.s[::4], "mink-diff", num_cpus=4) +nt2.export.makeMovie(framerate=10, input="mink-diff/", number=5, output="mink-diff.mp4") diff --git a/setups/tests/deposit/plot-sr.py b/setups/tests/deposit/plot-sr.py new file mode 100644 index 000000000..0ba5dff24 --- /dev/null +++ b/setups/tests/deposit/plot-sr.py @@ -0,0 +1,29 @@ +import nt2 +import matplotlib.pyplot as plt +import matplotlib as mpl + +data = nt2.Data(path=f"sr-np8") + + +def plot(ti): + fig = plt.figure(figsize=(9, 6), dpi=300) + gs = mpl.gridspec.GridSpec(1, 3, figure=fig) + axs = [fig.add_subplot(gs[0, i]) for i in range(3)] + for i, (ax, j) in enumerate(zip(axs, ["Jr", "Jth", "Jph"])): + data.fields.isel(t=ti)[j].polar.pcolor( + ax=ax, + cbar_position="top", + cbar_size="2%", + norm=mpl.colors.SymLogNorm(linthresh=1e-8, vmin=-1e-4, vmax=1e-4), + cmap="seismic", + ) + ax.set_title(None) + ax.add_artist(mpl.patches.Circle((0, 0), 1, color="k", alpha=0.2)) + ax.add_artist(mpl.patches.Circle((0, 0), 5, edgecolor="k", facecolor="none")) + if i > 0: + ax.set_yticklabels([]) + ax.set_ylabel(None) + + +nt2.export.makeFrames(plot, data.fields.s, "sr-dep", num_cpus=4) +nt2.export.makeMovie(framerate=10, input="sr-dep/", number=5, output="sr-dep.mp4") diff --git a/setups/tests/deposit/run-mink.sh b/setups/tests/deposit/run-mink.sh new file mode 100755 index 000000000..4b52d6642 --- /dev/null +++ b/setups/tests/deposit/run-mink.sh @@ -0,0 +1,14 @@ +#!/usr/bin/env bash + +for i in {1..8}; do + if [ $i -eq 7 ]; then + continue + fi + run=$(echo "np${i}") + cp deposit-mink.toml deposit-mink-${run}.toml && \ + sed -i 's/name[[:space:]]*=[[:space:]]*".*\?"/name = "mink-'${run}'"/g' deposit-mink-${run}.toml && \ + mpiexec -np ${i} ./entity.xc -input deposit-mink-${run}.toml && \ + rm deposit-mink-${run}.toml +done + +rm *.info *.err *.log *.csv From 1eccfa4a774dc59a7c6bdc2e13f45796a5f38b80 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludwig=20B=C3=B6ss?= Date: Thu, 24 Apr 2025 12:58:05 -0500 Subject: [PATCH 24/25] Revert "example settings for higher order field solver" This reverts commit 0515a9d9d7abdecc8ab0ccf8a8b0ba46f0d5bec8. --- input.example.toml | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/input.example.toml b/input.example.toml index 960d4b7b2..bd5d52b06 100644 --- a/input.example.toml +++ b/input.example.toml @@ -253,21 +253,6 @@ # @type: float: > 0 gamma_rad = "" - [algorithms.fieldsolver] - # Yee - all 0.0 - default - # 1D - deltax = -0.065 - # 2D - deltay = -0.065 - betaxy = -0.065 - betayx = -0.065 - # 3D - not yet tested - deltaz = 0.0 - betaxz = 0.0 - betazx = 0.0 - betayz = 0.0 - betazy = 0.0 - [particles] # Fiducial number of particles per cell: # @required From 82cf56b3fffe699836cb75b16de7eec0af2cc4ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludwig=20B=C3=B6ss?= Date: Thu, 24 Apr 2025 12:58:30 -0500 Subject: [PATCH 25/25] Revert "higher order Faraday kernel implemented by @vanthieg" This reverts commit b1ebdf17180d1aefcf92b257645d4e1b1ebf35b3. --- src/engines/srpic.hpp | 36 ++--------- src/framework/parameters.cpp | 63 ++----------------- src/global/defaults.h | 14 ----- src/kernels/faraday_mink.hpp | 119 +++++++---------------------------- 4 files changed, 32 insertions(+), 200 deletions(-) diff --git a/src/engines/srpic.hpp b/src/engines/srpic.hpp index 8f8215f1e..a10070c34 100644 --- a/src/engines/srpic.hpp +++ b/src/engines/srpic.hpp @@ -185,24 +185,6 @@ namespace ntt { if constexpr (M::CoordType == Coord::Cart) { // minkowski case const auto dx = math::sqrt(domain.mesh.metric.template h_<1, 1>({})); - const auto deltax = m_params.template get( - "algorithms.fieldsolver.deltax"); - const auto deltay = m_params.template get( - "algorithms.fieldsolver.deltay"); - const auto betaxy = m_params.template get( - "algorithms.fieldsolver.betaxy"); - const auto betayx = m_params.template get( - "algorithms.fieldsolver.betayx"); - const auto deltaz = m_params.template get( - "algorithms.fieldsolver.deltaz"); - const auto betaxz = m_params.template get( - "algorithms.fieldsolver.betaxz"); - const auto betazx = m_params.template get( - "algorithms.fieldsolver.betazx"); - const auto betayz = m_params.template get( - "algorithms.fieldsolver.betayz"); - const auto betazy = m_params.template get( - "algorithms.fieldsolver.betazy"); real_t coeff1, coeff2; if constexpr (M::Dim == Dim::_2D) { coeff1 = dT / SQR(dx); @@ -211,20 +193,10 @@ namespace ntt { coeff1 = dT / dx; coeff2 = ZERO; } - Kokkos::parallel_for("Faraday", - domain.mesh.rangeActiveCells(), - kernel::mink::Faraday_kernel(domain.fields.em, - coeff1, - coeff2, - deltax, - deltay, - betaxy, - betayx, - deltaz, - betaxz, - betazx, - betayz, - betazy)); + Kokkos::parallel_for( + "Faraday", + domain.mesh.rangeActiveCells(), + kernel::mink::Faraday_kernel(domain.fields.em, coeff1, coeff2)); } else { Kokkos::parallel_for("Faraday", domain.mesh.rangeActiveCells(), diff --git a/src/framework/parameters.cpp b/src/framework/parameters.cpp index 42655b207..8a30f52d9 100644 --- a/src/framework/parameters.cpp +++ b/src/framework/parameters.cpp @@ -31,10 +31,10 @@ namespace ntt { template - auto get_dx0_V0(const std::vector& resolution, - const boundaries_t& extent, - const std::map& params) - -> std::pair { + auto get_dx0_V0( + const std::vector& resolution, + const boundaries_t& extent, + const std::map& params) -> std::pair { const auto metric = M(resolution, extent, params); const auto dx0 = metric.dxMin(); coord_t x_corner { ZERO }; @@ -416,61 +416,6 @@ namespace ntt { set("algorithms.toggles.deposit", toml::find_or(toml_data, "algorithms", "toggles", "deposit", true)); - /* [algorithms.fieldsolver] --------------------------------------------- */ - set("algorithms.fieldsolver.deltax", - toml::find_or(toml_data, - "algorithms", - "fieldsolver", - "deltax", - defaults::fieldsolver::deltax)); - set("algorithms.fieldsolver.deltay", - toml::find_or(toml_data, - "algorithms", - "fieldsolver", - "deltay", - defaults::fieldsolver::deltay)); - set("algorithms.fieldsolver.deltaz", - toml::find_or(toml_data, - "algorithms", - "fieldsolver", - "deltaz", - defaults::fieldsolver::deltaz)); - set("algorithms.fieldsolver.betaxy", - toml::find_or(toml_data, - "algorithms", - "fieldsolver", - "betaxy", - defaults::fieldsolver::betaxy)); - set("algorithms.fieldsolver.betayx", - toml::find_or(toml_data, - "algorithms", - "fieldsolver", - "betayx", - defaults::fieldsolver::betayx)); - set("algorithms.fieldsolver.betaxz", - toml::find_or(toml_data, - "algorithms", - "fieldsolver", - "betaxz", - defaults::fieldsolver::betaxz)); - set("algorithms.fieldsolver.betazx", - toml::find_or(toml_data, - "algorithms", - "fieldsolver", - "betazx", - defaults::fieldsolver::betazx)); - set("algorithms.fieldsolver.betayz", - toml::find_or(toml_data, - "algorithms", - "fieldsolver", - "betayz", - defaults::fieldsolver::betayz)); - set("algorithms.fieldsolver.betazy", - toml::find_or(toml_data, - "algorithms", - "fieldsolver", - "betazy", - defaults::fieldsolver::betazy)); /* [algorithms.timestep] ------------------------------------------------ */ set("algorithms.timestep.CFL", toml::find_or(toml_data, "algorithms", "timestep", "CFL", defaults::cfl)); diff --git a/src/global/defaults.h b/src/global/defaults.h index b5120d6e4..e44103ed0 100644 --- a/src/global/defaults.h +++ b/src/global/defaults.h @@ -26,20 +26,6 @@ namespace ntt::defaults { const std::string ph_pusher = "Photon"; const timestep_t clear_interval = 100; - namespace fieldsolver { - const real_t deltax = 0.0; - - const real_t deltay = 0.0; - const real_t betaxy = 0.0; - const real_t betayx = 0.0; - - const real_t deltaz = 0.0; - const real_t betaxz = 0.0; - const real_t betazx = 0.0; - const real_t betayz = 0.0; - const real_t betazy = 0.0; - } // namespace fieldsolver - namespace qsph { const real_t r0 = 0.0; const real_t h = 0.0; diff --git a/src/kernels/faraday_mink.hpp b/src/kernels/faraday_mink.hpp index 1112e56e7..6d249b999 100644 --- a/src/kernels/faraday_mink.hpp +++ b/src/kernels/faraday_mink.hpp @@ -26,15 +26,6 @@ namespace kernel::mink { ndfield_t EB; const real_t coeff1; const real_t coeff2; - const real_t deltax; - const real_t deltay; - const real_t betaxy; - const real_t betayx; - const real_t deltaz; - const real_t betaxz; - const real_t betazx; - const real_t betayz; - const real_t betazy; public: /** @@ -42,34 +33,15 @@ namespace kernel::mink { * ! 2D: coeff1 = dt / dx^2, coeff2 = dt * ! 3D: coeff1 = dt / dx */ - Faraday_kernel(const ndfield_t& EB, real_t coeff1, real_t coeff2 - , real_t deltax, real_t deltay, real_t betaxy, real_t betayx - , real_t deltaz, real_t betaxz, real_t betazx, real_t betayz - , real_t betazy) + Faraday_kernel(const ndfield_t& EB, real_t coeff1, real_t coeff2) : EB { EB } , coeff1 { coeff1 } - , coeff2 { coeff2 } - , deltax { deltax } - , deltay { deltay } - , betaxy { betaxy } - , betayx { betayx } - , deltaz { deltaz } - , betaxz { betaxz } - , betazx { betazx } - , betayz { betayz } - , betazy { betazy } {} - - + , coeff2 { coeff2 } {} Inline void operator()(index_t i1) const { if constexpr (D == Dim::_1D) { - const auto alphax = ONE - THREE * deltax; - EB(i1, em::bx2) += coeff1 * ( - + alphax * (EB(i1 + 1, em::ex3) - EB(i1 , em::ex3)) - + deltax * (EB(i1 + 2, em::ex3) - EB(i1 - 1, em::ex3))); - EB(i1, em::bx3) += coeff1 * ( - - alphax * (EB(i1 + 1, em::ex2) - EB(i1 , em::ex2)) - - deltax * (EB(i1 + 2, em::ex2) - EB(i1 - 1, em::ex2))); + EB(i1, em::bx2) += coeff1 * (EB(i1 + 1, em::ex3) - EB(i1, em::ex3)); + EB(i1, em::bx3) += coeff1 * (EB(i1, em::ex2) - EB(i1 + 1, em::ex2)); } else { raise::KernelError(HERE, "Faraday_kernel: 1D implementation called for D != 1"); } @@ -77,27 +49,13 @@ namespace kernel::mink { Inline void operator()(index_t i1, index_t i2) const { if constexpr (D == Dim::_2D) { - const auto alphax = ONE - TWO * betaxy - THREE * deltax; - const auto alphay = ONE - TWO * betayx - THREE * deltay; - EB(i1, i2, em::bx1) += coeff1 * ( - - alphay * (EB(i1 , i2 + 1, em::ex3) - EB(i1 , i2 , em::ex3)) - - deltay * (EB(i1 , i2 + 2, em::ex3) - EB(i1 , i2 - 1, em::ex3)) - - betayx * (EB(i1 + 1, i2 + 1, em::ex3) - EB(i1 + 1, i2 , em::ex3)) - - betayx * (EB(i1 - 1, i2 + 1, em::ex3) - EB(i1 - 1, i2 , em::ex3))); - EB(i1, i2, em::bx2) += coeff1 * ( - + alphax * (EB(i1 + 1, i2 , em::ex3) - EB(i1 , i2 , em::ex3)) - + deltax * (EB(i1 + 2, i2 , em::ex3) - EB(i1 - 1, i2 , em::ex3)) - + betaxy * (EB(i1 + 1, i2 + 1, em::ex3) - EB(i1 , i2 + 1, em::ex3)) - + betaxy * (EB(i1 + 1, i2 - 1, em::ex3) - EB(i1 , i2 - 1, em::ex3))); - EB(i1, i2, em::bx3) += coeff2 * ( - + alphay * (EB(i1 , i2 + 1, em::ex1) - EB(i1 , i2 , em::ex1)) - + deltay * (EB(i1 , i2 + 2, em::ex1) - EB(i1 , i2 - 1, em::ex1)) - + betayx * (EB(i1 + 1, i2 + 1, em::ex1) - EB(i1 + 1, i2 , em::ex1)) - + betayx * (EB(i1 - 1, i2 + 1, em::ex1) - EB(i1 - 1, i2 , em::ex1)) - - alphax * (EB(i1 + 1, i2 , em::ex2) - EB(i1 , i2 , em::ex2)) - - deltax * (EB(i1 + 2, i2 , em::ex2) - EB(i1 - 1, i2 , em::ex2)) - - betaxy * (EB(i1 + 1, i2 + 1, em::ex2) - EB(i1 , i2 + 1, em::ex2)) - - betaxy * (EB(i1 + 1, i2 - 1, em::ex2) - EB(i1 , i2 - 1, em::ex2))); + EB(i1, i2, em::bx1) += coeff1 * + (EB(i1, i2, em::ex3) - EB(i1, i2 + 1, em::ex3)); + EB(i1, i2, em::bx2) += coeff1 * + (EB(i1 + 1, i2, em::ex3) - EB(i1, i2, em::ex3)); + EB(i1, i2, em::bx3) += coeff2 * + (EB(i1, i2 + 1, em::ex1) - EB(i1, i2, em::ex1) + + EB(i1, i2, em::ex2) - EB(i1 + 1, i2, em::ex2)); } else { raise::KernelError(HERE, "Faraday_kernel: 2D implementation called for D != 2"); @@ -106,48 +64,19 @@ namespace kernel::mink { Inline void operator()(index_t i1, index_t i2, index_t i3) const { if constexpr (D == Dim::_3D) { - const auto alphax = ONE - TWO * betaxy - TWO * betaxz - THREE * deltax; - const auto alphay = ONE - TWO * betayx - TWO * betayz - THREE * deltay; - const auto alphaz = ONE - TWO * betazx - TWO * betazy - THREE * deltaz; - EB(i1, i2, i3, em::bx1) += coeff1 * ( - + alphaz * (EB(i1 , i2 , i3 + 1, em::ex2) - EB(i1 , i2 , i3 , em::ex2)) - + deltaz * (EB(i1 , i2 , i3 + 2, em::ex2) - EB(i1 , i2 , i3 - 1, em::ex2)) - + betazx * (EB(i1 + 1, i2 , i3 + 1, em::ex2) - EB(i1 + 1, i2 , i3 , em::ex2)) - + betazx * (EB(i1 - 1, i2 , i3 + 1, em::ex2) - EB(i1 - 1, i2 , i3 , em::ex2)) - + betazy * (EB(i1 , i2 + 1, i3 + 1, em::ex2) - EB(i1 , i2 + 1, i3 , em::ex2)) - + betazy * (EB(i1 , i2 - 1, i3 + 1, em::ex2) - EB(i1 , i2 - 1, i3 , em::ex2)) - - alphay * (EB(i1 , i2 + 1, i3 , em::ex3) - EB(i1 , i2 , i3 , em::ex3)) - - deltay * (EB(i1 , i2 + 2, i3 , em::ex3) - EB(i1 , i2 - 1, i3 , em::ex3)) - - betayx * (EB(i1 + 1, i2 + 1, i3 , em::ex3) - EB(i1 + 1, i2 , i3 , em::ex3)) - - betayx * (EB(i1 - 1, i2 + 1, i3 , em::ex3) - EB(i1 - 1, i2 , i3 , em::ex3)) - - betayz * (EB(i1 , i2 + 1, i3 + 1, em::ex3) - EB(i1 , i2 , i3 + 1, em::ex3)) - - betayz * (EB(i1 , i2 + 1, i3 - 1, em::ex3) - EB(i1 , i2 , i3 - 1, em::ex3))); - EB(i1, i2, i3, em::bx2) += coeff1 * ( - + alphax * (EB(i1 + 1, i2 , i3 , em::ex3) - EB(i1 , i2 , i3 , em::ex3)) - + deltax * (EB(i1 + 2, i2 , i3 , em::ex3) - EB(i1 - 1, i2 , i3 , em::ex3)) - + betaxy * (EB(i1 + 1, i2 + 1, i3 , em::ex3) - EB(i1 , i2 + 1, i3 , em::ex3)) - + betaxy * (EB(i1 + 1, i2 - 1, i3 , em::ex3) - EB(i1 , i2 - 1, i3 , em::ex3)) - + betaxz * (EB(i1 + 1, i2 , i3 + 1, em::ex3) - EB(i1 , i2 , i3 + 1, em::ex3)) - + betaxz * (EB(i1 + 1, i2 , i3 - 1, em::ex3) - EB(i1 , i2 , i3 - 1, em::ex3)) - - alphaz * (EB(i1 , i2 , i3 + 1, em::ex1) - EB(i1 , i2 , i3 , em::ex1)) - - deltaz * (EB(i1 , i2 , i3 + 2, em::ex1) - EB(i1 , i2 , i3 - 1, em::ex1)) - - betazx * (EB(i1 + 1, i2 , i3 + 1, em::ex1) - EB(i1 + 1, i2 , i3 , em::ex1)) - - betazx * (EB(i1 - 1, i2 , i3 + 1, em::ex1) - EB(i1 - 1, i2 , i3 , em::ex1)) - - betazy * (EB(i1 , i2 + 1, i3 + 1, em::ex1) - EB(i1 , i2 + 1, i3 , em::ex1)) - - betazy * (EB(i1 , i2 - 1, i3 + 1, em::ex1) - EB(i1 , i2 - 1, i3 , em::ex1))); - EB(i1, i2, i3, em::bx3) += coeff1 * ( - + alphay * (EB(i1 , i2 + 1, i3 , em::ex1) - EB(i1 , i2 , i3 , em::ex1)) - + deltay * (EB(i1 , i2 + 2, i3 , em::ex1) - EB(i1 , i2 - 1, i3 , em::ex1)) - + betayx * (EB(i1 + 1, i2 + 1, i3 , em::ex1) - EB(i1 + 1, i2 , i3 , em::ex1)) - + betayx * (EB(i1 - 1, i2 + 1, i3 , em::ex1) - EB(i1 - 1, i2 , i3 , em::ex1)) - + betayz * (EB(i1 , i2 + 1, i3 + 1, em::ex1) - EB(i1 , i2 , i3 + 1, em::ex1)) - + betayz * (EB(i1 , i2 + 1, i3 - 1, em::ex1) - EB(i1 , i2 , i3 - 1, em::ex1)) - - alphax * (EB(i1 + 1, i2 , i3 , em::ex2) - EB(i1 , i2 , i3 , em::ex2)) - - deltax * (EB(i1 + 2, i2 , i3 , em::ex2) - EB(i1 - 1, i2 , i3 , em::ex2)) - - betaxy * (EB(i1 + 1, i2 + 1, i3 , em::ex2) - EB(i1 , i2 + 1, i3 , em::ex2)) - - betaxy * (EB(i1 + 1, i2 - 1, i3 , em::ex2) - EB(i1 , i2 - 1, i3 , em::ex2)) - - betaxz * (EB(i1 + 1, i2 , i3 + 1, em::ex2) - EB(i1 , i2 , i3 + 1, em::ex2)) - - betaxz * (EB(i1 + 1, i2 , i3 - 1, em::ex2) - EB(i1 , i2 , i3 - 1, em::ex2))); + EB(i1, i2, i3, em::bx1) += coeff1 * (EB(i1, i2, i3 + 1, em::ex2) - + EB(i1, i2, i3, em::ex2) + + EB(i1, i2, i3, em::ex3) - + EB(i1, i2 + 1, i3, em::ex3)); + EB(i1, i2, i3, em::bx2) += coeff1 * (EB(i1 + 1, i2, i3, em::ex3) - + EB(i1, i2, i3, em::ex3) + + EB(i1, i2, i3, em::ex1) - + EB(i1, i2, i3 + 1, em::ex1)); + EB(i1, i2, i3, em::bx3) += coeff1 * (EB(i1, i2 + 1, i3, em::ex1) - + EB(i1, i2, i3, em::ex1) + + EB(i1, i2, i3, em::ex2) - + EB(i1 + 1, i2, i3, em::ex2)); + } else { raise::KernelError(HERE, "Faraday_kernel: 3D implementation called for D != 3"); }