From a23cf8e2fa5b8d2d2a7be0e8c7e88285c5421c5e Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sun, 26 Mar 2023 09:48:40 -0700 Subject: [PATCH 1/2] Try different handling of low density Use Nelim to avoid dividing by small values. Force pressure to be consistent with the limited density and temperature. --- sd1d.cxx | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/sd1d.cxx b/sd1d.cxx index eb0fd70..547cb0e 100644 --- a/sd1d.cxx +++ b/sd1d.cxx @@ -382,6 +382,8 @@ class SD1D : public PhysicsModel { mesh->communicate(Ne, NVi, P); + Field3D P_solver = P; // Save the pressure in the solver + // Floor small values P = floor(P, 1e-10); Ne = floor(Ne, 1e-10); @@ -390,7 +392,7 @@ class SD1D : public PhysicsModel { Vi = NVi / Ne; - Field3D Te = 0.5 * P / Ne; // Assuming Te = Ti + Field3D Te = 0.5 * P / Nelim; // Assuming Te = Ti for (auto &i : Te.getRegion("RGN_NOBNDRY")) { if (Te[i] > 10.) @@ -1271,6 +1273,10 @@ class SD1D : public PhysicsModel { if (ADpar > 0.0) { ddt(P) += ADpar * AddedDissipation(1.0, P, P, true); } + + // Force P towards value in solver. + // Note: This only impacts solutions where Ne < floor + ddt(P) += 2. * Te * Ne - P_solver; } } From d9293838722e8305f955a94a43b1c75196266634 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sun, 26 Mar 2023 09:54:00 -0700 Subject: [PATCH 2/2] Allow numerical momentum flux through boundary Reduces oscillations close to the boundary. --- sd1d.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sd1d.cxx b/sd1d.cxx index 547cb0e..2afbebf 100644 --- a/sd1d.cxx +++ b/sd1d.cxx @@ -1159,7 +1159,7 @@ class SD1D : public PhysicsModel { if (rhs_explicit) { // Flux splitting with upwinding Field3D a = sqrt(gamma_sound * 2. * Te); // Local sound speed - ddt(NVi) = -FV::Div_par(NVi, Vi, a, bndry_flux_fix) // Momentum flow + ddt(NVi) = -FV::Div_par(NVi, Vi, a, false) // Momentum flow - Grad_par(P); if (atomic) {