From 3b9d47a5589785b875ca3d8330e92de6780abe80 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Thu, 2 Sep 2021 03:22:34 +0200 Subject: [PATCH 01/15] implement Parser from amrex/WarpX --- src/Hipace.H | 1 + src/Hipace.cpp | 32 +++++++ src/utils/Parser.H | 211 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 244 insertions(+) create mode 100755 src/utils/Parser.H diff --git a/src/Hipace.H b/src/Hipace.H index f90b0681be..422ff74ce3 100644 --- a/src/Hipace.H +++ b/src/Hipace.H @@ -9,6 +9,7 @@ #include "utils/AdaptiveTimeStep.H" #include "utils/GridCurrent.H" #include "utils/Constants.H" +#include "utils/Parser.H" #include "diagnostics/Diagnostic.H" #include "diagnostics/OpenPMDWriter.H" diff --git a/src/Hipace.cpp b/src/Hipace.cpp index b391b1daeb..6ec43e4948 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -83,6 +83,38 @@ Hipace::Hipace () : } else { m_phys_const = make_constants_SI(); } + + double test_d = 0; + bool test_d_b = queryWithParser(pph, "test_d", test_d); + amrex::Print() << "##### Parser Test double " << test_d_b << ": " << test_d << '\n'; + + float test_f = 0; + bool test_f_b = queryWithParser(pph, "test_f", test_f); + amrex::Print() << "##### Parser Test float " << test_f_b << ": " << test_f << '\n'; + + int test_i = 0; + bool test_i_b = queryWithParser(pph, "test_i", test_i); + amrex::Print() << "##### Parser Test int " << test_i_b << ": " << test_i << '\n'; + + std::vector test_v{}; + bool test_v_b = queryWithParser(pph, "test_v", test_v); + amrex::Print() << "##### Parser Test vector " << test_v_b << ": "; + for (auto t : test_v) { + amrex::Print() << t << " "; + } + amrex::Print() << '\n'; + + std::array test_a{}; + bool test_a_b = queryWithParser(pph, "test_a", test_a); + amrex::Print() << "##### Parser Test array " << test_a_b << ": "; + for (auto t : test_a) { + amrex::Print() << t << " "; + } + amrex::Print() << '\n'; + + auto test_func = getFunctionWithParser<3>(pph, "test_func", {"x","y","z"}); + amrex::Print() << "##### Parser Test function " << test_func(1.,2.,3.) << '\n'; + pph.query("dt", m_dt); pph.query("verbose", m_verbose); pph.query("numprocs_x", m_numprocs_x); diff --git a/src/utils/Parser.H b/src/utils/Parser.H new file mode 100755 index 0000000000..bdf40984d4 --- /dev/null +++ b/src/utils/Parser.H @@ -0,0 +1,211 @@ +#ifndef HIPACE_Parser_H_ +#define HIPACE_Parser_H_ + +#include "Constants.H" + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + + +template +inline void +getWithParser (const amrex::ParmParse& pp, char const * const str, T& val); + +template +inline bool +queryWithParser (const amrex::ParmParse& pp, char const * const str, T& val); + +inline int +safeCastToInt(const double x, const std::string& real_name) { + int result = 0; + bool error_detected = false; + std::string assert_msg; + // (2.0*(numeric_limits::max()/2+1)) converts numeric_limits::max()+1 to a real ensuring accuracy to all digits + // This accepts x = 2**31-1 but rejects 2**31. + if (x < (2.0*(std::numeric_limits::max()/2+1))) { + if (std::ceil(x) >= std::numeric_limits::min()) { + result = static_cast(x); + } else { + error_detected = true; + assert_msg = "Error: Negative overflow detected when casting " + real_name + " = " + std::to_string(x) + " to int"; + } + } else if (x > 0) { + error_detected = true; + assert_msg = "Error: Overflow detected when casting " + real_name + " = " + std::to_string(x) + " to int"; + } else { + error_detected = true; + assert_msg = "Error: NaN detected when casting " + real_name + " to int"; + } + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!error_detected, assert_msg); + return result; +} + +inline amrex::Parser +makeParser (std::string const& parse_function, amrex::Vector const& varnames) +{ + // Since queryWithParser recursively calls this routine, keep track of symbols + // in case an infinite recursion is found (a symbol's value depending on itself). + static std::set recursive_symbols{}; + + amrex::Parser parser(parse_function); + parser.registerVariables(varnames); + + std::set symbols = parser.symbols(); + for (auto const& v : varnames) symbols.erase(v); + + // User can provide inputs under this name, through which expressions + // can be provided for arbitrary variables. PICMI inputs are aware of + // this convention and use the same prefix as well. This potentially + // includes variable names that match physical or mathematical + // constants, in case the user wishes to enforce a different + // system of units or some form of quasi-physical behavior in the + // simulation. Thus, this needs to override any built-in + // constants. + amrex::ParmParse pp_my_constants("my_constants"); + + // Cache for evalutated constants + static std::map my_constants_cache{}; + + // Physical / Numerical Constants available to parsed expressions + const auto phys_const = make_constants_SI(); + + static std::map hipace_constants + { + {"clight", phys_const.c}, + {"epsilon0", phys_const.ep0}, + {"mu0", phys_const.mu0}, + {"q_e", phys_const.q_e}, + {"m_e", phys_const.m_e}, + {"m_p", phys_const.m_p}, + {"pi", MathConst::pi} + }; + + for (auto const& s : symbols) { + + auto usr_constant = my_constants_cache.find(s); + if (usr_constant != my_constants_cache.end()) { + parser.setConstant(s, usr_constant->second); + continue; + } + + //AmrexAlwaysAssert(recursive_symbols.count(*it)==0, "Expressions contains recursive symbol "+*it); + double v; + recursive_symbols.insert(s); + const bool is_input = queryWithParser(pp_my_constants, s.c_str(), v); + recursive_symbols.erase(s); + + if (is_input) { + my_constants_cache.insert({s, v}); + parser.setConstant(s, v); + continue; + } + + auto phy_constant = hipace_constants.find(s); + if (phy_constant != hipace_constants.end()) { + parser.setConstant(s, phy_constant->second); + continue; + } + + amrex::Abort("makeParser::Unknown symbol " + s); + } + parser.print(); + return parser; +} + +template +inline auto +getFunctionWithParser (const amrex::ParmParse& pp, + char const * const str, + amrex::Vector const& varnames) +{ + std::string func_str{}; + std::vector f; + pp.getarr(str, f); + for (auto const& s : f) { + func_str += s; + } + return makeParser(func_str, varnames).compile(); +} + +inline void +fillWithParser (std::string const& str, double& val) { + val = makeParser(str, {}).compileHost<0>()(); +} + +inline void +fillWithParser (std::string const& str, float& val) { + val = static_cast(makeParser(str, {}).compileHost<0>()()); +} + +inline void +fillWithParser (std::string const& str, int& val) { + val = safeCastToInt(makeParser(str, {}).compileHost<0>()(),str); +} + +template +inline void +fillWithParserArr (std::vector const& str_arr, T& val) { + std::string str{}; + for (auto const& s : str_arr) { + str += s; + } + fillWithParser(str, val); +} + +template +inline void +fillWithParserArr (std::vector const& str_arr, std::vector& val_arr) { + auto const n = str_arr.size(); + val_arr.resize(n); + for (auto i=0ul ; i != n ; ++i) { + fillWithParser(str_arr[i], val_arr[i]); + } +} + +template +inline void +fillWithParserArr (std::vector const& str_arr, std::array& val_arr) { + const auto n = str_arr.size(); + if (n != size) { + for( auto const& str : str_arr) { + amrex::ErrorStream() << str << ' '; + } + amrex::ErrorStream() << "has wrong length " << n << " should be " << size << '\n'; + } + AMREX_ALWAYS_ASSERT( n == size ); + for (auto i=0ul ; i != n ; ++i) { + fillWithParser(str_arr[i], val_arr[i]); + } +} + +template +inline void +getWithParser (const amrex::ParmParse& pp, char const * const str, T& val) +{ + std::vector f; + pp.getarr(str, f); + fillWithParserArr(f, val); +} + +template +inline bool +queryWithParser (const amrex::ParmParse& pp, char const * const str, T& val) +{ + std::vector f; + bool is_specified = pp.queryarr(str, f); + if (is_specified) + { + fillWithParserArr(f, val); + } + return is_specified; +} + +#endif From a8f89f2460ffbc71141249794aad9d4edf3b9224 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Wed, 8 Sep 2021 16:07:18 +0200 Subject: [PATCH 02/15] initialize PhysConst earlier --- src/Hipace.H | 19 +++++++++--------- src/Hipace.cpp | 24 +++++++++++------------ src/particles/BeamParticleContainer.cpp | 8 +------- src/particles/PlasmaParticleContainer.cpp | 10 ++-------- 4 files changed, 24 insertions(+), 37 deletions(-) diff --git a/src/Hipace.H b/src/Hipace.H index 6e62cb575c..6ffbab9e02 100644 --- a/src/Hipace.H +++ b/src/Hipace.H @@ -190,6 +190,15 @@ public: /** status of the physical time send request */ MPI_Request m_tsend_request = MPI_REQUEST_NULL; +private: + /** Pointer to current (and only) instance of class Hipace */ + static Hipace* m_instance; +public: + /** Whether to use normalized units */ + static bool m_normalized_units; + /** Struct containing physical constants (which values depends on the unit system, determined + * at runtime): SI or normalized units. */ + PhysConst m_phys_const; /** All field data (3D array, slices) and field methods */ Fields m_fields; /** Contains all beam species */ @@ -208,11 +217,6 @@ public: static amrex::Real m_physical_time; /** Level of verbosity */ static int m_verbose; - /** Whether to use normalized units */ - static bool m_normalized_units; - /** Struct containing physical constants (which values depends on the unit system, determined - * at runtime): SI or normalized units. */ - PhysConst m_phys_const; /** Order of the field gather and current deposition shape factor in the transverse directions */ static int m_depos_order_xy; @@ -292,9 +296,6 @@ public: amrex::RealVect patch_lo {0., 0., 0.}; /**< 3D array with lower ends of the refined grid */ amrex::RealVect patch_hi {0., 0., 0.}; /**< 3D array with upper ends of the refined grid */ private: - /** Pointer to current (and only) instance of class Hipace */ - static Hipace* m_instance; - amrex::Vector m_slice_geom; amrex::Vector m_slice_dm; amrex::Vector m_slice_ba; @@ -348,8 +349,6 @@ private: void PredictorCorrectorLoopToSolveBxBy (const int islice_local, const int lev, const amrex::Box bx, amrex::Vector bins, const int ibox); - void Ionisation (const int lev); - /** \brief define Geometry, DistributionMapping and BoxArray for the slice. * The slice MultiFAB and the plasma particles are defined on this GDB. * diff --git a/src/Hipace.cpp b/src/Hipace.cpp index bd3421ecd8..39dc9b349d 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -27,9 +27,9 @@ namespace { Hipace* Hipace::m_instance = nullptr; +bool Hipace::m_normalized_units = false; int Hipace::m_max_step = 0; amrex::Real Hipace::m_dt = 0.0; -bool Hipace::m_normalized_units = false; amrex::Real Hipace::m_physical_time = 0.0; int Hipace::m_verbose = 0; int Hipace::m_depos_order_xy = 2; @@ -54,20 +54,26 @@ bool Hipace::m_do_tiling = true; Hipace& Hipace::GetInstance () { - if (!m_instance) { - m_instance = new Hipace(); - } + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_instance, "instance has not been initialized yet"); return *m_instance; } Hipace::Hipace () : + m_phys_const([this](){ + m_instance = this; + amrex::ParmParse pph("hipace"); + pph.query("normalized_units", this->m_normalized_units); + if (this->m_normalized_units){ + return make_constants_normalized(); + } else { + return make_constants_SI(); + } + }()), m_fields(this), m_multi_beam(this), m_multi_plasma(this), m_diags(this->maxLevel()+1) { - m_instance = this; - amrex::ParmParse pp;// Traditionally, max_step and stop_time do not have prefix. pp.query("max_step", m_max_step); @@ -77,12 +83,6 @@ Hipace::Hipace () : #endif amrex::ParmParse pph("hipace"); - pph.query("normalized_units", m_normalized_units); - if (m_normalized_units){ - m_phys_const = make_constants_normalized(); - } else { - m_phys_const = make_constants_SI(); - } pph.query("dt", m_dt); pph.query("verbose", m_verbose); pph.query("numprocs_x", m_numprocs_x); diff --git a/src/particles/BeamParticleContainer.cpp b/src/particles/BeamParticleContainer.cpp index 56bb7294d2..707589c64a 100644 --- a/src/particles/BeamParticleContainer.cpp +++ b/src/particles/BeamParticleContainer.cpp @@ -7,13 +7,7 @@ namespace { void QueryElementSetChargeMass (amrex::ParmParse& pp, amrex::Real& charge, amrex::Real& mass) { - // normalized_units is directly queried here so we can defined the appropriate PhysConst - // locally. We cannot use Hipace::m_phys_const as it has not been initialized when the - // PlasmaParticleContainer constructor is called. - amrex::ParmParse pph("hipace"); - bool normalized_units = false; - pph.query("normalized_units", normalized_units); - PhysConst phys_const = normalized_units ? make_constants_normalized() : make_constants_SI(); + PhysConst phys_const = get_phys_const(); std::string element = "electron"; pp.query("element", element); diff --git a/src/particles/PlasmaParticleContainer.cpp b/src/particles/PlasmaParticleContainer.cpp index ef2ce19a0d..909faa82e1 100644 --- a/src/particles/PlasmaParticleContainer.cpp +++ b/src/particles/PlasmaParticleContainer.cpp @@ -11,13 +11,7 @@ void PlasmaParticleContainer::ReadParameters () { - // normalized_units is directly queried here so we can defined the appropriate PhysConst - // locally. We cannot use Hipace::m_phys_const as it has not been initialized when the - // PlasmaParticleContainer constructor is called. - amrex::ParmParse pph("hipace"); - bool normalized_units = false; - pph.query("normalized_units", normalized_units); - PhysConst phys_const = normalized_units ? make_constants_normalized() : make_constants_SI(); + PhysConst phys_const = get_phys_const(); amrex::ParmParse pp(m_name); std::string element = ""; @@ -51,7 +45,7 @@ PlasmaParticleContainer::ReadParameters () pp.query("can_ionize", m_can_ionize); if(m_can_ionize) { m_neutralize_background = false; // change default - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!normalized_units, + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!Hipace::GetInstance().m_normalized_units, "Cannot use Ionization Module in normalized units"); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_init_ion_lev >= 0, "The initial Ion level must be specified"); From 624fb0e5333c9d3a6de9dd352048168f2c6af190 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Wed, 8 Sep 2021 18:31:55 +0200 Subject: [PATCH 03/15] readd assert --- src/utils/Parser.H | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/src/utils/Parser.H b/src/utils/Parser.H index bdf40984d4..32776b1720 100755 --- a/src/utils/Parser.H +++ b/src/utils/Parser.H @@ -9,6 +9,7 @@ #include #include +#include #include #include #include @@ -28,18 +29,20 @@ safeCastToInt(const double x, const std::string& real_name) { int result = 0; bool error_detected = false; std::string assert_msg; - // (2.0*(numeric_limits::max()/2+1)) converts numeric_limits::max()+1 to a real ensuring accuracy to all digits - // This accepts x = 2**31-1 but rejects 2**31. + // (2.0*(numeric_limits::max()/2+1)) converts numeric_limits::max()+1 to a real + // ensuring accuracy to all digits. This accepts x = 2**31-1 but rejects 2**31. if (x < (2.0*(std::numeric_limits::max()/2+1))) { if (std::ceil(x) >= std::numeric_limits::min()) { result = static_cast(x); } else { error_detected = true; - assert_msg = "Error: Negative overflow detected when casting " + real_name + " = " + std::to_string(x) + " to int"; + assert_msg = "Error: Negative overflow detected when casting " + + real_name + " = " + std::to_string(x) + " to int"; } } else if (x > 0) { error_detected = true; - assert_msg = "Error: Overflow detected when casting " + real_name + " = " + std::to_string(x) + " to int"; + assert_msg = "Error: Overflow detected when casting " + + real_name + " = " + std::to_string(x) + " to int"; } else { error_detected = true; assert_msg = "Error: NaN detected when casting " + real_name + " to int"; @@ -48,7 +51,7 @@ safeCastToInt(const double x, const std::string& real_name) { return result; } -inline amrex::Parser +static amrex::Parser makeParser (std::string const& parse_function, amrex::Vector const& varnames) { // Since queryWithParser recursively calls this routine, keep track of symbols @@ -62,8 +65,7 @@ makeParser (std::string const& parse_function, amrex::Vector const& for (auto const& v : varnames) symbols.erase(v); // User can provide inputs under this name, through which expressions - // can be provided for arbitrary variables. PICMI inputs are aware of - // this convention and use the same prefix as well. This potentially + // can be provided for arbitrary variables. This potentially // includes variable names that match physical or mathematical // constants, in case the user wishes to enforce a different // system of units or some form of quasi-physical behavior in the @@ -75,7 +77,7 @@ makeParser (std::string const& parse_function, amrex::Vector const& static std::map my_constants_cache{}; // Physical / Numerical Constants available to parsed expressions - const auto phys_const = make_constants_SI(); + const auto phys_const = get_phys_const(); static std::map hipace_constants { @@ -96,7 +98,8 @@ makeParser (std::string const& parse_function, amrex::Vector const& continue; } - //AmrexAlwaysAssert(recursive_symbols.count(*it)==0, "Expressions contains recursive symbol "+*it); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(recursive_symbols.count(*it)==0, + "Expressions contains recursive symbol " + *it); double v; recursive_symbols.insert(s); const bool is_input = queryWithParser(pp_my_constants, s.c_str(), v); @@ -123,8 +126,8 @@ makeParser (std::string const& parse_function, amrex::Vector const& template inline auto getFunctionWithParser (const amrex::ParmParse& pp, - char const * const str, - amrex::Vector const& varnames) + char const * const str, + amrex::Vector const& varnames) { std::string func_str{}; std::vector f; @@ -147,7 +150,7 @@ fillWithParser (std::string const& str, float& val) { inline void fillWithParser (std::string const& str, int& val) { - val = safeCastToInt(makeParser(str, {}).compileHost<0>()(),str); + val = safeCastToInt(std::round(makeParser(str, {}).compileHost<0>()()),str); } template From 76b10a30a3ea2000ef810c5f4424beeb8c7ab09d Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Wed, 8 Sep 2021 18:57:07 +0200 Subject: [PATCH 04/15] start to add Parser --- src/Hipace.cpp | 48 +++++----- src/fields/Fields.cpp | 201 ++++++++++++++++-------------------------- 2 files changed, 100 insertions(+), 149 deletions(-) diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 3755051fa5..a39ec89e17 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -75,11 +75,11 @@ Hipace::Hipace () : m_diags(this->maxLevel()+1) { amrex::ParmParse pp;// Traditionally, max_step and stop_time do not have prefix. - pp.query("max_step", m_max_step); + queryWithParser(pp, "max_step", m_max_step); #ifndef AMREX_USE_GPU int seed; - if (pp.query("random_seed", seed)) amrex::ResetRandomSeed(seed); + if (queryWithParser(pp, "random_seed", seed)) amrex::ResetRandomSeed(seed); #endif amrex::ParmParse pph("hipace"); @@ -115,35 +115,35 @@ Hipace::Hipace () : auto test_func = getFunctionWithParser<3>(pph, "test_func", {"x","y","z"}); amrex::Print() << "##### Parser Test function " << test_func(1.,2.,3.) << '\n'; - pph.query("dt", m_dt); - pph.query("verbose", m_verbose); - pph.query("numprocs_x", m_numprocs_x); - pph.query("numprocs_y", m_numprocs_y); + queryWithParser(pph, "dt", m_dt); + queryWithParser(pph, "verbose", m_verbose); + queryWithParser(pph, "numprocs_x", m_numprocs_x); + queryWithParser(pph, "numprocs_y", m_numprocs_y); m_numprocs_z = amrex::ParallelDescriptor::NProcs() / (m_numprocs_x*m_numprocs_y); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_numprocs_z <= m_max_step+1, "Please use more or equal time steps than number of ranks"); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_numprocs_x*m_numprocs_y*m_numprocs_z == amrex::ParallelDescriptor::NProcs(), "Check hipace.numprocs_x and hipace.numprocs_y"); - pph.query("boxes_in_z", m_boxes_in_z); + queryWithParser(pph, "boxes_in_z", m_boxes_in_z); if (m_boxes_in_z > 1) AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_numprocs_z == 1, "Multiple boxes per rank only implemented for one rank."); - pph.query("depos_order_xy", m_depos_order_xy); - pph.query("depos_order_z", m_depos_order_z); - pph.query("predcorr_B_error_tolerance", m_predcorr_B_error_tolerance); - pph.query("predcorr_max_iterations", m_predcorr_max_iterations); - pph.query("predcorr_B_mixing_factor", m_predcorr_B_mixing_factor); - pph.query("output_period", m_output_period); + queryWithParser(pph, "depos_order_xy", m_depos_order_xy); + queryWithParser(pph, "depos_order_z", m_depos_order_z); + queryWithParser(pph, "predcorr_B_error_tolerance", m_predcorr_B_error_tolerance); + queryWithParser(pph, "predcorr_max_iterations", m_predcorr_max_iterations); + queryWithParser(pph, "predcorr_B_mixing_factor", m_predcorr_B_mixing_factor); + queryWithParser(pph, "output_period", m_output_period); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_output_period != 0, "To avoid output, please use output_period = -1."); - pph.query("beam_injection_cr", m_beam_injection_cr); - pph.query("do_beam_jx_jy_deposition", m_do_beam_jx_jy_deposition); - pph.query("do_device_synchronize", m_do_device_synchronize); - pph.query("external_ExmBy_slope", m_external_ExmBy_slope); - pph.query("external_Ez_slope", m_external_Ez_slope); - pph.query("external_Ez_uniform", m_external_Ez_uniform); + queryWithParser(pph, "beam_injection_cr", m_beam_injection_cr); + queryWithParser(pph, "do_beam_jx_jy_deposition", m_do_beam_jx_jy_deposition); + queryWithParser(pph, "do_device_synchronize", m_do_device_synchronize); + queryWithParser(pph, "external_ExmBy_slope", m_external_ExmBy_slope); + queryWithParser(pph, "external_Ez_slope", m_external_Ez_slope); + queryWithParser(pph, "external_Ez_uniform", m_external_Ez_uniform); std::string solver = "predictor-corrector"; - pph.query("bxby_solver", solver); + queryWithParser(pph, "bxby_solver", solver); AMREX_ALWAYS_ASSERT_WITH_MESSAGE( solver == "predictor-corrector" || solver == "explicit", @@ -153,9 +153,9 @@ Hipace::Hipace () : !(m_explicit && !m_multi_plasma.AllSpeciesNeutralizeBackground()), "Ion motion with explicit solver is not implemented, need to use neutralize_background"); - pph.query("MG_tolerance_rel", m_MG_tolerance_rel); - pph.query("MG_tolerance_abs", m_MG_tolerance_abs); - pph.query("do_tiling", m_do_tiling); + queryWithParser(pph, "MG_tolerance_rel", m_MG_tolerance_rel); + queryWithParser(pph, "MG_tolerance_abs", m_MG_tolerance_abs); + queryWithParser(pph, "do_tiling", m_do_tiling); #ifdef AMREX_USE_GPU AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_do_tiling==0, "Tiling must be turned off to run on GPU."); #endif @@ -172,7 +172,7 @@ Hipace::Hipace () : } #ifdef AMREX_USE_MPI - pph.query("skip_empty_comms", m_skip_empty_comms); + queryWithParser(pph, "skip_empty_comms", m_skip_empty_comms); int myproc = amrex::ParallelDescriptor::MyProc(); m_rank_z = myproc/(m_numprocs_x*m_numprocs_y); MPI_Comm_split(amrex::ParallelDescriptor::Communicator(), m_rank_z, myproc, &m_comm_xy); diff --git a/src/fields/Fields.cpp b/src/fields/Fields.cpp index c4f99dda2f..0a59b61298 100644 --- a/src/fields/Fields.cpp +++ b/src/fields/Fields.cpp @@ -12,7 +12,7 @@ Fields::Fields (Hipace const* a_hipace) : m_slices(a_hipace->maxLevel()+1) { amrex::ParmParse ppf("fields"); - ppf.query("do_dirichlet_poisson", m_do_dirichlet_poisson); + queryWithParser(ppf, "do_dirichlet_poisson", m_do_dirichlet_poisson); } void @@ -233,14 +233,12 @@ Fields::LongitudinalDerivative (const amrex::MultiFab& src1, const amrex::MultiF void -Fields::Copy (const int lev, const int i_slice, const int slice_comp, const int full_comp, - const amrex::Gpu::DeviceVector& diag_comps_vect, - const int ncomp, amrex::FArrayBox& fab, const int slice_dir, - const amrex::IntVect diag_coarsen, const amrex::Geometry geom) +Fields::Copy (int lev, int i_slice, FieldCopyType copy_type, int slice_comp, int full_comp, + int ncomp, amrex::FArrayBox& fab, int slice_dir, amrex::Geometry geom) { using namespace amrex::literals; HIPACE_PROFILE("Fields::Copy()"); - auto& slice_mf = m_slices[lev][WhichSlice::This]; // copy from the current slice + auto& slice_mf = m_slices[lev][WhichSlice::This]; // copy from/to the current slice amrex::Array4 slice_array; // There is only one Box. for (amrex::MFIter mfi(slice_mf); mfi.isValid(); ++mfi) { auto& slice_fab = slice_mf[mfi]; @@ -251,92 +249,64 @@ Fields::Copy (const int lev, const int i_slice, const int slice_comp, const int // slice_array's longitude index is i_slice. } - const int full_array_z = i_slice / diag_coarsen[2]; - const amrex::IntVect ncells_global = geom.Domain().length(); - amrex::Box const& vbx = fab.box(); - if (vbx.smallEnd(Direction::z) <= full_array_z and - vbx.bigEnd (Direction::z) >= full_array_z and - ( i_slice % diag_coarsen[2] == diag_coarsen[2]/2 or - ( i_slice == ncells_global[2] - 1 and - ( ncells_global[2] - 1 ) % diag_coarsen[2] < diag_coarsen[2]/2 ))) + if (vbx.smallEnd(Direction::z) <= i_slice and + vbx.bigEnd (Direction::z) >= i_slice) { amrex::Box copy_box = vbx; - copy_box.setSmall(Direction::z, full_array_z); - copy_box.setBig (Direction::z, full_array_z); + copy_box.setSmall(Direction::z, i_slice); + copy_box.setBig (Direction::z, i_slice); amrex::Array4 const& full_array = fab.array(); - const int even_slice_x = ncells_global[0] % 2 == 0 and slice_dir == 0; - const int even_slice_y = ncells_global[1] % 2 == 0 and slice_dir == 1; - - const int coarse_x = diag_coarsen[0]; - const int coarse_y = diag_coarsen[1]; - - const int ncells_x = ncells_global[0]; - const int ncells_y = ncells_global[1]; - - const int *diag_comps = diag_comps_vect.data(); - - amrex::ParallelFor(copy_box, ncomp, - [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - const int m = n[diag_comps]; - - // coarsening in slice direction is always 1 - const int i_c_start = amrex::min(i*coarse_x +(coarse_x-1)/2 -even_slice_x, ncells_x-1); - const int i_c_stop = amrex::min(i*coarse_x +coarse_x/2+1, ncells_x); - const int j_c_start = amrex::min(j*coarse_y +(coarse_y-1)/2 -even_slice_y, ncells_y-1); - const int j_c_stop = amrex::min(j*coarse_y +coarse_y/2+1, ncells_y); + const amrex::IntVect ncells_global = geom.Domain().length(); + const bool nx_even = ncells_global[0] % 2 == 0; + const bool ny_even = ncells_global[1] % 2 == 0; - amrex::Real field_value = 0._rt; - int n_values = 0; - - for (int j_c = j_c_start; j_c != j_c_stop; ++j_c) { - for (int i_c = i_c_start; i_c != i_c_stop; ++i_c) { - field_value += slice_array(i_c, j_c, i_slice, m+slice_comp); - ++n_values; + if (copy_type == FieldCopyType::FtoS) { + amrex::ParallelFor(copy_box, ncomp, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + slice_array(i,j,k,n+slice_comp) = full_array(i,j,k,n+full_comp); + }); + } else { + amrex::ParallelFor(copy_box, ncomp, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + if (slice_dir ==-1 /* 3D data */){ + full_array(i,j,k,n+full_comp) = slice_array(i,j,k,n+slice_comp); + } else if (slice_dir == 0 /* yz slice */){ + full_array(i,j,k,n+full_comp) = + nx_even ? 0.5_rt * (slice_array(i-1,j,k,n+slice_comp) + + slice_array(i,j,k,n+slice_comp)) + : slice_array(i,j,k,n+slice_comp); + } else /* slice_dir == 1, xz slice */{ + full_array(i,j,k,n+full_comp) = + ny_even ? 0.5_rt * ( slice_array(i,j-1,k,n+slice_comp) + + slice_array(i,j,k,n+slice_comp)) + : slice_array(i,j,k,n+slice_comp); } - } - - full_array(i,j,k,n+full_comp) = field_value / amrex::max(n_values,1); - }); + }); + } } } void -Fields::ShiftSlices (int nlev, int islice, amrex::Geometry geom, amrex::Real patch_lo, - amrex::Real patch_hi) +Fields::ShiftSlices (int lev) { HIPACE_PROFILE("Fields::ShiftSlices()"); - - for (int lev = 0; lev < nlev; ++lev) { - - // skip all slices which are not existing on level 1 - if (lev == 1) { - // use geometry of coarse grid to determine whether slice is to be solved - const amrex::Real* problo = geom.ProbLo(); - const amrex::Real* dx = geom.CellSize(); - const amrex::Real pos = (islice+0.5)*dx[2]+problo[2]; - if (pos < patch_lo || pos > patch_hi) continue; - } - - // shift Bx, By amrex::MultiFab::Copy( getSlices(lev, WhichSlice::Previous2), getSlices(lev, WhichSlice::Previous1), Comps[WhichSlice::Previous1]["Bx"], Comps[WhichSlice::Previous2]["Bx"], 2, m_slices_nguards); - // shift Ez, Bx, By, Bz, jx, jx_beam, jy, jy_beam amrex::MultiFab::Copy( getSlices(lev, WhichSlice::Previous1), getSlices(lev, WhichSlice::This), - Comps[WhichSlice::This]["Ez"], Comps[WhichSlice::Previous1]["Ez"], - 8, m_slices_nguards); - // shift rho, Psi + Comps[WhichSlice::This]["Bx"], Comps[WhichSlice::Previous1]["Bx"], + 2, m_slices_nguards); amrex::MultiFab::Copy( getSlices(lev, WhichSlice::Previous1), getSlices(lev, WhichSlice::This), - Comps[WhichSlice::This]["rho"], Comps[WhichSlice::Previous1]["rho"], - 2, m_slices_nguards); - } + Comps[WhichSlice::This]["jx"], Comps[WhichSlice::Previous1]["jx"], + 4, m_slices_nguards); } void @@ -370,14 +340,13 @@ Fields::AddBeamCurrents (const int lev, const int which_slice) void Fields::InterpolateBoundaries (amrex::Vector const& geom, const int lev, - std::string component, const int islice) + std::string component) { // To solve a Poisson equation with non-zero Dirichlet boundary conditions, the source term // must be corrected at the outmost grid points in x by -field_value_at_guard_cell / dx^2 and // in y by -field_value_at_guard_cell / dy^2, where dx and dy are those of the fine grid // This follows Van Loan, C. (1992). Computational frameworks for the fast Fourier transform. // Page 254 ff. - // The interpolation is done in second order transversely and linearly in longitudinal direction HIPACE_PROFILE("Fields::InterpolateBoundaries()"); if (lev == 0) return; // only interpolate boundaries to lev 1 @@ -386,18 +355,11 @@ Fields::InterpolateBoundaries (amrex::Vector const& geom, const const auto dx = geom[lev].CellSizeArray(); const auto plo_coarse = geom[lev-1].ProbLoArray(); const auto dx_coarse = geom[lev-1].CellSizeArray(); -constexpr int interp_order = 2; - - // get relative position of fine grid slice between coarse grids for longitudinal lin. interpol. - const amrex::Real z = plo_coarse[2] + (islice+0.5_rt)*dx[2]; - const int idz_coarse = (z-plo_coarse[2])/dx_coarse[2]; - const amrex::Real rel_z = (z - (plo_coarse[2] + (idz_coarse)*dx_coarse[2])) / dx_coarse[2]; + constexpr int interp_order = 2; // get level 0 for interpolation to source term of level 1 amrex::MultiFab lhs_coarse(getSlices(lev-1, WhichSlice::This), amrex::make_alias, Comps[WhichSlice::This][component], 1); - amrex::MultiFab lhs_coarse_prev(getSlices(lev-1, WhichSlice::Previous1), amrex::make_alias, - Comps[WhichSlice::Previous1][component], 1); amrex::FArrayBox& lhs_fab = lhs_coarse[0]; amrex::Box lhs_bx = lhs_fab.box(); lhs_bx.grow({-m_slices_nguards[0], -m_slices_nguards[1], 0}); @@ -422,8 +384,7 @@ constexpr int interp_order = 2; const int nx_fine_high = big[0]; const int ny_fine_high = big[1]; amrex::Array4 data_array = m_poisson_solver[lev]->StagingArea().array(mfi); - amrex::Array4 arr_coarse = lhs_coarse.array(mfi); - amrex::Array4 arr_coarse_prev = lhs_coarse_prev.array(mfi); + amrex::Array4 data_array_coarse = lhs_coarse.array(mfi); // Loop over the valid indices on the fine grid and interpolate the value of the coarse grid // at the location of the guard cell on the fine grid to the first/last valid grid point on @@ -453,22 +414,23 @@ constexpr int interp_order = 2; // sx_cell shape factor along x const amrex::Real xmid = (x - plo_coarse[0])/dx_coarse[0]; amrex::Real sx_cell[interp_order + 1]; - const int j_cell = compute_shape_factor(sx_cell, xmid-0.5_rt); + const int j_cell = compute_shape_factor(sx_cell, + xmid-0.5_rt); // y direction const amrex::Real ymid = (y - plo_coarse[1])/dx_coarse[1]; amrex::Real sy_cell[interp_order + 1]; - const int k_cell = compute_shape_factor(sy_cell, ymid-0.5_rt); + const int k_cell = compute_shape_factor(sy_cell, + ymid-0.5_rt); amrex::Real boundary_value = 0.0_rt; // add interpolated contribution to boundary value for (int iy=0; iy<=interp_order; iy++){ for (int ix=0; ix<=interp_order; ix++){ - boundary_value += sx_cell[ix]*sy_cell[iy]* - ((1.0_rt-rel_z)*arr_coarse(lo_coarse[0]+j_cell+ix, - lo_coarse[1]+k_cell+iy, lo_coarse[2]) - + rel_z*arr_coarse_prev(lo_coarse[0]+j_cell+ix, - lo_coarse[1]+k_cell+iy, lo_coarse[2])); + boundary_value += data_array_coarse(lo_coarse[0]+j_cell+ix, + lo_coarse[1]+k_cell+iy, + lo_coarse[2]) + *sx_cell[ix]*sy_cell[iy]; } } @@ -493,22 +455,23 @@ constexpr int interp_order = 2; // sx_cell shape factor along x const amrex::Real xmid = (x - plo_coarse[0])/dx_coarse[0]; amrex::Real sx_cell[interp_order + 1]; - const int j_cell = compute_shape_factor(sx_cell, xmid-0.5_rt); + const int j_cell = compute_shape_factor(sx_cell, + xmid-0.5_rt); // y direction const amrex::Real ymid = (y - plo_coarse[1])/dx_coarse[1]; amrex::Real sy_cell[interp_order + 1]; - const int k_cell = compute_shape_factor(sy_cell, ymid-0.5_rt); + const int k_cell = compute_shape_factor(sy_cell, + ymid-0.5_rt); amrex::Real boundary_value = 0.0_rt; // add interpolated contribution to boundary value for (int iy=0; iy<=interp_order; iy++){ for (int ix=0; ix<=interp_order; ix++){ - boundary_value += sx_cell[ix]*sy_cell[iy]* - ((1.0_rt-rel_z)*arr_coarse(lo_coarse[0]+j_cell+ix, - lo_coarse[1]+k_cell+iy, lo_coarse[2]) - + rel_z*arr_coarse_prev(lo_coarse[0]+j_cell+ix, - lo_coarse[1]+k_cell+iy, lo_coarse[2])); + boundary_value += data_array_coarse(lo_coarse[0]+j_cell+ix, + lo_coarse[1]+k_cell+iy, + lo_coarse[2]) + *sx_cell[ix]*sy_cell[iy]; } } @@ -522,12 +485,11 @@ constexpr int interp_order = 2; void Fields::InterpolateFromLev0toLev1 (amrex::Vector const& geom, const int lev, - std::string component, const int islice) + std::string component) { // This function interpolates values from the coarse to the fine grid with second order. // This is required for rho to fix the incomplete deposition close to the boundary and for Psi // to fill the guard cell, which is needed for the transverse derivative - // The interpolation is done in second order transversely and linearly in longitudinal direction HIPACE_PROFILE("Fields::InterpolateFromLev0toLev1()"); if (lev == 0) return; // only interpolate boundaries to lev 1 @@ -538,16 +500,9 @@ Fields::InterpolateFromLev0toLev1 (amrex::Vector const& geom, c const auto dx_coarse = geom[lev-1].CellSizeArray(); constexpr int interp_order = 2; - // get relative position of fine grid slice between coarse grids for longitudinal lin. interpol. - const amrex::Real z = plo_coarse[2] + (islice+0.5_rt)*dx[2]; - const int idz_coarse = (z-plo_coarse[2])/dx_coarse[2]; - const amrex::Real rel_z = (z - (plo_coarse[2] + (idz_coarse)*dx_coarse[2])) / dx_coarse[2]; - // get level 0 array amrex::MultiFab lhs_coarse(getSlices(lev-1, WhichSlice::This), amrex::make_alias, Comps[WhichSlice::This][component], 1); - amrex::MultiFab lhs_coarse_prev(getSlices(lev-1, WhichSlice::Previous1), amrex::make_alias, - Comps[WhichSlice::Previous1][component], 1); amrex::FArrayBox& lhs_fab = lhs_coarse[0]; amrex::Box lhs_bx = lhs_fab.box(); // lhs_bx should only have valid cells @@ -581,8 +536,7 @@ Fields::InterpolateFromLev0toLev1 (amrex::Vector const& geom, c const int y_range = (component == "rho") ? m_slices_nguards[1] : 1; amrex::Array4 data_array = lhs_fine.array(mfi); - amrex::Array4 arr_coarse = lhs_coarse.array(mfi); - amrex::Array4 arr_coarse_prev = lhs_coarse_prev.array(mfi); + amrex::Array4 data_array_coarse = lhs_coarse.array(mfi); // Loop over the valid indices on the fine grid and interpolate the value of the coarse grid amrex::ParallelFor( @@ -614,11 +568,9 @@ Fields::InterpolateFromLev0toLev1 (amrex::Vector const& geom, c // sum interpolated contributions for (int iy=0; iy<=interp_order; iy++){ for (int ix=0; ix<=interp_order; ix++){ - coarse_value += sx_cell[ix]*sy_cell[iy]* - ((1.0_rt-rel_z)*arr_coarse(lo_coarse[0]+j_cell+ix, - lo_coarse[1]+k_cell+iy, lo_coarse[2]) - + rel_z*arr_coarse_prev(lo_coarse[0]+j_cell+ix, - lo_coarse[1]+k_cell+iy, lo_coarse[2])); + coarse_value += data_array_coarse(lo_coarse[0]+j_cell+ix, + lo_coarse[1]+k_cell+iy, + lo_coarse[2])*sx_cell[ix]*sy_cell[iy]; } } @@ -629,10 +581,9 @@ Fields::InterpolateFromLev0toLev1 (amrex::Vector const& geom, c } } - void Fields::SolvePoissonExmByAndEypBx (amrex::Vector const& geom, - const MPI_Comm& m_comm_xy, const int lev, const int islice) + const MPI_Comm& m_comm_xy, const int lev) { /* Solves Laplacian(Psi) = 1/episilon0 * -(rho-Jz/c) and * calculates Ex-c By, Ey + c Bx from grad(-Psi) @@ -645,7 +596,7 @@ Fields::SolvePoissonExmByAndEypBx (amrex::Vector const& geom, amrex::MultiFab lhs(getSlices(lev, WhichSlice::This), amrex::make_alias, Comps[WhichSlice::This]["Psi"], 1); - InterpolateFromLev0toLev1(geom, lev, "rho", islice); + InterpolateFromLev0toLev1(geom, lev, "rho"); // calculating the right-hand side 1/episilon0 * -(rho-Jz/c) CopyToStagingArea(getSlices(lev,WhichSlice::This), SliceOperatorType::Assign, @@ -655,7 +606,7 @@ Fields::SolvePoissonExmByAndEypBx (amrex::Vector const& geom, Comps[WhichSlice::This]["rho"], lev); m_poisson_solver[lev]->StagingArea().mult(-1./phys_const.ep0); - InterpolateBoundaries(geom, lev, "Psi", islice); + InterpolateBoundaries(geom, lev, "Psi"); m_poisson_solver[lev]->SolvePoissonEquation(lhs); /* ---------- Transverse FillBoundary Psi ---------- */ @@ -663,7 +614,7 @@ Fields::SolvePoissonExmByAndEypBx (amrex::Vector const& geom, lhs.FillBoundary(geom[lev].periodicity()); amrex::ParallelContext::pop(); - InterpolateFromLev0toLev1(geom, lev, "Psi", islice); + InterpolateFromLev0toLev1(geom, lev, "Psi"); /* Compute ExmBy and Eypbx from grad(-psi) */ TransverseDerivative( @@ -689,7 +640,7 @@ Fields::SolvePoissonExmByAndEypBx (amrex::Vector const& geom, void -Fields::SolvePoissonEz (amrex::Vector const& geom, const int lev, const int islice) +Fields::SolvePoissonEz (amrex::Vector const& geom, const int lev) { /* Solves Laplacian(Ez) = 1/(episilon0 *c0 )*(d_x(jx) + d_y(jy)) */ HIPACE_PROFILE("Fields::SolvePoissonEz()"); @@ -718,7 +669,7 @@ Fields::SolvePoissonEz (amrex::Vector const& geom, const int le SliceOperatorType::Add, Comps[WhichSlice::This]["jy"], 0, 1); - InterpolateBoundaries(geom, lev, "Ez", islice); + InterpolateBoundaries(geom, lev, "Ez"); // Solve Poisson equation. // The RHS is in the staging area of poisson_solver. // The LHS will be returned as lhs. @@ -727,7 +678,7 @@ Fields::SolvePoissonEz (amrex::Vector const& geom, const int le void Fields::SolvePoissonBx (amrex::MultiFab& Bx_iter, amrex::Vector const& geom, - const int lev, const int islice) + const int lev) { /* Solves Laplacian(Bx) = mu_0*(- d_y(jz) + d_z(jy) ) */ HIPACE_PROFILE("Fields::SolvePoissonBx()"); @@ -754,7 +705,7 @@ Fields::SolvePoissonBx (amrex::MultiFab& Bx_iter, amrex::Vector Comps[WhichSlice::Previous1]["jy"], Comps[WhichSlice::Next]["jy"]); - InterpolateBoundaries(geom, lev, "Bx", islice); + InterpolateBoundaries(geom, lev, "Bx"); // Solve Poisson equation. // The RHS is in the staging area of poisson_solver. // The LHS will be returned as lhs. @@ -763,7 +714,7 @@ Fields::SolvePoissonBx (amrex::MultiFab& Bx_iter, amrex::Vector void Fields::SolvePoissonBy (amrex::MultiFab& By_iter, amrex::Vector const& geom, - const int lev, const int islice) + const int lev) { /* Solves Laplacian(By) = mu_0*(d_x(jz) - d_z(jx) ) */ HIPACE_PROFILE("Fields::SolvePoissonBy()"); @@ -790,7 +741,7 @@ Fields::SolvePoissonBy (amrex::MultiFab& By_iter, amrex::Vector Comps[WhichSlice::Previous1]["jx"], Comps[WhichSlice::Next]["jx"]); - InterpolateBoundaries(geom, lev, "By", islice); + InterpolateBoundaries(geom, lev, "By"); // Solve Poisson equation. // The RHS is in the staging area of poisson_solver. // The LHS will be returned as lhs. @@ -798,7 +749,7 @@ Fields::SolvePoissonBy (amrex::MultiFab& By_iter, amrex::Vector } void -Fields::SolvePoissonBz (amrex::Vector const& geom, const int lev, const int islice) +Fields::SolvePoissonBz (amrex::Vector const& geom, const int lev) { /* Solves Laplacian(Bz) = mu_0*(d_y(jx) - d_x(jy)) */ HIPACE_PROFILE("Fields::SolvePoissonBz()"); @@ -827,7 +778,7 @@ Fields::SolvePoissonBz (amrex::Vector const& geom, const int le SliceOperatorType::Add, Comps[WhichSlice::This]["jy"], 0, 1); - InterpolateBoundaries(geom, lev, "Bz", islice); + InterpolateBoundaries(geom, lev, "Bz"); // Solve Poisson equation. // The RHS is in the staging area of m_poisson_solver. // The LHS will be returned as lhs. From fe1dcf4226cd59b58bc6b8becff86dee1d46933f Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Wed, 8 Sep 2021 19:50:33 +0200 Subject: [PATCH 05/15] use Parser --- src/Hipace.cpp | 4 +- src/diagnostics/Diagnostic.cpp | 10 ++-- src/diagnostics/OpenPMDWriter.cpp | 11 ++-- .../fft_poisson_solver/fft/WrapCuDST.cpp | 2 +- src/particles/BeamParticleContainer.cpp | 60 +++++++++---------- src/particles/MultiBeam.cpp | 12 ++-- src/particles/MultiPlasma.cpp | 6 +- src/particles/PlasmaParticleContainer.cpp | 34 +++++------ src/particles/PlasmaParticleContainerInit.cpp | 2 +- src/particles/profiles/GetInitialDensity.cpp | 8 +-- src/particles/profiles/GetInitialMomentum.cpp | 6 +- src/utils/AdaptiveTimeStep.cpp | 4 +- src/utils/GridCurrent.cpp | 10 ++-- src/utils/Parser.H | 28 +++++++-- 14 files changed, 107 insertions(+), 90 deletions(-) diff --git a/src/Hipace.cpp b/src/Hipace.cpp index a39ec89e17..18c1df1a02 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -165,9 +165,9 @@ Hipace::Hipace () : AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!m_explicit, "Mesh refinement + explicit solver is not yet" " supported! Please use hipace.bxby_solver = predictor-corrector"); amrex::Array loc_array; - pph.get("patch_lo", loc_array); + getWithParser(pph, "patch_lo", loc_array); for (int idim=0; idim diag_coarsen_arr{1,1,1}; // set all levels the same for now - ppd.query("coarsening", diag_coarsen_arr); + queryWithParser(ppd, "coarsening", diag_coarsen_arr); if(m_slice_dir == 0 || m_slice_dir == 1) { diag_coarsen_arr[m_slice_dir] = 1; } @@ -35,7 +35,7 @@ Diagnostic::Diagnostic (int nlev) "Coarsening ratio must be >= 1"); } - ppd.queryarr("field_data", m_comps_output); + queryWithParser(ppd, "field_data", m_comps_output); const amrex::Vector all_field_comps {"ExmBy", "EypBx", "Ez", "Bx", "By", "Bz", "jx", "jx_beam", "jy", "jy_beam", "jz", "jz_beam", "rho", "Psi"}; @@ -68,9 +68,9 @@ Diagnostic::Diagnostic (int nlev) amrex::ParmParse ppb("beams"); // read in all beam names amrex::Vector all_beam_names; - ppb.queryarr("names", all_beam_names); + queryWithParser(ppb, "names", all_beam_names); // read in which beam should be written to file - ppd.queryarr("beam_data", m_output_beam_names); + queryWithParser(ppd, "beam_data", m_output_beam_names); if(m_output_beam_names.empty()) { m_output_beam_names = all_beam_names; diff --git a/src/diagnostics/OpenPMDWriter.cpp b/src/diagnostics/OpenPMDWriter.cpp index 2115fdca1e..90c9866042 100644 --- a/src/diagnostics/OpenPMDWriter.cpp +++ b/src/diagnostics/OpenPMDWriter.cpp @@ -12,12 +12,12 @@ OpenPMDWriter::OpenPMDWriter () AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_real_names.size() == BeamIdx::nattribs, "List of real names in openPMD Writer class do not match BeamIdx::nattribs"); amrex::ParmParse pp("hipace"); - pp.query("file_prefix", m_file_prefix); - pp.query("openpmd_backend", m_openpmd_backend); + queryWithParser(pp, "file_prefix", m_file_prefix); + queryWithParser(pp, "openpmd_backend", m_openpmd_backend); // temporary workaround until openPMD-viewer gets fixed amrex::ParmParse ppd("diagnostic"); - ppd.query("openpmd_viewer_u_workaround", m_openpmd_viewer_workaround); + queryWithParser(ppd, "openpmd_viewer_u_workaround", m_openpmd_viewer_workaround); } void @@ -101,11 +101,12 @@ OpenPMDWriter::WriteFieldData ( auto meshes = iteration.meshes; // loop over field components - for ( int icomp = 0; icomp < varnames.size(); ++icomp ) + for (std::string fieldname : varnames) { + int icomp = Comps[WhichSlice::This][fieldname]; // "B" "x" (todo) // "Bx" "" (just for now) - openPMD::Mesh field = meshes[varnames[icomp]]; + openPMD::Mesh field = meshes[fieldname]; openPMD::MeshRecordComponent field_comp = field[openPMD::MeshRecordComponent::SCALAR]; // meta-data diff --git a/src/fields/fft_poisson_solver/fft/WrapCuDST.cpp b/src/fields/fft_poisson_solver/fft/WrapCuDST.cpp index 36fba3454d..123a4556f2 100644 --- a/src/fields/fft_poisson_solver/fft/WrapCuDST.cpp +++ b/src/fields/fft_poisson_solver/fft/WrapCuDST.cpp @@ -233,7 +233,7 @@ namespace AnyDST amrex::ParmParse pp("hipace"); dst_plan.use_small_dst = (std::max(real_size[0], real_size[1]) >= 511); - pp.query("use_small_dst", dst_plan.use_small_dst); + queryWithParser(pp, "use_small_dst", dst_plan.use_small_dst); if(!dst_plan.use_small_dst) { const int nx = real_size[0]; diff --git a/src/particles/BeamParticleContainer.cpp b/src/particles/BeamParticleContainer.cpp index 707589c64a..f7e70f2767 100644 --- a/src/particles/BeamParticleContainer.cpp +++ b/src/particles/BeamParticleContainer.cpp @@ -10,7 +10,7 @@ namespace PhysConst phys_const = get_phys_const(); std::string element = "electron"; - pp.query("element", element); + queryWithParser(pp, "element", element); if (element == "electron"){ charge = -phys_const.q_e; mass = phys_const.m_e; @@ -31,22 +31,22 @@ BeamParticleContainer::ReadParameters () amrex::ParmParse pp(m_name); QueryElementSetChargeMass(pp, m_charge, m_mass); // Overwrite element's charge and mass if user specifies them explicitly - pp.query("charge", m_charge); - pp.query("mass", m_mass); + queryWithParser(pp, "charge", m_charge); + queryWithParser(pp, "mass", m_mass); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_mass != 0, "The beam particle mass must not be 0"); - pp.get("injection_type", m_injection_type); + getWithParser(pp, "injection_type", m_injection_type); amrex::Vector tmp_vector; - if (pp.queryarr("ppc", tmp_vector)){ + if (queryWithParser(pp, "ppc", tmp_vector)){ AMREX_ALWAYS_ASSERT(tmp_vector.size() == AMREX_SPACEDIM); for (int i=0; i= 1, "n_subcycles must be >= 1"); if (m_injection_type == "fixed_ppc" || m_injection_type == "from_file"){ AMREX_ALWAYS_ASSERT_WITH_MESSAGE( (m_dx_per_dzeta == 0.) && (m_dy_per_dzeta == 0.) @@ -63,12 +63,12 @@ BeamParticleContainer::InitData (const amrex::Geometry& geom) if (m_injection_type == "fixed_ppc") { amrex::ParmParse pp(m_name); - pp.get("zmin", m_zmin); - pp.get("zmax", m_zmax); - pp.get("radius", m_radius); - pp.query("min_density", m_min_density); + getWithParser(pp, "zmin", m_zmin); + getWithParser(pp, "zmax", m_zmax); + getWithParser(pp, "radius", m_radius); + queryWithParser(pp, "min_density", m_min_density); amrex::Vector random_ppc {false, false, false}; - pp.queryarr("random_ppc", random_ppc); + queryWithParser(pp, "random_ppc", random_ppc); const GetInitialDensity get_density(m_name); const GetInitialMomentum get_momentum(m_name); InitBeamFixedPPC(m_ppc, get_density, get_momentum, geom, m_zmin, @@ -81,25 +81,25 @@ BeamParticleContainer::InitData (const amrex::Geometry& geom) bool can = false; amrex::Real zmin = 0, zmax = 0; std::string profile = "gaussian"; - pp.query("profile", profile); + queryWithParser(pp, "profile", profile); if (profile == "can") { can = true; - pp.get("zmin", zmin); - pp.get("zmax", zmax); + getWithParser(pp, "zmin", zmin); + getWithParser(pp, "zmax", zmax); } else if (profile == "gaussian") { } else { amrex::Abort("Only gaussian and can are supported with fixed_weight beam injection"); } - pp.get("position_mean", loc_array); + getWithParser(pp, "position_mean", loc_array); for (int idim=0; idim beam_names) { amrex::ParmParse pp("beams"); std::string all_input_file = ""; - if(!pp.query("all_from_file", all_input_file)) { + if(!queryWithParser(pp, "all_from_file", all_input_file)) { return; } int iteration = 0; - const bool multi_iteration = pp.query("iteration", iteration); + const bool multi_iteration = queryWithParser(pp, "iteration", iteration); amrex::Real plasma_density = 0; - const bool multi_plasma_density = pp.query("plasma_density", plasma_density); + const bool multi_plasma_density = queryWithParser(pp, "plasma_density", plasma_density); std::vector file_coordinates_xyz; - const bool multi_file_coordinates_xyz = pp.queryarr("file_coordinates_xyz", - file_coordinates_xyz); + const bool multi_file_coordinates_xyz = queryWithParser(pp, "file_coordinates_xyz", + file_coordinates_xyz); for( std::string name : beam_names ) { amrex::ParmParse pp_beam(name); diff --git a/src/particles/MultiPlasma.cpp b/src/particles/MultiPlasma.cpp index 1f94db540b..ebb6e17fff 100644 --- a/src/particles/MultiPlasma.cpp +++ b/src/particles/MultiPlasma.cpp @@ -8,9 +8,9 @@ MultiPlasma::MultiPlasma (amrex::AmrCore* amr_core) { amrex::ParmParse pp("plasmas"); - pp.getarr("names", m_names); - pp.query("adaptive_density", m_adaptive_density); - pp.query("sort_bin_size", m_sort_bin_size); + getWithParser(pp, "names", m_names); + queryWithParser(pp, "adaptive_density", m_adaptive_density); + queryWithParser(pp, "sort_bin_size", m_sort_bin_size); if (m_names[0] == "no_plasma") return; m_nplasmas = m_names.size(); diff --git a/src/particles/PlasmaParticleContainer.cpp b/src/particles/PlasmaParticleContainer.cpp index 909faa82e1..7e1bbc0b68 100644 --- a/src/particles/PlasmaParticleContainer.cpp +++ b/src/particles/PlasmaParticleContainer.cpp @@ -16,7 +16,7 @@ PlasmaParticleContainer::ReadParameters () amrex::ParmParse pp(m_name); std::string element = ""; amrex::Real mass_Da = 0; - pp.query("element", element); + queryWithParser(pp, "element", element); if (element == "electron") { m_charge = -phys_const.q_e; m_mass = phys_const.m_e; @@ -32,17 +32,17 @@ PlasmaParticleContainer::ReadParameters () AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mass_Da != 0, "Unknown Element"); } - pp.query("mass_Da", mass_Da); + queryWithParser(pp, "mass_Da", mass_Da); if(mass_Da != 0) { m_mass = phys_const.m_p * mass_Da / 1.007276466621; } - pp.query("mass", m_mass); + queryWithParser(pp, "mass", m_mass); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_mass != 0, "The plasma particle mass must be specified"); - bool ion_lev_specified = pp.query("initial_ion_level", m_init_ion_lev); + bool ion_lev_specified = queryWithParser(pp, "initial_ion_level", m_init_ion_lev); m_can_ionize = pp.contains("ionization_product"); - pp.query("can_ionize", m_can_ionize); + queryWithParser(pp, "can_ionize", m_can_ionize); if(m_can_ionize) { m_neutralize_background = false; // change default AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!Hipace::GetInstance().m_normalized_units, @@ -50,11 +50,11 @@ PlasmaParticleContainer::ReadParameters () AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_init_ion_lev >= 0, "The initial Ion level must be specified"); } - pp.query("neutralize_background", m_neutralize_background); + queryWithParser(pp, "neutralize_background", m_neutralize_background); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!m_can_ionize || !m_neutralize_background, "Cannot use neutralize_background for Ion plasma"); - if(!pp.query("charge", m_charge)) { + if(!queryWithParser(pp, "charge", m_charge)) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_charge != 0, "The plasma particle charge must be specified"); } @@ -62,18 +62,18 @@ PlasmaParticleContainer::ReadParameters () if(ion_lev_specified && !m_can_ionize) { m_charge *= m_init_ion_lev; } - pp.query("ionization_product", m_product_name); - pp.query("density", m_density); - pp.query("level", m_level); - pp.query("radius", m_radius); - pp.query("hollow_core_radius", m_hollow_core_radius); + queryWithParser(pp, "ionization_product", m_product_name); + queryWithParser(pp, "density", m_density); + queryWithParser(pp, "level", m_level); + queryWithParser(pp, "radius", m_radius); + queryWithParser(pp, "hollow_core_radius", m_hollow_core_radius); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_hollow_core_radius < m_radius, "The hollow core plasma radius must not be smaller than the " "plasma radius itself"); - pp.query("parabolic_curvature", m_parabolic_curvature); - pp.query("max_qsa_weighting_factor", m_max_qsa_weighting_factor); + queryWithParser(pp, "parabolic_curvature", m_parabolic_curvature); + queryWithParser(pp, "max_qsa_weighting_factor", m_max_qsa_weighting_factor); amrex::Vector tmp_vector; - if (pp.queryarr("ppc", tmp_vector)){ + if (queryWithParser(pp, "ppc", tmp_vector)){ AMREX_ALWAYS_ASSERT_WITH_MESSAGE(tmp_vector.size() == AMREX_SPACEDIM-1, "ppc is only specified in transverse directions for plasma particles, " "it is 1 in the longitudinal direction z. " @@ -81,12 +81,12 @@ PlasmaParticleContainer::ReadParameters () for (int i=0; i loc_array; - if (pp.query("u_mean", loc_array)) { + if (queryWithParser(pp, "u_mean", loc_array)) { for (int idim=0; idim < AMREX_SPACEDIM; ++idim) { m_u_mean[idim] = loc_array[idim]; } } - if (pp.query("u_std", loc_array)) { + if (queryWithParser(pp, "u_std", loc_array)) { for (int idim=0; idim < AMREX_SPACEDIM; ++idim) { m_u_std[idim] = loc_array[idim]; } diff --git a/src/particles/PlasmaParticleContainerInit.cpp b/src/particles/PlasmaParticleContainerInit.cpp index 85f827d168..c295a6ae27 100644 --- a/src/particles/PlasmaParticleContainerInit.cpp +++ b/src/particles/PlasmaParticleContainerInit.cpp @@ -207,7 +207,7 @@ InitIonizationModule (const amrex::Geometry& geom, m_product_pc = product_pc; amrex::ParmParse pp(m_name); std::string physical_element; - pp.get("element", physical_element); + getWithParser(pp, "element", physical_element); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(ion_map_ids.count(physical_element) != 0, "There are no ionization energies available for this element. " "Please update src/utils/IonizationEnergiesTable.H using write_atomic_data_cpp.py"); diff --git a/src/particles/profiles/GetInitialDensity.cpp b/src/particles/profiles/GetInitialDensity.cpp index 2a3b9d843a..21b17e0ce7 100644 --- a/src/particles/profiles/GetInitialDensity.cpp +++ b/src/particles/profiles/GetInitialDensity.cpp @@ -5,8 +5,8 @@ GetInitialDensity::GetInitialDensity (const std::string& name) { amrex::ParmParse pp(name); std::string profile; - pp.get("density", m_density); - pp.get("profile", profile); + getWithParser(pp, "density", m_density); + getWithParser(pp, "profile", profile); if (profile == "gaussian") { m_profile = BeamProfileType::Gaussian; @@ -18,12 +18,12 @@ GetInitialDensity::GetInitialDensity (const std::string& name) if (m_profile == BeamProfileType::Gaussian) { amrex::Array loc_array; - if (pp.query("position_mean", loc_array)) { + if (queryWithParser(pp, "position_mean", loc_array)) { for (int idim=0; idim < AMREX_SPACEDIM; ++idim) { m_position_mean[idim] = loc_array[idim]; } } - if (pp.query("position_std", loc_array)) { + if (queryWithParser(pp, "position_std", loc_array)) { for (int idim=0; idim < AMREX_SPACEDIM; ++idim) { m_position_std[idim] = loc_array[idim]; } diff --git a/src/particles/profiles/GetInitialMomentum.cpp b/src/particles/profiles/GetInitialMomentum.cpp index 9cf4347fab..b04b4247b2 100644 --- a/src/particles/profiles/GetInitialMomentum.cpp +++ b/src/particles/profiles/GetInitialMomentum.cpp @@ -9,18 +9,18 @@ GetInitialMomentum::GetInitialMomentum (const std::string& name) if (m_momentum_profile == BeamMomentumType::Gaussian) { amrex::Array loc_array; - if (pp.query("u_mean", loc_array)) { + if (queryWithParser(pp, "u_mean", loc_array)) { for (int idim=0; idim < AMREX_SPACEDIM; ++idim) { m_u_mean[idim] = loc_array[idim]; } } - if (pp.query("u_std", loc_array)) { + if (queryWithParser(pp, "u_std", loc_array)) { for (int idim=0; idim < AMREX_SPACEDIM; ++idim) { m_u_std[idim] = loc_array[idim]; } } bool do_symmetrize = false; - pp.query("do_symmetrize", do_symmetrize); + queryWithParser(pp, "do_symmetrize", do_symmetrize); if (do_symmetrize) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE( std::fabs(m_u_mean[0]) +std::fabs(m_u_mean[1]) < std::numeric_limits::epsilon(), diff --git a/src/utils/AdaptiveTimeStep.cpp b/src/utils/AdaptiveTimeStep.cpp index d601f08bbd..7796861381 100644 --- a/src/utils/AdaptiveTimeStep.cpp +++ b/src/utils/AdaptiveTimeStep.cpp @@ -17,8 +17,8 @@ struct WhichDouble { AdaptiveTimeStep::AdaptiveTimeStep () { amrex::ParmParse ppa("hipace"); - ppa.query("do_adaptive_time_step", m_do_adaptive_time_step); - ppa.query("nt_per_omega_betatron", m_nt_per_omega_betatron); + queryWithParser(ppa, "do_adaptive_time_step", m_do_adaptive_time_step); + queryWithParser(ppa, "nt_per_omega_betatron", m_nt_per_omega_betatron); } #ifdef AMREX_USE_MPI diff --git a/src/utils/GridCurrent.cpp b/src/utils/GridCurrent.cpp index eb8f0e1440..125c2ba4d8 100644 --- a/src/utils/GridCurrent.cpp +++ b/src/utils/GridCurrent.cpp @@ -7,14 +7,14 @@ GridCurrent::GridCurrent () { amrex::ParmParse pp("grid_current"); - if (pp.query("use_grid_current", m_use_grid_current) ) { - pp.get("peak_current_density", m_peak_current_density); + if (queryWithParser(pp, "use_grid_current", m_use_grid_current) ) { + getWithParser(pp, "peak_current_density", m_peak_current_density); amrex::Array loc_array; - pp.get("position_mean", loc_array); + getWithParser(pp, "position_mean", loc_array); for (int idim=0; idim < AMREX_SPACEDIM; ++idim) m_position_mean[idim] = loc_array[idim]; - pp.get("position_std", loc_array); + getWithParser(pp, "position_std", loc_array); for (int idim=0; idim < AMREX_SPACEDIM; ++idim) m_position_std[idim] = loc_array[idim]; - pp.query("finest_level", m_finest_level); + queryWithParser(pp, "finest_level", m_finest_level); } } diff --git a/src/utils/Parser.H b/src/utils/Parser.H index 32776b1720..a9348aae92 100755 --- a/src/utils/Parser.H +++ b/src/utils/Parser.H @@ -21,7 +21,7 @@ inline void getWithParser (const amrex::ParmParse& pp, char const * const str, T& val); template -inline bool +inline int queryWithParser (const amrex::ParmParse& pp, char const * const str, T& val); inline int @@ -97,9 +97,8 @@ makeParser (std::string const& parse_function, amrex::Vector const& parser.setConstant(s, usr_constant->second); continue; } - - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(recursive_symbols.count(*it)==0, - "Expressions contains recursive symbol " + *it); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(recursive_symbols.count(s)==0, + "Expressions contains recursive symbol " + s); double v; recursive_symbols.insert(s); const bool is_input = queryWithParser(pp_my_constants, s.c_str(), v); @@ -153,6 +152,16 @@ fillWithParser (std::string const& str, int& val) { val = safeCastToInt(std::round(makeParser(str, {}).compileHost<0>()()),str); } +inline void +fillWithParser (std::string const& str, bool& val) { + val = makeParser(str, {}).compileHost<0>()(); +} + +inline void +fillWithParser (std::string const& str, std::string& val) { + val = str; +} + template inline void fillWithParserArr (std::vector const& str_arr, T& val) { @@ -161,6 +170,7 @@ fillWithParserArr (std::vector const& str_arr, T& val) { str += s; } fillWithParser(str, val); + amrex::Print() << val; } template @@ -170,6 +180,7 @@ fillWithParserArr (std::vector const& str_arr, std::vector& val_ val_arr.resize(n); for (auto i=0ul ; i != n ; ++i) { fillWithParser(str_arr[i], val_arr[i]); + amrex::Print() << val_arr[i] << " "; } } @@ -186,6 +197,7 @@ fillWithParserArr (std::vector const& str_arr, std::array& AMREX_ALWAYS_ASSERT( n == size ); for (auto i=0ul ; i != n ; ++i) { fillWithParser(str_arr[i], val_arr[i]); + amrex::Print() << val_arr[i] << " "; } } @@ -195,19 +207,23 @@ getWithParser (const amrex::ParmParse& pp, char const * const str, T& val) { std::vector f; pp.getarr(str, f); + amrex::Print() << str << " "; fillWithParserArr(f, val); + amrex::Print() << '\n'; } template -inline bool +inline int queryWithParser (const amrex::ParmParse& pp, char const * const str, T& val) { std::vector f; - bool is_specified = pp.queryarr(str, f); + const int is_specified = pp.queryarr(str, f); + amrex::Print() << str << " "; if (is_specified) { fillWithParserArr(f, val); } + amrex::Print() << '\n'; return is_specified; } From abdc747e64b38e0d7f125b535f0a9e15ffa33094 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Wed, 8 Sep 2021 19:56:06 +0200 Subject: [PATCH 06/15] fix Fields.cpp --- src/fields/Fields.cpp | 199 ++++++++++++++++++++++++++---------------- 1 file changed, 124 insertions(+), 75 deletions(-) diff --git a/src/fields/Fields.cpp b/src/fields/Fields.cpp index 0a59b61298..a334093482 100644 --- a/src/fields/Fields.cpp +++ b/src/fields/Fields.cpp @@ -233,12 +233,14 @@ Fields::LongitudinalDerivative (const amrex::MultiFab& src1, const amrex::MultiF void -Fields::Copy (int lev, int i_slice, FieldCopyType copy_type, int slice_comp, int full_comp, - int ncomp, amrex::FArrayBox& fab, int slice_dir, amrex::Geometry geom) +Fields::Copy (const int lev, const int i_slice, const int slice_comp, const int full_comp, + const amrex::Gpu::DeviceVector& diag_comps_vect, + const int ncomp, amrex::FArrayBox& fab, const int slice_dir, + const amrex::IntVect diag_coarsen, const amrex::Geometry geom) { using namespace amrex::literals; HIPACE_PROFILE("Fields::Copy()"); - auto& slice_mf = m_slices[lev][WhichSlice::This]; // copy from/to the current slice + auto& slice_mf = m_slices[lev][WhichSlice::This]; // copy from the current slice amrex::Array4 slice_array; // There is only one Box. for (amrex::MFIter mfi(slice_mf); mfi.isValid(); ++mfi) { auto& slice_fab = slice_mf[mfi]; @@ -249,64 +251,92 @@ Fields::Copy (int lev, int i_slice, FieldCopyType copy_type, int slice_comp, int // slice_array's longitude index is i_slice. } + const int full_array_z = i_slice / diag_coarsen[2]; + const amrex::IntVect ncells_global = geom.Domain().length(); + amrex::Box const& vbx = fab.box(); - if (vbx.smallEnd(Direction::z) <= i_slice and - vbx.bigEnd (Direction::z) >= i_slice) + if (vbx.smallEnd(Direction::z) <= full_array_z and + vbx.bigEnd (Direction::z) >= full_array_z and + ( i_slice % diag_coarsen[2] == diag_coarsen[2]/2 or + ( i_slice == ncells_global[2] - 1 and + ( ncells_global[2] - 1 ) % diag_coarsen[2] < diag_coarsen[2]/2 ))) { amrex::Box copy_box = vbx; - copy_box.setSmall(Direction::z, i_slice); - copy_box.setBig (Direction::z, i_slice); + copy_box.setSmall(Direction::z, full_array_z); + copy_box.setBig (Direction::z, full_array_z); amrex::Array4 const& full_array = fab.array(); - const amrex::IntVect ncells_global = geom.Domain().length(); - const bool nx_even = ncells_global[0] % 2 == 0; - const bool ny_even = ncells_global[1] % 2 == 0; + const int even_slice_x = ncells_global[0] % 2 == 0 and slice_dir == 0; + const int even_slice_y = ncells_global[1] % 2 == 0 and slice_dir == 1; - if (copy_type == FieldCopyType::FtoS) { - amrex::ParallelFor(copy_box, ncomp, - [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - slice_array(i,j,k,n+slice_comp) = full_array(i,j,k,n+full_comp); - }); - } else { - amrex::ParallelFor(copy_box, ncomp, - [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - if (slice_dir ==-1 /* 3D data */){ - full_array(i,j,k,n+full_comp) = slice_array(i,j,k,n+slice_comp); - } else if (slice_dir == 0 /* yz slice */){ - full_array(i,j,k,n+full_comp) = - nx_even ? 0.5_rt * (slice_array(i-1,j,k,n+slice_comp) + - slice_array(i,j,k,n+slice_comp)) - : slice_array(i,j,k,n+slice_comp); - } else /* slice_dir == 1, xz slice */{ - full_array(i,j,k,n+full_comp) = - ny_even ? 0.5_rt * ( slice_array(i,j-1,k,n+slice_comp) + - slice_array(i,j,k,n+slice_comp)) - : slice_array(i,j,k,n+slice_comp); + const int coarse_x = diag_coarsen[0]; + const int coarse_y = diag_coarsen[1]; + + const int ncells_x = ncells_global[0]; + const int ncells_y = ncells_global[1]; + + const int *diag_comps = diag_comps_vect.data(); + + amrex::ParallelFor(copy_box, ncomp, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + const int m = n[diag_comps]; + + // coarsening in slice direction is always 1 + const int i_c_start = amrex::min(i*coarse_x +(coarse_x-1)/2 -even_slice_x, ncells_x-1); + const int i_c_stop = amrex::min(i*coarse_x +coarse_x/2+1, ncells_x); + const int j_c_start = amrex::min(j*coarse_y +(coarse_y-1)/2 -even_slice_y, ncells_y-1); + const int j_c_stop = amrex::min(j*coarse_y +coarse_y/2+1, ncells_y); + + amrex::Real field_value = 0._rt; + int n_values = 0; + + for (int j_c = j_c_start; j_c != j_c_stop; ++j_c) { + for (int i_c = i_c_start; i_c != i_c_stop; ++i_c) { + field_value += slice_array(i_c, j_c, i_slice, m+slice_comp); + ++n_values; } - }); - } + } + + full_array(i,j,k,n+full_comp) = field_value / amrex::max(n_values,1); + }); } } void -Fields::ShiftSlices (int lev) +Fields::ShiftSlices (int nlev, int islice, amrex::Geometry geom, amrex::Real patch_lo, + amrex::Real patch_hi) { HIPACE_PROFILE("Fields::ShiftSlices()"); + + for (int lev = 0; lev < nlev; ++lev) { + + // skip all slices which are not existing on level 1 + if (lev == 1) { + // use geometry of coarse grid to determine whether slice is to be solved + const amrex::Real* problo = geom.ProbLo(); + const amrex::Real* dx = geom.CellSize(); + const amrex::Real pos = (islice+0.5)*dx[2]+problo[2]; + if (pos < patch_lo || pos > patch_hi) continue; + } + + // shift Bx, By amrex::MultiFab::Copy( getSlices(lev, WhichSlice::Previous2), getSlices(lev, WhichSlice::Previous1), Comps[WhichSlice::Previous1]["Bx"], Comps[WhichSlice::Previous2]["Bx"], 2, m_slices_nguards); + // shift Ez, Bx, By, Bz, jx, jx_beam, jy, jy_beam amrex::MultiFab::Copy( getSlices(lev, WhichSlice::Previous1), getSlices(lev, WhichSlice::This), - Comps[WhichSlice::This]["Bx"], Comps[WhichSlice::Previous1]["Bx"], - 2, m_slices_nguards); + Comps[WhichSlice::This]["Ez"], Comps[WhichSlice::Previous1]["Ez"], + 8, m_slices_nguards); + // shift rho, Psi amrex::MultiFab::Copy( getSlices(lev, WhichSlice::Previous1), getSlices(lev, WhichSlice::This), - Comps[WhichSlice::This]["jx"], Comps[WhichSlice::Previous1]["jx"], - 4, m_slices_nguards); + Comps[WhichSlice::This]["rho"], Comps[WhichSlice::Previous1]["rho"], + 2, m_slices_nguards); + } } void @@ -340,13 +370,14 @@ Fields::AddBeamCurrents (const int lev, const int which_slice) void Fields::InterpolateBoundaries (amrex::Vector const& geom, const int lev, - std::string component) + std::string component, const int islice) { // To solve a Poisson equation with non-zero Dirichlet boundary conditions, the source term // must be corrected at the outmost grid points in x by -field_value_at_guard_cell / dx^2 and // in y by -field_value_at_guard_cell / dy^2, where dx and dy are those of the fine grid // This follows Van Loan, C. (1992). Computational frameworks for the fast Fourier transform. // Page 254 ff. + // The interpolation is done in second order transversely and linearly in longitudinal direction HIPACE_PROFILE("Fields::InterpolateBoundaries()"); if (lev == 0) return; // only interpolate boundaries to lev 1 @@ -355,11 +386,18 @@ Fields::InterpolateBoundaries (amrex::Vector const& geom, const const auto dx = geom[lev].CellSizeArray(); const auto plo_coarse = geom[lev-1].ProbLoArray(); const auto dx_coarse = geom[lev-1].CellSizeArray(); - constexpr int interp_order = 2; +constexpr int interp_order = 2; + + // get relative position of fine grid slice between coarse grids for longitudinal lin. interpol. + const amrex::Real z = plo_coarse[2] + (islice+0.5_rt)*dx[2]; + const int idz_coarse = (z-plo_coarse[2])/dx_coarse[2]; + const amrex::Real rel_z = (z - (plo_coarse[2] + (idz_coarse)*dx_coarse[2])) / dx_coarse[2]; // get level 0 for interpolation to source term of level 1 amrex::MultiFab lhs_coarse(getSlices(lev-1, WhichSlice::This), amrex::make_alias, Comps[WhichSlice::This][component], 1); + amrex::MultiFab lhs_coarse_prev(getSlices(lev-1, WhichSlice::Previous1), amrex::make_alias, + Comps[WhichSlice::Previous1][component], 1); amrex::FArrayBox& lhs_fab = lhs_coarse[0]; amrex::Box lhs_bx = lhs_fab.box(); lhs_bx.grow({-m_slices_nguards[0], -m_slices_nguards[1], 0}); @@ -384,7 +422,8 @@ Fields::InterpolateBoundaries (amrex::Vector const& geom, const const int nx_fine_high = big[0]; const int ny_fine_high = big[1]; amrex::Array4 data_array = m_poisson_solver[lev]->StagingArea().array(mfi); - amrex::Array4 data_array_coarse = lhs_coarse.array(mfi); + amrex::Array4 arr_coarse = lhs_coarse.array(mfi); + amrex::Array4 arr_coarse_prev = lhs_coarse_prev.array(mfi); // Loop over the valid indices on the fine grid and interpolate the value of the coarse grid // at the location of the guard cell on the fine grid to the first/last valid grid point on @@ -414,23 +453,22 @@ Fields::InterpolateBoundaries (amrex::Vector const& geom, const // sx_cell shape factor along x const amrex::Real xmid = (x - plo_coarse[0])/dx_coarse[0]; amrex::Real sx_cell[interp_order + 1]; - const int j_cell = compute_shape_factor(sx_cell, - xmid-0.5_rt); + const int j_cell = compute_shape_factor(sx_cell, xmid-0.5_rt); // y direction const amrex::Real ymid = (y - plo_coarse[1])/dx_coarse[1]; amrex::Real sy_cell[interp_order + 1]; - const int k_cell = compute_shape_factor(sy_cell, - ymid-0.5_rt); + const int k_cell = compute_shape_factor(sy_cell, ymid-0.5_rt); amrex::Real boundary_value = 0.0_rt; // add interpolated contribution to boundary value for (int iy=0; iy<=interp_order; iy++){ for (int ix=0; ix<=interp_order; ix++){ - boundary_value += data_array_coarse(lo_coarse[0]+j_cell+ix, - lo_coarse[1]+k_cell+iy, - lo_coarse[2]) - *sx_cell[ix]*sy_cell[iy]; + boundary_value += sx_cell[ix]*sy_cell[iy]* + ((1.0_rt-rel_z)*arr_coarse(lo_coarse[0]+j_cell+ix, + lo_coarse[1]+k_cell+iy, lo_coarse[2]) + + rel_z*arr_coarse_prev(lo_coarse[0]+j_cell+ix, + lo_coarse[1]+k_cell+iy, lo_coarse[2])); } } @@ -455,23 +493,22 @@ Fields::InterpolateBoundaries (amrex::Vector const& geom, const // sx_cell shape factor along x const amrex::Real xmid = (x - plo_coarse[0])/dx_coarse[0]; amrex::Real sx_cell[interp_order + 1]; - const int j_cell = compute_shape_factor(sx_cell, - xmid-0.5_rt); + const int j_cell = compute_shape_factor(sx_cell, xmid-0.5_rt); // y direction const amrex::Real ymid = (y - plo_coarse[1])/dx_coarse[1]; amrex::Real sy_cell[interp_order + 1]; - const int k_cell = compute_shape_factor(sy_cell, - ymid-0.5_rt); + const int k_cell = compute_shape_factor(sy_cell, ymid-0.5_rt); amrex::Real boundary_value = 0.0_rt; // add interpolated contribution to boundary value for (int iy=0; iy<=interp_order; iy++){ for (int ix=0; ix<=interp_order; ix++){ - boundary_value += data_array_coarse(lo_coarse[0]+j_cell+ix, - lo_coarse[1]+k_cell+iy, - lo_coarse[2]) - *sx_cell[ix]*sy_cell[iy]; + boundary_value += sx_cell[ix]*sy_cell[iy]* + ((1.0_rt-rel_z)*arr_coarse(lo_coarse[0]+j_cell+ix, + lo_coarse[1]+k_cell+iy, lo_coarse[2]) + + rel_z*arr_coarse_prev(lo_coarse[0]+j_cell+ix, + lo_coarse[1]+k_cell+iy, lo_coarse[2])); } } @@ -485,11 +522,12 @@ Fields::InterpolateBoundaries (amrex::Vector const& geom, const void Fields::InterpolateFromLev0toLev1 (amrex::Vector const& geom, const int lev, - std::string component) + std::string component, const int islice) { // This function interpolates values from the coarse to the fine grid with second order. // This is required for rho to fix the incomplete deposition close to the boundary and for Psi // to fill the guard cell, which is needed for the transverse derivative + // The interpolation is done in second order transversely and linearly in longitudinal direction HIPACE_PROFILE("Fields::InterpolateFromLev0toLev1()"); if (lev == 0) return; // only interpolate boundaries to lev 1 @@ -500,9 +538,16 @@ Fields::InterpolateFromLev0toLev1 (amrex::Vector const& geom, c const auto dx_coarse = geom[lev-1].CellSizeArray(); constexpr int interp_order = 2; + // get relative position of fine grid slice between coarse grids for longitudinal lin. interpol. + const amrex::Real z = plo_coarse[2] + (islice+0.5_rt)*dx[2]; + const int idz_coarse = (z-plo_coarse[2])/dx_coarse[2]; + const amrex::Real rel_z = (z - (plo_coarse[2] + (idz_coarse)*dx_coarse[2])) / dx_coarse[2]; + // get level 0 array amrex::MultiFab lhs_coarse(getSlices(lev-1, WhichSlice::This), amrex::make_alias, Comps[WhichSlice::This][component], 1); + amrex::MultiFab lhs_coarse_prev(getSlices(lev-1, WhichSlice::Previous1), amrex::make_alias, + Comps[WhichSlice::Previous1][component], 1); amrex::FArrayBox& lhs_fab = lhs_coarse[0]; amrex::Box lhs_bx = lhs_fab.box(); // lhs_bx should only have valid cells @@ -536,7 +581,8 @@ Fields::InterpolateFromLev0toLev1 (amrex::Vector const& geom, c const int y_range = (component == "rho") ? m_slices_nguards[1] : 1; amrex::Array4 data_array = lhs_fine.array(mfi); - amrex::Array4 data_array_coarse = lhs_coarse.array(mfi); + amrex::Array4 arr_coarse = lhs_coarse.array(mfi); + amrex::Array4 arr_coarse_prev = lhs_coarse_prev.array(mfi); // Loop over the valid indices on the fine grid and interpolate the value of the coarse grid amrex::ParallelFor( @@ -568,9 +614,11 @@ Fields::InterpolateFromLev0toLev1 (amrex::Vector const& geom, c // sum interpolated contributions for (int iy=0; iy<=interp_order; iy++){ for (int ix=0; ix<=interp_order; ix++){ - coarse_value += data_array_coarse(lo_coarse[0]+j_cell+ix, - lo_coarse[1]+k_cell+iy, - lo_coarse[2])*sx_cell[ix]*sy_cell[iy]; + coarse_value += sx_cell[ix]*sy_cell[iy]* + ((1.0_rt-rel_z)*arr_coarse(lo_coarse[0]+j_cell+ix, + lo_coarse[1]+k_cell+iy, lo_coarse[2]) + + rel_z*arr_coarse_prev(lo_coarse[0]+j_cell+ix, + lo_coarse[1]+k_cell+iy, lo_coarse[2])); } } @@ -581,9 +629,10 @@ Fields::InterpolateFromLev0toLev1 (amrex::Vector const& geom, c } } + void Fields::SolvePoissonExmByAndEypBx (amrex::Vector const& geom, - const MPI_Comm& m_comm_xy, const int lev) + const MPI_Comm& m_comm_xy, const int lev, const int islice) { /* Solves Laplacian(Psi) = 1/episilon0 * -(rho-Jz/c) and * calculates Ex-c By, Ey + c Bx from grad(-Psi) @@ -596,7 +645,7 @@ Fields::SolvePoissonExmByAndEypBx (amrex::Vector const& geom, amrex::MultiFab lhs(getSlices(lev, WhichSlice::This), amrex::make_alias, Comps[WhichSlice::This]["Psi"], 1); - InterpolateFromLev0toLev1(geom, lev, "rho"); + InterpolateFromLev0toLev1(geom, lev, "rho", islice); // calculating the right-hand side 1/episilon0 * -(rho-Jz/c) CopyToStagingArea(getSlices(lev,WhichSlice::This), SliceOperatorType::Assign, @@ -606,7 +655,7 @@ Fields::SolvePoissonExmByAndEypBx (amrex::Vector const& geom, Comps[WhichSlice::This]["rho"], lev); m_poisson_solver[lev]->StagingArea().mult(-1./phys_const.ep0); - InterpolateBoundaries(geom, lev, "Psi"); + InterpolateBoundaries(geom, lev, "Psi", islice); m_poisson_solver[lev]->SolvePoissonEquation(lhs); /* ---------- Transverse FillBoundary Psi ---------- */ @@ -614,7 +663,7 @@ Fields::SolvePoissonExmByAndEypBx (amrex::Vector const& geom, lhs.FillBoundary(geom[lev].periodicity()); amrex::ParallelContext::pop(); - InterpolateFromLev0toLev1(geom, lev, "Psi"); + InterpolateFromLev0toLev1(geom, lev, "Psi", islice); /* Compute ExmBy and Eypbx from grad(-psi) */ TransverseDerivative( @@ -640,7 +689,7 @@ Fields::SolvePoissonExmByAndEypBx (amrex::Vector const& geom, void -Fields::SolvePoissonEz (amrex::Vector const& geom, const int lev) +Fields::SolvePoissonEz (amrex::Vector const& geom, const int lev, const int islice) { /* Solves Laplacian(Ez) = 1/(episilon0 *c0 )*(d_x(jx) + d_y(jy)) */ HIPACE_PROFILE("Fields::SolvePoissonEz()"); @@ -669,7 +718,7 @@ Fields::SolvePoissonEz (amrex::Vector const& geom, const int le SliceOperatorType::Add, Comps[WhichSlice::This]["jy"], 0, 1); - InterpolateBoundaries(geom, lev, "Ez"); + InterpolateBoundaries(geom, lev, "Ez", islice); // Solve Poisson equation. // The RHS is in the staging area of poisson_solver. // The LHS will be returned as lhs. @@ -678,7 +727,7 @@ Fields::SolvePoissonEz (amrex::Vector const& geom, const int le void Fields::SolvePoissonBx (amrex::MultiFab& Bx_iter, amrex::Vector const& geom, - const int lev) + const int lev, const int islice) { /* Solves Laplacian(Bx) = mu_0*(- d_y(jz) + d_z(jy) ) */ HIPACE_PROFILE("Fields::SolvePoissonBx()"); @@ -705,7 +754,7 @@ Fields::SolvePoissonBx (amrex::MultiFab& Bx_iter, amrex::Vector Comps[WhichSlice::Previous1]["jy"], Comps[WhichSlice::Next]["jy"]); - InterpolateBoundaries(geom, lev, "Bx"); + InterpolateBoundaries(geom, lev, "Bx", islice); // Solve Poisson equation. // The RHS is in the staging area of poisson_solver. // The LHS will be returned as lhs. @@ -714,7 +763,7 @@ Fields::SolvePoissonBx (amrex::MultiFab& Bx_iter, amrex::Vector void Fields::SolvePoissonBy (amrex::MultiFab& By_iter, amrex::Vector const& geom, - const int lev) + const int lev, const int islice) { /* Solves Laplacian(By) = mu_0*(d_x(jz) - d_z(jx) ) */ HIPACE_PROFILE("Fields::SolvePoissonBy()"); @@ -741,7 +790,7 @@ Fields::SolvePoissonBy (amrex::MultiFab& By_iter, amrex::Vector Comps[WhichSlice::Previous1]["jx"], Comps[WhichSlice::Next]["jx"]); - InterpolateBoundaries(geom, lev, "By"); + InterpolateBoundaries(geom, lev, "By", islice); // Solve Poisson equation. // The RHS is in the staging area of poisson_solver. // The LHS will be returned as lhs. @@ -749,7 +798,7 @@ Fields::SolvePoissonBy (amrex::MultiFab& By_iter, amrex::Vector } void -Fields::SolvePoissonBz (amrex::Vector const& geom, const int lev) +Fields::SolvePoissonBz (amrex::Vector const& geom, const int lev, const int islice) { /* Solves Laplacian(Bz) = mu_0*(d_y(jx) - d_x(jy)) */ HIPACE_PROFILE("Fields::SolvePoissonBz()"); @@ -778,7 +827,7 @@ Fields::SolvePoissonBz (amrex::Vector const& geom, const int le SliceOperatorType::Add, Comps[WhichSlice::This]["jy"], 0, 1); - InterpolateBoundaries(geom, lev, "Bz"); + InterpolateBoundaries(geom, lev, "Bz", islice); // Solve Poisson equation. // The RHS is in the staging area of m_poisson_solver. // The LHS will be returned as lhs. From c099d45836232d3d0c269278d5c6432f6fe11a30 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Wed, 8 Sep 2021 20:05:34 +0200 Subject: [PATCH 07/15] fix OpenPMDWriter.cpp --- src/diagnostics/OpenPMDWriter.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/diagnostics/OpenPMDWriter.cpp b/src/diagnostics/OpenPMDWriter.cpp index 90c9866042..534108ac13 100644 --- a/src/diagnostics/OpenPMDWriter.cpp +++ b/src/diagnostics/OpenPMDWriter.cpp @@ -101,12 +101,11 @@ OpenPMDWriter::WriteFieldData ( auto meshes = iteration.meshes; // loop over field components - for (std::string fieldname : varnames) + for ( int icomp = 0; icomp < varnames.size(); ++icomp ) { - int icomp = Comps[WhichSlice::This][fieldname]; // "B" "x" (todo) // "Bx" "" (just for now) - openPMD::Mesh field = meshes[fieldname]; + openPMD::Mesh field = meshes[varnames[icomp]]; openPMD::MeshRecordComponent field_comp = field[openPMD::MeshRecordComponent::SCALAR]; // meta-data From 9ac9910a0dcb19355960314a451e77ae71008977 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Wed, 8 Sep 2021 23:28:41 +0200 Subject: [PATCH 08/15] fix UB --- src/Hipace.cpp | 33 +---------- src/particles/profiles/GetInitialDensity.cpp | 2 +- src/particles/profiles/GetInitialMomentum.cpp | 2 +- src/utils/Parser.H | 59 +++++++++++++------ 4 files changed, 46 insertions(+), 50 deletions(-) diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 18c1df1a02..d47d7d049b 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -84,36 +84,9 @@ Hipace::Hipace () : amrex::ParmParse pph("hipace"); - double test_d = 0; - bool test_d_b = queryWithParser(pph, "test_d", test_d); - amrex::Print() << "##### Parser Test double " << test_d_b << ": " << test_d << '\n'; - - float test_f = 0; - bool test_f_b = queryWithParser(pph, "test_f", test_f); - amrex::Print() << "##### Parser Test float " << test_f_b << ": " << test_f << '\n'; - - int test_i = 0; - bool test_i_b = queryWithParser(pph, "test_i", test_i); - amrex::Print() << "##### Parser Test int " << test_i_b << ": " << test_i << '\n'; - - std::vector test_v{}; - bool test_v_b = queryWithParser(pph, "test_v", test_v); - amrex::Print() << "##### Parser Test vector " << test_v_b << ": "; - for (auto t : test_v) { - amrex::Print() << t << " "; - } - amrex::Print() << '\n'; - - std::array test_a{}; - bool test_a_b = queryWithParser(pph, "test_a", test_a); - amrex::Print() << "##### Parser Test array " << test_a_b << ": "; - for (auto t : test_a) { - amrex::Print() << t << " "; - } - amrex::Print() << '\n'; - - auto test_func = getFunctionWithParser<3>(pph, "test_func", {"x","y","z"}); - amrex::Print() << "##### Parser Test function " << test_func(1.,2.,3.) << '\n'; + amrex::Parser parser; + auto test_func = getFunctionWithParser<3>(pph, "function", parser, {"x","y","z"}); + amrex::Print() << "function : " << test_func(1.,2.,3.) << '\n'; queryWithParser(pph, "dt", m_dt); queryWithParser(pph, "verbose", m_verbose); diff --git a/src/particles/profiles/GetInitialDensity.cpp b/src/particles/profiles/GetInitialDensity.cpp index 21b17e0ce7..58cf2e7fa4 100644 --- a/src/particles/profiles/GetInitialDensity.cpp +++ b/src/particles/profiles/GetInitialDensity.cpp @@ -1,5 +1,5 @@ #include "GetInitialDensity.H" -#include +#include "utils/Parser.H" GetInitialDensity::GetInitialDensity (const std::string& name) { diff --git a/src/particles/profiles/GetInitialMomentum.cpp b/src/particles/profiles/GetInitialMomentum.cpp index b04b4247b2..41a8d6bb64 100644 --- a/src/particles/profiles/GetInitialMomentum.cpp +++ b/src/particles/profiles/GetInitialMomentum.cpp @@ -1,5 +1,5 @@ #include "GetInitialMomentum.H" -#include +#include "utils/Parser.H" GetInitialMomentum::GetInitialMomentum (const std::string& name) { diff --git a/src/utils/Parser.H b/src/utils/Parser.H index a9348aae92..bc0e6a940f 100755 --- a/src/utils/Parser.H +++ b/src/utils/Parser.H @@ -8,6 +8,7 @@ #include #include +#include #include #include #include @@ -51,14 +52,14 @@ safeCastToInt(const double x, const std::string& real_name) { return result; } -static amrex::Parser -makeParser (std::string const& parse_function, amrex::Vector const& varnames) +static amrex::Parser& +initParser (amrex::Parser& parser, amrex::Vector const& varnames) { // Since queryWithParser recursively calls this routine, keep track of symbols // in case an infinite recursion is found (a symbol's value depending on itself). static std::set recursive_symbols{}; - amrex::Parser parser(parse_function); + //amrex::Parser parser(parse_function); parser.registerVariables(varnames); std::set symbols = parser.symbols(); @@ -92,13 +93,15 @@ makeParser (std::string const& parse_function, amrex::Vector const& for (auto const& s : symbols) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(recursive_symbols.count(s)==0, + "Expressions contains recursive symbol " + s); + auto usr_constant = my_constants_cache.find(s); if (usr_constant != my_constants_cache.end()) { parser.setConstant(s, usr_constant->second); continue; } - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(recursive_symbols.count(s)==0, - "Expressions contains recursive symbol " + s); + double v; recursive_symbols.insert(s); const bool is_input = queryWithParser(pp_my_constants, s.c_str(), v); @@ -118,7 +121,6 @@ makeParser (std::string const& parse_function, amrex::Vector const& amrex::Abort("makeParser::Unknown symbol " + s); } - parser.print(); return parser; } @@ -126,40 +128,50 @@ template inline auto getFunctionWithParser (const amrex::ParmParse& pp, char const * const str, + amrex::Parser& parser, amrex::Vector const& varnames) { std::string func_str{}; std::vector f; pp.getarr(str, f); for (auto const& s : f) { - func_str += s; + func_str += s + ' '; } - return makeParser(func_str, varnames).compile(); + parser.define(func_str); + initParser(parser, varnames); + parser.print(); + return parser.compile(); } inline void fillWithParser (std::string const& str, double& val) { - val = makeParser(str, {}).compileHost<0>()(); + amrex::Parser parser(str); + val = initParser(parser, {}).compileHost<0>()(); } inline void fillWithParser (std::string const& str, float& val) { - val = static_cast(makeParser(str, {}).compileHost<0>()()); + amrex::Parser parser(str); + val = static_cast(initParser(parser, {}).compileHost<0>()()); } inline void fillWithParser (std::string const& str, int& val) { - val = safeCastToInt(std::round(makeParser(str, {}).compileHost<0>()()),str); + amrex::Parser parser(str); + val = safeCastToInt(std::round(initParser(parser, {}).compileHost<0>()()),str); } inline void fillWithParser (std::string const& str, bool& val) { - val = makeParser(str, {}).compileHost<0>()(); + amrex::Parser parser(str); + val = initParser(parser, {}).compileHost<0>()(); } inline void fillWithParser (std::string const& str, std::string& val) { val = str; + auto end_pos = std::remove(val.begin(), val.end(), ' '); + val.erase(end_pos, val.end()); } template @@ -167,10 +179,10 @@ inline void fillWithParserArr (std::vector const& str_arr, T& val) { std::string str{}; for (auto const& s : str_arr) { - str += s; + str += s + ' '; } fillWithParser(str, val); - amrex::Print() << val; + amrex::Print() << " " << val; } template @@ -180,7 +192,18 @@ fillWithParserArr (std::vector const& str_arr, std::vector& val_ val_arr.resize(n); for (auto i=0ul ; i != n ; ++i) { fillWithParser(str_arr[i], val_arr[i]); - amrex::Print() << val_arr[i] << " "; + amrex::Print() << " " << val_arr[i]; + } +} + +template +inline void +fillWithParserArr (std::vector const& str_arr, amrex::Vector& val_arr) { + auto const n = str_arr.size(); + val_arr.resize(n); + for (auto i=0ul ; i != n ; ++i) { + fillWithParser(str_arr[i], val_arr[i]); + amrex::Print() << " " << val_arr[i]; } } @@ -197,7 +220,7 @@ fillWithParserArr (std::vector const& str_arr, std::array& AMREX_ALWAYS_ASSERT( n == size ); for (auto i=0ul ; i != n ; ++i) { fillWithParser(str_arr[i], val_arr[i]); - amrex::Print() << val_arr[i] << " "; + amrex::Print() << " " << val_arr[i]; } } @@ -207,7 +230,7 @@ getWithParser (const amrex::ParmParse& pp, char const * const str, T& val) { std::vector f; pp.getarr(str, f); - amrex::Print() << str << " "; + amrex::Print() << str << " :"; fillWithParserArr(f, val); amrex::Print() << '\n'; } @@ -218,7 +241,7 @@ queryWithParser (const amrex::ParmParse& pp, char const * const str, T& val) { std::vector f; const int is_specified = pp.queryarr(str, f); - amrex::Print() << str << " "; + amrex::Print() << str << " :"; if (is_specified) { fillWithParserArr(f, val); From 931ac3e2c9005dacb19cb355558ff3b219d59190 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Wed, 15 Sep 2021 16:13:06 +0200 Subject: [PATCH 09/15] change BeamParticleContainer to ParticleTile --- src/particles/BoxSort.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/particles/BoxSort.cpp b/src/particles/BoxSort.cpp index 0e40b813aa..5eac417b9e 100644 --- a/src/particles/BoxSort.cpp +++ b/src/particles/BoxSort.cpp @@ -38,7 +38,7 @@ void BoxSorter::sortParticlesByBox (BeamParticleContainer& a_beam, amrex::Gpu::exclusive_scan(m_box_counts.begin(), m_box_counts.end(), m_box_offsets.begin()); - BeamParticleContainer tmp(a_beam.get_name()); + amrex::ParticleTile<0, 0, BeamIdx::nattribs, 0> tmp{}; tmp.resize(np); auto p_box_offsets = m_box_offsets.dataPtr(); @@ -49,7 +49,8 @@ void BoxSorter::sortParticlesByBox (BeamParticleContainer& a_beam, p_dst_indices[i] += p_box_offsets[dst_box]; }); - amrex::scatterParticles(tmp, a_beam, np, dst_indices.dataPtr()); + amrex::scatterParticles>(tmp, a_beam, np, + dst_indices.dataPtr()); a_beam.swap(tmp); } From f3d461227dcccd4b7c05f014aac2ac5beddd09a8 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Thu, 16 Sep 2021 00:56:27 +0200 Subject: [PATCH 10/15] Parse AMReX parameters --- src/Hipace.H | 20 +- src/Hipace.cpp | 30 +-- src/utils/Parser.H | 447 ++++++++++++++++++++++++++------------------- 3 files changed, 294 insertions(+), 203 deletions(-) diff --git a/src/Hipace.H b/src/Hipace.H index c50ee0130c..ce53d76bde 100644 --- a/src/Hipace.H +++ b/src/Hipace.H @@ -25,8 +25,21 @@ #include + +/** \brief Helper struct to initialize m_phys_const and Parser before amrex::AmrCore */ +struct Hipace_early_init +{ + /** Constructor of m_phys_const, initialize Hipace::m_instance, Hipace::m_normalized_units + * and Parser Constants */ + Hipace_early_init (Hipace* instance); + + /** Struct containing physical constants (which values depends on the unit system, determined + * at runtime): SI or normalized units. */ + PhysConst m_phys_const; +}; + /** \brief Singleton class, that intialize, runs and finalizes the simulation */ -class Hipace final : public amrex::AmrCore +class Hipace final : public Hipace_early_init, public amrex::AmrCore { public: /** Ctor: read general input parameters, call constructors of main member variables @@ -191,15 +204,10 @@ public: /** status of the physical time send request */ MPI_Request m_tsend_request = MPI_REQUEST_NULL; -private: /** Pointer to current (and only) instance of class Hipace */ static Hipace* m_instance; -public: /** Whether to use normalized units */ static bool m_normalized_units; - /** Struct containing physical constants (which values depends on the unit system, determined - * at runtime): SI or normalized units. */ - PhysConst m_phys_const; /** All field data (3D array, slices) and field methods */ Fields m_fields; /** Contains all beam species */ diff --git a/src/Hipace.cpp b/src/Hipace.cpp index d47d7d049b..241aeeebba 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -51,6 +51,20 @@ bool Hipace::m_do_tiling = false; bool Hipace::m_do_tiling = true; #endif +Hipace_early_init::Hipace_early_init (Hipace* instance) +{ + Hipace::m_instance = instance; + amrex::ParmParse pph("hipace"); + queryWithParser(pph ,"normalized_units", Hipace::m_normalized_units); + if (Hipace::m_normalized_units) { + m_phys_const = make_constants_normalized(); + } else { + m_phys_const = make_constants_SI(); + } + Parser::addConstantsToParser(m_phys_const); + Parser::replaceAmrexParamsWithParser(); +} + Hipace& Hipace::GetInstance () { @@ -59,16 +73,8 @@ Hipace::GetInstance () } Hipace::Hipace () : - m_phys_const([this](){ - m_instance = this; - amrex::ParmParse pph("hipace"); - pph.query("normalized_units", this->m_normalized_units); - if (this->m_normalized_units){ - return make_constants_normalized(); - } else { - return make_constants_SI(); - } - }()), + Hipace_early_init(this), + amrex::AmrCore(), m_fields(this), m_multi_beam(this), m_multi_plasma(this), @@ -84,10 +90,6 @@ Hipace::Hipace () : amrex::ParmParse pph("hipace"); - amrex::Parser parser; - auto test_func = getFunctionWithParser<3>(pph, "function", parser, {"x","y","z"}); - amrex::Print() << "function : " << test_func(1.,2.,3.) << '\n'; - queryWithParser(pph, "dt", m_dt); queryWithParser(pph, "verbose", m_verbose); queryWithParser(pph, "numprocs_x", m_numprocs_x); diff --git a/src/utils/Parser.H b/src/utils/Parser.H index bc0e6a940f..420fdcb7af 100755 --- a/src/utils/Parser.H +++ b/src/utils/Parser.H @@ -3,8 +3,6 @@ #include "Constants.H" -#include -#include #include #include @@ -17,237 +15,320 @@ #include -template -inline void -getWithParser (const amrex::ParmParse& pp, char const * const str, T& val); - template inline int queryWithParser (const amrex::ParmParse& pp, char const * const str, T& val); -inline int -safeCastToInt(const double x, const std::string& real_name) { - int result = 0; - bool error_detected = false; - std::string assert_msg; - // (2.0*(numeric_limits::max()/2+1)) converts numeric_limits::max()+1 to a real - // ensuring accuracy to all digits. This accepts x = 2**31-1 but rejects 2**31. - if (x < (2.0*(std::numeric_limits::max()/2+1))) { - if (std::ceil(x) >= std::numeric_limits::min()) { - result = static_cast(x); - } else { - error_detected = true; - assert_msg = "Error: Negative overflow detected when casting " + - real_name + " = " + std::to_string(x) + " to int"; - } - } else if (x > 0) { - error_detected = true; - assert_msg = "Error: Overflow detected when casting " + - real_name + " = " + std::to_string(x) + " to int"; - } else { - error_detected = true; - assert_msg = "Error: NaN detected when casting " + real_name + " to int"; - } - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!error_detected, assert_msg); - return result; -} - -static amrex::Parser& -initParser (amrex::Parser& parser, amrex::Vector const& varnames) +namespace Parser { - // Since queryWithParser recursively calls this routine, keep track of symbols - // in case an infinite recursion is found (a symbol's value depending on itself). - static std::set recursive_symbols{}; - - //amrex::Parser parser(parse_function); - parser.registerVariables(varnames); - - std::set symbols = parser.symbols(); - for (auto const& v : varnames) symbols.erase(v); - - // User can provide inputs under this name, through which expressions - // can be provided for arbitrary variables. This potentially - // includes variable names that match physical or mathematical - // constants, in case the user wishes to enforce a different - // system of units or some form of quasi-physical behavior in the - // simulation. Thus, this needs to override any built-in - // constants. - amrex::ParmParse pp_my_constants("my_constants"); - - // Cache for evalutated constants + // Cache for evaluated constants static std::map my_constants_cache{}; // Physical / Numerical Constants available to parsed expressions - const auto phys_const = get_phys_const(); - static std::map hipace_constants { - {"clight", phys_const.c}, - {"epsilon0", phys_const.ep0}, - {"mu0", phys_const.mu0}, - {"q_e", phys_const.q_e}, - {"m_e", phys_const.m_e}, - {"m_p", phys_const.m_p}, - {"pi", MathConst::pi} + {"pi", MathConst::pi}, + {"true", 1.}, + {"false", 0.} }; - for (auto const& s : symbols) { + /** \brief add Physical constants to Parser constants + * + * \param[in] phys_const the constants to add + */ + inline void + addConstantsToParser (const PhysConst& phys_const) { + hipace_constants.insert({"clight", phys_const.c}); + hipace_constants.insert({"epsilon0", phys_const.ep0}); + hipace_constants.insert({"mu0", phys_const.mu0}); + hipace_constants.insert({"q_e", phys_const.q_e}); + hipace_constants.insert({"m_e", phys_const.m_e}); + hipace_constants.insert({"m_p", phys_const.m_p}); + } - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(recursive_symbols.count(s)==0, - "Expressions contains recursive symbol " + s); + /** \brief replace ParmParse input with a Parsed version + * + * \param[in] pp ParmParse object + * \param[in] str name of input value to replace + */ + template + inline void + replaceWithParser (amrex::ParmParse& pp, char const * const str) { + T val{}; + if(queryWithParser(pp, str, val)) { + pp.remove(str); + pp.add(str, val); + } + } - auto usr_constant = my_constants_cache.find(s); - if (usr_constant != my_constants_cache.end()) { - parser.setConstant(s, usr_constant->second); - continue; + /** \brief array version of replaceWithParser + * + * \param[in] pp ParmParse object + * \param[in] str name of input value to replace + */ + template + inline void + replaceArrWithParser (amrex::ParmParse& pp, char const * const str) { + std::vector val{}; + if(queryWithParser(pp, str, val)) { + pp.remove(str); + pp.addarr(str, val); } + } - double v; - recursive_symbols.insert(s); - const bool is_input = queryWithParser(pp_my_constants, s.c_str(), v); - recursive_symbols.erase(s); + /** \brief replace AMReX input parameters with Parsed version + * + * AMReX only uses ParmParse when reading in these parameter, not a Parser + */ + inline void + replaceAmrexParamsWithParser () { + amrex::ParmParse pp_amr("amr"); + replaceArrWithParser(pp_amr, "n_cell"); + replaceWithParser(pp_amr, "blocking_factor"); + replaceWithParser(pp_amr, "max_level"); + amrex::ParmParse pp_geom("geometry"); + replaceWithParser(pp_geom, "coord_sys"); + replaceArrWithParser(pp_geom, "is_periodic"); + replaceArrWithParser(pp_geom, "prob_lo"); + replaceArrWithParser(pp_geom, "prob_hi"); + } - if (is_input) { - my_constants_cache.insert({s, v}); - parser.setConstant(s, v); - continue; + /** \brief return valid int, asserts if inf or NaN + * + * \param[in] x value to cast + * \param[in] real_name name of value for error message + */ + inline int + safeCastToInt (const double x, const std::string& real_name) { + int result = 0; + bool error_detected = false; + std::string assert_msg; + // (2.0*(numeric_limits::max()/2+1)) converts numeric_limits::max()+1 to a real + // ensuring accuracy to all digits. This accepts x = 2**31-1 but rejects 2**31. + if (x < (2.0*(std::numeric_limits::max()/2+1))) { + if (std::ceil(x) >= std::numeric_limits::min()) { + result = static_cast(x); + } else { + error_detected = true; + assert_msg = "Error: Negative overflow detected when casting " + + real_name + " = " + std::to_string(x) + " to int"; + } + } else if (x > 0) { + error_detected = true; + assert_msg = "Error: Overflow detected when casting " + + real_name + " = " + std::to_string(x) + " to int"; + } else { + error_detected = true; + assert_msg = "Error: NaN detected when casting " + real_name + " to int"; } + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!error_detected, assert_msg); + return result; + } - auto phy_constant = hipace_constants.find(s); - if (phy_constant != hipace_constants.end()) { - parser.setConstant(s, phy_constant->second); - continue; + /** \brief return Parser ready to compile + * + * \param[in,out] parser Parser that has been defined + * \param[in] varnames names of input variables if a function is Parsed + */ + static amrex::Parser& + initParser (amrex::Parser& parser, amrex::Vector const& varnames) { + // Since queryWithParser recursively calls this routine, keep track of symbols + // in case an infinite recursion is found (a symbol's value depending on itself). + static std::set recursive_symbols{}; + + //amrex::Parser parser(parse_function); + parser.registerVariables(varnames); + + std::set symbols = parser.symbols(); + for (auto const& v : varnames) symbols.erase(v); + + // User can provide inputs under this name, through which expressions + // can be provided for arbitrary variables. This potentially + // includes variable names that match physical or mathematical + // constants, in case the user wishes to enforce a different + // system of units or some form of quasi-physical behavior in the + // simulation. Thus, this needs to override any built-in + // constants. + amrex::ParmParse pp_my_constants("my_constants"); + + for (auto const& s : symbols) { + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(recursive_symbols.count(s)==0, + "Expressions contains recursive symbol " + s); + + auto usr_constant = my_constants_cache.find(s); + if (usr_constant != my_constants_cache.end()) { + parser.setConstant(s, usr_constant->second); + continue; + } + + double v; + recursive_symbols.insert(s); + const bool is_input = queryWithParser(pp_my_constants, s.c_str(), v); + recursive_symbols.erase(s); + + if (is_input) { + my_constants_cache.insert({s, v}); + parser.setConstant(s, v); + continue; + } + + auto phy_constant = hipace_constants.find(s); + if (phy_constant != hipace_constants.end()) { + parser.setConstant(s, phy_constant->second); + continue; + } + + amrex::Abort("makeParser::Unknown symbol " + s); } - - amrex::Abort("makeParser::Unknown symbol " + s); + return parser; } - return parser; -} -template -inline auto -getFunctionWithParser (const amrex::ParmParse& pp, - char const * const str, - amrex::Parser& parser, - amrex::Vector const& varnames) -{ - std::string func_str{}; - std::vector f; - pp.getarr(str, f); - for (auto const& s : f) { - func_str += s + ' '; + /** \brief fill second argument val with a value obtained through Parsing str + * for std::string: val is same as str without spaces + * + * \param[in] str expression in string form + * \param[out] val value parsed, can be double, float, int, bool or std::string + */ + inline void + fillWithParser (std::string const& str, double& val) { + amrex::Parser parser(str); + val = initParser(parser, {}).compileHost<0>()(); } - parser.define(func_str); - initParser(parser, varnames); - parser.print(); - return parser.compile(); -} - -inline void -fillWithParser (std::string const& str, double& val) { - amrex::Parser parser(str); - val = initParser(parser, {}).compileHost<0>()(); -} - -inline void -fillWithParser (std::string const& str, float& val) { - amrex::Parser parser(str); - val = static_cast(initParser(parser, {}).compileHost<0>()()); -} -inline void -fillWithParser (std::string const& str, int& val) { - amrex::Parser parser(str); - val = safeCastToInt(std::round(initParser(parser, {}).compileHost<0>()()),str); -} + inline void + fillWithParser (std::string const& str, float& val) { + amrex::Parser parser(str); + val = static_cast(initParser(parser, {}).compileHost<0>()()); + } -inline void -fillWithParser (std::string const& str, bool& val) { - amrex::Parser parser(str); - val = initParser(parser, {}).compileHost<0>()(); -} + inline void + fillWithParser (std::string const& str, int& val) { + amrex::Parser parser(str); + val = safeCastToInt(std::round(initParser(parser, {}).compileHost<0>()()),str); + } -inline void -fillWithParser (std::string const& str, std::string& val) { - val = str; - auto end_pos = std::remove(val.begin(), val.end(), ' '); - val.erase(end_pos, val.end()); -} + inline void + fillWithParser (std::string const& str, bool& val) { + amrex::Parser parser(str); + val = initParser(parser, {}).compileHost<0>()(); + } -template -inline void -fillWithParserArr (std::vector const& str_arr, T& val) { - std::string str{}; - for (auto const& s : str_arr) { - str += s + ' '; + inline void + fillWithParser (std::string const& str, std::string& val) { + val = str; + auto end_pos = std::remove(val.begin(), val.end(), ' '); + val.erase(end_pos, val.end()); } - fillWithParser(str, val); - amrex::Print() << " " << val; -} -template -inline void -fillWithParserArr (std::vector const& str_arr, std::vector& val_arr) { - auto const n = str_arr.size(); - val_arr.resize(n); - for (auto i=0ul ; i != n ; ++i) { - fillWithParser(str_arr[i], val_arr[i]); - amrex::Print() << " " << val_arr[i]; + /** \brief fill second argument val of array type with a value obtained through Parsing str_arr + * if val is just a single value, str_arr is reduced to a single string with spaces as + * separators + * + * \param[in] str_arr vector of expressions to be parsed + * \param[out] val value parsed, can be scalar, std::vector, amrex::Vector or std::array + */ + template + inline void + fillWithParserArr (std::vector const& str_arr, T& val) { + std::string str{}; + for (auto const& s : str_arr) { + str += s + ' '; + } + fillWithParser(str, val); } -} -template -inline void -fillWithParserArr (std::vector const& str_arr, amrex::Vector& val_arr) { - auto const n = str_arr.size(); - val_arr.resize(n); - for (auto i=0ul ; i != n ; ++i) { - fillWithParser(str_arr[i], val_arr[i]); - amrex::Print() << " " << val_arr[i]; + template + inline void + fillWithParserArr (std::vector const& str_arr, std::vector& val_arr) { + auto const n = str_arr.size(); + val_arr.resize(n); + for (auto i=0ul ; i != n ; ++i) { + fillWithParser(str_arr[i], val_arr[i]); + } } -} -template -inline void -fillWithParserArr (std::vector const& str_arr, std::array& val_arr) { - const auto n = str_arr.size(); - if (n != size) { - for( auto const& str : str_arr) { - amrex::ErrorStream() << str << ' '; + template + inline void + fillWithParserArr (std::vector const& str_arr, amrex::Vector& val_arr) { + auto const n = str_arr.size(); + val_arr.resize(n); + for (auto i=0ul ; i != n ; ++i) { + fillWithParser(str_arr[i], val_arr[i]); } - amrex::ErrorStream() << "has wrong length " << n << " should be " << size << '\n'; } - AMREX_ALWAYS_ASSERT( n == size ); - for (auto i=0ul ; i != n ; ++i) { - fillWithParser(str_arr[i], val_arr[i]); - amrex::Print() << " " << val_arr[i]; + + template + inline void + fillWithParserArr (std::vector const& str_arr, std::array& val_arr) { + const auto n = str_arr.size(); + if (n != size) { + for( auto const& str : str_arr) { + amrex::ErrorStream() << str << ' '; + } + amrex::ErrorStream() << "has wrong length " << n << " should be " << size << '\n'; + } + AMREX_ALWAYS_ASSERT( n == size ); + for (auto i=0ul ; i != n ; ++i) { + fillWithParser(str_arr[i], val_arr[i]); + } } } +/** \brief fill val with the evaluated expression from the input file + * + * \param[in] pp ParmParse that is searched for the expression + * \param[in] str name of expression + * \param[out] val value to be filled, see fillWithParserArr and fillWithParser for supported types + */ template inline void -getWithParser (const amrex::ParmParse& pp, char const * const str, T& val) -{ +getWithParser (const amrex::ParmParse& pp, char const * const str, T& val) { std::vector f; pp.getarr(str, f); - amrex::Print() << str << " :"; - fillWithParserArr(f, val); - amrex::Print() << '\n'; + Parser::fillWithParserArr(f, val); } +/** \brief return if input file contains the expression, if so it is parsed into val + * + * \param[in] pp ParmParse that is searched for the expression + * \param[in] str name of expression + * \param[out] val value to be filled, see fillWithParserArr and fillWithParser for supported types + */ template inline int -queryWithParser (const amrex::ParmParse& pp, char const * const str, T& val) -{ +queryWithParser (const amrex::ParmParse& pp, char const * const str, T& val) { std::vector f; const int is_specified = pp.queryarr(str, f); - amrex::Print() << str << " :"; - if (is_specified) - { - fillWithParserArr(f, val); + if (is_specified) { + Parser::fillWithParserArr(f, val); } - amrex::Print() << '\n'; return is_specified; } +/** \brief return function object for Host and Device from the input file + * + * \param[in] pp ParmParse that is searched for the function + * \param[in] str name of the function in pp + * \param[out] parser Parser which owns the data of the returned function + * \param[in] varnemes names of the N arguments used in the parsed function + */ +template +inline auto +getFunctionWithParser (const amrex::ParmParse& pp, + char const * const str, + amrex::Parser& parser, + amrex::Vector const& varnames) +{ + std::string func_str{}; + std::vector f; + pp.getarr(str, f); + for (auto const& s : f) { + func_str += s + ' '; + } + parser.define(func_str); + Parser::initParser(parser, varnames); + return parser.compile(); +} + #endif From 68ccf8503254dbf5698412a262d31dc2f5a65fee Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Thu, 16 Sep 2021 01:23:44 +0200 Subject: [PATCH 11/15] add Parser to some examples --- examples/blowout_wake/inputs_SI | 15 ++++++++++----- examples/blowout_wake/inputs_normalized | 4 ++-- examples/linear_wake/inputs_SI | 11 ++++++++--- 3 files changed, 20 insertions(+), 10 deletions(-) diff --git a/examples/blowout_wake/inputs_SI b/examples/blowout_wake/inputs_SI index c8a194a461..5ebe8dac58 100644 --- a/examples/blowout_wake/inputs_SI +++ b/examples/blowout_wake/inputs_SI @@ -4,6 +4,11 @@ hipace.predcorr_max_iterations = 1 hipace.predcorr_B_mixing_factor = 0.12 hipace.predcorr_B_error_tolerance = -1 +my_constants.kp_inv = 10e-6 +my_constants.kp = 1 / kp_inv +my_constants.wp = clight * kp +my_constants.ne = wp^2 * m_e * epsilon0 / q_e^2 + amr.blocking_factor = 2 amr.max_level = 0 @@ -16,9 +21,9 @@ hipace.numprocs_y = 1 hipace.depos_order_xy = 2 geometry.coord_sys = 0 # 0: Cartesian -geometry.is_periodic = 1 1 0 # Is periodic? -geometry.prob_lo = -80.e-6 -80.e-6 -60.e-6 # physical domain -geometry.prob_hi = 80.e-6 80.e-6 60.e-6 +geometry.is_periodic = true true false # Is periodic? +geometry.prob_lo = -8*kp_inv -8*kp_inv -6*kp_inv # physical domain +geometry.prob_hi = 8*kp_inv 8*kp_inv 6*kp_inv beams.names = beam beam.injection_type = fixed_ppc @@ -26,7 +31,7 @@ beam.profile = gaussian beam.zmin = -59.e-6 beam.zmax = 59.e-6 beam.radius = 12.e-6 -beam.density = 8.47187610257747e23 # 3*ne +beam.density = 3*ne beam.u_mean = 0. 0. 2000 beam.u_std = 0. 0. 0. beam.position_mean = 0. 0. 0 @@ -34,7 +39,7 @@ beam.position_std = 3.e-6 3.e-6 14.1e-6 beam.ppc = 1 1 1 plasmas.names = plasma -plasma.density = 2.8239587008591567e23 # at this density, 1/kp = 10um +plasma.density = ne plasma.ppc = 1 1 plasma.u_mean = 0.0 0.0 0. plasma.element = electron diff --git a/examples/blowout_wake/inputs_normalized b/examples/blowout_wake/inputs_normalized index a2dee695f8..34a508b533 100644 --- a/examples/blowout_wake/inputs_normalized +++ b/examples/blowout_wake/inputs_normalized @@ -17,8 +17,8 @@ hipace.numprocs_y = 1 hipace.depos_order_xy = 2 geometry.coord_sys = 0 # 0: Cartesian -geometry.is_periodic = 1 1 0 # Is periodic? -geometry.prob_lo = -8. -8. -6 # physical domain +geometry.is_periodic = true true false # Is periodic? +geometry.prob_lo = -8. -8. -6 # physical domain geometry.prob_hi = 8. 8. 6 beams.names = beam diff --git a/examples/linear_wake/inputs_SI b/examples/linear_wake/inputs_SI index cd2aabaa8e..dc5ca49b67 100644 --- a/examples/linear_wake/inputs_SI +++ b/examples/linear_wake/inputs_SI @@ -1,5 +1,10 @@ amr.n_cell = 32 32 200 +my_constants.kp_inv = 10e-6 +my_constants.kp = 1 / kp_inv +my_constants.wp = clight * kp +my_constants.ne = wp^2 * m_e * epsilon0 / q_e^2 + amr.blocking_factor = 2 amr.max_level = 0 @@ -12,7 +17,7 @@ hipace.numprocs_y = 1 hipace.depos_order_xy = 2 geometry.coord_sys = 0 # 0: Cartesian -geometry.is_periodic = 1 1 0 # Is periodic? +geometry.is_periodic = true true false # Is periodic? geometry.prob_lo = -100.e-6 -100.e-6 -75.e-6 # physical domain geometry.prob_hi = 100.e-6 100.e-6 20.e-6 @@ -22,13 +27,13 @@ beam.profile = flattop beam.zmin = -10.e-6 beam.zmax = 10.e-6 beam.radius = 30.e-6 -beam.density = 2.8239587008591567e21 +beam.density = 0.01 * ne beam.u_mean = 0. 0. 2000 beam.u_std = 0. 0. 0. beam.ppc = 1 1 1 plasmas.names = plasma -plasma.density = 2.8239587008591567e23 # at this density, 1/kp = 10um +plasma.density = ne plasma.ppc = 1 1 plasma.u_mean = 0.0 0.0 0. plasma.element = electron From a034d97796b926f2e00a6864bc8ee364b9432659 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Thu, 16 Sep 2021 01:56:26 +0200 Subject: [PATCH 12/15] fix doc --- src/utils/Parser.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/Parser.H b/src/utils/Parser.H index 420fdcb7af..38fab2e5bd 100755 --- a/src/utils/Parser.H +++ b/src/utils/Parser.H @@ -311,7 +311,7 @@ queryWithParser (const amrex::ParmParse& pp, char const * const str, T& val) { * \param[in] pp ParmParse that is searched for the function * \param[in] str name of the function in pp * \param[out] parser Parser which owns the data of the returned function - * \param[in] varnemes names of the N arguments used in the parsed function + * \param[in] varnames names of the N arguments used in the parsed function */ template inline auto From e17a7c1d35b57e0b978b9647619599a9ad953f69 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Thu, 16 Sep 2021 02:36:13 +0200 Subject: [PATCH 13/15] fix small rounding error --- .../gaussian_linear_wake.SI.1Rank.json | 26 +++++++++---------- .../benchmarks_json/linear_wake.SI.1Rank.json | 26 +++++++++---------- 2 files changed, 26 insertions(+), 26 deletions(-) diff --git a/tests/checksum/benchmarks_json/gaussian_linear_wake.SI.1Rank.json b/tests/checksum/benchmarks_json/gaussian_linear_wake.SI.1Rank.json index e6b8d0ffac..dd91b5eebe 100644 --- a/tests/checksum/benchmarks_json/gaussian_linear_wake.SI.1Rank.json +++ b/tests/checksum/benchmarks_json/gaussian_linear_wake.SI.1Rank.json @@ -6,25 +6,25 @@ "ux": 0.0, "uy": 0.0, "uz": 318304000.0, - "w": 250837608.05724, + "w": 250837609.69339, "x": 6.78895, "y": 6.78895, "z": 4.6790688 }, "lev=0": { - "Bx": 2204.9807052641, - "By": 2204.9807052646, - "Bz": 8.4174863114613, - "ExmBy": 1390604996052.0, - "EypBx": 1390604996051.9, - "Ez": 3002632583254.2, - "Psi": 40644624.184591, - "jx": 268464858918650.0, + "Bx": 2204.9807071426, + "By": 2204.9807071426, + "Bz": 8.4174862449737, + "ExmBy": 1390604994672.3, + "EypBx": 1390604994672.3, + "Ez": 3002632592545.8, + "Psi": 40644624.093825, + "jx": 268464859936910.0, "jx_beam": 0.0, - "jy": 268464858918640.0, + "jy": 268464859936900.0, "jy_beam": 0.0, - "jz": 798344094734230.0, - "jz_beam": 514058134649270.0, - "rho": 3631659.7301984 + "jz": 798344096425600.0, + "jz_beam": 514058138002350.0, + "rho": 3631659.7458766 } } \ No newline at end of file diff --git a/tests/checksum/benchmarks_json/linear_wake.SI.1Rank.json b/tests/checksum/benchmarks_json/linear_wake.SI.1Rank.json index b4ccd903ed..fc1f42c5e0 100644 --- a/tests/checksum/benchmarks_json/linear_wake.SI.1Rank.json +++ b/tests/checksum/benchmarks_json/linear_wake.SI.1Rank.json @@ -6,25 +6,25 @@ "ux": 0.0, "uy": 0.0, "uz": 6384000.0, - "w": 167253366.49385, + "w": 167253367.58481, "x": 0.041475, "y": 0.041475, "z": 0.0159201 }, "lev=0": { - "Bx": 2216.4133520235, - "By": 2216.4133520235, - "Bz": 11.07576944688, - "ExmBy": 2500407278489.8, - "EypBx": 2500407278489.8, - "Ez": 6817523944889.0, - "Psi": 67231365.551324, - "jx": 671867841626580.0, + "Bx": 2216.413355091, + "By": 2216.4133550911, + "Bz": 11.075769499901, + "ExmBy": 2500407290516.7, + "EypBx": 2500407290516.7, + "Ez": 6817523976026.3, + "Psi": 67231365.830154, + "jx": 671867845710570.0, "jx_beam": 0.0, - "jy": 671867841626580.0, + "jy": 671867845710570.0, "jy_beam": 0.0, - "jz": 1890631652716800.0, - "jz_beam": 432964477422490.0, - "rho": 6545516.9772153 + "jz": 1890631670399600.0, + "jz_beam": 432964480246620.0, + "rho": 6545517.04493 } } \ No newline at end of file From f650a228370dbcc7d23a871dcf9edcdccc5d47a6 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Tue, 21 Sep 2021 13:24:15 +0200 Subject: [PATCH 14/15] add ionization Parser example --- examples/blowout_wake/inputs_ionization_SI | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/examples/blowout_wake/inputs_ionization_SI b/examples/blowout_wake/inputs_ionization_SI index 0001cad2e4..a9884608e2 100644 --- a/examples/blowout_wake/inputs_ionization_SI +++ b/examples/blowout_wake/inputs_ionization_SI @@ -5,6 +5,11 @@ hipace.predcorr_max_iterations = 1 hipace.predcorr_B_mixing_factor = 0.12 hipace.predcorr_B_error_tolerance = -1 +my_constants.ne = 1.25e24 +my_constants.wp = sqrt(ne * q_e^2 / (epsilon0 * m_e)) +my_constants.kp = wp / clight +my_constants.kp_inv = 1 / kp # =4.753069541um + amr.blocking_factor = 2 amr.max_level = 0 @@ -24,23 +29,23 @@ geometry.prob_hi = 20.e-6 20.e-6 30.e-6 beams.names = beam beam.injection_type = fixed_ppc beam.profile = flattop -beam.zmin = 15.493860918e-6 +beam.zmin = 25e-6 - 2 * kp_inv beam.zmax = 25e-6 -beam.radius = 2.376534771e-6 -beam.density = 5e24 +beam.radius = kp_inv / 2 +beam.density = 4 * ne beam.u_mean = 0. 0. 2000 beam.u_std = 0. 0. 0. beam.ppc = 1 1 1 plasmas.names = elec ion -elec.density = 1.25e24 # at this density, 1/kp = 4.753069541um +elec.density = ne elec.ppc = 0 0 elec.u_mean = 0.0 0.0 0.0 elec.element = electron -elec.neutralize_background = 0 +elec.neutralize_background = false -ion.density = 1.25e24 # at this density, 1/kp = 4.753069541um +ion.density = ne ion.ppc = 1 1 ion.u_mean = 0.0 0.0 0.0 ion.element = H From cd1e77374efe18cc883979a832fafc28f81270a7 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Tue, 21 Sep 2021 19:28:55 +0200 Subject: [PATCH 15/15] add suggestions --- examples/blowout_wake/inputs_SI | 10 +++++----- examples/blowout_wake/inputs_ionization_SI | 8 ++++---- examples/linear_wake/inputs_SI | 4 ++-- src/utils/Parser.H | 22 ++++++++++++---------- 4 files changed, 23 insertions(+), 21 deletions(-) diff --git a/examples/blowout_wake/inputs_SI b/examples/blowout_wake/inputs_SI index 5ebe8dac58..7b01df19ee 100644 --- a/examples/blowout_wake/inputs_SI +++ b/examples/blowout_wake/inputs_SI @@ -4,8 +4,8 @@ hipace.predcorr_max_iterations = 1 hipace.predcorr_B_mixing_factor = 0.12 hipace.predcorr_B_error_tolerance = -1 -my_constants.kp_inv = 10e-6 -my_constants.kp = 1 / kp_inv +my_constants.kp_inv = 10.e-6 +my_constants.kp = 1. / kp_inv my_constants.wp = clight * kp my_constants.ne = wp^2 * m_e * epsilon0 / q_e^2 @@ -22,8 +22,8 @@ hipace.depos_order_xy = 2 geometry.coord_sys = 0 # 0: Cartesian geometry.is_periodic = true true false # Is periodic? -geometry.prob_lo = -8*kp_inv -8*kp_inv -6*kp_inv # physical domain -geometry.prob_hi = 8*kp_inv 8*kp_inv 6*kp_inv +geometry.prob_lo = -8.*kp_inv -8.*kp_inv -6.*kp_inv # physical domain +geometry.prob_hi = 8.*kp_inv 8.*kp_inv 6.*kp_inv beams.names = beam beam.injection_type = fixed_ppc @@ -31,7 +31,7 @@ beam.profile = gaussian beam.zmin = -59.e-6 beam.zmax = 59.e-6 beam.radius = 12.e-6 -beam.density = 3*ne +beam.density = 3.*ne beam.u_mean = 0. 0. 2000 beam.u_std = 0. 0. 0. beam.position_mean = 0. 0. 0 diff --git a/examples/blowout_wake/inputs_ionization_SI b/examples/blowout_wake/inputs_ionization_SI index a9884608e2..a082aaca0a 100644 --- a/examples/blowout_wake/inputs_ionization_SI +++ b/examples/blowout_wake/inputs_ionization_SI @@ -8,7 +8,7 @@ hipace.predcorr_B_error_tolerance = -1 my_constants.ne = 1.25e24 my_constants.wp = sqrt(ne * q_e^2 / (epsilon0 * m_e)) my_constants.kp = wp / clight -my_constants.kp_inv = 1 / kp # =4.753069541um +my_constants.kp_inv = 1. / kp # =4.753069541um amr.blocking_factor = 2 amr.max_level = 0 @@ -29,10 +29,10 @@ geometry.prob_hi = 20.e-6 20.e-6 30.e-6 beams.names = beam beam.injection_type = fixed_ppc beam.profile = flattop -beam.zmin = 25e-6 - 2 * kp_inv -beam.zmax = 25e-6 +beam.zmin = 25.e-6 - 2. * kp_inv +beam.zmax = 25.e-6 beam.radius = kp_inv / 2 -beam.density = 4 * ne +beam.density = 4. * ne beam.u_mean = 0. 0. 2000 beam.u_std = 0. 0. 0. beam.ppc = 1 1 1 diff --git a/examples/linear_wake/inputs_SI b/examples/linear_wake/inputs_SI index dc5ca49b67..1d2814ade9 100644 --- a/examples/linear_wake/inputs_SI +++ b/examples/linear_wake/inputs_SI @@ -1,7 +1,7 @@ amr.n_cell = 32 32 200 -my_constants.kp_inv = 10e-6 -my_constants.kp = 1 / kp_inv +my_constants.kp_inv = 10.e-6 +my_constants.kp = 1. / kp_inv my_constants.wp = clight * kp my_constants.ne = wp^2 * m_e * epsilon0 / q_e^2 diff --git a/src/utils/Parser.H b/src/utils/Parser.H index 38fab2e5bd..60d911f370 100755 --- a/src/utils/Parser.H +++ b/src/utils/Parser.H @@ -28,8 +28,8 @@ namespace Parser static std::map hipace_constants { {"pi", MathConst::pi}, - {"true", 1.}, - {"false", 0.} + {"true", 1}, + {"false", 0} }; /** \brief add Physical constants to Parser constants @@ -125,18 +125,17 @@ namespace Parser return result; } - /** \brief return Parser ready to compile + /** \brief init Parser ready to compile * * \param[in,out] parser Parser that has been defined * \param[in] varnames names of input variables if a function is Parsed */ - static amrex::Parser& + static void initParser (amrex::Parser& parser, amrex::Vector const& varnames) { // Since queryWithParser recursively calls this routine, keep track of symbols // in case an infinite recursion is found (a symbol's value depending on itself). static std::set recursive_symbols{}; - //amrex::Parser parser(parse_function); parser.registerVariables(varnames); std::set symbols = parser.symbols(); @@ -181,7 +180,6 @@ namespace Parser amrex::Abort("makeParser::Unknown symbol " + s); } - return parser; } /** \brief fill second argument val with a value obtained through Parsing str @@ -193,25 +191,29 @@ namespace Parser inline void fillWithParser (std::string const& str, double& val) { amrex::Parser parser(str); - val = initParser(parser, {}).compileHost<0>()(); + initParser(parser, {}); + val = parser.compileHost<0>()(); } inline void fillWithParser (std::string const& str, float& val) { amrex::Parser parser(str); - val = static_cast(initParser(parser, {}).compileHost<0>()()); + initParser(parser, {}); + val = static_cast(parser.compileHost<0>()()); } inline void fillWithParser (std::string const& str, int& val) { amrex::Parser parser(str); - val = safeCastToInt(std::round(initParser(parser, {}).compileHost<0>()()),str); + initParser(parser, {}); + val = safeCastToInt(std::round(parser.compileHost<0>()()),str); } inline void fillWithParser (std::string const& str, bool& val) { amrex::Parser parser(str); - val = initParser(parser, {}).compileHost<0>()(); + initParser(parser, {}); + val = parser.compileHost<0>()(); } inline void