From cfb27333d6f7829e74c26dfcb9244234f7a59680 Mon Sep 17 00:00:00 2001 From: James Date: Tue, 16 Jan 2024 13:13:01 -0500 Subject: [PATCH 1/4] bgk_params: load h0, k from file --- src/psc_bgk_util/bgk_params.hxx | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/psc_bgk_util/bgk_params.hxx b/src/psc_bgk_util/bgk_params.hxx index b6d35cacd..ec35d418c 100644 --- a/src/psc_bgk_util/bgk_params.hxx +++ b/src/psc_bgk_util/bgk_params.hxx @@ -76,8 +76,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 +133,8 @@ struct PscBgkParams T_e_coef = parsedParams.getOrDefault("T_e_coef", 1); maxwellian = parsedParams.getOrDefault("maxwellian", false); + k = parsedParams.getOrDefault("k", .1); + h0 = parsedParams.getOrDefault("h0", .9); n_grid_3 = parsedParams.getOrDefault("n_grid_3", 1); box_size_3 = parsedParams.getAndWarnOrDefault("box_size_3", -1); From bd5609bd1770c6cc18a5a4e1dd4d2f115dc93c23 Mon Sep 17 00:00:00 2001 From: James Date: Tue, 16 Jan 2024 13:13:50 -0500 Subject: [PATCH 2/4] sample_bgk_params: +h0, k --- src/psc_bgk_util/sample_bgk_params.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/psc_bgk_util/sample_bgk_params.txt b/src/psc_bgk_util/sample_bgk_params.txt index 99e046131..1934c9f67 100644 --- 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 From e40b9f0417e454908a773457f96c99d388974cbf Mon Sep 17 00:00:00 2001 From: James Date: Tue, 30 Jan 2024 11:19:04 -0500 Subject: [PATCH 3/4] bgk, bgk_params: mv get_beta --- src/psc_bgk.cxx | 21 ++++----------------- src/psc_bgk_util/bgk_params.hxx | 14 ++++++++++++++ 2 files changed, 18 insertions(+), 17 deletions(-) diff --git a/src/psc_bgk.cxx b/src/psc_bgk.cxx index 9237d958b..67b04ceb8 100644 --- a/src/psc_bgk.cxx +++ b/src/psc_bgk.cxx @@ -208,19 +208,6 @@ 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 @@ -228,8 +215,8 @@ double get_beta() 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 index ec35d418c..746d8b6e8 100644 --- 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 From 9470f8876cf6f7b02935e39f7808595174a889b0 Mon Sep 17 00:00:00 2001 From: James Date: Tue, 30 Jan 2024 11:19:52 -0500 Subject: [PATCH 4/4] bgk_params: don't hardcode beta --- src/psc_bgk_util/bgk_params.hxx | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/psc_bgk_util/bgk_params.hxx b/src/psc_bgk_util/bgk_params.hxx index 746d8b6e8..acec7a8b8 100644 --- a/src/psc_bgk_util/bgk_params.hxx +++ b/src/psc_bgk_util/bgk_params.hxx @@ -31,7 +31,7 @@ double get_beta(ParsedData& data) // 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 @@ -44,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; } @@ -66,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);