More parser-enabled ParmParse.query and ParmParse.get (#1501)

* Parser can be used for most input parameters, including charge and mass

* update documentation for the math parser

* eol

* remove typo and update doc for parser

* Apply suggestions from code review

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

8 files changed, 104 insertions, 112 deletions

diff --git a/Source/Diagnostics/BTDiagnostics.cpp b/Source/Diagnostics/BTDiagnostics.cpp
index b651874e4..2f6a2dfbc 100644
--- a/Source/Diagnostics/BTDiagnostics.cpp
+++ b/Source/Diagnostics/BTDiagnostics.cpp
@@ -109,8 +109,8 @@ BTDiagnostics::ReadParameters ()
 
     // Read either dz_snapshots_lab or dt_snapshots_lab
     bool snapshot_interval_is_specified = false;
-    snapshot_interval_is_specified = pp.query("dt_snapshots_lab", m_dt_snapshots_lab);
-    if ( pp.query("dz_snapshots_lab", m_dz_snapshots_lab) ) {
+    snapshot_interval_is_specified = queryWithParser(pp, "dt_snapshots_lab", m_dt_snapshots_lab);
+    if ( queryWithParser(pp, "dz_snapshots_lab", m_dz_snapshots_lab) ) {
         m_dt_snapshots_lab = m_dz_snapshots_lab/PhysConst::c;
         snapshot_interval_is_specified = true;
     }
diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp
index d4c9e9ba5..6c9f66bc0 100644
--- a/Source/Initialization/PlasmaInjector.cpp
+++ b/Source/Initialization/PlasmaInjector.cpp
@@ -22,7 +22,6 @@
 #include <string>
 #include <memory>
 
-
 using namespace amrex;
 
 namespace {
@@ -34,40 +33,6 @@ namespace {
         string = stringstream.str();
         amrex::Abort(string.c_str());
     }
-
-    Real parseChargeName(const ParmParse& pp, const std::string& name) {
-        Real result;
-        if (name == "q_e") {
-            return PhysConst::q_e;
-        } else if (pp.query("charge", result)) {
-            return result;
-        } else {
-            StringParseAbortMessage("Charge", name);
-            return 0.0;
-        }
-    }
-
-    Real parseChargeString(const ParmParse& pp, const std::string& name) {
-        if(name.substr(0, 1) == "-")
-            return -1.0 * parseChargeName(pp, name.substr(1, name.size() - 1));
-        return parseChargeName(pp, name);
-    }
-
-    Real parseMassString(const ParmParse& pp, const std::string& name) {
-        Real result;
-        if (name == "m_e") {
-            return PhysConst::m_e;
-        } else if (name == "m_p"){
-            return PhysConst::m_p;
-        } else if (name == "inf"){
-            return std::numeric_limits<double>::infinity();
-        } else if (pp.query("mass", result)) {
-            return result;
-        } else {
-            StringParseAbortMessage("Mass", name);
-            return 0.0;
-        }
-    }
 }
 
 PlasmaInjector::PlasmaInjector () {}
@@ -136,9 +101,7 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name)
     bool species_is_specified = pp.query("species_type", physical_species_s);
     if (species_is_specified){
         physical_species = species::from_string( physical_species_s );
-        // charge = SpeciesCharge[physical_species];
         charge = species::get_charge( physical_species );
-        // mass = SpeciesMass[physical_species];
         mass = species::get_mass( physical_species );
     }
 
@@ -146,18 +109,9 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name)
     pp.query("injection_style", s_inj_style);
 
     // parse charge and mass
-    std::string charge_s;
-    std::string mass_s;
-    bool charge_is_specified = pp.query("charge", charge_s);
-    bool mass_is_specified = pp.query("mass", mass_s);
-
-    if (charge_is_specified){
-        std::transform(charge_s.begin(),
-                       charge_s.end(),
-                       charge_s.begin(),
-                       ::tolower);
-        charge = parseChargeString(pp, charge_s);
-    }
+    bool charge_is_specified = queryWithParser(pp, "charge", charge);
+    bool mass_is_specified = queryWithParser(pp, "mass", mass);
+
     if ( charge_is_specified && species_is_specified ){
         Print() << "WARNING: Both '" << species_name << ".charge' and "
                 << species_name << ".species_type' are specified\n'"
@@ -168,13 +122,6 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name)
         amrex::Abort("Need to specify at least one of species_type or charge");
     }
 
-    if (mass_is_specified){
-        std::transform(mass_s.begin(),
-                       mass_s.end(),
-                       mass_s.begin(),
-                       ::tolower);
-        mass = parseMassString(pp, mass_s);
-    }
     if ( mass_is_specified && species_is_specified ){
         Print() << "WARNING: Both '" << species_name << ".mass' and "
                 << species_name << ".species_type' are specified\n'"
@@ -205,16 +152,16 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name)
         add_single_particle = true;
         return;
     } else if (part_pos_s == "gaussian_beam") {
-        pp.get("x_m", x_m);
-        pp.get("y_m", y_m);
-        pp.get("z_m", z_m);
-        pp.get("x_rms", x_rms);
-        pp.get("y_rms", y_rms);
-        pp.get("z_rms", z_rms);
-        pp.query("x_cut", x_cut);
-        pp.query("y_cut", y_cut);
-        pp.query("z_cut", z_cut);
-        pp.get("q_tot", q_tot);
+        getWithParser(pp, "x_m", x_m);
+        getWithParser(pp, "y_m", y_m);
+        getWithParser(pp, "z_m", z_m);
+        getWithParser(pp, "x_rms", x_rms);
+        getWithParser(pp, "y_rms", y_rms);
+        getWithParser(pp, "z_rms", z_rms);
+        queryWithParser(pp, "x_cut", x_cut);
+        queryWithParser(pp, "y_cut", y_cut);
+        queryWithParser(pp, "z_cut", z_cut);
+        getWithParser(pp, "q_tot", q_tot);
         pp.get("npart", npart);
         pp.query("do_symmetrize", do_symmetrize);
         gaussian_beam = true;
@@ -274,8 +221,8 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name)
         std::string str_injection_file;
         pp.get("injection_file", str_injection_file);
         // optional parameters
-        pp.query("q_tot", q_tot);
-        pp.query("z_shift",z_shift);
+        queryWithParser(pp, "q_tot", q_tot);
+        queryWithParser(pp, "z_shift",z_shift);
 
 #ifdef WARPX_USE_OPENPMD
         if (ParallelDescriptor::IOProcessor()) {
@@ -454,9 +401,9 @@ void PlasmaInjector::parseMomentum (ParmParse& pp)
         Real ux = 0.;
         Real uy = 0.;
         Real uz = 0.;
-        pp.query("ux", ux);
-        pp.query("uy", uy);
-        pp.query("uz", uz);
+        queryWithParser(pp, "ux", ux);
+        queryWithParser(pp, "uy", uy);
+        queryWithParser(pp, "uz", uz);
         // Construct InjectorMomentum with InjectorMomentumConstant.
         h_inj_mom.reset(new InjectorMomentum((InjectorMomentumConstant*)nullptr, ux,uy, uz));
     } else if (mom_dist_s == "custom") {
@@ -469,12 +416,12 @@ void PlasmaInjector::parseMomentum (ParmParse& pp)
         Real ux_th = 0.;
         Real uy_th = 0.;
         Real uz_th = 0.;
-        pp.query("ux_m", ux_m);
-        pp.query("uy_m", uy_m);
-        pp.query("uz_m", uz_m);
-        pp.query("ux_th", ux_th);
-        pp.query("uy_th", uy_th);
-        pp.query("uz_th", uz_th);
+        queryWithParser(pp, "ux_m", ux_m);
+        queryWithParser(pp, "uy_m", uy_m);
+        queryWithParser(pp, "uz_m", uz_m);
+        queryWithParser(pp, "ux_th", ux_th);
+        queryWithParser(pp, "uy_th", uy_th);
+        queryWithParser(pp, "uz_th", uz_th);
         // Construct InjectorMomentum with InjectorMomentumGaussian.
         h_inj_mom.reset(new InjectorMomentum((InjectorMomentumGaussian*)nullptr,
                                              ux_m, uy_m, uz_m, ux_th, uy_th, uz_th));
@@ -483,11 +430,11 @@ void PlasmaInjector::parseMomentum (ParmParse& pp)
         Real theta = 10.;
         int dir = 0;
         std::string direction = "x";
-        pp.query("beta", beta);
+        queryWithParser(pp, "beta", beta);
         if(beta < 0){
             amrex::Abort("Please enter a positive beta value. Drift direction is set with <s_name>.bulk_vel_dir = 'x' or '+x', '-x', 'y' or '+y', etc.");
         }
-        pp.query("theta", theta);
+        queryWithParser(pp, "theta", theta);
         pp.query("bulk_vel_dir", direction);
         if(direction[0] == '-'){
             beta = -beta;
@@ -514,11 +461,11 @@ void PlasmaInjector::parseMomentum (ParmParse& pp)
         Real theta = 10.;
         int dir = 0;
         std::string direction = "x";
-        pp.query("beta", beta);
+        queryWithParser(pp, "beta", beta);
         if(beta < 0){
             amrex::Abort("Please enter a positive beta value. Drift direction is set with <s_name>.bulk_vel_dir = 'x' or '+x', '-x', 'y' or '+y', etc.");
         }
-        pp.query("theta", theta);
+        queryWithParser(pp, "theta", theta);
         pp.query("bulk_vel_dir", direction);
         if(direction[0] == '-'){
             beta = -beta;
@@ -542,7 +489,7 @@ void PlasmaInjector::parseMomentum (ParmParse& pp)
         h_inj_mom.reset(new InjectorMomentum((InjectorMomentumJuttner*)nullptr, theta, beta, dir));
     } else if (mom_dist_s == "radial_expansion") {
         Real u_over_r = 0.;
-        pp.query("u_over_r", u_over_r);
+        queryWithParser(pp, "u_over_r", u_over_r);
         // Construct InjectorMomentum with InjectorMomentumRadialExpansion.
         h_inj_mom.reset(new InjectorMomentum
                         ((InjectorMomentumRadialExpansion*)nullptr, u_over_r));
diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp
index 05c6c6f29..8f74a9252 100644
--- a/Source/Laser/LaserParticleContainer.cpp
+++ b/Source/Laser/LaserParticleContainer.cpp
@@ -51,8 +51,8 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies,
     pp.getarr("polarization", m_p_X);
 
     pp.query("pusher_algo", m_pusher_algo);
-    pp.get("wavelength", m_wavelength);
-    pp.get("e_max", m_e_max);
+    getWithParser(pp, "wavelength", m_wavelength);
+    getWithParser(pp, "e_max", m_e_max);
     pp.query("do_continuous_injection", do_continuous_injection);
     pp.query("min_particles_per_mode", m_min_particles_per_mode);
 
diff --git a/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp b/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp
index e7dc795a5..b795174d2 100644
--- a/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp
+++ b/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp
@@ -8,6 +8,7 @@
 #include "Laser/LaserProfiles.H"
 #include "Utils/WarpX_Complex.H"
 #include "Utils/WarpXConst.H"
+#include "Utils/WarpXUtil.H"
 
 #include <cmath>
 
@@ -24,14 +25,14 @@ WarpXLaserProfiles::GaussianLaserProfile::init (
     m_common_params = params;
 
     // Parse the properties of the Gaussian profile
-    ppl.get("profile_waist", m_params.waist);
-    ppl.get("profile_duration", m_params.duration);
-    ppl.get("profile_t_peak", m_params.t_peak);
-    ppl.get("profile_focal_distance", m_params.focal_distance);
-    ppl.query("zeta", m_params.zeta);
-    ppl.query("beta", m_params.beta);
-    ppl.query("phi2", m_params.phi2);
-    ppl.query("phi0", m_params.phi0);
+    getWithParser(ppl, "profile_waist", m_params.waist);
+    getWithParser(ppl, "profile_duration", m_params.duration);
+    getWithParser(ppl, "profile_t_peak", m_params.t_peak);
+    getWithParser(ppl, "profile_focal_distance", m_params.focal_distance);
+    queryWithParser(ppl, "zeta", m_params.zeta);
+    queryWithParser(ppl, "beta", m_params.beta);
+    queryWithParser(ppl, "phi2", m_params.phi2);
+    queryWithParser(ppl, "phi0", m_params.phi0);
 
     m_params.stc_direction = m_common_params.p_X;
     ppl.queryarr("stc_direction", m_params.stc_direction);
diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp
index b27b67dc9..8e91b12bb 100644
--- a/Source/Particles/RigidInjectedParticleContainer.cpp
+++ b/Source/Particles/RigidInjectedParticleContainer.cpp
@@ -37,7 +37,7 @@ RigidInjectedParticleContainer::RigidInjectedParticleContainer (AmrCore* amr_cor
 
     ParmParse pp(species_name);
 
-    pp.get("zinject_plane", zinject_plane);
+    getWithParser(pp, "zinject_plane", zinject_plane);
     pp.query("projected", projected);
     pp.query("focused", focused);
     pp.query("rigid_advance", rigid_advance);
diff --git a/Source/Utils/WarpXUtil.H b/Source/Utils/WarpXUtil.H
index e5a01a44a..d72ed683b 100644
--- a/Source/Utils/WarpXUtil.H
+++ b/Source/Utils/WarpXUtil.H
@@ -42,7 +42,7 @@ void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin,
  * \param query_string ParmParse.query will look for this string
  * \param stored_string variable in which the string to parse is stored
  */
-void Store_parserString(amrex::ParmParse &pp, std::string query_string,
+void Store_parserString(const amrex::ParmParse &pp, std::string query_string,
                         std::string& stored_string);
 
 namespace WarpXUtilIO{
@@ -173,7 +173,21 @@ WarpXParser makeParser (std::string const& parse_function, std::vector<std::stri
  * \param[in] str name of the parameter to read
  * \param[out] val where the value queried and parsed is stored
  */
-int queryWithParser (amrex::ParmParse& a_pp, char const * const str, amrex::Real& val);
+int queryWithParser (const amrex::ParmParse& a_pp, char const * const str, amrex::Real& val);
+
+/**
+ * \brief Similar to amrex::ParmParse::get, but also supports math expressions for the value.
+ *
+ * amrex::ParmParse::get reads a name and a value from the input file. This function does the
+ * same, and applies the WarpXParser to the value, so the user has the choice to specify a value or
+ * a math expression (including user-defined constants).
+ * Only works for amrex::Real numbers, one would need another version for integers etc.
+ *
+ * \param[in] a_pp amrex::ParmParse object
+ * \param[in] str name of the parameter to read
+ * \param[out] val where the value queried and parsed is stored
+ */
+void getWithParser (const amrex::ParmParse& a_pp, char const * const str, amrex::Real& val);
 
 namespace WarpXUtilMsg{
 
diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp
index 4372c132e..7b66b416e 100644
--- a/Source/Utils/WarpXUtil.cpp
+++ b/Source/Utils/WarpXUtil.cpp
@@ -21,7 +21,7 @@ void ReadBoostedFrameParameters(Real& gamma_boost, Real& beta_boost,
                                 Vector<int>& boost_direction)
 {
     ParmParse pp("warpx");
-    pp.query("gamma_boost", gamma_boost);
+    queryWithParser(pp, "gamma_boost", gamma_boost);
     if( gamma_boost > 1. ) {
         beta_boost = std::sqrt(1.-1./pow(gamma_boost,2));
         std::string s;
@@ -178,7 +178,7 @@ namespace WarpXUtilIO{
     }
 }
 
-void Store_parserString(amrex::ParmParse& pp, std::string query_string,
+void Store_parserString(const amrex::ParmParse& pp, std::string query_string,
                         std::string& stored_string)
 {
     std::vector<std::string> f;
@@ -200,10 +200,28 @@ WarpXParser makeParser (std::string const& parse_function, std::vector<std::stri
     for (auto it = symbols.begin(); it != symbols.end(); ) {
         Real v;
         if (pp.query(it->c_str(), v)) {
-           parser.setConstant(*it, v);
-           it = symbols.erase(it);
+            parser.setConstant(*it, v);
+            it = symbols.erase(it);
+        } else if (std::strcmp(it->c_str(), "q_e") == 0) {
+            parser.setConstant(*it, PhysConst::q_e);
+            it = symbols.erase(it);
+        } else if (std::strcmp(it->c_str(), "m_e") == 0) {
+            parser.setConstant(*it, PhysConst::m_e);
+            it = symbols.erase(it);
+        } else if (std::strcmp(it->c_str(), "m_p") == 0) {
+            parser.setConstant(*it, PhysConst::m_p);
+            it = symbols.erase(it);
+        } else if (std::strcmp(it->c_str(), "epsilon0") == 0) {
+            parser.setConstant(*it, PhysConst::ep0);
+            it = symbols.erase(it);
+        } else if (std::strcmp(it->c_str(), "clight") == 0) {
+            parser.setConstant(*it, PhysConst::c);
+            it = symbols.erase(it);
+        } else if (std::strcmp(it->c_str(), "pi") == 0) {
+            parser.setConstant(*it, MathConst::pi);
+            it = symbols.erase(it);
         } else {
-           ++it;
+            ++it;
         }
     }
     for (auto const& s : symbols) {
@@ -213,7 +231,7 @@ WarpXParser makeParser (std::string const& parse_function, std::vector<std::stri
 }
 
 int
-queryWithParser (amrex::ParmParse& a_pp, char const * const str, amrex::Real& val)
+queryWithParser (const amrex::ParmParse& a_pp, char const * const str, amrex::Real& val)
 {
     // call amrex::ParmParse::query, check if the user specified str.
     std::string tmp_str;
@@ -231,6 +249,17 @@ queryWithParser (amrex::ParmParse& a_pp, char const * const str, amrex::Real& va
     return is_specified;
 }
 
+void
+getWithParser (const amrex::ParmParse& a_pp, char const * const str, amrex::Real& val)
+{
+    // If so, create a parser object and apply it to the value provided by the user.
+    std::string str_val;
+    Store_parserString(a_pp, str, str_val);
+
+    auto parser = makeParser(str_val, {});
+    val = parser.eval();
+}
+
 /**
  * \brief Ensures that the blocks are setup correctly for the RZ spectral solver
  * When using the RZ spectral solver, the Hankel transform cannot be
@@ -345,4 +374,3 @@ namespace WarpXUtilStr
     }
 
 }
-
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 99d360291..a8afa17ec 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -321,7 +321,7 @@ WarpX::ReadParameters ()
     {
         ParmParse pp;// Traditionally, max_step and stop_time do not have prefix.
         pp.query("max_step", max_step);
-        pp.query("stop_time", stop_time);
+        queryWithParser(pp, "stop_time", stop_time);
         pp.query("authors", authors);
     }
 
@@ -368,7 +368,7 @@ WarpX::ReadParameters ()
             }
         }
 
-        pp.query("cfl", cfl);
+        queryWithParser(pp, "cfl", cfl);
         pp.query("verbose", verbose);
         pp.query("regrid_int", regrid_int);
         pp.query("do_subcycling", do_subcycling);
@@ -388,7 +388,7 @@ WarpX::ReadParameters ()
         // pp.query returns 1 if argument zmax_plasma_to_compute_max_step is
         // specified by the user, 0 otherwise.
         do_compute_max_step_from_zmax =
-            pp.query("zmax_plasma_to_compute_max_step",
+            queryWithParser(pp, "zmax_plasma_to_compute_max_step",
                       zmax_plasma_to_compute_max_step);
 
         pp.query("do_moving_window", do_moving_window);
@@ -439,8 +439,8 @@ WarpX::ReadParameters ()
             // Read either dz_snapshots_lab or dt_snapshots_lab
             bool snapshot_interval_is_specified = 0;
             Real dz_snapshots_lab = 0;
-            snapshot_interval_is_specified += pp.query("dt_snapshots_lab", dt_snapshots_lab);
-            if ( pp.query("dz_snapshots_lab", dz_snapshots_lab) ){
+            snapshot_interval_is_specified += queryWithParser(pp, "dt_snapshots_lab", dt_snapshots_lab);
+            if ( queryWithParser(pp, "dz_snapshots_lab", dz_snapshots_lab) ){
                 dt_snapshots_lab = dz_snapshots_lab/PhysConst::c;
                 snapshot_interval_is_specified = 1;
             }
@@ -520,10 +520,12 @@ WarpX::ReadParameters ()
                 sort_bin_size[i] = vect_sort_bin_size[i];
         }
 
-        double quantum_xi;
-        int quantum_xi_is_specified = pp.query("quantum_xi", quantum_xi);
-        if (quantum_xi_is_specified)
+        amrex::Real quantum_xi_tmp;
+        int quantum_xi_is_specified = queryWithParser(pp, "quantum_xi", quantum_xi_tmp);
+        if (quantum_xi_is_specified) {
+            double const quantum_xi = quantum_xi_tmp;
             quantum_xi_c2 = quantum_xi * PhysConst::c * PhysConst::c;
+        }
 
         pp.query("do_pml", do_pml);
         pp.query("pml_ncell", pml_ncell);