Merge pull request #328 from JamesMcClung/pr/bgk-configure-h0

BGK: configure h0, k
psc-code · Apr 3, 2024 · 62a4ece · 62a4ece
2 parents 875b7eb + 9470f88
commit 62a4ece
Show file tree

Hide file tree

Showing 3 changed files with 26 additions and 22 deletions.
diff --git a/src/psc_bgk.cxx b/src/psc_bgk.cxx
@@ -208,28 +208,15 @@ inline double getIonDensity(double rho)
   return std::exp(-potential / g.T_i);
 }
 
-// ======================================================================
-// get_beta
-// Return the conversion factor from paper units to psc units.
-
-double get_beta()
-{
-  // PSC is normalized as c=1, but the paper has electron thermal velocity
-  // v_e=1. Beta is v_e/c = sqrt(Te_paper) / sqrt(Te_psc)
-  const double PAPER_ELECTRON_TEMPERATURE = 1.;
-  const double pscElectronTemperature = parsedData->get_interpolated(COL_TE, 0);
-  return std::sqrt(pscElectronTemperature / PAPER_ELECTRON_TEMPERATURE);
-}
-
 // ======================================================================
 // v_phi_cdf
 // Cumulative distribution function for azimuthal electron velocity
 
 double v_phi_cdf(double v_phi, double rho)
 {
   // convert units from psc to paper
-  v_phi /= get_beta();
-  rho /= get_beta();
+  v_phi /= get_beta(*parsedData);
+  rho /= get_beta(*parsedData);
 
   double B = g.Hx;
   double sqrt2 = std::sqrt(2);
@@ -259,8 +246,8 @@ struct pdist
       z{z},
       rho{rho},
       v_phi_dist{[=](double v_phi) { return v_phi_cdf(v_phi, rho); }},
-      v_rho_dist{0, get_beta()},
-      v_x_dist{0, get_beta()}
+      v_rho_dist{0, get_beta(*parsedData)},
+      v_x_dist{0, get_beta(*parsedData)}
   {}
 
   Double3 operator()()

diff --git a/src/psc_bgk_util/bgk_params.hxx b/src/psc_bgk_util/bgk_params.hxx
@@ -4,6 +4,20 @@
 #include "params_parser.hxx"
 #include "input_parser.hxx"
 
+// ======================================================================
+// get_beta
+// Return the conversion factor from paper units to psc units.
+
+double get_beta(ParsedData& data)
+{
+  // PSC is normalized as c=1, but the paper has electron thermal velocity
+  // v_e=1. Beta is v_e/c = sqrt(Te_paper) / sqrt(Te_psc)
+  const double PAPER_ELECTRON_TEMPERATURE = 1.;
+  const int COL_TE = 3;
+  const double pscElectronTemperature = data.get_interpolated(COL_TE, 0);
+  return std::sqrt(pscElectronTemperature / PAPER_ELECTRON_TEMPERATURE);
+}
+
 // ======================================================================
 // getCalculatedBoxSize
 
@@ -17,7 +31,7 @@
 // The negative Gaussian has a mean given by the RHS of the equation. It drifts
 // up, approaching a line with slope B/2. The negative Gaussian is the source
 // of the spike.
-double getSpikeSize(double B, double k)
+double getSpikeSize(double B, double k, double beta)
 {
   double max_v = 4.5;
   // solve cubic with linear coefficient = 0
@@ -30,7 +44,6 @@ double getSpikeSize(double B, double k)
   double q = p * p * p + t;
   double s = sqrt(t * (2. * q - t));
 
-  double beta = .001; // FIXME don't hardcode this (see psc_bgk.cxx get_beta())
   double r = (std::cbrt(q + s) + std::cbrt(q - s) + p) * beta;
   return r;
 }
@@ -52,7 +65,7 @@ double getHoleSize(ParsedData& data)
 // Calculates a box size big enough to resolve the spike and the hole.
 double getCalculatedBoxSize(double B, double k, ParsedData& data)
 {
-  double spike_size = getSpikeSize(B, k);
+  double spike_size = getSpikeSize(B, k, get_beta(data));
   double hole_size = getHoleSize(data);
   LOG_INFO("spike radius: %f\thole radius: %f\n", spike_size, hole_size);
   return 2 * std::max(spike_size, hole_size);
@@ -76,8 +89,8 @@ struct PscBgkParams
   int n_patches;       // number of patches
   int nicell;          // number of particles per gripdoint when density=1
 
-  double k = .1;  // a parameter for BGK solutions
-  double h0 = .9; // a parameter for BGK solutions
+  double k;  // a parameter for BGK solutions
+  double h0; // a parameter for BGK solutions
 
   int fields_every;    // interval for pfd output
   int moments_every;   // interval for pfd_moments output
@@ -133,6 +146,8 @@ struct PscBgkParams
     T_e_coef = parsedParams.getOrDefault<double>("T_e_coef", 1);
 
     maxwellian = parsedParams.getOrDefault<bool>("maxwellian", false);
+    k = parsedParams.getOrDefault<double>("k", .1);
+    h0 = parsedParams.getOrDefault<double>("h0", .9);
 
     n_grid_3 = parsedParams.getOrDefault<int>("n_grid_3", 1);
     box_size_3 = parsedParams.getAndWarnOrDefault<double>("box_size_3", -1);

diff --git a/src/psc_bgk_util/sample_bgk_params.txt b/src/psc_bgk_util/sample_bgk_params.txt
@@ -22,6 +22,8 @@ T_i 0 normally 0, or 1e-6 for boltzmann
 # Run-specific parameters
 path_to_data /CHANGE/ME
 H_x 1
+h0 .9
+k .1
 
 n_grid 16
 n_grid_3 1 number of gridpoints on x axis; 1 for 2D, 2+ for 3D