From 64dfe2885d6c9221c1106a605c0139a60492ccfc Mon Sep 17 00:00:00 2001 From: James Date: Wed, 3 Jul 2024 15:00:37 -0400 Subject: [PATCH 1/5] table: +has_column --- src/psc_bgk_util/table.hxx | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/psc_bgk_util/table.hxx b/src/psc_bgk_util/table.hxx index 47da585bd..5c535bab4 100644 --- a/src/psc_bgk_util/table.hxx +++ b/src/psc_bgk_util/table.hxx @@ -148,6 +148,11 @@ public: int n_rows() { return data.size(); } int n_cols() { return header.size(); } + bool has_column(std::string column_name) + { + return header.count(column_name) != 0; + } + double get(std::string column_name, int row) { int col = header[column_name]; From efc2aeeca13ed09d292d049dedf1574a9f78f8e6 Mon Sep 17 00:00:00 2001 From: James Date: Wed, 3 Jul 2024 15:01:18 -0400 Subject: [PATCH 2/5] bgk: use B_x from ic_table when present --- src/psc_bgk.cxx | 27 +++++++++++++++++++++------ 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/src/psc_bgk.cxx b/src/psc_bgk.cxx index 75d049287..e13bba3a6 100644 --- a/src/psc_bgk.cxx +++ b/src/psc_bgk.cxx @@ -207,6 +207,9 @@ double v_phi_cdf(double v_phi, double rho) rho /= g.beta; double B = g.Hx; + if (ic_table->has_column("B_x")) { + B = ic_table->get_interpolated("B_x", "rho", rho); + } double sqrt2 = std::sqrt(2); double rho_sqr = sqr(rho); @@ -444,12 +447,24 @@ void initializeDivGradPhi(BgkMfields& gradPhi, BgkMfields& divGradPhi) void initializeFields(MfieldsState& mflds, BgkMfields& gradPhi) { - setupFields(mflds, [&](int m, double crd[3]) { - switch (m) { - case HX: return g.Hx; - default: return 0.; - } - }); + if (ic_table->has_column("B_x")) { + setupFields(mflds, [&](int m, double crd[3]) { + double y = getCoord(crd[1]); + double z = getCoord(crd[2]); + double rho = sqrt(sqr(y) + sqr(z)); + switch (m) { + case HX: return ic_table->get_interpolated("B_x", "rho", rho); + default: return 0.; + } + }); + } else { + setupFields(mflds, [&](int m, double crd[3]) { + switch (m) { + case HX: return g.Hx; + default: return 0.; + } + }); + } // initialize E separately mflds.storage().view(_all, _all, _all, _s(EX, EX + 3)) = -gradPhi.gt(); From 438f79686729dc1b15d1d71d26065bd60de09f31 Mon Sep 17 00:00:00 2001 From: James Date: Fri, 5 Jul 2024 11:38:51 -0400 Subject: [PATCH 3/5] bgk: use A_phi when present --- src/psc_bgk.cxx | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/psc_bgk.cxx b/src/psc_bgk.cxx index e13bba3a6..981ef1832 100644 --- a/src/psc_bgk.cxx +++ b/src/psc_bgk.cxx @@ -270,7 +270,10 @@ struct pdist_case4 z{z}, rho{rho}, v_phi_dist{g.beta * 8.0 * g.k * sqr(rho / g.beta) * - (0.5 * g.Hx * rho / g.beta) / + ((ic_table->has_column("A_phi") + ? ic_table->get_interpolated("A_phi", "rho", rho) + : 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}, From 925c4158acadb5bcf1078348ea4a4816fa785b1e Mon Sep 17 00:00:00 2001 From: James Date: Fri, 5 Jul 2024 11:58:26 -0400 Subject: [PATCH 4/5] bgk: use A_phi instead of B --- src/psc_bgk.cxx | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/psc_bgk.cxx b/src/psc_bgk.cxx index 981ef1832..cb3f30ceb 100644 --- a/src/psc_bgk.cxx +++ b/src/psc_bgk.cxx @@ -202,26 +202,27 @@ inline double getIonDensity(double rho) double v_phi_cdf(double v_phi, double rho) { + double A_phi = ic_table->has_column("A_phi") + ? ic_table->get_interpolated("A_phi", "rho", rho) + : .5 * g.Hx * rho; + // convert units from psc to paper v_phi /= g.beta; rho /= g.beta; + A_phi /= g.beta; - double B = g.Hx; - if (ic_table->has_column("B_x")) { - B = ic_table->get_interpolated("B_x", "rho", rho); - } 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 alpha = + 1 - g.h0 / std::sqrt(gamma) * std::exp(-4 * g.k * sqr(A_phi * rho) / gamma); double mean0 = 0; double stdev0 = 1; - double mean1 = 4 * g.k * B * rho_sqr * rho / gamma; + double mean1 = 8 * g.k * A_phi * rho_sqr / gamma; double stdev1 = 1 / std::sqrt(gamma); double m0 = (1 + std::erf((v_phi - mean0) / (stdev0 * sqrt2))) / 2; From 40af70086460c321c9b669159d66fb32d72af12f Mon Sep 17 00:00:00 2001 From: James Date: Fri, 5 Jul 2024 11:58:54 -0400 Subject: [PATCH 5/5] bgk: inline some vars --- src/psc_bgk.cxx | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/src/psc_bgk.cxx b/src/psc_bgk.cxx index cb3f30ceb..208aedbdb 100644 --- a/src/psc_bgk.cxx +++ b/src/psc_bgk.cxx @@ -211,22 +211,18 @@ double v_phi_cdf(double v_phi, double rho) rho /= g.beta; A_phi /= g.beta; - double sqrt2 = std::sqrt(2); - - double rho_sqr = sqr(rho); - - double gamma = 1 + 8 * g.k * rho_sqr; + double gamma = 1 + 8 * g.k * sqr(rho); double alpha = 1 - g.h0 / std::sqrt(gamma) * std::exp(-4 * g.k * sqr(A_phi * rho) / gamma); double mean0 = 0; double stdev0 = 1; - double mean1 = 8 * g.k * A_phi * rho_sqr / gamma; + double mean1 = 8 * g.k * A_phi * 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; + double m0 = (1 + std::erf((v_phi - mean0) / (stdev0 * std::sqrt(2)))) / 2; + double m1 = (1 + std::erf((v_phi - mean1) / (stdev1 * std::sqrt(2)))) / 2; return m0 / alpha + m1 * (1 - 1 / alpha); }