From 3c398a56009d960952dd290201ca0e6b8fa8d5eb Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 16 Dec 2025 10:44:58 -0500 Subject: [PATCH 1/8] psc_bnd_fields_impl: use std::min,max --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 98 +++++++++---------- 1 file changed, 49 insertions(+), 49 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 5bec697cf..7a368419b 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -145,8 +145,8 @@ struct BndFields_ : BndFieldsBase if (d == 1) { #ifdef DEBUG for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = MAX(-2, F.ib_[0]); - ix < MIN(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { + 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)); @@ -158,8 +158,8 @@ struct BndFields_ : BndFieldsBase #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 = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]); - ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { F(EX, ix, 0, iz) = 0.; F(EX, ix, -1, iz) = F(EX, ix, 1, iz); @@ -183,8 +183,8 @@ struct BndFields_ : BndFieldsBase } #endif for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]); - ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { F(EX, ix, iy, 0) = 0.; F(EX, ix, iy, -1) = F(EX, ix, iy, 1); @@ -209,8 +209,8 @@ struct BndFields_ : BndFieldsBase int my _mrc_unused = ldims[1]; #ifdef DEBUG for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = MAX(-2, F.ib_[0]); - ix < MIN(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { + 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)); @@ -221,8 +221,8 @@ struct BndFields_ : BndFieldsBase } #endif for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]); - ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { F(EX, ix, my, iz) = 0.; F(EX, ix, my + 1, iz) = F(EX, ix, my - 1, iz); @@ -247,8 +247,8 @@ struct BndFields_ : BndFieldsBase } #endif for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]); - ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { F(EX, ix, iy, mz) = 0.; F(EX, ix, iy, mz + 1) = F(EX, ix, iy, mz - 1); @@ -272,8 +272,8 @@ struct BndFields_ : BndFieldsBase if (d == 1) { #ifdef DEBUG for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = MAX(-2, F.ib_[0]); - ix < MIN(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { + 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)); @@ -284,8 +284,8 @@ struct BndFields_ : BndFieldsBase } #endif for (int iz = -1; iz < ldims[2] + 2; iz++) { - for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]); - ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { F(HX, ix, -1, iz) = -F(HX, ix, 0, iz); F(HY, ix, -1, iz) = F(HY, ix, 1, iz); @@ -307,8 +307,8 @@ struct BndFields_ : BndFieldsBase } #endif for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]); - ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { F(HX, ix, iy, -1) = -F(HX, ix, iy, 0); F(HY, ix, iy, -1) = -F(HY, ix, iy, 0); @@ -332,8 +332,8 @@ struct BndFields_ : BndFieldsBase int my _mrc_unused = ldims[1]; #ifdef DEBUG for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = MAX(-2, F.ib_[0]); - ix < MIN(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { + 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)); @@ -343,8 +343,8 @@ struct BndFields_ : BndFieldsBase } #endif for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]); - ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { F(HX, ix, my, iz) = -F(HX, ix, my - 1, iz); F(HY, ix, my + 1, iz) = F(HY, ix, my - 1, iz); @@ -366,8 +366,8 @@ struct BndFields_ : BndFieldsBase } #endif for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]); - ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { F(HX, ix, iy, mz) = -F(HX, ix, iy, mz - 1); F(HY, ix, iy, mz) = -F(HY, ix, iy, mz - 1); @@ -388,8 +388,8 @@ struct BndFields_ : BndFieldsBase if (d == 1) { for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]); - ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { F(JXI, ix, 1, iz) += F(JXI, ix, -1, iz); F(JXI, ix, -1, iz) = 0.; @@ -402,8 +402,8 @@ struct BndFields_ : BndFieldsBase } } else if (d == 2) { for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]); - ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { F(JXI, ix, iy, 1) += F(JXI, ix, iy, -1); F(JXI, ix, iy, -1) = 0.; @@ -428,8 +428,8 @@ struct BndFields_ : BndFieldsBase if (d == 1) { int my _mrc_unused = ldims[1]; for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]); - ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { F(JXI, ix, my - 1, iz) += F(JXI, ix, my + 1, iz); F(JXI, ix, my + 1, iz) = 0.; @@ -443,8 +443,8 @@ struct BndFields_ : BndFieldsBase } else if (d == 2) { int mz = ldims[2]; for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]); - ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { F(JXI, ix, iy, mz - 1) += F(JXI, ix, iy, mz + 1); F(JXI, ix, iy, mz + 1) = 0.; @@ -479,8 +479,8 @@ struct BndFields_ : BndFieldsBase if (d == 1) { #ifdef DEBUG for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = MAX(-2, F.ib_[0]); - ix < MIN(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { + 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)); @@ -493,8 +493,8 @@ struct BndFields_ : BndFieldsBase real_t dt = F.grid().dt, dy = F.grid().domain.dx[1], dz = F.grid().domain.dx[2]; for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = MAX(-2, F.ib_[0]); - ix < MIN(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { + for (int ix = std::max(-2, F.ib_[0]); + ix < std::min(ldims[0] + 2, F.ib_[0] + F.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) @@ -512,8 +512,8 @@ struct BndFields_ : BndFieldsBase } else if (d == 2) { #ifdef DEBUG for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = MAX(-2, F.ib_[0]); - ix < MIN(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { + 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, iy, -1)); fields_t_set_nan(&F(HX, ix, iy, -2)); fields_t_set_nan(&F(HY, ix, iy, -1)); @@ -526,8 +526,8 @@ struct BndFields_ : BndFieldsBase real_t dt = F.grid().dt, dy = F.grid().domain.dx[1], dz = F.grid().domain.dx[2]; for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = MAX(-2, F.ib_[0]); - ix < MIN(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { + for (int ix = std::max(-2, F.ib_[0]); + ix < std::min(ldims[0] + 2, F.ib_[0] + F.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) - @@ -562,8 +562,8 @@ struct BndFields_ : BndFieldsBase int my _mrc_unused = ldims[1]; #ifdef DEBUG for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = MAX(-2, F.ib_[0]); - ix < MIN(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { + 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, iz)); @@ -575,8 +575,8 @@ struct BndFields_ : BndFieldsBase real_t dt = F.grid().dt, dy = F.grid().domain.dx[1], dz = F.grid().domain.dx[2]; for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = MAX(-2, F.ib_[0]); - ix < MIN(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { + for (int ix = std::max(-2, F.ib_[0]); + ix < std::min(ldims[0] + 2, F.ib_[0] + F.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) @@ -609,8 +609,8 @@ struct BndFields_ : BndFieldsBase real_t dt = F.grid().dt, dy = F.grid().domain.dx[1], dz = F.grid().domain.dx[2]; for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = MAX(-2, F.ib_[0]); - ix < MIN(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { + for (int ix = std::max(-2, F.ib_[0]); + ix < std::min(ldims[0] + 2, F.ib_[0] + F.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) + @@ -644,7 +644,7 @@ struct BndFieldsNone : BndFieldsBase { using Mfields = MF; - void fill_ghosts_E(Mfields& mflds){}; - void fill_ghosts_H(Mfields& mflds){}; - void add_ghosts_J(Mfields& mflds){}; + void fill_ghosts_E(Mfields& mflds) {}; + void fill_ghosts_H(Mfields& mflds) {}; + void add_ghosts_J(Mfields& mflds) {}; }; From e9d022c7403c8bd63b7f82b07c7648ad62d5e1ff Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 16 Dec 2025 11:20:14 -0500 Subject: [PATCH 2/8] fields_item: replace DomainBoundary with func --- src/include/fields_item.hxx | 41 ++++++++----------------------------- 1 file changed, 9 insertions(+), 32 deletions(-) diff --git a/src/include/fields_item.hxx b/src/include/fields_item.hxx index 2a4cbe857..5c89117fc 100644 --- a/src/include/fields_item.hxx +++ b/src/include/fields_item.hxx @@ -34,23 +34,11 @@ inline std::vector addKindSuffix( // ====================================================================== // ItemMomentBnd -template -class DomainBoundary +template +void add_ghosts_boundary(const Grid_t& grid, FE& mres_gt, const Int3& ib, int p, + int mb, int me, centering::Centering c) { -public: - template - static void add_ghosts_boundary(const Grid_t& grid, FE& mres_gt, - const Int3& ib, int p, int mb, int me); -}; - -template <> -class DomainBoundary -{ -public: - template - static void add_ghosts_boundary(const Grid_t& grid, FE& mres_gt, - const Int3& ib, int p, int mb, int me) - { + if (c == centering::CC) { // lo for (int d = 0; d < 3; d++) { // FIXME why reflect for open BCs? @@ -67,17 +55,7 @@ public: add_ghosts_reflecting_hi_cc(grid.ldims, mres_gt, ib, p, d, mb, me); } } - } -}; - -template <> -class DomainBoundary -{ -public: - template - static void add_ghosts_boundary(const Grid_t& grid, FE& mres_gt, - const Int3& ib, int p, int mb, int me) - { + } else if (c == centering::NC) { // lo for (int d = 0; d < 3; d++) { // FIXME why reflect for open BCs? @@ -95,7 +73,7 @@ public: } } } -}; +} // ====================================================================== // ItemMomentBnd @@ -111,8 +89,7 @@ public: void add_ghosts(const Grid_t& grid, storage_type& mres_gt, const Int3& ib) { for (int p = 0; p < mres_gt.shape(4); p++) { - DomainBoundary::add_ghosts_boundary(grid, mres_gt, ib, p, 0, - mres_gt.shape(3)); + add_ghosts_boundary(grid, mres_gt, ib, p, 0, mres_gt.shape(3), C); } bnd_.add_ghosts(grid, mres_gt, ib, 0, mres_gt.shape(3)); @@ -152,8 +129,8 @@ public: { Int3 ib = -mprts.grid().ibn; // FIXME, gt::gtensor and psc::gtensor are slightly different, and ideally - // zeros() shouldn't actually allocate, but probably it does, so this wastes - // memory and a copy + // zeros() shouldn't actually allocate, but probably it does, so this + // wastes memory and a copy storage_type mres = psc::mflds::zeros(mprts.grid(), n_comps(), ib); moment_type{}(mres, ib, mprts); From dbcb77c70b9ea77a11e7dfd26fc5c086a22682a8 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 16 Dec 2025 11:21:47 -0500 Subject: [PATCH 3/8] fields_item: don't reflect for open bcs --- src/include/fields_item.hxx | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/src/include/fields_item.hxx b/src/include/fields_item.hxx index 5c89117fc..75c70ec3a 100644 --- a/src/include/fields_item.hxx +++ b/src/include/fields_item.hxx @@ -41,34 +41,26 @@ void add_ghosts_boundary(const Grid_t& grid, FE& mres_gt, const Int3& ib, int p, if (c == centering::CC) { // lo for (int d = 0; d < 3; d++) { - // FIXME why reflect for open BCs? - if (grid.atBoundaryLo(p, d) && (grid.bc.prt_lo[d] == BND_PRT_REFLECTING || - grid.bc.prt_lo[d] == BND_PRT_OPEN)) { + if (grid.atBoundaryLo(p, d) && grid.bc.prt_lo[d] == BND_PRT_REFLECTING) { add_ghosts_reflecting_lo_cc(grid.ldims, mres_gt, ib, p, d, mb, me); } } // hi for (int d = 0; d < 3; d++) { - // FIXME why reflect for open BCs? - if (grid.atBoundaryHi(p, d) && (grid.bc.prt_hi[d] == BND_PRT_REFLECTING || - grid.bc.prt_hi[d] == BND_PRT_OPEN)) { + if (grid.atBoundaryHi(p, d) && grid.bc.prt_hi[d] == BND_PRT_REFLECTING) { add_ghosts_reflecting_hi_cc(grid.ldims, mres_gt, ib, p, d, mb, me); } } } else if (c == centering::NC) { // lo for (int d = 0; d < 3; d++) { - // FIXME why reflect for open BCs? - if (grid.atBoundaryLo(p, d) && (grid.bc.prt_lo[d] == BND_PRT_REFLECTING || - grid.bc.prt_lo[d] == BND_PRT_OPEN)) { + if (grid.atBoundaryLo(p, d) && grid.bc.prt_lo[d] == BND_PRT_REFLECTING) { add_ghosts_reflecting_lo_nc(grid.ldims, mres_gt, ib, p, d, mb, me); } } // hi for (int d = 0; d < 3; d++) { - // FIXME why reflect for open BCs? - if (grid.atBoundaryHi(p, d) && (grid.bc.prt_hi[d] == BND_PRT_REFLECTING || - grid.bc.prt_hi[d] == BND_PRT_OPEN)) { + if (grid.atBoundaryHi(p, d) && grid.bc.prt_hi[d] == BND_PRT_REFLECTING) { add_ghosts_reflecting_hi_nc(grid.ldims, mres_gt, ib, p, d, mb, me); } } From 5d7bfcb17bcd5800bce7c151009bcb761e752ba4 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 16 Dec 2025 11:42:34 -0500 Subject: [PATCH 4/8] fields_item: rename func to maybe_... --- src/include/fields_item.hxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/include/fields_item.hxx b/src/include/fields_item.hxx index 75c70ec3a..83546a28c 100644 --- a/src/include/fields_item.hxx +++ b/src/include/fields_item.hxx @@ -35,8 +35,8 @@ inline std::vector addKindSuffix( // ItemMomentBnd template -void add_ghosts_boundary(const Grid_t& grid, FE& mres_gt, const Int3& ib, int p, - int mb, int me, centering::Centering c) +void maybe_add_ghosts_boundary(const Grid_t& grid, FE& mres_gt, const Int3& ib, + int p, int mb, int me, centering::Centering c) { if (c == centering::CC) { // lo @@ -81,7 +81,7 @@ public: void add_ghosts(const Grid_t& grid, storage_type& mres_gt, const Int3& ib) { for (int p = 0; p < mres_gt.shape(4); p++) { - add_ghosts_boundary(grid, mres_gt, ib, p, 0, mres_gt.shape(3), C); + maybe_add_ghosts_boundary(grid, mres_gt, ib, p, 0, mres_gt.shape(3), C); } bnd_.add_ghosts(grid, mres_gt, ib, 0, mres_gt.shape(3)); From d91a96aab4639d7688f04d3e16c94d4be9e7a6c9 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 16 Dec 2025 14:13:17 -0500 Subject: [PATCH 5/8] checks_*: +exit_on_failure option --- src/include/checks_params.hxx | 1 + src/libpsc/psc_checks/checks_impl.hxx | 11 ++++++++--- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/src/include/checks_params.hxx b/src/include/checks_params.hxx index 8b328f6bc..3cc26f8d0 100644 --- a/src/include/checks_params.hxx +++ b/src/include/checks_params.hxx @@ -7,6 +7,7 @@ struct CheckParams double err_threshold = 1e-13; // maximum acceptable error bool print_max_err_always = true; // always print error, even if acceptable bool dump_always = false; // always dump compared fields, even if acceptable + bool exit_on_failure = false; // exit the program if check fails bool enabled() { return check_interval > 0; } diff --git a/src/libpsc/psc_checks/checks_impl.hxx b/src/libpsc/psc_checks/checks_impl.hxx index 2eb70759a..071a6a9f1 100644 --- a/src/libpsc/psc_checks/checks_impl.hxx +++ b/src/libpsc/psc_checks/checks_impl.hxx @@ -1,4 +1,3 @@ - #pragma once #include "fields.hxx" @@ -102,7 +101,10 @@ public: MPI_Barrier(grid.comm()); } - assert(max_err < err_threshold); + if (exit_on_failure && max_err >= err_threshold) { + LOG_ERROR("exiting due to failed continuity check\n"); + } + last_max_err = max_err; } @@ -185,7 +187,10 @@ public: writer_.end_step(); } - assert(max_err < err_threshold); + if (exit_on_failure && max_err >= err_threshold) { + LOG_ERROR("exiting due to failed gauss check\n"); + } + last_max_err = max_err; } From f6931ade6f910c8ae77de2417e74091d20169992 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 16 Dec 2025 14:48:55 -0500 Subject: [PATCH 6/8] checks_impl: set rho to dive at conductive surface this fixes 2 things: 1. postprocessing no longer has to manually ignore surface charges, which required knowing BCs 2. perhaps more importantly, the print_diffs helper now prints correct indices --- src/libpsc/psc_checks/checks_impl.hxx | 30 ++++++++++++--------------- 1 file changed, 13 insertions(+), 17 deletions(-) diff --git a/src/libpsc/psc_checks/checks_impl.hxx b/src/libpsc/psc_checks/checks_impl.hxx index 071a6a9f1..0271bf126 100644 --- a/src/libpsc/psc_checks/checks_impl.hxx +++ b/src/libpsc/psc_checks/checks_impl.hxx @@ -144,34 +144,30 @@ public: auto rho = psc::mflds::interior(grid, item_rho(mprts)); auto dive = psc::mflds::interior(grid, item_dive(mflds)); - double max_err = 0.; for (int p = 0; p < grid.n_patches(); p++) { - int l[3] = {0, 0, 0}, r[3] = {0, 0, 0}; for (int d = 0; d < 3; d++) { + 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)) { - l[d] = 1; + 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); } } + } - auto patch_rho = - rho.view(_s(l[0], -r[0]), _s(l[1], -r[1]), _s(l[2], -r[2]), 0, p); - auto patch_dive = - dive.view(_s(l[0], -r[0]), _s(l[1], -r[1]), _s(l[2], -r[2]), 0, p); + double local_err = gt::norm_linf(dive - rho); - double patch_err = gt::norm_linf(patch_dive - patch_rho); - max_err = std::max(max_err, patch_err); + double max_err; + MPI_Allreduce(&local_err, &max_err, 1, MPI_DOUBLE, MPI_MAX, grid.comm()); - if (should_print_diffs(patch_err)) { - mpi_printf(grid.comm(), "gauss: rho -- div E\n"); - psc::helper::print_diff(patch_rho, patch_dive, err_threshold); - } + if (should_print_diffs(max_err)) { + mpi_printf(grid.comm(), "gauss: rho -- div E\n"); + psc::helper::print_diff(rho, dive, err_threshold); } - // find global max - double tmp = max_err; - MPI_Allreduce(&tmp, &max_err, 1, MPI_DOUBLE, MPI_MAX, grid.comm()); - if (should_print_max_err(max_err)) { mpi_printf(grid.comm(), "gauss: max_err = %g (thres %g)\n", max_err, err_threshold); From c3e0e51debcbd5e22e6e8554f0ae61f3e69bb259 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 16 Dec 2025 15:32:14 -0500 Subject: [PATCH 7/8] checks_impl: improve print timing --- src/libpsc/psc_checks/checks_impl.hxx | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/src/libpsc/psc_checks/checks_impl.hxx b/src/libpsc/psc_checks/checks_impl.hxx index 0271bf126..1b098fec0 100644 --- a/src/libpsc/psc_checks/checks_impl.hxx +++ b/src/libpsc/psc_checks/checks_impl.hxx @@ -85,11 +85,6 @@ public: psc::helper::print_diff(d_rho, -dt_divj, err_threshold); } - if (should_print_max_err(max_err)) { - mpi_printf(grid.comm(), "continuity: max_err = %g (thres %g)\n", max_err, - err_threshold); - } - if (should_dump(max_err)) { if (!writer_) { writer_.open("continuity"); @@ -101,6 +96,11 @@ public: MPI_Barrier(grid.comm()); } + if (should_print_max_err(max_err)) { + mpi_printf(grid.comm(), "continuity: max_err = %g (thres %g)\n", max_err, + err_threshold); + } + if (exit_on_failure && max_err >= err_threshold) { LOG_ERROR("exiting due to failed continuity check\n"); } @@ -168,11 +168,6 @@ public: psc::helper::print_diff(rho, dive, err_threshold); } - if (should_print_max_err(max_err)) { - mpi_printf(grid.comm(), "gauss: max_err = %g (thres %g)\n", max_err, - err_threshold); - } - if (should_dump(max_err)) { if (!writer_) { writer_.open("gauss"); @@ -181,6 +176,12 @@ public: writer_.write(rho, grid, "rho", {"rho"}); writer_.write(dive, grid, "dive", {"dive"}); writer_.end_step(); + MPI_Barrier(grid.comm()); + } + + if (should_print_max_err(max_err)) { + mpi_printf(grid.comm(), "gauss: max_err = %g (thres %g)\n", max_err, + err_threshold); } if (exit_on_failure && max_err >= err_threshold) { From a7e31e95e43a66a81129a1f60e8a4982d0f8b678 Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 22 Dec 2025 10:00:10 -0500 Subject: [PATCH 8/8] psc_bnd_fields_impl: fix format --- 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 7a368419b..1c670e5aa 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -644,7 +644,7 @@ struct BndFieldsNone : BndFieldsBase { using Mfields = MF; - void fill_ghosts_E(Mfields& mflds) {}; - void fill_ghosts_H(Mfields& mflds) {}; - void add_ghosts_J(Mfields& mflds) {}; + void fill_ghosts_E(Mfields& mflds){}; + void fill_ghosts_H(Mfields& mflds){}; + void add_ghosts_J(Mfields& mflds){}; };