Merge pull request #335 from JamesMcClung/pr/bgk-variable-input-cols

BGK: generalize initial conditions parser

src/psc_bgk.cxx

@@ -8,7 +8,7 @@
 #include "psc_config.hxx"
 
 #include "psc_bgk_util/bgk_params.hxx"
-#include "psc_bgk_util/input_parser.hxx"
+#include "psc_bgk_util/table.hxx"
 #include "psc_bgk_util/params_parser.hxx"
 
 // ======================================================================
@@ -37,24 +37,13 @@ using Checks = PscConfig::Checks;
 using Marder = PscConfig::Marder;
 using OutputParticles = PscConfig::OutputParticles;
 
-// for parser
-enum DATA_COL
-{
-  COL_RHO,
-  COL_NE,
-  COL_V_PHI,
-  COL_TE,
-  COL_E_RHO,
-  COL_PHI,
-  n_cols
-};
-
 // ======================================================================
 // Global parameters
 
 namespace
 {
-ParsedData* parsedData;
+// table of initial conditions
+Table* ic_table;
 
 std::string read_checkpoint_filename;
 
@@ -78,9 +67,8 @@ void setupParameters(int argc, char** argv)
   }
   std::string path_to_params(argv[1]);
   ParsedParams parsedParams(path_to_params);
-  parsedData = new ParsedData(parsedParams.get<std::string>("path_to_data"),
-                              n_cols, COL_RHO, 1);
-  g.loadParams(parsedParams, *parsedData);
+  ic_table = new Table(parsedParams.get<std::string>("path_to_data"));
+  g.loadParams(parsedParams, *ic_table);
 
   psc_params.nmax = parsedParams.get<int>("nmax");
   psc_params.stats_every = parsedParams.get<int>("stats_every");
@@ -204,7 +192,7 @@ inline double getIonDensity(double rho)
 {
   if (!g.do_ion)
     return g.n_i;
-  double potential = parsedData->get_interpolated(COL_PHI, rho);
+  double potential = ic_table->get_interpolated("Psi", "rho", rho);
   return std::exp(-potential / g.T_i);
 }
 
@@ -360,7 +348,7 @@ void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts,
           // other's moments are given in input file
 
           double psi_cs =
-            parsedData->get_interpolated(COL_PHI, rho) / sqr(g.beta);
+            ic_table->get_interpolated("Psi", "rho", rho) / sqr(g.beta);
           double n_background = exp(psi_cs);
           double p_background = n_background / np.n;
           np.p = pdist_case4(p_background, y, z, rho);
@@ -369,12 +357,12 @@ void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts,
         }
 
         if (rho == 0) {
-          double Te = parsedData->get_interpolated(COL_TE, rho);
+          double Te = ic_table->get_interpolated("Te", "rho", rho);
           np.p = setup_particles.createMaxwellian(
             {np.kind, np.n, {0, 0, 0}, {Te, Te, Te}, np.tag});
         } else if (g.maxwellian) {
-          double Te = parsedData->get_interpolated(COL_TE, rho);
-          double vphi = parsedData->get_interpolated(COL_V_PHI, rho);
+          double Te = ic_table->get_interpolated("Te", "rho", rho);
+          double vphi = ic_table->get_interpolated("v_phi", "rho", rho);
           double coef = g.v_e_coef * (g.reverse_v ? -1 : 1) *
                         (g.reverse_v_half && y < 0 ? -1 : 1);
 
@@ -419,7 +407,7 @@ void initializePhi(BgkMfields& phi)
   setupScalarField(
     phi, Centering::Centerer(Centering::NC), [&](int m, double crd[3]) {
       double rho = sqrt(sqr(getCoord(crd[1])) + sqr(getCoord(crd[2])));
-      return parsedData->get_interpolated(COL_PHI, rho);
+      return ic_table->get_interpolated("Psi", "rho", rho);
     });
 
   writeMF(phi, "phi", {"phi"});
@@ -556,7 +544,7 @@ static void run(int argc, char** argv)
     read_checkpoint(read_checkpoint_filename, *grid_ptr, mprts, mflds);
   }
 
-  delete parsedData;
+  delete ic_table;
 
   // ----------------------------------------------------------------------
   // Hand off to PscIntegrator to run the simulation

src/psc_bgk_util/bgk_params.hxx

@@ -2,7 +2,7 @@
 
 #include <cmath>
 #include "params_parser.hxx"
-#include "input_parser.hxx"
+#include "table.hxx"
 
 // ======================================================================
 // getCalculatedBoxSize
@@ -36,23 +36,23 @@ double getSpikeSize(double B, double k, double beta)
 
 // The hole size is determined empirically from the input data. This function
 // finds where |electron density - 1| <= epsilon, where epsilon is small.
-double getHoleSize(ParsedData& data)
+double getHoleSize(Table& ic_table)
 {
   double epsilon = 1e-4;
   const int COL_RHO = 0;
   const int COL_NE = 1;
-  for (int row = data.get_nrows() - 1; row >= 0; row--) {
-    if (std::abs(data[row][COL_NE] - 1) > epsilon)
-      return data[row][COL_RHO];
+  for (int row = ic_table.n_rows() - 1; row >= 0; row--) {
+    if (std::abs(ic_table.get("ne", row) - 1) > epsilon)
+      return ic_table.get("rho", row);
   }
   throw "Unable to determine hole size.";
 }
 
 // Calculates a box size big enough to resolve the spike and the hole.
-double getCalculatedBoxSize(double B, double k, double beta, ParsedData& data)
+double getCalculatedBoxSize(double B, double k, double beta, Table& ic_table)
 {
   double spike_size = getSpikeSize(B, k, beta);
-  double hole_size = getHoleSize(data);
+  double hole_size = getHoleSize(ic_table);
   LOG_INFO("spike radius: %f\thole radius: %f\n", spike_size, hole_size);
   return 2 * std::max(spike_size, hole_size);
 }
@@ -106,7 +106,7 @@ struct PscBgkParams
   double rel_box_size_3; // length of 3rd dimension in calculated units
   int n_patches_3;       // number of patches in 3rd dimension
 
-  void loadParams(ParsedParams parsedParams, ParsedData& data)
+  void loadParams(ParsedParams parsedParams, Table& ic_table)
   {
     box_size = parsedParams.getAndWarnOrDefault<double>("box_size", -1);
     rel_box_size = parsedParams.getOrDefault<double>("rel_box_size", 1);
@@ -162,8 +162,8 @@ struct PscBgkParams
     }
 
     if (box_size <= 0)
-      box_size = rel_box_size * getCalculatedBoxSize(Hx, k, beta, data);
+      box_size = rel_box_size * getCalculatedBoxSize(Hx, k, beta, ic_table);
     if (box_size_3 <= 0)
-      box_size_3 = rel_box_size_3 * getCalculatedBoxSize(Hx, k, beta, data);
+      box_size_3 = rel_box_size_3 * getCalculatedBoxSize(Hx, k, beta, ic_table);
   }
 };