Refactor gauss and continuity checks API

src/include/checks_params.hxx

@@ -0,0 +1,41 @@
+
+#pragma once
+
+struct CheckParams
+{
+  int check_interval = 0;           // number of steps per check
+  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 enabled() { return check_interval > 0; }
+
+  bool should_do_check(int timestep)
+  {
+    return enabled() && timestep % check_interval == 0;
+  }
+
+  bool should_print_diffs(double err) { return err > err_threshold; }
+
+  bool should_print_max_err(double max_err)
+  {
+    return print_max_err_always || max_err > err_threshold;
+  }
+
+  bool should_dump(double max_err)
+  {
+    return dump_always || max_err > err_threshold;
+  }
+};
+
+struct ContinuityCheckParams : CheckParams
+{};
+
+struct GaussCheckParams : CheckParams
+{};
+
+struct ChecksParams
+{
+  ContinuityCheckParams continuity;
+  GaussCheckParams gauss;
+};

src/include/psc.hxx

@@ -8,7 +8,7 @@
 #include <setup_particles.hxx>
 
 #include "../libpsc/vpic/fields_item_vpic.hxx"
-#include <checks.hxx>
+#include <checks_params.hxx>
 #include <output_particles.hxx>
 #include <push_particles.hxx>
 
@@ -489,11 +489,10 @@ struct Psc
     inject_particles();
     prof_stop(pr_inject_prts);
 
-    if (checks_.continuity_every_step > 0 &&
-        timestep % checks_.continuity_every_step == 0) {
+    if (checks_.continuity.should_do_check(timestep)) {
       mpi_printf(comm, "***** Checking continuity...\n");
       prof_start(pr_checks);
-      checks_.continuity_before_particle_push(mprts_);
+      checks_.continuity.before_particle_push(mprts_);
       prof_stop(pr_checks);
     }
 
@@ -558,10 +557,9 @@ struct Psc
     // state is now: x^{n+3/2}, p^{n+1}, E^{n+3/2}, B^{n+3/2}
 #endif
 
-    if (checks_.continuity_every_step > 0 &&
-        timestep % checks_.continuity_every_step == 0) {
+    if (checks_.continuity.should_do_check(timestep)) {
       prof_restart(pr_checks);
-      checks_.continuity_after_particle_push(mprts_, mflds_);
+      checks_.continuity.after_particle_push(mprts_, mflds_);
       prof_stop(pr_checks);
     }
 
@@ -575,8 +573,7 @@ struct Psc
       prof_stop(pr_marder);
     }
 
-    if (checks_.gauss_every_step > 0 &&
-        timestep % checks_.gauss_every_step == 0) {
+    if (checks_.gauss.should_do_check(timestep)) {
       prof_restart(pr_checks);
       checks_.gauss(mprts_, mflds_);
       prof_stop(pr_checks);

src/libpsc/integrate.cxx

@@ -6,7 +6,7 @@
 #include "sort.hxx"
 #include "collision.hxx"
 #include "bnd_particles.hxx"
-#include "checks.hxx"
+#include "checks_params.hxx"
 
 #include <mrc_common.h>
 #include <mrc_profile.h>

src/libpsc/psc_checks/checks_impl.hxx

@@ -3,7 +3,7 @@
 
 #include "fields.hxx"
 #include "fields_item.hxx"
-#include "checks.hxx"
+#include "checks_params.hxx"
 #include "../libpsc/psc_output_fields/fields_item_fields.hxx"
 #include "../libpsc/psc_output_fields/fields_item_moments_1st.hxx"
 #include <psc/helper.hxx>
@@ -32,13 +32,15 @@ namespace checks
 // psc::checks::continuity: Charge Continuity
 
 template <typename S, typename ITEM_RHO>
-class continuity : ChecksParams
+class continuity : public ContinuityCheckParams
 {
 public:
   using storage_type = S;
   using Item_rho = ITEM_RHO;
 
-  continuity(const ChecksParams& params) : ChecksParams(params) {}
+  continuity(const ContinuityCheckParams& params)
+    : ContinuityCheckParams(params)
+  {}
 
   // ----------------------------------------------------------------------
   // before_particle_push
@@ -47,8 +49,7 @@ public:
   void before_particle_push(const Mparticles& mprts)
   {
     const auto& grid = mprts.grid();
-    if (continuity_every_step <= 0 ||
-        grid.timestep() % continuity_every_step != 0) {
+    if (!should_do_check(grid.timestep())) {
       return;
     }
 
@@ -63,8 +64,7 @@ public:
   void after_particle_push(const Mparticles& mprts, MfieldsState& mflds)
   {
     const Grid_t& grid = mprts.grid();
-    if (continuity_every_step <= 0 ||
-        grid.timestep() % continuity_every_step != 0) {
+    if (!should_do_check(grid.timestep())) {
       return;
     }
 
@@ -81,16 +81,16 @@ public:
     double max_err;
     MPI_Allreduce(&local_err, &max_err, 1, MPI_DOUBLE, MPI_MAX, grid.comm());
 
-    if (max_err >= continuity_threshold) {
-      psc::helper::print_diff(d_rho, -dt_divj, continuity_threshold);
+    if (should_print_diffs(max_err)) {
+      psc::helper::print_diff(d_rho, -dt_divj, err_threshold);
     }
 
-    if (continuity_verbose || max_err >= continuity_threshold) {
+    if (should_print_max_err(max_err)) {
       mpi_printf(grid.comm(), "continuity: max_err = %g (thres %g)\n", max_err,
-                 continuity_threshold);
+                 err_threshold);
     }
 
-    if (continuity_dump_always || max_err >= continuity_threshold) {
+    if (should_dump(max_err)) {
       if (!writer_) {
         writer_.open("continuity");
       }
@@ -101,9 +101,13 @@ public:
       MPI_Barrier(grid.comm());
     }
 
-    assert(max_err < continuity_threshold);
+    assert(max_err < err_threshold);
+    last_max_err = max_err;
   }
 
+public:
+  double last_max_err;
+
 private:
   storage_type rho_m_;
   WriterDefault writer_;
@@ -113,13 +117,13 @@ private:
 // psc::checks::gauss: Gauss's Law div E = rho
 
 template <typename S, typename ITEM_RHO>
-class gauss : ChecksParams
+class gauss : public GaussCheckParams
 {
 public:
   using storage_type = S;
   using Item_rho = ITEM_RHO;
 
-  gauss(const ChecksParams& params) : ChecksParams(params) {}
+  gauss(const GaussCheckParams params) : GaussCheckParams(params) {}
 
   // ----------------------------------------------------------------------
   // operator()
@@ -128,7 +132,7 @@ public:
   void operator()(Mparticles& mprts, MfieldsState& mflds)
   {
     const auto& grid = mprts.grid();
-    if (gauss_every_step <= 0 || grid.timestep() % gauss_every_step != 0) {
+    if (!should_do_check(grid.timestep())) {
       return;
     }
 
@@ -155,21 +159,21 @@ public:
       double patch_err = gt::norm_linf(patch_dive - patch_rho);
       max_err = std::max(max_err, patch_err);
 
-      if (patch_err > gauss_threshold) {
-        psc::helper::print_diff(patch_rho, patch_dive, gauss_threshold);
+      if (should_print_diffs(patch_err)) {
+        psc::helper::print_diff(patch_rho, patch_dive, err_threshold);
       }
     }
 
     // find global max
     double tmp = max_err;
     MPI_Allreduce(&tmp, &max_err, 1, MPI_DOUBLE, MPI_MAX, grid.comm());
 
-    if (gauss_verbose || max_err >= gauss_threshold) {
+    if (should_print_max_err(max_err)) {
       mpi_printf(grid.comm(), "gauss: max_err = %g (thres %g)\n", max_err,
-                 gauss_threshold);
+                 err_threshold);
     }
 
-    if (gauss_dump_always || max_err >= gauss_threshold) {
+    if (should_dump(max_err)) {
       if (!writer_) {
         writer_.open("gauss");
       }
@@ -179,9 +183,13 @@ public:
       writer_.end_step();
     }
 
-    assert(max_err < gauss_threshold);
+    assert(max_err < err_threshold);
+    last_max_err = max_err;
   }
 
+public:
+  double last_max_err;
+
 private:
   WriterDefault writer_;
 };
@@ -202,37 +210,20 @@ struct checks_order_2nd
 };
 
 template <typename MP, typename S, typename ITEM_RHO>
-class ChecksCommon : public ChecksParams
+class ChecksCommon
 {
 public:
   using Mparticles = MP;
   using storage_type = S;
   using item_rho_type = ITEM_RHO;
 
   ChecksCommon(const Grid_t& grid, MPI_Comm comm, const ChecksParams& params)
-    : ChecksParams{params}, continuity_{params}, gauss_{params}
+    : continuity{params.continuity}, gauss{params.gauss}
   {}
 
-  void continuity_before_particle_push(Mparticles& mprts)
-  {
-    continuity_.before_particle_push(mprts);
-  }
-
-  template <typename MfieldsState>
-  void continuity_after_particle_push(Mparticles& mprts, MfieldsState& mflds)
-  {
-    continuity_.after_particle_push(mprts, mflds);
-  }
-
-  template <typename MfieldsState>
-  void gauss(Mparticles& mprts, MfieldsState& mflds)
-  {
-    gauss_(mprts, mflds);
-  }
-
-private:
-  psc::checks::continuity<storage_type, item_rho_type> continuity_;
-  psc::checks::gauss<storage_type, item_rho_type> gauss_;
+public:
+  psc::checks::continuity<storage_type, item_rho_type> continuity;
+  psc::checks::gauss<storage_type, item_rho_type> gauss;
 };
 
 template <typename MP, typename MF, typename ORDER, typename D>

src/libpsc/tests/test_push_particles_2.cxx

@@ -71,17 +71,17 @@ TYPED_TEST(PushParticlesTest, Accel)
   BndParticles bndp_{grid};
   Bnd bnd_{};
   ChecksParams checks_params{};
-  checks_params.continuity_every_step = 1;
-  checks_params.continuity_threshold = 1e-7;
-  checks_params.continuity_verbose = false;
+  checks_params.continuity.check_interval = 1;
+  checks_params.continuity.err_threshold = 1e-7;
+  checks_params.continuity.print_max_err_always = false;
   Checks checks_{grid, MPI_COMM_WORLD, checks_params};
   for (int n = 0; n < n_steps; n++) {
-    checks_.continuity_before_particle_push(mprts);
+    checks_.continuity.before_particle_push(mprts);
     pushp_.push_mprts(mprts, mflds);
     bndp_(mprts);
     bnd_.add_ghosts(mflds, JXI, JXI + 3);
     bnd_.fill_ghosts(mflds, JXI, JXI + 3);
-    checks_.continuity_after_particle_push(mprts, mflds);
+    checks_.continuity.after_particle_push(mprts, mflds);
 
     auto accessor = mprts.accessor();
     for (auto prt : accessor[0]) {
@@ -149,17 +149,17 @@ TYPED_TEST(PushParticlesTest, Cyclo)
   BndParticles bndp_{grid};
   Bnd bnd_{};
   ChecksParams checks_params{};
-  checks_params.continuity_every_step = 1;
-  checks_params.continuity_threshold = 1e-7;
-  checks_params.continuity_verbose = false;
+  checks_params.continuity.check_interval = 1;
+  checks_params.continuity.err_threshold = 1e-7;
+  checks_params.continuity.print_max_err_always = false;
   Checks checks_{grid, MPI_COMM_WORLD, checks_params};
   for (int n = 0; n < n_steps; n++) {
-    checks_.continuity_before_particle_push(mprts);
+    checks_.continuity.before_particle_push(mprts);
     pushp_.push_mprts(mprts, mflds);
     bndp_(mprts);
     bnd_.add_ghosts(mflds, JXI, JXI + 3);
     bnd_.fill_ghosts(mflds, JXI, JXI + 3);
-    checks_.continuity_after_particle_push(mprts, mflds);
+    checks_.continuity.after_particle_push(mprts, mflds);
 
     double ux = (cos(2 * M_PI * (0.125 * n_steps - (n + 1)) / (double)n_steps) /
                  cos(2 * M_PI * (0.125 * n_steps) / (double)n_steps));

src/psc_2d_shock.cxx

@@ -439,14 +439,14 @@ void run()
 
   // -- Checks
   ChecksParams checks_params{};
-  checks_params.continuity_every_step = 100;
-  checks_params.continuity_threshold = 1e-4;
-  checks_params.continuity_verbose = true;
-  checks_params.continuity_dump_always = false;
-  checks_params.gauss_every_step = 100;
-  checks_params.gauss_threshold = 1e-4;
-  checks_params.gauss_verbose = true;
-  checks_params.gauss_dump_always = false;
+  checks_params.continuity.check_interval = 100;
+  checks_params.continuity.err_threshold = 1e-4;
+  checks_params.continuity.print_max_err_always = true;
+  checks_params.continuity.dump_always = false;
+  checks_params.gauss.check_interval = 100;
+  checks_params.gauss.err_threshold = 1e-4;
+  checks_params.gauss.print_max_err_always = true;
+  checks_params.gauss.dump_always = false;
   Checks checks{grid, MPI_COMM_WORLD, checks_params};
 
   // -- Marder correction

src/psc_bgk.cxx

@@ -510,9 +510,9 @@ static void run(int argc, char** argv)
 
   // -- Checks
   ChecksParams checks_params{};
-  checks_params.gauss_every_step = g.gauss_every;
-  // checks_params.gauss_dump_always = true;
-  checks_params.gauss_threshold = 1e-5;
+  checks_params.gauss.check_interval = g.gauss_every;
+  // checks_params.gauss.dump_always = true;
+  checks_params.gauss.err_threshold = 1e-5;
 
   Checks checks{grid, MPI_COMM_WORLD, checks_params};