From c212e3a4b2f5ced7f99766007ff7f3d2d3047b7d Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 3 Feb 2026 14:48:16 -0500 Subject: [PATCH 01/36] psc_bnd_fields_impl: rm old fixme --- src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx | 1 - 1 file changed, 1 deletion(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 2e8c18516..3a3d7b67a 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -7,7 +7,6 @@ #include -// FIXME, needs public access to Fields::ib, im // #define DEBUG template From 86bb9ffb9dad71ef49947ece2898d29eac8a8f86 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 3 Feb 2026 14:48:49 -0500 Subject: [PATCH 02/36] psc_bnd_fields_impl: rename to MfieldsState --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 42 +++++++++---------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 3a3d7b67a..0889183f4 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -9,19 +9,19 @@ // #define DEBUG -template +template struct BndFields_ : BndFieldsBase { - using Self = BndFields_; - using Mfields = MF; - using real_t = typename Mfields::real_t; - using fields_view_t = typename Mfields::fields_view_t; + using Self = BndFields_; + using MfieldsState = MFIELDS_STATE; + using real_t = typename MfieldsState::real_t; + using fields_view_t = typename MfieldsState::fields_view_t; using dim_t = Dim; // ---------------------------------------------------------------------- // fill_ghosts_E - void fill_ghosts_E(Mfields& mflds) + void fill_ghosts_E(MfieldsState& mflds) { const auto& grid = mflds.grid(); @@ -59,7 +59,7 @@ struct BndFields_ : BndFieldsBase // ---------------------------------------------------------------------- // fill_ghosts_H - void fill_ghosts_H(Mfields& mflds) + void fill_ghosts_H(MfieldsState& mflds) { const auto& grid = mflds.grid(); @@ -96,7 +96,7 @@ struct BndFields_ : BndFieldsBase // ---------------------------------------------------------------------- // add_ghosts_J - void add_ghosts_J(Mfields& mflds) + void add_ghosts_J(MfieldsState& mflds) { const auto& grid = mflds.grid(); @@ -135,7 +135,7 @@ struct BndFields_ : BndFieldsBase *f = std::numeric_limits::quiet_NaN(); } - void conducting_wall_E_lo(Mfields& mflds, int p, int d) + void conducting_wall_E_lo(MfieldsState& mflds, int p, int d) { auto F = make_Fields3d(mflds[p]); const int* ldims = mflds.grid().ldims; @@ -198,7 +198,7 @@ struct BndFields_ : BndFieldsBase } } - void conducting_wall_E_hi(Mfields& mflds, int p, int d) + void conducting_wall_E_hi(MfieldsState& mflds, int p, int d) { auto F = make_Fields3d(mflds[p]); const int* ldims = mflds.grid().ldims; @@ -262,7 +262,7 @@ struct BndFields_ : BndFieldsBase } } - void conducting_wall_H_lo(Mfields& mflds, int p, int d) + void conducting_wall_H_lo(MfieldsState& mflds, int p, int d) { auto F = make_Fields3d(mflds[p]); const int* ldims = mflds.grid().ldims; @@ -320,7 +320,7 @@ struct BndFields_ : BndFieldsBase } } - void conducting_wall_H_hi(Mfields& mflds, int p, int d) + void conducting_wall_H_hi(MfieldsState& mflds, int p, int d) { auto F = make_Fields3d(mflds[p]); @@ -379,7 +379,7 @@ struct BndFields_ : BndFieldsBase } } - void conducting_wall_J_lo(Mfields& mflds, int p, int d) + void conducting_wall_J_lo(MfieldsState& mflds, int p, int d) { auto F = make_Fields3d(mflds[p]); const int* ldims = mflds.grid().ldims; @@ -418,7 +418,7 @@ struct BndFields_ : BndFieldsBase } } - void conducting_wall_J_hi(Mfields& mflds, int p, int d) + void conducting_wall_J_hi(MfieldsState& mflds, int p, int d) { auto F = make_Fields3d(mflds[p]); const int* ldims = mflds.grid().ldims; @@ -465,7 +465,7 @@ struct BndFields_ : BndFieldsBase // ---------------------------------------------------------------------- // open_H_lo - void open_H_lo(Mfields& mflds, int p, int d) + void open_H_lo(MfieldsState& mflds, int p, int d) { auto F = make_Fields3d(mflds[p]); const Grid_t& grid = mflds.grid(); @@ -549,7 +549,7 @@ struct BndFields_ : BndFieldsBase #endif } - void open_H_hi(Mfields& mflds, int p, int d) + void open_H_hi(MfieldsState& mflds, int p, int d) { auto F = make_Fields3d(mflds[p]); const Grid_t& grid = mflds.grid(); @@ -641,14 +641,14 @@ struct BndFields_ : BndFieldsBase // used by CUDA -template +template struct BndFieldsNone : BndFieldsBase { - using Mfields = MF; + using MfieldsState = MFIELDS_STATE; // clang-format off - void fill_ghosts_E(Mfields& mflds) {}; - void fill_ghosts_H(Mfields& mflds) {}; - void add_ghosts_J(Mfields& mflds) {}; + void fill_ghosts_E(MfieldsState& mflds) {}; + void fill_ghosts_H(MfieldsState& mflds) {}; + void add_ghosts_J(MfieldsState& mflds) {}; // clang-format on }; From 0f15520b67b1cf7378487ace3c62213a1f90ca57 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 3 Feb 2026 14:51:56 -0500 Subject: [PATCH 03/36] psc_bnd_fields_impl: rename open -> radiative --- src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 0889183f4..75159ce7a 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -72,7 +72,7 @@ struct BndFields_ : BndFieldsBase case BND_FLD_CONDUCTING_WALL: conducting_wall_H_lo(mflds, p, d); break; - case BND_FLD_OPEN: open_H_lo(mflds, p, d); break; + case BND_FLD_OPEN: radiative_H_lo(mflds, p, d); break; default: assert(0); } } @@ -85,7 +85,7 @@ struct BndFields_ : BndFieldsBase case BND_FLD_CONDUCTING_WALL: conducting_wall_H_hi(mflds, p, d); break; - case BND_FLD_OPEN: open_H_hi(mflds, p, d); break; + case BND_FLD_OPEN: radiative_H_hi(mflds, p, d); break; default: assert(0); } } @@ -459,13 +459,7 @@ struct BndFields_ : BndFieldsBase } } - // ====================================================================== - // open - - // ---------------------------------------------------------------------- - // open_H_lo - - void open_H_lo(MfieldsState& mflds, int p, int d) + void radiative_H_lo(MfieldsState& mflds, int p, int d) { auto F = make_Fields3d(mflds[p]); const Grid_t& grid = mflds.grid(); @@ -549,7 +543,7 @@ struct BndFields_ : BndFieldsBase #endif } - void open_H_hi(MfieldsState& mflds, int p, int d) + void radiative_H_hi(MfieldsState& mflds, int p, int d) { auto F = make_Fields3d(mflds[p]); const Grid_t& grid = mflds.grid(); From d82a7cf5f4d5d9ae5a58013ac3d0064da7629c1b Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 3 Feb 2026 14:55:06 -0500 Subject: [PATCH 04/36] psc_bnd_fields_impl: no fallthrough --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 92 ++++++++++++++----- 1 file changed, 68 insertions(+), 24 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 75159ce7a..bca20ce26 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -30,12 +30,19 @@ struct BndFields_ : BndFieldsBase for (int d = 0; d < 3; d++) { if (grid.atBoundaryLo(p, d)) { switch (grid.bc.fld_lo[d]) { - case BND_FLD_PERIODIC: break; - case BND_FLD_CONDUCTING_WALL: + case BND_FLD_PERIODIC: { + break; + } + case BND_FLD_CONDUCTING_WALL: { conducting_wall_E_lo(mflds, p, d); break; - case BND_FLD_OPEN: break; - default: assert(0); + } + case BND_FLD_OPEN: { + break; + } + default: { + assert(0); + } } } } @@ -44,12 +51,19 @@ struct BndFields_ : BndFieldsBase for (int d = 0; d < 3; d++) { if (grid.atBoundaryHi(p, d)) { switch (grid.bc.fld_hi[d]) { - case BND_FLD_PERIODIC: break; - case BND_FLD_CONDUCTING_WALL: + case BND_FLD_PERIODIC: { + break; + } + case BND_FLD_CONDUCTING_WALL: { conducting_wall_E_hi(mflds, p, d); break; - case BND_FLD_OPEN: break; - default: assert(0); + } + case BND_FLD_OPEN: { + break; + } + default: { + assert(0); + } } } } @@ -68,12 +82,20 @@ struct BndFields_ : BndFieldsBase for (int d = 0; d < 3; d++) { if (grid.atBoundaryLo(p, d)) { switch (grid.bc.fld_lo[d]) { - case BND_FLD_PERIODIC: break; - case BND_FLD_CONDUCTING_WALL: + case BND_FLD_PERIODIC: { + break; + } + case BND_FLD_CONDUCTING_WALL: { conducting_wall_H_lo(mflds, p, d); break; - case BND_FLD_OPEN: radiative_H_lo(mflds, p, d); break; - default: assert(0); + } + case BND_FLD_OPEN: { + radiative_H_lo(mflds, p, d); + break; + } + default: { + assert(0); + } } } } @@ -81,12 +103,20 @@ struct BndFields_ : BndFieldsBase for (int d = 0; d < 3; d++) { if (grid.atBoundaryHi(p, d)) { switch (grid.bc.fld_hi[d]) { - case BND_FLD_PERIODIC: break; - case BND_FLD_CONDUCTING_WALL: + case BND_FLD_PERIODIC: { + break; + } + case BND_FLD_CONDUCTING_WALL: { conducting_wall_H_hi(mflds, p, d); break; - case BND_FLD_OPEN: radiative_H_hi(mflds, p, d); break; - default: assert(0); + } + case BND_FLD_OPEN: { + radiative_H_hi(mflds, p, d); + break; + } + default: { + assert(0); + } } } } @@ -105,12 +135,19 @@ struct BndFields_ : BndFieldsBase for (int d = 0; d < 3; d++) { if (grid.atBoundaryLo(p, d)) { switch (grid.bc.fld_lo[d]) { - case BND_FLD_PERIODIC: - case BND_FLD_OPEN: break; - case BND_FLD_CONDUCTING_WALL: + case BND_FLD_PERIODIC: { + break; + } + case BND_FLD_CONDUCTING_WALL: { conducting_wall_J_lo(mflds, p, d); break; - default: assert(0); + } + case BND_FLD_OPEN: { + break; + } + default: { + assert(0); + } } } } @@ -118,12 +155,19 @@ struct BndFields_ : BndFieldsBase for (int d = 0; d < 3; d++) { if (grid.atBoundaryHi(p, d)) { switch (grid.bc.fld_hi[d]) { - case BND_FLD_PERIODIC: - case BND_FLD_OPEN: break; - case BND_FLD_CONDUCTING_WALL: + case BND_FLD_PERIODIC: { + break; + } + case BND_FLD_CONDUCTING_WALL: { conducting_wall_J_hi(mflds, p, d); break; - default: assert(0); + } + case BND_FLD_OPEN: { + break; + } + default: { + assert(0); + } } } } From d6c73d49378b819ea801efc916b83e5dda3be754 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 3 Feb 2026 15:38:53 -0500 Subject: [PATCH 05/36] fields: index by Int3 --- src/include/fields.hxx | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/include/fields.hxx b/src/include/fields.hxx index 7696df04a..1194c7b7d 100644 --- a/src/include/fields.hxx +++ b/src/include/fields.hxx @@ -65,6 +65,16 @@ public: return type_traits::get(e_(i, j, k, m)); } + GT_INLINE const_reference operator()(int m, const Int3& i3) const + { + return (*this)(m, i3[0], i3[1], i3[2]); + } + + GT_INLINE reference operator()(int m, const Int3& i3) + { + return (*this)(m, i3[0], i3[1], i3[2]); + } + private: fields_t e_; Int3 ib_; From dcf3c2f43457a98a7e4349d5b253e1e9e0f9e2dc Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 4 Feb 2026 10:22:52 -0500 Subject: [PATCH 06/36] psc_bnd_fields_impl: +background_e,h and bcs --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 77 +++++++++++++++++++ 1 file changed, 77 insertions(+) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index bca20ce26..6a8c01a0f 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -1,5 +1,6 @@ #include "psc.h" +#include "kg/VecRange.hxx" #include "fields.hxx" #include "bnd_fields.hxx" @@ -38,6 +39,7 @@ struct BndFields_ : BndFieldsBase break; } case BND_FLD_OPEN: { + set_background_E_lo(mflds, p, d); break; } default: { @@ -59,6 +61,7 @@ struct BndFields_ : BndFieldsBase break; } case BND_FLD_OPEN: { + set_background_E_hi(mflds, p, d); break; } default: { @@ -91,6 +94,7 @@ struct BndFields_ : BndFieldsBase } case BND_FLD_OPEN: { radiative_H_lo(mflds, p, d); + add_background_H_lo(mflds, p, d); break; } default: { @@ -112,6 +116,7 @@ struct BndFields_ : BndFieldsBase } case BND_FLD_OPEN: { radiative_H_hi(mflds, p, d); + add_background_H_hi(mflds, p, d); break; } default: { @@ -672,6 +677,78 @@ struct BndFields_ : BndFieldsBase } #endif } + + void set_background_E_lo(MfieldsState& mflds, int p, int d) + { + auto F = make_Fields3d(mflds[p]); + Int3 start = mflds.ib(); + Int3 stop = mflds.im(); + stop[d] = 0; + for (Int3 i3 : VecRange(start, stop)) { + F(EX, i3) = background_e[0]; + F(EY, i3) = background_e[1]; + F(EZ, i3) = background_e[2]; + } + } + + void set_background_E_hi(MfieldsState& mflds, int p, int d) + { + auto F = make_Fields3d(mflds[p]); + Int3 start = mflds.ib(); + Int3 stop = mflds.im(); + start[d] = mflds.grid().ldims[d] + 1; + + Int3 neg1 = {0, 0, 0}; + neg1[d] = -1; + + for (Int3 i3 : VecRange(start, stop)) { + F(EX, i3) = background_e[0]; + F(EY, i3) = background_e[1]; + F(EZ, i3) = background_e[2]; + + // the other two components of e are in the domain at this index + F(EX + d, i3 + neg1) = background_e[d]; + } + } + + void add_background_H_lo(MfieldsState& mflds, int p, int d) + { + auto F = make_Fields3d(mflds[p]); + Int3 start = mflds.ib(); + Int3 stop = mflds.im(); + stop[d] = 0; + for (Int3 i3 : VecRange(start, stop)) { + F(HX, i3) += background_h[0]; + F(HY, i3) += background_h[1]; + F(HZ, i3) += background_h[2]; + } + } + + void add_background_H_hi(MfieldsState& mflds, int p, int d) + { + auto F = make_Fields3d(mflds[p]); + Int3 start = mflds.ib(); + Int3 stop = mflds.im(); + start[d] = mflds.grid().ldims[d] + 1; + + Int3 neg1 = {0, 0, 0}; + neg1[d] = -1; + + for (Int3 i3 : VecRange(start, stop)) { + F(HX, i3) += background_h[0]; + F(HY, i3) += background_h[1]; + F(HZ, i3) += background_h[2]; + + // the third component of h is in the domain at this index + int d1 = (d + 1) % 3; + int d2 = (d + 2) % 3; + F(HX + d1, i3 + neg1) += background_h[d1]; + F(HX + d2, i3 + neg1) += background_h[d2]; + } + } + + Vec3 background_e = {0.0, 0.0, 0.0}; + Vec3 background_h = {0.0, 0.0, 0.0}; }; // ====================================================================== From 2a51f971f16764678cf68d82b9578b051627a3fc Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 6 Mar 2026 14:55:05 -0500 Subject: [PATCH 07/36] psc_bnd_fields_impl: rm radiative z impls for now --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 64 ------------------- 1 file changed, 64 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 6a8c01a0f..0b3d2b4ff 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -554,38 +554,6 @@ struct BndFields_ : BndFieldsBase (1.f + dt / dy); } } - } else if (d == 2) { -#ifdef DEBUG - for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - fields_t_set_nan(&F(HX, ix, iy, -1)); - fields_t_set_nan(&F(HX, ix, iy, -2)); - fields_t_set_nan(&F(HY, ix, iy, -1)); - fields_t_set_nan(&F(HY, ix, iy, -2)); - fields_t_set_nan(&F(HZ, ix, iy, -1)); - fields_t_set_nan(&F(HZ, ix, iy, -2)); - } - } -#endif - for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - F(HY, ix, iy, -1) = - (/* + 4.f * C_s_pulse_z1(x+0.5*dx,y,z,t) */ - -2.f * F(EX, ix, iy, 0) - - dt / dy * (F(HZ, ix, iy, 0) - F(HZ, ix, iy - 1, 0)) - - (1.f - dt / dz) * F(HY, ix, iy, 0) + dt * F(JXI, ix, iy, 0)) / - (1.f + dt / dz); - F(HX, ix, iy, -1) = - (/* - 4.f * C_p_pulse_z1(x+0.5*dx,y,z,t) */ - +2.f * F(EY, ix, iy, 0) - /*- dt/dx * (F(HZ, ix,iy,0) - F(HZ, ix-1,iy,0)) FIXME not in yz 2d - */ - - (1.f - dt / dz) * F(HY, ix, iy, 0) - dt * F(JYI, ix, iy, 0)) / - (1.f + dt / dz); - } - } } else { assert(0); } @@ -640,38 +608,6 @@ struct BndFields_ : BndFieldsBase (1.f + dt / dy); } } - } else if (d == 2) { - int mz = ldims[2]; -#ifdef DEBUG - for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = -2; ix < ldims[0] + 2; ix++) { - fields_t_set_nan(&F(HX, ix, iy, mz)); - fields_t_set_nan(&F(HX, ix, iy, mz + 1)); - fields_t_set_nan(&F(HY, ix, iy, mz)); - fields_t_set_nan(&F(HY, ix, iy, mz + 1)); - fields_t_set_nan(&F(HZ, ix, iy, mz + 1)); - } - } -#endif - for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - F(HY, ix, iy, mz) = - (/* - 4.f * C_s_pulse_z2(x+0.5*dx,y,z,t) */ - +2.f * F(EX, ix, iy, mz) + - dt / dy * (F(HZ, ix, iy, mz) - F(HZ, ix, iy - 1, mz)) - - (1.f - dt / dz) * F(HY, ix, iy, mz - 1) - - dt * F(JXI, ix, iy, mz)) / - (1.f + dt / dz); - F(HX, ix, iy, mz) = (/* + 4.f * C_p_pulse_z2(x+0.5*dx,y,z,t) */ - -2.f * F(EY, ix, iy, mz) - /*+ dt/dx * (F(HZ, ix,iy,mz) - F(HZ, ix-1,iy,mz)) - FIXME not in yz 2d*/ - - (1.f - dt / dz) * F(HX, ix, iy, mz - 1) + - dt * F(JYI, ix, iy, mz)) / - (1.f + dt / dz); - } - } } else { assert(0); } From 1e94ede46c9496af016e28fa121e5b2c04621b42 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 6 Mar 2026 14:56:03 -0500 Subject: [PATCH 08/36] psc_bnd_fields_impl: rm old ifdefs --- src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 0b3d2b4ff..99e53df97 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -520,8 +520,6 @@ struct BndFields_ : BndFieldsBase real_t dy = grid.domain.dx[1]; real_t dz = grid.domain.dx[2]; -#if 1 - if (d == 1) { #ifdef DEBUG for (int iz = -2; iz < ldims[2] + 2; iz++) { @@ -536,7 +534,6 @@ struct BndFields_ : BndFieldsBase } } #endif - for (int iz = -2; iz < ldims[2] + 2; iz++) { for (int ix = std::max(-2, ib[0]); ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { @@ -557,7 +554,6 @@ struct BndFields_ : BndFieldsBase } else { assert(0); } -#endif } void radiative_H_hi(MfieldsState& mflds, int p, int d) @@ -572,9 +568,6 @@ struct BndFields_ : BndFieldsBase real_t dy = grid.domain.dx[1]; real_t dz = grid.domain.dx[2]; - // TODO - // assert(0); -#if 1 if (d == 1) { int my _mrc_unused = ldims[1]; #ifdef DEBUG @@ -611,7 +604,6 @@ struct BndFields_ : BndFieldsBase } else { assert(0); } -#endif } void set_background_E_lo(MfieldsState& mflds, int p, int d) From bf8e8a8f72be28a5ed0ed8b9e5c9a408e8931957 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 6 Mar 2026 15:09:26 -0500 Subject: [PATCH 09/36] psc_bnd_fields_impl: use VecRange in radiative_H_lo --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 48 ++++++++++++------- 1 file changed, 32 insertions(+), 16 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 99e53df97..5dd289d02 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -534,23 +534,39 @@ struct BndFields_ : BndFieldsBase } } #endif - for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - F(HX, ix, -1, iz) = - (/* + 4.f * C_s_pulse_y1(x,y,z+0.5*dz,t) */ - -2.f * F(EZ, ix, 0, iz) - /*- dt/dx * (F(HY, ix,0,iz) - F(HY, ix-1,0,iz)) */ - - (1.f - dt / dy) * F(HX, ix, 0, iz) + dt * F(JZI, ix, 0, iz)) / - (1.f + dt / dy); - F(HZ, ix, -1, iz) = - (/* + 4.f * C_p_pulse_y1(x+.5*dx,y,z,t) */ - +2.f * F(EX, ix, 0, iz) - - dt / dz * (F(HY, ix, 0, iz) - F(HY, ix, 0, iz - 1)) - - (1.f - dt / dy) * F(HZ, ix, 0, iz) + dt * F(JXI, ix, 0, iz)) / - (1.f + dt / dy); - } + int d0 = (d + 0) % 3; + int d1 = (d + 1) % 3; + int d2 = (d + 2) % 3; + + Int3 start = mflds.ib(); + Int3 stop = mflds.im(); + start[d0] = -1; + stop[d0] = 0; + + for (Int3 i3 : VecRange(start, stop)) { + Int3 edge_idx = i3; + edge_idx[d0] = 0; + + Int3 neg1_d1{0, 0, 0}; + neg1_d1[d1] = -1; + + Int3 neg1_d2{0, 0, 0}; + neg1_d2[d2] = -1; + + F(HX, i3) = + (/* + 4.f * C_s_pulse_y1(x,y,z+0.5*dz,t) */ + -2.f * F(EZ, edge_idx) + /*- dt/dx * (F(HY, edge_idx) - F(HY, edge_idx + neg1_d2)) */ + - (1.f - dt / dy) * F(HX, edge_idx) + dt * F(JZI, edge_idx)) / + (1.f + dt / dy); + F(HZ, i3) = + (/* + 4.f * C_p_pulse_y1(x+.5*dx,y,z,t) */ + +2.f * F(EX, edge_idx) - + dt / dz * (F(HY, edge_idx) - F(HY, edge_idx + neg1_d1)) - + (1.f - dt / dy) * F(HZ, edge_idx) + dt * F(JXI, edge_idx)) / + (1.f + dt / dy); } + } else { assert(0); } From aa32d74bbef9bfd8f3df3b4e046064e0b755332a Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 6 Mar 2026 15:11:11 -0500 Subject: [PATCH 10/36] psc_bnd_fields_impl: uncomment a grad term --- src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 5dd289d02..0d27a99e5 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -555,9 +555,9 @@ struct BndFields_ : BndFieldsBase F(HX, i3) = (/* + 4.f * C_s_pulse_y1(x,y,z+0.5*dz,t) */ - -2.f * F(EZ, edge_idx) - /*- dt/dx * (F(HY, edge_idx) - F(HY, edge_idx + neg1_d2)) */ - - (1.f - dt / dy) * F(HX, edge_idx) + dt * F(JZI, edge_idx)) / + -2.f * F(EZ, edge_idx) - + dt / dx * (F(HY, edge_idx) - F(HY, edge_idx + neg1_d2)) - + (1.f - dt / dy) * F(HX, edge_idx) + dt * F(JZI, edge_idx)) / (1.f + dt / dy); F(HZ, i3) = (/* + 4.f * C_p_pulse_y1(x+.5*dx,y,z,t) */ From 43d7a265dce08a18a70724780110c8230aa5bf23 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 6 Mar 2026 15:14:45 -0500 Subject: [PATCH 11/36] psc_bnd_fields_impl: don't hardcode comp --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 0d27a99e5..6471d9c97 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -553,17 +553,19 @@ struct BndFields_ : BndFieldsBase Int3 neg1_d2{0, 0, 0}; neg1_d2[d2] = -1; - F(HX, i3) = + F(HX + d2, i3) = (/* + 4.f * C_s_pulse_y1(x,y,z+0.5*dz,t) */ - -2.f * F(EZ, edge_idx) - - dt / dx * (F(HY, edge_idx) - F(HY, edge_idx + neg1_d2)) - - (1.f - dt / dy) * F(HX, edge_idx) + dt * F(JZI, edge_idx)) / + -2.f * F(EX + d1, edge_idx) - + dt / dx * (F(HX + d0, edge_idx) - F(HX + d0, edge_idx + neg1_d2)) - + (1.f - dt / dy) * F(HX + d2, edge_idx) + + dt * F(JXI + d1, edge_idx)) / (1.f + dt / dy); - F(HZ, i3) = + F(HX + d1, i3) = (/* + 4.f * C_p_pulse_y1(x+.5*dx,y,z,t) */ - +2.f * F(EX, edge_idx) - - dt / dz * (F(HY, edge_idx) - F(HY, edge_idx + neg1_d1)) - - (1.f - dt / dy) * F(HZ, edge_idx) + dt * F(JXI, edge_idx)) / + +2.f * F(EX + d2, edge_idx) - + dt / dz * (F(HX + d0, edge_idx) - F(HX + d0, edge_idx + neg1_d1)) - + (1.f - dt / dy) * F(HX + d1, edge_idx) + + dt * F(JXI + d2, edge_idx)) / (1.f + dt / dy); } From db17b81d2d39a286d89e27aa08cf1a1de21cb844 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 6 Mar 2026 15:19:09 -0500 Subject: [PATCH 12/36] psc_bnd_fields_impl: don't hardcode deltas --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 6471d9c97..3bb064646 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -516,9 +516,7 @@ struct BndFields_ : BndFieldsBase Int3 ib = mflds.ib(); Int3 im = mflds.im(); real_t dt = grid.dt; - real_t dx = grid.domain.dx[0]; - real_t dy = grid.domain.dx[1]; - real_t dz = grid.domain.dx[2]; + auto dxi = grid.domain.dx_inv; if (d == 1) { #ifdef DEBUG @@ -556,17 +554,19 @@ struct BndFields_ : BndFieldsBase F(HX + d2, i3) = (/* + 4.f * C_s_pulse_y1(x,y,z+0.5*dz,t) */ -2.f * F(EX + d1, edge_idx) - - dt / dx * (F(HX + d0, edge_idx) - F(HX + d0, edge_idx + neg1_d2)) - - (1.f - dt / dy) * F(HX + d2, edge_idx) + + dt * dxi[d2] * + (F(HX + d0, edge_idx) - F(HX + d0, edge_idx + neg1_d2)) - + (1.f - dt * dxi[d0]) * F(HX + d2, edge_idx) + dt * F(JXI + d1, edge_idx)) / - (1.f + dt / dy); + (1.f + dt * dxi[d0]); F(HX + d1, i3) = (/* + 4.f * C_p_pulse_y1(x+.5*dx,y,z,t) */ +2.f * F(EX + d2, edge_idx) - - dt / dz * (F(HX + d0, edge_idx) - F(HX + d0, edge_idx + neg1_d1)) - - (1.f - dt / dy) * F(HX + d1, edge_idx) + + dt * dxi[d1] * + (F(HX + d0, edge_idx) - F(HX + d0, edge_idx + neg1_d1)) - + (1.f - dt * dxi[d0]) * F(HX + d1, edge_idx) + dt * F(JXI + d2, edge_idx)) / - (1.f + dt / dy); + (1.f + dt * dxi[d0]); } } else { From 276d851470ff6f1cc84080f3a68e60906d5c72a9 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 6 Mar 2026 16:05:29 -0500 Subject: [PATCH 13/36] psc_bnd_fields_impl: +nan_lower --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 87 +++++-------------- 1 file changed, 24 insertions(+), 63 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 3bb064646..777cf7711 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -184,26 +184,33 @@ struct BndFields_ : BndFieldsBase *f = std::numeric_limits::quiet_NaN(); } + void nan_lower(MfieldsState& mflds, int p, int d, int mb, int me) + { +#ifndef DEBUG + return; +#endif + + auto F = make_Fields3d(mflds[p]); + Int3 start = mflds.ib(); + Int3 stop = mflds.im(); + stop[d] = 0; + + for (int m = mb; m < me; m++) { + for (Int3 i3 : VecRange(start, stop)) { + fields_t_set_nan(&F(m, i3)); + } + } + } + void conducting_wall_E_lo(MfieldsState& mflds, int p, int d) { + nan_lower(mflds, p, d, EX, EX + 3); + auto F = make_Fields3d(mflds[p]); const int* ldims = mflds.grid().ldims; Int3 ib = mflds.ib(), im = mflds.im(); if (d == 1) { -#ifdef DEBUG - for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = std::max(-2, F.ib_[0]); - ix < std::min(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { - fields_t_set_nan(&F(EX, ix, -1, iz)); - fields_t_set_nan(&F(EX, ix, -2, iz)); - fields_t_set_nan(&F(EY, ix, -1, iz)); - fields_t_set_nan(&F(EY, ix, -2, iz)); - fields_t_set_nan(&F(EZ, ix, -1, iz)); - fields_t_set_nan(&F(EZ, ix, -2, iz)); - } - } -#endif for (int iz = -2; iz < ldims[2] + 2; iz++) { // FIXME, needs to be for other dir, too, and it's ugly for (int ix = std::max(-2, ib[0]); @@ -218,18 +225,6 @@ struct BndFields_ : BndFieldsBase } } } else if (d == 2) { -#ifdef DEBUG - for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = -2; ix < ldims[0] + 2; ix++) { - fields_t_set_nan(&F(EX, ix, iy, -1)); - fields_t_set_nan(&F(EX, ix, iy, -2)); - fields_t_set_nan(&F(EY, ix, iy, -1)); - fields_t_set_nan(&F(EY, ix, iy, -2)); - fields_t_set_nan(&F(EZ, ix, iy, -1)); - fields_t_set_nan(&F(EZ, ix, iy, -2)); - } - } -#endif for (int iy = -2; iy < ldims[1] + 2; iy++) { for (int ix = std::max(-2, ib[0]); ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { @@ -313,24 +308,13 @@ struct BndFields_ : BndFieldsBase void conducting_wall_H_lo(MfieldsState& mflds, int p, int d) { + nan_lower(mflds, p, d, HX, HX + 3); + auto F = make_Fields3d(mflds[p]); const int* ldims = mflds.grid().ldims; Int3 ib = mflds.ib(), im = mflds.im(); if (d == 1) { -#ifdef DEBUG - for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = std::max(-2, F.ib_[0]); - ix < std::min(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { - fields_t_set_nan(&F(HX, ix, -1, iz)); - fields_t_set_nan(&F(HX, ix, -2, iz)); - fields_t_set_nan(&F(HY, ix, -1, iz)); - fields_t_set_nan(&F(HY, ix, -2, iz)); - fields_t_set_nan(&F(HZ, ix, -1, iz)); - fields_t_set_nan(&F(HZ, ix, -2, iz)); - } - } -#endif for (int iz = -1; iz < ldims[2] + 2; iz++) { for (int ix = std::max(-2, ib[0]); ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { @@ -342,18 +326,6 @@ struct BndFields_ : BndFieldsBase } } } else if (d == 2) { -#ifdef DEBUG - for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = -2; ix < ldims[0] + 2; ix++) { - fields_t_set_nan(&F(HX, ix, iy, -1)); - fields_t_set_nan(&F(HX, ix, iy, -2)); - fields_t_set_nan(&F(HY, ix, iy, -1)); - fields_t_set_nan(&F(HY, ix, iy, -2)); - fields_t_set_nan(&F(HZ, ix, iy, -1)); - fields_t_set_nan(&F(HZ, ix, iy, -2)); - } - } -#endif for (int iy = -2; iy < ldims[1] + 2; iy++) { for (int ix = std::max(-2, ib[0]); ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { @@ -510,6 +482,8 @@ struct BndFields_ : BndFieldsBase void radiative_H_lo(MfieldsState& mflds, int p, int d) { + nan_lower(mflds, p, d, HX, HX + 3); + auto F = make_Fields3d(mflds[p]); const Grid_t& grid = mflds.grid(); Int3 ldims = grid.ldims; @@ -519,19 +493,6 @@ struct BndFields_ : BndFieldsBase auto dxi = grid.domain.dx_inv; if (d == 1) { -#ifdef DEBUG - for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - fields_t_set_nan(&F(HX, ix, -1, iz)); - fields_t_set_nan(&F(HX, ix, -2, iz)); - fields_t_set_nan(&F(HY, ix, -1, iz)); - fields_t_set_nan(&F(HY, ix, -2, iz)); - fields_t_set_nan(&F(HZ, ix, -1, iz)); - fields_t_set_nan(&F(HZ, ix, -2, iz)); - } - } -#endif int d0 = (d + 0) % 3; int d1 = (d + 1) % 3; int d2 = (d + 2) % 3; From a538820d4b6a1eaf50bfacc80cdab0bce3c06728 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 6 Mar 2026 16:11:49 -0500 Subject: [PATCH 14/36] psc_bnd_fields_impl: +Real3 --- src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx | 1 + 1 file changed, 1 insertion(+) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 777cf7711..4a0b1321f 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -16,6 +16,7 @@ struct BndFields_ : BndFieldsBase using Self = BndFields_; using MfieldsState = MFIELDS_STATE; using real_t = typename MfieldsState::real_t; + using Real3 = Vec3; using fields_view_t = typename MfieldsState::fields_view_t; using dim_t = Dim; From bad960ada69932d78357e9aeb1613ee1aa92bec2 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 6 Mar 2026 16:16:53 -0500 Subject: [PATCH 15/36] psc_bnd_fields_impl: +set_lower --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 4a0b1321f..a2d332ff3 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -185,27 +185,32 @@ struct BndFields_ : BndFieldsBase *f = std::numeric_limits::quiet_NaN(); } - void nan_lower(MfieldsState& mflds, int p, int d, int mb, int me) + void nan_lower(MfieldsState& mflds, int p, int d, int mb) { #ifndef DEBUG return; #endif + real_t nan = std::numeric_limits::quiet_NaN(); + set_lower(mflds, p, d, mb, {nan, nan, nan}); + } + void set_lower(MfieldsState& mflds, int p, int d, int mb, Real3 val) + { auto F = make_Fields3d(mflds[p]); Int3 start = mflds.ib(); Int3 stop = mflds.im(); stop[d] = 0; - for (int m = mb; m < me; m++) { + for (int m = mb; m < mb + 3; m++) { for (Int3 i3 : VecRange(start, stop)) { - fields_t_set_nan(&F(m, i3)); + F(m, i3) = val[m - mb]; } } } void conducting_wall_E_lo(MfieldsState& mflds, int p, int d) { - nan_lower(mflds, p, d, EX, EX + 3); + nan_lower(mflds, p, d, EX); auto F = make_Fields3d(mflds[p]); const int* ldims = mflds.grid().ldims; @@ -309,7 +314,7 @@ struct BndFields_ : BndFieldsBase void conducting_wall_H_lo(MfieldsState& mflds, int p, int d) { - nan_lower(mflds, p, d, HX, HX + 3); + nan_lower(mflds, p, d, HX); auto F = make_Fields3d(mflds[p]); const int* ldims = mflds.grid().ldims; @@ -483,7 +488,7 @@ struct BndFields_ : BndFieldsBase void radiative_H_lo(MfieldsState& mflds, int p, int d) { - nan_lower(mflds, p, d, HX, HX + 3); + nan_lower(mflds, p, d, HX); auto F = make_Fields3d(mflds[p]); const Grid_t& grid = mflds.grid(); From 187d1b9e770345b6c32f2332bd00eb3ffce00940 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 6 Mar 2026 16:22:12 -0500 Subject: [PATCH 16/36] psc_bnd_fields_impl: use set_lower for open E --- src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx | 15 +-------------- 1 file changed, 1 insertion(+), 14 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index a2d332ff3..51c5fa1a3 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -40,7 +40,7 @@ struct BndFields_ : BndFieldsBase break; } case BND_FLD_OPEN: { - set_background_E_lo(mflds, p, d); + set_lower(mflds, p, d, EX, background_e); break; } default: { @@ -591,19 +591,6 @@ struct BndFields_ : BndFieldsBase } } - void set_background_E_lo(MfieldsState& mflds, int p, int d) - { - auto F = make_Fields3d(mflds[p]); - Int3 start = mflds.ib(); - Int3 stop = mflds.im(); - stop[d] = 0; - for (Int3 i3 : VecRange(start, stop)) { - F(EX, i3) = background_e[0]; - F(EY, i3) = background_e[1]; - F(EZ, i3) = background_e[2]; - } - } - void set_background_E_hi(MfieldsState& mflds, int p, int d) { auto F = make_Fields3d(mflds[p]); From eed585448d73d7d3f6f3f26558ea67db0c8db7d7 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 6 Mar 2026 16:24:18 -0500 Subject: [PATCH 17/36] psc_bnd_fields_impl: radiative_h_lo is generic --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 74 +++++++++---------- 1 file changed, 33 insertions(+), 41 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 51c5fa1a3..c42c30cb3 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -492,52 +492,44 @@ struct BndFields_ : BndFieldsBase auto F = make_Fields3d(mflds[p]); const Grid_t& grid = mflds.grid(); - Int3 ldims = grid.ldims; - Int3 ib = mflds.ib(); - Int3 im = mflds.im(); real_t dt = grid.dt; auto dxi = grid.domain.dx_inv; - if (d == 1) { - int d0 = (d + 0) % 3; - int d1 = (d + 1) % 3; - int d2 = (d + 2) % 3; - - Int3 start = mflds.ib(); - Int3 stop = mflds.im(); - start[d0] = -1; - stop[d0] = 0; + int d0 = (d + 0) % 3; + int d1 = (d + 1) % 3; + int d2 = (d + 2) % 3; - for (Int3 i3 : VecRange(start, stop)) { - Int3 edge_idx = i3; - edge_idx[d0] = 0; - - Int3 neg1_d1{0, 0, 0}; - neg1_d1[d1] = -1; - - Int3 neg1_d2{0, 0, 0}; - neg1_d2[d2] = -1; - - F(HX + d2, i3) = - (/* + 4.f * C_s_pulse_y1(x,y,z+0.5*dz,t) */ - -2.f * F(EX + d1, edge_idx) - - dt * dxi[d2] * - (F(HX + d0, edge_idx) - F(HX + d0, edge_idx + neg1_d2)) - - (1.f - dt * dxi[d0]) * F(HX + d2, edge_idx) + - dt * F(JXI + d1, edge_idx)) / - (1.f + dt * dxi[d0]); - F(HX + d1, i3) = - (/* + 4.f * C_p_pulse_y1(x+.5*dx,y,z,t) */ - +2.f * F(EX + d2, edge_idx) - - dt * dxi[d1] * - (F(HX + d0, edge_idx) - F(HX + d0, edge_idx + neg1_d1)) - - (1.f - dt * dxi[d0]) * F(HX + d1, edge_idx) + - dt * F(JXI + d2, edge_idx)) / - (1.f + dt * dxi[d0]); - } + Int3 start = mflds.ib(); + Int3 stop = mflds.im(); + start[d0] = -1; + stop[d0] = 0; - } else { - assert(0); + for (Int3 i3 : VecRange(start, stop)) { + Int3 edge_idx = i3; + edge_idx[d0] = 0; + + Int3 neg1_d1{0, 0, 0}; + neg1_d1[d1] = -1; + + Int3 neg1_d2{0, 0, 0}; + neg1_d2[d2] = -1; + + F(HX + d2, i3) = + (/* + 4.f * C_s_pulse_y1(x,y,z+0.5*dz,t), where d0=y */ + -2.f * F(EX + d1, edge_idx) - + dt * dxi[d2] * + (F(HX + d0, edge_idx) - F(HX + d0, edge_idx + neg1_d2)) - + (1.f - dt * dxi[d0]) * F(HX + d2, edge_idx) + + dt * F(JXI + d1, edge_idx)) / + (1.f + dt * dxi[d0]); + F(HX + d1, i3) = + (/* + 4.f * C_p_pulse_y1(x+.5*dx,y,z,t), where d0=y */ + +2.f * F(EX + d2, edge_idx) - + dt * dxi[d1] * + (F(HX + d0, edge_idx) - F(HX + d0, edge_idx + neg1_d1)) - + (1.f - dt * dxi[d0]) * F(HX + d1, edge_idx) + + dt * F(JXI + d2, edge_idx)) / + (1.f + dt * dxi[d0]); } } From fed8412f378edfd22d304fee7e3bdff1a34a890a Mon Sep 17 00:00:00 2001 From: James McClung Date: Sat, 7 Mar 2026 11:18:27 -0500 Subject: [PATCH 18/36] psc_bnd_fields_impl: +nan_upper --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 116 +++++++++--------- 1 file changed, 56 insertions(+), 60 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index c42c30cb3..19dec6830 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -194,6 +194,15 @@ struct BndFields_ : BndFieldsBase set_lower(mflds, p, d, mb, {nan, nan, nan}); } + void nan_upper(MfieldsState& mflds, int p, int d, int mb, bool include_edge) + { +#ifndef DEBUG + return; +#endif + real_t nan = std::numeric_limits::quiet_NaN(); + set_upper(mflds, p, d, mb, {nan, nan, nan}, include_edge); + } + void set_lower(MfieldsState& mflds, int p, int d, int mb, Real3 val) { auto F = make_Fields3d(mflds[p]); @@ -208,6 +217,47 @@ struct BndFields_ : BndFieldsBase } } + /** + * @brief Set E or B upper ghosts to the given constants (each component has + * its own constant). + * @param mflds mflds + * @param p patch index + * @param d which dimension to set the ghosts of + * @param mb `EX` or `HX`; note that `mb+1` and `mb+2` are also set + * @param val the constants + * @param include_edge whether or not values located on exact domain edges + * should be considered "ghosts" + */ + void set_upper(MfieldsState& mflds, int p, int d, int mb, Real3 val, + bool include_edge) + { + auto F = make_Fields3d(mflds[p]); + Int3 start = mflds.ib(); + Int3 stop = mflds.im(); + start[d] = mflds.grid().ldims[d] + 1; + + for (int m = mb; m < mb + 3; m++) { + for (Int3 i3 : VecRange(start, stop)) { + F(m, i3) = val[m - mb]; + } + } + + Int3 edge_start = mflds.ib(); + Int3 edge_stop = mflds.im(); + edge_start[d] = mflds.grid().ldims[d]; + edge_stop[d] = mflds.grid().ldims[d] + 1; + + for (int m = mb; m < mb + 3; m++) { + bool not_edge_ec = mb == EX && m - mb == d; + bool not_edge_fc = mb == HX && m - mb != d; + + if (not_edge_ec || not_edge_fc || include_edge) + for (Int3 i3 : VecRange(edge_start, edge_stop)) { + F(m, i3) = val[m - mb]; + } + } + } + void conducting_wall_E_lo(MfieldsState& mflds, int p, int d) { nan_lower(mflds, p, d, EX); @@ -250,25 +300,14 @@ struct BndFields_ : BndFieldsBase void conducting_wall_E_hi(MfieldsState& mflds, int p, int d) { + nan_upper(mflds, p, d, EX, true); + auto F = make_Fields3d(mflds[p]); const int* ldims = mflds.grid().ldims; Int3 ib = mflds.ib(), im = mflds.im(); if (d == 1) { int my _mrc_unused = ldims[1]; -#ifdef DEBUG - for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = std::max(-2, F.ib_[0]); - ix < std::min(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { - fields_t_set_nan(&F(EX, ix, my, iz)); - fields_t_set_nan(&F(EX, ix, my + 1, iz)); - fields_t_set_nan(&F(EY, ix, my, iz)); - fields_t_set_nan(&F(EY, ix, my + 1, iz)); - fields_t_set_nan(&F(EZ, ix, my, iz)); - fields_t_set_nan(&F(EZ, ix, my + 1, iz)); - } - } -#endif for (int iz = -2; iz < ldims[2] + 2; iz++) { for (int ix = std::max(-2, ib[0]); ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { @@ -283,18 +322,6 @@ struct BndFields_ : BndFieldsBase } } else if (d == 2) { int mz = ldims[2]; -#ifdef DEBUG - for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = -2; ix < ldims[0] + 2; ix++) { - fields_t_set_nan(&F(EX, ix, iy, mz)); - fields_t_set_nan(&F(EX, ix, iy, mz + 1)); - fields_t_set_nan(&F(EY, ix, iy, mz)); - fields_t_set_nan(&F(EY, ix, iy, mz + 1)); - fields_t_set_nan(&F(EZ, ix, iy, mz)); - fields_t_set_nan(&F(EZ, ix, iy, mz + 1)); - } - } -#endif for (int iy = -2; iy < ldims[1] + 2; iy++) { for (int ix = std::max(-2, ib[0]); ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { @@ -349,6 +376,8 @@ struct BndFields_ : BndFieldsBase void conducting_wall_H_hi(MfieldsState& mflds, int p, int d) { + nan_upper(mflds, p, d, HX, false); + auto F = make_Fields3d(mflds[p]); const int* ldims = mflds.grid().ldims; @@ -356,18 +385,6 @@ struct BndFields_ : BndFieldsBase if (d == 1) { int my _mrc_unused = ldims[1]; -#ifdef DEBUG - for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = std::max(-2, F.ib_[0]); - ix < std::min(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { - fields_t_set_nan(&F(HX, ix, my, iz)); - fields_t_set_nan(&F(HX, ix, my + 1, iz)); - fields_t_set_nan(&F(HY, ix, my + 1, iz)); - fields_t_set_nan(&F(HZ, ix, my, iz)); - fields_t_set_nan(&F(HZ, ix, my + 1, iz)); - } - } -#endif for (int iz = -2; iz < ldims[2] + 2; iz++) { for (int ix = std::max(-2, ib[0]); ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { @@ -380,17 +397,6 @@ struct BndFields_ : BndFieldsBase } } else if (d == 2) { int mz = ldims[2]; -#ifdef DEBUG - for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = -2; ix < ldims[0] + 2; ix++) { - fields_t_set_nan(&F(HX, ix, iy, mz)); - fields_t_set_nan(&F(HX, ix, iy, mz + 1)); - fields_t_set_nan(&F(HY, ix, iy, mz)); - fields_t_set_nan(&F(HY, ix, iy, mz + 1)); - fields_t_set_nan(&F(HZ, ix, iy, mz + 1)); - } - } -#endif for (int iy = -2; iy < ldims[1] + 2; iy++) { for (int ix = std::max(-2, ib[0]); ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { @@ -535,6 +541,8 @@ struct BndFields_ : BndFieldsBase void radiative_H_hi(MfieldsState& mflds, int p, int d) { + nan_upper(mflds, p, d, HX, false); + auto F = make_Fields3d(mflds[p]); const Grid_t& grid = mflds.grid(); Int3 ldims = grid.ldims; @@ -547,18 +555,6 @@ struct BndFields_ : BndFieldsBase if (d == 1) { int my _mrc_unused = ldims[1]; -#ifdef DEBUG - for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - fields_t_set_nan(&F(HX, ix, my, iz)); - fields_t_set_nan(&F(HX, ix, my + 1, iz)); - fields_t_set_nan(&F(HY, ix, my, iz)); - fields_t_set_nan(&F(HY, ix, my + 1, iz)); - fields_t_set_nan(&F(HZ, ix, my + 1, iz)); - } - } -#endif for (int iz = -2; iz < ldims[2] + 2; iz++) { for (int ix = std::max(-2, ib[0]); ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { From 0672079ea2fa74208ca539480ce903ceb87b4ffd Mon Sep 17 00:00:00 2001 From: James McClung Date: Sat, 7 Mar 2026 11:29:52 -0500 Subject: [PATCH 19/36] psc_bnd_fields_impl: +include_edge to set_lower --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 46 ++++++++++++++++--- 1 file changed, 39 insertions(+), 7 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 19dec6830..e960b083d 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -40,7 +40,7 @@ struct BndFields_ : BndFieldsBase break; } case BND_FLD_OPEN: { - set_lower(mflds, p, d, EX, background_e); + set_lower(mflds, p, d, EX, background_e, false); break; } default: { @@ -185,13 +185,13 @@ struct BndFields_ : BndFieldsBase *f = std::numeric_limits::quiet_NaN(); } - void nan_lower(MfieldsState& mflds, int p, int d, int mb) + void nan_lower(MfieldsState& mflds, int p, int d, int mb, bool include_edge) { #ifndef DEBUG return; #endif real_t nan = std::numeric_limits::quiet_NaN(); - set_lower(mflds, p, d, mb, {nan, nan, nan}); + set_lower(mflds, p, d, mb, {nan, nan, nan}, include_edge); } void nan_upper(MfieldsState& mflds, int p, int d, int mb, bool include_edge) @@ -203,7 +203,19 @@ struct BndFields_ : BndFieldsBase set_upper(mflds, p, d, mb, {nan, nan, nan}, include_edge); } - void set_lower(MfieldsState& mflds, int p, int d, int mb, Real3 val) + /** + * @brief Set E or B lower ghosts to the given constants (each component has + * its own constant). + * @param mflds mflds + * @param p patch index + * @param d which dimension to set the ghosts of + * @param mb `EX` or `HX`; note that `mb+1` and `mb+2` are also set + * @param val the constants + * @param include_edge whether or not values located on exact domain edges + * should be considered "ghosts" + */ + void set_lower(MfieldsState& mflds, int p, int d, int mb, Real3 val, + bool include_edge) { auto F = make_Fields3d(mflds[p]); Int3 start = mflds.ib(); @@ -215,6 +227,26 @@ struct BndFields_ : BndFieldsBase F(m, i3) = val[m - mb]; } } + + if (!include_edge) { + return; + } + + Int3 edge_start = mflds.ib(); + Int3 edge_stop = mflds.im(); + edge_start[d] = 0; + edge_stop[d] = 1; + + for (int m = mb; m < mb + 3; m++) { + bool edge_ec = mb == EX && m - mb != d; + bool edge_fc = mb == HX && m - mb == d; + + if (edge_ec || edge_fc) { + for (Int3 i3 : VecRange(edge_start, edge_stop)) { + F(m, i3) = val[m - mb]; + } + } + } } /** @@ -260,7 +292,7 @@ struct BndFields_ : BndFieldsBase void conducting_wall_E_lo(MfieldsState& mflds, int p, int d) { - nan_lower(mflds, p, d, EX); + nan_lower(mflds, p, d, EX, true); auto F = make_Fields3d(mflds[p]); const int* ldims = mflds.grid().ldims; @@ -341,7 +373,7 @@ struct BndFields_ : BndFieldsBase void conducting_wall_H_lo(MfieldsState& mflds, int p, int d) { - nan_lower(mflds, p, d, HX); + nan_lower(mflds, p, d, HX, false); auto F = make_Fields3d(mflds[p]); const int* ldims = mflds.grid().ldims; @@ -494,7 +526,7 @@ struct BndFields_ : BndFieldsBase void radiative_H_lo(MfieldsState& mflds, int p, int d) { - nan_lower(mflds, p, d, HX); + nan_lower(mflds, p, d, HX, false); auto F = make_Fields3d(mflds[p]); const Grid_t& grid = mflds.grid(); From f7af23bc0688740a8f35c1e5a56c971aa68a5efd Mon Sep 17 00:00:00 2001 From: James McClung Date: Sat, 7 Mar 2026 11:31:55 -0500 Subject: [PATCH 20/36] psc_bnd_fields_impl: -set_background_E_hi --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 22 +------------------ 1 file changed, 1 insertion(+), 21 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index e960b083d..309f9b9ba 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -62,7 +62,7 @@ struct BndFields_ : BndFieldsBase break; } case BND_FLD_OPEN: { - set_background_E_hi(mflds, p, d); + set_upper(mflds, p, d, EX, background_e, false); break; } default: { @@ -611,26 +611,6 @@ struct BndFields_ : BndFieldsBase } } - void set_background_E_hi(MfieldsState& mflds, int p, int d) - { - auto F = make_Fields3d(mflds[p]); - Int3 start = mflds.ib(); - Int3 stop = mflds.im(); - start[d] = mflds.grid().ldims[d] + 1; - - Int3 neg1 = {0, 0, 0}; - neg1[d] = -1; - - for (Int3 i3 : VecRange(start, stop)) { - F(EX, i3) = background_e[0]; - F(EY, i3) = background_e[1]; - F(EZ, i3) = background_e[2]; - - // the other two components of e are in the domain at this index - F(EX + d, i3 + neg1) = background_e[d]; - } - } - void add_background_H_lo(MfieldsState& mflds, int p, int d) { auto F = make_Fields3d(mflds[p]); From 1d039b778d407a9ed6ef99864cfd96db20fae073 Mon Sep 17 00:00:00 2001 From: James McClung Date: Sat, 7 Mar 2026 11:32:50 -0500 Subject: [PATCH 21/36] psc_bnd_fields_impl: +TODOs about gtensor views --- src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 309f9b9ba..ee949ae36 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -222,6 +222,8 @@ struct BndFields_ : BndFieldsBase Int3 stop = mflds.im(); stop[d] = 0; + // TODO use gtensor views instead of VecRange + for (int m = mb; m < mb + 3; m++) { for (Int3 i3 : VecRange(start, stop)) { F(m, i3) = val[m - mb]; @@ -268,6 +270,8 @@ struct BndFields_ : BndFieldsBase Int3 stop = mflds.im(); start[d] = mflds.grid().ldims[d] + 1; + // TODO use gtensor views instead of VecRange + for (int m = mb; m < mb + 3; m++) { for (Int3 i3 : VecRange(start, stop)) { F(m, i3) = val[m - mb]; From 9b46b64856aa48231dc25a600e8de31fa8411b6a Mon Sep 17 00:00:00 2001 From: James McClung Date: Sat, 7 Mar 2026 11:34:26 -0500 Subject: [PATCH 22/36] psc_bnd_fields_impl: rename new setters --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 34 ++++++++++--------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index ee949ae36..80b871549 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -40,7 +40,7 @@ struct BndFields_ : BndFieldsBase break; } case BND_FLD_OPEN: { - set_lower(mflds, p, d, EX, background_e, false); + set_lower_ghosts(mflds, p, d, EX, background_e, false); break; } default: { @@ -62,7 +62,7 @@ struct BndFields_ : BndFieldsBase break; } case BND_FLD_OPEN: { - set_upper(mflds, p, d, EX, background_e, false); + set_upper_ghosts(mflds, p, d, EX, background_e, false); break; } default: { @@ -185,22 +185,24 @@ struct BndFields_ : BndFieldsBase *f = std::numeric_limits::quiet_NaN(); } - void nan_lower(MfieldsState& mflds, int p, int d, int mb, bool include_edge) + void set_lower_ghosts_to_nan(MfieldsState& mflds, int p, int d, int mb, + bool include_edge) { #ifndef DEBUG return; #endif real_t nan = std::numeric_limits::quiet_NaN(); - set_lower(mflds, p, d, mb, {nan, nan, nan}, include_edge); + set_lower_ghosts(mflds, p, d, mb, {nan, nan, nan}, include_edge); } - void nan_upper(MfieldsState& mflds, int p, int d, int mb, bool include_edge) + void set_upper_ghosts_to_nan(MfieldsState& mflds, int p, int d, int mb, + bool include_edge) { #ifndef DEBUG return; #endif real_t nan = std::numeric_limits::quiet_NaN(); - set_upper(mflds, p, d, mb, {nan, nan, nan}, include_edge); + set_upper_ghosts(mflds, p, d, mb, {nan, nan, nan}, include_edge); } /** @@ -214,8 +216,8 @@ struct BndFields_ : BndFieldsBase * @param include_edge whether or not values located on exact domain edges * should be considered "ghosts" */ - void set_lower(MfieldsState& mflds, int p, int d, int mb, Real3 val, - bool include_edge) + void set_lower_ghosts(MfieldsState& mflds, int p, int d, int mb, Real3 val, + bool include_edge) { auto F = make_Fields3d(mflds[p]); Int3 start = mflds.ib(); @@ -262,8 +264,8 @@ struct BndFields_ : BndFieldsBase * @param include_edge whether or not values located on exact domain edges * should be considered "ghosts" */ - void set_upper(MfieldsState& mflds, int p, int d, int mb, Real3 val, - bool include_edge) + void set_upper_ghosts(MfieldsState& mflds, int p, int d, int mb, Real3 val, + bool include_edge) { auto F = make_Fields3d(mflds[p]); Int3 start = mflds.ib(); @@ -296,7 +298,7 @@ struct BndFields_ : BndFieldsBase void conducting_wall_E_lo(MfieldsState& mflds, int p, int d) { - nan_lower(mflds, p, d, EX, true); + set_lower_ghosts_to_nan(mflds, p, d, EX, true); auto F = make_Fields3d(mflds[p]); const int* ldims = mflds.grid().ldims; @@ -336,7 +338,7 @@ struct BndFields_ : BndFieldsBase void conducting_wall_E_hi(MfieldsState& mflds, int p, int d) { - nan_upper(mflds, p, d, EX, true); + set_upper_ghosts_to_nan(mflds, p, d, EX, true); auto F = make_Fields3d(mflds[p]); const int* ldims = mflds.grid().ldims; @@ -377,7 +379,7 @@ struct BndFields_ : BndFieldsBase void conducting_wall_H_lo(MfieldsState& mflds, int p, int d) { - nan_lower(mflds, p, d, HX, false); + set_lower_ghosts_to_nan(mflds, p, d, HX, false); auto F = make_Fields3d(mflds[p]); const int* ldims = mflds.grid().ldims; @@ -412,7 +414,7 @@ struct BndFields_ : BndFieldsBase void conducting_wall_H_hi(MfieldsState& mflds, int p, int d) { - nan_upper(mflds, p, d, HX, false); + set_upper_ghosts_to_nan(mflds, p, d, HX, false); auto F = make_Fields3d(mflds[p]); @@ -530,7 +532,7 @@ struct BndFields_ : BndFieldsBase void radiative_H_lo(MfieldsState& mflds, int p, int d) { - nan_lower(mflds, p, d, HX, false); + set_lower_ghosts_to_nan(mflds, p, d, HX, false); auto F = make_Fields3d(mflds[p]); const Grid_t& grid = mflds.grid(); @@ -577,7 +579,7 @@ struct BndFields_ : BndFieldsBase void radiative_H_hi(MfieldsState& mflds, int p, int d) { - nan_upper(mflds, p, d, HX, false); + set_upper_ghosts_to_nan(mflds, p, d, HX, false); auto F = make_Fields3d(mflds[p]); const Grid_t& grid = mflds.grid(); From f8e6276fcbfa48439906f3e68467deba810efbd6 Mon Sep 17 00:00:00 2001 From: James McClung Date: Sat, 7 Mar 2026 11:44:35 -0500 Subject: [PATCH 23/36] psc_bnd_fields_impl: use VecRange in radiative hi --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 49 +++++++++++-------- 1 file changed, 29 insertions(+), 20 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 80b871549..a40f10b7f 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -584,33 +584,42 @@ struct BndFields_ : BndFieldsBase auto F = make_Fields3d(mflds[p]); const Grid_t& grid = mflds.grid(); Int3 ldims = grid.ldims; - Int3 ib = mflds.ib(); - Int3 im = mflds.im(); real_t dt = grid.dt; real_t dx = grid.domain.dx[0]; real_t dy = grid.domain.dx[1]; real_t dz = grid.domain.dx[2]; + int d0 = (d + 0) % 3; + int d1 = (d + 1) % 3; + int d2 = (d + 2) % 3; + + Int3 start = mflds.ib(); + Int3 stop = mflds.im(); + start[d0] = grid.ldims[d0]; + stop[d0] = grid.ldims[d0] + 1; + if (d == 1) { int my _mrc_unused = ldims[1]; - for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - F(HX, ix, my, iz) = - (/* + 4.f * C_s_pulse_y2(x,y,z+0.5*dz,t) */ - +2.f * F(EZ, ix, my, iz) - /*+ dt/dx * (F(HY, ix,my,iz) - F(HY, ix-1,my,iz)) */ - - (1.f - dt / dy) * F(HX, ix, my - 1, iz) - - dt * F(JZI, ix, my, iz)) / - (1.f + dt / dy); - F(HZ, ix, my, iz) = - (/* + 4.f * C_p_pulse_y2(x+.5*dx,y,z,t) */ - -2.f * F(EX, ix, my, iz) + - dt / dz * (F(HY, ix, my, iz) - F(HY, ix, my, iz - 1)) - - (1.f - dt / dy) * F(HZ, ix, my - 1, iz) - - dt * F(JXI, ix, my, iz)) / - (1.f + dt / dy); - } + for (Int3 i3 : VecRange(start, stop)) { + Int3 edge_idx = i3; + edge_idx[d0] -= 1; + + Int3 neg1_d1{0, 0, 0}; + neg1_d1[d1] = -1; + + Int3 neg1_d2{0, 0, 0}; + neg1_d2[d2] = -1; + + F(HX, i3) = + (/* + 4.f * C_s_pulse_y2(x,y,z+0.5*dz,t) */ + +2.f * F(EZ, i3) + dt / dx * (F(HY, i3) - F(HY, i3 + neg1_d2)) - + (1.f - dt / dy) * F(HX, edge_idx) - dt * F(JZI, i3)) / + (1.f + dt / dy); + F(HZ, i3) = + (/* + 4.f * C_p_pulse_y2(x+.5*dx,y,z,t) */ + -2.f * F(EX, i3) + dt / dz * (F(HY, i3) - F(HY, i3 + neg1_d1)) - + (1.f - dt / dy) * F(HZ, edge_idx) - dt * F(JXI, i3)) / + (1.f + dt / dy); } } else { assert(0); From 7f90ea847a947aa7c21f2744277ba4d626b9cd50 Mon Sep 17 00:00:00 2001 From: James McClung Date: Sat, 7 Mar 2026 11:47:47 -0500 Subject: [PATCH 24/36] psc_bnd_fields_impl: comp from d --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index a40f10b7f..d9e9e226a 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -610,15 +610,17 @@ struct BndFields_ : BndFieldsBase Int3 neg1_d2{0, 0, 0}; neg1_d2[d2] = -1; - F(HX, i3) = - (/* + 4.f * C_s_pulse_y2(x,y,z+0.5*dz,t) */ - +2.f * F(EZ, i3) + dt / dx * (F(HY, i3) - F(HY, i3 + neg1_d2)) - - (1.f - dt / dy) * F(HX, edge_idx) - dt * F(JZI, i3)) / + F(HX + d2, i3) = + (/* + 4.f * C_s_pulse_y2(x,y,z+0.5*dz,t), where d0=y */ + +2.f * F(EX + d1, i3) + + dt / dx * (F(HX + d0, i3) - F(HX + d0, i3 + neg1_d2)) - + (1.f - dt / dy) * F(HX + d2, edge_idx) - dt * F(JXI + d1, i3)) / (1.f + dt / dy); - F(HZ, i3) = - (/* + 4.f * C_p_pulse_y2(x+.5*dx,y,z,t) */ - -2.f * F(EX, i3) + dt / dz * (F(HY, i3) - F(HY, i3 + neg1_d1)) - - (1.f - dt / dy) * F(HZ, edge_idx) - dt * F(JXI, i3)) / + F(HX + d1, i3) = + (/* + 4.f * C_p_pulse_y2(x+.5*dx,y,z,t), where d0=y */ + -2.f * F(EX + d2, i3) + + dt / dz * (F(HX + d0, i3) - F(HX + d0, i3 + neg1_d1)) - + (1.f - dt / dy) * F(HX + d1, edge_idx) - dt * F(JXI + d2, i3)) / (1.f + dt / dy); } } else { From 6796ed6721fb4a3fafdb9c4603c3688cd8417a88 Mon Sep 17 00:00:00 2001 From: James McClung Date: Sat, 7 Mar 2026 11:49:01 -0500 Subject: [PATCH 25/36] psc_bnd_fields_impl: delta from d --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index d9e9e226a..a2c1f1ffd 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -585,9 +585,7 @@ struct BndFields_ : BndFieldsBase const Grid_t& grid = mflds.grid(); Int3 ldims = grid.ldims; real_t dt = grid.dt; - real_t dx = grid.domain.dx[0]; - real_t dy = grid.domain.dx[1]; - real_t dz = grid.domain.dx[2]; + auto dxi = grid.domain.dx_inv; int d0 = (d + 0) % 3; int d1 = (d + 1) % 3; @@ -613,15 +611,15 @@ struct BndFields_ : BndFieldsBase F(HX + d2, i3) = (/* + 4.f * C_s_pulse_y2(x,y,z+0.5*dz,t), where d0=y */ +2.f * F(EX + d1, i3) + - dt / dx * (F(HX + d0, i3) - F(HX + d0, i3 + neg1_d2)) - - (1.f - dt / dy) * F(HX + d2, edge_idx) - dt * F(JXI + d1, i3)) / - (1.f + dt / dy); + dt * dxi[d2] * (F(HX + d0, i3) - F(HX + d0, i3 + neg1_d2)) - + (1.f - dt * dxi[d0]) * F(HX + d2, edge_idx) - dt * F(JXI + d1, i3)) / + (1.f + dt * dxi[d0]); F(HX + d1, i3) = (/* + 4.f * C_p_pulse_y2(x+.5*dx,y,z,t), where d0=y */ -2.f * F(EX + d2, i3) + - dt / dz * (F(HX + d0, i3) - F(HX + d0, i3 + neg1_d1)) - - (1.f - dt / dy) * F(HX + d1, edge_idx) - dt * F(JXI + d2, i3)) / - (1.f + dt / dy); + dt * dxi[d1] * (F(HX + d0, i3) - F(HX + d0, i3 + neg1_d1)) - + (1.f - dt * dxi[d0]) * F(HX + d1, edge_idx) - dt * F(JXI + d2, i3)) / + (1.f + dt * dxi[d0]); } } else { assert(0); From decbdad234f073b8b09e828bbe3b6207ee37a505 Mon Sep 17 00:00:00 2001 From: James McClung Date: Sat, 7 Mar 2026 11:50:14 -0500 Subject: [PATCH 26/36] psc_bnd_fields_impl: upper radiative is generic --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 49 +++++++++---------- 1 file changed, 22 insertions(+), 27 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index a2c1f1ffd..059256385 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -596,33 +596,28 @@ struct BndFields_ : BndFieldsBase start[d0] = grid.ldims[d0]; stop[d0] = grid.ldims[d0] + 1; - if (d == 1) { - int my _mrc_unused = ldims[1]; - for (Int3 i3 : VecRange(start, stop)) { - Int3 edge_idx = i3; - edge_idx[d0] -= 1; - - Int3 neg1_d1{0, 0, 0}; - neg1_d1[d1] = -1; - - Int3 neg1_d2{0, 0, 0}; - neg1_d2[d2] = -1; - - F(HX + d2, i3) = - (/* + 4.f * C_s_pulse_y2(x,y,z+0.5*dz,t), where d0=y */ - +2.f * F(EX + d1, i3) + - dt * dxi[d2] * (F(HX + d0, i3) - F(HX + d0, i3 + neg1_d2)) - - (1.f - dt * dxi[d0]) * F(HX + d2, edge_idx) - dt * F(JXI + d1, i3)) / - (1.f + dt * dxi[d0]); - F(HX + d1, i3) = - (/* + 4.f * C_p_pulse_y2(x+.5*dx,y,z,t), where d0=y */ - -2.f * F(EX + d2, i3) + - dt * dxi[d1] * (F(HX + d0, i3) - F(HX + d0, i3 + neg1_d1)) - - (1.f - dt * dxi[d0]) * F(HX + d1, edge_idx) - dt * F(JXI + d2, i3)) / - (1.f + dt * dxi[d0]); - } - } else { - assert(0); + for (Int3 i3 : VecRange(start, stop)) { + Int3 edge_idx = i3; + edge_idx[d0] -= 1; + + Int3 neg1_d1{0, 0, 0}; + neg1_d1[d1] = -1; + + Int3 neg1_d2{0, 0, 0}; + neg1_d2[d2] = -1; + + F(HX + d2, i3) = + (/* + 4.f * C_s_pulse_y2(x,y,z+0.5*dz,t), where d0=y */ + +2.f * F(EX + d1, i3) + + dt * dxi[d2] * (F(HX + d0, i3) - F(HX + d0, i3 + neg1_d2)) - + (1.f - dt * dxi[d0]) * F(HX + d2, edge_idx) - dt * F(JXI + d1, i3)) / + (1.f + dt * dxi[d0]); + F(HX + d1, i3) = + (/* + 4.f * C_p_pulse_y2(x+.5*dx,y,z,t), where d0=y */ + -2.f * F(EX + d2, i3) + + dt * dxi[d1] * (F(HX + d0, i3) - F(HX + d0, i3 + neg1_d1)) - + (1.f - dt * dxi[d0]) * F(HX + d1, edge_idx) - dt * F(JXI + d2, i3)) / + (1.f + dt * dxi[d0]); } } From cfae91cabc1528fd41280927c0e7001d41115baf Mon Sep 17 00:00:00 2001 From: James McClung Date: Sat, 7 Mar 2026 12:01:17 -0500 Subject: [PATCH 27/36] psc_bnd_fields_impl: +H0, H1, etc. --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 48 +++++++++---------- 1 file changed, 22 insertions(+), 26 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 059256385..eb508f2aa 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -539,9 +539,10 @@ struct BndFields_ : BndFieldsBase real_t dt = grid.dt; auto dxi = grid.domain.dx_inv; - int d0 = (d + 0) % 3; - int d1 = (d + 1) % 3; - int d2 = (d + 2) % 3; + int d0 = (d + 0) % 3, d1 = (d + 1) % 3, d2 = (d + 2) % 3; + int H0 = HX + d0, H1 = HX + d1, H2 = HX + d2; + int E1 = EX + d1, E2 = EX + d2; + int J1 = JXI + d1, J2 = JXI + d2; Int3 start = mflds.ib(); Int3 stop = mflds.im(); @@ -558,21 +559,17 @@ struct BndFields_ : BndFieldsBase Int3 neg1_d2{0, 0, 0}; neg1_d2[d2] = -1; - F(HX + d2, i3) = + F(H2, i3) = (/* + 4.f * C_s_pulse_y1(x,y,z+0.5*dz,t), where d0=y */ - -2.f * F(EX + d1, edge_idx) - - dt * dxi[d2] * - (F(HX + d0, edge_idx) - F(HX + d0, edge_idx + neg1_d2)) - - (1.f - dt * dxi[d0]) * F(HX + d2, edge_idx) + - dt * F(JXI + d1, edge_idx)) / + -2.f * F(E1, edge_idx) - + dt * dxi[d2] * (F(H0, edge_idx) - F(H0, edge_idx + neg1_d2)) - + (1.f - dt * dxi[d0]) * F(H2, edge_idx) + dt * F(J1, edge_idx)) / (1.f + dt * dxi[d0]); - F(HX + d1, i3) = + F(H1, i3) = (/* + 4.f * C_p_pulse_y1(x+.5*dx,y,z,t), where d0=y */ - +2.f * F(EX + d2, edge_idx) - - dt * dxi[d1] * - (F(HX + d0, edge_idx) - F(HX + d0, edge_idx + neg1_d1)) - - (1.f - dt * dxi[d0]) * F(HX + d1, edge_idx) + - dt * F(JXI + d2, edge_idx)) / + +2.f * F(E2, edge_idx) - + dt * dxi[d1] * (F(H0, edge_idx) - F(H0, edge_idx + neg1_d1)) - + (1.f - dt * dxi[d0]) * F(H1, edge_idx) + dt * F(J2, edge_idx)) / (1.f + dt * dxi[d0]); } } @@ -587,9 +584,10 @@ struct BndFields_ : BndFieldsBase real_t dt = grid.dt; auto dxi = grid.domain.dx_inv; - int d0 = (d + 0) % 3; - int d1 = (d + 1) % 3; - int d2 = (d + 2) % 3; + int d0 = (d + 0) % 3, d1 = (d + 1) % 3, d2 = (d + 2) % 3; + int H0 = HX + d0, H1 = HX + d1, H2 = HX + d2; + int E1 = EX + d1, E2 = EX + d2; + int J1 = JXI + d1, J2 = JXI + d2; Int3 start = mflds.ib(); Int3 stop = mflds.im(); @@ -606,17 +604,15 @@ struct BndFields_ : BndFieldsBase Int3 neg1_d2{0, 0, 0}; neg1_d2[d2] = -1; - F(HX + d2, i3) = + F(H2, i3) = (/* + 4.f * C_s_pulse_y2(x,y,z+0.5*dz,t), where d0=y */ - +2.f * F(EX + d1, i3) + - dt * dxi[d2] * (F(HX + d0, i3) - F(HX + d0, i3 + neg1_d2)) - - (1.f - dt * dxi[d0]) * F(HX + d2, edge_idx) - dt * F(JXI + d1, i3)) / + +2.f * F(E1, i3) + dt * dxi[d2] * (F(H0, i3) - F(H0, i3 + neg1_d2)) - + (1.f - dt * dxi[d0]) * F(H2, edge_idx) - dt * F(J1, i3)) / (1.f + dt * dxi[d0]); - F(HX + d1, i3) = + F(H1, i3) = (/* + 4.f * C_p_pulse_y2(x+.5*dx,y,z,t), where d0=y */ - -2.f * F(EX + d2, i3) + - dt * dxi[d1] * (F(HX + d0, i3) - F(HX + d0, i3 + neg1_d1)) - - (1.f - dt * dxi[d0]) * F(HX + d1, edge_idx) - dt * F(JXI + d2, i3)) / + -2.f * F(E2, i3) + dt * dxi[d1] * (F(H0, i3) - F(H0, i3 + neg1_d1)) - + (1.f - dt * dxi[d0]) * F(H1, edge_idx) - dt * F(J2, i3)) / (1.f + dt * dxi[d0]); } } From e1099d63cf8f9619fabf124e6851643060df5121 Mon Sep 17 00:00:00 2001 From: James McClung Date: Sat, 7 Mar 2026 12:05:03 -0500 Subject: [PATCH 28/36] psc_bnd_fields_impl: +dtdx --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 38 +++++++++---------- 1 file changed, 18 insertions(+), 20 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index eb508f2aa..b29219958 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -537,7 +537,7 @@ struct BndFields_ : BndFieldsBase auto F = make_Fields3d(mflds[p]); const Grid_t& grid = mflds.grid(); real_t dt = grid.dt; - auto dxi = grid.domain.dx_inv; + Real3 dtdx = dt * Real3(grid.domain.dx_inv); int d0 = (d + 0) % 3, d1 = (d + 1) % 3, d2 = (d + 2) % 3; int H0 = HX + d0, H1 = HX + d1, H2 = HX + d2; @@ -559,18 +559,16 @@ struct BndFields_ : BndFieldsBase Int3 neg1_d2{0, 0, 0}; neg1_d2[d2] = -1; - F(H2, i3) = - (/* + 4.f * C_s_pulse_y1(x,y,z+0.5*dz,t), where d0=y */ - -2.f * F(E1, edge_idx) - - dt * dxi[d2] * (F(H0, edge_idx) - F(H0, edge_idx + neg1_d2)) - - (1.f - dt * dxi[d0]) * F(H2, edge_idx) + dt * F(J1, edge_idx)) / - (1.f + dt * dxi[d0]); - F(H1, i3) = - (/* + 4.f * C_p_pulse_y1(x+.5*dx,y,z,t), where d0=y */ - +2.f * F(E2, edge_idx) - - dt * dxi[d1] * (F(H0, edge_idx) - F(H0, edge_idx + neg1_d1)) - - (1.f - dt * dxi[d0]) * F(H1, edge_idx) + dt * F(J2, edge_idx)) / - (1.f + dt * dxi[d0]); + F(H2, i3) = (/* + 4.f * C_s_pulse_y1(x,y,z+0.5*dz,t), where d0=y */ + -2.f * F(E1, edge_idx) - + dtdx[d2] * (F(H0, edge_idx) - F(H0, edge_idx + neg1_d2)) - + (1.f - dtdx[d0]) * F(H2, edge_idx) + dt * F(J1, edge_idx)) / + (1.f + dtdx[d0]); + F(H1, i3) = (/* + 4.f * C_p_pulse_y1(x+.5*dx,y,z,t), where d0=y */ + +2.f * F(E2, edge_idx) - + dtdx[d1] * (F(H0, edge_idx) - F(H0, edge_idx + neg1_d1)) - + (1.f - dtdx[d0]) * F(H1, edge_idx) + dt * F(J2, edge_idx)) / + (1.f + dtdx[d0]); } } @@ -582,7 +580,7 @@ struct BndFields_ : BndFieldsBase const Grid_t& grid = mflds.grid(); Int3 ldims = grid.ldims; real_t dt = grid.dt; - auto dxi = grid.domain.dx_inv; + Real3 dtdx = dt * Real3(grid.domain.dx_inv); int d0 = (d + 0) % 3, d1 = (d + 1) % 3, d2 = (d + 2) % 3; int H0 = HX + d0, H1 = HX + d1, H2 = HX + d2; @@ -606,14 +604,14 @@ struct BndFields_ : BndFieldsBase F(H2, i3) = (/* + 4.f * C_s_pulse_y2(x,y,z+0.5*dz,t), where d0=y */ - +2.f * F(E1, i3) + dt * dxi[d2] * (F(H0, i3) - F(H0, i3 + neg1_d2)) - - (1.f - dt * dxi[d0]) * F(H2, edge_idx) - dt * F(J1, i3)) / - (1.f + dt * dxi[d0]); + +2.f * F(E1, i3) + dtdx[d2] * (F(H0, i3) - F(H0, i3 + neg1_d2)) - + (1.f - dtdx[d0]) * F(H2, edge_idx) - dt * F(J1, i3)) / + (1.f + dtdx[d0]); F(H1, i3) = (/* + 4.f * C_p_pulse_y2(x+.5*dx,y,z,t), where d0=y */ - -2.f * F(E2, i3) + dt * dxi[d1] * (F(H0, i3) - F(H0, i3 + neg1_d1)) - - (1.f - dt * dxi[d0]) * F(H1, edge_idx) - dt * F(J2, i3)) / - (1.f + dt * dxi[d0]); + -2.f * F(E2, i3) + dtdx[d1] * (F(H0, i3) - F(H0, i3 + neg1_d1)) - + (1.f - dtdx[d0]) * F(H1, edge_idx) - dt * F(J2, i3)) / + (1.f + dtdx[d0]); } } From cac9d9d528372127ae874316b2fd4e469eda6fb8 Mon Sep 17 00:00:00 2001 From: James McClung Date: Sat, 7 Mar 2026 12:06:54 -0500 Subject: [PATCH 29/36] psc_bnd_fields_impl: -add_background_H --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 38 ------------------- 1 file changed, 38 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index b29219958..efc088fd4 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -95,7 +95,6 @@ struct BndFields_ : BndFieldsBase } case BND_FLD_OPEN: { radiative_H_lo(mflds, p, d); - add_background_H_lo(mflds, p, d); break; } default: { @@ -117,7 +116,6 @@ struct BndFields_ : BndFieldsBase } case BND_FLD_OPEN: { radiative_H_hi(mflds, p, d); - add_background_H_hi(mflds, p, d); break; } default: { @@ -615,42 +613,6 @@ struct BndFields_ : BndFieldsBase } } - void add_background_H_lo(MfieldsState& mflds, int p, int d) - { - auto F = make_Fields3d(mflds[p]); - Int3 start = mflds.ib(); - Int3 stop = mflds.im(); - stop[d] = 0; - for (Int3 i3 : VecRange(start, stop)) { - F(HX, i3) += background_h[0]; - F(HY, i3) += background_h[1]; - F(HZ, i3) += background_h[2]; - } - } - - void add_background_H_hi(MfieldsState& mflds, int p, int d) - { - auto F = make_Fields3d(mflds[p]); - Int3 start = mflds.ib(); - Int3 stop = mflds.im(); - start[d] = mflds.grid().ldims[d] + 1; - - Int3 neg1 = {0, 0, 0}; - neg1[d] = -1; - - for (Int3 i3 : VecRange(start, stop)) { - F(HX, i3) += background_h[0]; - F(HY, i3) += background_h[1]; - F(HZ, i3) += background_h[2]; - - // the third component of h is in the domain at this index - int d1 = (d + 1) % 3; - int d2 = (d + 2) % 3; - F(HX + d1, i3 + neg1) += background_h[d1]; - F(HX + d2, i3 + neg1) += background_h[d2]; - } - } - Vec3 background_e = {0.0, 0.0, 0.0}; Vec3 background_h = {0.0, 0.0, 0.0}; }; From 6e50ef65a121e86ec2921f42d69a48ec03309428 Mon Sep 17 00:00:00 2001 From: James McClung Date: Sat, 7 Mar 2026 12:09:48 -0500 Subject: [PATCH 30/36] psc_bnd_fields_impl: rad takes bg into account --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 40 +++++++++++-------- 1 file changed, 24 insertions(+), 16 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index efc088fd4..d41ef3cf5 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -558,15 +558,19 @@ struct BndFields_ : BndFieldsBase neg1_d2[d2] = -1; F(H2, i3) = (/* + 4.f * C_s_pulse_y1(x,y,z+0.5*dz,t), where d0=y */ - -2.f * F(E1, edge_idx) - + -2.f * (F(E1, edge_idx) - background_e[d1]) - dtdx[d2] * (F(H0, edge_idx) - F(H0, edge_idx + neg1_d2)) - - (1.f - dtdx[d0]) * F(H2, edge_idx) + dt * F(J1, edge_idx)) / - (1.f + dtdx[d0]); + (1.f - dtdx[d0]) * (F(H2, edge_idx) - background_h[d2]) + + dt * F(J1, edge_idx)) / + (1.f + dtdx[d0]) + + background_h[d2]; F(H1, i3) = (/* + 4.f * C_p_pulse_y1(x+.5*dx,y,z,t), where d0=y */ - +2.f * F(E2, edge_idx) - + +2.f * (F(E2, edge_idx) - background_e[d2]) - dtdx[d1] * (F(H0, edge_idx) - F(H0, edge_idx + neg1_d1)) - - (1.f - dtdx[d0]) * F(H1, edge_idx) + dt * F(J2, edge_idx)) / - (1.f + dtdx[d0]); + (1.f - dtdx[d0]) * (F(H1, edge_idx) - background_h[d1]) + + dt * F(J2, edge_idx)) / + (1.f + dtdx[d0]) + + background_h[d1]; } } @@ -600,16 +604,20 @@ struct BndFields_ : BndFieldsBase Int3 neg1_d2{0, 0, 0}; neg1_d2[d2] = -1; - F(H2, i3) = - (/* + 4.f * C_s_pulse_y2(x,y,z+0.5*dz,t), where d0=y */ - +2.f * F(E1, i3) + dtdx[d2] * (F(H0, i3) - F(H0, i3 + neg1_d2)) - - (1.f - dtdx[d0]) * F(H2, edge_idx) - dt * F(J1, i3)) / - (1.f + dtdx[d0]); - F(H1, i3) = - (/* + 4.f * C_p_pulse_y2(x+.5*dx,y,z,t), where d0=y */ - -2.f * F(E2, i3) + dtdx[d1] * (F(H0, i3) - F(H0, i3 + neg1_d1)) - - (1.f - dtdx[d0]) * F(H1, edge_idx) - dt * F(J2, i3)) / - (1.f + dtdx[d0]); + F(H2, i3) = (/* + 4.f * C_s_pulse_y2(x,y,z+0.5*dz,t), where d0=y */ + +2.f * (F(E1, i3) - background_e[d1]) + + dtdx[d2] * (F(H0, i3) - F(H0, i3 + neg1_d2)) - + (1.f - dtdx[d0]) * (F(H2, edge_idx) - background_h[d2]) - + dt * F(J1, i3)) / + (1.f + dtdx[d0]) + + background_h[d2]; + F(H1, i3) = (/* + 4.f * C_p_pulse_y2(x+.5*dx,y,z,t), where d0=y */ + -2.f * (F(E2, i3) - background_e[d2]) + + dtdx[d1] * (F(H0, i3) - F(H0, i3 + neg1_d1)) - + (1.f - dtdx[d0]) * (F(H1, edge_idx) - background_h[d1]) - + dt * F(J2, i3)) / + (1.f + dtdx[d0]) + + background_h[d1]; } } From 3c241ee2aedd8af559f4140049c258440fadd5ce Mon Sep 17 00:00:00 2001 From: James McClung Date: Sat, 7 Mar 2026 12:11:12 -0500 Subject: [PATCH 31/36] psc_bnd_fields_impl: -fields_t_set_nan --- src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index d41ef3cf5..4ba023fac 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -178,11 +178,6 @@ struct BndFields_ : BndFieldsBase } } - static void fields_t_set_nan(real_t* f) - { - *f = std::numeric_limits::quiet_NaN(); - } - void set_lower_ghosts_to_nan(MfieldsState& mflds, int p, int d, int mb, bool include_edge) { From a5e3debb29a6ab83343373a1e819edc79d186683 Mon Sep 17 00:00:00 2001 From: James McClung Date: Sat, 7 Mar 2026 12:12:10 -0500 Subject: [PATCH 32/36] psc_bnd_fields_impl: make static --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 4ba023fac..f596d1fda 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -178,8 +178,8 @@ struct BndFields_ : BndFieldsBase } } - void set_lower_ghosts_to_nan(MfieldsState& mflds, int p, int d, int mb, - bool include_edge) + static void set_lower_ghosts_to_nan(MfieldsState& mflds, int p, int d, int mb, + bool include_edge) { #ifndef DEBUG return; @@ -188,8 +188,8 @@ struct BndFields_ : BndFieldsBase set_lower_ghosts(mflds, p, d, mb, {nan, nan, nan}, include_edge); } - void set_upper_ghosts_to_nan(MfieldsState& mflds, int p, int d, int mb, - bool include_edge) + static void set_upper_ghosts_to_nan(MfieldsState& mflds, int p, int d, int mb, + bool include_edge) { #ifndef DEBUG return; @@ -209,8 +209,8 @@ struct BndFields_ : BndFieldsBase * @param include_edge whether or not values located on exact domain edges * should be considered "ghosts" */ - void set_lower_ghosts(MfieldsState& mflds, int p, int d, int mb, Real3 val, - bool include_edge) + static void set_lower_ghosts(MfieldsState& mflds, int p, int d, int mb, + Real3 val, bool include_edge) { auto F = make_Fields3d(mflds[p]); Int3 start = mflds.ib(); @@ -257,8 +257,8 @@ struct BndFields_ : BndFieldsBase * @param include_edge whether or not values located on exact domain edges * should be considered "ghosts" */ - void set_upper_ghosts(MfieldsState& mflds, int p, int d, int mb, Real3 val, - bool include_edge) + static void set_upper_ghosts(MfieldsState& mflds, int p, int d, int mb, + Real3 val, bool include_edge) { auto F = make_Fields3d(mflds[p]); Int3 start = mflds.ib(); From 6907dc57fd89118047f41c493c0ec4fc74336b90 Mon Sep 17 00:00:00 2001 From: James McClung Date: Sat, 7 Mar 2026 12:38:47 -0500 Subject: [PATCH 33/36] psc: make bndf public necessary to set background e, h --- src/include/psc.hxx | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/include/psc.hxx b/src/include/psc.hxx index 93a8feed7..3adc7fd66 100644 --- a/src/include/psc.hxx +++ b/src/include/psc.hxx @@ -215,12 +215,12 @@ struct Psc void pre_first_step() { - bndf_.fill_ghosts_H(mflds_); + bndf.fill_ghosts_H(mflds_); bnd_.fill_ghosts(mflds_, HX, HX + 3); bnd_.fill_ghosts(mflds_, JXI, JXI + 3); - bndf_.fill_ghosts_E(mflds_); + bndf.fill_ghosts_E(mflds_); bnd_.fill_ghosts(mflds_, EX, EX + 3); // FIXME: do a half-step on p to bring it to its natural time, @@ -399,7 +399,7 @@ struct Psc mpi_printf(comm, "***** Bnd fields J...\n"); prof_start(pr_bndf); - bndf_.add_ghosts_J(mflds_); + bndf.add_ghosts_J(mflds_); bnd_.add_ghosts(mflds_, JXI, JXI + 3); bnd_.fill_ghosts(mflds_, JXI, JXI + 3); prof_stop(pr_bndf); @@ -413,7 +413,7 @@ struct Psc mpi_printf(comm, "***** Bnd fields B (1 of 2)...\n"); prof_restart(pr_bndf); - bndf_.fill_ghosts_H(mflds_); + bndf.fill_ghosts_H(mflds_); bnd_.fill_ghosts(mflds_, HX, HX + 3); prof_stop(pr_bndf); // state is now: x^{n+3/2}, p^{n+1}, E^{n+1/2}, B^{n+1}, j^{n+1} @@ -426,7 +426,7 @@ struct Psc mpi_printf(comm, "***** Bnd fields E...\n"); prof_restart(pr_bndf); - bndf_.fill_ghosts_E(mflds_); + bndf.fill_ghosts_E(mflds_); bnd_.fill_ghosts(mflds_, EX, EX + 3); prof_stop(pr_bndf); // state is now: x^{n+3/2}, p^{n+1}, E^{n+3/2}, B^{n+1} @@ -439,7 +439,7 @@ struct Psc mpi_printf(comm, "***** Bnd fields B (2 of 2)...\n"); prof_restart(pr_bndf); - bndf_.fill_ghosts_H(mflds_); + bndf.fill_ghosts_H(mflds_); bnd_.fill_ghosts(mflds_, HX, HX + 3); prof_stop(pr_bndf); // state is now: x^{n+3/2}, p^{n+1}, E^{n+3/2}, B^{n+3/2} @@ -517,6 +517,7 @@ private: public: const Grid_t& grid() { return *grid_; } + BndFields bndf; private: double time_start_; @@ -540,7 +541,6 @@ protected: PushParticles pushp_; PushFields pushf_; Bnd bnd_; - BndFields bndf_; BndParticles bndp_; Checkpointing checkpointing_; From 2c652ec48fdfcacd41ab339523af3f0cd3dabbf8 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 13 Mar 2026 12:05:15 -0400 Subject: [PATCH 34/36] checks_impl: disable lower check for open BCs --- src/libpsc/psc_checks/checks_impl.hxx | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/libpsc/psc_checks/checks_impl.hxx b/src/libpsc/psc_checks/checks_impl.hxx index 1b098fec0..d9d2fea75 100644 --- a/src/libpsc/psc_checks/checks_impl.hxx +++ b/src/libpsc/psc_checks/checks_impl.hxx @@ -146,12 +146,15 @@ public: for (int p = 0; p < grid.n_patches(); p++) { for (int d = 0; d < 3; d++) { - Int3 r = grid.ldims; - r[d] = 1; + if (grid.atBoundaryLo(p, d) && + (grid.bc.fld_lo[d] == BND_FLD_CONDUCTING_WALL || + grid.bc.fld_lo[d] == BND_FLD_OPEN)) { + + // account for implicit surface charges + + Int3 r = grid.ldims; + r[d] = 1; - // account for implicit surface charges - if (grid.bc.fld_lo[d] == BND_FLD_CONDUCTING_WALL && - grid.atBoundaryLo(p, d)) { rho.view(_s(0, r[0]), _s(0, r[1]), _s(0, r[2]), 0, p) = dive.view(_s(0, r[0]), _s(0, r[1]), _s(0, r[2]), 0, p); } From e8892ff3953c90831358b2aa63865d9602ca0a42 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 13 Mar 2026 12:08:02 -0400 Subject: [PATCH 35/36] psc_bnd_fields_impl: fix stop --- src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index f596d1fda..41cab1876 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -214,7 +214,7 @@ struct BndFields_ : BndFieldsBase { auto F = make_Fields3d(mflds[p]); Int3 start = mflds.ib(); - Int3 stop = mflds.im(); + Int3 stop = mflds.ib() + mflds.im(); stop[d] = 0; // TODO use gtensor views instead of VecRange @@ -262,7 +262,7 @@ struct BndFields_ : BndFieldsBase { auto F = make_Fields3d(mflds[p]); Int3 start = mflds.ib(); - Int3 stop = mflds.im(); + Int3 stop = mflds.ib() + mflds.im(); start[d] = mflds.grid().ldims[d] + 1; // TODO use gtensor views instead of VecRange @@ -538,7 +538,7 @@ struct BndFields_ : BndFieldsBase int J1 = JXI + d1, J2 = JXI + d2; Int3 start = mflds.ib(); - Int3 stop = mflds.im(); + Int3 stop = mflds.ib() + mflds.im(); start[d0] = -1; stop[d0] = 0; @@ -585,7 +585,7 @@ struct BndFields_ : BndFieldsBase int J1 = JXI + d1, J2 = JXI + d2; Int3 start = mflds.ib(); - Int3 stop = mflds.im(); + Int3 stop = mflds.ib() + mflds.im(); start[d0] = grid.ldims[d0]; stop[d0] = grid.ldims[d0] + 1; From f37c514690b13bddca9412f0b9157cfcf1708079 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 13 Mar 2026 12:12:04 -0400 Subject: [PATCH 36/36] psc_bnd_fields_impl: missing some stops --- src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 41cab1876..afe33fae5 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -230,7 +230,7 @@ struct BndFields_ : BndFieldsBase } Int3 edge_start = mflds.ib(); - Int3 edge_stop = mflds.im(); + Int3 edge_stop = mflds.ib() + mflds.im(); edge_start[d] = 0; edge_stop[d] = 1; @@ -274,7 +274,7 @@ struct BndFields_ : BndFieldsBase } Int3 edge_start = mflds.ib(); - Int3 edge_stop = mflds.im(); + Int3 edge_stop = mflds.ib() + mflds.im(); edge_start[d] = mflds.grid().ldims[d]; edge_stop[d] = mflds.grid().ldims[d] + 1;