From a59d6ea14d62c56f160f41a5fddad8a4dfa21b5e Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 11 Feb 2026 16:12:44 -0800 Subject: [PATCH 01/29] BOUT-dev: Use branch with SNES solver updates --- external/BOUT-dev | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/BOUT-dev b/external/BOUT-dev index 7d28d67c3..b9b038b02 160000 --- a/external/BOUT-dev +++ b/external/BOUT-dev @@ -1 +1 @@ -Subproject commit 7d28d67c3f12c24ec281c0982e870f5369c65a6f +Subproject commit b9b038b0243d3a54e8f89448fb0eb6664f7821d9 From 40cdfbdc73b702a2767f4936c053e08e99a02194 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 12 Feb 2026 11:02:02 -0800 Subject: [PATCH 02/29] BOUT-dev: More SNES changes Added jacobian_persists option --- external/BOUT-dev | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/BOUT-dev b/external/BOUT-dev index b9b038b02..74ce9aad8 160000 --- a/external/BOUT-dev +++ b/external/BOUT-dev @@ -1 +1 @@ -Subproject commit b9b038b0243d3a54e8f89448fb0eb6664f7821d9 +Subproject commit 74ce9aad802cf355ce7fdabab55e7c90cd0f5c9b From ffcb5c85d048f6e3331287de87b28eef5a3e8d52 Mon Sep 17 00:00:00 2001 From: malamast Date: Mon, 29 Dec 2025 15:57:22 -0800 Subject: [PATCH 03/29] evolve_density and pressure: Added a functionality to initialize number density and pressure from grid variables with names varName_init. fixed_density and isothermal: Added a functionality to read the fixed denstiy or temperature from grid variables with names varName_init --- include/evolve_density.hxx | 3 +++ include/evolve_pressure.hxx | 2 ++ include/fixed_density.hxx | 17 +++++++++++++++++ include/isothermal.hxx | 12 ++++++++++++ src/evolve_density.cxx | 13 +++++++++++++ src/evolve_pressure.cxx | 15 +++++++++++++-- src/isothermal.cxx | 17 ++++++++++++++++- 7 files changed, 76 insertions(+), 3 deletions(-) diff --git a/include/evolve_density.hxx b/include/evolve_density.hxx index 273d1999d..3459e0b03 100644 --- a/include/evolve_density.hxx +++ b/include/evolve_density.hxx @@ -86,6 +86,9 @@ private: bool diagnose; ///< Output additional diagnostics? Field3D flow_xlow, flow_ylow; ///< Particle flow diagnostics + bool initialize_from_mesh; ///< Initilize the Field3D N from 2D profiles stored in the mesh file. + + /// This sets in the state /// - species /// - diff --git a/include/evolve_pressure.hxx b/include/evolve_pressure.hxx index 32b2ef9ba..4bd1a2ef3 100644 --- a/include/evolve_pressure.hxx +++ b/include/evolve_pressure.hxx @@ -100,6 +100,8 @@ private: bool fix_momentum_boundary_flux; ///< Fix momentum flux to boundary condition? Field3D Sp_nvh; ///< Pressure source due to artificial viscosity Field3D E_PdivV, E_VgradP; ///< Diagnostic energy source terms for p*Div(V) and V*Grad(P) + bool initialize_from_mesh; ///< Initilize the Field3D P from 2D profiles stored in the mesh file. + /// Inputs /// - species diff --git a/include/fixed_density.hxx b/include/fixed_density.hxx index c0e910281..1ca183ff5 100644 --- a/include/fixed_density.hxx +++ b/include/fixed_density.hxx @@ -4,6 +4,8 @@ #include "component.hxx" +using bout::globals::mesh; + /// Set ion density to a fixed value /// struct FixedDensity : public Component { @@ -26,6 +28,18 @@ struct FixedDensity : public Component { // Get the density and normalise N = options["density"].as() / Nnorm; + + // Initilize the Field3D N from 2D profiles stored in the mesh file. + initialize_from_mesh = options["initialize_from_mesh"] + .doc("Initilize field from 2D profiles stored in the mesh file?") + .withDefault(false); + + if (initialize_from_mesh) { + // Try to read the initial field from the mesh + mesh->get(N, std::string("N") + name + "_init"); // Units: [m^-3/s] + N /= Nnorm; // Normalization + } + substitutePermissions("name", {name}); substitutePermissions("vars", {"AA", "charge", "density"}); } @@ -51,6 +65,9 @@ private: Field3D N; ///< Species density (normalised) + bool initialize_from_mesh; ///< Initilize the Field3D N from 2D profiles stored in the mesh file. + + /// Sets in the state the density, mass and charge of the species /// /// - species diff --git a/include/isothermal.hxx b/include/isothermal.hxx index 1b385c5bd..54d322235 100644 --- a/include/isothermal.hxx +++ b/include/isothermal.hxx @@ -31,6 +31,18 @@ private: /// - pressure (if density is set) /// void transform_impl(GuardedOptions& state) override; + + void outputVars(Options &state) override; +private: + std::string name; // Species name + + // BoutReal T; ///< The normalised temperature + Field3D T; ///< The normalised temperature + Field3D P; ///< The normalised pressure + + bool initialize_from_mesh; ///< Initilize the Field3D T from 2D profiles stored in the mesh file. + + bool diagnose; ///< Output additional diagnostics? }; namespace { diff --git a/src/evolve_density.cxx b/src/evolve_density.cxx index 6efd3a8f3..f22f5f5f6 100644 --- a/src/evolve_density.cxx +++ b/src/evolve_density.cxx @@ -143,6 +143,19 @@ EvolveDensity::EvolveDensity(std::string name, Options& alloptions, Solver* solv } substitutePermissions("name", {name}); substitutePermissions("outputs", outputs); + + // Initilize the Field3D N from 2D profiles stored in the mesh file. + initialize_from_mesh = n_options["initialize_from_mesh"] + .doc("Initilize field from 2D profiles stored in the mesh file?") + .withDefault(false); + + if (initialize_from_mesh) { + // Try to read the initial field from the mesh + mesh->get(N, std::string("N") + name + "_init"); // Units: [m^-3/s] + N /= Nnorm; // Normalization + + } + } void EvolveDensity::transform_impl(GuardedOptions& state) { diff --git a/src/evolve_pressure.cxx b/src/evolve_pressure.cxx index ebfb606f8..5b81e705a 100644 --- a/src/evolve_pressure.cxx +++ b/src/evolve_pressure.cxx @@ -100,10 +100,10 @@ EvolvePressure::EvolvePressure(std::string name, Options& alloptions, Solver* so const BoutReal Nnorm = units["inv_meters_cubed"]; const BoutReal Tnorm = units["eV"]; const BoutReal Omega_ci = 1. / units["seconds"].as(); + const BoutReal Pnorm = SI::qe * Tnorm * Nnorm; // Pressure normalisation auto& p_options = alloptions[std::string("P") + name]; - source_normalisation = - SI::qe * Nnorm * Tnorm * Omega_ci; // [Pa/s] or [W/m^3] if converted to energy + source_normalisation = Pnorm * Omega_ci; // [Pa/s] or [W/m^3] if converted to energy time_normalisation = 1. / Omega_ci; // [s] // Try to read the pressure source from the mesh @@ -164,6 +164,17 @@ EvolvePressure::EvolvePressure(std::string name, Options& alloptions, Solver* so .doc("Include parallel heat conduction?") .withDefault(true); + // Initilize the Field3D P from 2D profiles stored in the mesh file. + initialize_from_mesh = p_options["initialize_from_mesh"] + .doc("Initilize field from 2D profiles stored in the mesh file?") + .withDefault(false); + + if (initialize_from_mesh) { + // Try to read the initial field from the mesh + mesh->get(P, std::string("P") + name + "_init"); // Units: Pascal [Pa] + P /= Pnorm; // Normalization + } + if (source_time_dependent) { setPermissions(readOnly("time")); } diff --git a/src/isothermal.cxx b/src/isothermal.cxx index eaa47c480..7889e596e 100644 --- a/src/isothermal.cxx +++ b/src/isothermal.cxx @@ -1,8 +1,12 @@ #include - +#include #include "../include/isothermal.hxx" +using bout::globals::mesh; + +using bout::globals::mesh; + Isothermal::Isothermal(std::string name, Options& alloptions, Solver* UNUSED(solver)) : Component({readIfSet(fmt::format("species:{}:density", name), Regions::Interior), readWrite(fmt::format("species:{}:temperature", name)), @@ -15,6 +19,17 @@ Isothermal::Isothermal(std::string name, Options& alloptions, Solver* UNUSED(sol T = options["temperature"].doc("Constant temperature [eV]").as() / Tnorm; // Normalise + // Initilize the Field3D T from 2D profiles stored in the mesh file. + initialize_from_mesh = options["initialize_from_mesh"] + .doc("Initilize field from 2D profiles stored in the mesh file?") + .withDefault(false); + + if (initialize_from_mesh) { + // Try to read the initial field from the mesh + mesh->get(T, std::string("T") + name + "_init"); // Units: [eV] + T /= Tnorm; // Normalization + } + diagnose = options["diagnose"] .doc("Save additional output diagnostics") .withDefault(false); From a165e565d5825cda62bc7943976906045525ce63 Mon Sep 17 00:00:00 2001 From: malamast Date: Mon, 29 Dec 2025 16:49:04 -0800 Subject: [PATCH 04/29] zero_current: Set species momentum and added momentum to outputVars --- include/zero_current.hxx | 16 ++++++++++++++++ src/zero_current.cxx | 21 ++++++++++++++++++++- 2 files changed, 36 insertions(+), 1 deletion(-) diff --git a/include/zero_current.hxx b/include/zero_current.hxx index 5a00eb758..1d6279e11 100644 --- a/include/zero_current.hxx +++ b/include/zero_current.hxx @@ -49,6 +49,22 @@ private: /// - velocity /// void transform_impl(GuardedOptions& state) override; + + void finally(const Options &state) override { + // Get the velocity with boundary condition applied. + // This is for output only + velocity = get(state["species"][name]["velocity"]); + } + + void outputVars(Options &state) override; +private: + std::string name; ///< Name of this species + BoutReal charge; ///< The charge of this species + BoutReal atomic_mass; // Atomic mass + + Field3D velocity; ///< Species velocity (for writing to output) + Field3D momentum; ///< Species momentum (for writing to output) + }; namespace { diff --git a/src/zero_current.cxx b/src/zero_current.cxx index 4d1a9e887..12242af54 100644 --- a/src/zero_current.cxx +++ b/src/zero_current.cxx @@ -2,6 +2,8 @@ #include #include "../include/zero_current.hxx" +#include +#include "../include/hermes_utils.hxx" ZeroCurrent::ZeroCurrent(std::string name, Options& alloptions, Solver*) : Component({readIfSet("species:{all_species}:charge"), @@ -14,6 +16,9 @@ ZeroCurrent::ZeroCurrent(std::string name, Options& alloptions, Solver*) ASSERT0(charge != 0.0); + atomic_mass = options["AA"].doc("Particle atomic number.").withDefault(0.0); // Atomic mass + // ASSERT0(atomic_mass != 0.0); + substitutePermissions("inputs", {"density", "velocity"}); } @@ -64,12 +69,16 @@ void ZeroCurrent::transform_impl(GuardedOptions& state) { } Field3D N = getNoBoundary(species["density"]); - velocity = current / (-charge * floor(N, 1e-5)); + velocity = current / (-charge * softFloor(N, 1e-7)); set(species["velocity"], velocity); + + momentum = atomic_mass * N * velocity; + set(species["momentum"], momentum); } void ZeroCurrent::outputVars(Options& state) { auto Cs0 = get(state["Cs0"]); + auto Nnorm = get(state["Nnorm"]); // Save the velocity set_with_attrs(state[std::string("V") + name], velocity, @@ -80,4 +89,14 @@ void ZeroCurrent::outputVars(Options& state) { {"standard_name", "velocity"}, {"species", name}, {"source", "zero_current"}}); + + set_with_attrs(state[std::string("NV") + name], momentum, + {{"time_dimension", "t"}, + {"units", "kg / m^2 / s"}, + {"conversion", SI::Mp * Nnorm * Cs0}, + {"long_name", name + " parallel momentum"}, + {"standard_name", "momentum"}, + {"species", name}, + {"source", "zero_current"}}); + } From 1334bb6dadda35951e59be9a7eabfc473f52f830 Mon Sep 17 00:00:00 2001 From: malamast Date: Mon, 29 Dec 2025 16:55:15 -0800 Subject: [PATCH 05/29] anomalous_diffusion: make anomalous_D variable accessible to other components. This willbe needed in neutral_mixed for calculating the contribution of the ion velocity. --- src/anomalous_diffusion.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/anomalous_diffusion.cxx b/src/anomalous_diffusion.cxx index ef2b8d450..602c41bc2 100644 --- a/src/anomalous_diffusion.cxx +++ b/src/anomalous_diffusion.cxx @@ -137,6 +137,8 @@ void AnomalousDiffusion::transform_impl(GuardedOptions& state) { flow_xlow, flow_ylow)); add(species["energy_flow_xlow"], flow_xlow); add(species["energy_flow_ylow"], flow_ylow); + + set(species["anomalous_D"], anomalous_D); } if (include_chi) { From 904e071903b93f610a21edfcbb661a827cb5573e Mon Sep 17 00:00:00 2001 From: malamast Date: Mon, 29 Dec 2025 18:39:29 -0800 Subject: [PATCH 06/29] neutral_mixed: Added the contribution of ion velocity to the neutral perpendicular velocity to make it consistent with AFN models. --- include/neutral_mixed.hxx | 9 ++++ src/neutral_mixed.cxx | 94 ++++++++++++++++++++++++++++++++++++++- 2 files changed, 101 insertions(+), 2 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index ca4f8b7c7..69b89d54a 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -37,6 +37,9 @@ private: Field3D Tn; ///< Neutral temperature Field3D Nnlim, Pnlim, logPnlim; // Limited in regions of low density + Field3D NVn_err; ///< Difference from momentum as input from solver + Field3D NVn_solver; ///< Momentum as calculated in the solver + BoutReal AA; ///< Atomic mass (proton = 1) std::vector collision_names; ///< Collisions used for collisionality @@ -66,6 +69,10 @@ private: bool normalise_sources; ///< Normalise input sources? Field3D kappa_n, eta_n; ///< Neutral conduction and viscosity + Field3D kappa_n_perp, eta_n_perp; ///< Neutral conduction and viscosity + Field3D kappa_n_par, eta_n_par; ///< Neutral conduction and viscosity + + BoutReal neutral_lmax; bool nonorthogonal_operators; ///< Use nonorthogonal operators for radial transport? bool precondition{true}; ///< Enable preconditioner? @@ -77,6 +84,8 @@ private: Field3D Sn, Sp, Snv; ///< Particle, pressure and momentum source Field3D sound_speed; ///< Sound speed for use with Lax flux + bool zero_timederivs; ///< Set the time derivatives to zero? + bool output_ddt; ///< Save time derivatives? bool diagnose; ///< Save additional diagnostics? diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 6f632c8f9..f57ca437a 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -93,6 +93,9 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Limit diffusive fluxes to fraction of thermal speed. <0 means off.") .withDefault(0.2); + neutral_lmax = options["neutral_lmax"].doc("Maximum length scale due to the present of walls.") + .withDefault(0.1) / get(alloptions["units"]["meters"]); // Normalised length + diffusion_limit = options["diffusion_limit"] .doc("Upper limit on diffusion coefficient [m^2/s]. <0 means off") .withDefault(-1.0) @@ -127,6 +130,10 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* inv->setOuterBoundaryFlags(INVERT_DC_GRAD | INVERT_AC_GRAD); } + zero_timederivs = options["zero_timederivs"] + .doc("Set the time derivatives to zero?") + .withDefault(false); + // Optionally output time derivatives output_ddt = options["output_ddt"].doc("Save derivatives to output?").withDefault(false); @@ -230,6 +237,16 @@ void NeutralMixed::transform_impl(GuardedOptions& state) { Vn = NVn / (AA * Nnlim); Vn.applyBoundary("neumann"); + // NVn.applyBoundary(); + NVn_solver = NVn; // Save the momentum as calculated by the solver + NVn = AA * Nn * Vn; // Re-calculate consistent with V and N + // Note: Now NV and NV_solver will differ when N < density_floor + NVn_err = NVn - NVn_solver; // This is used in the finally() function + NVn.applyBoundary(); + + Pnlim = softFloor(Pn, pressure_floor); + Pnlim.applyBoundary(); + ///////////////////////////////////////////////////// // Parallel boundary conditions TRACE("Neutral boundary conditions"); @@ -409,10 +426,12 @@ void NeutralMixed::finally(const Options& state) { } } if (flux_limit > 0.0) { + // Thermal velocity of neutrals + Field3D Vnth = sqrt(Tnlim / AA); + // Apply flux limit to diffusion, // using the local thermal speed and pressure gradient magnitude - Field3D Dmax = - flux_limit * sqrt(Tnlim / AA) / (abs(Grad(logPnlim)) + 1. / neutral_lmax); + Field3D Dmax = flux_limit * Vnth / (abs(Grad(logPnlim)) + 1. / neutral_lmax); BOUT_FOR(i, Dnn.getRegion("RGN_NOBNDRY")) { Dnn[i] = Dnn[i] * Dmax[i] / (Dnn[i] + Dmax[i]); } @@ -616,6 +635,70 @@ void NeutralMixed::finally(const Options& state) { Snv = 0; } + // Add the contribution of ion perp velocity (i.e. anomalous transport) + // See eq 20 and 21 by Horsten et al., (2017) + const Options& allspecies = state["species"]; + + for (auto& kv : allspecies.getChildren()) { + // NOTE:: This is only true for d+ ions. How do we generalize? + // How do we include the perpendicular ion velocity from other drifts? + + const Options& species = kv.second; + + if ((kv.first == "e") or !species.isSet("charge") + or (fabs(get(species["charge"])) < 1e-5)) { + continue; // Skip electrons and non-charged ions + } + + // sources/sinks due to anomalous transport + if (species.isSet("anomalous_D")) { + const Field2D anomalous_D = get(species["anomalous_D"]); + + const Field3D Ni = get(species["density"]); + Field2D Ni2D = DC(Ni); + + // Apply Neumann Y boundary condition, so no additional flux into boundary + // Note: Not setting radial (X) boundaries since those set radial fluxes + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + Ni2D(r.ind, mesh->ystart - 1) = Ni2D(r.ind, mesh->ystart); + } + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + Ni2D(r.ind, mesh->yend + 1) = Ni2D(r.ind, mesh->yend); + } + + ddt(Nn) += Div_a_Grad_perp_upwind (Nn * anomalous_D / softFloor(Ni,density_floor), Ni2D); + //NOTE: Here, we used Nn as is done in UEDGE but it supposted to be the equilibrium value of Nn. + + ddt(Pn) += (5. / 3) * Div_a_Grad_perp_upwind ( Pn * anomalous_D / softFloor(Ni,density_floor), Ni2D); + + if (evolve_momentum) { + ddt(NVn) += Div_a_Grad_perp_upwind (NVn * anomalous_D / softFloor(Ni,density_floor), Ni2D); + } + + } + } + + // If N < density_floor then NV and NV_solver may differ + // -> Add term to force NV_solver towards NV + // Note: This correction is calculated in transform() + ddt(NVn) += NVn_err; + + // Ste time derivatives to zero + if (zero_timederivs) { + + Field3D zero {0.0}; + zero.splitParallelSlices(); + zero.yup() = 0.0; + zero.ydown() = 0.0; + + ddt(Nn) = zero; + ddt(Pn) = zero; + if (evolve_momentum) { + ddt(NVn) = zero; + } + return; + } + // Scale time derivatives if (state.isSet("scale_timederivs")) { Field3D scale_timederivs = get(state["scale_timederivs"]); @@ -659,6 +742,13 @@ void NeutralMixed::finally(const Options& state) { } } + // NOTE: Do we need to do that? + // Restore NV to the value returned by the solver + // so that restart files contain the correct values + // Note: Copy boundary condition so dump file has correct boundary. + NVn_solver.setBoundaryTo(NVn); + NVn = NVn_solver; + #if CHECKLEVEL >= 1 for (auto& i : Nn.getRegion("RGN_NOBNDRY")) { if (!std::isfinite(ddt(Nn)[i])) { From ed2e2105f520c4b01b374c43140db98d5371fc7b Mon Sep 17 00:00:00 2001 From: malamast Date: Tue, 30 Dec 2025 10:00:19 -0800 Subject: [PATCH 07/29] neutral_mixed: I separated the way we calculate the flux limiters in the perpendicular and parallel directions. This is what they do in UEDGE as well. -Added the contibution of ion velocities to the neural diffusion perpendicular velocity. This term comes form momentum source terms due to collisions of neutrals with ions. See eq 20 amd 21 of Horsten 2017. --- src/neutral_mixed.cxx | 109 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 88 insertions(+), 21 deletions(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index f57ca437a..aa0fc1068 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -213,6 +213,9 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* DnnPn.setBoundary(std::string("Dnn") + name); DnnNVn.setBoundary(std::string("Dnn") + name); + kappa_n_perp = 0.0, eta_n_perp = 0.0; + kappa_n_par = 0.0, eta_n_par = 0.0; + substitutePermissions("name", {name}); substitutePermissions( "outputs", {"AA", "density", "pressure", "temperature", "momentum", "velocity"}); @@ -419,22 +422,61 @@ void NeutralMixed::finally(const Options& state) { nu += GET_VALUE(Field3D, localstate["collision_frequencies"][collision_name]); } + // Dnn = Vth^2 / sigma Dnn = (Tnlim / AA) / (nu + Rnn); } else { Dnn = (Tnlim / AA) / Rnn; } } + // Dnn = Vth^2 / sigma + Dnn = (Tnlim / AA) / (nu + Rnn); + } else { + Dnn = (Tnlim / AA) / Rnn; + } + + + // Heat conductivity + // Note: This is kappa_n = (5/2) * Pn / (m * nu) + // where nu is the collision frequency used in Dnn + kappa_n = (5. / 2) * Dnn * Nnlim; + + // Viscosity + // Relationship between heat conduction and viscosity for neutral + // gas Chapman, Cowling "The Mathematical Theory of Non-Uniform + // Gases", CUP 1952 Ferziger, Kaper "Mathematical Theory of + // Transport Processes in Gases", 1972 + // eta_n = (2. / 5) * m_n * kappa_n; + // + eta_n = AA * (2. / 5) * kappa_n; + + if (flux_limit > 0.0) { // Thermal velocity of neutrals Field3D Vnth = sqrt(Tnlim / AA); // Apply flux limit to diffusion, // using the local thermal speed and pressure gradient magnitude - Field3D Dmax = flux_limit * Vnth / (abs(Grad(logPnlim)) + 1. / neutral_lmax); + Field3D Dmax = flux_limit * Vnth / (abs(Grad_perp(logPnlim)) + 1. / neutral_lmax); BOUT_FOR(i, Dnn.getRegion("RGN_NOBNDRY")) { Dnn[i] = Dnn[i] * Dmax[i] / (Dnn[i] + Dmax[i]); } + + Field3D kappa_n_max_perp = flux_limit * (3.0 / 2.0 * Vnth * Nnlim) / (abs(Grad_perp(Tn))/Tnlim + 1. / neutral_lmax); + Field3D kappa_n_max_par = flux_limit * (3.0 / 2.0 * Vnth * Nnlim) / (abs(Grad_par(Tn))/Tnlim + 1. / neutral_lmax); + + BOUT_FOR(i, kappa_n.getRegion("RGN_NOBNDRY")) { + kappa_n_perp[i] = kappa_n[i] * kappa_n_max_perp[i] / (kappa_n[i] + kappa_n_max_perp[i]); + kappa_n_par[i] = kappa_n[i] * kappa_n_max_par[i] / (kappa_n[i] + kappa_n_max_par[i]); + } + + Field3D viscosity_factor_perp = 1.0 / (1.0 + eta_n * abs(Grad_perp(Vn)) / (flux_limit * Pnlim)); + Field3D viscosity_factor_par = 1.0 / (1.0 + eta_n * abs(Grad_par(Vn)) / (flux_limit * Pnlim)); + BOUT_FOR(i, eta_n.getRegion("RGN_NOBNDRY")) { + eta_n_perp[i] = eta_n[i] * viscosity_factor_perp[i]; + eta_n_par[i] = eta_n[i] * viscosity_factor_par[i]; + } + } if (diffusion_limit > 0.0) { @@ -448,6 +490,31 @@ void NeutralMixed::finally(const Options& state) { Dnn.clearParallelSlices(); Dnn.applyBoundary(); + mesh->communicate(kappa_n); + kappa_n.clearParallelSlices(); + kappa_n.applyBoundary("neumann"); + + mesh->communicate(kappa_n_perp); + kappa_n_perp.clearParallelSlices(); + kappa_n_perp.applyBoundary("neumann"); + + mesh->communicate(kappa_n_par); + kappa_n_par.clearParallelSlices(); + kappa_n_par.applyBoundary("neumann"); + + mesh->communicate(eta_n); + eta_n.clearParallelSlices(); + eta_n.applyBoundary("neumann"); + + mesh->communicate(eta_n_perp); + eta_n_perp.clearParallelSlices(); + eta_n_perp.applyBoundary("neumann"); + + mesh->communicate(eta_n_par); + eta_n_par.clearParallelSlices(); + eta_n_par.applyBoundary("neumann"); + + // Neutral diffusion parameters have the same boundary condition as Dnn DnnNn = Dnn * Nnlim; DnnPn = Dnn * Pnlim; @@ -464,6 +531,13 @@ void NeutralMixed::finally(const Options& state) { DnnNn(r.ind, mesh->ystart - 1, jz) = -DnnNn(r.ind, mesh->ystart, jz); DnnPn(r.ind, mesh->ystart - 1, jz) = -DnnPn(r.ind, mesh->ystart, jz); DnnNVn(r.ind, mesh->ystart - 1, jz) = -DnnNVn(r.ind, mesh->ystart, jz); + + // Neumann BC for the transport coef + // NOTE: Do we need that? + auto i = indexAt(kappa_n_par, r.ind, mesh->ystart, jz); + auto im = i.ym(); + kappa_n_par[im] = kappa_n_par[i]; + eta_n_par[im] = eta_n_par[i]; } } } @@ -475,6 +549,13 @@ void NeutralMixed::finally(const Options& state) { DnnNn(r.ind, mesh->yend + 1, jz) = -DnnNn(r.ind, mesh->yend, jz); DnnPn(r.ind, mesh->yend + 1, jz) = -DnnPn(r.ind, mesh->yend, jz); DnnNVn(r.ind, mesh->yend + 1, jz) = -DnnNVn(r.ind, mesh->yend, jz); + + // Neumann BC for the transport coef + // NOTE: Do we need that? + auto i = indexAt(kappa_n_par, r.ind, mesh->yend, jz); + auto ip = i.yp(); + kappa_n_par[ip] = kappa_n_par[i]; + eta_n_par[ip] = eta_n_par[i]; } } } @@ -489,20 +570,6 @@ void NeutralMixed::finally(const Options& state) { } } - // Heat conductivity - // Note: This is kappa_n = (5/2) * Pn / (m * nu) - // where nu is the collision frequency used in Dnn - kappa_n = (5. / 2) * DnnNn; - - // Viscosity - // Relationship between heat conduction and viscosity for neutral - // gas Chapman, Cowling "The Mathematical Theory of Non-Uniform - // Gases", CUP 1952 Ferziger, Kaper "Mathematical Theory of - // Transport Processes in Gases", 1972 - // eta_n = (2. / 5) * m_n * kappa_n; - // - eta_n = AA * (2. / 5) * kappa_n; - ///////////////////////////////////////////////////// // Neutral density TRACE("Neutral density"); @@ -550,7 +617,7 @@ void NeutralMixed::finally(const Options& state) { if (neutral_conduction) { ddt(Pn) += (2. / 3) - * Div_par_K_Grad_par_mod(kappa_n, Tn, // Parallel conduction + * Div_par_K_Grad_par_mod(kappa_n_par, Tn, // Parallel conduction ef_cond_par_ylow, false); // No conduction through target boundary @@ -558,11 +625,11 @@ void NeutralMixed::finally(const Options& state) { if (nonorthogonal_operators) { ddt(Pn) += (2. / 3) - * Div_a_Grad_perp_nonorthog(kappa_n, Tn, ef_cond_perp_xlow, ef_cond_perp_ylow); + * Div_a_Grad_perp_nonorthog(kappa_n_perp, Tn, ef_cond_perp_xlow, ef_cond_perp_ylow); } else { ddt(Pn) += (2. / 3) - * Div_a_Grad_perp_flows(kappa_n, Tn, ef_cond_perp_xlow, ef_cond_perp_ylow); + * Div_a_Grad_perp_flows(kappa_n_perp, Tn, ef_cond_perp_xlow, ef_cond_perp_ylow); } // The factor here is likely 3/2 as this is pure energy flow, but needs checking. @@ -608,17 +675,17 @@ void NeutralMixed::finally(const Options& state) { // eta_n = (2. / 5) * kappa_n; Field3D viscosity_source = Div_par_K_Grad_par_mod( // Parallel viscosity - eta_n, Vn, mf_visc_par_ylow, + eta_n_par, Vn, mf_visc_par_ylow, false) // No viscosity through target boundary ; // Perpendicular viscosity if (nonorthogonal_operators) { viscosity_source += - Div_a_Grad_perp_nonorthog(eta_n, Vn, mf_visc_perp_xlow, mf_visc_perp_ylow); + Div_a_Grad_perp_nonorthog(eta_n_perp, Vn, mf_visc_perp_xlow, mf_visc_perp_ylow); } else { viscosity_source += - Div_a_Grad_perp_flows(eta_n, Vn, mf_visc_perp_xlow, mf_visc_perp_ylow); + Div_a_Grad_perp_flows(eta_n_perp, Vn, mf_visc_perp_xlow, mf_visc_perp_ylow); } ddt(NVn) += viscosity_source; From adfae6a2d0fc0b8dd56aea6f3cf991902c435126 Mon Sep 17 00:00:00 2001 From: malamast Date: Tue, 30 Dec 2025 10:20:25 -0800 Subject: [PATCH 08/29] braginskii_collisions: added a user-defined density_floor --- include/braginskii_collisions.hxx | 1 + src/braginskii_collisions.cxx | 11 ++++++++--- 2 files changed, 9 insertions(+), 3 deletions(-) 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/src/braginskii_collisions.cxx b/src/braginskii_collisions.cxx index d66716e69..a9e4cb03e 100644 --- a/src/braginskii_collisions.cxx +++ b/src/braginskii_collisions.cxx @@ -59,8 +59,13 @@ BraginskiiCollisions::BraginskiiCollisions(const std::string& name, Options& all .doc("User-set arbitrary multiplier on electron-ion collision rate") .withDefault(1.0); - diagnose = - options["diagnose"].doc("Output additional diagnostics?").withDefault(false); + density_floor = options["density_floor"] + .doc("Minimum density floor") + .withDefault(1e-8); + + diagnose = options["diagnose"] + .doc("Output additional diagnostics?") + .withDefault(false); setPermissions(readOnly("species:{electrons}:temperature", Regions::Interior)); setPermissions(readOnly("species:{electrons}:density", Regions::Interior)); @@ -144,7 +149,7 @@ 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) * density1[i] / softFloor(density2[i], density_floor); }); add(species2["collision_frequency"], nu); // Total collision frequency From 8f2081f6de19dd9118dccb90cff7424b89a1f899 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Wed, 25 Feb 2026 14:18:01 +0000 Subject: [PATCH 09/29] Fix rebase --- include/isothermal.hxx | 14 ++------------ include/neutral_mixed.hxx | 1 - include/zero_current.hxx | 18 ++---------------- src/neutral_mixed.cxx | 5 ----- 4 files changed, 4 insertions(+), 34 deletions(-) diff --git a/include/isothermal.hxx b/include/isothermal.hxx index 54d322235..00b2e3d3d 100644 --- a/include/isothermal.hxx +++ b/include/isothermal.hxx @@ -16,6 +16,8 @@ private: BoutReal T; ///< The normalised temperature Field3D P; ///< The normalised pressure + bool initialize_from_mesh; ///< Initilize the Field3D T from 2D profiles stored in the mesh file. + bool diagnose; ///< Output additional diagnostics? /// Inputs @@ -31,18 +33,6 @@ private: /// - pressure (if density is set) /// void transform_impl(GuardedOptions& state) override; - - void outputVars(Options &state) override; -private: - std::string name; // Species name - - // BoutReal T; ///< The normalised temperature - Field3D T; ///< The normalised temperature - Field3D P; ///< The normalised pressure - - bool initialize_from_mesh; ///< Initilize the Field3D T from 2D profiles stored in the mesh file. - - bool diagnose; ///< Output additional diagnostics? }; namespace { diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 69b89d54a..2bdbe9efd 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -72,7 +72,6 @@ private: Field3D kappa_n_perp, eta_n_perp; ///< Neutral conduction and viscosity Field3D kappa_n_par, eta_n_par; ///< Neutral conduction and viscosity - BoutReal neutral_lmax; bool nonorthogonal_operators; ///< Use nonorthogonal operators for radial transport? bool precondition{true}; ///< Enable preconditioner? diff --git a/include/zero_current.hxx b/include/zero_current.hxx index 1d6279e11..260ce5730 100644 --- a/include/zero_current.hxx +++ b/include/zero_current.hxx @@ -30,8 +30,10 @@ struct ZeroCurrent : public Component { private: std::string name; ///< Name of this species BoutReal charge; ///< The charge of this species + BoutReal atomic_mass; // Atomic mass Field3D velocity; ///< Species velocity (for writing to output) + Field3D momentum; ///< Species momentum (for writing to output) /// Required inputs /// - species @@ -49,22 +51,6 @@ private: /// - velocity /// void transform_impl(GuardedOptions& state) override; - - void finally(const Options &state) override { - // Get the velocity with boundary condition applied. - // This is for output only - velocity = get(state["species"][name]["velocity"]); - } - - void outputVars(Options &state) override; -private: - std::string name; ///< Name of this species - BoutReal charge; ///< The charge of this species - BoutReal atomic_mass; // Atomic mass - - Field3D velocity; ///< Species velocity (for writing to output) - Field3D momentum; ///< Species momentum (for writing to output) - }; namespace { diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index aa0fc1068..59bab42c9 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -429,11 +429,6 @@ void NeutralMixed::finally(const Options& state) { Dnn = (Tnlim / AA) / Rnn; } } - // Dnn = Vth^2 / sigma - Dnn = (Tnlim / AA) / (nu + Rnn); - } else { - Dnn = (Tnlim / AA) / Rnn; - } // Heat conductivity From d90cf0d71d72b54b69cc597904c1419ab7488ce0 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Thu, 26 Feb 2026 11:31:09 +0000 Subject: [PATCH 10/29] Add new diagnostics --- include/neutral_mixed.hxx | 3 ++- src/neutral_mixed.cxx | 42 +++++++++++++++++++++++++++++++-------- 2 files changed, 36 insertions(+), 9 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 2bdbe9efd..84d90c751 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -46,7 +46,8 @@ private: std::string diffusion_collisions_mode; ///< Collision selection, either afn or multispecies Field3D nu; ///< Collisionality to use for diffusion - Field3D Dnn; ///< Diffusion coefficient + Field3D nu_pseudo_mfp; ///< Pseudo-collision frequency based on mean free path + Field3D Dnn, Dnn_unlimited, Dmax; ///< Diffusion coefficient Field3D DnnNn, DnnPn, DnnTn, DnnNVn; ///< Used for operators BoutReal flux_limit; ///< Diffusive flux limit BoutReal diffusion_limit; ///< Maximum diffusion coefficient diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 59bab42c9..79a470f94 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -354,11 +354,11 @@ void NeutralMixed::finally(const Options& state) { // // - Field3D Rnn = sqrt(Tnlim / AA) - / neutral_lmax; // Neutral-neutral collisions [normalised frequency] + // Pseudo-collisionality representing domain size based neutral MFP limit + nu_pseudo_mfp = sqrt(Tnlim / AA) / neutral_lmax; if (collisionality_override > 0.0) { // user has set an override for collision frequency - Dnn = (Tn / AA) / collisionality_override; + Dnn_unlimited = (Tn / AA) / collisionality_override; } else { if (localstate.isSet("collision_frequency")) { // Collisionality @@ -424,9 +424,9 @@ void NeutralMixed::finally(const Options& state) { // Dnn = Vth^2 / sigma - Dnn = (Tnlim / AA) / (nu + Rnn); + Dnn_unlimited = (Tnlim / AA) / (nu + nu_pseudo_mfp); } else { - Dnn = (Tnlim / AA) / Rnn; + Dnn_unlimited = (Tnlim / AA) / nu_pseudo_mfp; } } @@ -434,7 +434,7 @@ void NeutralMixed::finally(const Options& state) { // Heat conductivity // Note: This is kappa_n = (5/2) * Pn / (m * nu) // where nu is the collision frequency used in Dnn - kappa_n = (5. / 2) * Dnn * Nnlim; + kappa_n = (5. / 2) * Dnn_unlimited * Nnlim; // Viscosity // Relationship between heat conduction and viscosity for neutral @@ -445,6 +445,8 @@ void NeutralMixed::finally(const Options& state) { // eta_n = AA * (2. / 5) * kappa_n; + Dnn = emptyFrom(Dnn_unlimited); + Dmax = emptyFrom(Dnn_unlimited); if (flux_limit > 0.0) { // Thermal velocity of neutrals @@ -452,9 +454,9 @@ void NeutralMixed::finally(const Options& state) { // Apply flux limit to diffusion, // using the local thermal speed and pressure gradient magnitude - Field3D Dmax = flux_limit * Vnth / (abs(Grad_perp(logPnlim)) + 1. / neutral_lmax); + Dmax = flux_limit * Vnth / (abs(Grad_perp(logPnlim)) + 1. / neutral_lmax); BOUT_FOR(i, Dnn.getRegion("RGN_NOBNDRY")) { - Dnn[i] = Dnn[i] * Dmax[i] / (Dnn[i] + Dmax[i]); + Dnn[i] = Dnn_unlimited[i] * Dmax[i] / (Dnn_unlimited[i] + Dmax[i]); } Field3D kappa_n_max_perp = flux_limit * (3.0 / 2.0 * Vnth * Nnlim) / (abs(Grad_perp(Tn))/Tnlim + 1. / neutral_lmax); @@ -472,6 +474,9 @@ void NeutralMixed::finally(const Options& state) { eta_n_par[i] = eta_n[i] * viscosity_factor_par[i]; } + } else { + Dnn = Dnn_unlimited; + Dmax = -1.0; } if (diffusion_limit > 0.0) { @@ -902,6 +907,27 @@ void NeutralMixed::outputVars(Options& state) { {"standard_name", "diffusion coefficient"}, {"long_name", name + " diffusion coefficient"}, {"source", "neutral_mixed"}}); + set_with_attrs(state[fmt::format("Dnn{}_unlimited", name)], Dnn_unlimited, + {{"time_dimension", "t"}, + {"units", "m^2/s"}, + {"conversion", Cs0 * Cs0 / Omega_ci}, + {"standard_name", "diffusion coefficient"}, + {"long_name", name + " unlimited diffusion coefficient"}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[fmt::format("Dnn{}_max", name)], Dmax, + {{"time_dimension", "t"}, + {"units", "m^2/s"}, + {"conversion", Cs0 * Cs0 / Omega_ci}, + {"standard_name", "diffusion coefficient"}, + {"long_name", name + " maximum diffusion coefficient"}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[fmt::format("K{}_mfp_pseudo_coll", name)], nu_pseudo_mfp, + {{"time_dimension", "t"}, + {"units", "s^-1"}, + {"conversion", Omega_ci}, + {"standard_name", "collision frequency"}, + {"long_name", name + " MFP limit pseudo-collisionality"}, + {"source", "neutral_mixed"}}); set_with_attrs(state[std::string("SN") + name], Sn, {{"time_dimension", "t"}, {"units", "m^-3 s^-1"}, From 54c69c60742f6f5621a6d8234de99aa96971f826 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Thu, 26 Feb 2026 11:56:09 +0000 Subject: [PATCH 11/29] Improve Dnn logic flow --- src/neutral_mixed.cxx | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 79a470f94..013b607c7 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -356,10 +356,7 @@ void NeutralMixed::finally(const Options& state) { // Pseudo-collisionality representing domain size based neutral MFP limit nu_pseudo_mfp = sqrt(Tnlim / AA) / neutral_lmax; - if (collisionality_override > 0.0) { - // user has set an override for collision frequency - Dnn_unlimited = (Tn / AA) / collisionality_override; - } else { + if (localstate.isSet("collision_frequency")) { // Collisionality // Braginskii mode: plasma - self collisions and ei, neutrals - CX, IZ @@ -422,14 +419,20 @@ void NeutralMixed::finally(const Options& state) { nu += GET_VALUE(Field3D, localstate["collision_frequencies"][collision_name]); } + } else { + nu = 0.0; + } - // Dnn = Vth^2 / sigma - Dnn_unlimited = (Tnlim / AA) / (nu + nu_pseudo_mfp); + Field3D nu_total; + if (collisionality_override > 0.0) { + nu_total = collisionality_override; } else { - Dnn_unlimited = (Tnlim / AA) / nu_pseudo_mfp; - } + nu_total = nu + nu_pseudo_mfp; } + // Dnn = Vth^2 / nu_total + Dnn_unlimited = (Tnlim / AA) / nu_total; + // Heat conductivity // Note: This is kappa_n = (5/2) * Pn / (m * nu) From c5a582e6b0a52de874fffe56852acd8727e56e62 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Thu, 26 Feb 2026 11:56:23 +0000 Subject: [PATCH 12/29] Add temporary debug variable --- include/neutral_mixed.hxx | 2 + src/neutral_mixed.cxx | 123 +++++++++++++++++++++----------------- 2 files changed, 69 insertions(+), 56 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 84d90c751..85033a967 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -53,6 +53,8 @@ private: BoutReal diffusion_limit; ///< Maximum diffusion coefficient BoutReal neutral_lmax; + Field3D debug; ///< Debug variable FIXME: remove + bool sheath_ydown, sheath_yup; BoutReal density_floor; ///< Minimum Nn used when dividing NVn by Nn to get Vn. diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 013b607c7..d4270b06c 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -250,6 +250,8 @@ void NeutralMixed::transform_impl(GuardedOptions& state) { Pnlim = softFloor(Pn, pressure_floor); Pnlim.applyBoundary(); + debug = 0.0; + ///////////////////////////////////////////////////// // Parallel boundary conditions TRACE("Neutral boundary conditions"); @@ -357,83 +359,85 @@ void NeutralMixed::finally(const Options& state) { // Pseudo-collisionality representing domain size based neutral MFP limit nu_pseudo_mfp = sqrt(Tnlim / AA) / neutral_lmax; - if (localstate.isSet("collision_frequency")) { - // Collisionality - // Braginskii mode: plasma - self collisions and ei, neutrals - CX, IZ - if (collision_names.empty()) { // Calculate only once - at the beginning - - if (diffusion_collisions_mode == "afn") { - for (const auto& collision : - localstate["collision_frequencies"].getChildren()) { - - std::string collision_name = collision.second.name(); - - if ( // Charge exchange - (collisionSpeciesMatch(collision_name, name, "+", "cx", "partial")) or - // Ionisation - (collisionSpeciesMatch(collision_name, name, "+", "iz", "partial")) or - // Neutral-neutral collisions - (collisionSpeciesMatch(collision_name, name, name, "coll", "exact"))) { - collision_names.push_back(collision_name); - } + if (localstate.isSet("collision_frequency")) { + // Collisionality + // Braginskii mode: plasma - self collisions and ei, neutrals - CX, IZ + if (collision_names.empty()) { // Calculate only once - at the beginning + + if (diffusion_collisions_mode == "afn") { + for (const auto& collision : + localstate["collision_frequencies"].getChildren()) { + + std::string collision_name = collision.second.name(); + + if ( // Charge exchange + (collisionSpeciesMatch(collision_name, name, "+", "cx", "partial")) or + // Ionisation + (collisionSpeciesMatch(collision_name, name, "+", "iz", "partial")) or + // Neutral-neutral collisions + (collisionSpeciesMatch(collision_name, name, name, "coll", "exact"))) { + collision_names.push_back(collision_name); } - // Multispecies mode: all collisions and CX are included - } else if (diffusion_collisions_mode == "multispecies") { - for (const auto& collision : - localstate["collision_frequencies"].getChildren()) { - - std::string collision_name = collision.second.name(); - - if ( // Charge exchange - (collisionSpeciesMatch(collision_name, name, "", "cx", "partial")) or - // Any collision (en, in, ee, ii, nn) - (collisionSpeciesMatch(collision_name, name, "", "coll", "partial"))) { - collision_names.push_back(collision_name); - } + } + // Multispecies mode: all collisions and CX are included + } else if (diffusion_collisions_mode == "multispecies") { + for (const auto& collision : + localstate["collision_frequencies"].getChildren()) { + + std::string collision_name = collision.second.name(); + + if ( // Charge exchange + (collisionSpeciesMatch(collision_name, name, "", "cx", "partial")) or + // Any collision (en, in, ee, ii, nn) + (collisionSpeciesMatch(collision_name, name, "", "coll", "partial"))) { + collision_names.push_back(collision_name); } - - } else { - throw BoutException("\ndiffusion_collisions_mode for {:s} must be either " - "multispecies or braginskii", - name); } - if (collision_names.empty()) { - throw BoutException("\tNo collisions found for {:s} in neutral_mixed for " - "selected collisions mode", - name); - } + } else { + throw BoutException("\ndiffusion_collisions_mode for {:s} must be either " + "multispecies or braginskii", + name); + } - // Write chosen collisions to log file - output_info.write("\t{:s} neutral collisionality mode: '{:s}' using ", name, - diffusion_collisions_mode); - for (const auto& collision : collision_names) { - output_info.write("{:s} ", collision); - } - output_info.write("\n"); + if (collision_names.empty()) { + throw BoutException("\tNo collisions found for {:s} in neutral_mixed for " + "selected collisions mode", + name); } - // Collect the collisionalities based on list of names - nu = 0; - for (const auto& collision_name : collision_names) { - nu += GET_VALUE(Field3D, localstate["collision_frequencies"][collision_name]); + // Write chosen collisions to log file + output_info.write("\t{:s} neutral collisionality mode: '{:s}' using ", name, + diffusion_collisions_mode); + for (const auto& collision : collision_names) { + output_info.write("{:s} ", collision); } + output_info.write("\n"); + } + + // Collect the collisionalities based on list of names + nu = 0; + for (const auto& collision_name : collision_names) { + nu += GET_VALUE(Field3D, localstate["collision_frequencies"][collision_name]); + } } else { nu = 0.0; } - + Field3D nu_total; if (collisionality_override > 0.0) { nu_total = collisionality_override; - } else { + } else { nu_total = nu + nu_pseudo_mfp; } - + // Dnn = Vth^2 / nu_total Dnn_unlimited = (Tnlim / AA) / nu_total; + debug = nu; + // Heat conductivity // Note: This is kappa_n = (5/2) * Pn / (m * nu) // where nu is the collision frequency used in Dnn @@ -896,6 +900,13 @@ void NeutralMixed::outputVars(Options& state) { {"source", "neutral_mixed"}}); } if (diagnose) { + set_with_attrs(state["debug"], debug, + {{"time_dimension", "t"}, + {"units", "-"}, + {"conversion", 1}, + {"standard_name", "debug"}, + {"long_name", "debug"}, + {"source", "neutral_mixed"}}); set_with_attrs(state[std::string("T") + name], Tn, {{"time_dimension", "t"}, {"units", "eV"}, From 62cd8c74cc37bf18df94f740d2415f44b61160c0 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Thu, 26 Feb 2026 12:00:59 +0000 Subject: [PATCH 13/29] Fix segfault/wasted allocation --- include/neutral_mixed.hxx | 1 + src/neutral_mixed.cxx | 16 +++++++++++----- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 85033a967..65354ff51 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -47,6 +47,7 @@ private: diffusion_collisions_mode; ///< Collision selection, either afn or multispecies Field3D nu; ///< Collisionality to use for diffusion Field3D nu_pseudo_mfp; ///< Pseudo-collision frequency based on mean free path + Field3D nu_total; ///< Total collision frequency used for diffusion, including pseudo-collisions Field3D Dnn, Dnn_unlimited, Dmax; ///< Diffusion coefficient Field3D DnnNn, DnnPn, DnnTn, DnnNVn; ///< Used for operators BoutReal flux_limit; ///< Diffusive flux limit diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index d4270b06c..ed9b4672c 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -424,14 +424,13 @@ void NeutralMixed::finally(const Options& state) { } else { nu = 0.0; } - - Field3D nu_total; + + + nu_total = nu + nu_pseudo_mfp; if (collisionality_override > 0.0) { nu_total = collisionality_override; - } else { - nu_total = nu + nu_pseudo_mfp; } - + // Dnn = Vth^2 / nu_total Dnn_unlimited = (Tnlim / AA) / nu_total; @@ -942,6 +941,13 @@ void NeutralMixed::outputVars(Options& state) { {"standard_name", "collision frequency"}, {"long_name", name + " MFP limit pseudo-collisionality"}, {"source", "neutral_mixed"}}); + set_with_attrs(state[fmt::format("K{}_Dnn_coll", name)], nu_total, + {{"time_dimension", "t"}, + {"units", "s^-1"}, + {"conversion", Omega_ci}, + {"standard_name", "collision frequency"}, + {"long_name", name + " total collision frequency used in neutral diffusion coefficient"}, + {"source", "neutral_mixed"}}); set_with_attrs(state[std::string("SN") + name], Sn, {{"time_dimension", "t"}, {"units", "m^-3 s^-1"}, From f4b37f0753b107047df02e623d0fdd9080f6ee61 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Thu, 26 Feb 2026 13:53:22 +0000 Subject: [PATCH 14/29] Add diagnostics to allow reproduction in test --- include/neutral_mixed.hxx | 4 +- src/neutral_mixed.cxx | 113 +++++++++++++++++++++++++++++++------- 2 files changed, 97 insertions(+), 20 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 65354ff51..54f2bbbaa 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -50,6 +50,7 @@ private: Field3D nu_total; ///< Total collision frequency used for diffusion, including pseudo-collisions Field3D Dnn, Dnn_unlimited, Dmax; ///< Diffusion coefficient Field3D DnnNn, DnnPn, DnnTn, DnnNVn; ///< Used for operators + Field3D kappa_n_unlimited, kappa_n_max_par, kappa_n_max_perp; BoutReal flux_limit; ///< Diffusive flux limit BoutReal diffusion_limit; ///< Maximum diffusion coefficient BoutReal neutral_lmax; @@ -72,9 +73,10 @@ private: bool evolve_momentum; ///< Evolve parallel momentum? bool normalise_sources; ///< Normalise input sources? - Field3D kappa_n, eta_n; ///< Neutral conduction and viscosity + Field3D kappa_n, eta_n_unlimited; ///< Neutral conduction and viscosity Field3D kappa_n_perp, eta_n_perp; ///< Neutral conduction and viscosity Field3D kappa_n_par, eta_n_par; ///< Neutral conduction and viscosity + Field3D eta_factor_par, eta_factor_perp; ///< Viscosity reduction factor for flux-limiting bool nonorthogonal_operators; ///< Use nonorthogonal operators for radial transport? diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index ed9b4672c..279f77d82 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -250,7 +250,6 @@ void NeutralMixed::transform_impl(GuardedOptions& state) { Pnlim = softFloor(Pn, pressure_floor); Pnlim.applyBoundary(); - debug = 0.0; ///////////////////////////////////////////////////// // Parallel boundary conditions @@ -440,7 +439,7 @@ void NeutralMixed::finally(const Options& state) { // Heat conductivity // Note: This is kappa_n = (5/2) * Pn / (m * nu) // where nu is the collision frequency used in Dnn - kappa_n = (5. / 2) * Dnn_unlimited * Nnlim; + kappa_n_unlimited = (5. / 2) * Dnn_unlimited * Nnlim; // Viscosity // Relationship between heat conduction and viscosity for neutral @@ -449,7 +448,7 @@ void NeutralMixed::finally(const Options& state) { // Transport Processes in Gases", 1972 // eta_n = (2. / 5) * m_n * kappa_n; // - eta_n = AA * (2. / 5) * kappa_n; + eta_n_unlimited = AA * (2. / 5) * kappa_n_unlimited; Dnn = emptyFrom(Dnn_unlimited); Dmax = emptyFrom(Dnn_unlimited); @@ -465,24 +464,30 @@ void NeutralMixed::finally(const Options& state) { Dnn[i] = Dnn_unlimited[i] * Dmax[i] / (Dnn_unlimited[i] + Dmax[i]); } - Field3D kappa_n_max_perp = flux_limit * (3.0 / 2.0 * Vnth * Nnlim) / (abs(Grad_perp(Tn))/Tnlim + 1. / neutral_lmax); - Field3D kappa_n_max_par = flux_limit * (3.0 / 2.0 * Vnth * Nnlim) / (abs(Grad_par(Tn))/Tnlim + 1. / neutral_lmax); + kappa_n_max_perp = flux_limit * (3.0 / 2.0 * Vnth * Nnlim) / (abs(Grad_perp(Tn))/Tnlim + 1. / neutral_lmax); + kappa_n_max_par = flux_limit * (3.0 / 2.0 * Vnth * Nnlim) / (abs(Grad_par(Tn))/Tnlim + 1. / neutral_lmax); - BOUT_FOR(i, kappa_n.getRegion("RGN_NOBNDRY")) { - kappa_n_perp[i] = kappa_n[i] * kappa_n_max_perp[i] / (kappa_n[i] + kappa_n_max_perp[i]); - kappa_n_par[i] = kappa_n[i] * kappa_n_max_par[i] / (kappa_n[i] + kappa_n_max_par[i]); + debug = abs(Grad_par(Tn)); + + BOUT_FOR(i, kappa_n_unlimited.getRegion("RGN_NOBNDRY")) { + kappa_n_perp[i] = kappa_n_unlimited[i] * kappa_n_max_perp[i] / (kappa_n_unlimited[i] + kappa_n_max_perp[i]); + kappa_n_par[i] = kappa_n_unlimited[i] * kappa_n_max_par[i] / (kappa_n_unlimited[i] + kappa_n_max_par[i]); } - Field3D viscosity_factor_perp = 1.0 / (1.0 + eta_n * abs(Grad_perp(Vn)) / (flux_limit * Pnlim)); - Field3D viscosity_factor_par = 1.0 / (1.0 + eta_n * abs(Grad_par(Vn)) / (flux_limit * Pnlim)); - BOUT_FOR(i, eta_n.getRegion("RGN_NOBNDRY")) { - eta_n_perp[i] = eta_n[i] * viscosity_factor_perp[i]; - eta_n_par[i] = eta_n[i] * viscosity_factor_par[i]; + eta_factor_perp = 1.0 / (1.0 + eta_n_unlimited * abs(Grad_perp(Vn)) / (flux_limit * Pnlim)); + eta_factor_par = 1.0 / (1.0 + eta_n_unlimited * abs(Grad_par(Vn)) / (flux_limit * Pnlim)); + BOUT_FOR(i, eta_n_unlimited.getRegion("RGN_NOBNDRY")) { + eta_n_perp[i] = eta_n_unlimited[i] * eta_factor_perp[i]; + eta_n_par[i] = eta_n_unlimited[i] * eta_factor_par[i]; } } else { Dnn = Dnn_unlimited; Dmax = -1.0; + kappa_n_perp = kappa_n_unlimited; + kappa_n_par = kappa_n_unlimited; + kappa_n_max_perp = -1.0; + kappa_n_max_par = -1.0; } if (diffusion_limit > 0.0) { @@ -496,9 +501,9 @@ void NeutralMixed::finally(const Options& state) { Dnn.clearParallelSlices(); Dnn.applyBoundary(); - mesh->communicate(kappa_n); - kappa_n.clearParallelSlices(); - kappa_n.applyBoundary("neumann"); + mesh->communicate(kappa_n_unlimited); + kappa_n_unlimited.clearParallelSlices(); + kappa_n_unlimited.applyBoundary("neumann"); mesh->communicate(kappa_n_perp); kappa_n_perp.clearParallelSlices(); @@ -508,9 +513,9 @@ void NeutralMixed::finally(const Options& state) { kappa_n_par.clearParallelSlices(); kappa_n_par.applyBoundary("neumann"); - mesh->communicate(eta_n); - eta_n.clearParallelSlices(); - eta_n.applyBoundary("neumann"); + mesh->communicate(eta_n_unlimited); + eta_n_unlimited.clearParallelSlices(); + eta_n_unlimited.applyBoundary("neumann"); mesh->communicate(eta_n_perp); eta_n_perp.clearParallelSlices(); @@ -948,6 +953,76 @@ void NeutralMixed::outputVars(Options& state) { {"standard_name", "collision frequency"}, {"long_name", name + " total collision frequency used in neutral diffusion coefficient"}, {"source", "neutral_mixed"}}); + set_with_attrs(state[fmt::format("kappa_{}_unlimited", name)], kappa_n_unlimited, + {{"time_dimension", "t"}, + {"units", "W / m / eV"}, + {"conversion", (Pnorm * Omega_ci * SQ(rho_s0)) / Tnorm}, + {"standard_name", "conductivity"}, + {"long_name", name + " unlimited conductivity"}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[fmt::format("kappamax_{}_perp", name)], kappa_n_max_perp, + {{"time_dimension", "t"}, + {"units", "W / m / eV"}, + {"conversion", (Pnorm * Omega_ci * SQ(rho_s0)) / Tnorm}, + {"standard_name", "conductivity"}, + {"long_name", name + " maximum perpendicular conductivity"}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[fmt::format("kappamax_{}_par", name)], kappa_n_max_par, + {{"time_dimension", "t"}, + {"units", "W / m / eV"}, + {"conversion", (Pnorm * Omega_ci * SQ(rho_s0)) / Tnorm}, + {"standard_name", "conductivity"}, + {"long_name", name + " maximum parallel conductivity"}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[fmt::format("kappa_{}_perp", name)], kappa_n_perp, + {{"time_dimension", "t"}, + {"units", "W / m / eV"}, + {"conversion", (Pnorm * Omega_ci * SQ(rho_s0)) / Tnorm}, + {"standard_name", "conductivity"}, + {"long_name", name + " perpendicular conductivity"}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[fmt::format("kappa_{}_par", name)], kappa_n_par, + {{"time_dimension", "t"}, + {"units", "W / m / eV"}, + {"conversion", (Pnorm * Omega_ci * SQ(rho_s0)) / Tnorm}, + {"standard_name", "conductivity"}, + {"long_name", name + " parallel conductivity"}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[fmt::format("eta_{}_unlimited", name)], eta_n_unlimited, + {{"time_dimension", "t"}, + {"units", "Pa s"}, + {"conversion", Pnorm / Omega_ci}, + {"standard_name", "viscosity"}, + {"long_name", name + " unlimited viscosity"}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[fmt::format("eta_{}_perp", name)], eta_n_perp, + {{"time_dimension", "t"}, + {"units", "Pa s"}, + {"conversion", Pnorm / Omega_ci}, + {"standard_name", "viscosity"}, + {"long_name", name + " perpendicular viscosity"}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[fmt::format("eta_{}_par", name)], eta_n_par, + {{"time_dimension", "t"}, + {"units", "Pa s"}, + {"conversion", Pnorm / Omega_ci}, + {"standard_name", "viscosity"}, + {"long_name", name + " parallel viscosity"}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[fmt::format("eta_factor_{}_perp", name)], eta_factor_perp, + {{"time_dimension", "t"}, + {"units", "-"}, + {"conversion", 1}, + {"standard_name", "viscosity factor"}, + {"long_name", name + " perpendicular viscosity factor"}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[fmt::format("eta_factor_{}_par", name)], eta_factor_par, + {{"time_dimension", "t"}, + {"units", "-"}, + {"conversion", 1}, + {"standard_name", "viscosity factor"}, + {"long_name", name + " parallel viscosity factor"}, + {"source", "neutral_mixed"}}); set_with_attrs(state[std::string("SN") + name], Sn, {{"time_dimension", "t"}, {"units", "m^-3 s^-1"}, From b183c15f95f2690a29f82c3cdac3fd5dba787d37 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Thu, 26 Feb 2026 14:13:02 +0000 Subject: [PATCH 15/29] Refactor eta for consistency with kappa --- include/neutral_mixed.hxx | 2 +- src/neutral_mixed.cxx | 29 +++++++++++++++-------------- 2 files changed, 16 insertions(+), 15 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 54f2bbbaa..f3bde18c8 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -76,7 +76,7 @@ private: Field3D kappa_n, eta_n_unlimited; ///< Neutral conduction and viscosity Field3D kappa_n_perp, eta_n_perp; ///< Neutral conduction and viscosity Field3D kappa_n_par, eta_n_par; ///< Neutral conduction and viscosity - Field3D eta_factor_par, eta_factor_perp; ///< Viscosity reduction factor for flux-limiting + Field3D eta_n_max_par, eta_n_max_perp; ///< Viscosity reduction factor for flux-limiting bool nonorthogonal_operators; ///< Use nonorthogonal operators for radial transport? diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 279f77d82..64bfc9a62 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -474,11 +474,12 @@ void NeutralMixed::finally(const Options& state) { kappa_n_par[i] = kappa_n_unlimited[i] * kappa_n_max_par[i] / (kappa_n_unlimited[i] + kappa_n_max_par[i]); } - eta_factor_perp = 1.0 / (1.0 + eta_n_unlimited * abs(Grad_perp(Vn)) / (flux_limit * Pnlim)); - eta_factor_par = 1.0 / (1.0 + eta_n_unlimited * abs(Grad_par(Vn)) / (flux_limit * Pnlim)); + eta_n_max_perp = flux_limit * Pnlim / abs(Grad_perp(Vn)); + eta_n_max_par = flux_limit * Pnlim / abs(Grad_par(Vn)); + BOUT_FOR(i, eta_n_unlimited.getRegion("RGN_NOBNDRY")) { - eta_n_perp[i] = eta_n_unlimited[i] * eta_factor_perp[i]; - eta_n_par[i] = eta_n_unlimited[i] * eta_factor_par[i]; + eta_n_perp[i] = eta_n_unlimited[i] * eta_n_max_perp[i] / (eta_n_unlimited[i] + eta_n_max_perp[i]); + eta_n_par[i] = eta_n_unlimited[i] * eta_n_max_par[i] / (eta_n_unlimited[i] + eta_n_max_par[i]); } } else { @@ -1009,19 +1010,19 @@ void NeutralMixed::outputVars(Options& state) { {"standard_name", "viscosity"}, {"long_name", name + " parallel viscosity"}, {"source", "neutral_mixed"}}); - set_with_attrs(state[fmt::format("eta_factor_{}_perp", name)], eta_factor_perp, + set_with_attrs(state[fmt::format("eta_max_{}_perp", name)], eta_n_max_perp, {{"time_dimension", "t"}, - {"units", "-"}, - {"conversion", 1}, - {"standard_name", "viscosity factor"}, - {"long_name", name + " perpendicular viscosity factor"}, + {"units", "Pa s"}, + {"conversion", Pnorm / Omega_ci}, + {"standard_name", "viscosity"}, + {"long_name", name + " maximum perpendicular viscosity"}, {"source", "neutral_mixed"}}); - set_with_attrs(state[fmt::format("eta_factor_{}_par", name)], eta_factor_par, + set_with_attrs(state[fmt::format("eta_max_{}_par", name)], eta_n_max_par, {{"time_dimension", "t"}, - {"units", "-"}, - {"conversion", 1}, - {"standard_name", "viscosity factor"}, - {"long_name", name + " parallel viscosity factor"}, + {"units", "Pa s"}, + {"conversion", Pnorm / Omega_ci}, + {"standard_name", "viscosity"}, + {"long_name", name + " maximum parallel viscosity"}, {"source", "neutral_mixed"}}); set_with_attrs(state[std::string("SN") + name], Sn, {{"time_dimension", "t"}, From 4369ef347208d96055480cc806bdcda5e113be49 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Thu, 26 Feb 2026 14:14:39 +0000 Subject: [PATCH 16/29] Improve diagnostic naming --- src/neutral_mixed.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 64bfc9a62..4831ebc9d 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -933,7 +933,7 @@ void NeutralMixed::outputVars(Options& state) { {"standard_name", "diffusion coefficient"}, {"long_name", name + " unlimited diffusion coefficient"}, {"source", "neutral_mixed"}}); - set_with_attrs(state[fmt::format("Dnn{}_max", name)], Dmax, + set_with_attrs(state[fmt::format("Dnn_max_{}", name)], Dmax, {{"time_dimension", "t"}, {"units", "m^2/s"}, {"conversion", Cs0 * Cs0 / Omega_ci}, @@ -961,14 +961,14 @@ void NeutralMixed::outputVars(Options& state) { {"standard_name", "conductivity"}, {"long_name", name + " unlimited conductivity"}, {"source", "neutral_mixed"}}); - set_with_attrs(state[fmt::format("kappamax_{}_perp", name)], kappa_n_max_perp, + set_with_attrs(state[fmt::format("kappa_max_{}_perp", name)], kappa_n_max_perp, {{"time_dimension", "t"}, {"units", "W / m / eV"}, {"conversion", (Pnorm * Omega_ci * SQ(rho_s0)) / Tnorm}, {"standard_name", "conductivity"}, {"long_name", name + " maximum perpendicular conductivity"}, {"source", "neutral_mixed"}}); - set_with_attrs(state[fmt::format("kappamax_{}_par", name)], kappa_n_max_par, + set_with_attrs(state[fmt::format("kappa_max_{}_par", name)], kappa_n_max_par, {{"time_dimension", "t"}, {"units", "W / m / eV"}, {"conversion", (Pnorm * Omega_ci * SQ(rho_s0)) / Tnorm}, From 4bc4a46a6f96251316d901aac02e9c3fd65c70ca Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Thu, 26 Feb 2026 14:57:55 +0000 Subject: [PATCH 17/29] Formatting --- include/neutral_mixed.hxx | 20 ++++---- src/neutral_mixed.cxx | 98 +++++++++++++++++++++------------------ 2 files changed, 63 insertions(+), 55 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index f3bde18c8..ded96914d 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -37,8 +37,8 @@ private: Field3D Tn; ///< Neutral temperature Field3D Nnlim, Pnlim, logPnlim; // Limited in regions of low density - Field3D NVn_err; ///< Difference from momentum as input from solver - Field3D NVn_solver; ///< Momentum as calculated in the solver + Field3D NVn_err; ///< Difference from momentum as input from solver + Field3D NVn_solver; ///< Momentum as calculated in the solver BoutReal AA; ///< Atomic mass (proton = 1) @@ -47,12 +47,13 @@ private: diffusion_collisions_mode; ///< Collision selection, either afn or multispecies Field3D nu; ///< Collisionality to use for diffusion Field3D nu_pseudo_mfp; ///< Pseudo-collision frequency based on mean free path - Field3D nu_total; ///< Total collision frequency used for diffusion, including pseudo-collisions - Field3D Dnn, Dnn_unlimited, Dmax; ///< Diffusion coefficient + Field3D nu_total; ///< Total collision frequency used for diffusion, including + ///< pseudo-collisions + Field3D Dnn, Dnn_unlimited, Dmax; ///< Diffusion coefficient Field3D DnnNn, DnnPn, DnnTn, DnnNVn; ///< Used for operators Field3D kappa_n_unlimited, kappa_n_max_par, kappa_n_max_perp; - BoutReal flux_limit; ///< Diffusive flux limit - BoutReal diffusion_limit; ///< Maximum diffusion coefficient + BoutReal flux_limit; ///< Diffusive flux limit + BoutReal diffusion_limit; ///< Maximum diffusion coefficient BoutReal neutral_lmax; Field3D debug; ///< Debug variable FIXME: remove @@ -73,12 +74,11 @@ private: bool evolve_momentum; ///< Evolve parallel momentum? bool normalise_sources; ///< Normalise input sources? - Field3D kappa_n, eta_n_unlimited; ///< Neutral conduction and viscosity - Field3D kappa_n_perp, eta_n_perp; ///< Neutral conduction and viscosity - Field3D kappa_n_par, eta_n_par; ///< Neutral conduction and viscosity + Field3D kappa_n, eta_n_unlimited; ///< Neutral conduction and viscosity + Field3D kappa_n_perp, eta_n_perp; ///< Neutral conduction and viscosity + Field3D kappa_n_par, eta_n_par; ///< Neutral conduction and viscosity Field3D eta_n_max_par, eta_n_max_perp; ///< Viscosity reduction factor for flux-limiting - bool nonorthogonal_operators; ///< Use nonorthogonal operators for radial transport? bool precondition{true}; ///< Enable preconditioner? bool precon_laplacexy{false}; ///< Use LaplaceXY? diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 4831ebc9d..3df122b5e 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -93,8 +93,10 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Limit diffusive fluxes to fraction of thermal speed. <0 means off.") .withDefault(0.2); - neutral_lmax = options["neutral_lmax"].doc("Maximum length scale due to the present of walls.") - .withDefault(0.1) / get(alloptions["units"]["meters"]); // Normalised length + neutral_lmax = options["neutral_lmax"] + .doc("Maximum length scale due to the present of walls.") + .withDefault(0.1) + / get(alloptions["units"]["meters"]); // Normalised length diffusion_limit = options["diffusion_limit"] .doc("Upper limit on diffusion coefficient [m^2/s]. <0 means off") @@ -131,8 +133,8 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* } zero_timederivs = options["zero_timederivs"] - .doc("Set the time derivatives to zero?") - .withDefault(false); + .doc("Set the time derivatives to zero?") + .withDefault(false); // Optionally output time derivatives output_ddt = @@ -213,8 +215,8 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* DnnPn.setBoundary(std::string("Dnn") + name); DnnNVn.setBoundary(std::string("Dnn") + name); - kappa_n_perp = 0.0, eta_n_perp = 0.0; - kappa_n_par = 0.0, eta_n_par = 0.0; + kappa_n_perp = 0.0, eta_n_perp = 0.0; + kappa_n_par = 0.0, eta_n_par = 0.0; substitutePermissions("name", {name}); substitutePermissions( @@ -241,7 +243,7 @@ void NeutralMixed::transform_impl(GuardedOptions& state) { Vn.applyBoundary("neumann"); // NVn.applyBoundary(); - NVn_solver = NVn; // Save the momentum as calculated by the solver + NVn_solver = NVn; // Save the momentum as calculated by the solver NVn = AA * Nn * Vn; // Re-calculate consistent with V and N // Note: Now NV and NV_solver will differ when N < density_floor NVn_err = NVn - NVn_solver; // This is used in the finally() function @@ -250,7 +252,6 @@ void NeutralMixed::transform_impl(GuardedOptions& state) { Pnlim = softFloor(Pn, pressure_floor); Pnlim.applyBoundary(); - ///////////////////////////////////////////////////// // Parallel boundary conditions TRACE("Neutral boundary conditions"); @@ -364,8 +365,7 @@ void NeutralMixed::finally(const Options& state) { if (collision_names.empty()) { // Calculate only once - at the beginning if (diffusion_collisions_mode == "afn") { - for (const auto& collision : - localstate["collision_frequencies"].getChildren()) { + for (const auto& collision : localstate["collision_frequencies"].getChildren()) { std::string collision_name = collision.second.name(); @@ -380,8 +380,7 @@ void NeutralMixed::finally(const Options& state) { } // Multispecies mode: all collisions and CX are included } else if (diffusion_collisions_mode == "multispecies") { - for (const auto& collision : - localstate["collision_frequencies"].getChildren()) { + for (const auto& collision : localstate["collision_frequencies"].getChildren()) { std::string collision_name = collision.second.name(); @@ -424,7 +423,6 @@ void NeutralMixed::finally(const Options& state) { nu = 0.0; } - nu_total = nu + nu_pseudo_mfp; if (collisionality_override > 0.0) { nu_total = collisionality_override; @@ -433,10 +431,9 @@ void NeutralMixed::finally(const Options& state) { // Dnn = Vth^2 / nu_total Dnn_unlimited = (Tnlim / AA) / nu_total; - debug = nu; - // Heat conductivity + // Heat conductivity // Note: This is kappa_n = (5/2) * Pn / (m * nu) // where nu is the collision frequency used in Dnn kappa_n_unlimited = (5. / 2) * Dnn_unlimited * Nnlim; @@ -455,8 +452,8 @@ void NeutralMixed::finally(const Options& state) { if (flux_limit > 0.0) { // Thermal velocity of neutrals - Field3D Vnth = sqrt(Tnlim / AA); - + Field3D Vnth = sqrt(Tnlim / AA); + // Apply flux limit to diffusion, // using the local thermal speed and pressure gradient magnitude Dmax = flux_limit * Vnth / (abs(Grad_perp(logPnlim)) + 1. / neutral_lmax); @@ -464,22 +461,28 @@ void NeutralMixed::finally(const Options& state) { Dnn[i] = Dnn_unlimited[i] * Dmax[i] / (Dnn_unlimited[i] + Dmax[i]); } - kappa_n_max_perp = flux_limit * (3.0 / 2.0 * Vnth * Nnlim) / (abs(Grad_perp(Tn))/Tnlim + 1. / neutral_lmax); - kappa_n_max_par = flux_limit * (3.0 / 2.0 * Vnth * Nnlim) / (abs(Grad_par(Tn))/Tnlim + 1. / neutral_lmax); + kappa_n_max_perp = flux_limit * (3.0 / 2.0 * Vnth * Nnlim) + / (abs(Grad_perp(Tn)) / Tnlim + 1. / neutral_lmax); + kappa_n_max_par = flux_limit * (3.0 / 2.0 * Vnth * Nnlim) + / (abs(Grad_par(Tn)) / Tnlim + 1. / neutral_lmax); debug = abs(Grad_par(Tn)); BOUT_FOR(i, kappa_n_unlimited.getRegion("RGN_NOBNDRY")) { - kappa_n_perp[i] = kappa_n_unlimited[i] * kappa_n_max_perp[i] / (kappa_n_unlimited[i] + kappa_n_max_perp[i]); - kappa_n_par[i] = kappa_n_unlimited[i] * kappa_n_max_par[i] / (kappa_n_unlimited[i] + kappa_n_max_par[i]); + kappa_n_perp[i] = kappa_n_unlimited[i] * kappa_n_max_perp[i] + / (kappa_n_unlimited[i] + kappa_n_max_perp[i]); + kappa_n_par[i] = kappa_n_unlimited[i] * kappa_n_max_par[i] + / (kappa_n_unlimited[i] + kappa_n_max_par[i]); } eta_n_max_perp = flux_limit * Pnlim / abs(Grad_perp(Vn)); eta_n_max_par = flux_limit * Pnlim / abs(Grad_par(Vn)); - + BOUT_FOR(i, eta_n_unlimited.getRegion("RGN_NOBNDRY")) { - eta_n_perp[i] = eta_n_unlimited[i] * eta_n_max_perp[i] / (eta_n_unlimited[i] + eta_n_max_perp[i]); - eta_n_par[i] = eta_n_unlimited[i] * eta_n_max_par[i] / (eta_n_unlimited[i] + eta_n_max_par[i]); + eta_n_perp[i] = eta_n_unlimited[i] * eta_n_max_perp[i] + / (eta_n_unlimited[i] + eta_n_max_perp[i]); + eta_n_par[i] = + eta_n_unlimited[i] * eta_n_max_par[i] / (eta_n_unlimited[i] + eta_n_max_par[i]); } } else { @@ -526,7 +529,6 @@ void NeutralMixed::finally(const Options& state) { eta_n_par.clearParallelSlices(); eta_n_par.applyBoundary("neumann"); - // Neutral diffusion parameters have the same boundary condition as Dnn DnnNn = Dnn * Nnlim; DnnPn = Dnn * Pnlim; @@ -543,7 +545,7 @@ void NeutralMixed::finally(const Options& state) { DnnNn(r.ind, mesh->ystart - 1, jz) = -DnnNn(r.ind, mesh->ystart, jz); DnnPn(r.ind, mesh->ystart - 1, jz) = -DnnPn(r.ind, mesh->ystart, jz); DnnNVn(r.ind, mesh->ystart - 1, jz) = -DnnNVn(r.ind, mesh->ystart, jz); - + // Neumann BC for the transport coef // NOTE: Do we need that? auto i = indexAt(kappa_n_par, r.ind, mesh->ystart, jz); @@ -635,9 +637,9 @@ void NeutralMixed::finally(const Options& state) { // Perpendicular conduction if (nonorthogonal_operators) { - ddt(Pn) += - (2. / 3) - * Div_a_Grad_perp_nonorthog(kappa_n_perp, Tn, ef_cond_perp_xlow, ef_cond_perp_ylow); + ddt(Pn) += (2. / 3) + * Div_a_Grad_perp_nonorthog(kappa_n_perp, Tn, ef_cond_perp_xlow, + ef_cond_perp_ylow); } else { ddt(Pn) += (2. / 3) @@ -693,8 +695,8 @@ void NeutralMixed::finally(const Options& state) { // Perpendicular viscosity if (nonorthogonal_operators) { - viscosity_source += - Div_a_Grad_perp_nonorthog(eta_n_perp, Vn, mf_visc_perp_xlow, mf_visc_perp_ylow); + viscosity_source += Div_a_Grad_perp_nonorthog(eta_n_perp, Vn, mf_visc_perp_xlow, + mf_visc_perp_ylow); } else { viscosity_source += Div_a_Grad_perp_flows(eta_n_perp, Vn, mf_visc_perp_xlow, mf_visc_perp_ylow); @@ -719,7 +721,7 @@ void NeutralMixed::finally(const Options& state) { const Options& allspecies = state["species"]; for (auto& kv : allspecies.getChildren()) { - // NOTE:: This is only true for d+ ions. How do we generalize? + // NOTE:: This is only true for d+ ions. How do we generalize? // How do we include the perpendicular ion velocity from other drifts? const Options& species = kv.second; @@ -745,15 +747,19 @@ void NeutralMixed::finally(const Options& state) { Ni2D(r.ind, mesh->yend + 1) = Ni2D(r.ind, mesh->yend); } - ddt(Nn) += Div_a_Grad_perp_upwind (Nn * anomalous_D / softFloor(Ni,density_floor), Ni2D); - //NOTE: Here, we used Nn as is done in UEDGE but it supposted to be the equilibrium value of Nn. + ddt(Nn) += + Div_a_Grad_perp_upwind(Nn * anomalous_D / softFloor(Ni, density_floor), Ni2D); + // NOTE: Here, we used Nn as is done in UEDGE but it supposted to be the equilibrium + // value of Nn. - ddt(Pn) += (5. / 3) * Div_a_Grad_perp_upwind ( Pn * anomalous_D / softFloor(Ni,density_floor), Ni2D); + ddt(Pn) += + (5. / 3) + * Div_a_Grad_perp_upwind(Pn * anomalous_D / softFloor(Ni, density_floor), Ni2D); if (evolve_momentum) { - ddt(NVn) += Div_a_Grad_perp_upwind (NVn * anomalous_D / softFloor(Ni,density_floor), Ni2D); + ddt(NVn) += Div_a_Grad_perp_upwind( + NVn * anomalous_D / softFloor(Ni, density_floor), Ni2D); } - } } @@ -765,7 +771,7 @@ void NeutralMixed::finally(const Options& state) { // Ste time derivatives to zero if (zero_timederivs) { - Field3D zero {0.0}; + Field3D zero{0.0}; zero.splitParallelSlices(); zero.yup() = 0.0; zero.ydown() = 0.0; @@ -947,13 +953,15 @@ void NeutralMixed::outputVars(Options& state) { {"standard_name", "collision frequency"}, {"long_name", name + " MFP limit pseudo-collisionality"}, {"source", "neutral_mixed"}}); - set_with_attrs(state[fmt::format("K{}_Dnn_coll", name)], nu_total, - {{"time_dimension", "t"}, - {"units", "s^-1"}, - {"conversion", Omega_ci}, - {"standard_name", "collision frequency"}, - {"long_name", name + " total collision frequency used in neutral diffusion coefficient"}, - {"source", "neutral_mixed"}}); + set_with_attrs( + state[fmt::format("K{}_Dnn_coll", name)], nu_total, + {{"time_dimension", "t"}, + {"units", "s^-1"}, + {"conversion", Omega_ci}, + {"standard_name", "collision frequency"}, + {"long_name", + name + " total collision frequency used in neutral diffusion coefficient"}, + {"source", "neutral_mixed"}}); set_with_attrs(state[fmt::format("kappa_{}_unlimited", name)], kappa_n_unlimited, {{"time_dimension", "t"}, {"units", "W / m / eV"}, From 79615af6cb0aa162484bdfb8aec0501da2813818 Mon Sep 17 00:00:00 2001 From: mikekryjak <62797494+mikekryjak@users.noreply.github.com> Date: Thu, 26 Feb 2026 14:57:45 +0000 Subject: [PATCH 18/29] [bot] Apply format changes --- include/evolve_density.hxx | 4 ++-- include/evolve_pressure.hxx | 4 ++-- include/fixed_density.hxx | 15 ++++++++------- include/isothermal.hxx | 3 ++- include/zero_current.hxx | 2 +- src/braginskii_collisions.cxx | 11 ++++------- src/evolve_density.cxx | 13 ++++++------- src/evolve_pressure.cxx | 11 ++++++----- src/isothermal.cxx | 13 +++++++------ src/zero_current.cxx | 8 ++++---- 10 files changed, 42 insertions(+), 42 deletions(-) diff --git a/include/evolve_density.hxx b/include/evolve_density.hxx index 3459e0b03..6d8f0df33 100644 --- a/include/evolve_density.hxx +++ b/include/evolve_density.hxx @@ -86,8 +86,8 @@ private: bool diagnose; ///< Output additional diagnostics? Field3D flow_xlow, flow_ylow; ///< Particle flow diagnostics - bool initialize_from_mesh; ///< Initilize the Field3D N from 2D profiles stored in the mesh file. - + bool initialize_from_mesh; ///< Initilize the Field3D N from 2D profiles stored in the + ///< mesh file. /// This sets in the state /// - species diff --git a/include/evolve_pressure.hxx b/include/evolve_pressure.hxx index 4bd1a2ef3..2f06f3f69 100644 --- a/include/evolve_pressure.hxx +++ b/include/evolve_pressure.hxx @@ -100,8 +100,8 @@ private: bool fix_momentum_boundary_flux; ///< Fix momentum flux to boundary condition? Field3D Sp_nvh; ///< Pressure source due to artificial viscosity Field3D E_PdivV, E_VgradP; ///< Diagnostic energy source terms for p*Div(V) and V*Grad(P) - bool initialize_from_mesh; ///< Initilize the Field3D P from 2D profiles stored in the mesh file. - + bool initialize_from_mesh; ///< Initilize the Field3D P from 2D profiles stored in the + ///< mesh file. /// Inputs /// - species diff --git a/include/fixed_density.hxx b/include/fixed_density.hxx index 1ca183ff5..da20e56b5 100644 --- a/include/fixed_density.hxx +++ b/include/fixed_density.hxx @@ -30,15 +30,16 @@ struct FixedDensity : public Component { N = options["density"].as() / Nnorm; // Initilize the Field3D N from 2D profiles stored in the mesh file. - initialize_from_mesh = options["initialize_from_mesh"] - .doc("Initilize field from 2D profiles stored in the mesh file?") - .withDefault(false); + initialize_from_mesh = + options["initialize_from_mesh"] + .doc("Initilize field from 2D profiles stored in the mesh file?") + .withDefault(false); if (initialize_from_mesh) { // Try to read the initial field from the mesh mesh->get(N, std::string("N") + name + "_init"); // Units: [m^-3/s] - N /= Nnorm; // Normalization - } + N /= Nnorm; // Normalization + } substitutePermissions("name", {name}); substitutePermissions("vars", {"AA", "charge", "density"}); @@ -65,8 +66,8 @@ private: Field3D N; ///< Species density (normalised) - bool initialize_from_mesh; ///< Initilize the Field3D N from 2D profiles stored in the mesh file. - + bool initialize_from_mesh; ///< Initilize the Field3D N from 2D profiles stored in the + ///< mesh file. /// Sets in the state the density, mass and charge of the species /// diff --git a/include/isothermal.hxx b/include/isothermal.hxx index 00b2e3d3d..0c1eead80 100644 --- a/include/isothermal.hxx +++ b/include/isothermal.hxx @@ -16,7 +16,8 @@ private: BoutReal T; ///< The normalised temperature Field3D P; ///< The normalised pressure - bool initialize_from_mesh; ///< Initilize the Field3D T from 2D profiles stored in the mesh file. + bool initialize_from_mesh; ///< Initilize the Field3D T from 2D profiles stored in the + ///< mesh file. bool diagnose; ///< Output additional diagnostics? diff --git a/include/zero_current.hxx b/include/zero_current.hxx index 260ce5730..5649c506d 100644 --- a/include/zero_current.hxx +++ b/include/zero_current.hxx @@ -34,7 +34,7 @@ private: Field3D velocity; ///< Species velocity (for writing to output) Field3D momentum; ///< Species momentum (for writing to output) - + /// Required inputs /// - species /// - diff --git a/src/braginskii_collisions.cxx b/src/braginskii_collisions.cxx index a9e4cb03e..d788cacac 100644 --- a/src/braginskii_collisions.cxx +++ b/src/braginskii_collisions.cxx @@ -59,13 +59,10 @@ 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); + density_floor = options["density_floor"].doc("Minimum density floor").withDefault(1e-8); + + diagnose = + options["diagnose"].doc("Output additional diagnostics?").withDefault(false); setPermissions(readOnly("species:{electrons}:temperature", Regions::Interior)); setPermissions(readOnly("species:{electrons}:density", Regions::Interior)); diff --git a/src/evolve_density.cxx b/src/evolve_density.cxx index f22f5f5f6..21216a22a 100644 --- a/src/evolve_density.cxx +++ b/src/evolve_density.cxx @@ -145,17 +145,16 @@ EvolveDensity::EvolveDensity(std::string name, Options& alloptions, Solver* solv substitutePermissions("outputs", outputs); // Initilize the Field3D N from 2D profiles stored in the mesh file. - initialize_from_mesh = n_options["initialize_from_mesh"] - .doc("Initilize field from 2D profiles stored in the mesh file?") - .withDefault(false); + initialize_from_mesh = + n_options["initialize_from_mesh"] + .doc("Initilize field from 2D profiles stored in the mesh file?") + .withDefault(false); if (initialize_from_mesh) { - // Try to read the initial field from the mesh + // Try to read the initial field from the mesh mesh->get(N, std::string("N") + name + "_init"); // Units: [m^-3/s] - N /= Nnorm; // Normalization - + N /= Nnorm; // Normalization } - } void EvolveDensity::transform_impl(GuardedOptions& state) { diff --git a/src/evolve_pressure.cxx b/src/evolve_pressure.cxx index 5b81e705a..604b71c12 100644 --- a/src/evolve_pressure.cxx +++ b/src/evolve_pressure.cxx @@ -165,14 +165,15 @@ EvolvePressure::EvolvePressure(std::string name, Options& alloptions, Solver* so .withDefault(true); // Initilize the Field3D P from 2D profiles stored in the mesh file. - initialize_from_mesh = p_options["initialize_from_mesh"] - .doc("Initilize field from 2D profiles stored in the mesh file?") - .withDefault(false); + initialize_from_mesh = + p_options["initialize_from_mesh"] + .doc("Initilize field from 2D profiles stored in the mesh file?") + .withDefault(false); if (initialize_from_mesh) { - // Try to read the initial field from the mesh + // Try to read the initial field from the mesh mesh->get(P, std::string("P") + name + "_init"); // Units: Pascal [Pa] - P /= Pnorm; // Normalization + P /= Pnorm; // Normalization } if (source_time_dependent) { diff --git a/src/isothermal.cxx b/src/isothermal.cxx index 7889e596e..075b07d4c 100644 --- a/src/isothermal.cxx +++ b/src/isothermal.cxx @@ -1,7 +1,7 @@ +#include "../include/isothermal.hxx" #include #include -#include "../include/isothermal.hxx" using bout::globals::mesh; @@ -20,14 +20,15 @@ Isothermal::Isothermal(std::string name, Options& alloptions, Solver* UNUSED(sol / Tnorm; // Normalise // Initilize the Field3D T from 2D profiles stored in the mesh file. - initialize_from_mesh = options["initialize_from_mesh"] - .doc("Initilize field from 2D profiles stored in the mesh file?") - .withDefault(false); + initialize_from_mesh = + options["initialize_from_mesh"] + .doc("Initilize field from 2D profiles stored in the mesh file?") + .withDefault(false); if (initialize_from_mesh) { - // Try to read the initial field from the mesh + // Try to read the initial field from the mesh mesh->get(T, std::string("T") + name + "_init"); // Units: [eV] - T /= Tnorm; // Normalization + T /= Tnorm; // Normalization } diagnose = options["diagnose"] diff --git a/src/zero_current.cxx b/src/zero_current.cxx index 12242af54..5c281aed7 100644 --- a/src/zero_current.cxx +++ b/src/zero_current.cxx @@ -1,9 +1,9 @@ #include +#include "../include/hermes_utils.hxx" #include "../include/zero_current.hxx" #include -#include "../include/hermes_utils.hxx" ZeroCurrent::ZeroCurrent(std::string name, Options& alloptions, Solver*) : Component({readIfSet("species:{all_species}:charge"), @@ -16,7 +16,8 @@ ZeroCurrent::ZeroCurrent(std::string name, Options& alloptions, Solver*) ASSERT0(charge != 0.0); - atomic_mass = options["AA"].doc("Particle atomic number.").withDefault(0.0); // Atomic mass + atomic_mass = + options["AA"].doc("Particle atomic number.").withDefault(0.0); // Atomic mass // ASSERT0(atomic_mass != 0.0); substitutePermissions("inputs", {"density", "velocity"}); @@ -91,12 +92,11 @@ void ZeroCurrent::outputVars(Options& state) { {"source", "zero_current"}}); set_with_attrs(state[std::string("NV") + name], momentum, - {{"time_dimension", "t"}, + {{"time_dimension", "t"}, {"units", "kg / m^2 / s"}, {"conversion", SI::Mp * Nnorm * Cs0}, {"long_name", name + " parallel momentum"}, {"standard_name", "momentum"}, {"species", name}, {"source", "zero_current"}}); - } From 33088f92cebd6472fd0137cdface107aaa06cc0f Mon Sep 17 00:00:00 2001 From: mikekryjak <62797494+mikekryjak@users.noreply.github.com> Date: Thu, 26 Feb 2026 14:57:45 +0000 Subject: [PATCH 19/29] [bot] Add last format changes commit to ignore file --- .git-blame-ignore-revs | 1 + 1 file changed, 1 insertion(+) diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs index eef5ab924..a59015dd6 100644 --- a/.git-blame-ignore-revs +++ b/.git-blame-ignore-revs @@ -1,3 +1,4 @@ f5de620399e88f139e236e40675facfd0cfd6b13 07616298519738fad0c1bb259e4d021cd6f51d58 37efa12abe0961add8715b797d3926664977d453 +79615af6cb0aa162484bdfb8aec0501da2813818 From 06e73589eac8c20ec192bb94dfa898f90315502b Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 2 Mar 2026 16:01:27 +0000 Subject: [PATCH 20/29] Make double counting lmax a flag, refactor neutral_lmax was already used in Dnn, it shouldn't be also put into Dmax as this limits it twice! The flag enables legacy behaviour by default. --- include/neutral_mixed.hxx | 1 + src/neutral_mixed.cxx | 60 ++++++++++++++++++++++++--------------- 2 files changed, 38 insertions(+), 23 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index ded96914d..fcd004f52 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -73,6 +73,7 @@ private: bool neutral_conduction; ///< Include heat conduction? bool evolve_momentum; ///< Evolve parallel momentum? bool normalise_sources; ///< Normalise input sources? + bool double_count_lmax; ///< Include neutral_lmax in Dmax and kappa_max as well as Dnn? Field3D kappa_n, eta_n_unlimited; ///< Neutral conduction and viscosity Field3D kappa_n_perp, eta_n_perp; ///< Neutral conduction and viscosity diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 3df122b5e..749b3213d 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -136,6 +136,12 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Set the time derivatives to zero?") .withDefault(false); + // FIXME: Temporary options to enable legacy behaviour. Will be removed. + + double_count_lmax = options["double_count_lmax"] + .doc("Include neutral_lmax in Dmax and kappa_max as well as Dnn?") + .withDefault(true); + // Optionally output time derivatives output_ddt = options["output_ddt"].doc("Save derivatives to output?").withDefault(false); @@ -454,37 +460,45 @@ void NeutralMixed::finally(const Options& state) { // Thermal velocity of neutrals Field3D Vnth = sqrt(Tnlim / AA); - // Apply flux limit to diffusion, - // using the local thermal speed and pressure gradient magnitude - Dmax = flux_limit * Vnth / (abs(Grad_perp(logPnlim)) + 1. / neutral_lmax); - BOUT_FOR(i, Dnn.getRegion("RGN_NOBNDRY")) { - Dnn[i] = Dnn_unlimited[i] * Dmax[i] / (Dnn_unlimited[i] + Dmax[i]); - } - - kappa_n_max_perp = flux_limit * (3.0 / 2.0 * Vnth * Nnlim) - / (abs(Grad_perp(Tn)) / Tnlim + 1. / neutral_lmax); - kappa_n_max_par = flux_limit * (3.0 / 2.0 * Vnth * Nnlim) - / (abs(Grad_par(Tn)) / Tnlim + 1. / neutral_lmax); - - debug = abs(Grad_par(Tn)); + + // Calculate maximum diffusivities + if (double_count_lmax) { + Dmax = flux_limit * Vnth / (abs(Grad_perp(logPnlim)) + 1. / neutral_lmax); + kappa_n_max_perp = flux_limit * (3.0 / 2.0 * Vnth * Nnlim) + / (abs(Grad_perp(Tn)) / Tnlim + 1. / neutral_lmax); + kappa_n_max_par = flux_limit * (3.0 / 2.0 * Vnth * Nnlim) + / (abs(Grad_par(Tn)) / Tnlim + 1. / neutral_lmax); - BOUT_FOR(i, kappa_n_unlimited.getRegion("RGN_NOBNDRY")) { - kappa_n_perp[i] = kappa_n_unlimited[i] * kappa_n_max_perp[i] - / (kappa_n_unlimited[i] + kappa_n_max_perp[i]); - kappa_n_par[i] = kappa_n_unlimited[i] * kappa_n_max_par[i] - / (kappa_n_unlimited[i] + kappa_n_max_par[i]); + } else { + Dmax = flux_limit * Vnth / (abs(Grad_perp(logPnlim))); + kappa_n_max_perp = + flux_limit * (3.0 / 2.0 * Vnth * Nnlim) / (abs(Grad_perp(Tn)) / Tnlim); + kappa_n_max_par = + flux_limit * (3.0 / 2.0 * Vnth * Nnlim) / (abs(Grad_par(Tn)) / Tnlim); } + + // Eta_max never had neutral_lmax added eta_n_max_perp = flux_limit * Pnlim / abs(Grad_perp(Vn)); eta_n_max_par = flux_limit * Pnlim / abs(Grad_par(Vn)); - BOUT_FOR(i, eta_n_unlimited.getRegion("RGN_NOBNDRY")) { - eta_n_perp[i] = eta_n_unlimited[i] * eta_n_max_perp[i] - / (eta_n_unlimited[i] + eta_n_max_perp[i]); - eta_n_par[i] = - eta_n_unlimited[i] * eta_n_max_par[i] / (eta_n_unlimited[i] + eta_n_max_par[i]); + + // Apply limits + static auto apply_limiter = [](BoutReal unlimited, BoutReal max) -> BoutReal { + return unlimited * max / (unlimited + max); + }; + + BOUT_FOR(i, Dnn.getRegion("RGN_NOBNDRY")) { + Dnn[i] = apply_limiter(Dnn_unlimited[i], Dmax[i]); + kappa_n_perp[i] = apply_limiter(kappa_n_unlimited[i], kappa_n_max_perp[i]); + kappa_n_par[i] = apply_limiter(kappa_n_unlimited[i], kappa_n_max_par[i]); + eta_n_perp[i] = apply_limiter(eta_n_unlimited[i], eta_n_max_perp[i]); + eta_n_par[i] = apply_limiter(eta_n_unlimited[i], eta_n_max_par[i]); } + debug = abs(Grad_par(Tn)); + + } else { Dnn = Dnn_unlimited; Dmax = -1.0; From 846f3341246c99c6715714693bde49bacb13a74c Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 2 Mar 2026 16:16:04 +0000 Subject: [PATCH 21/29] Improve formulation of thermal speed Left a flag for legacy behaviour. --- include/neutral_mixed.hxx | 6 +++++- src/neutral_mixed.cxx | 38 +++++++++++++++++++++++++++----------- 2 files changed, 32 insertions(+), 12 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index fcd004f52..fddaeb777 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -56,7 +56,7 @@ private: BoutReal diffusion_limit; ///< Maximum diffusion coefficient BoutReal neutral_lmax; - Field3D debug; ///< Debug variable FIXME: remove + bool sheath_ydown, sheath_yup; @@ -73,7 +73,11 @@ private: bool neutral_conduction; ///< Include heat conduction? bool evolve_momentum; ///< Evolve parallel momentum? bool normalise_sources; ///< Normalise input sources? + + // Temporary variables + Field3D debug; ///< Debug variable FIXME: remove bool double_count_lmax; ///< Include neutral_lmax in Dmax and kappa_max as well as Dnn? + bool legacy_thermal_speed; ///< Use legacy definition of thermal speed in flux limiter? Field3D kappa_n, eta_n_unlimited; ///< Neutral conduction and viscosity Field3D kappa_n_perp, eta_n_perp; ///< Neutral conduction and viscosity diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 749b3213d..8ff58640d 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -142,6 +142,10 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Include neutral_lmax in Dmax and kappa_max as well as Dnn?") .withDefault(true); + legacy_thermal_speed = options["legacy_thermal_speed"] + .doc("Use legacy definition of thermal speed in flux limiter?") + .withDefault(true); + // Optionally output time derivatives output_ddt = options["output_ddt"].doc("Save derivatives to output?").withDefault(false); @@ -458,31 +462,43 @@ void NeutralMixed::finally(const Options& state) { if (flux_limit > 0.0) { // Thermal velocity of neutrals - Field3D Vnth = sqrt(Tnlim / AA); + // Old formulation: 3D thermal speed for advection, that times 3/2 for heat flux + // New formulation: 1D particle flow in 3D maxwellian, 1D heat flow in 3D Maxwellian + // See Stangeby + + Field3D Vnth_pf = 0.0; + Field3D Vnth_hf = 0.0; + + if (legacy_thermal_speed) { + Vnth_pf = sqrt(Tnlim / AA); + Vnth_hf = 3.0 / 2.0 * Vnth_pf; + } else { + Vnth_pf = 0.25 * sqrt(8.0 * Tnlim / (PI * AA)); + Vnth_hf = sqrt(2.0 * Tnlim / (PI * AA)); + } - // Calculate maximum diffusivities + // double_count_lmax includes neutral_lmax also in the maximum Dnn and kappa, + // which is legacy behaviour and double counting. + // eta_max never had neutral_lmax added so is omitted if (double_count_lmax) { - Dmax = flux_limit * Vnth / (abs(Grad_perp(logPnlim)) + 1. / neutral_lmax); - kappa_n_max_perp = flux_limit * (3.0 / 2.0 * Vnth * Nnlim) + Dmax = flux_limit * Vnth_pf / (abs(Grad_perp(logPnlim)) + 1. / neutral_lmax); + kappa_n_max_perp = flux_limit * (Vnth_hf * Nnlim) / (abs(Grad_perp(Tn)) / Tnlim + 1. / neutral_lmax); - kappa_n_max_par = flux_limit * (3.0 / 2.0 * Vnth * Nnlim) + kappa_n_max_par = flux_limit * (Vnth_hf * Nnlim) / (abs(Grad_par(Tn)) / Tnlim + 1. / neutral_lmax); } else { - Dmax = flux_limit * Vnth / (abs(Grad_perp(logPnlim))); + Dmax = flux_limit * Vnth_pf / (abs(Grad_perp(logPnlim))); kappa_n_max_perp = - flux_limit * (3.0 / 2.0 * Vnth * Nnlim) / (abs(Grad_perp(Tn)) / Tnlim); + flux_limit * (Vnth_hf * Nnlim) / (abs(Grad_perp(Tn)) / Tnlim); kappa_n_max_par = - flux_limit * (3.0 / 2.0 * Vnth * Nnlim) / (abs(Grad_par(Tn)) / Tnlim); + flux_limit * (Vnth_hf * Nnlim) / (abs(Grad_par(Tn)) / Tnlim); } - - // Eta_max never had neutral_lmax added eta_n_max_perp = flux_limit * Pnlim / abs(Grad_perp(Vn)); eta_n_max_par = flux_limit * Pnlim / abs(Grad_par(Vn)); - // Apply limits static auto apply_limiter = [](BoutReal unlimited, BoutReal max) -> BoutReal { return unlimited * max / (unlimited + max); From d6e6a09b83b2c1f15a89c5287d968feafde22aa2 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 2 Mar 2026 17:28:12 +0000 Subject: [PATCH 22/29] Add catch for boolean flux_limit setting This was boolean before. Now it's the numerical setting for the limiter. --- src/neutral_mixed.cxx | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 8ff58640d..350df5877 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -88,6 +88,16 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* neutral_lmax = 0.1 / meters; // Normalised length + // Guard against accidentally setting flux_limit = true (parses as 1.0) + if (options.isSet("flux_limit")) { + const std::string raw = options["flux_limit"].as(); + if (raw == "true" || raw == "false") { + throw BoutException( + "flux_limit is no longer a boolean setting. Please use a numeric value.", + raw); + } + } + flux_limit = options["flux_limit"] .doc("Limit diffusive fluxes to fraction of thermal speed. <0 means off.") From 90ae2e3167600766b78da4faf0f7aec96430c508 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 2 Mar 2026 17:41:02 +0000 Subject: [PATCH 23/29] Separate neutral limiters With intuitive defaults. --- include/neutral_mixed.hxx | 7 ++++- src/neutral_mixed.cxx | 64 ++++++++++++++++++++++++++++++--------- 2 files changed, 56 insertions(+), 15 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index fddaeb777..6d7b41465 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -52,7 +52,12 @@ private: Field3D Dnn, Dnn_unlimited, Dmax; ///< Diffusion coefficient Field3D DnnNn, DnnPn, DnnTn, DnnNVn; ///< Used for operators Field3D kappa_n_unlimited, kappa_n_max_par, kappa_n_max_perp; - BoutReal flux_limit; ///< Diffusive flux limit + BoutReal flux_limit_adv; ///< Diffusive flux limit + BoutReal flux_limit_cond_par; ///< Limit for parallel conductive flux + BoutReal flux_limit_cond_perp; ///< Limit for perpendicular conductive flux + BoutReal flux_limit_visc_par; ///< Limit for parallel viscous flux + BoutReal flux_limit_visc_perp; ///< Limit for perpendicular viscous flux + BoutReal diffusion_limit; ///< Maximum diffusion coefficient BoutReal neutral_lmax; diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 350df5877..7300c8482 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -98,11 +98,45 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* } } - flux_limit = + // Guard against accidentally setting flux_limit = true (parses as 1.0) + if (options.isSet("flux_limit")) { + const std::string raw = options["flux_limit"].as(); + if (raw == "true" || raw == "false") { + throw BoutException( + "flux_limit is no longer a boolean setting. Please use a numeric value.", + raw); + } + } + + flux_limit_adv = options["flux_limit"] - .doc("Limit diffusive fluxes to fraction of thermal speed. <0 means off.") + .doc("Limit advective fluxes to fraction of thermal speed. <=0 means no limits applied.") .withDefault(0.2); + flux_limit_cond_perp = + options["flux_limit_cond_perp"] + .doc("Limit perpendicular conductive fluxes to fraction of free-streaming. " + "Defaults to same as advection limiter.") + .withDefault(flux_limit_adv); + + flux_limit_cond_par = + options["flux_limit_cond_par"] + .doc("Limit parallel conductive fluxes to fraction of free-streaming. Defaults " + "to same as perpendicular conduction limiter value.") + .withDefault(flux_limit_cond_perp); + + flux_limit_visc_perp = + options["flux_limit_visc_perp"] + .doc("Limit perpendicular viscous fluxes to fraction of free-streaming. " + "Defaults to same as advection limiter value.") + .withDefault(flux_limit_adv); + + flux_limit_visc_par = + options["flux_limit_visc_par"] + .doc("Limit parallel viscous fluxes to fraction of free-streaming. Defaults to " + "same as perpendicular viscous limiter value.") + .withDefault(flux_limit_visc_perp); + neutral_lmax = options["neutral_lmax"] .doc("Maximum length scale due to the present of walls.") .withDefault(0.1) @@ -470,7 +504,7 @@ void NeutralMixed::finally(const Options& state) { Dnn = emptyFrom(Dnn_unlimited); Dmax = emptyFrom(Dnn_unlimited); - if (flux_limit > 0.0) { + if (flux_limit_adv > 0.0) { // Thermal velocity of neutrals // Old formulation: 3D thermal speed for advection, that times 3/2 for heat flux // New formulation: 1D particle flow in 3D maxwellian, 1D heat flow in 3D Maxwellian @@ -492,22 +526,22 @@ void NeutralMixed::finally(const Options& state) { // which is legacy behaviour and double counting. // eta_max never had neutral_lmax added so is omitted if (double_count_lmax) { - Dmax = flux_limit * Vnth_pf / (abs(Grad_perp(logPnlim)) + 1. / neutral_lmax); - kappa_n_max_perp = flux_limit * (Vnth_hf * Nnlim) + Dmax = flux_limit_adv * Vnth_pf / (abs(Grad_perp(logPnlim)) + 1. / neutral_lmax); + kappa_n_max_perp = flux_limit_cond_perp * (Vnth_hf * Nnlim) / (abs(Grad_perp(Tn)) / Tnlim + 1. / neutral_lmax); - kappa_n_max_par = flux_limit * (Vnth_hf * Nnlim) + kappa_n_max_par = flux_limit_cond_par * (Vnth_hf * Nnlim) / (abs(Grad_par(Tn)) / Tnlim + 1. / neutral_lmax); } else { - Dmax = flux_limit * Vnth_pf / (abs(Grad_perp(logPnlim))); + Dmax = flux_limit_adv * Vnth_pf / (abs(Grad_perp(logPnlim))); kappa_n_max_perp = - flux_limit * (Vnth_hf * Nnlim) / (abs(Grad_perp(Tn)) / Tnlim); + flux_limit_cond_perp * (Vnth_hf * Nnlim) / (abs(Grad_perp(Tn)) / Tnlim); kappa_n_max_par = - flux_limit * (Vnth_hf * Nnlim) / (abs(Grad_par(Tn)) / Tnlim); + flux_limit_cond_par * (Vnth_hf * Nnlim) / (abs(Grad_par(Tn)) / Tnlim); } - eta_n_max_perp = flux_limit * Pnlim / abs(Grad_perp(Vn)); - eta_n_max_par = flux_limit * Pnlim / abs(Grad_par(Vn)); + eta_n_max_perp = flux_limit_visc_perp * Pnlim / abs(Grad_perp(Vn)); + eta_n_max_par = flux_limit_visc_par * Pnlim / abs(Grad_par(Vn)); // Apply limits static auto apply_limiter = [](BoutReal unlimited, BoutReal max) -> BoutReal { @@ -527,11 +561,13 @@ void NeutralMixed::finally(const Options& state) { } else { Dnn = Dnn_unlimited; - Dmax = -1.0; + Dmax = Dnn_unlimited; kappa_n_perp = kappa_n_unlimited; kappa_n_par = kappa_n_unlimited; - kappa_n_max_perp = -1.0; - kappa_n_max_par = -1.0; + kappa_n_max_perp = kappa_n_unlimited; + kappa_n_max_par = kappa_n_unlimited; + eta_n_max_perp = eta_n_unlimited; + eta_n_max_par = eta_n_unlimited; } if (diffusion_limit > 0.0) { From 26208639b6d4d5b88e833bc19b51a0726a5e3f03 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 2 Mar 2026 17:41:41 +0000 Subject: [PATCH 24/29] Formatting --- include/neutral_mixed.hxx | 16 +++++++-------- src/neutral_mixed.cxx | 41 +++++++++++++++++++-------------------- 2 files changed, 27 insertions(+), 30 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 6d7b41465..16891447c 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -52,17 +52,15 @@ private: Field3D Dnn, Dnn_unlimited, Dmax; ///< Diffusion coefficient Field3D DnnNn, DnnPn, DnnTn, DnnNVn; ///< Used for operators Field3D kappa_n_unlimited, kappa_n_max_par, kappa_n_max_perp; - BoutReal flux_limit_adv; ///< Diffusive flux limit - BoutReal flux_limit_cond_par; ///< Limit for parallel conductive flux - BoutReal flux_limit_cond_perp; ///< Limit for perpendicular conductive flux - BoutReal flux_limit_visc_par; ///< Limit for parallel viscous flux - BoutReal flux_limit_visc_perp; ///< Limit for perpendicular viscous flux + BoutReal flux_limit_adv; ///< Diffusive flux limit + BoutReal flux_limit_cond_par; ///< Limit for parallel conductive flux + BoutReal flux_limit_cond_perp; ///< Limit for perpendicular conductive flux + BoutReal flux_limit_visc_par; ///< Limit for parallel viscous flux + BoutReal flux_limit_visc_perp; ///< Limit for perpendicular viscous flux BoutReal diffusion_limit; ///< Maximum diffusion coefficient BoutReal neutral_lmax; - - bool sheath_ydown, sheath_yup; BoutReal density_floor; ///< Minimum Nn used when dividing NVn by Nn to get Vn. @@ -80,8 +78,8 @@ private: bool normalise_sources; ///< Normalise input sources? // Temporary variables - Field3D debug; ///< Debug variable FIXME: remove - bool double_count_lmax; ///< Include neutral_lmax in Dmax and kappa_max as well as Dnn? + Field3D debug; ///< Debug variable FIXME: remove + bool double_count_lmax; ///< Include neutral_lmax in Dmax and kappa_max as well as Dnn? bool legacy_thermal_speed; ///< Use legacy definition of thermal speed in flux limiter? Field3D kappa_n, eta_n_unlimited; ///< Neutral conduction and viscosity diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 7300c8482..ff66f4c18 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -93,8 +93,7 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* const std::string raw = options["flux_limit"].as(); if (raw == "true" || raw == "false") { throw BoutException( - "flux_limit is no longer a boolean setting. Please use a numeric value.", - raw); + "flux_limit is no longer a boolean setting. Please use a numeric value.", raw); } } @@ -103,15 +102,14 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* const std::string raw = options["flux_limit"].as(); if (raw == "true" || raw == "false") { throw BoutException( - "flux_limit is no longer a boolean setting. Please use a numeric value.", - raw); + "flux_limit is no longer a boolean setting. Please use a numeric value.", raw); } } - flux_limit_adv = - options["flux_limit"] - .doc("Limit advective fluxes to fraction of thermal speed. <=0 means no limits applied.") - .withDefault(0.2); + flux_limit_adv = options["flux_limit"] + .doc("Limit advective fluxes to fraction of thermal speed. <=0 " + "means no limits applied.") + .withDefault(0.2); flux_limit_cond_perp = options["flux_limit_cond_perp"] @@ -182,13 +180,15 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* // FIXME: Temporary options to enable legacy behaviour. Will be removed. - double_count_lmax = options["double_count_lmax"] - .doc("Include neutral_lmax in Dmax and kappa_max as well as Dnn?") - .withDefault(true); + double_count_lmax = + options["double_count_lmax"] + .doc("Include neutral_lmax in Dmax and kappa_max as well as Dnn?") + .withDefault(true); - legacy_thermal_speed = options["legacy_thermal_speed"] - .doc("Use legacy definition of thermal speed in flux limiter?") - .withDefault(true); + legacy_thermal_speed = + options["legacy_thermal_speed"] + .doc("Use legacy definition of thermal speed in flux limiter?") + .withDefault(true); // Optionally output time derivatives output_ddt = @@ -539,7 +539,7 @@ void NeutralMixed::finally(const Options& state) { kappa_n_max_par = flux_limit_cond_par * (Vnth_hf * Nnlim) / (abs(Grad_par(Tn)) / Tnlim); } - + eta_n_max_perp = flux_limit_visc_perp * Pnlim / abs(Grad_perp(Vn)); eta_n_max_par = flux_limit_visc_par * Pnlim / abs(Grad_par(Vn)); @@ -549,16 +549,15 @@ void NeutralMixed::finally(const Options& state) { }; BOUT_FOR(i, Dnn.getRegion("RGN_NOBNDRY")) { - Dnn[i] = apply_limiter(Dnn_unlimited[i], Dmax[i]); - kappa_n_perp[i] = apply_limiter(kappa_n_unlimited[i], kappa_n_max_perp[i]); - kappa_n_par[i] = apply_limiter(kappa_n_unlimited[i], kappa_n_max_par[i]); - eta_n_perp[i] = apply_limiter(eta_n_unlimited[i], eta_n_max_perp[i]); - eta_n_par[i] = apply_limiter(eta_n_unlimited[i], eta_n_max_par[i]); + Dnn[i] = apply_limiter(Dnn_unlimited[i], Dmax[i]); + kappa_n_perp[i] = apply_limiter(kappa_n_unlimited[i], kappa_n_max_perp[i]); + kappa_n_par[i] = apply_limiter(kappa_n_unlimited[i], kappa_n_max_par[i]); + eta_n_perp[i] = apply_limiter(eta_n_unlimited[i], eta_n_max_perp[i]); + eta_n_par[i] = apply_limiter(eta_n_unlimited[i], eta_n_max_par[i]); } debug = abs(Grad_par(Tn)); - } else { Dnn = Dnn_unlimited; Dmax = Dnn_unlimited; From 23c7d7cbedbda80ab0cfe4c54862f45093e897a1 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Tue, 3 Mar 2026 10:28:45 +0000 Subject: [PATCH 25/29] Remove spurious neutral_lmax re-definition --- src/neutral_mixed.cxx | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index ff66f4c18..308989353 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -86,8 +86,6 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* lax_flux = options["lax_flux"].doc("Enable stabilising lax flux?").withDefault(true); - neutral_lmax = 0.1 / meters; // Normalised length - // Guard against accidentally setting flux_limit = true (parses as 1.0) if (options.isSet("flux_limit")) { const std::string raw = options["flux_limit"].as(); From 7c3b10b1f0e497d3531b6e2b6a5dba2126f26eb4 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Tue, 3 Mar 2026 10:42:27 +0000 Subject: [PATCH 26/29] Remove flux limit exception duplicate due to poor merge --- src/neutral_mixed.cxx | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 308989353..0eb8a8b3d 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -95,15 +95,6 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* } } - // Guard against accidentally setting flux_limit = true (parses as 1.0) - if (options.isSet("flux_limit")) { - const std::string raw = options["flux_limit"].as(); - if (raw == "true" || raw == "false") { - throw BoutException( - "flux_limit is no longer a boolean setting. Please use a numeric value.", raw); - } - } - flux_limit_adv = options["flux_limit"] .doc("Limit advective fluxes to fraction of thermal speed. <=0 " "means no limits applied.") From b62906bc8b7a38e74657e47c9604c88cd6370c82 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Tue, 3 Mar 2026 14:33:07 +0000 Subject: [PATCH 27/29] Implement AFN style limiter formulation Unlike the previous form, there is a sharpness setting. This is referred to as gamma in the AFN parlance. --- include/neutral_mixed.hxx | 2 ++ src/neutral_mixed.cxx | 45 +++++++++++++++++++++++++++++++-------- 2 files changed, 38 insertions(+), 9 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 16891447c..b4932d4f3 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -60,6 +60,8 @@ private: BoutReal diffusion_limit; ///< Maximum diffusion coefficient BoutReal neutral_lmax; + BoutReal flux_limiter_sharpness; ///< Sharpness of flux limiter transition + bool legacy_limiter_form; ///< Use legacy form of flux limiter rather than SOLPS-style bool sheath_ydown, sheath_yup; diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 0eb8a8b3d..ccfcd24f4 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -124,8 +124,13 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* "same as perpendicular viscous limiter value.") .withDefault(flux_limit_visc_perp); + flux_limiter_sharpness = options["flux_limiter_sharpness"] + .doc("Sharpness of flux limiter transition. Only used if " + "legacy_limiter_form is false. Default is 2.0") + .withDefault(2.0); + neutral_lmax = options["neutral_lmax"] - .doc("Maximum length scale due to the present of walls.") + .doc("Maximum length scale due to the presence of walls.") .withDefault(0.1) / get(alloptions["units"]["meters"]); // Normalised length @@ -145,7 +150,7 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* collisionality_override = options["collisionality_override"] .doc( - "Paramter for overriding the neutral collision frequency in Dn for testing") + "Parameter for overriding the neutral collision frequency in Dn for testing") .withDefault(-1.0); normalise_sources = options["normalise_sources"] @@ -179,6 +184,12 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Use legacy definition of thermal speed in flux limiter?") .withDefault(true); + legacy_limiter_form = + options["legacy_limiter_form"] + .doc("Use legacy form of flux limiter rather than SOLPS-style with sharpness parameter?") + .withDefault(true); + + // Optionally output time derivatives output_ddt = options["output_ddt"].doc("Save derivatives to output?").withDefault(false); @@ -533,16 +544,32 @@ void NeutralMixed::finally(const Options& state) { eta_n_max_par = flux_limit_visc_par * Pnlim / abs(Grad_par(Vn)); // Apply limits - static auto apply_limiter = [](BoutReal unlimited, BoutReal max) -> BoutReal { - return unlimited * max / (unlimited + max); + static auto apply_limiter = [](BoutReal unlimited, BoutReal max, + BoutReal flux_limiter_sharpness, + bool legacy_limiter_form) -> BoutReal { + if (legacy_limiter_form) { + + // Simple harmonic average + return unlimited * max / (unlimited + max); + } else { + // SOLPS style limiter + return unlimited + * pow(1.0 + pow(unlimited / max, flux_limiter_sharpness), + -1.0 / flux_limiter_sharpness); + } }; BOUT_FOR(i, Dnn.getRegion("RGN_NOBNDRY")) { - Dnn[i] = apply_limiter(Dnn_unlimited[i], Dmax[i]); - kappa_n_perp[i] = apply_limiter(kappa_n_unlimited[i], kappa_n_max_perp[i]); - kappa_n_par[i] = apply_limiter(kappa_n_unlimited[i], kappa_n_max_par[i]); - eta_n_perp[i] = apply_limiter(eta_n_unlimited[i], eta_n_max_perp[i]); - eta_n_par[i] = apply_limiter(eta_n_unlimited[i], eta_n_max_par[i]); + Dnn[i] = apply_limiter(Dnn_unlimited[i], Dmax[i], flux_limiter_sharpness, + legacy_limiter_form); + kappa_n_perp[i] = apply_limiter(kappa_n_unlimited[i], kappa_n_max_perp[i], + flux_limiter_sharpness, legacy_limiter_form); + kappa_n_par[i] = apply_limiter(kappa_n_unlimited[i], kappa_n_max_par[i], + flux_limiter_sharpness, legacy_limiter_form); + eta_n_perp[i] = apply_limiter(eta_n_unlimited[i], eta_n_max_perp[i], + flux_limiter_sharpness, legacy_limiter_form); + eta_n_par[i] = apply_limiter(eta_n_unlimited[i], eta_n_max_par[i], + flux_limiter_sharpness, legacy_limiter_form); } debug = abs(Grad_par(Tn)); From d75bac3a0a3f2735e13600ed34647cbd65da3563 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Tue, 3 Mar 2026 15:47:06 +0000 Subject: [PATCH 28/29] Add setting for perp ion velocity coupling Previously always on --- include/neutral_mixed.hxx | 5 ++- src/neutral_mixed.cxx | 72 ++++++++++++++++++++++----------------- 2 files changed, 44 insertions(+), 33 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index b4932d4f3..0558fbc91 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -61,7 +61,7 @@ private: BoutReal diffusion_limit; ///< Maximum diffusion coefficient BoutReal neutral_lmax; BoutReal flux_limiter_sharpness; ///< Sharpness of flux limiter transition - bool legacy_limiter_form; ///< Use legacy form of flux limiter rather than SOLPS-style + bool sheath_ydown, sheath_yup; @@ -78,11 +78,14 @@ private: bool neutral_conduction; ///< Include heat conduction? bool evolve_momentum; ///< Evolve parallel momentum? bool normalise_sources; ///< Normalise input sources? + bool perp_ion_coupling; ///< Include coupling to ion perpendicular velocity? // Temporary variables Field3D debug; ///< Debug variable FIXME: remove bool double_count_lmax; ///< Include neutral_lmax in Dmax and kappa_max as well as Dnn? bool legacy_thermal_speed; ///< Use legacy definition of thermal speed in flux limiter? + bool legacy_limiter_form; ///< Use legacy form of flux limiter rather than SOLPS-style + Field3D kappa_n, eta_n_unlimited; ///< Neutral conduction and viscosity Field3D kappa_n_perp, eta_n_perp; ///< Neutral conduction and viscosity diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index ccfcd24f4..4b4353b69 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -147,6 +147,10 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Include neutral gas heat conduction?") .withDefault(true); + perp_ion_coupling = options["perp_ion_coupling"] + .doc("Include coupling to ion perpendicular velocity?") + .withDefault(true); + collisionality_override = options["collisionality_override"] .doc( @@ -809,47 +813,51 @@ void NeutralMixed::finally(const Options& state) { // Add the contribution of ion perp velocity (i.e. anomalous transport) // See eq 20 and 21 by Horsten et al., (2017) - const Options& allspecies = state["species"]; + if (perp_ion_coupling) { - for (auto& kv : allspecies.getChildren()) { - // NOTE:: This is only true for d+ ions. How do we generalize? - // How do we include the perpendicular ion velocity from other drifts? + + const Options& allspecies = state["species"]; - const Options& species = kv.second; + for (auto& kv : allspecies.getChildren()) { + // NOTE:: This is only true for d+ ions. How do we generalize? + // How do we include the perpendicular ion velocity from other drifts? - if ((kv.first == "e") or !species.isSet("charge") - or (fabs(get(species["charge"])) < 1e-5)) { - continue; // Skip electrons and non-charged ions - } + const Options& species = kv.second; - // sources/sinks due to anomalous transport - if (species.isSet("anomalous_D")) { - const Field2D anomalous_D = get(species["anomalous_D"]); + if ((kv.first == "e") or !species.isSet("charge") + or (fabs(get(species["charge"])) < 1e-5)) { + continue; // Skip electrons and non-charged ions + } - const Field3D Ni = get(species["density"]); - Field2D Ni2D = DC(Ni); + // sources/sinks due to anomalous transport + if (species.isSet("anomalous_D")) { + const Field2D anomalous_D = get(species["anomalous_D"]); - // Apply Neumann Y boundary condition, so no additional flux into boundary - // Note: Not setting radial (X) boundaries since those set radial fluxes - for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { - Ni2D(r.ind, mesh->ystart - 1) = Ni2D(r.ind, mesh->ystart); - } - for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { - Ni2D(r.ind, mesh->yend + 1) = Ni2D(r.ind, mesh->yend); - } + const Field3D Ni = get(species["density"]); + Field2D Ni2D = DC(Ni); - ddt(Nn) += - Div_a_Grad_perp_upwind(Nn * anomalous_D / softFloor(Ni, density_floor), Ni2D); - // NOTE: Here, we used Nn as is done in UEDGE but it supposted to be the equilibrium - // value of Nn. + // Apply Neumann Y boundary condition, so no additional flux into boundary + // Note: Not setting radial (X) boundaries since those set radial fluxes + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + Ni2D(r.ind, mesh->ystart - 1) = Ni2D(r.ind, mesh->ystart); + } + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + Ni2D(r.ind, mesh->yend + 1) = Ni2D(r.ind, mesh->yend); + } - ddt(Pn) += - (5. / 3) - * Div_a_Grad_perp_upwind(Pn * anomalous_D / softFloor(Ni, density_floor), Ni2D); + ddt(Nn) += + Div_a_Grad_perp_upwind(Nn * anomalous_D / softFloor(Ni, density_floor), Ni2D); + // NOTE: Here, we used Nn as is done in UEDGE but it supposted to be the equilibrium + // value of Nn. - if (evolve_momentum) { - ddt(NVn) += Div_a_Grad_perp_upwind( - NVn * anomalous_D / softFloor(Ni, density_floor), Ni2D); + ddt(Pn) += + (5. / 3) + * Div_a_Grad_perp_upwind(Pn * anomalous_D / softFloor(Ni, density_floor), Ni2D); + + if (evolve_momentum) { + ddt(NVn) += Div_a_Grad_perp_upwind( + NVn * anomalous_D / softFloor(Ni, density_floor), Ni2D); + } } } } From 42085c85d1661a28a370262f705a51ec88c14012 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Tue, 3 Mar 2026 15:47:33 +0000 Subject: [PATCH 29/29] Formatting --- include/neutral_mixed.hxx | 6 ++---- src/neutral_mixed.cxx | 39 ++++++++++++++++++--------------------- 2 files changed, 20 insertions(+), 25 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 0558fbc91..d78dc050a 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -61,7 +61,6 @@ private: BoutReal diffusion_limit; ///< Maximum diffusion coefficient BoutReal neutral_lmax; BoutReal flux_limiter_sharpness; ///< Sharpness of flux limiter transition - bool sheath_ydown, sheath_yup; @@ -78,14 +77,13 @@ private: bool neutral_conduction; ///< Include heat conduction? bool evolve_momentum; ///< Evolve parallel momentum? bool normalise_sources; ///< Normalise input sources? - bool perp_ion_coupling; ///< Include coupling to ion perpendicular velocity? + bool perp_ion_coupling; ///< Include coupling to ion perpendicular velocity? // Temporary variables Field3D debug; ///< Debug variable FIXME: remove bool double_count_lmax; ///< Include neutral_lmax in Dmax and kappa_max as well as Dnn? bool legacy_thermal_speed; ///< Use legacy definition of thermal speed in flux limiter? - bool legacy_limiter_form; ///< Use legacy form of flux limiter rather than SOLPS-style - + bool legacy_limiter_form; ///< Use legacy form of flux limiter rather than SOLPS-style Field3D kappa_n, eta_n_unlimited; ///< Neutral conduction and viscosity Field3D kappa_n_perp, eta_n_perp; ///< Neutral conduction and viscosity diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 4b4353b69..8715e2b88 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -125,9 +125,9 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .withDefault(flux_limit_visc_perp); flux_limiter_sharpness = options["flux_limiter_sharpness"] - .doc("Sharpness of flux limiter transition. Only used if " - "legacy_limiter_form is false. Default is 2.0") - .withDefault(2.0); + .doc("Sharpness of flux limiter transition. Only used if " + "legacy_limiter_form is false. Default is 2.0") + .withDefault(2.0); neutral_lmax = options["neutral_lmax"] .doc("Maximum length scale due to the presence of walls.") @@ -148,14 +148,13 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .withDefault(true); perp_ion_coupling = options["perp_ion_coupling"] - .doc("Include coupling to ion perpendicular velocity?") - .withDefault(true); + .doc("Include coupling to ion perpendicular velocity?") + .withDefault(true); - collisionality_override = - options["collisionality_override"] - .doc( - "Parameter for overriding the neutral collision frequency in Dn for testing") - .withDefault(-1.0); + collisionality_override = options["collisionality_override"] + .doc("Parameter for overriding the neutral collision " + "frequency in Dn for testing") + .withDefault(-1.0); normalise_sources = options["normalise_sources"] .doc("Normalise input sources?") @@ -188,11 +187,10 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Use legacy definition of thermal speed in flux limiter?") .withDefault(true); - legacy_limiter_form = - options["legacy_limiter_form"] - .doc("Use legacy form of flux limiter rather than SOLPS-style with sharpness parameter?") - .withDefault(true); - + legacy_limiter_form = options["legacy_limiter_form"] + .doc("Use legacy form of flux limiter rather than " + "SOLPS-style with sharpness parameter?") + .withDefault(true); // Optionally output time derivatives output_ddt = @@ -815,7 +813,6 @@ void NeutralMixed::finally(const Options& state) { // See eq 20 and 21 by Horsten et al., (2017) if (perp_ion_coupling) { - const Options& allspecies = state["species"]; for (auto& kv : allspecies.getChildren()) { @@ -847,12 +844,12 @@ void NeutralMixed::finally(const Options& state) { ddt(Nn) += Div_a_Grad_perp_upwind(Nn * anomalous_D / softFloor(Ni, density_floor), Ni2D); - // NOTE: Here, we used Nn as is done in UEDGE but it supposted to be the equilibrium - // value of Nn. + // NOTE: Here, we used Nn as is done in UEDGE but it supposted to be the + // equilibrium value of Nn. - ddt(Pn) += - (5. / 3) - * Div_a_Grad_perp_upwind(Pn * anomalous_D / softFloor(Ni, density_floor), Ni2D); + ddt(Pn) += (5. / 3) + * Div_a_Grad_perp_upwind( + Pn * anomalous_D / softFloor(Ni, density_floor), Ni2D); if (evolve_momentum) { ddt(NVn) += Div_a_Grad_perp_upwind(