Merge pull request #318 from JamesMcClung/pr/bgk-exact

Use exact distribution for BGK case

src/psc_bgk.cxx

@@ -209,6 +209,82 @@ 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();
+
+  double B = g.Hx;
+  double sqrt2 = std::sqrt(2);
+
+  double rho_sqr = sqr(rho);
+
+  double gamma = 1 + 8 * g.k * rho_sqr;
+  double alpha = 1 - g.h0 / std::sqrt(gamma) *
+                       std::exp(-g.k * sqr(B) * sqr(rho_sqr) / gamma);
+
+  double mean0 = 0;
+  double stdev0 = 1;
+
+  double mean1 = 4 * g.k * B * rho_sqr * rho / gamma;
+  double stdev1 = 1 / std::sqrt(gamma);
+
+  double m0 = (1 + std::erf((v_phi - mean0) / (stdev0 * sqrt2))) / 2;
+  double m1 = (1 + std::erf((v_phi - mean1) / (stdev1 * sqrt2))) / 2;
+
+  return m0 / alpha + m1 * (1 - 1 / alpha);
+}
+
+struct pdist
+{
+  pdist(double y, double z, double rho)
+    : y{y},
+      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()}
+  {}
+
+  Double3 operator()()
+  {
+    double v_phi = v_phi_dist.get();
+    double v_rho = v_rho_dist.get();
+    double v_x = v_x_dist.get();
+
+    double coef = g.v_e_coef * (g.reverse_v ? -1 : 1) *
+                  (g.reverse_v_half && y < 0 ? -1 : 1);
+    double p_x = coef * g.m_e * v_x;
+    double p_y = coef * g.m_e * (v_phi * -z + v_rho * y) / rho;
+    double p_z = coef * g.m_e * (v_phi * y + v_rho * z) / rho;
+    return Double3{p_x, p_y, p_z};
+  }
+
+private:
+  double y, z, rho;
+  rng::InvertedCdf<double> v_phi_dist;
+  rng::Normal<double> v_rho_dist;
+  rng::Normal<double> v_x_dist;
+};
+
 // ======================================================================
 // initializeParticles
 
@@ -220,43 +296,50 @@ void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts,
 
   auto&& qDensity = -psc::mflds::interior(divGradPhi.grid(), divGradPhi.gt());
 
-  auto npt_init = [&](int kind, Double3 crd, int p, Int3 idx,
-                      psc_particle_npt& npt) {
+  auto init_np = [&](int kind, Double3 crd, int p, Int3 idx,
+                     psc_particle_np& np) {
     double y = getCoord(crd[1]);
     double z = getCoord(crd[2]);
     double rho = sqrt(sqr(y) + sqr(z));
 
+    double Ti = g.T_i;
     switch (kind) {
 
       case KIND_ELECTRON:
-        npt.n = (qDensity(idx[0], idx[1], idx[2], 0, p) -
-                 getIonDensity(rho) * g.q_i) /
-                g.q_e;
-        if (rho != 0) {
-          double v_phi = parsedData->get_interpolated(COL_V_PHI, rho);
+        np.n = (qDensity(idx[0], idx[1], idx[2], 0, p) -
+                getIonDensity(rho) * g.q_i) /
+               g.q_e;
+        if (rho == 0) {
+          double Te = parsedData->get_interpolated(COL_TE, 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 coef = g.v_e_coef * (g.reverse_v ? -1 : 1) *
                         (g.reverse_v_half && y < 0 ? -1 : 1);
-          npt.p[0] = 0;
-          npt.p[1] = coef * g.m_e * v_phi * -z / rho;
-          npt.p[2] = coef * g.m_e * v_phi * y / rho;
+
+          double pz = coef * g.m_e * vphi * y / rho;
+          double py = coef * g.m_e * -vphi * z / rho;
+          np.p = setup_particles.createMaxwellian(
+            {np.kind, np.n, {0, py, pz}, {Te, Te, Te}, np.tag});
         } else {
-          setAll(npt.p, 0);
+          np.p = pdist(y, z, rho);
         }
-        setAll(npt.T, g.T_e_coef * parsedData->get_interpolated(COL_TE, rho));
         break;
 
       case KIND_ION:
-        npt.n = getIonDensity(rho);
-        setAll(npt.p, 0);
-        setAll(npt.T, g.T_i);
+        np.n = getIonDensity(rho);
+        np.p = setup_particles.createMaxwellian(
+          {np.kind, np.n, {0, 0, 0}, {Ti, Ti, Ti}, np.tag});
         break;
 
       default: assert(false);
     }
   };
 
   partitionAndSetupParticles(setup_particles, balance, grid_ptr, mprts,
-                             npt_init);
+                             init_np);
 }
 
 // ======================================================================

src/psc_bgk_util/bgk_params.hxx

@@ -1,13 +1,47 @@
 #pragma once
 
+#include <cmath>
 #include "params_parser.hxx"
 
+// ======================================================================
+// getCalculatedBoxSize
+//
+// Calculate the radius where the spike in the exact distribution function ends,
+// according to the equation (in paper units):
+//    v_phi = max_v = 4*k*B*r^3 / (1 + 8*k*r^2)
+// The exact solution can be decomposed into the difference of two Gaussians.
+// The positive Gaussian has a mean of 0 and stdev of 1, independently of
+// all parameters. "max_v" represents an upper limit for v_phi according
+// to this term.
+// 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 getCalculatedBoxSize(double B, double k)
+{
+  double max_v = 3.;
+  // solve cubic with linear coefficient = 0
+  double a = 4. * k * B;
+  double b = -8. * max_v * k;
+  double d = -max_v;
+
+  double p = -b / (3. * a);
+  double t = -d / (2. * a);
+  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 extra_multiplier = 1.5;
+  double r = (std::cbrt(q + s) + std::cbrt(q - s) + p) * beta;
+  return extra_multiplier * 2 * r;
+}
+
 // ======================================================================
 // PscBgkParams
 
 struct PscBgkParams
 {
-  double box_size; // physical length of region along y and z
+  double box_size; // physical length of region along y and z; -1 -> auto
   double Hx;       // strength of transverse magnetic field
   double q_i;      // ion charge
   double n_i;      // ion number density
@@ -18,6 +52,9 @@ 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
+
   int fields_every;    // interval for pfd output
   int moments_every;   // interval for pfd_moments output
   int gauss_every;     // interval for gauss output/checking
@@ -32,6 +69,7 @@ struct PscBgkParams
   double T_e_coef;     // multiplier for electron temperature
   bool reverse_v;      // whether or not to reverse electron velocity
   bool reverse_v_half; // whether or not to reverse electron velocity for y<0
+  bool maxwellian;     // whether or not to use Maxwellian instead of exact sol
 
   // For 3D cases
   int n_grid_3;      // number of grid points in 3rd dimension
@@ -67,6 +105,8 @@ struct PscBgkParams
     v_e_coef = parsedParams.getOrDefault<double>("v_e_coef", 1);
     T_e_coef = parsedParams.getOrDefault<double>("T_e_coef", 1);
 
+    maxwellian = parsedParams.getOrDefault<bool>("maxwellian", false);
+
     n_grid_3 = parsedParams.getOrDefault<int>("n_grid_3", 1);
     box_size_3 = parsedParams.getOrDefault<double>("box_size_3", 1);
     if (n_grid_3 < parsedParams.get<int>("n_cells_per_patch")) {
@@ -76,5 +116,8 @@ struct PscBgkParams
       if (n_patches_3 <= 0)
         n_patches_3 = n_grid_3 / parsedParams.get<int>("n_cells_per_patch");
     }
+
+    if (box_size <= 0)
+      box_size = getCalculatedBoxSize(Hx, k);
   }
-};
+};