diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs index 5ea4e3c5d..356a81dfd 100644 --- a/.git-blame-ignore-revs +++ b/.git-blame-ignore-revs @@ -2,3 +2,4 @@ f5de620399e88f139e236e40675facfd0cfd6b13 07616298519738fad0c1bb259e4d021cd6f51d58 37efa12abe0961add8715b797d3926664977d453 8e49d5933131ddcb2f97bb2ac7095a1be80ad7b7 +eaf52851d149ca10f69371e3c78e413f2e3e60cf diff --git a/include/braginskii_collisions.hxx b/include/braginskii_collisions.hxx index 5303d72c6..bc3e9d0db 100644 --- a/include/braginskii_collisions.hxx +++ b/include/braginskii_collisions.hxx @@ -57,6 +57,7 @@ private: neutral_neutral; BoutReal ei_multiplier; // Arbitrary user-set multiplier on electron-ion collisions + BoutReal density_floor; /// Calculated collision rates saved for post-processing and use by other components /// Saved in options, the BOUT++ dictionary-like object diff --git a/include/braginskii_electron_viscosity.hxx b/include/braginskii_electron_viscosity.hxx index ccb57feaa..19300cfa7 100644 --- a/include/braginskii_electron_viscosity.hxx +++ b/include/braginskii_electron_viscosity.hxx @@ -37,6 +37,9 @@ private: bool diagnose; ///< Output viscosity diagnostic? Field3D viscosity; ///< The viscosity momentum source + BoutReal density_floor; + BoutReal pressure_floor; ///< When non-zero pressure is needed + /// Inputs /// - species /// - e diff --git a/include/braginskii_ion_viscosity.hxx b/include/braginskii_ion_viscosity.hxx index 242759959..4ea894942 100644 --- a/include/braginskii_ion_viscosity.hxx +++ b/include/braginskii_ion_viscosity.hxx @@ -61,6 +61,7 @@ private: ///< frequency change BoutReal bounce_frequency_R; ///< Input major radius BoutReal density_floor; ///< Minimum density used in calculating Pi_ciperp + BoutReal pressure_floor; ///< When non-zero pressure is needed bool diagnose; ///< Output additional diagnostics? /// Per-species diagnostics diff --git a/src/braginskii_collisions.cxx b/src/braginskii_collisions.cxx index d66716e69..9f8d4d971 100644 --- a/src/braginskii_collisions.cxx +++ b/src/braginskii_collisions.cxx @@ -59,6 +59,8 @@ BraginskiiCollisions::BraginskiiCollisions(const std::string& name, Options& all .doc("User-set arbitrary multiplier on electron-ion collision rate") .withDefault(1.0); + density_floor = options["density_floor"].doc("Minimum density floor").withDefault(1e-8); + diagnose = options["diagnose"].doc("Output additional diagnostics?").withDefault(false); @@ -144,7 +146,8 @@ void BraginskiiCollisions::collide(GuardedOptions& species1, GuardedOptions& spe const Field3D density2 = GET_NOBOUNDARY(Field3D, species2["density"]); const Field3D nu = filledFrom(nu_12, [&](auto& i) { - return nu_12[i] * (A1 / A2) * density1[i] / softFloor(density2[i], 1e-5); + return nu_12[i] * (A1 / A2) * softFloor(density1[i], density_floor) + / softFloor(density2[i], density_floor); }); add(species2["collision_frequency"], nu); // Total collision frequency diff --git a/src/braginskii_electron_viscosity.cxx b/src/braginskii_electron_viscosity.cxx index ac9269782..9331803e1 100644 --- a/src/braginskii_electron_viscosity.cxx +++ b/src/braginskii_electron_viscosity.cxx @@ -29,6 +29,14 @@ BraginskiiElectronViscosity::BraginskiiElectronViscosity(const std::string& name .doc("Viscosity flux limiter coefficient. <0 = turned off") .withDefault(-1.0); + density_floor = options["density_floor"].doc("Minimum density floor").withDefault(1e-8); + + BoutReal temperature_floor = + options["temperature_floor"].doc("Low temperature scale").withDefault(0.1) + / get(alloptions["units"]["eV"]); + + pressure_floor = density_floor * temperature_floor; + diagnose = options["diagnose"].doc("Output diagnostics?").withDefault(false); } @@ -44,8 +52,8 @@ void BraginskiiElectronViscosity::transform_impl(GuardedOptions& state) { throw BoutException("No electron velocity => Can't calculate electron viscosity"); } - const Field3D tau = 1. / get(species["collision_frequency"]); - const Field3D P = get(species["pressure"]); + const Field3D tau = 1. / softFloor(get(species["collision_frequency"]), 1e-15); + const Field3D P = floor(get(species["pressure"]), pressure_floor); const Field3D V = get(species["velocity"]); Coordinates* coord = P.getCoordinates(); @@ -59,10 +67,10 @@ void BraginskiiElectronViscosity::transform_impl(GuardedOptions& state) { // SOLPS-style flux limiter // Values of alpha ~ 0.5 typically - const Field3D q_cl = eta * Grad_par(V); // Collisional value + const Field3D q_cl = eta * abs(Grad_par(V)); // Collisional value const Field3D q_fl = eta_limit_alpha * P; // Flux limit - eta = eta / (1. + abs(q_cl / q_fl)); + eta = eta / (1. + softFloor(q_cl, 1e-15) / softFloor(q_fl, 1e-15)); eta.getMesh()->communicate(eta); eta.applyBoundary("neumann"); diff --git a/src/braginskii_friction.cxx b/src/braginskii_friction.cxx index fc76f00a5..0f261a502 100644 --- a/src/braginskii_friction.cxx +++ b/src/braginskii_friction.cxx @@ -105,7 +105,6 @@ void BraginskiiFriction::transform_impl(GuardedOptions& state) { } const Field3D nu = GET_VALUE(Field3D, species1["collision_frequencies"][coll_name]); - const Field3D density2 = GET_VALUE(Field3D, species2["density"]); const Field3D velocity2 = species2.isSet("velocity") ? GET_NOBOUNDARY(Field3D, species2["velocity"]) : 0.0; diff --git a/src/braginskii_ion_viscosity.cxx b/src/braginskii_ion_viscosity.cxx index 9578d1af5..6a66e38bf 100644 --- a/src/braginskii_ion_viscosity.cxx +++ b/src/braginskii_ion_viscosity.cxx @@ -48,6 +48,14 @@ BraginskiiIonViscosity::BraginskiiIonViscosity(const std::string& name, .doc("Viscosity flux limiter coefficient. <0 = turned off") .withDefault(-1.0); + density_floor = options["density_floor"].doc("Minimum density floor").withDefault(1e-8); + + BoutReal temperature_floor = + options["temperature_floor"].doc("Low temperature scale").withDefault(0.1) + / get(alloptions["units"]["eV"]); + + pressure_floor = density_floor * temperature_floor; + diagnose = options["diagnose"].doc("Output additional diagnostics?").withDefault(false); @@ -84,8 +92,6 @@ BraginskiiIonViscosity::BraginskiiIonViscosity(const std::string& name, "frequency modification to viscosity") .withDefault(2.0); - density_floor = options["density_floor"].doc("Minimum density floor").withDefault(1e-7); - if (perpendicular) { // Read curvature vector try { @@ -230,7 +236,7 @@ void BraginskiiIonViscosity::transform_impl(GuardedOptions& state) { } const Field3D tau = 1. / nu; - const Field3D P = get(species["pressure"]); + const Field3D P = floor(get(species["pressure"]), pressure_floor); const Field3D T = get(species["temperature"]); const Field3D N = get(species["density"]); @@ -272,10 +278,10 @@ void BraginskiiIonViscosity::transform_impl(GuardedOptions& state) { // SOLPS-style flux limiter // Values of alpha ~ 0.5 typically - const Field3D q_cl = eta * Grad_par(V); // Collisional value + const Field3D q_cl = eta * abs(Grad_par(V)); // Collisional value const Field3D q_fl = eta_limit_alpha * P; // Flux limit - eta = eta / (1. + abs(q_cl / q_fl)); + eta = eta / (1. + (q_cl / q_fl)); eta.getMesh()->communicate(eta); eta.applyBoundary("neumann"); diff --git a/src/braginskii_thermal_force.cxx b/src/braginskii_thermal_force.cxx index 99f1c238b..98910c006 100644 --- a/src/braginskii_thermal_force.cxx +++ b/src/braginskii_thermal_force.cxx @@ -38,7 +38,7 @@ void BraginskiiThermalForce::transform_impl(GuardedOptions& state) { const BoutReal Z = get(species["charge"]); // Don't need density boundary - const Field3D nz = GET_NOBOUNDARY(Field3D, species["density"]); + const Field3D nz = floor(GET_NOBOUNDARY(Field3D, species["density"]), 0.0); const Field3D ion_force = nz * (0.71 * SQ(Z)) * Grad_Te; @@ -123,7 +123,8 @@ void BraginskiiThermalForce::transform_impl(GuardedOptions& state) { const BoutReal mz = get((*heavy)["AA"]); const BoutReal Z = get((*heavy)["charge"]); - const Field3D nz = GET_NOBOUNDARY(Field3D, (*heavy)["density"]); + const Field3D nz = floor(GET_NOBOUNDARY(Field3D, (*heavy)["density"]), 0.0); + ; if (Z == 0.0) { continue; // Check that the charge is not zero diff --git a/src/hydrogen_charge_exchange.cxx b/src/hydrogen_charge_exchange.cxx index 44cca5f8d..e7ed7bfdb 100644 --- a/src/hydrogen_charge_exchange.cxx +++ b/src/hydrogen_charge_exchange.cxx @@ -20,7 +20,7 @@ void HydrogenChargeExchange::calculate_rates(GuardedOptions&& atom1, GuardedOpti // Calculate effective temperature in eV Field3D Teff = (Tatom / Aatom + Tion / Aion) * Tnorm; - for (auto& i : Teff.getRegion("RGN_NOBNDRY")) { + for (auto& i : Teff.getRegion("RGN_ALL")) { if (Teff[i] < 0.01) { Teff[i] = 0.01; } else if (Teff[i] > 10000) { @@ -44,8 +44,8 @@ void HydrogenChargeExchange::calculate_rates(GuardedOptions&& atom1, GuardedOpti // Optionally multiply by arbitrary multiplier const Field3D sigmav = exp(ln_sigmav) * (1e-6 * Nnorm / FreqNorm) * rate_multiplier; - const Field3D Natom = get(atom1["density"]); - const Field3D Nion = get(ion1["density"]); + const Field3D Natom = softFloor(get(atom1["density"]), 0.0); + const Field3D Nion = softFloor(get(ion1["density"]), 0.0); R = Natom * Nion * sigmav; // Rate coefficient in [m^-3 s^-1] diff --git a/src/relax_potential.cxx b/src/relax_potential.cxx index 2eb404d6c..1ba8658fe 100644 --- a/src/relax_potential.cxx +++ b/src/relax_potential.cxx @@ -7,6 +7,7 @@ #include #include #include +#include using bout::globals::mesh; @@ -843,6 +844,20 @@ void RelaxPotential::finally(const Options& state) { ddt(phi1) = lambda_1 * (phi_vort - Vort); } + +#if CHECKLEVEL >= 1 + for (auto& i : Vort.getRegion("RGN_NOBNDRY")) { + if (!std::isfinite(ddt(Vort)[i])) { + throw BoutException("ddt(Vort) non-finite at {}\n", i); + } + } + + for (auto& i : phi1.getRegion("RGN_NOBNDRY")) { + if (!std::isfinite(ddt(phi1)[i])) { + throw BoutException("ddt(phi1) non-finite at {}\n", i); + } + } +#endif } void RelaxPotential::outputVars(Options& state) { diff --git a/src/sheath_boundary_simple.cxx b/src/sheath_boundary_simple.cxx index ebe666925..83047b54a 100644 --- a/src/sheath_boundary_simple.cxx +++ b/src/sheath_boundary_simple.cxx @@ -43,9 +43,6 @@ BoutReal limitFree(BoutReal fm, BoutReal fc, BoutReal mode) { throw BoutException("Unknown boundary mode"); } - return fp; // Extrapolation - - #if CHECKLEVEL >= 2 if (!std::isfinite(fp)) { throw BoutException("SheathBoundary limitFree: {}, {} -> {}", fm, fc, fp); diff --git a/src/vorticity.cxx b/src/vorticity.cxx index 116fb4c22..fb28ecdc8 100644 --- a/src/vorticity.cxx +++ b/src/vorticity.cxx @@ -4,11 +4,12 @@ #include "../include/hermes_utils.hxx" #include -#include -#include #include #include +#include +#include #include +#include using bout::globals::mesh; @@ -901,6 +902,14 @@ void Vorticity::finally(const Options& state) { } } } + +#if CHECKLEVEL >= 1 + for (auto& i : Vort.getRegion("RGN_NOBNDRY")) { + if (!std::isfinite(ddt(Vort)[i])) { + throw BoutException("ddt(Vort) non-finite at {}\n", i); + } + } +#endif } void Vorticity::outputVars(Options& state) {