diff --git a/src/include/checks.hxx b/src/include/checks.hxx deleted file mode 100644 index 50185dc47a..0000000000 --- a/src/include/checks.hxx +++ /dev/null @@ -1,19 +0,0 @@ - -#pragma once - -struct ChecksParams -{ - int continuity_every_step = - 0; // check charge continuity eqn every so many steps - double continuity_threshold = 1e-13; // acceptable error in continuity eqn - bool continuity_verbose = - true; // always print continuity error, even if acceptable - bool continuity_dump_always = - false; // always dump d_rho, div_j, even if acceptable - - int gauss_every_step = 0; // check Gauss's Law every so many steps - double gauss_threshold = 1e-13; // acceptable error in Gauss's Law - bool gauss_verbose = - true; // always print Gauss's Law error, even if acceptable - bool gauss_dump_always = false; // always dump E, div_rho, even if acceptable -}; diff --git a/src/include/checks_params.hxx b/src/include/checks_params.hxx new file mode 100644 index 0000000000..8b328f6bc6 --- /dev/null +++ b/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; +}; \ No newline at end of file diff --git a/src/include/psc.hxx b/src/include/psc.hxx index 4a25af1eb8..c96ba36c4c 100644 --- a/src/include/psc.hxx +++ b/src/include/psc.hxx @@ -8,7 +8,7 @@ #include #include "../libpsc/vpic/fields_item_vpic.hxx" -#include +#include #include #include @@ -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); diff --git a/src/libpsc/integrate.cxx b/src/libpsc/integrate.cxx index 2a279638f9..6f67ba3b91 100644 --- a/src/libpsc/integrate.cxx +++ b/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 #include diff --git a/src/libpsc/psc_checks/checks_impl.hxx b/src/libpsc/psc_checks/checks_impl.hxx index c513712909..33ddbe55eb 100644 --- a/src/libpsc/psc_checks/checks_impl.hxx +++ b/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 @@ -32,13 +32,15 @@ namespace checks // psc::checks::continuity: Charge Continuity template -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 -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,8 +159,8 @@ 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); } } @@ -164,12 +168,12 @@ public: 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,7 +210,7 @@ struct checks_order_2nd }; template -class ChecksCommon : public ChecksParams +class ChecksCommon { public: using Mparticles = MP; @@ -210,29 +218,12 @@ public: 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 - void continuity_after_particle_push(Mparticles& mprts, MfieldsState& mflds) - { - continuity_.after_particle_push(mprts, mflds); - } - - template - void gauss(Mparticles& mprts, MfieldsState& mflds) - { - gauss_(mprts, mflds); - } - -private: - psc::checks::continuity continuity_; - psc::checks::gauss gauss_; +public: + psc::checks::continuity continuity; + psc::checks::gauss gauss; }; template diff --git a/src/libpsc/tests/test_push_particles_2.cxx b/src/libpsc/tests/test_push_particles_2.cxx index 7dd682d250..ac5d300cb1 100644 --- a/src/libpsc/tests/test_push_particles_2.cxx +++ b/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)); diff --git a/src/psc_2d_shock.cxx b/src/psc_2d_shock.cxx index 443697208c..58eaeef73e 100644 --- a/src/psc_2d_shock.cxx +++ b/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 diff --git a/src/psc_bgk.cxx b/src/psc_bgk.cxx index 208aedbdb9..d6f901ca3d 100644 --- a/src/psc_bgk.cxx +++ b/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}; diff --git a/src/psc_flatfoil_yz.cxx b/src/psc_flatfoil_yz.cxx index 342e520100..1f2469e9da 100644 --- a/src/psc_flatfoil_yz.cxx +++ b/src/psc_flatfoil_yz.cxx @@ -478,24 +478,24 @@ void run() // -- Checks ChecksParams checks_params{}; #if CASE == CASE_2D_SMALL - checks_params.continuity_every_step = 1; - checks_params.continuity_dump_always = true; + checks_params.continuity.check_interval = 1; + checks_params.continuity.dump_always = true; #else - checks_params.continuity_every_step = 0; - checks_params.continuity_dump_always = false; + checks_params.continuity.check_interval = 0; + checks_params.continuity.dump_always = false; #endif - checks_params.continuity_threshold = 1e-4; - checks_params.continuity_verbose = true; + checks_params.continuity.err_threshold = 1e-4; + checks_params.continuity.print_max_err_always = true; #if CASE == CASE_2D_SMALL - checks_params.gauss_every_step = 1; - checks_params.gauss_dump_always = true; + checks_params.gauss.check_interval = 1; + checks_params.gauss.dump_always = true; #else - checks_params.gauss_every_step = 100; - checks_params.gauss_dump_always = false; + checks_params.gauss.check_interval = 100; + checks_params.gauss.dump_always = false; #endif - checks_params.gauss_threshold = 1e-4; - checks_params.gauss_verbose = true; + checks_params.gauss.err_threshold = 1e-4; + checks_params.gauss.print_max_err_always = true; Checks checks{grid, MPI_COMM_WORLD, checks_params}; diff --git a/src/psc_harris_yz.cxx b/src/psc_harris_yz.cxx index 669c83d9a2..0d1f27a126 100644 --- a/src/psc_harris_yz.cxx +++ b/src/psc_harris_yz.cxx @@ -385,15 +385,15 @@ void run() // -- Checks ChecksParams checks_params{}; - checks_params.continuity_every_step = 0; - checks_params.continuity_dump_always = false; - checks_params.continuity_threshold = 1e-4; - checks_params.continuity_verbose = true; - - checks_params.gauss_every_step = -100; - checks_params.gauss_dump_always = false; - checks_params.gauss_threshold = 1e-4; - checks_params.gauss_verbose = true; + checks_params.continuity.check_interval = 0; + checks_params.continuity.dump_always = false; + checks_params.continuity.err_threshold = 1e-4; + checks_params.continuity.print_max_err_always = true; + + checks_params.gauss.check_interval = -100; + checks_params.gauss.dump_always = false; + checks_params.gauss.err_threshold = 1e-4; + checks_params.gauss.print_max_err_always = true; Checks checks{grid, MPI_COMM_WORLD, checks_params}; diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index d90053dade..6a8dc48f9c 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -278,9 +278,9 @@ static void run(int argc, char** argv) // -- Checks ChecksParams checks_params{}; - checks_params.gauss_every_step = out_interval; - // checks_params.gauss_dump_always = true; - checks_params.gauss_threshold = 1e-5; + checks_params.gauss.check_interval = out_interval; + // checks_params.gauss.dump_always = true; + checks_params.gauss.err_threshold = 1e-5; Checks checks{grid, MPI_COMM_WORLD, checks_params}; diff --git a/src/psc_whistler.cxx b/src/psc_whistler.cxx index 04a5500009..793cd2f89e 100644 --- a/src/psc_whistler.cxx +++ b/src/psc_whistler.cxx @@ -292,9 +292,9 @@ void run() // -- Checks ChecksParams checks_params{}; - checks_params.continuity_every_step = 50; - checks_params.continuity_threshold = 1e-5; - checks_params.continuity_verbose = false; + checks_params.continuity.check_interval = 50; + checks_params.continuity.err_threshold = 1e-5; + checks_params.continuity.print_max_err_always = false; Checks checks{grid, MPI_COMM_WORLD, checks_params}; // -- Marder correction