Merge pull request #333 from JamesMcClung/pr/bgk-case4

BGK case 4

src/psc_bgk.cxx

@@ -131,8 +131,8 @@ Grid_t* setupGrid()
   kinds[KIND_ELECTRON] = {g.q_e, g.m_e, "e"};
   kinds[KIND_ION] = {g.q_i, g.m_i, "i"};
 
-  mpi_printf(MPI_COMM_WORLD, "lambda_D = %g\n",
-             sqrt(parsedData->get_interpolated(COL_TE, 0)));
+  // lambda_D = v_e / omega_pe = v_e = beta
+  mpi_printf(MPI_COMM_WORLD, "lambda_D = %g\n", g.beta);
 
   // --- generic setup
   auto norm_params = Grid_t::NormalizationParams::dimensionless();
@@ -215,8 +215,8 @@ inline double getIonDensity(double rho)
 double v_phi_cdf(double v_phi, double rho)
 {
   // convert units from psc to paper
-  v_phi /= get_beta(*parsedData);
-  rho /= get_beta(*parsedData);
+  v_phi /= g.beta;
+  rho /= g.beta;
 
   double B = g.Hx;
   double sqrt2 = std::sqrt(2);
@@ -246,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(*parsedData)},
-      v_x_dist{0, get_beta(*parsedData)}
+      v_rho_dist{0, g.beta},
+      v_x_dist{0, g.beta}
   {}
 
   Double3 operator()()
@@ -271,6 +271,65 @@ struct pdist
   rng::Normal<double> v_x_dist;
 };
 
+struct pdist_case4
+{
+  pdist_case4(double p_background, double y, double z, double rho)
+    : p_background{p_background},
+      y{y},
+      z{z},
+      rho{rho},
+      v_phi_dist{g.beta * 8.0 * g.k * sqr(rho / g.beta) *
+                   (0.5 * g.Hx * rho / g.beta) /
+                   (1.0 + 8.0 * g.k * sqr(rho / g.beta)),
+                 g.beta / sqrt(1.0 + 8.0 * g.k * sqr(rho / g.beta))},
+      v_rho_dist{0, g.beta},
+      v_x_dist{g.beta * 2.0 * g.xi * g.A_x0 / (1.0 + 2.0 * g.xi),
+               g.beta / sqrt(1.0 + 2.0 * g.xi)},
+      simple_dist{0.0, g.beta},
+      uniform{0.0, 1.0}
+  {}
+
+  Double3 operator()()
+  {
+    double coef = g.v_e_coef * (g.reverse_v ? -1 : 1) *
+                  (g.reverse_v_half && y < 0 ? -1 : 1);
+    if (rho == 0) {
+      //  p_y and p_z reduce to same case at rho=0, but not p_x
+      double p_x =
+        coef * g.m_e *
+        (uniform.get() < p_background ? simple_dist : v_x_dist).get();
+      double p_y = coef * g.m_e * simple_dist.get();
+      double p_z = coef * g.m_e * simple_dist.get();
+      return Double3{p_x, p_y, p_z};
+    }
+
+    double v_phi, v_rho, v_x;
+
+    if (uniform.get() < p_background) {
+      v_phi = simple_dist.get();
+      v_rho = simple_dist.get();
+      v_x = simple_dist.get();
+    } else {
+      v_phi = v_phi_dist.get();
+      v_rho = v_rho_dist.get();
+      v_x = v_x_dist.get();
+    }
+
+    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 p_background, y, z, rho;
+  rng::Normal<double> v_phi_dist;
+  rng::Normal<double> v_rho_dist;
+  rng::Normal<double> v_x_dist;
+  rng::Normal<double> simple_dist;
+  rng::Uniform<double> uniform;
+};
+
 // ======================================================================
 // initializeParticles
 
@@ -295,6 +354,20 @@ void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts,
         np.n = (qDensity(idx[0], idx[1], idx[2], 0, p) -
                 getIonDensity(rho) * g.q_i) /
                g.q_e;
+
+        if (g.xi != 0) {
+          // case 4: sum of two maxwellians, where one has easy moments and the
+          // other's moments are given in input file
+
+          double psi_cs =
+            parsedData->get_interpolated(COL_PHI, 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);
+
+          break;
+        }
+
         if (rho == 0) {
           double Te = parsedData->get_interpolated(COL_TE, rho);
           np.p = setup_particles.createMaxwellian(

src/psc_bgk_util/bgk_params.hxx

@@ -4,20 +4,6 @@
 #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
 
@@ -49,23 +35,23 @@ double getSpikeSize(double B, double k, double beta)
 }
 
 // The hole size is determined empirically from the input data. This function
-// finds where the electron density <= 1 + epsilon, where epsilon is small.
+// finds where |electron density - 1| <= epsilon, where epsilon is small.
 double getHoleSize(ParsedData& data)
 {
   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 (data[row][COL_NE] > 1 + epsilon)
+    if (std::abs(data[row][COL_NE] - 1) > epsilon)
       return data[row][COL_RHO];
   }
   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, ParsedData& data)
+double getCalculatedBoxSize(double B, double k, double beta, ParsedData& data)
 {
-  double spike_size = getSpikeSize(B, k, get_beta(data));
+  double spike_size = getSpikeSize(B, k, beta);
   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);
@@ -91,6 +77,11 @@ struct PscBgkParams
 
   double k;  // a parameter for BGK solutions
   double h0; // a parameter for BGK solutions
+  double xi; // a parameter for BGK solutions
+
+  double beta; // electron thermal velocity / c
+
+  double A_x0; // vector potential, used when xi != 0
 
   int fields_every;    // interval for pfd output
   int moments_every;   // interval for pfd_moments output
@@ -145,9 +136,19 @@ struct PscBgkParams
     v_e_coef = parsedParams.getOrDefault<double>("v_e_coef", 1);
     T_e_coef = parsedParams.getOrDefault<double>("T_e_coef", 1);
 
+    beta = parsedParams.getOrDefault<double>("beta", 1e-3);
+
     maxwellian = parsedParams.getOrDefault<bool>("maxwellian", false);
     k = parsedParams.getOrDefault<double>("k", .1);
     h0 = parsedParams.getOrDefault<double>("h0", .9);
+    xi = parsedParams.getOrDefault<double>("xi", 0.0);
+
+    if (xi == 0.0) {
+      parsedParams.warnIfPresent("A_x0",
+                                 "A_x0 is only used in cases when xi != 0");
+    } else {
+      A_x0 = parsedParams.get<double>("A_x0");
+    }
 
     n_grid_3 = parsedParams.getOrDefault<int>("n_grid_3", 1);
     box_size_3 = parsedParams.getAndWarnOrDefault<double>("box_size_3", -1);
@@ -161,8 +162,8 @@ struct PscBgkParams
     }
 
     if (box_size <= 0)
-      box_size = rel_box_size * getCalculatedBoxSize(Hx, k, data);
+      box_size = rel_box_size * getCalculatedBoxSize(Hx, k, beta, data);
     if (box_size_3 <= 0)
-      box_size_3 = rel_box_size_3 * getCalculatedBoxSize(Hx, k, data);
+      box_size_3 = rel_box_size_3 * getCalculatedBoxSize(Hx, k, beta, data);
   }
 };