From 6053c28b8d67388f9d16a0fb5518f6b69141bc41 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Fri, 17 Jul 2026 22:47:35 -0500 Subject: [PATCH 01/14] add pinning tests for simplified gas kinetics model Co-Authored-By: Claude Fable 5 --- tests/plume/test_simplified_gaskinetics.py | 161 +++++++++++++++++++++ 1 file changed, 161 insertions(+) create mode 100644 tests/plume/test_simplified_gaskinetics.py diff --git a/tests/plume/test_simplified_gaskinetics.py b/tests/plume/test_simplified_gaskinetics.py new file mode 100644 index 0000000..30c9ac0 --- /dev/null +++ b/tests/plume/test_simplified_gaskinetics.py @@ -0,0 +1,161 @@ +# ======================== +# PyRPOD: tests/plume/test_simplified_gaskinetics.py +# ======================== +# Pinning tests for the SimplifiedGasKinetics class (Cai & Wang 2012, +# "Numerical Validations for a Set of Collisionless Rocket Plume Solutions", +# JSR 49(1), DOI 10.2514/1.A32046). +# +# Reference values were computed independently from the paper's equations +# with mpmath at 40 significant digits (Eqs. 13-19 and 22-24 implemented +# directly from the printed formulas, not from pyrpod code). They pin the +# known-correct behavior of the simplified far-field model so that later +# refactors (sympy removal, overflow-safe exponentials) can be verified +# to preserve behavior. + +import numpy as np +import pytest + +from pyrpod.plume.RarefiedPlumeGasKinetics import SimplifiedGasKinetics + +pytestmark = pytest.mark.plume + +# Fixed thruster gas properties used for all cases (argon-like). +R_SPECIFIC = 208.13 # J / (kg K) +T_0 = 500.0 # K +GAMMA = 5.0 / 3.0 +N_0 = 1.0e20 # m^-3 +D_NOZZLE = 0.2 # m (R_0 = 0.1 m) +T_W = 300.0 # K +SIGMA = 1.0 + +# Relative tolerance for pinned closed-form values: the implementation +# evaluates the same closed-form expressions in float64 (machine epsilon +# ~1e-16, tens of flops => accumulated rounding well below 1e-12), while +# the references carry 17 correct digits. 1e-9 leaves ample headroom for +# rounding-level refactors (e.g. sympy.erf -> scipy.special.erf, combined +# exponentials) yet fails loudly on any actual formula change. +RTOL_PINNED = 1e-9 + + +def make_plume(distance, theta, S_0): + """Build a SimplifiedGasKinetics instance with exit speed ratio S_0.""" + ve = S_0 * np.sqrt(2 * R_SPECIFIC * T_0) + thruster_characteristics = { + 'd': D_NOZZLE, 've': ve, 'R': R_SPECIFIC, + 'gamma': GAMMA, 'Te': T_0, 'n': N_0, + } + return SimplifiedGasKinetics(distance, theta, thruster_characteristics, + T_W, SIGMA) + + +# (distance, theta, S_0) -> (n/n0, U*sqrt(beta0), W*sqrt(beta0), T/T0) +# per Eqs. 14-17 with Q' = X^2/(X^2+Z^2) (Eq. 13). +FIELD_CASES = [ + ((1.0, 0.2, 2.0), (0.036338635060152093, 2.3636427215997389, + 0.47913410002529834, 0.27393936134733628)), + ((2.0, 0.7, 2.0), (0.0010317934247266036, 1.5832262838600719, + 1.3335331025390796, 0.2543577108315116)), + ((5.0, 1.2, 2.0), (4.4418364832646119e-6, 0.54854957955069592, + 1.4109526908580327, 0.20360334802441553)), + ((1.5, 0.4, 1.0), (0.0046772248683472122, 1.5085528987977848, + 0.63780593571949983, 0.21734975912503903)), + ((3.0, 0.9, 3.0), (1.0977476844642729e-5, 1.4507999313815155, + 1.8282374555518849, 0.27008477291007084)), +] + +# (X, S_0) -> (n/n0, U*sqrt(beta0)) per Eqs. 18-19 (centerline). +CENTERLINE_CASES = [ + ((0.5, 2.0), (0.15923868831612018, 2.4062627987321269)), + ((1.0, 2.0), (0.043598168447392599, 2.434502010083779)), + ((5.0, 2.0), (0.0017976259018792059, 2.4441625633761652)), + ((1.0, 1.0), (0.014624592798873514, 1.6841200399767633)), + ((1.0, 3.0), (0.089793762869630585, 3.3010882334220839)), +] + +# S_0 -> (U_inf*sqrt(beta0), T_inf/T0) per Eqs. 22-24, evaluated with +# mpmath. (Eq. 23's denominator is 3*[S0 + (1/2+S0^2)*sqrt(pi)* +# (1+erf(S0))*exp(S0^2)], i.e. 3*K(Q=1), consistent with +# -2/3*G^2 + 4N(1)/(3K(1)).) +ASYMPTOTE_CASES = { + 1.0: (1.689948557881606, 0.22268161973032418), + 2.0: (2.444572023854802, 0.2754744452640276), + 3.0: (3.3157896665554836, 0.30193859122021722), +} + + +@pytest.mark.parametrize("point,expected", FIELD_CASES) +def test_field_solution_pinned(point, expected): + """Pin Eqs. 14-17 (off-centerline simplified field solutions).""" + distance, theta, S_0 = point + n_ref, U_ref, W_ref, T_ref = expected + plume = make_plume(distance, theta, S_0) + assert plume.get_num_density_ratio() == pytest.approx(n_ref, rel=RTOL_PINNED) + assert plume.get_U_normalized() == pytest.approx(U_ref, rel=RTOL_PINNED) + assert plume.get_W_normalized() == pytest.approx(W_ref, rel=RTOL_PINNED) + assert plume.get_temp_ratio() == pytest.approx(T_ref, rel=RTOL_PINNED) + + +@pytest.mark.parametrize("point,expected", CENTERLINE_CASES) +def test_centerline_solution_pinned(point, expected): + """Pin Eqs. 18-19 (centerline density and velocity closed forms).""" + X, S_0 = point + n_ref, U_ref = expected + plume = make_plume(X, 0.0, S_0) + assert plume.get_num_density_centerline() == pytest.approx(n_ref, + rel=RTOL_PINNED) + assert plume.get_velocity_centerline() == pytest.approx(U_ref, + rel=RTOL_PINNED) + + +@pytest.mark.parametrize("S_0", [1.0, 2.0, 3.0]) +def test_centerline_velocity_approaches_asymptote(S_0): + """Centerline U converges to the Eq. 22/24 far-field asymptote. + + Convergence is O((R_0/X)^2): at X/R_0 = 1000 the residual is ~4e-7 + relative, so rel=1e-5 passes with an order-of-magnitude margin while + still requiring genuine convergence. + """ + U_inf = ASYMPTOTE_CASES[S_0][0] + X = 1000 * (D_NOZZLE / 2) + plume = make_plume(X, 0.0, S_0) + assert plume.get_velocity_centerline() == pytest.approx(U_inf, rel=1e-5) + + +@pytest.mark.parametrize("S_0", [1.0, 2.0, 3.0]) +def test_centerline_temperature_approaches_asymptote(S_0): + """Centerline T converges to the Eq. 23 far-field asymptote. + + T = -2/3*U^2 + 4/3*/ subtracts two O(S_0^2) quantities, so the + O((R_0/X)^2 * S_0^2) residual of each term is amplified by ~U^2/T + (a factor ~25 at S_0 = 3) in the relative error of T. At + X/R_0 = 5000 the amplified residual is <~1e-5 for S_0 <= 3, so + rel=1e-4 passes with margin for both the constant-N approximation + and the exact Eq. 21 quadrature, while still requiring convergence. + """ + T_inf = ASYMPTOTE_CASES[S_0][1] + X = 5000 * (D_NOZZLE / 2) + plume = make_plume(X, 0.0, S_0) + assert plume.get_temp_centerline() == pytest.approx(T_inf, rel=1e-4) + + +@pytest.mark.parametrize("point,expected", [ + ((1.0, 0.2, 2.0), (0.36855461276327334, -0.056815859301478394, + 143.3923778236429)), + ((2.0, 0.7, 2.0), (0.005299929030826453, -0.003007472717182616, + 1.9141355894449432)), +]) +def test_surface_interaction_pinned(point, expected): + """Pin the Maxwell gas-surface interface outputs (off-centerline). + + Unlike the field-solution pins above, these reference values were + generated from the current implementation (not independently), so + they pin implementation behavior for refactor safety rather than + correctness against the paper. Tolerance rationale as RTOL_PINNED. + """ + distance, theta, S_0 = point + p_ref, tau_ref, q_ref = expected + plume = make_plume(distance, theta, S_0) + assert float(plume.get_pressure()) == pytest.approx(p_ref, rel=RTOL_PINNED) + assert float(plume.get_shear_pressure()) == pytest.approx(tau_ref, + rel=RTOL_PINNED) + assert float(plume.get_heat_flux()) == pytest.approx(q_ref, rel=RTOL_PINNED) From 720f675a0239eaf1e01ea4e528c04b0695c73e3c Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Fri, 17 Jul 2026 22:50:58 -0500 Subject: [PATCH 02/14] replace sympy with scipy and remove dead collisionless plume code Co-Authored-By: Claude Fable 5 --- pyrpod/plume/RarefiedPlumeGasKinetics.py | 382 ++++++++--------------- 1 file changed, 122 insertions(+), 260 deletions(-) diff --git a/pyrpod/plume/RarefiedPlumeGasKinetics.py b/pyrpod/plume/RarefiedPlumeGasKinetics.py index e571e27..d02282d 100644 --- a/pyrpod/plume/RarefiedPlumeGasKinetics.py +++ b/pyrpod/plume/RarefiedPlumeGasKinetics.py @@ -35,14 +35,91 @@ """ import numpy as np -import sympy as sp from scipy import integrate -import matplotlib.pyplot as plt +from scipy.special import erf #define constants AVOGADROS_NUMBER = 6.0221e23 GAS_CONSTANT = 8.314 + +def get_K_factor(Q, S_0): + ''' + Scaled special factor exp(-S_0^2) * K [Cai & Wang 2012, Eq. 10]. + + The exp(-S_0^2) prefactor of the field solutions (Eqs. 5-8, 14) + is combined analytically with Eq. 10's exp(Q * S_0^2) into + exp(-S_0^2 * (1 - Q)). Since Q <= 1 both exponentials are <= 1, + so this never overflows for large speed ratios, while the plain + exp(Q * S_0^2) would. The [1 + erf(...)] factor is bounded by 2 + and needs no scaling. + + Parameters + ---------- + Q : float + special factor Q (full or simplified), 0 < Q <= 1 + S_0 : float + molecular speed ratio at the nozzle exit + + Returns + ------- + float + exp(-S_0^2) * K(Q, S_0) + ''' + erf_term = (1 + erf(S_0 * np.sqrt(Q))) * np.exp(-S_0 ** 2 * (1 - Q)) + term1 = Q * S_0 * np.exp(-S_0 ** 2) + term2 = (0.5 + Q * S_0 ** 2) * np.sqrt(np.pi * Q) + return Q * (term1 + term2 * erf_term) + + +def get_M_factor(Q, S_0): + ''' + Scaled special factor exp(-S_0^2) * M [Cai & Wang 2012, Eq. 11]. + + Overflow-safe exponential combination as in get_K_factor. + + Parameters + ---------- + Q : float + special factor Q (full or simplified), 0 < Q <= 1 + S_0 : float + molecular speed ratio at the nozzle exit + + Returns + ------- + float + exp(-S_0^2) * M(Q, S_0) + ''' + erf_term = (1 + erf(S_0 * np.sqrt(Q))) * np.exp(-S_0 ** 2 * (1 - Q)) + term1 = (1 + Q * S_0 ** 2) * np.exp(-S_0 ** 2) + term2 = S_0 * (1.5 + Q * S_0 ** 2) * np.sqrt(np.pi * Q) + return Q ** 2 * (term1 + term2 * erf_term) + + +def get_N_factor(Q, S_0): + ''' + Scaled special factor exp(-S_0^2) * N [Cai & Wang 2012, Eq. 12]. + + Overflow-safe exponential combination as in get_K_factor. + + Parameters + ---------- + Q : float + special factor Q (full or simplified), 0 < Q <= 1 + S_0 : float + molecular speed ratio at the nozzle exit + + Returns + ------- + float + exp(-S_0^2) * N(Q, S_0) + ''' + erf_term = (1 + erf(S_0 * np.sqrt(Q))) * np.exp(-S_0 ** 2 * (1 - Q)) + term1 = S_0 * Q ** 2 * (1.25 + Q * S_0 ** 2 / 2) * np.exp(-S_0 ** 2) + term2 = 0.5 * np.sqrt(np.pi * Q ** 3) + term3 = 0.75 + 3 * Q * S_0 ** 2 + Q ** 2 * S_0 ** 4 + return term1 + term2 * term3 * erf_term + def get_maxwellian_pressure(rho_inf, U, S, sigma, theta, T, T_w): ''' Rarefied Gas Dynamics - Shen - eq. 4.19 @@ -79,7 +156,7 @@ def get_maxwellian_pressure(rho_inf, U, S, sigma, theta, T, T_w): p1 *= np.exp(- (S * np.cos(theta)) ** 2) p2 = (2 - sigma) * ((S * np.cos(theta)) ** 2 + 0.5) p2 += (S * np.cos(theta) * (sigma / 2) * np.sqrt(np.pi * T_w / T)) - p2 *= 1 + sp.erf(S * np.cos(theta)) + p2 *= 1 + erf(S * np.cos(theta)) p = p1 + p2 p *= (rho_inf * U ** 2) / (2 * S ** 2) return p @@ -111,7 +188,7 @@ def get_maxwellian_shear_pressure(rho_inf, U, S, sigma, theta): tau1 = np.exp(- (S * np.cos(theta)) ** 2) tau2 = np.sqrt(np.pi) * S * np.cos(theta) - tau2 *= (1 + sp.erf(S * np.cos(theta))) + tau2 *= (1 + erf(S * np.cos(theta))) tau = tau1 + tau2 tau *= -(sigma * rho_inf * np.sin(theta) * U ** 2) / (2 * np.sqrt(np.pi) * S) return tau @@ -150,7 +227,7 @@ def get_maxwellian_heat_transfer(rho_inf, S, sigma, theta, T, T_r, R, gamma): the pressure exerted on the surface element (N / m^2) ''' q = (S ** 2) + (gamma / (gamma - 1)) - (((gamma + 1) * T_r) / (2 * (gamma - 1) * T)) - q *= np.exp(- (S * np.cos(theta)) ** 2) + (np.sqrt(np.pi) * (S * np.cos(theta)) * (1 + sp.erf(S * np.cos(theta)))) + q *= np.exp(- (S * np.cos(theta)) ** 2) + (np.sqrt(np.pi) * (S * np.cos(theta)) * (1 + erf(S * np.cos(theta)))) q -= 0.5 * np.exp(- (S * np.cos(theta)) ** 2) q *= sigma * rho_inf * R * T * np.sqrt(R * T / (2 * np.pi)) return q @@ -302,7 +379,10 @@ def get_plume_angular_density_decay_function(self, theta): def set_normalization_constant(self): ''' From Lumpkin 1999. Setter for normalization constant. - Integrates theta from 0 to the limiting turn angle. + Numerically integrates sin(theta) * cos^kappa(pi*theta/(2*theta_max)) + from 0 to the limiting turn angle. The integrand is real and + non-negative on [0, theta_max] for any kappa, so no complex + arithmetic can arise. Parameters ---------- @@ -313,17 +393,13 @@ def set_normalization_constant(self): None. ''' theta_max = self.get_limiting_turn_angle() - theta = sp.symbols('theta') - f = (sp.cos((sp.pi / 2) * (theta / theta_max))) ** (2 / (self.gamma - 1)) - integrand = sp.sin(theta) * f - integral = sp.integrate(integrand, (theta, 0, theta_max)) #integrate from 0 to max turning angle - A = 0.5 * np.sqrt((self.gamma - 1) / (self.gamma + 1)) / (integral.evalf()) - - # if A is complex, take the real part - if type(A) == sp.Add: - A = sp.re(A) + kappa = 2 / (self.gamma - 1) + + def integrand(theta): + return np.sin(theta) * np.cos((np.pi / 2) * (theta / theta_max)) ** kappa - self.A = A + integral, _ = integrate.quad(integrand, 0, theta_max) + self.A = 0.5 * np.sqrt((self.gamma - 1) / (self.gamma + 1)) / integral return def get_sonic_velocity(self): @@ -600,8 +676,11 @@ def set_Q_simple(self): def set_K_simple(self): ''' - Setter for simplified special factor K. This simplification is just the - substitution of Q for Q'. + Setter for simplified special factor K [Cai & Wang 2012, Eq. 10] + with Q substituted by Q'. Stored scaled by exp(-S_0^2) for + overflow safety (see get_K_factor); ratio methods (U, W, T) + are unaffected since the scaling cancels, and the density + method uses the scaled factor directly. Parameters ---------- @@ -611,21 +690,15 @@ def set_K_simple(self): ------- None. ''' - #K_simple = Q_simple * ((Q_simple * S_0) + ((0.5 + (Q_simple * S_0 ** 2)) * np.sqrt(np.pi * Q_simple) * - #(1 + sp.erf(S_0 * np.sqrt(Q_simple))) ** (Q_simple * S_0 ** 2))) - term1 = self.Q_simple * self.S_0 - term2 = 0.5 + self.Q_simple * self.S_0 ** 2 - term3 = np.sqrt(np.pi * self.Q_simple) - term4 = (1 + sp.erf(self.S_0 * np.sqrt(self.Q_simple))) * np.exp(self.Q_simple * self.S_0 ** 2) - K_simple = self.Q_simple * (term1 + term2 * term3 * term4) - self.K_simple = K_simple + self.K_simple = get_K_factor(self.Q_simple, self.S_0) return - + def set_M_simple(self): ''' - Setter for simplified special factor M. This simplification is just the - substitution of Q for Q'. + Setter for simplified special factor M [Cai & Wang 2012, Eq. 11] + with Q substituted by Q'. Stored scaled by exp(-S_0^2) for + overflow safety (see get_K_factor). Paramters --------- @@ -635,21 +708,15 @@ def set_M_simple(self): ------- None. ''' - #M_simple = (Q_simple ** 2) * ((Q_simple * S_0 ** 2) + 1 + (S_0 * (1.5 + (Q_simple * S_0 ** 2)) * - #np.sqrt(np.pi * Q_simple)) * (1 + sp.erf(S_0 * np.sqrt(Q_simple))) ** (Q_simple *S_0 ** 2)) - term1 = 1 + self.Q_simple * self.S_0 ** 2 - term2 = self.S_0 * (1.5 + self.Q_simple * self.S_0 ** 2) - term3 = np.sqrt(np.pi * self.Q_simple) - term4 = (1 + sp.erf(self.S_0 * np.sqrt(self.Q_simple))) * np.exp(self.Q_simple * self.S_0 ** 2) - M_simple = self.Q_simple ** 2 * (term1 + term2 * term3 * term4) - self.M_simple = M_simple + self.M_simple = get_M_factor(self.Q_simple, self.S_0) return - + def set_N_simple(self): ''' - Setter for simplified special factor N. This simplification is just the - substitution of Q for Q'. + Setter for simplified special factor N [Cai & Wang 2012, Eq. 12] + with Q substituted by Q'. Stored scaled by exp(-S_0^2) for + overflow safety (see get_K_factor). Parameters ---------- @@ -659,14 +726,7 @@ def set_N_simple(self): ------- None. ''' - #N_simple = S_0 * (Q_simple ** 2) * (1.25 + (Q_simple * S_0 ** 2) / 2) - #N_simple += (0.5 * np.sqrt(np.pi * Q_simple ** 3)) * (0.75 + 3 * Q_simple * S_0 **2 + Q_simple ** 2 * S_0 ** 4) * (1 + sp.erf(S_0 * np.sqrt(Q_simple))) ** (Q_simple * S_0 ** 2) - term1 = self.S_0 * self.Q_simple ** 2 * (1.25 + self.Q_simple * self.S_0 ** 2 / 2) - term2 = 0.5 * np.sqrt(np.pi * self.Q_simple ** 3) - term3 = 0.75 + 3 * self.Q_simple * self.S_0 ** 2 + self.Q_simple ** 2 * self.S_0 ** 4 - term4 = (1 + sp.erf(self.S_0 * np.sqrt(self.Q_simple))) * np.exp(self.Q_simple * self.S_0 ** 2) - N_simple = term1 + term2 * term3 * term4 - self.N_simple = N_simple + self.N_simple = get_N_factor(self.Q_simple, self.S_0) return @@ -684,10 +744,9 @@ def get_num_density_ratio(self): float number density at a point (X, 0, Z) vs number density at the nozzle exit ''' - # num_density_ratio = n_1s(X, 0, Z) / n_0 + # num_density_ratio = n_1s(X, 0, Z) / n_0 [Cai & Wang 2012, Eq. 14] + # K_simple already carries the exp(-S_0^2) prefactor (see set_K_simple) num_density_ratio = (self.K_simple / (2 * np.sqrt(np.pi)) * (self.R_0 / self.X) ** 2) - num_density_ratio *= np.exp(-(self.S_0 ** 2)) - num_density_ratio = float(num_density_ratio) return num_density_ratio def get_U_normalized(self): @@ -705,9 +764,8 @@ def get_U_normalized(self): float returns U normalized ''' - # U_normalized = U_1s (X, 0, Z) * sqrt(beta) + # U_normalized = U_1s (X, 0, Z) * sqrt(beta) [Cai & Wang 2012, Eq. 15] U_normalized = self.M_simple / self.K_simple - U_normalized = float(U_normalized) return U_normalized def get_W_normalized(self): @@ -725,9 +783,8 @@ def get_W_normalized(self): float returns W normalized ''' - # W_normalized = W_1s (X, 0, Z) * sqrt(beta) + # W_normalized = W_1s (X, 0, Z) * sqrt(beta) [Cai & Wang 2012, Eq. 16] W_normalized = (self.M_simple / self.K_simple) * (self.Z / self.X) - W_normalized = float(W_normalized) return W_normalized def get_temp_ratio(self): @@ -744,11 +801,9 @@ def get_temp_ratio(self): float ratio of temperature at a point (X, 0, Z) to the temperature at the nozzle exit ''' - # T_ratio = T_1s / T_0 + # T_ratio = T_1s / T_0 [Cai & Wang 2012, Eq. 17] T_ratio = ((-2 * self.M_simple ** 2) / (3 * self.Q_simple * self.K_simple ** 2)) T_ratio += (4 * self.N_simple / (3 * self.K_simple)) - T_ratio = float(T_ratio) - #print(f'S0 = {S_0}, Qs = {Q_simple}, Ks = {K_simple}, Ms = {M_simple}, Ns = {N_simple}, T_ratio = {T_ratio}') return T_ratio def get_num_density_centerline(self): @@ -766,10 +821,10 @@ def get_num_density_centerline(self): number density at a point on the centerline vs the number density at the nozzle exit ''' - p1 = self.X / np.sqrt(self.X ** 2 + self.R_0 ** 2) + # n_1(X, 0, 0) / n_0 [Cai & Wang 2012, Eq. 18] + p1 = self.X / np.sqrt(self.X ** 2 + self.R_0 ** 2) p2 = self.R_0 / np.sqrt(self.X ** 2 + self.R_0 ** 2) - n_ratio = 0.5 + 0.5 * sp.erf(self.S_0) - (p1 * np.exp(-self.S_0 ** 2 * p2 ** 2) / 2) * (1 + sp.erf(p1 * self.S_0)) - n_ratio = float(n_ratio) + n_ratio = 0.5 + 0.5 * erf(self.S_0) - (p1 * np.exp(-self.S_0 ** 2 * p2 ** 2) / 2) * (1 + erf(p1 * self.S_0)) return n_ratio @@ -791,11 +846,11 @@ def get_velocity_centerline(self): velocity at a point on the centerline (X, 0, 0) normalized with the parameter beta at the exit ''' - p1 = self.X / np.sqrt(self.X ** 2 + self.R_0 ** 2) + # U_1(X, 0, 0) * sqrt(beta_0) [Cai & Wang 2012, Eq. 19] + p1 = self.X / np.sqrt(self.X ** 2 + self.R_0 ** 2) p2 = self.R_0 / np.sqrt(self.X ** 2 + self.R_0 ** 2) n_ratio = self.get_num_density_centerline() - U_ratio = 1 / (2 * n_ratio) * ((p2 ** 2 * np.exp(- self.S_0 ** 2) / np.sqrt(np.pi)) + (self.S_0 * (1 + sp.erf(self.S_0))) - (np.exp(- p2 ** 2 * self.S_0 ** 2) * p1 ** 3 * self.S_0 * (1 + sp.erf(p1 * self.S_0)))) - U_ratio = float(U_ratio) + U_ratio = 1 / (2 * n_ratio) * ((p2 ** 2 * np.exp(- self.S_0 ** 2) / np.sqrt(np.pi)) + (self.S_0 * (1 + erf(self.S_0))) - (np.exp(- p2 ** 2 * self.S_0 ** 2) * p1 ** 3 * self.S_0 * (1 + erf(p1 * self.S_0)))) return U_ratio def get_temp_centerline(self): @@ -813,17 +868,12 @@ def get_temp_centerline(self): temperature on the point on the centerline (X, 0, 0) vs temperature at the nozzle exit ''' + # T_1(X, 0, 0) / T_0 [Cai & Wang 2012, Eq. 21] + # N_simple already carries the exp(-S_0^2) prefactor (see set_N_simple) n_ratio = self.get_num_density_centerline() U1 = self.get_velocity_centerline() - ''' - r = sp.symbols("r") - f = N * r - integral = sp.integrate(f, (r, 0, R_0)) - temp_ratio = (4 * np.exp(- S_0 ** 2)) / (3 * n_ratio * np.sqrt(np.pi) * X ** 2) * integral.evalf() - (U1 ** 2 / (3/2)) - ''' integral = 0.5 * self.N_simple * self.R_0 ** 2 - temp_ratio = (4 * np.exp(- self.S_0 ** 2)) / (3 * n_ratio * np.sqrt(np.pi) * self.X ** 2) * integral - (U1 ** 2 / (3/2)) - temp_ratio = float(temp_ratio) + temp_ratio = 4 / (3 * n_ratio * np.sqrt(np.pi) * self.X ** 2) * integral - (U1 ** 2 / (3/2)) return temp_ratio def get_pressure(self): @@ -946,191 +996,3 @@ def get_heat_flux(self): heat_flux = get_maxwellian_heat_transfer(rho_inf, S, self.sigma, self.theta, T, self.T_w, self.R, self.gamma) return heat_flux -''' -import math -from scipy.special import legendre - -class CollisionlessPlume: - - def __init__(self, U_0, R, T_0, n_iters, conv_tol): - self.n_iters = n_iters - self.conv_tol = conv_tol - self.U_0 = U_0 - self.R = R - self.T_0 = T_0 - self.S_0 = self.get_speed_ratio() - return - - def get_speed_ratio(self): - S_0 = self.U_0 / np.sqrt(2 * self.R * self.T_0) - return S_0 - - #solves a summation series of legengre polynomials of the first kind - #from degree 0 to degree n or until the convergance tolerance, tol is met - #n: int, x: float (value to evaluate function over), tol: float - def get_Q(self, X, Z, r, epsilon): - - # TEST TODO TEST TODO TEST - - psi = np.arctan(Z/X) - leg_x = np.sin(psi) * np.sin(epsilon) - - sum = 0 - - P_n_minus_2 = 1 #P_0 = 1 - P_n_minus_1 = leg_x #P_1 = x - sum += P_n_minus_2 + (P_n_minus_1 * r / np.sqrt(X ** 2 + Z ** 2)) - - - for degree in range (2, self.n_iters + 1, 1): - - P_n = legendre(degree)(leg_x) - #print(f'degree: {degree}; x: {leg_x}; P_n: {P_n}') - sum += P_n - if abs(P_n) < self.conv_tol and degree % 2 == 0: - break - - #my recurrence legendre solver - #solve for polynomial of degree of current iter - - P_n = (2 * (degree - 1) + 1) * leg_x * (P_n_minus_1) - P_n -= (degree - 1) * P_n_minus_2 - P_n /= degree - - sum += P_n * (r / np.sqrt(X ** 2 + Z ** 2)) ** degree - - P_n_minus_2 = P_n_minus_1 - P_n_minus_1 = P_n - - if abs(P_n) < self.conv_tol: - break - - sum = 0.712 - Q = (np.cos(psi) ** 2) * (sum ** 2) - - return Q - - def get_K(self, X, Z, r, epsilon): - #K = Q * ((Q * S_0) + ((0.5 + (Q * S_0 ** 2)) * sqrt(pi * Q) * - #(1 + erf(S_0 * sqrt(Q))) ** (Q * S_0 ** 2))) - Q = self.get_Q(X, Z, r, epsilon) - term1 = Q * self.S_0 - term2 = 0.5 + Q * self.S_0 ** 2 - term3 = np.sqrt(np.pi * Q) - term4 = (1 + sp.erf(self.S_0 * np.sqrt(Q))) * np.exp(Q * self.S_0 ** 2) - K = Q * (term1 + term2 * term3 * term4) - return K - - def get_M(self, X, Z, r, epsilon): - #M = (Q ** 2) * ((Q * S_0 ** 2) + 1 + (S_0 * (1.5 + (Q * S_0 ** 2)) * - #sqrt(pi * Q_simple)) * (1 + erf(S_0 * sqrt(Q))) ** (Q *S_0 ** 2)) - Q = self.get_Q(X, Z, r, epsilon) - term1 = Q * self.S_0 ** 2 - term2 = 1 + self.S_0 * (1.5 + Q * self.S_0 ** 2) - term3 = np.sqrt(np.pi * Q) - term4 = (1 + sp.erf(self.S_0 * np.sqrt(Q))) * np.exp(Q * self.S_0 ** 2) - M = Q ** 2 * (term1 + term2 * term3 * term4) - return M - - def get_N(self, X, Z, r, epsilon): - #N = S_0 * (Q ** 2) * (1.25 + (Q * S_0 ** 2) / 2) + - #(0.5 * sqrt(pi * Q ** 3)) * (0.75 + 3 * Q * S_0 **2 + Q ** 2 * S_0 ** 4) * - #(1 + erf(S_0 * sqrt(Q))) ** (Q * S_0 ** 2) - Q = self.get_Q(X, Z, r, epsilon) - term1 = self.S_0 * Q ** 2 * (1.25 + Q * self.S_0 ** 2 / 2) - term2 = 0.5 * np.sqrt(np.pi * Q ** 3) - term3 = 0.75 + 3 * Q * self.S_0 ** 2 + Q ** 2 * self.S_0 ** 4 - term4 = (1 + sp.erf(self.S_0 * np.sqrt(Q))) * np.exp(Q * self.S_0 ** 2) - N = term1 + term2 * term3 * term4 - return N - - def get_num_density_ratio(self, X, Z, R_0): - n_ratio = np.exp(-self.S_0 ** 2) / (X ** 2 * np.pi ** (3/2)) - n = 15 - e_a, e_b = -np.pi, np.pi - e_h = (e_b - e_a) / n - - e_vals = np.arange(e_a, e_b, e_h) - e_vals = np.append(e_vals, e_b) - dbl_integral = 0 - - epsilon = e_a - for idx, epsilon in enumerate(e_vals): - if epsilon == e_a or idx == len(e_vals) - 1: - integral = self.get_integral(X, Z, R_0, epsilon) - dbl_integral += integral - elif ((epsilon / e_h) % 3 == 0): - integral = self.get_integral(X, Z, R_0, epsilon) - dbl_integral += 2 * integral - else: - integral = self.get_integral(X, Z, R_0, epsilon) - dbl_integral += 3 * integral - dbl_integral *= 3 * e_h / 8 - n_ratio *= dbl_integral - print(f'S0: {self.S_0}; X: {X}; Z: {Z}; {n_ratio}') - return n_ratio - - def get_integral(self, X, Z, R_0, epsilon): - - r_a, r_b = 0, R_0 - n = 15 - r_h = (r_b - r_a) / n - - r_vals = np.arange(r_a, r_b, r_h) - r_vals = np.append(r_vals, r_b) - integral = 0 - r = r_a - for idx, r in enumerate(r_vals): - if r == r_a or idx == len(r_vals) - 1: - integral += r * self.get_K(X, Z, r, epsilon) - elif ((r / r_h) % 3 == 0): - integral += 2 * r * self.get_K(X, Z, r, epsilon) - else: - integral += 3 * r * self.get_K(X, Z, r, epsilon) - integral *= (3 * r_h / 8) - - return integral - -import numpy as np -import matplotlib.pyplot as plt - -# Assuming you have an instance of CollisionlessPlume named plume - -# Define the range of X / (2 * R_0) values -x_values = np.linspace(0.3, 10, 100) # Adjust the range as needed - -# Values of S_0 to plot -u0_values = [2] -R_0 = 5 -# Create subplots -plt.figure(figsize=(10, 6)) -plt.title('Normalized analytical number density along centerline') -plt.xlabel('X / D') -plt.ylabel('num_density_ratio') - -# Plot for each S_0 value -for u0 in u0_values: - num_density_ratios = [] - plume = CollisionlessPlume(u0 * 15.15255, 0.287, 400, 999, 0.00001) - # Calculate num_density_ratio for each X / (2 * R_0) - for x_ratio in x_values: - X = x_ratio * 2 * R_0 # Calculate X from the ratio - Z = 0 # Fixed Z value - num_density_ratio = plume.get_num_density_ratio(X, Z, R_0) - num_density_ratios.append(num_density_ratio) - - # Plot the results - plt.plot(x_values, num_density_ratios, label=f'S_0 = {u0}') - -# Show legend -plt.legend() - -# Show the plot -plt.show() - -#test = CollisionlessPlume.get_num_density_ratio() - #def get_U_normalized(): - #def get_W_normalized(): - #def get_temp_ratio(): - -''' \ No newline at end of file From c2c7e4abad3fbb945a64b7d62664dd091d20bb79 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Fri, 17 Jul 2026 22:52:57 -0500 Subject: [PATCH 03/14] fix simons model gamma dependence and limiting angle guard Co-Authored-By: Claude Fable 5 --- pyrpod/plume/RarefiedPlumeGasKinetics.py | 99 +++++++++++++---- tests/plume/test_simons.py | 136 +++++++++++++++++++++++ 2 files changed, 215 insertions(+), 20 deletions(-) create mode 100644 tests/plume/test_simons.py diff --git a/pyrpod/plume/RarefiedPlumeGasKinetics.py b/pyrpod/plume/RarefiedPlumeGasKinetics.py index d02282d..fb33ebb 100644 --- a/pyrpod/plume/RarefiedPlumeGasKinetics.py +++ b/pyrpod/plume/RarefiedPlumeGasKinetics.py @@ -283,7 +283,7 @@ class Simons: Number density from continuity equation with constant mass flux across different spherical surfaces. Return the ratio of number density at an analyzed point outside of the exit vs at the exit. ''' - def __init__(self, gamma, R, T_c, P_c, R_0, r): + def __init__(self, gamma, R, T_c, P_c, R_0, r, kappa=None): ''' Simple constructor, saves parameters to self. @@ -300,8 +300,15 @@ def __init__(self, gamma, R, T_c, P_c, R_0, r): R_0 : float nozzle exit radius (m) r : float - distance from the evaluated point + distance from the evaluated point to the nozzle exit center (m) + kappa : float, optional + plume beaming exponent for the cosine-law decay + function [Cai & Wang 2012, Eq. 26]. Defaults to + Boyton's kappa = 2/(gamma - 1) [Cai & Wang 2012, + Sec. II.C, refs. 22-23]. Paper-cited alternatives: + kappa = 2 (Ashkenas & Sherman, ref. 20) and + kappa = 1/(gamma - 1) (Albini, ref. 21). Returns ------- @@ -314,12 +321,11 @@ def __init__(self, gamma, R, T_c, P_c, R_0, r): self.P_c = P_c self.r = r self.R_0 = R_0 + if kappa is None: + kappa = 2 / (gamma - 1) #Boyton 1967/68 from Cai2012 [22][23] + self.kappa = kappa self.set_normalization_constant() - #self.rho_throat = self.get_nozzle_throat_density() - #self.theta_max = self.get_limiting_turn_angle() - #self.A = self.get_normalization_constant() - #self.U_t = self.get_limiting_velocity() def get_nozzle_throat_density(self): ''' @@ -336,8 +342,10 @@ def get_nozzle_throat_density(self): ''' #ideal gas law P_throat = rho_throat * R * T_throat #therefore: rho_throat = P_throat / (R * T_throat) - P_throat = self.P_c * 0.5283 #from Isentropic Flow Tables @ M = 1 - T_throat = self.T_c * 0.8333 #from Isentropic Flow Tables @ M = 1 + #isentropic ratios at M = 1: T*/T_c = 2/(gamma+1), + #P*/P_c = (2/(gamma+1))^(gamma/(gamma-1)) + T_throat = self.T_c * (2 / (self.gamma + 1)) + P_throat = self.P_c * (2 / (self.gamma + 1)) ** (self.gamma / (self.gamma - 1)) rho_throat = P_throat / (self.R * T_throat) return rho_throat @@ -358,22 +366,30 @@ def get_limiting_turn_angle(self): def get_plume_angular_density_decay_function(self, theta): ''' - From Cai 2012. Solve for the density decay function at a given off-centerline angle. - Expression for kappa is chosen from Boyton 1967/68. + Solve for the density decay function at a given off-centerline + angle [Cai & Wang 2012, Eq. 26], f(theta) = + cos^kappa(pi*theta/(2*theta_max)) with the beaming exponent + kappa chosen at construction (default Boyton 2/(gamma-1)). + + For theta >= theta_max the plume model region is empty (the + gas cannot turn past the limiting angle), so the decay + function is 0. Without this guard the cosine goes negative + and fractional kappa would produce complex numbers. Parameters ---------- theta : float plume centerline off-angle of current position (rad) - + Returns ------- float evaluation of plume angular density decay function at theta ''' - kappa = 2 / (self.gamma - 1) #Boyton 1967/68 from Cai2012 [22][23] theta_max = self.get_limiting_turn_angle() - f = (np.cos((np.pi / 2) * (theta / theta_max))) ** kappa + if theta >= theta_max: + return 0.0 + f = (np.cos((np.pi / 2) * (theta / theta_max))) ** self.kappa return f def set_normalization_constant(self): @@ -393,7 +409,7 @@ def set_normalization_constant(self): None. ''' theta_max = self.get_limiting_turn_angle() - kappa = 2 / (self.gamma - 1) + kappa = self.kappa def integrand(theta): return np.sin(theta) * np.cos((np.pi / 2) * (theta / theta_max)) ** kappa @@ -416,7 +432,7 @@ def get_sonic_velocity(self): float sonic velocity (m/s) ''' - T_throat = self.T_c * 0.8333 #from Isentropic Flow Tables @ Mach = 1 + T_throat = self.T_c * (2 / (self.gamma + 1)) #isentropic ratio at M = 1 sonic_velocity = np.sqrt(self.gamma * self.R * T_throat) return sonic_velocity @@ -450,25 +466,68 @@ def get_static_pressure(self, rho_ratio): def get_num_density_ratio(self, theta): ''' - Number density from continuity equation with constant mass flux across different spherical surfaces. - Return the ratio of number density at an analyzed point outside of the exit vs at the exit. + Number density from continuity equation with constant mass flux + across different spherical surfaces [Cai & Wang 2012, Eq. 25]. + + NOTE: this ratio is THROAT-referenced, n/n_s: the returned + value normalizes by the number density at the nozzle throat + (rho_s in the paper), not at the nozzle exit. For a ratio + comparable to the gas-kinetic classes (which normalize by the + exit density n_0), use get_num_density_ratio_exit. + + Returns 0.0 for theta >= theta_max (empty plume region). Parameters ---------- theta : float Angle off centerline of the current point being analyzed (rad). - + Returns ------- float - the ratio of number density at an analyzed point outside of the exit vs at the exit + number density at the analyzed point normalized by the + nozzle THROAT number density, n/n_s ''' + theta_max = self.get_limiting_turn_angle() + if theta >= theta_max: + return 0.0 f = self.get_plume_angular_density_decay_function(theta) - #??? make own function for rho_ratio??? rho_ratio = self.A * ((self.R_0/self.r) ** 2) * f #rho_ratio = density / nozzle throat denisty aka rho / rho_s n_ratio = rho_ratio return n_ratio + def get_num_density_ratio_exit(self, theta, exit_mach): + ''' + Exit-referenced number density ratio n/n_0. + + Rescales the throat-referenced Eq. 25 result by the isentropic + density ratio between the throat (M = 1) and the nozzle exit + (M = exit_mach): + + n_s / n_0 = [(1 + (gamma-1)/2 * M_e^2) / ((gamma+1)/2)]^(1/(gamma-1)) + + This makes the Simons model directly comparable to the + gas-kinetic classes (SimplifiedGasKinetics, + CollisionlessGasKinetics), which normalize by the exit + number density n_0. + + Parameters + ---------- + theta : float + Angle off centerline of the current point being analyzed (rad). + exit_mach : float + Mach number at the nozzle exit plane. + + Returns + ------- + float + number density at the analyzed point normalized by the + nozzle EXIT number density, n/n_0 + ''' + throat_to_exit = ((1 + (self.gamma - 1) / 2 * exit_mach ** 2) + / ((self.gamma + 1) / 2)) ** (1 / (self.gamma - 1)) + return self.get_num_density_ratio(theta) * throat_to_exit + # TODO save plume constants into self class SimplifiedGasKinetics: ''' diff --git a/tests/plume/test_simons.py b/tests/plume/test_simons.py new file mode 100644 index 0000000..445d1ba --- /dev/null +++ b/tests/plume/test_simons.py @@ -0,0 +1,136 @@ +# ======================== +# PyRPOD: tests/plume/test_simons.py +# ======================== +# Tests for the Simons cosine-law plume model (Cai & Wang 2012, Sec. II.C, +# Eqs. 25-26) after the gamma-generalization fixes: +# - isentropic throat ratios computed from gamma instead of the +# hardcoded gamma=1.4 table values 0.5283 / 0.8333, +# - theta >= theta_max returns an empty-region density of 0.0, +# - parameterizable beaming exponent kappa. + +import numpy as np +import pytest + +from pyrpod.plume.RarefiedPlumeGasKinetics import Simons + +pytestmark = pytest.mark.plume + +R_SPECIFIC = 287.0 # J / (kg K) +T_C = 500.0 # K +P_C = 1.0e5 # N / m^2 +R_0 = 0.1 # m +R_FIELD = 2.0 # m + + +def make_simons(gamma, kappa=None): + return Simons(gamma, R_SPECIFIC, T_C, P_C, R_0, R_FIELD, kappa=kappa) + + +def test_throat_ratios_reproduce_gamma_1_4_table_values(): + """The closed forms T*/T_c = 2/(gamma+1), P*/P_c = (2/(gamma+1))^ + (gamma/(gamma-1)) must reproduce the isentropic-table values + 0.8333 / 0.5283 for gamma = 1.4 to 4 decimal places.""" + gamma = 1.4 + t_ratio = 2 / (gamma + 1) + p_ratio = t_ratio ** (gamma / (gamma - 1)) + assert t_ratio == pytest.approx(0.8333, abs=5e-5) + assert p_ratio == pytest.approx(0.5283, abs=5e-5) + + # And the class must use them: density/sonic velocity at gamma=1.4 + # match the old hardcoded-table implementation to table precision. + simons = make_simons(gamma) + rho_old_table = (P_C * 0.5283) / (R_SPECIFIC * T_C * 0.8333) + assert simons.get_nozzle_throat_density() == pytest.approx(rho_old_table, + rel=1e-4) + a_old_table = np.sqrt(gamma * R_SPECIFIC * T_C * 0.8333) + assert simons.get_sonic_velocity() == pytest.approx(a_old_table, rel=1e-4) + + +def test_gamma_2_yields_real_finite_normalization(): + """Regression for the old 'broken when gamma = 2' failure: the + numerically integrated normalization constant is a plain real, + positive, finite float for fractional/even kappa alike.""" + simons = make_simons(2.0) + assert isinstance(simons.A, float) + assert np.isfinite(simons.A) + assert simons.A > 0.0 + + +@pytest.mark.parametrize("gamma", [1.4, 5.0 / 3.0, 2.0]) +def test_density_zero_beyond_limiting_angle(gamma): + """theta >= theta_max is a physically empty region: both the decay + function and the density ratio must return exactly 0.0 (previously + the cosine went negative and fractional kappa produced complex + numbers).""" + simons = make_simons(gamma) + theta_max = simons.get_limiting_turn_angle() + for theta in [theta_max, 1.01 * theta_max, 2.0 * theta_max]: + f = simons.get_plume_angular_density_decay_function(theta) + n = simons.get_num_density_ratio(theta) + assert f == 0.0 + assert n == 0.0 + assert not isinstance(f, complex) + assert not isinstance(n, complex) + + +def test_density_positive_inside_limiting_angle(): + simons = make_simons(2.0) + theta_max = simons.get_limiting_turn_angle() + n = simons.get_num_density_ratio(0.99 * theta_max) + assert np.isreal(n) + assert 0.0 < n < simons.get_num_density_ratio(0.0) + + +def test_default_kappa_is_boyton(): + """Default beaming exponent is Boyton's kappa = 2/(gamma-1) + [Cai & Wang 2012, Sec. II.C, refs. 22-23].""" + gamma = 1.4 + simons = make_simons(gamma) + assert simons.kappa == pytest.approx(2 / (gamma - 1)) + + +@pytest.mark.parametrize("kappa", [2.0, 2.5, 1 / 0.4]) +def test_kappa_parameterization(kappa): + """Custom kappa (e.g. Ashkenas & Sherman's 2, Albini's 1/(gamma-1)) + is used by the decay function: f = cos^kappa at the half angle, + and f(0) = 1 for every kappa (paper Sec. II.C).""" + simons = make_simons(1.4, kappa=kappa) + assert simons.kappa == kappa + assert simons.get_plume_angular_density_decay_function(0.0) == 1.0 + theta_max = simons.get_limiting_turn_angle() + theta = 0.5 * theta_max + expected = np.cos((np.pi / 2) * (theta / theta_max)) ** kappa + f = simons.get_plume_angular_density_decay_function(theta) + assert f == pytest.approx(expected, rel=1e-12) + + +def test_kappa_changes_normalization_constant(): + """A comes from mass-flow continuity over the chosen decay function, + so it must respond to kappa: a larger kappa concentrates the plume, + requiring a larger A.""" + a_kappa_2 = make_simons(1.4, kappa=2.0).A + a_kappa_5 = make_simons(1.4, kappa=5.0).A + assert a_kappa_5 > a_kappa_2 > 0.0 + + +def test_exit_referenced_density_at_mach_1_equals_throat_referenced(): + """With exit Mach 1 the exit IS the throat, so n/n_0 == n/n_s.""" + simons = make_simons(1.4) + theta = 0.3 + assert simons.get_num_density_ratio_exit(theta, 1.0) == pytest.approx( + simons.get_num_density_ratio(theta), rel=1e-12) + + +def test_exit_referenced_density_scaling(): + """n/n_0 = (n/n_s) * (n_s/n_0) with the isentropic throat-to-exit + density ratio [(1+(gamma-1)/2*Me^2)/((gamma+1)/2)]^(1/(gamma-1)). + The throat is denser than a supersonic exit, so the factor is > 1.""" + gamma = 1.4 + exit_mach = 5.0 + simons = make_simons(gamma) + theta = 0.3 + factor = ((1 + (gamma - 1) / 2 * exit_mach ** 2) + / ((gamma + 1) / 2)) ** (1 / (gamma - 1)) + assert factor > 1.0 + assert simons.get_num_density_ratio_exit(theta, exit_mach) == pytest.approx( + simons.get_num_density_ratio(theta) * factor, rel=1e-12) From 248b727e2f5cafd365c9f7d5bba2784eb3e12a43 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Fri, 17 Jul 2026 22:57:31 -0500 Subject: [PATCH 04/14] add full collisionless analytical plume model with closed-form q Co-Authored-By: Claude Fable 5 --- pyrpod/plume/RarefiedPlumeGasKinetics.py | 345 +++++++++++++++++++++++ tests/plume/test_q_closed_form.py | 74 +++++ 2 files changed, 419 insertions(+) create mode 100644 tests/plume/test_q_closed_form.py diff --git a/pyrpod/plume/RarefiedPlumeGasKinetics.py b/pyrpod/plume/RarefiedPlumeGasKinetics.py index fb33ebb..120b6a2 100644 --- a/pyrpod/plume/RarefiedPlumeGasKinetics.py +++ b/pyrpod/plume/RarefiedPlumeGasKinetics.py @@ -34,6 +34,8 @@ ' = simplified analytical results """ +import warnings + import numpy as np from scipy import integrate from scipy.special import erf @@ -120,6 +122,45 @@ def get_N_factor(Q, S_0): term3 = 0.75 + 3 * Q * S_0 ** 2 + Q ** 2 * S_0 ** 4 return term1 + term2 * term3 * erf_term + +def get_Q_full(r, epsilon, X, Z): + ''' + Full special factor Q [Cai & Wang 2012, Eq. 9] in closed form. + + Eq. 9 defines Q = cos^2(psi) * [sum_n P_n(sin(psi) sin(epsilon)) + * (r/sqrt(X^2+Z^2))^n]^2 with P_n the Legendre polynomials. The + series is the Legendre generating function sum_n P_n(x) t^n = + 1/sqrt(1 - 2*x*t + t^2) evaluated at x = sin(psi)sin(epsilon), + t = r/sqrt(X^2+Z^2), which collapses (with cos^2(psi) = + X^2/(X^2+Z^2)) to + + Q = X^2 / (X^2 + Z^2 - 2*Z*r*sin(epsilon) + r^2) + + so no series truncation or convergence loop is needed. The + denominator equals X^2 + (Z - r*sin(epsilon))^2 + + (r*cos(epsilon))^2 >= X^2 > 0, hence 0 < Q <= 1 always, which + also guarantees the overflow-safe exponential combination in + get_K_factor. On the centerline (Z = 0) this reduces to + Q(r) = X^2/(X^2 + r^2). + + Parameters + ---------- + r : float or ndarray + radial integration variable over the exit disk, 0 <= r <= R_0 (m) + epsilon : float or ndarray + angular integration variable, -pi/2 <= epsilon <= pi/2 (rad) + X : float + axial coordinate of the field point, X > 0 (m) + Z : float + transverse coordinate of the field point (m) + + Returns + ------- + float or ndarray + special factor Q at (r, epsilon) for field point (X, 0, Z) + ''' + return X ** 2 / (X ** 2 + Z ** 2 - 2 * Z * r * np.sin(epsilon) + r ** 2) + def get_maxwellian_pressure(rho_inf, U, S, sigma, theta, T, T_w): ''' Rarefied Gas Dynamics - Shen - eq. 4.19 @@ -1055,3 +1096,307 @@ def get_heat_flux(self): heat_flux = get_maxwellian_heat_transfer(rho_inf, S, self.sigma, self.theta, T, self.T_w, self.R, self.gamma) return heat_flux + + +class CollisionlessGasKinetics(SimplifiedGasKinetics): + ''' + Full collisionless analytical plume model [Cai & Wang 2012, + Sec. II.A, Eqs. 5-12]: a free jet expanding from a round exit + into vacuum, evaluated at a point (X, 0, Z) in front of the + nozzle (X > 0). Unlike the parent SimplifiedGasKinetics (which + substitutes the far-field Q' of Eq. 13), this class integrates + the exact special factor Q of Eq. 9 -- via its closed form, see + get_Q_full -- over the finite exit disk, so it remains valid in + the near field. + + The field solutions (Eqs. 5-8) are double integrals over + r in [0, R_0] and epsilon in [-pi/2, pi/2]. They are evaluated + with tensor-product Gauss-Legendre quadrature: the integrands + are analytic on the compact rectangle (the denominator of Q is + bounded below by X^2 > 0), so Gauss-Legendre converges + geometrically. The order is doubled (40 -> 80 -> 160) until two + successive orders agree to QUAD_RTOL; the finer result is kept. + In practice order 40 already reaches machine precision except + very near the nozzle lip (X -> 0, Z ~ R_0), where the density + field has a singularity [Cai & Wang 2012, Sec. III]. + + Centerline closed forms (Eqs. 18-21) and the Maxwell gas-surface + interface are inherited from SimplifiedGasKinetics; the exact + analytical solutions reduce to Eqs. 18-21 on the centerline, so + the inherited methods are exact there. The inherited surface + methods (get_pressure, get_shear_pressure, get_heat_flux) + dispatch to this class's overridden field getters, so they are + fed by the full-model n, U, W, T. + + Attributes + ---------- + (all of SimplifiedGasKinetics, plus) + + I_K : float + integral of r * exp(-S_0^2) * K over the exit disk [Eq. 5] + + I_M : float + integral of r * exp(-S_0^2) * M over the exit disk [Eq. 6] + + I_W : float + integral of (Z - r sin(epsilon)) * r * exp(-S_0^2) * M over + the exit disk [Eq. 7] + + I_N : float + integral of r * exp(-S_0^2) * N over the exit disk [Eq. 8] + + Methods + ------- + get_num_density_ratio() + + get_U_normalized() + + get_W_normalized() + + get_Vr_normalized() + + get_temp_ratio() + + get_pressure_ratio() + + (get_pressure / get_shear_pressure / get_heat_flux and the + centerline closed forms are inherited from SimplifiedGasKinetics) + ''' + + QUAD_ORDERS = (40, 80, 160) + QUAD_RTOL = 1e-9 + + def __init__(self, distance, theta, thruster_characteristics, T_w, sigma): + ''' + Mirrors SimplifiedGasKinetics(distance, theta, + thruster_characteristics, T_w, sigma) exactly; X = d*cos(theta) + and Z = d*sin(theta) are derived internally. Additionally + precomputes the four field integrals of Eqs. 5-8. + ''' + super().__init__(distance, theta, thruster_characteristics, T_w, sigma) + self.set_field_integrals() + + return + + def _compute_field_integrals(self, order): + ''' + Evaluate the four exit-disk integrals of Eqs. 5-8 with a + tensor-product Gauss-Legendre rule of the given order per + axis, over r in [0, R_0] and epsilon in [-pi/2, pi/2]. + + The K, M, N factors carry the exp(-S_0^2) prefactor of the + field solutions (see get_K_factor), keeping every term + bounded for large speed ratios. + + Parameters + ---------- + order : int + number of Gauss-Legendre nodes per axis + + Returns + ------- + tuple of float + (I_K, I_M, I_W, I_N) + ''' + nodes, weights = np.polynomial.legendre.leggauss(order) + # map [-1, 1] to [0, R_0] (radial) and [-pi/2, pi/2] (angular) + r = 0.5 * self.R_0 * (nodes + 1) + w_r = 0.5 * self.R_0 * weights + eps = 0.5 * np.pi * nodes + w_eps = 0.5 * np.pi * weights + + R, E = np.meshgrid(r, eps, indexing='ij') + W2D = np.outer(w_r, w_eps) + + Q = get_Q_full(R, E, self.X, self.Z) + K = get_K_factor(Q, self.S_0) + M = get_M_factor(Q, self.S_0) + N = get_N_factor(Q, self.S_0) + + I_K = np.sum(W2D * R * K) + I_M = np.sum(W2D * R * M) + I_W = np.sum(W2D * (self.Z - R * np.sin(E)) * R * M) + I_N = np.sum(W2D * R * N) + return I_K, I_M, I_W, I_N + + def set_field_integrals(self): + ''' + Setter for the field integrals I_K, I_M, I_W, I_N of + Eqs. 5-8, with quadrature-order doubling until two + successive orders agree to QUAD_RTOL (see class docstring + for the accuracy rationale). Warns if the finest order is + reached without convergence (only possible extremely close + to the nozzle-lip singularity). + + Parameters + ---------- + None. + + Returns + ------- + None. + ''' + previous = None + for order in self.QUAD_ORDERS: + current = self._compute_field_integrals(order) + if previous is not None and self._integrals_converged(previous, current): + break + previous = current + else: + warnings.warn( + 'CollisionlessGasKinetics quadrature did not converge to ' + 'rtol={} at order {} for point (X={}, Z={}); using finest ' + 'result.'.format(self.QUAD_RTOL, self.QUAD_ORDERS[-1], + self.X, self.Z), + RuntimeWarning) + + self.I_K, self.I_M, self.I_W, self.I_N = current + + return + + def _integrals_converged(self, previous, current): + ''' + Convergence test between two quadrature orders. I_K, I_M and + I_N are strictly positive, so a relative test applies; I_W + can legitimately vanish (centerline, Eq. 20), so its + difference is measured against the natural magnitude of its + integrand, (|Z| + R_0) * I_M. + + Parameters + ---------- + previous : tuple of float + integrals from the coarser rule + current : tuple of float + integrals from the finer rule + + Returns + ------- + bool + True when every integral has converged to QUAD_RTOL + ''' + rtol = self.QUAD_RTOL + I_K0, I_M0, I_W0, I_N0 = previous + I_K1, I_M1, I_W1, I_N1 = current + scale_W = (abs(self.Z) + self.R_0) * abs(I_M1) + return (abs(I_K1 - I_K0) <= rtol * abs(I_K1) + and abs(I_M1 - I_M0) <= rtol * abs(I_M1) + and abs(I_N1 - I_N0) <= rtol * abs(I_N1) + and abs(I_W1 - I_W0) <= rtol * scale_W) + + def get_num_density_ratio(self): + ''' + Number density at (X, 0, Z) normalized by the exit number + density, n_1/n_0 [Cai & Wang 2012, Eq. 5]. The exp(-S_0^2) + prefactor is carried inside I_K (see get_K_factor). + + Parameters + ---------- + None. + + Returns + ------- + float + n_1(X, 0, Z) / n_0 + ''' + return self.I_K / (np.pi ** 1.5 * self.X ** 2) + + def get_U_normalized(self): + ''' + Macroscopic x-velocity at (X, 0, Z) normalized by + sqrt(beta_0), i.e. U_1 * sqrt(beta_0) [Cai & Wang 2012, + Eq. 6]. The shared prefactor and n_0/n_1 reduce the ratio to + I_M / I_K. + + Parameters + ---------- + None. + + Returns + ------- + float + U_1(X, 0, Z) * sqrt(beta_0) + ''' + return self.I_M / self.I_K + + def get_W_normalized(self): + ''' + Macroscopic z-velocity at (X, 0, Z) normalized by + sqrt(beta_0), i.e. W_1 * sqrt(beta_0) [Cai & Wang 2012, + Eq. 7]. Eq. 7's integrand factor is read as + (Z - r sin(epsilon)); the printed "(Z - r sin(theta))" is a + typo -- sin(epsilon) is the convention consistent with Q + (Eq. 9) and it makes W vanish on the centerline as Eq. 20 + requires. Eq. 7's extra 1/X (prefactor 1/X^3 vs 1/X^2) + reduces the ratio to I_W / (X * I_K). + + Parameters + ---------- + None. + + Returns + ------- + float + W_1(X, 0, Z) * sqrt(beta_0) + ''' + return self.I_W / (self.X * self.I_K) + + def get_Vr_normalized(self): + ''' + Radial (spherical, from the nozzle exit center) velocity + component in the Y = 0 plane, normalized by sqrt(beta_0): + V_r = (X*U + Z*W)/sqrt(X^2 + Z^2) [Cai & Wang 2012, + Figs. 16-18]. + + Parameters + ---------- + None. + + Returns + ------- + float + V_r(X, 0, Z) * sqrt(beta_0) + ''' + U = self.get_U_normalized() + W = self.get_W_normalized() + return (self.X * U + self.Z * W) / np.sqrt(self.X ** 2 + self.Z ** 2) + + def get_temp_ratio(self): + ''' + Temperature at (X, 0, Z) normalized by the exit temperature, + T_1/T_0 [Cai & Wang 2012, Eq. 8]. With beta_0 = 1/(2*R*T_0), + the kinetic term -(U_1^2 + W_1^2)/(3*R*T_0) equals + -(2/3) * [(U_1 sqrt(beta_0))^2 + (W_1 sqrt(beta_0))^2]. + + Parameters + ---------- + None. + + Returns + ------- + float + T_1(X, 0, Z) / T_0 + ''' + U = self.get_U_normalized() + W = self.get_W_normalized() + return -(2 / 3) * (U ** 2 + W ** 2) + (4 / 3) * self.I_N / self.I_K + + def get_pressure_ratio(self): + ''' + Flowfield static pressure at (X, 0, Z) normalized by the exit + static pressure: p_1/p_0 = (n_1/n_0) * (T_1/T_0) from the + ideal gas law with the LOCAL temperature. Cai & Wang 2012 + (p. 64, Fig. 13 discussion) explicitly warn that computing + the local pressure as n(X, 0, Z) * k * T_0 is theoretically + invalid, because the flowfield temperature is always lower + than the exit temperature T_0. + + Parameters + ---------- + None. + + Returns + ------- + float + p_1(X, 0, Z) / p_0 + ''' + return self.get_num_density_ratio() * self.get_temp_ratio() diff --git a/tests/plume/test_q_closed_form.py b/tests/plume/test_q_closed_form.py new file mode 100644 index 0000000..c6505a2 --- /dev/null +++ b/tests/plume/test_q_closed_form.py @@ -0,0 +1,74 @@ +# ======================== +# PyRPOD: tests/plume/test_q_closed_form.py +# ======================== +# Verifies the closed form of the special factor Q (Cai & Wang 2012, +# Eq. 9) against a direct truncation of the printed Legendre series. +# +# Eq. 9: Q = cos^2(psi) * [sum_n P_n(sin(psi)sin(eps)) * t^n]^2 with +# t = r/sqrt(X^2+Z^2). The series is the Legendre generating function, +# so Q collapses to X^2/(X^2 + Z^2 - 2*Z*r*sin(eps) + r^2). The series +# converges like t^n, so with the test points restricted to t <= 0.6 a +# 50-term truncation is accurate to ~t^51 ~ 5e-12, and rel=1e-9 both +# passes robustly and would catch any error in the closed form. + +import numpy as np +import pytest +from scipy.special import eval_legendre + +from pyrpod.plume.RarefiedPlumeGasKinetics import get_Q_full + +pytestmark = pytest.mark.plume + +N_TERMS = 50 + + +def q_series(r, epsilon, X, Z, n_terms=N_TERMS): + """Truncated Eq. 9 exactly as printed in the paper.""" + psi = np.arctan2(Z, X) + x = np.sin(psi) * np.sin(epsilon) + t = r / np.sqrt(X ** 2 + Z ** 2) + series = sum(eval_legendre(n, x) * t ** n for n in range(n_terms)) + return np.cos(psi) ** 2 * series ** 2 + + +# (r, epsilon, X, Z) with t = r/sqrt(X^2+Z^2) <= 0.6 +SERIES_POINTS = [ + (0.05, 0.0, 1.0, 0.5), + (0.10, 0.7, 1.0, 0.5), + (0.10, -1.2, 1.0, 0.5), + (0.30, 1.5, 0.8, 0.6), + (0.60, -0.4, 1.0, 0.2), + (0.05, 0.3, 0.1, 0.05), + (0.50, 1.0, 2.0, 1.5), +] + + +@pytest.mark.parametrize("r,epsilon,X,Z", SERIES_POINTS) +def test_closed_form_matches_legendre_series(r, epsilon, X, Z): + """Closed form == 50-term truncation of the printed Eq. 9 series.""" + assert get_Q_full(r, epsilon, X, Z) == pytest.approx( + q_series(r, epsilon, X, Z), rel=1e-9) + + +def test_centerline_reduction(): + """On the centerline (Z = 0), Q reduces to X^2/(X^2 + r^2) for any + epsilon (needed by the Eq. 21 centerline temperature quadrature).""" + X = 0.7 + for r in [0.0, 0.05, 0.1]: + for epsilon in [-1.0, 0.0, 1.3]: + assert get_Q_full(r, epsilon, X, 0.0) == pytest.approx( + X ** 2 / (X ** 2 + r ** 2), rel=1e-14) + + +def test_q_bounded_between_zero_and_one(): + """0 < Q <= 1 everywhere in the integration domain: the denominator + is X^2 + (Z - r sin(eps))^2 + (r cos(eps))^2 >= X^2. This bound is + what makes the exp(-S_0^2*(1-Q)) combination overflow-safe, so it + is asserted over a grid that includes the worst case Z ~ r*sin(eps).""" + r = np.linspace(0.0, 0.1, 41) + epsilon = np.linspace(-np.pi / 2, np.pi / 2, 41) + R, E = np.meshgrid(r, epsilon) + for X, Z in [(1.0, 0.5), (0.05, 0.08), (0.01, 0.0), (0.3, 3.0)]: + Q = get_Q_full(R, E, X, Z) + assert np.all(Q > 0.0) + assert np.all(Q <= 1.0) From f21b75c03c6e2f9d461139cf8c10e4d6cc500183 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Fri, 17 Jul 2026 22:58:27 -0500 Subject: [PATCH 05/14] add far-field centerline asymptote functions Co-Authored-By: Claude Fable 5 --- pyrpod/plume/RarefiedPlumeGasKinetics.py | 60 ++++++++++++++++++++++ tests/plume/test_simplified_gaskinetics.py | 21 ++++++-- 2 files changed, 78 insertions(+), 3 deletions(-) diff --git a/pyrpod/plume/RarefiedPlumeGasKinetics.py b/pyrpod/plume/RarefiedPlumeGasKinetics.py index 120b6a2..366c20e 100644 --- a/pyrpod/plume/RarefiedPlumeGasKinetics.py +++ b/pyrpod/plume/RarefiedPlumeGasKinetics.py @@ -161,6 +161,66 @@ def get_Q_full(r, epsilon, X, Z): ''' return X ** 2 / (X ** 2 + Z ** 2 - 2 * Z * r * np.sin(epsilon) + r ** 2) + +def get_far_field_velocity_normalized(S_0): + ''' + Far-field centerline velocity asymptote, lim U_1 * sqrt(beta_0) + as X -> infinity [Cai & Wang 2012, Eqs. 22 and 24]. A function of + the exit speed ratio only: far enough from the exit the problem + degenerates to a small-hole effusion flow, so no nozzle geometry + factors remain. + + Evaluated by dividing numerator and denominator through by + [1 + erf(S_0)] * exp(S_0^2), which would itself overflow for + S_0 >~ 26; the residual factor exp(-S_0^2) / (1 + erf(S_0)) is + bounded for S_0 > 0. + + Parameters + ---------- + S_0 : float + molecular speed ratio at the nozzle exit, S_0 > 0 + + Returns + ------- + float + lim (X -> infinity) U_1(X, 0, 0) * sqrt(beta_0) + ''' + inv_E = np.exp(-S_0 ** 2) / (1 + erf(S_0)) + sqpi = np.sqrt(np.pi) + return S_0 + (inv_E + sqpi * S_0) / (S_0 * inv_E + (0.5 + S_0 ** 2) * sqpi) + + +def get_far_field_temp_ratio(S_0): + ''' + Far-field centerline temperature asymptote, lim T_1/T_0 as + X -> infinity [Cai & Wang 2012, Eqs. 23-24]. A function of the + exit speed ratio only (see get_far_field_velocity_normalized). + + Implemented as -2/3 * G^2 + 4*N(Q=1)/(3*K(Q=1)) with G of + Eq. 24, i.e. the Q -> 1 limit of Eq. 17; this reads Eq. 23's + printed denominator "3 S_0 + (1/2 + S_0^2) sqrt(pi) [...]" with + the 3 distributing over the whole denominator, 3*K(Q=1). The + equivalence was verified independently to 40 digits. Overflow + safety as in get_far_field_velocity_normalized. + + Parameters + ---------- + S_0 : float + molecular speed ratio at the nozzle exit, S_0 > 0 + + Returns + ------- + float + lim (X -> infinity) T_1(X, 0, 0) / T_0 + ''' + inv_E = np.exp(-S_0 ** 2) / (1 + erf(S_0)) + sqpi = np.sqrt(np.pi) + G = get_far_field_velocity_normalized(S_0) + num = S_0 * (5 + 2 * S_0 ** 2) * inv_E + 2 * sqpi * (0.75 + 3 * S_0 ** 2 + S_0 ** 4) + den = 3 * (S_0 * inv_E + (0.5 + S_0 ** 2) * sqpi) + return -(2 / 3) * G ** 2 + num / den + + def get_maxwellian_pressure(rho_inf, U, S, sigma, theta, T, T_w): ''' Rarefied Gas Dynamics - Shen - eq. 4.19 diff --git a/tests/plume/test_simplified_gaskinetics.py b/tests/plume/test_simplified_gaskinetics.py index 30c9ac0..0ab6aa6 100644 --- a/tests/plume/test_simplified_gaskinetics.py +++ b/tests/plume/test_simplified_gaskinetics.py @@ -15,7 +15,11 @@ import numpy as np import pytest -from pyrpod.plume.RarefiedPlumeGasKinetics import SimplifiedGasKinetics +from pyrpod.plume.RarefiedPlumeGasKinetics import ( + SimplifiedGasKinetics, + get_far_field_temp_ratio, + get_far_field_velocity_normalized, +) pytestmark = pytest.mark.plume @@ -107,6 +111,17 @@ def test_centerline_solution_pinned(point, expected): rel=RTOL_PINNED) +@pytest.mark.parametrize("S_0", [1.0, 2.0, 3.0]) +def test_far_field_asymptote_functions_pinned(S_0): + """Pin the Eq. 22-24 module-level asymptote functions against the + independent mpmath evaluation (tolerance rationale as RTOL_PINNED).""" + U_ref, T_ref = ASYMPTOTE_CASES[S_0] + assert get_far_field_velocity_normalized(S_0) == pytest.approx( + U_ref, rel=RTOL_PINNED) + assert get_far_field_temp_ratio(S_0) == pytest.approx( + T_ref, rel=RTOL_PINNED) + + @pytest.mark.parametrize("S_0", [1.0, 2.0, 3.0]) def test_centerline_velocity_approaches_asymptote(S_0): """Centerline U converges to the Eq. 22/24 far-field asymptote. @@ -115,7 +130,7 @@ def test_centerline_velocity_approaches_asymptote(S_0): relative, so rel=1e-5 passes with an order-of-magnitude margin while still requiring genuine convergence. """ - U_inf = ASYMPTOTE_CASES[S_0][0] + U_inf = get_far_field_velocity_normalized(S_0) X = 1000 * (D_NOZZLE / 2) plume = make_plume(X, 0.0, S_0) assert plume.get_velocity_centerline() == pytest.approx(U_inf, rel=1e-5) @@ -132,7 +147,7 @@ def test_centerline_temperature_approaches_asymptote(S_0): rel=1e-4 passes with margin for both the constant-N approximation and the exact Eq. 21 quadrature, while still requiring convergence. """ - T_inf = ASYMPTOTE_CASES[S_0][1] + T_inf = get_far_field_temp_ratio(S_0) X = 5000 * (D_NOZZLE / 2) plume = make_plume(X, 0.0, S_0) assert plume.get_temp_centerline() == pytest.approx(T_inf, rel=1e-4) From c9fac49b49e64c4eb43d6f7b269ca1ad955e454c Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Fri, 17 Jul 2026 22:59:44 -0500 Subject: [PATCH 06/14] fix centerline temperature integral in simplified model Co-Authored-By: Claude Fable 5 --- pyrpod/plume/RarefiedPlumeGasKinetics.py | 18 ++++++++-- tests/plume/test_simplified_gaskinetics.py | 38 ++++++++++++++++++++++ 2 files changed, 53 insertions(+), 3 deletions(-) diff --git a/pyrpod/plume/RarefiedPlumeGasKinetics.py b/pyrpod/plume/RarefiedPlumeGasKinetics.py index 366c20e..47be556 100644 --- a/pyrpod/plume/RarefiedPlumeGasKinetics.py +++ b/pyrpod/plume/RarefiedPlumeGasKinetics.py @@ -1018,6 +1018,13 @@ def get_temp_centerline(self): Method to calculate the temperature at a point (X, 0, 0) outside of the nozzle. This temperature is normalized over the temperature at the nozzle exit. + Evaluates Eq. 21's integral of N(Q(r)) * r over the exit + radius by numerical quadrature with the exact centerline + Q(r) = X^2/(X^2 + r^2) [Cai & Wang 2012, Eqs. 9, 21]. (An + earlier version treated N as constant across the exit, + integral = 0.5 * N(Q') * R_0^2, which is only valid in the + far field X >> R_0 where Q(r) ~ Q' = 1.) + Parameters ---------- None. @@ -1025,14 +1032,19 @@ def get_temp_centerline(self): Returns ------- float - temperature on the point on the centerline (X, 0, 0) + temperature on the point on the centerline (X, 0, 0) vs temperature at the nozzle exit ''' # T_1(X, 0, 0) / T_0 [Cai & Wang 2012, Eq. 21] - # N_simple already carries the exp(-S_0^2) prefactor (see set_N_simple) + # get_N_factor carries the exp(-S_0^2) prefactor of Eq. 21 n_ratio = self.get_num_density_centerline() U1 = self.get_velocity_centerline() - integral = 0.5 * self.N_simple * self.R_0 ** 2 + + def integrand(r): + Q = get_Q_full(r, 0.0, self.X, 0.0) + return get_N_factor(Q, self.S_0) * r + + integral, _ = integrate.quad(integrand, 0, self.R_0) temp_ratio = 4 / (3 * n_ratio * np.sqrt(np.pi) * self.X ** 2) * integral - (U1 ** 2 / (3/2)) return temp_ratio diff --git a/tests/plume/test_simplified_gaskinetics.py b/tests/plume/test_simplified_gaskinetics.py index 0ab6aa6..7eacfa4 100644 --- a/tests/plume/test_simplified_gaskinetics.py +++ b/tests/plume/test_simplified_gaskinetics.py @@ -153,6 +153,44 @@ def test_centerline_temperature_approaches_asymptote(S_0): assert plume.get_temp_centerline() == pytest.approx(T_inf, rel=1e-4) +def old_temp_centerline_approximation(plume): + """The pre-fix Eq. 21 evaluation: N treated as constant across the + exit radius, integral ~ 0.5 * N(Q'=1) * R_0^2 (N_simple carries the + exp(-S_0^2) prefactor).""" + integral = 0.5 * plume.N_simple * plume.R_0 ** 2 + return (4 / (3 * plume.get_num_density_centerline() * np.sqrt(np.pi) + * plume.X ** 2) * integral + - plume.get_velocity_centerline() ** 2 / 1.5) + + +def test_temp_centerline_quadrature_matches_old_in_far_field(): + """(a) Far field: with X >> R_0, Q(r) = X^2/(X^2+r^2) ~ 1 across the + whole exit, so the exact quadrature and the constant-N approximation + must agree. Measured relative difference at X/R_0 = 1000 is ~5e-5 + (it scales as (R_0/X)^2 amplified by the -2/3 U^2 cancellation, see + the asymptote test above); rel=1e-3 gives ~20x margin.""" + X = 1000 * (D_NOZZLE / 2) + plume = make_plume(X, 0.0, 2.0) + assert plume.get_temp_centerline() == pytest.approx( + old_temp_centerline_approximation(plume), rel=1e-3) + + +def test_temp_centerline_quadrature_diverges_from_old_in_near_field(): + """(b) Near field: at X = 2*R_0 the exact Q(r) varies strongly over + the exit (Q(R_0) = 0.8), and the old approximation errs by ~87%, + yielding an unphysical T > T_0 (collisionless expansion can only + cool the gas below the exit temperature -- paper Fig. 13 argument). + Require a divergence of at least 50% to prove the fix changes the + near field meaningfully.""" + X = 2 * (D_NOZZLE / 2) + plume = make_plume(X, 0.0, 2.0) + T_new = plume.get_temp_centerline() + T_old = old_temp_centerline_approximation(plume) + assert abs(T_new / T_old - 1) > 0.5 + assert T_new < 1.0 + assert T_old > 1.0 + + @pytest.mark.parametrize("point,expected", [ ((1.0, 0.2, 2.0), (0.36855461276327334, -0.056815859301478394, 143.3923778236429)), From ed9b0395f8aeac026e019c080d07f99af5e177ff Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Fri, 17 Jul 2026 23:01:16 -0500 Subject: [PATCH 07/14] add verification tests for full analytical plume model Co-Authored-By: Claude Fable 5 --- tests/plume/test_collisionless_gaskinetics.py | 166 ++++++++++++++++++ 1 file changed, 166 insertions(+) create mode 100644 tests/plume/test_collisionless_gaskinetics.py diff --git a/tests/plume/test_collisionless_gaskinetics.py b/tests/plume/test_collisionless_gaskinetics.py new file mode 100644 index 0000000..6086190 --- /dev/null +++ b/tests/plume/test_collisionless_gaskinetics.py @@ -0,0 +1,166 @@ +# ======================== +# PyRPOD: tests/plume/test_collisionless_gaskinetics.py +# ======================== +# Verification tests for the full collisionless analytical plume model +# CollisionlessGasKinetics (Cai & Wang 2012, Eqs. 5-12), checked against +# the paper's own internal anchors: +# (a, b) the exact field integrals must reduce to the Eq. 18/19 +# closed forms on the centerline, +# (c) W must vanish on the centerline (Eq. 20), +# (d) the full and simplified models must agree in the far field +# within a bound derived from Eq. 30, +# (e) qualitative physics of Figs. 3, 5, 13 (monotonic centerline +# density decay, cooling below T_0), +# (f) convergence to the Eq. 22-24 far-field asymptotes. + +import numpy as np +import pytest + +from pyrpod.plume.RarefiedPlumeGasKinetics import ( + CollisionlessGasKinetics, + SimplifiedGasKinetics, + get_far_field_temp_ratio, + get_far_field_velocity_normalized, +) + +pytestmark = pytest.mark.plume + +R_SPECIFIC = 208.13 # J / (kg K), argon-like +T_0 = 500.0 # K +GAMMA = 5.0 / 3.0 +N_0 = 1.0e20 # m^-3 +D_NOZZLE = 0.2 # m (R_0 = 0.1 m) +R_0 = D_NOZZLE / 2 +T_W = 300.0 # K +SIGMA = 1.0 + +# The class doubles the quadrature order until successive orders agree +# to QUAD_RTOL = 1e-9; comparisons against exact closed forms therefore +# use a slightly looser 1e-8. +RTOL_QUAD = 1e-8 + + +def make_plume(cls, distance, theta, S_0): + ve = S_0 * np.sqrt(2 * R_SPECIFIC * T_0) + thruster_characteristics = { + 'd': D_NOZZLE, 've': ve, 'R': R_SPECIFIC, + 'gamma': GAMMA, 'Te': T_0, 'n': N_0, + } + return cls(distance, theta, thruster_characteristics, T_W, SIGMA) + + +CENTERLINE_POINTS = [(0.2, 1.0), (0.5, 2.0), (1.0, 2.0), (5.0, 2.0), + (1.0, 3.0)] + + +@pytest.mark.parametrize("X,S_0", CENTERLINE_POINTS) +def test_centerline_density_matches_eq18(X, S_0): + """(a) The Eq. 5 double integral evaluated at (X, 0, 0) must equal + the Eq. 18 closed form (they were derived from the same solution).""" + plume = make_plume(CollisionlessGasKinetics, X, 0.0, S_0) + assert plume.get_num_density_ratio() == pytest.approx( + plume.get_num_density_centerline(), rel=RTOL_QUAD) + + +@pytest.mark.parametrize("X,S_0", CENTERLINE_POINTS) +def test_centerline_velocity_matches_eq19(X, S_0): + """(b) The Eq. 6 integral at (X, 0, 0) must equal the Eq. 19 closed + form.""" + plume = make_plume(CollisionlessGasKinetics, X, 0.0, S_0) + assert plume.get_U_normalized() == pytest.approx( + plume.get_velocity_centerline(), rel=RTOL_QUAD) + + +@pytest.mark.parametrize("X,S_0", CENTERLINE_POINTS) +def test_centerline_W_vanishes(X, S_0): + """(c) Eq. 20: W(X, 0, 0) = 0. The quadrature can only leave + rounding noise, so |W| is required to be negligible against U. On + the centerline V_r must also reduce to U.""" + plume = make_plume(CollisionlessGasKinetics, X, 0.0, S_0) + U = plume.get_U_normalized() + assert abs(plume.get_W_normalized()) < 1e-12 * U + assert plume.get_Vr_normalized() == pytest.approx(U, rel=1e-12) + + +# Far-field points (distance, theta) with S_0 = 2. +FAR_FIELD_POINTS = [(20.0, 0.4), (50.0, 0.8), (100.0, 0.2)] + + +@pytest.mark.parametrize("distance,theta", FAR_FIELD_POINTS) +def test_far_field_agrees_with_simplified_within_eq30_bound(distance, theta): + """(d) Eq. 30 bounds the pointwise error of the simplified Q': + |Q' - Q|/Q <= max(R_0^2, 2*Z*R_0)/(X^2 + Z^2) over the exit disk. + The K, M, N factors depend on Q through powers up to Q^(5/2), the + factor exp(S_0^2 * Q), and bounded erf terms, so their logarithmic + sensitivity |d ln K / d ln Q| is at most ~(S_0^2 + 4). The relative + difference between the full and simplified field quantities is + therefore bounded by (S_0^2 + 4) * max(R_0^2, 2*Z*R_0)/(X^2 + Z^2).""" + S_0 = 2.0 + full = make_plume(CollisionlessGasKinetics, distance, theta, S_0) + simple = make_plume(SimplifiedGasKinetics, distance, theta, S_0) + + bound = ((S_0 ** 2 + 4) + * max(R_0 ** 2, 2 * full.Z * R_0) / (full.X ** 2 + full.Z ** 2)) + + assert full.get_num_density_ratio() == pytest.approx( + simple.get_num_density_ratio(), rel=bound) + assert full.get_U_normalized() == pytest.approx( + simple.get_U_normalized(), rel=bound) + assert full.get_W_normalized() == pytest.approx( + simple.get_W_normalized(), rel=bound) + assert full.get_temp_ratio() == pytest.approx( + simple.get_temp_ratio(), rel=bound) + + +def test_centerline_density_decays_monotonically(): + """(e) Paper Fig. 3: for fixed S_0 the centerline density decreases + monotonically with X.""" + S_0 = 2.0 + X_values = [0.15, 0.3, 0.6, 1.2, 2.5, 5.0, 10.0, 20.0] + densities = [make_plume(CollisionlessGasKinetics, X, 0.0, S_0) + .get_num_density_ratio() for X in X_values] + assert all(a > b for a, b in zip(densities, densities[1:])) + assert all(n > 0 for n in densities) + + +@pytest.mark.parametrize("distance,theta", [(0.2, 0.0), (1.0, 0.0), + (1.0, 0.5), (5.0, 1.0)]) +def test_temperature_below_exit_temperature_downstream(distance, theta): + """(e) Paper Figs. 5, 13: the collisionless expansion only cools the + gas, so T_1 < T_0 everywhere downstream; consequently the static + pressure ratio p_1/p_0 = (n/n_0)(T/T_0) lies below the density + ratio (the paper's p. 64 argument against using n*k*T_0).""" + plume = make_plume(CollisionlessGasKinetics, distance, theta, 2.0) + T_ratio = plume.get_temp_ratio() + assert 0.0 < T_ratio < 1.0 + assert 0.0 < plume.get_pressure_ratio() < plume.get_num_density_ratio() + + +@pytest.mark.parametrize("S_0", [1.0, 2.0, 3.0]) +def test_centerline_converges_to_far_field_asymptotes(S_0): + """(f) Full-model centerline U and T converge to the Eq. 22-24 + asymptotes at large X. Tolerances follow the O((R_0/X)^2) + convergence-rate arguments in test_simplified_gaskinetics.py + (rel=1e-5 for U at X/R_0 = 1000; rel=1e-4 for T at X/R_0 = 5000 + due to the -2/3 U^2 cancellation amplification).""" + plume_U = make_plume(CollisionlessGasKinetics, 1000 * R_0, 0.0, S_0) + assert plume_U.get_U_normalized() == pytest.approx( + get_far_field_velocity_normalized(S_0), rel=1e-5) + + plume_T = make_plume(CollisionlessGasKinetics, 5000 * R_0, 0.0, S_0) + assert plume_T.get_temp_ratio() == pytest.approx( + get_far_field_temp_ratio(S_0), rel=1e-4) + + +def test_surface_interaction_uses_full_model_and_is_finite(): + """The inherited Maxwell gas-surface interface dispatches to the + overridden full-model field getters; off the centerline it must + return finite, physically-signed values (p > 0, q finite, shear + opposing the flow direction).""" + plume = make_plume(CollisionlessGasKinetics, 1.0, 0.4, 2.0) + p = float(plume.get_pressure()) + tau = float(plume.get_shear_pressure()) + q = float(plume.get_heat_flux()) + assert np.isfinite(p) and p > 0.0 + assert np.isfinite(tau) and tau < 0.0 + assert np.isfinite(q) From 8f56cfaf0304fb811039c8af339630154ce56554 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Fri, 17 Jul 2026 23:27:09 -0500 Subject: [PATCH 08/14] rename new plume tests to follow subsystem naming convention Co-Authored-By: Claude Fable 5 --- .../plume/{test_q_closed_form.py => plume_unit_test_02.py} | 4 ++-- tests/plume/{test_simons.py => plume_unit_test_03.py} | 4 ++-- ...plified_gaskinetics.py => plume_verification_test_02.py} | 4 ++-- ...ionless_gaskinetics.py => plume_verification_test_03.py} | 6 +++--- 4 files changed, 9 insertions(+), 9 deletions(-) rename tests/plume/{test_q_closed_form.py => plume_unit_test_02.py} (96%) rename tests/plume/{test_simons.py => plume_unit_test_03.py} (98%) rename tests/plume/{test_simplified_gaskinetics.py => plume_verification_test_02.py} (98%) rename tests/plume/{test_collisionless_gaskinetics.py => plume_verification_test_03.py} (97%) diff --git a/tests/plume/test_q_closed_form.py b/tests/plume/plume_unit_test_02.py similarity index 96% rename from tests/plume/test_q_closed_form.py rename to tests/plume/plume_unit_test_02.py index c6505a2..8022e3a 100644 --- a/tests/plume/test_q_closed_form.py +++ b/tests/plume/plume_unit_test_02.py @@ -1,5 +1,5 @@ # ======================== -# PyRPOD: tests/plume/test_q_closed_form.py +# PyRPOD: tests/plume/plume_unit_test_02.py # ======================== # Verifies the closed form of the special factor Q (Cai & Wang 2012, # Eq. 9) against a direct truncation of the printed Legendre series. @@ -17,7 +17,7 @@ from pyrpod.plume.RarefiedPlumeGasKinetics import get_Q_full -pytestmark = pytest.mark.plume +pytestmark = [pytest.mark.plume, pytest.mark.unit] N_TERMS = 50 diff --git a/tests/plume/test_simons.py b/tests/plume/plume_unit_test_03.py similarity index 98% rename from tests/plume/test_simons.py rename to tests/plume/plume_unit_test_03.py index 445d1ba..a9d48de 100644 --- a/tests/plume/test_simons.py +++ b/tests/plume/plume_unit_test_03.py @@ -1,5 +1,5 @@ # ======================== -# PyRPOD: tests/plume/test_simons.py +# PyRPOD: tests/plume/plume_unit_test_03.py # ======================== # Tests for the Simons cosine-law plume model (Cai & Wang 2012, Sec. II.C, # Eqs. 25-26) after the gamma-generalization fixes: @@ -13,7 +13,7 @@ from pyrpod.plume.RarefiedPlumeGasKinetics import Simons -pytestmark = pytest.mark.plume +pytestmark = [pytest.mark.plume, pytest.mark.unit] R_SPECIFIC = 287.0 # J / (kg K) T_C = 500.0 # K diff --git a/tests/plume/test_simplified_gaskinetics.py b/tests/plume/plume_verification_test_02.py similarity index 98% rename from tests/plume/test_simplified_gaskinetics.py rename to tests/plume/plume_verification_test_02.py index 7eacfa4..54a7d61 100644 --- a/tests/plume/test_simplified_gaskinetics.py +++ b/tests/plume/plume_verification_test_02.py @@ -1,5 +1,5 @@ # ======================== -# PyRPOD: tests/plume/test_simplified_gaskinetics.py +# PyRPOD: tests/plume/plume_verification_test_02.py # ======================== # Pinning tests for the SimplifiedGasKinetics class (Cai & Wang 2012, # "Numerical Validations for a Set of Collisionless Rocket Plume Solutions", @@ -21,7 +21,7 @@ get_far_field_velocity_normalized, ) -pytestmark = pytest.mark.plume +pytestmark = [pytest.mark.plume, pytest.mark.verification] # Fixed thruster gas properties used for all cases (argon-like). R_SPECIFIC = 208.13 # J / (kg K) diff --git a/tests/plume/test_collisionless_gaskinetics.py b/tests/plume/plume_verification_test_03.py similarity index 97% rename from tests/plume/test_collisionless_gaskinetics.py rename to tests/plume/plume_verification_test_03.py index 6086190..fa3f5ba 100644 --- a/tests/plume/test_collisionless_gaskinetics.py +++ b/tests/plume/plume_verification_test_03.py @@ -1,5 +1,5 @@ # ======================== -# PyRPOD: tests/plume/test_collisionless_gaskinetics.py +# PyRPOD: tests/plume/plume_verification_test_03.py # ======================== # Verification tests for the full collisionless analytical plume model # CollisionlessGasKinetics (Cai & Wang 2012, Eqs. 5-12), checked against @@ -23,7 +23,7 @@ get_far_field_velocity_normalized, ) -pytestmark = pytest.mark.plume +pytestmark = [pytest.mark.plume, pytest.mark.verification] R_SPECIFIC = 208.13 # J / (kg K), argon-like T_0 = 500.0 # K @@ -140,7 +140,7 @@ def test_temperature_below_exit_temperature_downstream(distance, theta): def test_centerline_converges_to_far_field_asymptotes(S_0): """(f) Full-model centerline U and T converge to the Eq. 22-24 asymptotes at large X. Tolerances follow the O((R_0/X)^2) - convergence-rate arguments in test_simplified_gaskinetics.py + convergence-rate arguments in plume_verification_test_02.py (rel=1e-5 for U at X/R_0 = 1000; rel=1e-4 for T at X/R_0 = 5000 due to the -2/3 U^2 cancellation amplification).""" plume_U = make_plume(CollisionlessGasKinetics, 1000 * R_0, 0.0, S_0) From bf1387e83863ec1f4ce001341d8d2709e086ecfe Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Sat, 18 Jul 2026 07:48:57 -0500 Subject: [PATCH 09/14] add shared utilities for plume verification figures Co-Authored-By: Claude Fable 5 --- tests/plume/data/digitized/README.md | 29 ++ tests/plume/output/.gitignore | 2 + tests/plume/plume_figure_utils.py | 519 +++++++++++++++++++++++++++ 3 files changed, 550 insertions(+) create mode 100644 tests/plume/data/digitized/README.md create mode 100644 tests/plume/output/.gitignore create mode 100644 tests/plume/plume_figure_utils.py diff --git a/tests/plume/data/digitized/README.md b/tests/plume/data/digitized/README.md new file mode 100644 index 0000000..9f43dd1 --- /dev/null +++ b/tests/plume/data/digitized/README.md @@ -0,0 +1,29 @@ +# Digitized reference data for the plume verification figures + +Drop digitized curves from Cai & Wang 2012 (JSR 49(1), DOI +10.2514/1.A32046) here and the manual-run figure scripts +(`tests/plume/plume_verification_test_04` ... `_27`) will overlay them +automatically on their next run. No code changes are needed. + +## File convention + +- Name: `_.csv`, where `` is the paper figure + (e.g. `fig19`) and `` names the dataset. Examples: + - `fig03_dsmc.csv` — DSMC centerline density (S0 = 2), Fig. 3 + - `fig19_dsmc.csv` — DSMC centerline density, Kn = 100, Fig. 19 + - `fig25_dsmc_kn100.csv`, `fig25_dsmc_kn0p1.csv`, + `fig25_dsmc_kn0p01.csv` — the three DSMC mass-flux curves, Fig. 25 + - contour figures: one CSV per digitized contour polyline, e.g. + `fig06_dsmc_0p001.csv` (the DSMC n/n0 = 0.001 line of Fig. 6) +- Content: two comma-separated columns `x,y` with one header row. + x/y are in the figure's plotted units (X/D, Z/D, theta in degrees, + normalized quantity values). +- Series suffixes are prettified for legends: `dsmc` -> `DSMC`, + `kn0p1` -> `Kn=0.1`. +- Lower-half-plane series in the split contour figures (Figs. 6-8, + 11-18) should be digitized with positive Z/D values taken from the + paper's lower half; enter them as printed (negative y) — overlays are + drawn exactly as given. + +Each figure script's docstring lists the exact filenames it looks for. +Files that are absent are silently skipped. diff --git a/tests/plume/output/.gitignore b/tests/plume/output/.gitignore new file mode 100644 index 0000000..d6b7ef3 --- /dev/null +++ b/tests/plume/output/.gitignore @@ -0,0 +1,2 @@ +* +!.gitignore diff --git a/tests/plume/plume_figure_utils.py b/tests/plume/plume_figure_utils.py new file mode 100644 index 0000000..eae7606 --- /dev/null +++ b/tests/plume/plume_figure_utils.py @@ -0,0 +1,519 @@ +# ======================== +# PyRPOD: tests/plume/plume_figure_utils.py +# ======================== +# Shared helpers for the manual-run verification figure scripts +# plume_verification_test_04 ... _27, which reproduce Figs. 2-25 of +# Cai, C. and Wang, L., "Numerical Validations for a Set of Collisionless +# Rocket Plume Solutions," JSR 49(1), 2012, DOI 10.2514/1.A32046. +# +# Design notes +# ------------ +# * All figures use the paper's validation conditions: argon, +# D = 0.2 m, gamma = 5/3, R = 208.13 J/(kg K). Every plotted quantity +# is normalized (n/n0, U*sqrt(beta0), T/T0, ...), so the absolute +# T_0 and n_0 values are inert; the ones below match the existing +# verification tests. ve is set from the exit speed ratio, +# ve = S_0 * sqrt(2*R*T_0). +# * Full-model contour fields are evaluated with a vectorized +# Gauss-Legendre quadrature (fixed order 80 per axis, the same order +# at which CollisionlessGasKinetics converges to 1e-9). Every call +# cross-checks a few sample points against the class and raises if +# they disagree beyond 1e-6, so the vectorized path cannot silently +# diverge from the physics module. +# * The models require X > 0; grids start at X_MIN_OVER_D and angular +# sweeps stop at THETA_MAX_DEG < 90 deg. +# * Digitized-data convention (see tests/plume/data/digitized/README.md): +# any CSV named '_.csv' (two columns x,y, one header +# row) is overlaid automatically by overlay_digitized(); absent files +# are silently skipped, so the DSMC series appear once the paper's +# curves are digitized. +# * Simons overlays are EXIT-referenced (n/n_0) via the isentropic +# throat-to-exit conversion with Me = S_0*sqrt(2/gamma) (from +# S_0 = U0/sqrt(2RT0) and M = U0/sqrt(gamma*R*T0)). Following the +# paper's remark that the cosine-law curves with different kappa +# coincide on the centerline (p. 65: "the cosine-law model with +# different kappa values actually provides the same value"), the +# normalization constant A is always taken from the Boyton default +# kappa while the plotted decay exponent may vary (Figs. 22-24). + +import sys +import time +from pathlib import Path + +import numpy as np +import matplotlib +matplotlib.use('Agg') # save-only figures; must precede pyplot import +import matplotlib.pyplot as plt # noqa: E402 (backend must be set first) + +_THIS_DIR = Path(__file__).resolve().parent +_REPO_ROOT = _THIS_DIR.parents[1] +if str(_REPO_ROOT) not in sys.path: # direct `python tests/plume/...py` runs + sys.path.insert(0, str(_REPO_ROOT)) + +from pyrpod.plume.RarefiedPlumeGasKinetics import ( # noqa: E402 + CollisionlessGasKinetics, + SimplifiedGasKinetics, + Simons, + get_K_factor, + get_M_factor, + get_N_factor, + get_Q_full, +) + +# Paper validation conditions (Sec. III: argon, D = 0.2 m) +R_SPECIFIC = 208.13 # J / (kg K) +T_0 = 500.0 # K (inert for normalized results) +GAMMA = 5.0 / 3.0 +N_0 = 1.0e20 # m^-3 (inert for normalized results) +D_NOZZLE = 0.2 # m +R_0 = D_NOZZLE / 2 # m +T_W = 300.0 # K (unused by plotted quantities) +SIGMA = 1.0 +BOLTZMANN = 1.380649e-23 # J / K + +OUTPUT_DIR = _THIS_DIR / 'output' +DIGITIZED_DIR = _THIS_DIR / 'data' / 'digitized' + +X_MIN_OVER_D = 0.05 # models require X > 0 +THETA_MAX_DEG = 89.5 # angular sweeps stop short of 90 deg (X -> 0) +QUAD_ORDER = 80 # per-axis Gauss-Legendre order for grid fields +GRID_N = 161 # nodes along the longer contour-grid axis + +_GL_NODES, _GL_WEIGHTS = np.polynomial.legendre.leggauss(QUAD_ORDER) + + +# --------------------------------------------------------------------------- +# model factories +# --------------------------------------------------------------------------- + +def thruster_characteristics(S_0): + """Thruster dict for the gas-kinetic classes at exit speed ratio S_0.""" + ve = S_0 * np.sqrt(2 * R_SPECIFIC * T_0) + return {'d': D_NOZZLE, 've': ve, 'R': R_SPECIFIC, + 'gamma': GAMMA, 'Te': T_0, 'n': N_0} + + +def make_plume(cls, distance, theta, S_0): + """Instantiate SimplifiedGasKinetics or CollisionlessGasKinetics.""" + return cls(distance, theta, thruster_characteristics(S_0), T_W, SIGMA) + + +def exit_mach(S_0): + """Me from S_0: S_0 = U0/sqrt(2RT0), M = U0/sqrt(gamma R T0).""" + return S_0 * np.sqrt(2.0 / GAMMA) + + +def throat_to_exit_density_ratio(S_0): + """n_s/n_0 used by Simons exit referencing [module docstring note].""" + Me = exit_mach(S_0) + return ((1 + (GAMMA - 1) / 2 * Me ** 2) + / ((GAMMA + 1) / 2)) ** (1 / (GAMMA - 1)) + + +def make_simons(distance, S_0, kappa=None): + """Simons instance with chamber conditions isentropically consistent + with the exit state at Me(S_0). T_c/P_c only rescale absolute + density; they are inert for the plotted ratios.""" + Me = exit_mach(S_0) + stag = 1 + (GAMMA - 1) / 2 * Me ** 2 + T_c = T_0 * stag + P_c = (N_0 * BOLTZMANN * T_0) * stag ** (GAMMA / (GAMMA - 1)) + return Simons(GAMMA, R_SPECIFIC, T_c, P_c, R_0, distance, kappa=kappa) + + +def simons_density_field(X, Z, S_0, kappa_f=None): + """Exit-referenced Simons density n/n_0 at points (X, 0, Z) [m]. + + A is always the Boyton-kappa normalization constant (paper p. 65: + the cosine-law curves coincide at theta = 0 for all kappa); kappa_f + optionally varies the plotted decay exponent f = cos^kappa_f + (Figs. 22-24). Vectorized; theta >= theta_max gives 0. + """ + X = np.asarray(X, dtype=float) + Z = np.asarray(Z, dtype=float) + boyton = make_simons(1.0, S_0) # A from default Boyton kappa + theta_max = boyton.get_limiting_turn_angle() + kappa = boyton.kappa if kappa_f is None else kappa_f + + r_sph = np.sqrt(X ** 2 + Z ** 2) + theta = np.arctan2(np.abs(Z), X) + f = np.where(theta < theta_max, + np.cos((np.pi / 2) * (np.minimum(theta, theta_max) + / theta_max)) ** kappa, + 0.0) + n_over_ns = boyton.A * (R_0 / r_sph) ** 2 * f + return n_over_ns * throat_to_exit_density_ratio(S_0) + + +# --------------------------------------------------------------------------- +# centerline profiles (exact closed forms, Eqs. 18/19/21) +# --------------------------------------------------------------------------- + +def centerline_profiles(x_over_D, S_0, quantities=('n',)): + """Exact analytical centerline values at X/D array; quantities from + {'n', 'U', 'T'} -> Eqs. 18 / 19 / 21.""" + getters = {'n': 'get_num_density_centerline', + 'U': 'get_velocity_centerline', + 'T': 'get_temp_centerline'} + out = {q: np.empty(len(x_over_D)) for q in quantities} + for i, xd in enumerate(x_over_D): + plume = make_plume(SimplifiedGasKinetics, xd * D_NOZZLE, 0.0, S_0) + for q in quantities: + out[q][i] = getattr(plume, getters[q])() + return out + + +def simplified_centerline_density(x_over_D, S_0): + """Simplified-model centerline density, Eq. 14 with Q' = 1.""" + X = np.asarray(x_over_D) * D_NOZZLE + return get_K_factor(1.0, S_0) / (2 * np.sqrt(np.pi)) * (R_0 / X) ** 2 + + +# --------------------------------------------------------------------------- +# vectorized full-model field evaluation (Eqs. 5-8) +# --------------------------------------------------------------------------- + +def full_field_values(quantity, X, Z, S_0, _chunk=256): + """Full-model field quantity at 1-D point arrays X, Z [m]. + + quantity: 'n', 'U', 'W', 'Vr', 'T', or 'p'. Same integrals as + CollisionlessGasKinetics at its converged order (80); each call + verifies sample points against the class to 1e-6 (see module + docstring). + """ + X = np.asarray(X, dtype=float).ravel() + Z = np.asarray(Z, dtype=float).ravel() + + r = 0.5 * R_0 * (_GL_NODES + 1) + w_r = 0.5 * R_0 * _GL_WEIGHTS + eps = 0.5 * np.pi * _GL_NODES + w_eps = 0.5 * np.pi * _GL_WEIGHTS + Rn, En = np.meshgrid(r, eps, indexing='ij') + W2D = np.outer(w_r, w_eps) + sinE = np.sin(En) + + n_pts = X.size + I_K = np.empty(n_pts) + I_M = np.empty(n_pts) + I_W = np.empty(n_pts) + I_N = np.empty(n_pts) + for start in range(0, n_pts, _chunk): + sl = slice(start, min(start + _chunk, n_pts)) + Xc = X[sl][:, None, None] + Zc = Z[sl][:, None, None] + Q = Xc ** 2 / (Xc ** 2 + Zc ** 2 - 2 * Zc * Rn * sinE + Rn ** 2) + K = get_K_factor(Q, S_0) + M = get_M_factor(Q, S_0) + N = get_N_factor(Q, S_0) + I_K[sl] = np.sum(W2D * Rn * K, axis=(1, 2)) + I_M[sl] = np.sum(W2D * Rn * M, axis=(1, 2)) + I_W[sl] = np.sum(W2D * (Zc - Rn * sinE) * Rn * M, axis=(1, 2)) + I_N[sl] = np.sum(W2D * Rn * N, axis=(1, 2)) + + n = I_K / (np.pi ** 1.5 * X ** 2) + U = I_M / I_K + W = I_W / (X * I_K) + T = -(2 / 3) * (U ** 2 + W ** 2) + (4 / 3) * I_N / I_K + values = {'n': n, 'U': U, 'W': W, 'T': T, + 'Vr': (X * U + Z * W) / np.sqrt(X ** 2 + Z ** 2), + 'p': n * T} + _verify_against_class(values, X, Z, S_0) + return values[quantity] + + +def _verify_against_class(values, X, Z, S_0, rtol=1e-6): + """Cross-check sample points against CollisionlessGasKinetics.""" + for idx in {0, X.size // 2, X.size - 1}: + d = float(np.hypot(X[idx], Z[idx])) + th = float(np.arctan2(Z[idx], X[idx])) + ref = make_plume(CollisionlessGasKinetics, d, th, S_0) + checks = {'n': ref.get_num_density_ratio(), + 'U': ref.get_U_normalized(), + 'T': ref.get_temp_ratio()} + for q, expected in checks.items(): + got = values[q][idx] + if abs(got - expected) > rtol * abs(expected): + raise AssertionError( + f'vectorized field {q} diverged from ' + f'CollisionlessGasKinetics at (X={X[idx]}, Z={Z[idx]}): ' + f'{got} vs {expected}') + + +def full_field_grid(quantity, x_over_D, z_over_D, S_0): + """Full-model quantity on the tensor grid; returns shape + (len(z_over_D), len(x_over_D)) for use with plt.contour.""" + Xg, Zg = np.meshgrid(np.asarray(x_over_D) * D_NOZZLE, + np.asarray(z_over_D) * D_NOZZLE) + vals = full_field_values(quantity, Xg.ravel(), Zg.ravel(), S_0) + return vals.reshape(Xg.shape) + + +def simplified_field_grid(quantity, x_over_D, z_over_D, S_0): + """Simplified-model quantity (Eqs. 13-17) on the tensor grid.""" + Xg, Zg = np.meshgrid(np.asarray(x_over_D) * D_NOZZLE, + np.asarray(z_over_D) * D_NOZZLE) + Q = Xg ** 2 / (Xg ** 2 + Zg ** 2) + K = get_K_factor(Q, S_0) + M = get_M_factor(Q, S_0) + N = get_N_factor(Q, S_0) + n = K / (2 * np.sqrt(np.pi)) * (R_0 / Xg) ** 2 + U = M / K + W = U * Zg / Xg + T = -2 * M ** 2 / (3 * Q * K ** 2) + 4 * N / (3 * K) + values = {'n': n, 'U': U, 'W': W, 'T': T, + 'Vr': (Xg * U + Zg * W) / np.sqrt(Xg ** 2 + Zg ** 2), + 'p': n * T} + return values[quantity] + + +def default_contour_axes(x_max=10.0, z_max=10.0): + """Paper-style contour grid: X/D in [X_MIN, x_max], Z/D in [0, z_max] + (fields are symmetric in Z; mirror for the lower half-plane).""" + x = np.linspace(X_MIN_OVER_D, x_max, GRID_N) + z = np.linspace(0.0, z_max, (GRID_N + 1) // 2) + return x, z + + +# --------------------------------------------------------------------------- +# digitized-data overlay (see data/digitized/README.md) +# --------------------------------------------------------------------------- + +def load_digitized(fig_stem): + """dict {series_suffix: (x, y)} from data/digitized/_*.csv.""" + datasets = {} + for path in sorted(DIGITIZED_DIR.glob(f'{fig_stem}_*.csv')): + data = np.genfromtxt(path, delimiter=',', skip_header=1) + data = np.atleast_2d(data) + datasets[path.stem[len(fig_stem) + 1:]] = (data[:, 0], data[:, 1]) + return datasets + + +def _pretty_series_label(suffix): + tokens = [] + for token in suffix.split('_'): + if token.lower() == 'dsmc': + tokens.append('DSMC') + elif token.lower().startswith('kn'): + tokens.append('Kn=' + token[2:].replace('p', '.')) + else: + tokens.append(token) + return ' '.join(tokens) + + +def overlay_digitized(ax, fig_stem, style='markers', mirror_z=False, **kw): + """Overlay every digitized series for fig_stem; returns #series. + + style 'markers' suits profile figures; 'line' suits digitized + contour polylines. mirror_z flips y (for lower-half-plane DSMC + slots in the split contour figures). + """ + count = 0 + for suffix, (x, y) in load_digitized(fig_stem).items(): + y = -y if mirror_z else y + label = _pretty_series_label(suffix) + if style == 'line': + ax.plot(x, y, color='k', lw=1.0, label=label, **kw) + else: + ax.plot(x, y, linestyle='none', marker='d', mfc='none', + color='k', ms=4, label=label, **kw) + count += 1 + return count + + +# --------------------------------------------------------------------------- +# figure assembly helpers +# --------------------------------------------------------------------------- + +def max_rel_diff(candidate, reference): + """max |candidate - reference| / |reference| over finite entries.""" + candidate = np.asarray(candidate, dtype=float) + reference = np.asarray(reference, dtype=float) + mask = np.isfinite(candidate) & np.isfinite(reference) & (reference != 0) + return float(np.max(np.abs(candidate[mask] - reference[mask]) + / np.abs(reference[mask]))) + + +def annotate_error(ax, text, loc='lower left'): + """Small model-difference annotation box on the axes.""" + positions = {'lower left': (0.02, 0.02, 'left', 'bottom'), + 'lower right': (0.98, 0.02, 'right', 'bottom'), + 'upper left': (0.02, 0.98, 'left', 'top')} + x, y, ha, va = positions[loc] + ax.text(x, y, text, transform=ax.transAxes, fontsize=7, ha=ha, va=va, + bbox=dict(boxstyle='round', fc='white', ec='0.6', alpha=0.85)) + + +def save_figure(fig, name): + """Save PNG to tests/plume/output and return the path.""" + OUTPUT_DIR.mkdir(parents=True, exist_ok=True) + path = OUTPUT_DIR / f'{name}.png' + fig.savefig(path, dpi=200, bbox_inches='tight') + plt.close(fig) + return path + + +def run_script(generate_figure): + """Manual-run entry point shared by all figure scripts.""" + t0 = time.perf_counter() + path = generate_figure() + print(f'saved {path} ({time.perf_counter() - t0:.1f} s)') + + +# --------------------------------------------------------------------------- +# shared builders for the Kn-triplet figures (identical analytic content; +# separate figures so per-Kn DSMC overlays land in the right file) +# --------------------------------------------------------------------------- + +def density_contour_figure(kn_label, fig_stem): + """Figs. 6-8: number-density contours, S_0 = 2. Upper half-plane: + analytical (solid) + simplified (dashed); lower half-plane: Simons + (dashed) + DSMC digitized slot (solid lines when data present).""" + S_0 = 2.0 + levels = [0.001, 0.01, 0.1] + x, z = default_contour_axes() + n_full = full_field_grid('n', x, z, S_0) + n_simp = simplified_field_grid('n', x, z, S_0) + Xg, Zg = np.meshgrid(x * D_NOZZLE, z * D_NOZZLE) + n_simons = simons_density_field(Xg, Zg, S_0) + + fig, ax = plt.subplots(figsize=(6, 7)) + cs = ax.contour(x, z, n_full, levels=levels, colors='k', + linewidths=1.2) + ax.clabel(cs, fmt='%g', fontsize=7) + cs = ax.contour(x, z, n_simp, levels=levels, colors='k', + linewidths=0.9, linestyles='dashed') + ax.clabel(cs, fmt='%g', fontsize=7) + cs = ax.contour(x, -z, n_simons, levels=levels, colors='k', + linewidths=0.9, linestyles='dashed') + ax.clabel(cs, fmt='%g', fontsize=7) + n_dig = overlay_digitized(ax, fig_stem, style='line', mirror_z=False) + + ax.axhline(0.0, color='k', lw=0.6) + ax.set_xlim(0, 10) + ax.set_ylim(-10, 10) + ax.set_xlabel('X/D') + ax.set_ylabel('Z/D') + ax.set_title(f'Normalized number density, Kn = {kn_label}, ' + f'$S_0$ = 2.0') + handles = [plt.Line2D([], [], color='k', lw=1.2, label='Analytical (top)'), + plt.Line2D([], [], color='k', lw=0.9, ls='--', + label='Simplified (top)'), + plt.Line2D([], [], color='k', lw=0.9, ls='--', + label='Simons (bottom)')] + if n_dig: + handles.append(plt.Line2D([], [], color='k', lw=1.0, + label='DSMC (bottom, digitized)')) + ax.legend(handles=handles, fontsize=7, loc='upper right') + annotate_error( + ax, + f'max |n_S/n_A - 1| = {max_rel_diff(n_simp, n_full):.2g}\n' + f'max |n_Sim/n_A - 1| = {max_rel_diff(n_simons, n_full):.2g}\n' + '(over plotted grid; largest in the lip/tail regions)') + return save_figure(fig, fig_stem) + + +def top_bottom_contour_figure(quantity, levels, kn_label, fig_stem, + title_quantity, x_max=10.0, z_max=10.0, + z_axis_label='Z/D', fmt='%g'): + """Figs. 11-18 pattern: analytical contours in the upper half-plane, + lower half-plane reserved for the DSMC digitized overlay.""" + S_0 = 2.0 + x, z = default_contour_axes(x_max, z_max) + field = full_field_grid(quantity, x, z, S_0) + + fig, ax = plt.subplots(figsize=(6, 7)) + cs = ax.contour(x, z, field, levels=levels, colors='k', linewidths=1.1) + ax.clabel(cs, fmt=fmt, fontsize=7) + n_dig = overlay_digitized(ax, fig_stem, style='line') + + ax.axhline(0.0, color='k', lw=0.6) + ax.set_xlim(0, x_max) + ax.set_ylim(-z_max, z_max) + ax.set_xlabel('X/D') + ax.set_ylabel(z_axis_label) + ax.set_title(f'{title_quantity}, Kn = {kn_label}, $S_0$ = 2.0') + ax.text(0.7 * x_max, 0.55 * z_max, 'Analytical', fontsize=9) + ax.text(0.7 * x_max, -0.55 * z_max, + 'DSMC' if n_dig else 'DSMC (pending digitized data)', + fontsize=9) + return save_figure(fig, fig_stem) + + +def centerline_density_profile_figure(kn_label, fig_stem): + """Figs. 19-21: centerline density -- analytical (Eq. 18, circles), + simplified (Eq. 14, triangles), Simons (solid line, exit-referenced), + DSMC digitized slot.""" + S_0 = 2.0 + x_over_D = np.linspace(X_MIN_OVER_D, 10.0, 200) + analytical = centerline_profiles(x_over_D, S_0, ('n',))['n'] + simplified = simplified_centerline_density(x_over_D, S_0) + r = x_over_D * D_NOZZLE + simons = simons_density_field(r, np.zeros_like(r), S_0) + + fig, ax = plt.subplots(figsize=(6, 4.5)) + ax.plot(x_over_D, simons, 'k-', lw=1.2, label='Simons') + ax.plot(x_over_D, analytical, linestyle='none', marker='o', mfc='none', + color='k', ms=4, markevery=4, label='Analytical') + ax.plot(x_over_D, simplified, linestyle='none', marker='>', color='k', + ms=4, markevery=(2, 4), label='Simplified') + overlay_digitized(ax, fig_stem) + + ax.set_xlim(0, 10) + ax.set_ylim(0, 1.2) + ax.set_xlabel('X/D') + ax.set_ylabel('$n/n_0$') + ax.set_title(f'Centerline density profiles, Kn = {kn_label}, ' + f'$S_0$ = 2.0') + ax.legend(fontsize=8) + tail = x_over_D >= 1.0 + annotate_error( + ax, + f'max |n_S/n_A - 1| = {max_rel_diff(simplified[tail], analytical[tail]):.2g} ' + f'(X/D $\\geq$ 1)\n' + f'max |n_Sim/n_A - 1| = {max_rel_diff(simons[tail], analytical[tail]):.2g} ' + f'(X/D $\\geq$ 1)', loc='lower right') + return save_figure(fig, fig_stem) + + +def angular_density_profile_figure(kn_label, fig_stem): + """Figs. 22-24: density along r/D = 10 -- analytical, simplified, + Simons kappa = 1.5/2/3 (single Boyton-A normalization, see module + docstring), DSMC digitized slot.""" + S_0 = 2.0 + r = 10.0 * D_NOZZLE + theta = np.deg2rad(np.linspace(0.0, THETA_MAX_DEG, 90)) + X = r * np.cos(theta) + Z = r * np.sin(theta) + analytical = full_field_values('n', X, Z, S_0) + Q = X ** 2 / (X ** 2 + Z ** 2) + simplified = (get_K_factor(Q, S_0) / (2 * np.sqrt(np.pi)) + * (R_0 / X) ** 2) + theta_deg = np.rad2deg(theta) + + fig, ax = plt.subplots(figsize=(6, 4.5)) + for kappa, ls in [(2.0, '-'), (1.5, '--'), (3.0, '-.')]: + simons = simons_density_field(X, Z, S_0, kappa_f=kappa) + ax.plot(theta_deg, simons, 'k', ls=ls, lw=1.1, + label=f'Simons $\\kappa$={kappa:g}') + ax.plot(theta_deg, analytical, linestyle='none', marker='o', mfc='none', + color='k', ms=4, markevery=2, label='Analytical') + ax.plot(theta_deg, simplified, linestyle='none', marker='>', color='k', + ms=4, markevery=(1, 2), label='Simplified') + overlay_digitized(ax, fig_stem) + + ax.set_xlim(0, 90) + ax.set_ylim(0, 0.015) + ax.set_xticks([0, 30, 60, 90]) + ax.set_xlabel(r'$\theta$') + ax.set_ylabel('$n/n_0$') + ax.set_title(f'Density profiles along r/D = 10, Kn = {kn_label}, ' + f'$S_0$ = 2.0') + ax.legend(fontsize=8) + simons_boyton = simons_density_field(X, Z, S_0) + annotate_error( + ax, + f'max |n_S/n_A - 1| = {max_rel_diff(simplified, analytical):.2g}\n' + f'max |n_Sim($\\kappa$=3)/n_A - 1| = ' + f'{max_rel_diff(simons_boyton, analytical):.2g}', + loc='lower left') + return save_figure(fig, fig_stem) From 32a2437ba2e462b144fae46347159ab5b86b5916 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Sat, 18 Jul 2026 07:51:04 -0500 Subject: [PATCH 10/14] add plume boundary and centerline profile figures Co-Authored-By: Claude Fable 5 --- tests/plume/plume_verification_test_04.py | 44 +++++++++++++++++++++++ tests/plume/plume_verification_test_05.py | 38 ++++++++++++++++++++ tests/plume/plume_verification_test_06.py | 37 +++++++++++++++++++ tests/plume/plume_verification_test_07.py | 38 ++++++++++++++++++++ 4 files changed, 157 insertions(+) create mode 100644 tests/plume/plume_verification_test_04.py create mode 100644 tests/plume/plume_verification_test_05.py create mode 100644 tests/plume/plume_verification_test_06.py create mode 100644 tests/plume/plume_verification_test_07.py diff --git a/tests/plume/plume_verification_test_04.py b/tests/plume/plume_verification_test_04.py new file mode 100644 index 0000000..5978a31 --- /dev/null +++ b/tests/plume/plume_verification_test_04.py @@ -0,0 +1,44 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_04.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 2: plume boundaries (n/n0 = 0.001 +# contour of the full analytical model) for exit speed ratios +# S0 = 1, 2, 3; X/D in [0, 30], Z/D in [0, 15]. +# +# Manual-run verification script (design decision D5): it defines no +# pytest tests and is executed directly -- +# python tests/plume/plume_verification_test_04.py +# The figure is saved to tests/plume/output/. +# Digitized overlays looked for: fig02_*.csv (none published -- Fig. 2 +# has no DSMC data in the paper). + +import numpy as np + +import plume_figure_utils as u + + +def generate_figure(): + x = np.linspace(u.X_MIN_OVER_D, 30.0, 181) + z = np.linspace(0.0, 15.0, 91) + + fig, ax = u.plt.subplots(figsize=(7, 4.5)) + for S_0, ls in [(1.0, '-'), (2.0, '--'), (3.0, '-.')]: + n = u.full_field_grid('n', x, z, S_0) + cs = ax.contour(x, z, n, levels=[0.001], colors='k', + linestyles=ls, linewidths=1.2) + ax.clabel(cs, fmt='0.001', fontsize=7) + ax.plot([], [], 'k', ls=ls, label=f'$S_0$={S_0:g}') + u.overlay_digitized(ax, 'fig02') + + ax.set_xlim(0, 30) + ax.set_ylim(0, 15) + ax.set_xlabel('X/D') + ax.set_ylabel('Z/D') + ax.set_title('Plume boundaries at different exit speed ratios ' + '(Fig. 2)') + ax.legend(fontsize=8) + return u.save_figure(fig, 'fig02_plume_boundaries') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_05.py b/tests/plume/plume_verification_test_05.py new file mode 100644 index 0000000..430f282 --- /dev/null +++ b/tests/plume/plume_verification_test_05.py @@ -0,0 +1,38 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_05.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 3: normalized analytical number +# density along the centerline (Eq. 18) vs X/D for S0 = 1, 2, 3. +# +# Manual-run verification script (design decision D5): it defines no +# pytest tests and is executed directly -- +# python tests/plume/plume_verification_test_05.py +# Digitized overlay looked for: fig03_dsmc.csv (the paper's S0 = 2 +# DSMC circles). + +import numpy as np + +import plume_figure_utils as u + + +def generate_figure(): + x_over_D = np.linspace(u.X_MIN_OVER_D, 10.0, 200) + + fig, ax = u.plt.subplots(figsize=(6, 4.5)) + for S_0, ls in [(1.0, '-.'), (2.0, '--'), (3.0, '-')]: + n = u.centerline_profiles(x_over_D, S_0, ('n',))['n'] + ax.plot(x_over_D, n, 'k', ls=ls, lw=1.2, label=f'$S_0$={S_0:g}') + u.overlay_digitized(ax, 'fig03') + + ax.set_xlim(0, 10) + ax.set_ylim(0, 1.2) + ax.set_xlabel('X/D') + ax.set_ylabel('$n_1/n_0$') + ax.set_title('Normalized analytical number density along centerline ' + '(Fig. 3)') + ax.legend(fontsize=8) + return u.save_figure(fig, 'fig03_centerline_density') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_06.py b/tests/plume/plume_verification_test_06.py new file mode 100644 index 0000000..829d5dd --- /dev/null +++ b/tests/plume/plume_verification_test_06.py @@ -0,0 +1,37 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_06.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 4: normalized analytical U-velocity +# along the centerline (Eq. 19), U1*sqrt(beta0) vs X/D, S0 = 1, 2, 3. +# +# Manual-run verification script (design decision D5): it defines no +# pytest tests and is executed directly -- +# python tests/plume/plume_verification_test_06.py +# Digitized overlay looked for: fig04_dsmc.csv (S0 = 2 DSMC circles). + +import numpy as np + +import plume_figure_utils as u + + +def generate_figure(): + x_over_D = np.linspace(u.X_MIN_OVER_D, 10.0, 200) + + fig, ax = u.plt.subplots(figsize=(6, 4.5)) + for S_0, ls in [(1.0, '-.'), (2.0, '--'), (3.0, '-')]: + U = u.centerline_profiles(x_over_D, S_0, ('U',))['U'] + ax.plot(x_over_D, U, 'k', ls=ls, lw=1.2, label=f'$S_0$={S_0:g}') + u.overlay_digitized(ax, 'fig04') + + ax.set_xlim(0, 10) + ax.set_ylim(1, 4) + ax.set_xlabel('X/D') + ax.set_ylabel(r'$U_1\sqrt{\beta_0}$') + ax.set_title('Normalized analytical U-velocity along centerline ' + '(Fig. 4)') + ax.legend(fontsize=8, loc='lower right') + return u.save_figure(fig, 'fig04_centerline_velocity') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_07.py b/tests/plume/plume_verification_test_07.py new file mode 100644 index 0000000..3ca4a96 --- /dev/null +++ b/tests/plume/plume_verification_test_07.py @@ -0,0 +1,38 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_07.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 5: normalized analytical temperature +# along the centerline (Eq. 21, exact quadrature), T1/T0 vs X/D for +# S0 = 1, 2, 3. +# +# Manual-run verification script (design decision D5): it defines no +# pytest tests and is executed directly -- +# python tests/plume/plume_verification_test_07.py +# Digitized overlay looked for: fig05_dsmc.csv (S0 = 2 DSMC circles). + +import numpy as np + +import plume_figure_utils as u + + +def generate_figure(): + x_over_D = np.linspace(u.X_MIN_OVER_D, 10.0, 200) + + fig, ax = u.plt.subplots(figsize=(6, 4.5)) + for S_0, ls in [(1.0, '-.'), (2.0, '--'), (3.0, '-')]: + T = u.centerline_profiles(x_over_D, S_0, ('T',))['T'] + ax.plot(x_over_D, T, 'k', ls=ls, lw=1.2, label=f'$S_0$={S_0:g}') + u.overlay_digitized(ax, 'fig05') + + ax.set_xlim(0, 10) + ax.set_ylim(0.2, 1.2) + ax.set_xlabel('X/D') + ax.set_ylabel('$T_1/T_0$') + ax.set_title('Normalized analytical temperature along centerline ' + '(Fig. 5)') + ax.legend(fontsize=8) + return u.save_figure(fig, 'fig05_centerline_temperature') + + +if __name__ == '__main__': + u.run_script(generate_figure) From 27423da818ab54a649e0a03316c075e44c4b4167 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Sat, 18 Jul 2026 07:53:55 -0500 Subject: [PATCH 11/14] add flowfield contour verification figures Co-Authored-By: Claude Fable 5 --- tests/plume/plume_figure_utils.py | 7 ++-- tests/plume/plume_verification_test_08.py | 23 +++++++++++ tests/plume/plume_verification_test_09.py | 22 ++++++++++ tests/plume/plume_verification_test_10.py | 22 ++++++++++ tests/plume/plume_verification_test_11.py | 50 +++++++++++++++++++++++ tests/plume/plume_verification_test_12.py | 47 +++++++++++++++++++++ tests/plume/plume_verification_test_13.py | 25 ++++++++++++ tests/plume/plume_verification_test_14.py | 22 ++++++++++ tests/plume/plume_verification_test_15.py | 25 ++++++++++++ tests/plume/plume_verification_test_16.py | 22 ++++++++++ tests/plume/plume_verification_test_17.py | 25 ++++++++++++ tests/plume/plume_verification_test_18.py | 26 ++++++++++++ tests/plume/plume_verification_test_19.py | 24 +++++++++++ tests/plume/plume_verification_test_20.py | 24 +++++++++++ 14 files changed, 361 insertions(+), 3 deletions(-) create mode 100644 tests/plume/plume_verification_test_08.py create mode 100644 tests/plume/plume_verification_test_09.py create mode 100644 tests/plume/plume_verification_test_10.py create mode 100644 tests/plume/plume_verification_test_11.py create mode 100644 tests/plume/plume_verification_test_12.py create mode 100644 tests/plume/plume_verification_test_13.py create mode 100644 tests/plume/plume_verification_test_14.py create mode 100644 tests/plume/plume_verification_test_15.py create mode 100644 tests/plume/plume_verification_test_16.py create mode 100644 tests/plume/plume_verification_test_17.py create mode 100644 tests/plume/plume_verification_test_18.py create mode 100644 tests/plume/plume_verification_test_19.py create mode 100644 tests/plume/plume_verification_test_20.py diff --git a/tests/plume/plume_figure_utils.py b/tests/plume/plume_figure_utils.py index eae7606..76120c7 100644 --- a/tests/plume/plume_figure_utils.py +++ b/tests/plume/plume_figure_utils.py @@ -404,11 +404,12 @@ def density_contour_figure(kn_label, fig_stem): handles.append(plt.Line2D([], [], color='k', lw=1.0, label='DSMC (bottom, digitized)')) ax.legend(handles=handles, fontsize=7, loc='upper right') + far = x >= 1.0 # near the exit the simplified (R_0/X)^2 form diverges annotate_error( ax, - f'max |n_S/n_A - 1| = {max_rel_diff(n_simp, n_full):.2g}\n' - f'max |n_Sim/n_A - 1| = {max_rel_diff(n_simons, n_full):.2g}\n' - '(over plotted grid; largest in the lip/tail regions)') + f'max |n_S/n_A - 1| = {max_rel_diff(n_simp[:, far], n_full[:, far]):.2g}\n' + f'max |n_Sim/n_A - 1| = {max_rel_diff(n_simons[:, far], n_full[:, far]):.2g}\n' + '(over grid with X/D $\\geq$ 1; largest in the plume tails)') return save_figure(fig, fig_stem) diff --git a/tests/plume/plume_verification_test_08.py b/tests/plume/plume_verification_test_08.py new file mode 100644 index 0000000..bf5baec --- /dev/null +++ b/tests/plume/plume_verification_test_08.py @@ -0,0 +1,23 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_08.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 6: normalized number-density contours, +# Kn = 100, S0 = 2.0. Upper half-plane: analytical (solid) + simplified +# (dashed); lower half-plane: Simons (dashed) + DSMC slot. The analytic +# curves are Kn-independent (collisionless); Figs. 6-8 differ only in +# their DSMC data, so each gets its own script/overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_08.py +# Digitized overlays looked for: fig06_*.csv (e.g. fig06_dsmc_0p001.csv, +# one CSV per digitized DSMC contour polyline, lower half-plane). + +import plume_figure_utils as u + + +def generate_figure(): + return u.density_contour_figure('100', 'fig06') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_09.py b/tests/plume/plume_verification_test_09.py new file mode 100644 index 0000000..5dae292 --- /dev/null +++ b/tests/plume/plume_verification_test_09.py @@ -0,0 +1,22 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_09.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 7: normalized number-density contours, +# Kn = 0.1, S0 = 2.0. Same analytic content as Fig. 6 (collisionless +# models are Kn-independent); the Kn distinction lives in the DSMC +# overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_09.py +# Digitized overlays looked for: fig07_*.csv (one CSV per digitized +# DSMC contour polyline, lower half-plane). + +import plume_figure_utils as u + + +def generate_figure(): + return u.density_contour_figure('0.1', 'fig07') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_10.py b/tests/plume/plume_verification_test_10.py new file mode 100644 index 0000000..f5ab6b3 --- /dev/null +++ b/tests/plume/plume_verification_test_10.py @@ -0,0 +1,22 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_10.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 8: normalized number-density contours, +# Kn = 0.01, S0 = 2.0. Same analytic content as Fig. 6 (collisionless +# models are Kn-independent); the Kn distinction lives in the DSMC +# overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_10.py +# Digitized overlays looked for: fig08_*.csv (one CSV per digitized +# DSMC contour polyline, lower half-plane). + +import plume_figure_utils as u + + +def generate_figure(): + return u.density_contour_figure('0.01', 'fig08') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_11.py b/tests/plume/plume_verification_test_11.py new file mode 100644 index 0000000..deff98a --- /dev/null +++ b/tests/plume/plume_verification_test_11.py @@ -0,0 +1,50 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_11.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 9: relative density error between the +# analytical and Simons (cosine-law) solutions, |n_A/n_Simons - 1| +# [Eq. 28], Kn = 100 (analytic content is Kn-independent), S0 = 2.0, +# X/D and Y/D in [0, 10]. The Simons field is exit-referenced with the +# Boyton kappa (see plume_figure_utils docstring); the paper's A +# normalization additionally depends on exit Mach number, so absolute +# error magnitudes may differ somewhat from the printed figure. +# +# The paper notes the upper-nozzle-lip region is a density singularity +# whose error pattern "shall be neglected". +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_11.py +# Digitized overlays looked for: fig09_*.csv. + +import numpy as np + +import plume_figure_utils as u + + +def generate_figure(): + S_0 = 2.0 + x, y = u.default_contour_axes() + n_full = u.full_field_grid('n', x, y, S_0) + Xg, Yg = np.meshgrid(x * u.D_NOZZLE, y * u.D_NOZZLE) + n_simons = u.simons_density_field(Xg, Yg, S_0) + error = np.abs(n_full / n_simons - 1) + + fig, ax = u.plt.subplots(figsize=(6, 5)) + levels = [0.1, 0.4, 1, 2, 4, 5] + cs = ax.contour(x, y, error, levels=levels, colors='k', linewidths=1.0) + ax.clabel(cs, fmt='%g', fontsize=7) + u.overlay_digitized(ax, 'fig09', style='line') + + ax.set_xlim(0, 10) + ax.set_ylim(0, 10) + ax.set_xlabel('X/D') + ax.set_ylabel('Y/D') + ax.set_title('Relative density error: analytical vs Simons (Fig. 9)') + u.annotate_error( + ax, f'max |n_A/n_Sim - 1| = {np.nanmax(error):.2g} ' + '(peak at the nozzle-lip singularity)', loc='lower right') + return u.save_figure(fig, 'fig09_error_analytical_vs_simons') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_12.py b/tests/plume/plume_verification_test_12.py new file mode 100644 index 0000000..1e2f62f --- /dev/null +++ b/tests/plume/plume_verification_test_12.py @@ -0,0 +1,47 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_12.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 10: relative density error between the +# analytical and simplified analytical solutions, |n_A/n_As - 1| +# [Eq. 29], Kn = 100 (analytic content is Kn-independent), S0 = 2.0, +# X/D and Y/D in [0, 10]. +# +# The paper notes the upper-nozzle-lip region is a density singularity +# whose error pattern "shall be neglected". +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_12.py +# Digitized overlays looked for: fig10_*.csv. + +import numpy as np + +import plume_figure_utils as u + + +def generate_figure(): + S_0 = 2.0 + x, y = u.default_contour_axes() + n_full = u.full_field_grid('n', x, y, S_0) + n_simp = u.simplified_field_grid('n', x, y, S_0) + error = np.abs(n_full / n_simp - 1) + + fig, ax = u.plt.subplots(figsize=(6, 5)) + levels = [0.005, 0.03, 0.05, 0.1, 0.2, 0.5, 1, 2] + cs = ax.contour(x, y, error, levels=levels, colors='k', linewidths=1.0) + ax.clabel(cs, fmt='%g', fontsize=7) + u.overlay_digitized(ax, 'fig10', style='line') + + ax.set_xlim(0, 10) + ax.set_ylim(0, 10) + ax.set_xlabel('X/D') + ax.set_ylabel('Y/D') + ax.set_title('Relative density error: analytical vs simplified ' + '(Fig. 10)') + u.annotate_error( + ax, f'max |n_A/n_As - 1| = {np.nanmax(error):.2g} ' + '(peak at the nozzle-lip singularity)', loc='lower right') + return u.save_figure(fig, 'fig10_error_analytical_vs_simplified') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_13.py b/tests/plume/plume_verification_test_13.py new file mode 100644 index 0000000..09c2385 --- /dev/null +++ b/tests/plume/plume_verification_test_13.py @@ -0,0 +1,25 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_13.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 11: normalized pressure contours +# p1/p0 = (n1/n0)(T1/T0), Kn = 100, S0 = 2.0. Analytical field in the +# upper half-plane; lower half-plane reserved for the DSMC digitized +# overlay. Analytic content is Kn-independent; Figs. 11-12 differ only +# in their DSMC data. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_13.py +# Digitized overlays looked for: fig11_*.csv (one CSV per digitized +# DSMC contour polyline, lower half-plane). + +import plume_figure_utils as u + + +def generate_figure(): + return u.top_bottom_contour_figure( + 'p', [1e-5, 1e-4, 1e-3, 1e-2], '100', 'fig11', + 'Normalized pressure contours') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_14.py b/tests/plume/plume_verification_test_14.py new file mode 100644 index 0000000..a2a2bfe --- /dev/null +++ b/tests/plume/plume_verification_test_14.py @@ -0,0 +1,22 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_14.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 12: normalized pressure contours +# p1/p0, Kn = 0.1, S0 = 2.0. Same analytic content as Fig. 11; the Kn +# distinction lives in the DSMC overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_14.py +# Digitized overlays looked for: fig12_*.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.top_bottom_contour_figure( + 'p', [1e-5, 1e-4, 1e-3, 1e-2], '0.1', 'fig12', + 'Normalized pressure contours') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_15.py b/tests/plume/plume_verification_test_15.py new file mode 100644 index 0000000..493cf1b --- /dev/null +++ b/tests/plume/plume_verification_test_15.py @@ -0,0 +1,25 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_15.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 13: normalized temperature contours +# T1/T0, Kn = 100, S0 = 2.0; X/D in [0, 2.5], Y/D in [-2.5, 2.5]. +# Analytical field in the upper half-plane; lower half-plane reserved +# for the DSMC digitized overlay. This is the figure the paper uses to +# argue p = n*k*T0 is invalid (T1 < T0 everywhere downstream). +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_15.py +# Digitized overlays looked for: fig13_*.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.top_bottom_contour_figure( + 'T', [0.3, 0.4, 0.6, 0.9], '100', 'fig13', + 'Normalized temperature distribution', + x_max=2.5, z_max=2.5, z_axis_label='Y/D') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_16.py b/tests/plume/plume_verification_test_16.py new file mode 100644 index 0000000..f65f68f --- /dev/null +++ b/tests/plume/plume_verification_test_16.py @@ -0,0 +1,22 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_16.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 14: normalized U-velocity contours +# U1*sqrt(beta0), Kn = 100, S0 = 2.0. Analytical field in the upper +# half-plane; lower half-plane reserved for the DSMC digitized overlay. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_16.py +# Digitized overlays looked for: fig14_*.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.top_bottom_contour_figure( + 'U', [0.6, 1.0, 1.8, 2.4], '100', 'fig14', + 'Normalized U-velocity distribution') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_17.py b/tests/plume/plume_verification_test_17.py new file mode 100644 index 0000000..2c9a9e6 --- /dev/null +++ b/tests/plume/plume_verification_test_17.py @@ -0,0 +1,25 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_17.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 15: normalized "V-velocity" contours, +# Kn = 100, S0 = 2.0. In the plotted XOZ plane (Y = 0) the y-component +# V is identically zero by axisymmetry; the transverse component shown +# in the paper's Fig. 15 corresponds to the model's W (Eq. 7), which is +# what is plotted here. Analytical field in the upper half-plane; lower +# half-plane reserved for the DSMC digitized overlay. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_17.py +# Digitized overlays looked for: fig15_*.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.top_bottom_contour_figure( + 'W', [0.1, 0.5, 1.0], '100', 'fig15', + 'Normalized V-velocity distribution (transverse W at Y = 0)') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_18.py b/tests/plume/plume_verification_test_18.py new file mode 100644 index 0000000..d42e5ee --- /dev/null +++ b/tests/plume/plume_verification_test_18.py @@ -0,0 +1,26 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_18.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 16: normalized radial velocity Vr +# contours, Kn = 100, S0 = 2.0. Analytical field in the upper +# half-plane; lower half-plane reserved for the DSMC digitized overlay. +# Analytic content is Kn-independent; Figs. 16-18 differ only in their +# DSMC data. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_18.py +# Digitized overlays looked for: fig16_*.csv. + +import plume_figure_utils as u + +VR_LEVELS = [1.0, 1.8, 2.0, 2.2, 2.4, 2.8] + + +def generate_figure(): + return u.top_bottom_contour_figure( + 'Vr', VR_LEVELS, '100', 'fig16', + 'Normalized velocity $V_r$ distribution') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_19.py b/tests/plume/plume_verification_test_19.py new file mode 100644 index 0000000..150635d --- /dev/null +++ b/tests/plume/plume_verification_test_19.py @@ -0,0 +1,24 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_19.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 17: normalized radial velocity Vr +# contours, Kn = 0.1, S0 = 2.0. Same analytic content as Fig. 16; the +# Kn distinction lives in the DSMC overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_19.py +# Digitized overlays looked for: fig17_*.csv. + +import plume_figure_utils as u + +from plume_verification_test_18 import VR_LEVELS + + +def generate_figure(): + return u.top_bottom_contour_figure( + 'Vr', VR_LEVELS, '0.1', 'fig17', + 'Normalized velocity $V_r$ distribution') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_20.py b/tests/plume/plume_verification_test_20.py new file mode 100644 index 0000000..5e8b6e7 --- /dev/null +++ b/tests/plume/plume_verification_test_20.py @@ -0,0 +1,24 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_20.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 18: normalized radial velocity Vr +# contours, Kn = 0.01, S0 = 2.0. Same analytic content as Fig. 16; the +# Kn distinction lives in the DSMC overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_20.py +# Digitized overlays looked for: fig18_*.csv. + +import plume_figure_utils as u + +from plume_verification_test_18 import VR_LEVELS + + +def generate_figure(): + return u.top_bottom_contour_figure( + 'Vr', VR_LEVELS, '0.01', 'fig18', + 'Normalized velocity $V_r$ distribution') + + +if __name__ == '__main__': + u.run_script(generate_figure) From 85f03f9bc3b1987a0d78837276d2dadad6dab529 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Sat, 18 Jul 2026 07:55:13 -0500 Subject: [PATCH 12/14] add centerline, angular and mass flux profile figures Co-Authored-By: Claude Fable 5 --- tests/plume/plume_verification_test_21.py | 22 +++++++++ tests/plume/plume_verification_test_22.py | 20 ++++++++ tests/plume/plume_verification_test_23.py | 20 ++++++++ tests/plume/plume_verification_test_24.py | 24 +++++++++ tests/plume/plume_verification_test_25.py | 20 ++++++++ tests/plume/plume_verification_test_26.py | 20 ++++++++ tests/plume/plume_verification_test_27.py | 59 +++++++++++++++++++++++ 7 files changed, 185 insertions(+) create mode 100644 tests/plume/plume_verification_test_21.py create mode 100644 tests/plume/plume_verification_test_22.py create mode 100644 tests/plume/plume_verification_test_23.py create mode 100644 tests/plume/plume_verification_test_24.py create mode 100644 tests/plume/plume_verification_test_25.py create mode 100644 tests/plume/plume_verification_test_26.py create mode 100644 tests/plume/plume_verification_test_27.py diff --git a/tests/plume/plume_verification_test_21.py b/tests/plume/plume_verification_test_21.py new file mode 100644 index 0000000..ec85e17 --- /dev/null +++ b/tests/plume/plume_verification_test_21.py @@ -0,0 +1,22 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_21.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 19: centerline density profiles, +# Kn = 100, S0 = 2.0 -- analytical (Eq. 18), simplified (Eq. 14), +# Simons (exit-referenced cosine law) and the DSMC digitized slot. +# Analytic content is Kn-independent; Figs. 19-21 differ only in their +# DSMC data. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_21.py +# Digitized overlay looked for: fig19_dsmc.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.centerline_density_profile_figure('100', 'fig19') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_22.py b/tests/plume/plume_verification_test_22.py new file mode 100644 index 0000000..3d428c7 --- /dev/null +++ b/tests/plume/plume_verification_test_22.py @@ -0,0 +1,20 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_22.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 20: centerline density profiles, +# Kn = 0.1, S0 = 2.0. Same analytic content as Fig. 19; the Kn +# distinction lives in the DSMC overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_22.py +# Digitized overlay looked for: fig20_dsmc.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.centerline_density_profile_figure('0.1', 'fig20') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_23.py b/tests/plume/plume_verification_test_23.py new file mode 100644 index 0000000..2f9ff55 --- /dev/null +++ b/tests/plume/plume_verification_test_23.py @@ -0,0 +1,20 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_23.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 21: centerline density profiles, +# Kn = 0.01, S0 = 2.0. Same analytic content as Fig. 19; the Kn +# distinction lives in the DSMC overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_23.py +# Digitized overlay looked for: fig21_dsmc.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.centerline_density_profile_figure('0.01', 'fig21') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_24.py b/tests/plume/plume_verification_test_24.py new file mode 100644 index 0000000..cc3858a --- /dev/null +++ b/tests/plume/plume_verification_test_24.py @@ -0,0 +1,24 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_24.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 22: density profiles along r/D = 10, +# Kn = 100, S0 = 2.0 -- analytical, simplified, Simons kappa = 1.5/2/3 +# and the DSMC digitized slot. Following the paper (p. 65), the Simons +# curves share one normalization constant A (Boyton kappa) so they +# coincide at theta = 0; only the plotted decay exponent varies. +# Analytic content is Kn-independent; Figs. 22-24 differ only in their +# DSMC data. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_24.py +# Digitized overlay looked for: fig22_dsmc.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.angular_density_profile_figure('100', 'fig22') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_25.py b/tests/plume/plume_verification_test_25.py new file mode 100644 index 0000000..3932df9 --- /dev/null +++ b/tests/plume/plume_verification_test_25.py @@ -0,0 +1,20 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_25.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 23: density profiles along r/D = 10, +# Kn = 0.1, S0 = 2.0. Same analytic content as Fig. 22; the Kn +# distinction lives in the DSMC overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_25.py +# Digitized overlay looked for: fig23_dsmc.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.angular_density_profile_figure('0.1', 'fig23') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_26.py b/tests/plume/plume_verification_test_26.py new file mode 100644 index 0000000..93a2a1d --- /dev/null +++ b/tests/plume/plume_verification_test_26.py @@ -0,0 +1,20 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_26.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 24: density profiles along r/D = 10, +# Kn = 0.01, S0 = 2.0. Same analytic content as Fig. 22; the Kn +# distinction lives in the DSMC overlay slot. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_26.py +# Digitized overlay looked for: fig24_dsmc.csv. + +import plume_figure_utils as u + + +def generate_figure(): + return u.angular_density_profile_figure('0.01', 'fig24') + + +if __name__ == '__main__': + u.run_script(generate_figure) diff --git a/tests/plume/plume_verification_test_27.py b/tests/plume/plume_verification_test_27.py new file mode 100644 index 0000000..cf7ee1a --- /dev/null +++ b/tests/plume/plume_verification_test_27.py @@ -0,0 +1,59 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_test_27.py +# ======================== +# Reproduces Cai & Wang 2012 Fig. 25: normalized mass flux along +# r/D = 10 vs theta, S0 = 2.0 -- the analytical curve plus three DSMC +# digitized slots (Kn = 100, 0.1, 0.01). +# +# Normalization note: the paper states the flux is "normalized by the +# inflow mass flux at the exit plane" (p. 66), which would suggest +# rho*Vr / (rho0*U0) = (n/n0)*(Vr*sqrt(beta0))/S0. That convention +# peaks at ~0.0137 at theta = 0, half the paper's ~0.027. The plotted +# magnitudes of Fig. 25 are instead reproduced exactly by +# (n/n0)*(Vr*sqrt(beta0)), i.e. rho*Vr normalized by +# rho0*sqrt(2*R*T0) = rho0*U0/S0. That convention is adopted here +# (verified: theta = 0 value 0.0273 vs the paper's ~0.027); it should +# be re-calibrated against the digitized DSMC data when available. +# +# Manual-run verification script (design decision D5): no pytest tests; +# run directly -- python tests/plume/plume_verification_test_27.py +# Digitized overlays looked for: fig25_dsmc_kn100.csv, +# fig25_dsmc_kn0p1.csv, fig25_dsmc_kn0p01.csv. + +import numpy as np + +import plume_figure_utils as u + + +def mass_flux_curve(theta_deg, S_0=2.0, r_over_D=10.0): + """(n/n0)*(Vr*sqrt(beta0)) along the angular curve (see header).""" + theta = np.deg2rad(theta_deg) + r = r_over_D * u.D_NOZZLE + X = r * np.cos(theta) + Z = r * np.sin(theta) + n = u.full_field_values('n', X, Z, S_0) + Vr = u.full_field_values('Vr', X, Z, S_0) + return n * Vr + + +def generate_figure(): + theta_deg = np.linspace(0.0, u.THETA_MAX_DEG, 90) + flux = mass_flux_curve(theta_deg) + + fig, ax = u.plt.subplots(figsize=(6, 4.5)) + ax.plot(theta_deg, flux, linestyle='none', marker='o', mfc='none', + color='k', ms=4, markevery=2, label='Analytical') + u.overlay_digitized(ax, 'fig25', style='line') + + ax.set_xlim(0, 90) + ax.set_ylim(0, 0.03) + ax.set_xticks([0, 30, 60, 90]) + ax.set_xlabel(r'$\theta$') + ax.set_ylabel('Mass flux') + ax.set_title('Normalized mass flux along r/D = 10 (Fig. 25)') + ax.legend(fontsize=8) + return u.save_figure(fig, 'fig25_mass_flux') + + +if __name__ == '__main__': + u.run_script(generate_figure) From 03cd6681b7fa21756d17a68ed882f69abd12d109 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Sat, 18 Jul 2026 07:57:59 -0500 Subject: [PATCH 13/14] add model error summary and align centerline series with paper Co-Authored-By: Claude Fable 5 --- tests/plume/plume_figure_utils.py | 18 +-- .../plume/plume_verification_error_summary.py | 117 ++++++++++++++++++ 2 files changed, 128 insertions(+), 7 deletions(-) create mode 100644 tests/plume/plume_verification_error_summary.py diff --git a/tests/plume/plume_figure_utils.py b/tests/plume/plume_figure_utils.py index 76120c7..a7004e4 100644 --- a/tests/plume/plume_figure_utils.py +++ b/tests/plume/plume_figure_utils.py @@ -441,14 +441,18 @@ def top_bottom_contour_figure(quantity, levels, kn_label, fig_stem, def centerline_density_profile_figure(kn_label, fig_stem): - """Figs. 19-21: centerline density -- analytical (Eq. 18, circles), - simplified (Eq. 14, triangles), Simons (solid line, exit-referenced), - DSMC digitized slot.""" + """Figs. 19-21: centerline density -- analytical (Eq. 5 integral, + circles), simplified (the Eq. 18 closed form from the paper's + simplified-solutions section II.B, triangles), Simons (solid line, + exit-referenced), DSMC digitized slot. As in the paper, the + analytical and simplified series coincide on the centerline (the + Eq. 5 integral reduces exactly to Eq. 18) and converge to 1 at + X = 0.""" S_0 = 2.0 x_over_D = np.linspace(X_MIN_OVER_D, 10.0, 200) - analytical = centerline_profiles(x_over_D, S_0, ('n',))['n'] - simplified = simplified_centerline_density(x_over_D, S_0) r = x_over_D * D_NOZZLE + analytical = full_field_values('n', r, np.zeros_like(r), S_0) + simplified = centerline_profiles(x_over_D, S_0, ('n',))['n'] simons = simons_density_field(r, np.zeros_like(r), S_0) fig, ax = plt.subplots(figsize=(6, 4.5)) @@ -469,8 +473,8 @@ def centerline_density_profile_figure(kn_label, fig_stem): tail = x_over_D >= 1.0 annotate_error( ax, - f'max |n_S/n_A - 1| = {max_rel_diff(simplified[tail], analytical[tail]):.2g} ' - f'(X/D $\\geq$ 1)\n' + f'Eq. 5 vs Eq. 18: max rel diff = ' + f'{max_rel_diff(simplified, analytical):.1g}\n' f'max |n_Sim/n_A - 1| = {max_rel_diff(simons[tail], analytical[tail]):.2g} ' f'(X/D $\\geq$ 1)', loc='lower right') return save_figure(fig, fig_stem) diff --git a/tests/plume/plume_verification_error_summary.py b/tests/plume/plume_verification_error_summary.py new file mode 100644 index 0000000..104d194 --- /dev/null +++ b/tests/plume/plume_verification_error_summary.py @@ -0,0 +1,117 @@ +# ======================== +# PyRPOD: tests/plume/plume_verification_error_summary.py +# ======================== +# Model-vs-model analog of Cai & Wang 2012 Table 1: maximum relative +# differences in density and mass flux between the analytical, +# simplified and Simons (exit-referenced, Boyton kappa) plume models +# along the centerline and along the angular curve r/D = 10, S0 = 2.0. +# +# Table 1's reference columns are DSMC results, which PyRPOD does not +# have yet; those columns are emitted as "pending digitized data" and +# will be filled from tests/plume/data/digitized/ in a follow-up. As +# the paper notes, "the largest density relative error actually +# happens within the curve tails with very small values" -- each row +# therefore also reports the restricted-range maximum. +# +# The filename intentionally avoids pytest's collection patterns +# (test_*.py / *_test_*.py): this is a manual-run generator -- +# python tests/plume/plume_verification_error_summary.py +# It writes tests/plume/output/model_error_summary.csv and .md. + +import numpy as np + +import plume_figure_utils as u + +S_0 = 2.0 + + +def _row(curve, quantity, comparison, values, reference, mask, mask_note): + full = u.max_rel_diff(values, reference) + restricted = u.max_rel_diff(values[mask], reference[mask]) + return {'curve': curve, 'quantity': quantity, 'comparison': comparison, + 'max_rel_diff': full, + 'max_rel_diff_restricted': restricted, + 'restriction': mask_note, + 'vs_DSMC': 'pending digitized data'} + + +def build_rows(): + rows = [] + + # centerline, X/D in [X_MIN, 10]. As in the paper's Figs. 19-21, + # "simplified" on the centerline is the Eq. 18 closed form (paper + # Sec. II.B) and "analytical" is the Eq. 5 integral; they coincide + # to quadrature accuracy. + x_over_D = np.linspace(u.X_MIN_OVER_D, 10.0, 200) + r = x_over_D * u.D_NOZZLE + n_analytical = u.full_field_values('n', r, np.zeros_like(r), S_0) + n_simplified = u.centerline_profiles(x_over_D, S_0, ('n',))['n'] + n_simons = u.simons_density_field(r, np.zeros_like(r), S_0) + core = x_over_D >= 1.0 + rows.append(_row('centerline', 'density', + 'simplified (Eq. 18) vs analytical (Eq. 5)', + n_simplified, n_analytical, core, 'X/D >= 1')) + rows.append(_row('centerline', 'density', 'Simons vs analytical', + n_simons, n_analytical, core, 'X/D >= 1')) + + # angular curve r/D = 10 + theta = np.deg2rad(np.linspace(0.0, u.THETA_MAX_DEG, 90)) + X = 10.0 * u.D_NOZZLE * np.cos(theta) + Z = 10.0 * u.D_NOZZLE * np.sin(theta) + nA = u.full_field_values('n', X, Z, S_0) + Q = X ** 2 / (X ** 2 + Z ** 2) + nS = u.get_K_factor(Q, S_0) / (2 * np.sqrt(np.pi)) * (u.R_0 / X) ** 2 + nSim = u.simons_density_field(X, Z, S_0) + VrA = u.full_field_values('Vr', X, Z, S_0) + fluxA = nA * VrA + UA = u.full_field_values('U', X, Z, S_0) + WA = u.full_field_values('W', X, Z, S_0) + # simplified-model flux with the same Fig. 25 normalization + US = u.get_M_factor(Q, S_0) / u.get_K_factor(Q, S_0) + WS = US * Z / X + VrS = (X * US + Z * WS) / np.sqrt(X ** 2 + Z ** 2) + fluxS = nS * VrS + core = theta <= np.deg2rad(60.0) + rows.append(_row('angular r/D=10', 'density', 'simplified vs analytical', + nS, nA, core, 'theta <= 60 deg')) + rows.append(_row('angular r/D=10', 'density', 'Simons vs analytical', + nSim, nA, core, 'theta <= 60 deg')) + rows.append(_row('angular r/D=10', 'mass flux', 'simplified vs analytical', + fluxS, fluxA, core, 'theta <= 60 deg')) + return rows + + +def write_summary(): + rows = build_rows() + u.OUTPUT_DIR.mkdir(parents=True, exist_ok=True) + columns = ['curve', 'quantity', 'comparison', 'max_rel_diff', + 'max_rel_diff_restricted', 'restriction', 'vs_DSMC'] + + csv_path = u.OUTPUT_DIR / 'model_error_summary.csv' + with open(csv_path, 'w', encoding='utf-8', newline='') as fh: + fh.write(','.join(columns) + '\n') + for row in rows: + fh.write(','.join( + f'{row[c]:.4g}' if isinstance(row[c], float) else str(row[c]) + for c in columns) + '\n') + + md_path = u.OUTPUT_DIR / 'model_error_summary.md' + with open(md_path, 'w', encoding='utf-8', newline='') as fh: + fh.write('# Model-vs-model error summary (Table 1 analog)\n\n') + fh.write(f'Cai & Wang 2012 conditions, S0 = {S_0}. DSMC reference ' + 'columns are pending digitized data.\n\n') + fh.write('| ' + ' | '.join(columns) + ' |\n') + fh.write('|' + '---|' * len(columns) + '\n') + for row in rows: + fh.write('| ' + ' | '.join( + f'{row[c]:.4g}' if isinstance(row[c], float) else str(row[c]) + for c in columns) + ' |\n') + fh.write('\nNote: unrestricted maxima occur in the plume tails ' + 'where the reference density is very small (the paper ' + 'makes the same observation about Table 1).\n') + return csv_path, md_path + + +if __name__ == '__main__': + for path in write_summary(): + print(f'wrote {path}') From c277c7022a05b931a3998340e2ea3ca45dab9ef9 Mon Sep 17 00:00:00 2001 From: andytorrestb Date: Sat, 18 Jul 2026 12:40:39 -0500 Subject: [PATCH 14/14] delegate figure overlays to plume model classes Co-Authored-By: Claude Fable 5 --- tests/plume/plume_figure_utils.py | 149 +++++++++++------- .../plume/plume_verification_error_summary.py | 18 +-- 2 files changed, 98 insertions(+), 69 deletions(-) diff --git a/tests/plume/plume_figure_utils.py b/tests/plume/plume_figure_utils.py index a7004e4..194772e 100644 --- a/tests/plume/plume_figure_utils.py +++ b/tests/plume/plume_figure_utils.py @@ -14,12 +14,17 @@ # T_0 and n_0 values are inert; the ones below match the existing # verification tests. ve is set from the exit speed ratio, # ve = S_0 * sqrt(2*R*T_0). -# * Full-model contour fields are evaluated with a vectorized -# Gauss-Legendre quadrature (fixed order 80 per axis, the same order -# at which CollisionlessGasKinetics converges to 1e-9). Every call -# cross-checks a few sample points against the class and raises if -# they disagree beyond 1e-6, so the vectorized path cannot silently -# diverge from the physics module. +# * Single-source-of-truth policy: the simplified and Simons overlays +# delegate point-by-point to SimplifiedGasKinetics / Simons (they are +# closed-form and cheap), so any model change is reflected in the +# figures automatically. Only the FULL model keeps a vectorized +# re-implementation of its Eq. 5-8 quadrature here (per-point class +# instantiation would cost hours per contour sweep instead of +# seconds); as a guard, every call cross-checks sample points against +# CollisionlessGasKinetics and raises if they disagree beyond 1e-6, +# so that path cannot silently diverge from the physics module +# either -- a model change makes the figure scripts fail loudly +# until this file is updated to match. # * The models require X > 0; grids start at X_MIN_OVER_D and angular # sweeps stop at THETA_MAX_DEG < 90 deg. # * Digitized-data convention (see tests/plume/data/digitized/README.md): @@ -103,13 +108,6 @@ def exit_mach(S_0): return S_0 * np.sqrt(2.0 / GAMMA) -def throat_to_exit_density_ratio(S_0): - """n_s/n_0 used by Simons exit referencing [module docstring note].""" - Me = exit_mach(S_0) - return ((1 + (GAMMA - 1) / 2 * Me ** 2) - / ((GAMMA + 1) / 2)) ** (1 / (GAMMA - 1)) - - def make_simons(distance, S_0, kappa=None): """Simons instance with chamber conditions isentropically consistent with the exit state at Me(S_0). T_c/P_c only rescale absolute @@ -122,27 +120,35 @@ def make_simons(distance, S_0, kappa=None): def simons_density_field(X, Z, S_0, kappa_f=None): - """Exit-referenced Simons density n/n_0 at points (X, 0, Z) [m]. - - A is always the Boyton-kappa normalization constant (paper p. 65: - the cosine-law curves coincide at theta = 0 for all kappa); kappa_f - optionally varies the plotted decay exponent f = cos^kappa_f - (Figs. 22-24). Vectorized; theta >= theta_max gives 0. + """Exit-referenced Simons density n/n_0 at points (X, 0, Z) [m], + delegated point-by-point to Simons.get_num_density_ratio_exit so + that any change to the Simons model is reflected here. + + Following the paper (p. 65: the cosine-law curves coincide at + theta = 0 for every kappa), the normalization constant A is always + the Boyton-kappa value; kappa_f optionally varies only the plotted + decay exponent f = cos^kappa_f (Figs. 22-24), implemented by + overriding A on a kappa_f-built instance with the Boyton value. + A single instance is reused with its evaluation radius updated per + point (reconstructing would re-run the A quadrature 10^4 times). """ X = np.asarray(X, dtype=float) Z = np.asarray(Z, dtype=float) - boyton = make_simons(1.0, S_0) # A from default Boyton kappa - theta_max = boyton.get_limiting_turn_angle() - kappa = boyton.kappa if kappa_f is None else kappa_f + shape = np.broadcast(X, Z).shape + Xf = np.broadcast_to(X, shape).ravel() + Zf = np.broadcast_to(Z, shape).ravel() + + simons = make_simons(1.0, S_0, kappa=kappa_f) + if kappa_f is not None: + simons.A = make_simons(1.0, S_0).A # paper's shared Boyton A + Me = exit_mach(S_0) - r_sph = np.sqrt(X ** 2 + Z ** 2) - theta = np.arctan2(np.abs(Z), X) - f = np.where(theta < theta_max, - np.cos((np.pi / 2) * (np.minimum(theta, theta_max) - / theta_max)) ** kappa, - 0.0) - n_over_ns = boyton.A * (R_0 / r_sph) ** 2 * f - return n_over_ns * throat_to_exit_density_ratio(S_0) + out = np.empty(Xf.size) + for i in range(Xf.size): + simons.r = float(np.hypot(Xf[i], Zf[i])) + theta = float(np.arctan2(np.abs(Zf[i]), Xf[i])) + out[i] = simons.get_num_density_ratio_exit(theta, Me) + return out.reshape(shape) # --------------------------------------------------------------------------- @@ -163,10 +169,40 @@ def centerline_profiles(x_over_D, S_0, quantities=('n',)): return out -def simplified_centerline_density(x_over_D, S_0): - """Simplified-model centerline density, Eq. 14 with Q' = 1.""" - X = np.asarray(x_over_D) * D_NOZZLE - return get_K_factor(1.0, S_0) / (2 * np.sqrt(np.pi)) * (R_0 / X) ** 2 +# --------------------------------------------------------------------------- +# simplified-model field evaluation (delegated to SimplifiedGasKinetics) +# --------------------------------------------------------------------------- + +def simplified_field_values(quantity, X, Z, S_0): + """Simplified-model quantity at 1-D point arrays X, Z [m], + delegated point-by-point to SimplifiedGasKinetics (Eqs. 13-17) so + that any change to the simplified model is reflected here. The + class is closed-form, so the per-point loop costs well under a + second even on full contour grids. + + quantity: 'n', 'U', 'W', 'T', plus the derived compositions + 'Vr' = (X*U + Z*W)/sqrt(X^2+Z^2) and 'p' = n*T (the class exposes + no methods for those). + """ + getters = {'n': 'get_num_density_ratio', 'U': 'get_U_normalized', + 'W': 'get_W_normalized', 'T': 'get_temp_ratio'} + if quantity not in getters and quantity not in ('Vr', 'p'): + raise KeyError(quantity) + X = np.asarray(X, dtype=float).ravel() + Z = np.asarray(Z, dtype=float).ravel() + out = np.empty(X.size) + for i in range(X.size): + d = float(np.hypot(X[i], Z[i])) + theta = float(np.arctan2(Z[i], X[i])) + plume = make_plume(SimplifiedGasKinetics, d, theta, S_0) + if quantity in getters: + out[i] = getattr(plume, getters[quantity])() + elif quantity == 'Vr': + out[i] = (X[i] * plume.get_U_normalized() + + Z[i] * plume.get_W_normalized()) / d + else: # 'p' + out[i] = plume.get_num_density_ratio() * plume.get_temp_ratio() + return out # --------------------------------------------------------------------------- @@ -201,7 +237,7 @@ def full_field_values(quantity, X, Z, S_0, _chunk=256): sl = slice(start, min(start + _chunk, n_pts)) Xc = X[sl][:, None, None] Zc = Z[sl][:, None, None] - Q = Xc ** 2 / (Xc ** 2 + Zc ** 2 - 2 * Zc * Rn * sinE + Rn ** 2) + Q = get_Q_full(Rn, En, Xc, Zc) K = get_K_factor(Q, S_0) M = get_M_factor(Q, S_0) N = get_N_factor(Q, S_0) @@ -222,17 +258,26 @@ def full_field_values(quantity, X, Z, S_0, _chunk=256): def _verify_against_class(values, X, Z, S_0, rtol=1e-6): - """Cross-check sample points against CollisionlessGasKinetics.""" + """Cross-check sample points against CollisionlessGasKinetics. + + W legitimately vanishes on the centerline, so its deviation is + measured against the natural velocity scale max(|W|, |U|) instead + of a purely relative test. + """ for idx in {0, X.size // 2, X.size - 1}: d = float(np.hypot(X[idx], Z[idx])) th = float(np.arctan2(Z[idx], X[idx])) ref = make_plume(CollisionlessGasKinetics, d, th, S_0) - checks = {'n': ref.get_num_density_ratio(), - 'U': ref.get_U_normalized(), - 'T': ref.get_temp_ratio()} - for q, expected in checks.items(): + U_ref = ref.get_U_normalized() + W_ref = ref.get_W_normalized() + checks = {'n': (ref.get_num_density_ratio(), None), + 'U': (U_ref, None), + 'W': (W_ref, max(abs(W_ref), abs(U_ref))), + 'T': (ref.get_temp_ratio(), None)} + for q, (expected, scale) in checks.items(): + scale = abs(expected) if scale is None else scale got = values[q][idx] - if abs(got - expected) > rtol * abs(expected): + if abs(got - expected) > rtol * scale: raise AssertionError( f'vectorized field {q} diverged from ' f'CollisionlessGasKinetics at (X={X[idx]}, Z={Z[idx]}): ' @@ -249,21 +294,13 @@ def full_field_grid(quantity, x_over_D, z_over_D, S_0): def simplified_field_grid(quantity, x_over_D, z_over_D, S_0): - """Simplified-model quantity (Eqs. 13-17) on the tensor grid.""" + """Simplified-model quantity on the tensor grid (delegates to + simplified_field_values); returns shape + (len(z_over_D), len(x_over_D)) for use with plt.contour.""" Xg, Zg = np.meshgrid(np.asarray(x_over_D) * D_NOZZLE, np.asarray(z_over_D) * D_NOZZLE) - Q = Xg ** 2 / (Xg ** 2 + Zg ** 2) - K = get_K_factor(Q, S_0) - M = get_M_factor(Q, S_0) - N = get_N_factor(Q, S_0) - n = K / (2 * np.sqrt(np.pi)) * (R_0 / Xg) ** 2 - U = M / K - W = U * Zg / Xg - T = -2 * M ** 2 / (3 * Q * K ** 2) + 4 * N / (3 * K) - values = {'n': n, 'U': U, 'W': W, 'T': T, - 'Vr': (Xg * U + Zg * W) / np.sqrt(Xg ** 2 + Zg ** 2), - 'p': n * T} - return values[quantity] + vals = simplified_field_values(quantity, Xg.ravel(), Zg.ravel(), S_0) + return vals.reshape(Xg.shape) def default_contour_axes(x_max=10.0, z_max=10.0): @@ -490,9 +527,7 @@ def angular_density_profile_figure(kn_label, fig_stem): X = r * np.cos(theta) Z = r * np.sin(theta) analytical = full_field_values('n', X, Z, S_0) - Q = X ** 2 / (X ** 2 + Z ** 2) - simplified = (get_K_factor(Q, S_0) / (2 * np.sqrt(np.pi)) - * (R_0 / X) ** 2) + simplified = simplified_field_values('n', X, Z, S_0) theta_deg = np.rad2deg(theta) fig, ax = plt.subplots(figsize=(6, 4.5)) diff --git a/tests/plume/plume_verification_error_summary.py b/tests/plume/plume_verification_error_summary.py index 104d194..b99b6d7 100644 --- a/tests/plume/plume_verification_error_summary.py +++ b/tests/plume/plume_verification_error_summary.py @@ -54,23 +54,17 @@ def build_rows(): rows.append(_row('centerline', 'density', 'Simons vs analytical', n_simons, n_analytical, core, 'X/D >= 1')) - # angular curve r/D = 10 + # angular curve r/D = 10 (simplified series delegates to + # SimplifiedGasKinetics via plume_figure_utils) theta = np.deg2rad(np.linspace(0.0, u.THETA_MAX_DEG, 90)) X = 10.0 * u.D_NOZZLE * np.cos(theta) Z = 10.0 * u.D_NOZZLE * np.sin(theta) nA = u.full_field_values('n', X, Z, S_0) - Q = X ** 2 / (X ** 2 + Z ** 2) - nS = u.get_K_factor(Q, S_0) / (2 * np.sqrt(np.pi)) * (u.R_0 / X) ** 2 + nS = u.simplified_field_values('n', X, Z, S_0) nSim = u.simons_density_field(X, Z, S_0) - VrA = u.full_field_values('Vr', X, Z, S_0) - fluxA = nA * VrA - UA = u.full_field_values('U', X, Z, S_0) - WA = u.full_field_values('W', X, Z, S_0) - # simplified-model flux with the same Fig. 25 normalization - US = u.get_M_factor(Q, S_0) / u.get_K_factor(Q, S_0) - WS = US * Z / X - VrS = (X * US + Z * WS) / np.sqrt(X ** 2 + Z ** 2) - fluxS = nS * VrS + # mass flux with the same Fig. 25 normalization, n * Vr*sqrt(beta0) + fluxA = nA * u.full_field_values('Vr', X, Z, S_0) + fluxS = nS * u.simplified_field_values('Vr', X, Z, S_0) core = theta <= np.deg2rad(60.0) rows.append(_row('angular r/D=10', 'density', 'simplified vs analytical', nS, nA, core, 'theta <= 60 deg'))