diff --git a/src/include/DiagEnergies.h b/src/include/DiagEnergies.h index 67a5037b62..bbd72dc6c8 100644 --- a/src/include/DiagEnergies.h +++ b/src/include/DiagEnergies.h @@ -3,23 +3,29 @@ #include "DiagEnergiesField.h" #include "DiagEnergiesParticle.h" +#include "diagnostic_base.hxx" #include "psc.h" -class DiagEnergies +template +class DiagEnergies : public DiagnosticBase { public: DiagEnergies(); DiagEnergies(MPI_Comm comm, int interval); - template + void perform_diagnostic(Mparticles& mprts, MfieldsState& mflds) + { + (*this)(mprts, mflds); + } + void operator()(Mparticles& mprts, MfieldsState& mflds); private: template static std::string legend(const Item& item); - template + template void write_one(const Item& item, Mparticles& mprts, MfieldsState& mflds); private: diff --git a/src/include/DiagEnergies.inl b/src/include/DiagEnergies.inl index 7a81a0ecd3..25c660ab85 100644 --- a/src/include/DiagEnergies.inl +++ b/src/include/DiagEnergies.inl @@ -2,19 +2,21 @@ namespace { -void fclose_helper(FILE* fp) -{ - ::fclose(fp); -} +void fclose_helper(FILE* fp) { ::fclose(fp); } } // namespace // ---------------------------------------------------------------------- // DiagEnergies ctors -inline DiagEnergies::DiagEnergies() : file_{nullptr, fclose_helper} {} +template +inline DiagEnergies::DiagEnergies() + : file_{nullptr, fclose_helper} +{} -inline DiagEnergies::DiagEnergies(MPI_Comm comm, int interval) +template +inline DiagEnergies::DiagEnergies(MPI_Comm comm, + int interval) : comm_{comm}, interval_{interval}, file_{nullptr, fclose_helper} { MPI_Comm_rank(comm_, &rank_); @@ -32,7 +34,8 @@ inline DiagEnergies::DiagEnergies(MPI_Comm comm, int interval) // DiagEnergies::operator() template -inline void DiagEnergies::operator()(Mparticles& mprts, MfieldsState& mflds) +inline void DiagEnergies::operator()( + Mparticles& mprts, MfieldsState& mflds) { const auto& grid = mprts.grid(); @@ -56,8 +59,10 @@ inline void DiagEnergies::operator()(Mparticles& mprts, MfieldsState& mflds) // ---------------------------------------------------------------------- // legend +template template -inline std::string DiagEnergies::legend(const Item& item) +inline std::string DiagEnergies::legend( + const Item& item) { std::string s; for (auto& name : item.names()) { @@ -69,9 +74,10 @@ inline std::string DiagEnergies::legend(const Item& item) // ---------------------------------------------------------------------- // write_one -template -inline void DiagEnergies::write_one(const Item& item, Mparticles& mprts, - MfieldsState& mflds) +template +template +inline void DiagEnergies::write_one( + const Item& item, Mparticles& mprts, MfieldsState& mflds) { auto vals = item(mprts, mflds); diff --git a/src/include/DiagnosticsDefault.h b/src/include/DiagnosticsDefault.h deleted file mode 100644 index b5071d01eb..0000000000 --- a/src/include/DiagnosticsDefault.h +++ /dev/null @@ -1,49 +0,0 @@ - -#pragma once - -#include "DiagEnergies.h" - -// ====================================================================== -// DiagnosticsDefault - -template -class DiagnosticsDefault -{ -public: - DiagnosticsDefault(OutputFields& outf, OutputParticles& outp, - OutputEnergies& oute) - : outf_{outf}, outp_{outp}, oute_{oute} - {} - - template - void operator()(Mparticles& mprts, MfieldsState& mflds) - { - psc_stats_start(st_time_output); -#ifdef VPIC -#if 0 - TIC user_diagnostics(); TOC(user_diagnostics, 1); -#endif -#else - // FIXME - outf_(mflds, mprts); -#endif - outp_(mprts); - oute_(mprts, mflds); - psc_stats_stop(st_time_output); - } - -private: - OutputFields& outf_; - OutputParticles& outp_; - OutputEnergies& oute_; -}; - -template -DiagnosticsDefault -makeDiagnosticsDefault(OutputFields& outf, OutputParticles& outp, - OutputEnergies& oute) -{ - return {outf, outp, oute}; -} diff --git a/src/include/OutputFieldsDefault.h b/src/include/OutputFieldsDefault.h index 0d76bd5377..eafc7b5883 100644 --- a/src/include/OutputFieldsDefault.h +++ b/src/include/OutputFieldsDefault.h @@ -1,6 +1,7 @@ #pragma once +#include "diagnostic_base.hxx" #include "../libpsc/psc_output_fields/fields_item_fields.hxx" #include "../libpsc/psc_output_fields/fields_item_moments_1st.hxx" #ifdef USE_CUDA @@ -242,7 +243,7 @@ struct OutputFieldsParams template -class OutputFields +class OutputFields : public DiagnosticBase { public: // ---------------------------------------------------------------------- @@ -252,6 +253,11 @@ class OutputFields : fields{grid, prm.fields, ""}, moments{grid, prm.moments, "_moments"} {} + void perform_diagnostic(Mparticles& mprts, MfieldsState& mflds) override + { + (*this)(mflds, mprts); + } + // ---------------------------------------------------------------------- // operator() diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx index 72f9dc4eb6..3f86324573 100644 --- a/src/include/boundary_injector.hxx +++ b/src/include/boundary_injector.hxx @@ -11,6 +11,7 @@ #include "setup_particles.hxx" #include "kg/VecRange.hxx" #include "../libpsc/psc_push_particles/inc_push.cxx" +#include "injector_base.hxx" /// @brief A particle generator for use with @ref BoundaryInjector. Samples /// particles from a (possibly shifted) Maxwellian distribution. @@ -54,7 +55,7 @@ private: /// @brief Injects particles on a given boundary, sampling from a given particle /// generator. For precise control over multiple particle species, use one -/// BoundaryInjector per species, combined with @ref CompositeInjector. +/// BoundaryInjector per species. /// @tparam PARTICLE_GENERATOR a type that defines `get(min_pos, pos_range)` and /// returns an injectable particle within that range of positions (usually a /// grid cell); see @ref ParticleGeneratorMaxwellian @@ -62,6 +63,8 @@ private: /// `MfieldsState`, `Current`, `real_t`, etc. template class BoundaryInjector + : public InjectorBase { static const int INJECT_DIM_IDX_ = 1; @@ -88,7 +91,7 @@ public: /// the given ParticleGenerator. /// /// Some of these limitations may be removed in the future. - void operator()(Mparticles& mprts, MfieldsState& mflds) + void inject(Mparticles& mprts, MfieldsState& mflds) override { static_assert(INJECT_DIM_IDX_ == 1, "only injection at lower bound of y is supported"); diff --git a/src/include/composite_injector.hxx b/src/include/composite_injector.hxx deleted file mode 100644 index 105c605bff..0000000000 --- a/src/include/composite_injector.hxx +++ /dev/null @@ -1,48 +0,0 @@ -#pragma once - -/// @brief An injector comprising two other injectors. Enables multiple -/// injectors in a single psc case. Useful for injecting multiple species with -/// @ref BoundaryInjector. -/// @tparam INJECTOR_1 first type of injector -/// @tparam INJECTOR_2 second type of injector -template -class CompositeInjector -{ -public: - using Mparticles = typename INJECTOR_1::Mparticles; - using MfieldsState = typename INJECTOR_1::MfieldsState; - - CompositeInjector(INJECTOR_1 injector_1, INJECTOR_2 injector_2) - : injector_1{injector_1}, injector_2{injector_2} - {} - - /// @brief Calls both member injectors in order. - /// @param mprts particles - /// @param mflds fields - void operator()(Mparticles& mprts, MfieldsState& mflds) - { - injector_1(mprts, mflds); - injector_2(mprts, mflds); - } - -private: - INJECTOR_1 injector_1; - INJECTOR_2 injector_2; -}; - -/// @brief Helper method to deduce the type of a composite injector, enabling -/// e.g. `auto composite_injector = make_composite(injector_1, injector_2);` -/// -/// Wouldn't be necessary in C++17; see -/// https://isocpp.github.io/CppCoreGuidelines/CppCoreGuidelines#Rt-deduce. -/// @tparam INJECTOR_1 first type of injector -/// @tparam INJECTOR_2 second type of injector -/// @param injector_1 first injector -/// @param injector_2 second injector -/// @return -template -CompositeInjector make_composite(INJECTOR_1 injector_1, - INJECTOR_2 injector_2) -{ - return CompositeInjector(injector_1, injector_2); -} diff --git a/src/include/diagnostic_base.hxx b/src/include/diagnostic_base.hxx new file mode 100644 index 0000000000..df213f66b9 --- /dev/null +++ b/src/include/diagnostic_base.hxx @@ -0,0 +1,29 @@ +#pragma once + +template +struct DiagnosticBase +{ + virtual void perform_diagnostic(Mparticles& mprts, MfieldsState& mflds) = 0; +}; + +template +struct ParticleDiagnosticBase +{ + virtual void perform_diagnostic(Mparticles& mprts) = 0; +}; + +template +struct DiagnosticFromLambda : public DiagnosticBase +{ + DiagnosticFromLambda(std::function lambda) + : lambda{lambda} + {} + + void perform_diagnostic(Mparticles& mprts, MfieldsState& mflds) override + { + return lambda(mprts, mflds); + } + +private: + std::function lambda; +}; diff --git a/src/include/external_current_base.hxx b/src/include/external_current_base.hxx new file mode 100644 index 0000000000..48242f79e1 --- /dev/null +++ b/src/include/external_current_base.hxx @@ -0,0 +1,24 @@ +#pragma once + +template +struct ExternalCurrentBase +{ + using MfieldsState = MFIELDS_STATE; + + virtual void inject_current(MfieldsState& mflds) = 0; +}; + +template +struct ExternalCurrentFromLambda : ExternalCurrentBase +{ + using MfieldsState = MFIELDS_STATE; + + ExternalCurrentFromLambda(std::function lambda) + : lambda{lambda} + {} + + void inject_current(MfieldsState& mflds) override { return lambda(mflds); } + +private: + std::function lambda; +}; diff --git a/src/include/injector_base.hxx b/src/include/injector_base.hxx new file mode 100644 index 0000000000..5eddfe9c2a --- /dev/null +++ b/src/include/injector_base.hxx @@ -0,0 +1,29 @@ +#pragma once + +template +struct InjectorBase +{ + using Mparticles = MPARTICLES; + using MfieldsState = MFIELDS_STATE; + + virtual void inject(Mparticles& mprts, MfieldsState& mflds) = 0; +}; + +template +struct InjectFromLambda : InjectorBase +{ + using Mparticles = MPARTICLES; + using MfieldsState = MFIELDS_STATE; + + InjectFromLambda(std::function lambda) + : lambda{lambda} + {} + + void inject(Mparticles& mprts, MfieldsState& mflds) override + { + return lambda(mprts, mflds); + } + +private: + std::function lambda; +}; diff --git a/src/include/psc.hxx b/src/include/psc.hxx index 4b8a567f07..e275ae7fe2 100644 --- a/src/include/psc.hxx +++ b/src/include/psc.hxx @@ -8,6 +8,9 @@ #include #include "../libpsc/vpic/fields_item_vpic.hxx" +#include "diagnostic_base.hxx" +#include "injector_base.hxx" +#include "external_current_base.hxx" #include #include #include @@ -99,8 +102,7 @@ inline double courant_length(const Grid_t::Domain& domain) // ====================================================================== // Psc -template +template struct Psc { using Mparticles = typename PscConfig::Mparticles; @@ -116,14 +118,17 @@ struct Psc using BndFields = typename PscConfig::BndFields; using BndParticles = typename PscConfig::BndParticles; using Dim = typename PscConfig::Dim; + using DiagnosticBaseT = DiagnosticBase; + using ParticleDiagnosticBaseT = ParticleDiagnosticBase; + using InjectorBaseT = InjectorBase; + using ExternalCurrentBaseT = ExternalCurrentBase; // ---------------------------------------------------------------------- // ctor Psc(const PscParams& params, Grid_t& grid, MfieldsState& mflds, Mparticles& mprts, Balance& balance, Collision& collision, Checks& checks, - Marder& marder, Diagnostics& diagnostics, - InjectParticles& inject_particles, ExtCurrent& ext_current) + Marder& marder) : p_{params}, grid_{&grid}, mflds_{mflds}, @@ -133,9 +138,6 @@ struct Psc checks_{checks}, marder_{marder}, bndp_{grid}, - diagnostics_{diagnostics}, - inject_particles_{inject_particles}, - ext_current_{ext_current}, checkpointing_{params.write_checkpoint_every_step} { time_start_ = MPI_Wtime(); @@ -156,7 +158,40 @@ struct Psc #endif initialize_stats(); - initialize(); + } + + // ---------------------------------------------------------------------- + // API for modifying various internal components + + void add_injector(InjectorBaseT* injector) + { + if (injector) { + injectors_.push_back(injector); + } + } + + void add_external_current(ExternalCurrentBaseT* external_current) + { + if (external_current) { + external_currents_.push_back(external_current); + } + } + + void add_diagnostic(DiagnosticBaseT* diagnostic) + { + if (diagnostic) { + diagnostics_.push_back(diagnostic); + } + } + + void add_diagnostic(ParticleDiagnosticBaseT* diagnostic) + { + if (diagnostic) { + diagnostics_.push_back(new DiagnosticFromLambda( + [=](Mparticles& mprts, MfieldsState& mflds) { + return diagnostic->perform_diagnostic(mprts); + })); + } } // ---------------------------------------------------------------------- @@ -178,9 +213,9 @@ struct Psc } // ---------------------------------------------------------------------- - // initialize + // pre_first_step - void initialize() + void pre_first_step() { bndf_.fill_ghosts_H(mflds_); bnd_.fill_ghosts(mflds_, HX, HX + 3); @@ -198,16 +233,6 @@ struct Psc mpi_printf(grid().comm(), "Checking gauss.\n"); checks_.gauss(mprts_, mflds_); } - - // initial output / stats - mpi_printf(grid().comm(), "Performing initial diagnostics.\n"); - diagnostics(); - - psc_stats_val[st_nr_particles] = mprts_.size(); - - print_status(); - - mpi_printf(grid().comm(), "Initialization complete.\n"); } // ---------------------------------------------------------------------- @@ -220,6 +245,18 @@ struct Psc pr = prof_register("psc_step", 1., 0, 0); } + pre_first_step(); + + // initial output / stats + mpi_printf(grid().comm(), "Performing initial diagnostics.\n"); + perform_diagnostics(); + + psc_stats_val[st_nr_particles] = mprts_.size(); + + print_status(); + + mpi_printf(grid().comm(), "Initialization complete.\n"); + mpi_printf(grid().comm(), "*** Advancing\n"); double elapsed = MPI_Wtime(); @@ -240,7 +277,7 @@ struct Psc step(); - diagnostics(); + perform_diagnostics(); psc_stats_stop(st_time_step); prof_stop(pr); @@ -343,12 +380,16 @@ struct Psc // === particle injection prof_start(pr_inject_prts); - inject_particles_(mprts_, mflds_); + for (auto injector : injectors_) { + injector->inject(mprts_, mflds_); + } prof_stop(pr_inject_prts); // === external current prof_start(pr_external_current); - this->ext_current_(grid(), mflds_); + for (auto external_current : external_currents_) { + external_current->inject_current(mflds_); + } prof_stop(pr_external_current); mpi_printf(comm, "***** Bnd particles...\n"); @@ -454,9 +495,16 @@ private: } // ---------------------------------------------------------------------- - // diagnostics + // perform_diagnostics - void diagnostics() { diagnostics_(mprts_, mflds_); } + void perform_diagnostics() + { + psc_stats_start(st_time_output); + for (auto diagnostic : diagnostics_) { + diagnostic->perform_diagnostic(mprts_, mflds_); + } + psc_stats_stop(st_time_output); + } // ---------------------------------------------------------------------- // print_status @@ -484,9 +532,9 @@ protected: Collision& collision_; Checks& checks_; Marder& marder_; - Diagnostics& diagnostics_; - InjectParticles& inject_particles_; - ExtCurrent& ext_current_; + std::vector diagnostics_; + std::vector injectors_; + std::vector external_currents_; Sort sort_; PushParticles pushp_; @@ -507,59 +555,18 @@ protected: int st_time_step; }; -// ====================================================================== -// InjectParticlesNone - -class InjectParticlesNone -{ -public: - template - void operator()(Mparticles& mprts, MfieldsState& mflds) - {} -}; - -namespace -{ - -InjectParticlesNone injectParticlesNone; - -} // namespace - -// ====================================================================== -// ExtCurrentNone - -class ExtCurrentNone -{ -public: - template - void operator()(const Grid_t& grid, MfieldsState& mflds) - {} -}; - -namespace -{ - -ExtCurrentNone extCurrentNone; - -} // namespace // ====================================================================== // makePscIntegrator template -Psc makePscIntegrator( - const PscParams& params, Grid_t& grid, MfieldsState& mflds, Mparticles& mprts, - Balance& balance, Collision& collision, Checks& checks, Marder& marder, - Diagnostics& diagnostics, - InjectParticles& inject_particles = injectParticlesNone, - ExtCurrent& ext_current = extCurrentNone) + typename Marder> +Psc makePscIntegrator(const PscParams& params, Grid_t& grid, + MfieldsState& mflds, Mparticles& mprts, + Balance& balance, Collision& collision, + Checks& checks, Marder& marder) { - return {params, grid, mflds, mprts, balance, - collision, checks, marder, diagnostics, inject_particles, - ext_current}; + return {params, grid, mflds, mprts, balance, collision, checks, marder}; } // ====================================================================== @@ -655,4 +662,4 @@ void vpic_create_diagnostics(int interval) {} // ---------------------------------------------------------------------- // vpic_setup_diagnostics -void vpic_setup_diagnostics() {} +void vpic_setup_perform_diagnostics() {} diff --git a/src/libpsc/psc_output_particles/output_particles_ascii_impl.hxx b/src/libpsc/psc_output_particles/output_particles_ascii_impl.hxx index 406825fab3..760005a930 100644 --- a/src/libpsc/psc_output_particles/output_particles_ascii_impl.hxx +++ b/src/libpsc/psc_output_particles/output_particles_ascii_impl.hxx @@ -1,18 +1,21 @@ - +#include "diagnostic_base.hxx" #include "output_particles.hxx" +template struct OutputParticlesAscii : OutputParticlesParams , OutputParticlesBase + , public ParticleDiagnosticBase { OutputParticlesAscii(const Grid_t& grid, const OutputParticlesParams& params) : OutputParticlesParams(params), comm_{grid.comm()} {} + void perform_diagnostic(Mparticles& mprts) override { (*this)(mprts); } + // ---------------------------------------------------------------------- // operator() - template void operator()(Mparticles& mprts) { const auto& grid = mprts.grid(); diff --git a/src/libpsc/psc_output_particles/output_particles_hdf5_impl.hxx b/src/libpsc/psc_output_particles/output_particles_hdf5_impl.hxx index 1417bfd676..a698a0b7da 100644 --- a/src/libpsc/psc_output_particles/output_particles_hdf5_impl.hxx +++ b/src/libpsc/psc_output_particles/output_particles_hdf5_impl.hxx @@ -8,6 +8,7 @@ #include "grid.hxx" #include "output_particles.hxx" +#include "diagnostic_base.hxx" #include "psc_particles_single.h" #include "../libpsc/vpic/mparticles_vpic.hxx" @@ -758,8 +759,10 @@ private: } // namespace detail -template -class OutputParticlesHdf5 : OutputParticlesBase +template +class OutputParticlesHdf5 + : OutputParticlesBase + , public ParticleDiagnosticBase { static OutputParticlesParams adjust_params( const OutputParticlesParams& params_in, const Grid_t& grid) @@ -781,7 +784,8 @@ public: writer_{params_, grid.kinds, grid.comm()} {} - template + void perform_diagnostic(Mparticles& mprts) override { (*this)(mprts); } + void operator()(Mparticles& mprts) { const auto& grid = mprts.grid(); diff --git a/src/libpsc/tests/test_boundary_injector.cxx b/src/libpsc/tests/test_boundary_injector.cxx index ef51fd70f1..2d8b1e93f2 100644 --- a/src/libpsc/tests/test_boundary_injector.cxx +++ b/src/libpsc/tests/test_boundary_injector.cxx @@ -3,10 +3,8 @@ #include "test_common.hxx" #include "boundary_injector.hxx" -#include "composite_injector.hxx" #include "psc.hxx" -#include "DiagnosticsDefault.h" #include "OutputFieldsDefault.h" #include "../psc_config.hxx" @@ -137,18 +135,12 @@ TEST(BoundaryInjectorTest, Integration1Particle) Collision collision{grid, 0, 0.1}; Marder marder(grid, 0.9, 3, false); - OutputFields outf{grid, {}}; - OutputParticles outp{grid, {}}; - DiagEnergies oute{grid.comm(), 0}; - auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); - - auto inject_particles = - BoundaryInjector{ - ParticleGenerator(1, 1), grid}; - auto psc = makePscIntegrator(psc_params, grid, mflds, mprts, - balance, collision, checks, marder, - diagnostics, inject_particles); + balance, collision, checks, marder); + + psc.add_injector( + new BoundaryInjector( + ParticleGenerator(1, 1), grid)); // ---------------------------------------------------------------------- // set up initial conditions @@ -164,6 +156,7 @@ TEST(BoundaryInjectorTest, Integration1Particle) ASSERT_EQ(prts.size(), 0); + psc.pre_first_step(); for (; grid.timestep_ < psc_params.nmax;) { psc.step(); @@ -200,18 +193,12 @@ TEST(BoundaryInjectorTest, IntegrationManyParticles) Collision collision{grid, 0, 0.1}; Marder marder(grid, 0.9, 3, false); - OutputFields outf{grid, {}}; - OutputParticles outp{grid, {}}; - DiagEnergies oute{grid.comm(), 0}; - auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); - - auto inject_particles = - BoundaryInjector{ - ParticleGenerator(-1, 1), grid}; - auto psc = makePscIntegrator(psc_params, grid, mflds, mprts, - balance, collision, checks, marder, - diagnostics, inject_particles); + balance, collision, checks, marder); + + psc.add_injector( + new BoundaryInjector( + ParticleGenerator(-1, 1), grid)); // ---------------------------------------------------------------------- // set up initial conditions @@ -227,6 +214,7 @@ TEST(BoundaryInjectorTest, IntegrationManyParticles) ASSERT_EQ(prts.size(), 0); + psc.pre_first_step(); for (; grid.timestep_ < psc_params.nmax;) { psc.step(); @@ -263,11 +251,6 @@ TEST(BoundaryInjectorTest, IntegrationManySpecies) Collision collision{grid, 0, 0.1}; Marder marder(grid, 0.9, 3, false); - OutputFields outf{grid, {}}; - OutputParticles outp{grid, {}}; - DiagEnergies oute{grid.comm(), 0}; - auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); - auto inject_electrons = BoundaryInjector{ ParticleGenerator(-1, 0), grid}; @@ -275,11 +258,11 @@ TEST(BoundaryInjectorTest, IntegrationManySpecies) BoundaryInjector{ ParticleGenerator(-1, 1), grid}; - auto inject = make_composite(inject_electrons, inject_ions); - auto psc = makePscIntegrator(psc_params, grid, mflds, mprts, - balance, collision, checks, marder, - diagnostics, inject); + balance, collision, checks, marder); + + psc.add_injector(&inject_ions); + psc.add_injector(&inject_electrons); // ---------------------------------------------------------------------- // set up initial conditions @@ -295,6 +278,7 @@ TEST(BoundaryInjectorTest, IntegrationManySpecies) ASSERT_EQ(prts.size(), 0); + psc.pre_first_step(); for (; grid.timestep_ < psc_params.nmax;) { psc.step(); diff --git a/src/libpsc/tests/test_open_bcs_integration.cxx b/src/libpsc/tests/test_open_bcs_integration.cxx index 92e36fcde3..e887775629 100644 --- a/src/libpsc/tests/test_open_bcs_integration.cxx +++ b/src/libpsc/tests/test_open_bcs_integration.cxx @@ -4,7 +4,6 @@ #include #include -#include "DiagnosticsDefault.h" #include "OutputFieldsDefault.h" #include "add_ghosts_reflecting.hxx" #include "../psc_config.hxx" @@ -107,14 +106,8 @@ TEST(OpenBcsTest, IntegrationY) Collision collision{grid, 0, 0.1}; Marder marder(grid, 0.9, 3, false); - OutputFields outf{grid, {}}; - OutputParticles outp{grid, {}}; - DiagEnergies oute{grid.comm(), 0}; - auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); - - auto psc = - makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, balance, - collision, checks, marder, diagnostics); + auto psc = makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, + balance, collision, checks, marder); // ---------------------------------------------------------------------- // set up initial conditions @@ -148,6 +141,7 @@ TEST(OpenBcsTest, IntegrationY) ASSERT_EQ(prts[0].q(), -1.0); ASSERT_EQ(prts[1].q(), 1.0); + psc.pre_first_step(); for (; grid.timestep_ < psc_params.nmax;) { psc.step(); ASSERT_LT(checks.continuity.last_max_err, checks.continuity.err_threshold); diff --git a/src/libpsc/tests/test_output_particles.cxx b/src/libpsc/tests/test_output_particles.cxx index ff68aa0777..19e5523c39 100644 --- a/src/libpsc/tests/test_output_particles.cxx +++ b/src/libpsc/tests/test_output_particles.cxx @@ -65,11 +65,13 @@ struct OutputParticlesTest : ::testing::Test }; using OutputParticlesTestTypes = ::testing::Types< - Config + Config> #ifdef H5_HAVE_PARALLEL , - Config>, - Config> + Config>, + Config> #endif >; diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index e16721eb1f..0da074977b 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -4,7 +4,6 @@ #include #include -#include "DiagnosticsDefault.h" #include "OutputFieldsDefault.h" #include "add_ghosts_reflecting.hxx" #include "../psc_config.hxx" diff --git a/src/libpsc/tests/test_reflective_bcs_integration.cxx b/src/libpsc/tests/test_reflective_bcs_integration.cxx index 5068be2563..c04c5d20cc 100644 --- a/src/libpsc/tests/test_reflective_bcs_integration.cxx +++ b/src/libpsc/tests/test_reflective_bcs_integration.cxx @@ -4,7 +4,6 @@ #include #include -#include "DiagnosticsDefault.h" #include "OutputFieldsDefault.h" #include "add_ghosts_reflecting.hxx" #include "../psc_config.hxx" @@ -97,14 +96,8 @@ TEST(ReflectiveBcsTest, IntegrationY) Collision collision{grid, 0, 0.1}; Marder marder(grid, 0.9, 3, false); - OutputFields outf{grid, {}}; - OutputParticles outp{grid, {}}; - DiagEnergies oute{grid.comm(), 0}; - auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); - - auto psc = - makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, balance, - collision, checks, marder, diagnostics); + auto psc = makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, + balance, collision, checks, marder); // ---------------------------------------------------------------------- // set up initial conditions @@ -141,6 +134,7 @@ TEST(ReflectiveBcsTest, IntegrationY) bool about_to_reflect = false; bool reflected = false; + psc.pre_first_step(); for (; grid.timestep_ < psc_params.nmax;) { about_to_reflect = prts[0].x()[1] < grid.domain.dx[1] && prts[0].u()[1] < 0.0; @@ -184,14 +178,8 @@ TEST(ReflectiveBcsTest, IntegrationZ) Collision collision{grid, 0, 0.1}; Marder marder(grid, 0.9, 3, false); - OutputFields outf{grid, {}}; - OutputParticles outp{grid, {}}; - DiagEnergies oute{grid.comm(), 0}; - auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); - - auto psc = - makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, balance, - collision, checks, marder, diagnostics); + auto psc = makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, + balance, collision, checks, marder); // ---------------------------------------------------------------------- // set up initial conditions @@ -228,6 +216,7 @@ TEST(ReflectiveBcsTest, IntegrationZ) bool about_to_reflect = false; bool reflected = false; + psc.pre_first_step(); for (; grid.timestep_ < psc_params.nmax;) { about_to_reflect = prts[0].x()[2] < grid.domain.dx[2] && prts[0].u()[2] < 0.0; diff --git a/src/psc_2d_shock.cxx b/src/psc_2d_shock.cxx index 97d72c4d1c..8a66e258a1 100644 --- a/src/psc_2d_shock.cxx +++ b/src/psc_2d_shock.cxx @@ -3,7 +3,6 @@ #include #include -#include "DiagnosticsDefault.h" #include "OutputFieldsDefault.h" #include "psc_config.hxx" @@ -480,9 +479,7 @@ void run() OutputParticles outp{grid, outp_params}; int oute_interval = -100; - DiagEnergies oute{grid.comm(), oute_interval}; - - auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); + DiagEnergies oute{grid.comm(), oute_interval}; // ---------------------------------------------------------------------- // Set up objects specific to the flatfoil case @@ -609,8 +606,13 @@ void run() // hand off to PscIntegrator to run the simulation auto psc = makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, - balance, collision, checks, marder, - diagnostics, lf_inject_heat); + balance, collision, checks, marder); + + psc.add_diagnostic(&outf); + psc.add_diagnostic(&outp); + psc.add_diagnostic(&oute); + psc.add_injector( + new InjectFromLambda(lf_inject_heat)); psc.integrate(); } diff --git a/src/psc_bgk.cxx b/src/psc_bgk.cxx index 0141ec7c24..4b39ce6b76 100644 --- a/src/psc_bgk.cxx +++ b/src/psc_bgk.cxx @@ -3,7 +3,6 @@ #include #include -#include "DiagnosticsDefault.h" #include "OutputFieldsDefault.h" #include "psc_config.hxx" #include "input_params.hxx" @@ -534,9 +533,7 @@ static void run(int argc, char** argv) OutputParticles outp{grid, outp_params}; int oute_interval = -100; - DiagEnergies oute{grid.comm(), oute_interval}; - - auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); + DiagEnergies oute{grid.comm(), oute_interval}; // ---------------------------------------------------------------------- // Set up initial conditions @@ -556,9 +553,12 @@ static void run(int argc, char** argv) // ---------------------------------------------------------------------- // Hand off to PscIntegrator to run the simulation - auto psc = - makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, balance, - collision, checks, marder, diagnostics); + auto psc = makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, + balance, collision, checks, marder); + + psc.add_diagnostic(&outf); + psc.add_diagnostic(&outp); + psc.add_diagnostic(&oute); psc.integrate(); } diff --git a/src/psc_bubble_yz.cxx b/src/psc_bubble_yz.cxx index 93c67ca7f3..720789aef3 100644 --- a/src/psc_bubble_yz.cxx +++ b/src/psc_bubble_yz.cxx @@ -3,7 +3,6 @@ #include #include -#include "DiagnosticsDefault.h" #include "OutputFieldsDefault.h" #include "psc_config.hxx" @@ -335,9 +334,7 @@ static void run() OutputParticles outp{grid, outp_params}; int oute_interval = 100; - DiagEnergies oute{grid.comm(), oute_interval}; - - auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); + DiagEnergies oute{grid.comm(), oute_interval}; // ---------------------------------------------------------------------- // setup initial conditions @@ -352,9 +349,12 @@ static void run() // ---------------------------------------------------------------------- // hand off to PscIntegrator to run the simulation - auto psc = - makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, balance, - collision, checks, marder, diagnostics); + auto psc = makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, + balance, collision, checks, marder); + + psc.add_diagnostic(&outf); + psc.add_diagnostic(&outp); + psc.add_diagnostic(&oute); psc.integrate(); } diff --git a/src/psc_config.hxx b/src/psc_config.hxx index 9f1e40ac2e..5e5ed4f882 100644 --- a/src/psc_config.hxx +++ b/src/psc_config.hxx @@ -31,7 +31,9 @@ #include "../libpsc/vpic/fields_item_vpic.hxx" -using OutputParticlesDefault = OutputParticlesHdf5; +template +using OutputParticlesDefault = + OutputParticlesHdf5; struct SimulationNone { @@ -94,7 +96,7 @@ struct PscConfig_ using Checks = Checks_; using Marder = Marder_; using Simulation = _Simulation; - using OutputParticles = OutputParticlesDefault; + using OutputParticles = OutputParticlesDefault; }; #ifdef USE_CUDA @@ -119,7 +121,7 @@ struct PscConfig_<_Dim, _Mparticles, _MfieldsState, _Mfields, using Balance = Balance_; using Checks = ChecksCuda; using Marder = MarderCuda; - using OutputParticles = OutputParticlesDefault; + using OutputParticles = OutputParticlesDefault; }; template @@ -143,7 +145,7 @@ struct PscConfig_; using Checks = ChecksCuda; using Marder = MarderCuda; - using OutputParticles = OutputParticlesDefault; + using OutputParticles = OutputParticlesDefault; }; #endif diff --git a/src/psc_flatfoil_yz.cxx b/src/psc_flatfoil_yz.cxx index f2de74e21a..34eb008d45 100644 --- a/src/psc_flatfoil_yz.cxx +++ b/src/psc_flatfoil_yz.cxx @@ -3,7 +3,6 @@ #include #include -#include "DiagnosticsDefault.h" #include "OutputFieldsDefault.h" #include "psc_config.hxx" @@ -227,7 +226,8 @@ using Collision = PscConfig::Collision; using Checks = PscConfig::Checks; using Marder = PscConfig::Marder; #if CASE == CASE_2D_SMALL -using OutputParticles = OutputParticlesHdf5>; +using OutputParticles = + OutputParticlesHdf5>; #else using OutputParticles = PscConfig::OutputParticles; #endif @@ -543,9 +543,7 @@ void run() OutputParticles outp{grid, outp_params}; int oute_interval = -100; - DiagEnergies oute{grid.comm(), oute_interval}; - - auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); + DiagEnergies oute{grid.comm(), oute_interval}; // ---------------------------------------------------------------------- // Set up objects specific to the flatfoil case @@ -680,8 +678,14 @@ void run() // hand off to PscIntegrator to run the simulation auto psc = makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, - balance, collision, checks, marder, - diagnostics, lf_inject_heat); + balance, collision, checks, marder); + + psc.add_diagnostic(&outf); + psc.add_diagnostic(&outp); + psc.add_diagnostic(&oute); + + psc.add_injector( + new InjectFromLambda(lf_inject_heat)); MEM_STATS(); psc.integrate(); diff --git a/src/psc_harris_yz.cxx b/src/psc_harris_yz.cxx index aa2df93a27..76ddd408d4 100644 --- a/src/psc_harris_yz.cxx +++ b/src/psc_harris_yz.cxx @@ -4,7 +4,6 @@ #include #include -#include "DiagnosticsDefault.h" #include "OutputFieldsDefault.h" #include "psc_config.hxx" @@ -420,9 +419,7 @@ void run() OutputParticles outp{grid, outp_params}; int oute_interval = -100; - DiagEnergies oute{grid.comm(), oute_interval}; - - auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); + DiagEnergies oute{grid.comm(), oute_interval}; // ---------------------------------------------------------------------- // Set up objects specific to the Harris case @@ -442,9 +439,12 @@ void run() // ---------------------------------------------------------------------- // hand off to PscIntegrator to run the simulation - auto psc = - makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, balance, - collision, checks, marder, diagnostics); + auto psc = makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, + balance, collision, checks, marder); + + psc.add_diagnostic(&outf); + psc.add_diagnostic(&outp); + psc.add_diagnostic(&oute); MEM_STATS(); psc.integrate(); diff --git a/src/psc_radiation.cxx b/src/psc_radiation.cxx index 1f94a27f2d..c717a7b8eb 100644 --- a/src/psc_radiation.cxx +++ b/src/psc_radiation.cxx @@ -3,7 +3,6 @@ #include #include -#include "DiagnosticsDefault.h" #include "OutputFieldsDefault.h" #include "psc_config.hxx" @@ -234,7 +233,8 @@ void run() psc_params.marder_interval = -1; Marder marder(grid, marder_diffusion, marder_loop, marder_dump); - auto lf_ext_current = [&](const Grid_t& grid, MfieldsState& mflds) { + auto lf_ext_current = [&](MfieldsState& mflds) { + const Grid_t& grid = mflds.grid(); double time = grid.time(); auto& gdims = grid.domain.gdims; for (int p = 0; p < mflds.n_patches(); ++p) { @@ -289,24 +289,28 @@ void run() OutputParticles outp{grid, outp_params}; int oute_interval = -100; - DiagEnergies oute{grid.comm(), oute_interval}; - - auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); + DiagEnergies oute{grid.comm(), oute_interval}; // ---------------------------------------------------------------------- // setup initial conditions if (read_checkpoint_filename.empty()) { initializeFields(mflds); - lf_ext_current(*grid_ptr, mflds); + lf_ext_current(mflds); } // ---------------------------------------------------------------------- // hand off to PscIntegrator to run the simulation - auto psc = makePscIntegrator( - psc_params, *grid_ptr, mflds, mprts, balance, collision, checks, marder, - diagnostics, injectParticlesNone, lf_ext_current); + auto psc = makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, + balance, collision, checks, marder); + + psc.add_diagnostic(&outf); + psc.add_diagnostic(&outp); + psc.add_diagnostic(&oute); + + psc.add_external_current( + new ExternalCurrentFromLambda(lf_ext_current)); MEM_STATS(); psc.integrate(); diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index f89761745f..8017ed6048 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -2,7 +2,6 @@ #include #include -#include "DiagnosticsDefault.h" #include "OutputFieldsDefault.h" #include "psc_config.hxx" #include "input_params.hxx" @@ -302,9 +301,7 @@ static void run(int argc, char** argv) OutputParticles outp{grid, outp_params}; int oute_interval = -100; - DiagEnergies oute{grid.comm(), oute_interval}; - - auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); + DiagEnergies oute{grid.comm(), oute_interval}; // ---------------------------------------------------------------------- // set up initial conditions @@ -315,9 +312,12 @@ static void run(int argc, char** argv) // ---------------------------------------------------------------------- // run the simulation - auto psc = - makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, balance, - collision, checks, marder, diagnostics); + auto psc = makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, + balance, collision, checks, marder); + + psc.add_diagnostic(&outf); + psc.add_diagnostic(&outp); + psc.add_diagnostic(&oute); psc.integrate(); } diff --git a/src/psc_whistler.cxx b/src/psc_whistler.cxx index 6437bdac15..ae54879732 100644 --- a/src/psc_whistler.cxx +++ b/src/psc_whistler.cxx @@ -3,7 +3,6 @@ #include #include -#include "DiagnosticsDefault.h" #include "OutputFieldsDefault.h" #include "psc_config.hxx" @@ -314,9 +313,7 @@ void run() OutputParticles outp{grid, outp_params}; int oute_interval = 100; - DiagEnergies oute{grid.comm(), oute_interval}; - - auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); + DiagEnergies oute{grid.comm(), oute_interval}; // ---------------------------------------------------------------------- // setup initial conditions @@ -331,9 +328,12 @@ void run() // ---------------------------------------------------------------------- // hand off to PscIntegrator to run the simulation - auto psc = - makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, balance, - collision, checks, marder, diagnostics); + auto psc = makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, + balance, collision, checks, marder); + + psc.add_diagnostic(&outf); + psc.add_diagnostic(&outp); + psc.add_diagnostic(&oute); psc.integrate(); }