From 5a13b21e1c4ca07e65b983b0a74e42cd3640a58a Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 23 Mar 2026 08:46:08 -0400 Subject: [PATCH 01/16] +radiating_bnd --- src/libpsc/psc_bnd_fields/radiating_bnd.hxx | 11 +++++++++++ 1 file changed, 11 insertions(+) create mode 100644 src/libpsc/psc_bnd_fields/radiating_bnd.hxx diff --git a/src/libpsc/psc_bnd_fields/radiating_bnd.hxx b/src/libpsc/psc_bnd_fields/radiating_bnd.hxx new file mode 100644 index 000000000..42d6e1229 --- /dev/null +++ b/src/libpsc/psc_bnd_fields/radiating_bnd.hxx @@ -0,0 +1,11 @@ +#include "kg/Vec3.h" + +template +class RadiatingBoundary +{ + virtual real_t pulse_s_lower(int d, int p, Int3 i3) = 0; + virtual real_t pulse_p_lower(int d, int p, Int3 i3) = 0; + + virtual real_t pulse_s_upper(int d, int p, Int3 i3) = 0; + virtual real_t pulse_p_upper(int d, int p, Int3 i3) = 0; +}; \ No newline at end of file From c15c728f6980da7af945e86a69b0a791a186c359 Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 23 Mar 2026 09:27:55 -0400 Subject: [PATCH 02/16] psc_bnd_fields_impl: +radiation --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 29 ++++++++++++++----- 1 file changed, 21 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 afe33fae5..b9a98150b 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -3,6 +3,7 @@ #include "kg/VecRange.hxx" #include "fields.hxx" #include "bnd_fields.hxx" +#include "radiating_bnd.hxx" #include @@ -552,15 +553,20 @@ 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) - background_e[d1]) - + real_t s = 0.0; + real_t p = 0.0; + if (radiation) { + s = radiation->pulse_s_lower(d0, p, i3); + p = radiation->pulse_p_lower(d0, p, i3); + } + + F(H2, i3) = (4.f * s - 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) - 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) - background_e[d2]) - + F(H1, i3) = (-4.f * p + 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) - background_h[d1]) + dt * F(J2, edge_idx)) / @@ -599,15 +605,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) - background_e[d1]) + + real_t s = 0.0; + real_t p = 0.0; + if (radiation) { + s = radiation->pulse_s_upper(d0, p, i3); + p = radiation->pulse_p_upper(d0, p, i3); + } + + F(H2, i3) = (-4.f * s + 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]) + + F(H1, i3) = (4.f * p - 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)) / @@ -618,6 +629,8 @@ struct BndFields_ : BndFieldsBase Vec3 background_e = {0.0, 0.0, 0.0}; Vec3 background_h = {0.0, 0.0, 0.0}; + + RadiatingBoundary* radiation; }; // ====================================================================== From 769442f5941bfeadb098050baf6fb57b1285dfd7 Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 23 Mar 2026 10:52:30 -0400 Subject: [PATCH 03/16] radiating_bnd; *: +time arg --- src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx | 8 ++++---- src/libpsc/psc_bnd_fields/radiating_bnd.hxx | 8 ++++---- 2 files 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 b9a98150b..666ffe2cc 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -556,8 +556,8 @@ struct BndFields_ : BndFieldsBase real_t s = 0.0; real_t p = 0.0; if (radiation) { - s = radiation->pulse_s_lower(d0, p, i3); - p = radiation->pulse_p_lower(d0, p, i3); + s = radiation->pulse_s_lower(grid.time(), d0, p, i3); + p = radiation->pulse_p_lower(grid.time(), d0, p, i3); } F(H2, i3) = (4.f * s - 2.f * (F(E1, edge_idx) - background_e[d1]) - @@ -608,8 +608,8 @@ struct BndFields_ : BndFieldsBase real_t s = 0.0; real_t p = 0.0; if (radiation) { - s = radiation->pulse_s_upper(d0, p, i3); - p = radiation->pulse_p_upper(d0, p, i3); + s = radiation->pulse_s_upper(grid.time(), d0, p, i3); + p = radiation->pulse_p_upper(grid.time(), d0, p, i3); } F(H2, i3) = (-4.f * s + 2.f * (F(E1, i3) - background_e[d1]) + diff --git a/src/libpsc/psc_bnd_fields/radiating_bnd.hxx b/src/libpsc/psc_bnd_fields/radiating_bnd.hxx index 42d6e1229..40a4c2f6b 100644 --- a/src/libpsc/psc_bnd_fields/radiating_bnd.hxx +++ b/src/libpsc/psc_bnd_fields/radiating_bnd.hxx @@ -3,9 +3,9 @@ template class RadiatingBoundary { - virtual real_t pulse_s_lower(int d, int p, Int3 i3) = 0; - virtual real_t pulse_p_lower(int d, int p, Int3 i3) = 0; + virtual real_t pulse_s_lower(double t, int d, int p, Int3 i3) = 0; + virtual real_t pulse_p_lower(double t, int d, int p, Int3 i3) = 0; - virtual real_t pulse_s_upper(int d, int p, Int3 i3) = 0; - virtual real_t pulse_p_upper(int d, int p, Int3 i3) = 0; + virtual real_t pulse_s_upper(double t, int d, int p, Int3 i3) = 0; + virtual real_t pulse_p_upper(double t, int d, int p, Int3 i3) = 0; }; \ No newline at end of file From b0266d74490f35d46663f0a7ad72099bc538f9fe Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 23 Mar 2026 11:09:10 -0400 Subject: [PATCH 04/16] Vec3: +unit ctor --- src/kg/include/kg/Vec3.h | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/kg/include/kg/Vec3.h b/src/kg/include/kg/Vec3.h index b4c4fc4ea..02c1bbcd8 100644 --- a/src/kg/include/kg/Vec3.h +++ b/src/kg/include/kg/Vec3.h @@ -38,6 +38,15 @@ struct Vec : gt::sarray return res; } + // ---------------------------------------------------------------------- + // unit vector ctors + KG_INLINE static Vec unit(int d) + { + Vec res; + res[d] = value_type(1); + return res; + } + // ---------------------------------------------------------------------- // converting to Vec of different type (e.g., float -> double) From 1baeeea1e060824457784241698a1572e464a2ca Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 23 Mar 2026 11:12:03 -0400 Subject: [PATCH 05/16] psc_bnd_fields_impl: use unit --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 48 +++++++------------ 1 file changed, 18 insertions(+), 30 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 666ffe2cc..6d18a7a7b 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -544,14 +544,7 @@ struct BndFields_ : BndFieldsBase 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; + Int3 edge_idx = i3 + Int3::unit(d0); real_t s = 0.0; real_t p = 0.0; @@ -560,18 +553,20 @@ struct BndFields_ : BndFieldsBase p = radiation->pulse_p_lower(grid.time(), d0, p, i3); } - F(H2, i3) = (4.f * s - 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) - background_h[d2]) + - dt * F(J1, edge_idx)) / - (1.f + dtdx[d0]) + - background_h[d2]; - F(H1, i3) = (-4.f * p + 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) - background_h[d1]) + - dt * F(J2, edge_idx)) / - (1.f + dtdx[d0]) + - background_h[d1]; + F(H2, i3) = + (4.f * s - 2.f * (F(E1, edge_idx) - background_e[d1]) - + dtdx[d2] * (F(H0, edge_idx) - F(H0, edge_idx - Int3::unit(d2))) - + (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 * p + 2.f * (F(E2, edge_idx) - background_e[d2]) - + dtdx[d1] * (F(H0, edge_idx) - F(H0, edge_idx - Int3::unit(d1))) - + (1.f - dtdx[d0]) * (F(H1, edge_idx) - background_h[d1]) + + dt * F(J2, edge_idx)) / + (1.f + dtdx[d0]) + + background_h[d1]; } } @@ -596,14 +591,7 @@ struct BndFields_ : BndFieldsBase stop[d0] = grid.ldims[d0] + 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; + Int3 edge_idx = i3 - Int3::unit(d0); real_t s = 0.0; real_t p = 0.0; @@ -613,13 +601,13 @@ struct BndFields_ : BndFieldsBase } F(H2, i3) = (-4.f * s + 2.f * (F(E1, i3) - background_e[d1]) + - dtdx[d2] * (F(H0, i3) - F(H0, i3 + neg1_d2)) - + dtdx[d2] * (F(H0, i3) - F(H0, i3 - Int3::unit(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 * p - 2.f * (F(E2, i3) - background_e[d2]) + - dtdx[d1] * (F(H0, i3) - F(H0, i3 + neg1_d1)) - + dtdx[d1] * (F(H0, i3) - F(H0, i3 - Int3::unit(d1))) - (1.f - dtdx[d0]) * (F(H1, edge_idx) - background_h[d1]) - dt * F(J2, i3)) / (1.f + dtdx[d0]) + From 7dd710dd6e0c3090932f0d1948472a2f6df3930c Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 23 Mar 2026 11:18:43 -0400 Subject: [PATCH 06/16] radiating_bnd; *: pass x3 instead of i3 --- src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx | 12 ++++++++---- src/libpsc/psc_bnd_fields/radiating_bnd.hxx | 10 ++++++---- 2 files changed, 14 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 6d18a7a7b..f3e256398 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -549,8 +549,10 @@ struct BndFields_ : BndFieldsBase real_t s = 0.0; real_t p = 0.0; if (radiation) { - s = radiation->pulse_s_lower(grid.time(), d0, p, i3); - p = radiation->pulse_p_lower(grid.time(), d0, p, i3); + Real3 x3_s = (Real3(i3) + Real3::unit(d1) * 0.5) * grid.domain.dx; + Real3 x3_p = (Real3(i3) + Real3::unit(d2) * 0.5) * grid.domain.dx; + s = radiation->pulse_s_lower(grid.time(), d0, p, x3_s); + p = radiation->pulse_p_lower(grid.time(), d0, p, x3_p); } F(H2, i3) = @@ -596,8 +598,10 @@ struct BndFields_ : BndFieldsBase real_t s = 0.0; real_t p = 0.0; if (radiation) { - s = radiation->pulse_s_upper(grid.time(), d0, p, i3); - p = radiation->pulse_p_upper(grid.time(), d0, p, i3); + Real3 x3_s = (Real3(i3) + Real3::unit(d1) * 0.5) * grid.domain.dx; + Real3 x3_p = (Real3(i3) + Real3::unit(d2) * 0.5) * grid.domain.dx; + s = radiation->pulse_s_upper(grid.time(), d0, p, x3_s); + p = radiation->pulse_p_upper(grid.time(), d0, p, x3_p); } F(H2, i3) = (-4.f * s + 2.f * (F(E1, i3) - background_e[d1]) + diff --git a/src/libpsc/psc_bnd_fields/radiating_bnd.hxx b/src/libpsc/psc_bnd_fields/radiating_bnd.hxx index 40a4c2f6b..adf528f8c 100644 --- a/src/libpsc/psc_bnd_fields/radiating_bnd.hxx +++ b/src/libpsc/psc_bnd_fields/radiating_bnd.hxx @@ -3,9 +3,11 @@ template class RadiatingBoundary { - virtual real_t pulse_s_lower(double t, int d, int p, Int3 i3) = 0; - virtual real_t pulse_p_lower(double t, int d, int p, Int3 i3) = 0; + using Real3 = Vec3; - virtual real_t pulse_s_upper(double t, int d, int p, Int3 i3) = 0; - virtual real_t pulse_p_upper(double t, int d, int p, Int3 i3) = 0; + virtual real_t pulse_s_lower(double t, int d, int p, Real3 x3) = 0; + virtual real_t pulse_p_lower(double t, int d, int p, Real3 x3) = 0; + + virtual real_t pulse_s_upper(double t, int d, int p, Real3 x3) = 0; + virtual real_t pulse_p_upper(double t, int d, int p, Real3 x3) = 0; }; \ No newline at end of file From 23a30f1e1894dbf9921137a43d69c7b2e338749b Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 23 Mar 2026 11:19:57 -0400 Subject: [PATCH 07/16] radiating_bnd: is struct --- src/libpsc/psc_bnd_fields/radiating_bnd.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libpsc/psc_bnd_fields/radiating_bnd.hxx b/src/libpsc/psc_bnd_fields/radiating_bnd.hxx index adf528f8c..144e58409 100644 --- a/src/libpsc/psc_bnd_fields/radiating_bnd.hxx +++ b/src/libpsc/psc_bnd_fields/radiating_bnd.hxx @@ -1,7 +1,7 @@ #include "kg/Vec3.h" template -class RadiatingBoundary +struct RadiatingBoundary { using Real3 = Vec3; From f4722731d08dfdc216b0c71aeb8a5bbd16375ad4 Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 23 Mar 2026 11:57:39 -0400 Subject: [PATCH 08/16] radiating_bnd; psc_bnd_fields_impl: +update_cache --- src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx | 8 ++++++++ src/libpsc/psc_bnd_fields/radiating_bnd.hxx | 4 ++++ 2 files changed, 12 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 f3e256398..8dc621508 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -85,6 +85,10 @@ struct BndFields_ : BndFieldsBase for (int p = 0; p < mflds.n_patches(); p++) { // lo for (int d = 0; d < 3; d++) { + if (grid.bc.fld_lo[d] == BND_FLD_OPEN && radiation) { + radiation->update_cache_lower(grid.time(), d); + } + if (grid.atBoundaryLo(p, d)) { switch (grid.bc.fld_lo[d]) { case BND_FLD_PERIODIC: { @@ -106,6 +110,10 @@ struct BndFields_ : BndFieldsBase } // hi for (int d = 0; d < 3; d++) { + if (grid.bc.fld_hi[d] == BND_FLD_OPEN && radiation) { + radiation->update_cache_upper(grid.time(), d); + } + if (grid.atBoundaryHi(p, d)) { switch (grid.bc.fld_hi[d]) { case BND_FLD_PERIODIC: { diff --git a/src/libpsc/psc_bnd_fields/radiating_bnd.hxx b/src/libpsc/psc_bnd_fields/radiating_bnd.hxx index 144e58409..c6a40de9f 100644 --- a/src/libpsc/psc_bnd_fields/radiating_bnd.hxx +++ b/src/libpsc/psc_bnd_fields/radiating_bnd.hxx @@ -10,4 +10,8 @@ struct RadiatingBoundary virtual real_t pulse_s_upper(double t, int d, int p, Real3 x3) = 0; virtual real_t pulse_p_upper(double t, int d, int p, Real3 x3) = 0; + + // FIXME these are a hack for a specific subclass to work + virtual void update_cache_lower(double t, int d) {} + virtual void update_cache_upper(double t, int d) {} }; \ No newline at end of file From e524f023d35698c4636bf17cfe30d1356484e26e Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 23 Mar 2026 16:23:30 -0400 Subject: [PATCH 09/16] Vec3: +reverse --- src/kg/include/kg/Vec3.h | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/kg/include/kg/Vec3.h b/src/kg/include/kg/Vec3.h index 02c1bbcd8..b74c7db0c 100644 --- a/src/kg/include/kg/Vec3.h +++ b/src/kg/include/kg/Vec3.h @@ -286,6 +286,15 @@ struct Vec : gt::sarray return res; } + KG_INLINE Vec reverse() const + { + Vec res; + for (int i = 0; i < N; i++) { + res[i] = (*this)[N - 1 - i]; + } + return res; + } + // conversion to pointer KG_INLINE operator const T*() const { return this->data(); } From 9f18f7bec6117b2e1b01ddf73806107de8627980 Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 23 Mar 2026 16:46:41 -0400 Subject: [PATCH 10/16] psc_bnd_fields_impl: fix pulse location I think... this warrants a double check --- 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 8dc621508..8544efa5e 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -557,8 +557,8 @@ struct BndFields_ : BndFieldsBase real_t s = 0.0; real_t p = 0.0; if (radiation) { - Real3 x3_s = (Real3(i3) + Real3::unit(d1) * 0.5) * grid.domain.dx; - Real3 x3_p = (Real3(i3) + Real3::unit(d2) * 0.5) * grid.domain.dx; + Real3 x3_s = (Real3(edge_idx) + Real3::unit(d1) * 0.5) * grid.domain.dx; + Real3 x3_p = (Real3(edge_idx) + Real3::unit(d2) * 0.5) * grid.domain.dx; s = radiation->pulse_s_lower(grid.time(), d0, p, x3_s); p = radiation->pulse_p_lower(grid.time(), d0, p, x3_p); } @@ -606,8 +606,8 @@ struct BndFields_ : BndFieldsBase real_t s = 0.0; real_t p = 0.0; if (radiation) { - Real3 x3_s = (Real3(i3) + Real3::unit(d1) * 0.5) * grid.domain.dx; - Real3 x3_p = (Real3(i3) + Real3::unit(d2) * 0.5) * grid.domain.dx; + Real3 x3_s = (Real3(edge_idx) + Real3::unit(d1) * 0.5) * grid.domain.dx; + Real3 x3_p = (Real3(edge_idx) + Real3::unit(d2) * 0.5) * grid.domain.dx; s = radiation->pulse_s_upper(grid.time(), d0, p, x3_s); p = radiation->pulse_p_upper(grid.time(), d0, p, x3_p); } From e03faa2b892feee89b97cdee1022847fa69a9905 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 17 Feb 2026 14:02:36 -0500 Subject: [PATCH 11/16] psc_shock: +AdvectedPeriodicFields --- src/psc_shock.cxx | 151 +++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 148 insertions(+), 3 deletions(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index dd3159a64..7022e7f8d 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -32,6 +32,8 @@ using Collision = PscConfig::Collision; using Checks = PscConfig::Checks; using Marder = PscConfig::Marder; using OutputParticles = PscConfig::OutputParticles; +using real_t = PscConfig::Mfields::real_t; +using Real3 = Vec3; // ====================================================================== // Global parameters @@ -662,6 +664,145 @@ void initializeFields(MfieldsState& mflds) boost_fields(mflds, -v_upstream_y); } +struct AdvectedPeriodicFields : RadiatingBoundary +{ + static const int DIM_Y = 1; + + AdvectedPeriodicFields(MfieldsState& mflds, real_t v_advect, + Real3 background_e, Real3 background_h) + : v_advect(v_advect), grid(mflds.grid()) + { + // FIXME would be better to exclude J, but the interpolator uses EX, etc. + auto&& e_b_fields = mflds.storage().view(_all, _all, _all, _all, _all); + // FIXME this probably isn't the best way to copy a gtensor array + cycled_fields = gt::zeros_like(e_b_fields); + cycled_fields.view(_all, _all, _all, _all, _all) = e_b_fields; + + for (int d = 0; d < 3; d++) { + cycled_fields.view(_all, _all, _all, EX + d, _all) = + cycled_fields.view(_all, _all, _all, EX + d, _all) - background_e[d]; + cycled_fields.view(_all, _all, _all, HX + d, _all) = + cycled_fields.view(_all, _all, _all, HX + d, _all) - background_h[d]; + } + } + + Real3 advect_x3(Real3 x3, double t) + { + return x3 - Real3::unit(DIM_Y) * real_t(t * v_advect); + } + + int shift_to_patch_local(Real3& x3_advected) + { + real_t patch_size = grid.domain.length[DIM_Y] / grid.domain.np[DIM_Y]; + int n_patches_to_the_left = 0; + while (x3_advected[DIM_Y] < grid.domain.corner[DIM_Y]) { + x3_advected[DIM_Y] += patch_size; + n_patches_to_the_left += 1; + } + return n_patches_to_the_left; + } + + void calc_e_h(double t, int p, Real3 x3, int d_e, real_t& e, int d_h, + real_t& h) + { + Real3 x3_advected = advect_x3(x3, t); + int n_patches_to_the_left = shift_to_patch_local(x3_advected); + assert(n_patches_to_the_left == n_patch_cycles); + + ip.set_coeffs(x3_advected * grid.domain.dx_inv); + auto em = decltype(ip)::fields_t( + cycled_fields.view(_all, _all, _all, _all, p), -grid.ibn); + + switch (d_e) { + case 0: e = ip.ex(em); break; + case 1: e = ip.ey(em); break; + case 2: e = ip.ez(em); break; + } + switch (d_h) { + case 0: h = ip.hx(em); break; + case 1: h = ip.hy(em); break; + case 2: h = ip.hz(em); break; + } + } + + real_t pulse_s_lower(double t, int d, int p, Real3 x3) override + { + int d1 = (d + 1) % 3; + int d2 = (d + 2) % 3; + + real_t e, h; + calc_e_h(t, p, x3, d1, e, d2, h); + + return (e + h) / 2.0; + } + + real_t pulse_p_lower(double t, int d, int p, Real3 x3) override + { + int d1 = (d + 1) % 3; + int d2 = (d + 2) % 3; + + real_t e, h; + calc_e_h(t, p, x3, d2, e, d1, h); + + return (e - h) / 2.0; + } + + real_t pulse_s_upper(double t, int d, int p, Real3 x3) override + { + return 0.0; + } + + real_t pulse_p_upper(double t, int d, int p, Real3 x3) override + { + return 0.0; + } + + void update_cache_lower(double t, int d) override + { + if (d != DIM_Y) { + return; + } + + // TODO make this work with >1 patch per process? + assert(grid.n_patches() == 1); + + Real3 x3_advected = advect_x3({0.0, 0.0, 0.0}, t); + int n_patches_to_the_left = shift_to_patch_local(x3_advected); + + if (n_patches_to_the_left == n_patch_cycles) { + // "cache hit"; no need to cycle + return; + } + + LOG_INFO("cycling turbulence...\n"); + + // hack: guess the rank based on how mrc does it for simple domains + // (can't use mrc, because it wouldn't apply periodicity) + Int3 np = grid.domain.np; + Int3 proc = grid.localPatchInfo(0).idx3; + Int3 dest_proc = (proc + Int3::unit(DIM_Y)) % np; + int dest_rank = flatten_index(dest_proc.reverse(), np.reverse()); + Int3 source_proc = (proc - Int3::unit(DIM_Y) + np) % np; + int source_rank = flatten_index(source_proc.reverse(), np.reverse()); + + MPI_Status status; + int err = MPI_Sendrecv_replace(cycled_fields.data(), cycled_fields.size(), + MpiDtypeTraits::value(), dest_rank, + 0, source_rank, 0, grid.comm(), &status); + + n_patch_cycles += 1; + } + + const Grid_t& grid; + gt::gtensor cycled_fields; + int n_patch_cycles = 0; + real_t v_advect; + InterpolateEM< + Fields3d, + opt_ip_1st_ec, Dim> + ip; +}; + // ====================================================================== // run @@ -761,9 +902,13 @@ static void run(int argc, char** argv) psc.add_gauss_corrector(&marder); double gamma = 1 / sqrt(1 - sqr(v_upstream_y)); - psc.bndf.background_e = - -gamma * Double3{0, v_upstream_y, 0}.cross({b_x, b_y, b_z}); - psc.bndf.background_h = {b_x * gamma, b_y, b_z * gamma}; + Real3 background_e = + -gamma * Real3{0, v_upstream_y, 0}.cross({b_x, b_y, b_z}); + Real3 background_h = {b_x * gamma, b_y, b_z * gamma}; + psc.bndf.background_e = background_e; + psc.bndf.background_h = background_h; + psc.bndf.radiation = + new AdvectedPeriodicFields{mflds, v_upstream_y, background_e, background_h}; psc.add_diagnostic(&outf); psc.add_diagnostic(&outp); From b1ce9f5a9d5321898fc37cff10165ea1ac80a7fb Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 25 Mar 2026 11:14:44 -0400 Subject: [PATCH 12/16] psc_shock: advect at speed 1 --- src/psc_shock.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index 7022e7f8d..c6f576620 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -908,7 +908,7 @@ static void run(int argc, char** argv) psc.bndf.background_e = background_e; psc.bndf.background_h = background_h; psc.bndf.radiation = - new AdvectedPeriodicFields{mflds, v_upstream_y, background_e, background_h}; + new AdvectedPeriodicFields{mflds, 1.0, background_e, background_h}; psc.add_diagnostic(&outf); psc.add_diagnostic(&outp); From 0b8e801b18dadd4052390376c00ade2eb8be42b5 Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 25 Mar 2026 11:53:54 -0400 Subject: [PATCH 13/16] psc_shock: mv bg calcs --- src/psc_shock.cxx | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index c6f576620..d3cb2de80 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -56,6 +56,9 @@ double b_x; double b_y; double b_z; +Real3 background_e; +Real3 background_h; + int nx; int ny; int nz; @@ -114,6 +117,12 @@ void setupParameters(int argc, char** argv) b_y = b_mag * cos(b_angle_y_to_x_rad); b_z = 0.0; + Real3 v_upstream = {v_upstream_x, v_upstream_y, v_upstream_z}; + double gamma = 1 / sqrt(1 - v_upstream.mag2()); + background_e = -gamma * v_upstream.cross({b_x, b_y, b_z}); + // note: this only holds for vx=vz=0 + background_h = {b_x * gamma, b_y, b_z * gamma}; + nx = inputParams.get("nx"); ny = inputParams.get("ny"); nz = inputParams.get("nz"); @@ -901,10 +910,6 @@ static void run(int argc, char** argv) psc.add_gauss_corrector(&marder); - double gamma = 1 / sqrt(1 - sqr(v_upstream_y)); - Real3 background_e = - -gamma * Real3{0, v_upstream_y, 0}.cross({b_x, b_y, b_z}); - Real3 background_h = {b_x * gamma, b_y, b_z * gamma}; psc.bndf.background_e = background_e; psc.bndf.background_h = background_h; psc.bndf.radiation = From e34760a2245c3c4075df6c5b6a508e3ac2b122df Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 26 Mar 2026 08:30:23 -0400 Subject: [PATCH 14/16] psc_bnd_fields_impl: fix scalar t --- src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx | 12 ++++++++---- 1 file changed, 8 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 8544efa5e..27bfa7533 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -557,8 +557,10 @@ struct BndFields_ : BndFieldsBase real_t s = 0.0; real_t p = 0.0; if (radiation) { - Real3 x3_s = (Real3(edge_idx) + Real3::unit(d1) * 0.5) * grid.domain.dx; - Real3 x3_p = (Real3(edge_idx) + Real3::unit(d2) * 0.5) * grid.domain.dx; + Real3 x3_s = + (Real3(edge_idx) + Real3::unit(d1) * real_t(0.5)) * grid.domain.dx; + Real3 x3_p = + (Real3(edge_idx) + Real3::unit(d2) * real_t(0.5)) * grid.domain.dx; s = radiation->pulse_s_lower(grid.time(), d0, p, x3_s); p = radiation->pulse_p_lower(grid.time(), d0, p, x3_p); } @@ -606,8 +608,10 @@ struct BndFields_ : BndFieldsBase real_t s = 0.0; real_t p = 0.0; if (radiation) { - Real3 x3_s = (Real3(edge_idx) + Real3::unit(d1) * 0.5) * grid.domain.dx; - Real3 x3_p = (Real3(edge_idx) + Real3::unit(d2) * 0.5) * grid.domain.dx; + Real3 x3_s = + (Real3(edge_idx) + Real3::unit(d1) * real_t(0.5)) * grid.domain.dx; + Real3 x3_p = + (Real3(edge_idx) + Real3::unit(d2) * real_t(0.5)) * grid.domain.dx; s = radiation->pulse_s_upper(grid.time(), d0, p, x3_s); p = radiation->pulse_p_upper(grid.time(), d0, p, x3_p); } From b098e7eccd46e022dc0e1d1e47abfed219f996ba Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 26 Mar 2026 08:33:51 -0400 Subject: [PATCH 15/16] psc_bnd_fields_impl: fix vector t --- .../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 27bfa7533..bf7dbb7c6 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -557,10 +557,10 @@ struct BndFields_ : BndFieldsBase real_t s = 0.0; real_t p = 0.0; if (radiation) { - Real3 x3_s = - (Real3(edge_idx) + Real3::unit(d1) * real_t(0.5)) * grid.domain.dx; - Real3 x3_p = - (Real3(edge_idx) + Real3::unit(d2) * real_t(0.5)) * grid.domain.dx; + Real3 x3_s = (Real3(edge_idx) + Real3::unit(d1) * real_t(0.5)) * + Real3(grid.domain.dx); + Real3 x3_p = (Real3(edge_idx) + Real3::unit(d2) * real_t(0.5)) * + Real3(grid.domain.dx); s = radiation->pulse_s_lower(grid.time(), d0, p, x3_s); p = radiation->pulse_p_lower(grid.time(), d0, p, x3_p); } @@ -608,10 +608,10 @@ struct BndFields_ : BndFieldsBase real_t s = 0.0; real_t p = 0.0; if (radiation) { - Real3 x3_s = - (Real3(edge_idx) + Real3::unit(d1) * real_t(0.5)) * grid.domain.dx; - Real3 x3_p = - (Real3(edge_idx) + Real3::unit(d2) * real_t(0.5)) * grid.domain.dx; + Real3 x3_s = (Real3(edge_idx) + Real3::unit(d1) * real_t(0.5)) * + Real3(grid.domain.dx); + Real3 x3_p = (Real3(edge_idx) + Real3::unit(d2) * real_t(0.5)) * + Real3(grid.domain.dx); s = radiation->pulse_s_upper(grid.time(), d0, p, x3_s); p = radiation->pulse_p_upper(grid.time(), d0, p, x3_p); } From 576198488c9cf2df2ad97c7ab3098abd97bb2a58 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 26 Mar 2026 08:46:49 -0400 Subject: [PATCH 16/16] psc_bnd_fields_impl: radiation=nullptr --- src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx | 2 +- 1 file changed, 1 insertion(+), 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 bf7dbb7c6..7697893c7 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -634,7 +634,7 @@ struct BndFields_ : BndFieldsBase Vec3 background_e = {0.0, 0.0, 0.0}; Vec3 background_h = {0.0, 0.0, 0.0}; - RadiatingBoundary* radiation; + RadiatingBoundary* radiation = nullptr; }; // ======================================================================