From 7e52f2c2b375e4ec91a8b7a5e4002458270a3f4f Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 8 Aug 2023 21:44:38 -0700 Subject: [PATCH 01/14] neutral_mixed: model improvements Trying to get model closer to Horsten thesis (2019). The net effect of these changes is to: - Increase heat conduction by 5/3 - Increase heat advection by 5/3 - Increase viscosity by 5/2 - Add viscous heating The atomic rates used in diffusion, viscosity, and heat conduction coefficients are still different from Horsten 2019. --- src/neutral_mixed.cxx | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 232ba75a..5034c7e4 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -312,6 +312,11 @@ void NeutralMixed::finally(const Options& state) { } } + // Cross-field diffusivity + // Note: This is kappa_n = (5/2) * Pn / (m * nu) + // where nu is the collision frequency used in Dnn + Field3D kappa_n = (5. / 2) * DnnNn; + // Sound speed appearing in Lax flux for advection terms Field3D sound_speed = sqrt(Tn * (5. / 3) / AA); @@ -338,22 +343,6 @@ void NeutralMixed::finally(const Options& state) { + FV::Div_a_Grad_perp(DnnNVn, logPnlim) // Perpendicular diffusion ; - if (neutral_viscosity) { - // NOTE: The following viscosity terms are are not (yet) balanced - // by a viscous heating term - - // 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) * kappa_n; - // - - ddt(NVn) += AA * FV::Div_a_Grad_perp((2. / 5) * DnnNn, Vn) // Perpendicular viscosity - + AA * FV::Div_par_K_Grad_par((2. / 5) * DnnNn, Vn) // Parallel viscosity - ; - } - if (localstate.isSet("momentum_source")) { Snv = get(localstate["momentum_source"]); ddt(NVn) += Snv; @@ -364,10 +353,10 @@ void NeutralMixed::finally(const Options& state) { TRACE("Neutral pressure"); ddt(Pn) = -FV::Div_par_mod(Pn, Vn, sound_speed) // Advection - - (2. / 3) * Pn * Div_par(Vn) // Compression - + FV::Div_a_Grad_perp(DnnPn, logPnlim) // Perpendicular diffusion - + FV::Div_a_Grad_perp(DnnNn, Tn) // Conduction - + FV::Div_par_K_Grad_par(DnnNn, Tn) // Parallel conduction + - (2. / 3) * Pn * Div_par(Vn) // Compression + + FV::Div_a_Grad_perp((5. / 3) * DnnPn, logPnlim) // Perpendicular advection: q = 5/2 p u_perp + + (2. / 3) * (FV::Div_a_Grad_perp(kappa_n, Tn) // Perpendicular Conduction + + FV::Div_par_K_Grad_par(kappa_n, Tn)) // Parallel conduction ; Sp = pressure_source; @@ -376,6 +365,17 @@ void NeutralMixed::finally(const Options& state) { } ddt(Pn) += Sp; + if (neutral_viscosity) { + Field3D eta_n = AA * (2. / 5) * kappa_n; + + Field3D momentum_source = FV::Div_a_Grad_perp(eta_n, Vn) // Perpendicular viscosity + + FV::Div_par_K_Grad_par(eta_n, Vn) // Parallel viscosity + ; + + ddt(NVn) += momentum_source; // Viscosity + ddt(Pn) -= Vn * momentum_source; // Viscous heating + } + BOUT_FOR(i, Pn.getRegion("RGN_ALL")) { if ((Pn[i] < 1e-9) && (ddt(Pn)[i] < 0.0)) { ddt(Pn)[i] = 0.0; From 0682b562da1c2b0a004416616d90303cf20f9657 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 8 Aug 2023 22:29:42 -0700 Subject: [PATCH 02/14] neutral_mixed: Missed factor of 2/3 Convert from viscous heating power to rate of change of pressure. --- src/neutral_mixed.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 5034c7e4..029a7c71 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -373,7 +373,7 @@ void NeutralMixed::finally(const Options& state) { ; ddt(NVn) += momentum_source; // Viscosity - ddt(Pn) -= Vn * momentum_source; // Viscous heating + ddt(Pn) -= (2. / 3) * Vn * momentum_source; // Viscous heating } BOUT_FOR(i, Pn.getRegion("RGN_ALL")) { From a6ada668b002c73f52dfd215ff066939408ff42a Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 10 Aug 2023 22:51:26 -0700 Subject: [PATCH 03/14] Adding isotropic flux limiters Intended to be similar to Wim Van Uytven et al "Assessment of advanced fluid neutral models for the neutral atoms in the plasma edge and application in ITER geometry" Nucl. Fusion 62 (2022) 086023 Significantly slows convergence; may be something wrong in implementation. --- include/neutral_mixed.hxx | 10 +- src/neutral_mixed.cxx | 192 +++++++++++++++++++++++++++++--------- 2 files changed, 158 insertions(+), 44 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 4e28790b..f4e83514 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -43,13 +43,19 @@ private: Field3D Dnn; ///< Diffusion coefficient Field3D DnnNn, DnnPn, DnnTn, DnnNVn; ///< Used for operators + Field3D eta_n; ///< Viscosity + Field3D kappa_n; ///< Thermal conductivity bool sheath_ydown, sheath_yup; BoutReal nn_floor; ///< Minimum Nn used when dividing NVn by Nn to get Vn. - BoutReal flux_limit; ///< Diffusive flux limit - BoutReal diffusion_limit; ///< Maximum diffusion coefficient + bool flux_limit; ///< Impose flux limiter? + BoutReal flux_limit_alpha; + BoutReal flux_limit_gamma; + Field3D particle_flux_factor; ///< Particle flux scaling factor + Field3D momentum_flux_factor; + Field3D heat_flux_factor; bool neutral_viscosity; ///< include viscosity? diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 029a7c71..13d19599 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -51,13 +51,16 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .withDefault(true); flux_limit = options["flux_limit"] - .doc("Limit diffusive fluxes to fraction of thermal speed. <0 means off.") - .withDefault(0.2); + .doc("Use isotropic flux limiters?") + .withDefault(false); - diffusion_limit = options["diffusion_limit"] - .doc("Upper limit on diffusion coefficient [m^2/s]. <0 means off") - .withDefault(-1.0) - / (meters * meters / seconds); // Normalise + flux_limit_alpha = options["flux_limit_alpha"] + .doc("Scale flux limits") + .withDefault(1.0); + + flux_limit_gamma = options["flux_limit_gamma"] + .doc("Higher values increase sharpness of flux limiting") + .withDefault(2.0); neutral_viscosity = options["neutral_viscosity"] .doc("Include neutral gas viscosity?") @@ -156,6 +159,8 @@ void NeutralMixed::transform(Options& state) { Pnlim = floor(Nnlim * Tn, 1e-8); Pnlim.applyBoundary(); + Tnlim = Pnlim / Nnlim; + ///////////////////////////////////////////////////// // Parallel boundary conditions TRACE("Neutral boundary conditions"); @@ -252,30 +257,13 @@ void NeutralMixed::finally(const Options& state) { BoutReal neutral_lmax = 0.1 / get(state["units"]["meters"]); // Normalised length - Field3D Rnn = sqrt(Tn / AA) / neutral_lmax; // Neutral-neutral collisions [normalised frequency] + Field3D Rnn = sqrt(Tnlim / AA) / neutral_lmax; // Neutral-neutral collisions [normalised frequency] if (localstate.isSet("collision_frequency")) { // Dnn = Vth^2 / sigma - Dnn = (Tn / AA) / (get(localstate["collision_frequency"]) + Rnn); + Dnn = (Tnlim / AA) / (get(localstate["collision_frequency"]) + Rnn); } else { - Dnn = (Tn / AA) / Rnn; - } - - if (flux_limit > 0.0) { - // Apply flux limit to diffusion, - // using the local thermal speed and pressure gradient magnitude - Field3D Dmax = flux_limit * sqrt(Tn / AA) / - (abs(Grad(logPnlim)) + 1. / neutral_lmax); - BOUT_FOR(i, Dmax.getRegion("RGN_NOBNDRY")) { - Dnn[i] = BOUTMIN(Dnn[i], Dmax[i]); - } - } - - if (diffusion_limit > 0.0) { - // Impose an upper limit on the diffusion coefficient - BOUT_FOR(i, Dnn.getRegion("RGN_NOBNDRY")) { - Dnn[i] = BOUTMIN(Dnn[i], diffusion_limit); - } + Dnn = (Tnlim / AA) / Rnn; } mesh->communicate(Dnn); @@ -312,19 +300,100 @@ void NeutralMixed::finally(const Options& state) { } } - // Cross-field diffusivity + // Heat conductivity // Note: This is kappa_n = (5/2) * Pn / (m * nu) // where nu is the collision frequency used in Dnn - Field3D kappa_n = (5. / 2) * DnnNn; + kappa_n = (5. / 2) * DnnNn; + + // Viscosity + eta_n = AA * (2. / 5) * kappa_n; // Sound speed appearing in Lax flux for advection terms Field3D sound_speed = sqrt(Tn * (5. / 3) / AA); + // Set factors that multiply the fluxes + particle_flux_factor = 1.0; + momentum_flux_factor = 1.0; + heat_flux_factor = 1.0; + if (flux_limit) { + // Apply flux limiters + // Note: Fluxes calculated here are cell centre, rather than cell edge + + // Cross-field velocity + Vector3D v_perp = -Dnn * Grad_perp(logPnlim); + + // Total velocity, perpendicular + parallel + Vector3D v_total = v_perp; + ASSERT2(v_total.covariant == true); + auto* coord = mesh->getCoordinates(); + v_total.y = Vn * (coord->J * coord->Bxy); // Set parallel component + + Field3D v_sq = v_perp * v_perp + SQ(Vn); // v dot v + Field3D v_abs = sqrt(v_sq); // |v| + + // Magnitude of the particle flux + Field3D particle_flux_abs = Nnlim * v_abs; + + // Normalised particle flux limit + Field3D particle_limit = Nnlim * 0.25 * sqrt(8 * Tnlim / (PI * AA)); + + // Particle flux reduction factor + particle_flux_factor = pow(1. + pow(particle_flux_abs / (flux_limit_alpha * particle_limit), + flux_limit_gamma), + -1./flux_limit_gamma); + + // BOUT_FOR(i, particle_flux_factor.getRegion("RGN_NOBNDRY")) { + // output.write("particle_flux_factor {}: {} {} {}\n", i, particle_flux_factor[i], particle_flux_abs[i], particle_limit[i]); + // } + + // Flux of parallel momentum + // Question: Should flux-limited particle flux be used here, or original flux? + // Note: momentum flux here doesn't include pressure gradient term + // + // The resulting momentum_flux_factor is applied to the momentum advection and viscosity terms, + // and so also scales the advection of kinetic energy + Vector3D momentum_flux = NVn * v_total; + if (neutral_viscosity) { + momentum_flux -= eta_n * Grad_perp(Vn); + } + Field3D momentum_flux_abs = sqrt(momentum_flux * momentum_flux); + Field3D momentum_limit = Pnlim; + + momentum_flux_factor = pow(1. + pow(momentum_flux_abs / (flux_limit_alpha * momentum_limit), + flux_limit_gamma), + -1./flux_limit_gamma); + + // Flux of heat + // Note: + // - Limiting the heat flux, not energy flux + // - Advection and heat conduction are both vectors, and e.g + // could be in opposite directions. + // - Flux doesn't include compression term, or kinetic energy transport + // that is in the momentum equation + Vector3D heat_flux = (3./2) * Pn * v_total + Pn * v_perp - kappa_n * Grad(Tn); + + Field3D heat_flux_abs = sqrt(heat_flux * heat_flux); + Field3D heat_limit = Pnlim * sqrt(2. * Tnlim / (PI * AA)); + + heat_flux_factor = pow(1. + pow(heat_flux_abs / (flux_limit_alpha * heat_limit), + flux_limit_gamma), + -1./flux_limit_gamma); + + // Communicate guard cells and apply boundary conditions + // because the flux factors will be differentiated + mesh->communicate(particle_flux_factor, momentum_flux_factor, heat_flux_factor); + particle_flux_factor.applyBoundary("neumann"); + momentum_flux_factor.applyBoundary("neumann"); + heat_flux_factor.applyBoundary("neumann"); + } + ///////////////////////////////////////////////////// // Neutral density TRACE("Neutral density"); - ddt(Nn) = -FV::Div_par_mod(Nn, Vn, sound_speed) // Advection - + FV::Div_a_Grad_perp(DnnNn, logPnlim) // Perpendicular diffusion + + // Note: Parallel and perpendicular flux scaled by limiter + ddt(Nn) = -FV::Div_par_mod(Nn, Vn * particle_flux_factor, sound_speed) // Advection + - FV::Div_a_Grad_perp(DnnNn * particle_flux_factor, logPnlim) // Perpendicular diffusion ; Sn = density_source; // Save for possible output @@ -338,9 +407,9 @@ void NeutralMixed::finally(const Options& state) { TRACE("Neutral momentum"); ddt(NVn) = - -AA * FV::Div_par_fvv(Nnlim, Vn, sound_speed) // Momentum flow - - Grad_par(Pn) // Pressure gradient - + FV::Div_a_Grad_perp(DnnNVn, logPnlim) // Perpendicular diffusion + -AA * FV::Div_par_fvv(Nnlim, Vn * momentum_flux_factor, sound_speed) // Momentum flow + - Grad_par(Pn) // Pressure gradient. Not included in flux limit. + + FV::Div_a_Grad_perp(DnnNVn * momentum_flux_factor, logPnlim) // Perpendicular diffusion ; if (localstate.isSet("momentum_source")) { @@ -352,11 +421,11 @@ void NeutralMixed::finally(const Options& state) { // Neutral pressure TRACE("Neutral pressure"); - ddt(Pn) = -FV::Div_par_mod(Pn, Vn, sound_speed) // Advection + ddt(Pn) = -FV::Div_par_mod(Pn, Vn * heat_flux_factor, sound_speed) // Advection - (2. / 3) * Pn * Div_par(Vn) // Compression - + FV::Div_a_Grad_perp((5. / 3) * DnnPn, logPnlim) // Perpendicular advection: q = 5/2 p u_perp - + (2. / 3) * (FV::Div_a_Grad_perp(kappa_n, Tn) // Perpendicular Conduction - + FV::Div_par_K_Grad_par(kappa_n, Tn)) // Parallel conduction + + FV::Div_a_Grad_perp((5. / 3) * DnnPn * heat_flux_factor, logPnlim) // Perpendicular advection: q = 5/2 p u_perp + + (2. / 3) * (FV::Div_a_Grad_perp(kappa_n * heat_flux_factor, Tn) // Perpendicular Conduction + + FV::Div_par_K_Grad_par(kappa_n * heat_flux_factor, Tn)) // Parallel conduction ; Sp = pressure_source; @@ -366,10 +435,8 @@ void NeutralMixed::finally(const Options& state) { ddt(Pn) += Sp; if (neutral_viscosity) { - Field3D eta_n = AA * (2. / 5) * kappa_n; - - Field3D momentum_source = FV::Div_a_Grad_perp(eta_n, Vn) // Perpendicular viscosity - + FV::Div_par_K_Grad_par(eta_n, Vn) // Parallel viscosity + Field3D momentum_source = FV::Div_a_Grad_perp(eta_n * momentum_flux_factor, Vn) // Perpendicular viscosity + + FV::Div_par_K_Grad_par(eta_n * momentum_flux_factor, Vn) // Parallel viscosity ; ddt(NVn) += momentum_source; // Viscosity @@ -414,6 +481,7 @@ void NeutralMixed::outputVars(Options& state) { auto Tnorm = get(state["Tnorm"]); auto Omega_ci = get(state["Omega_ci"]); auto Cs0 = get(state["Cs0"]); + auto rho_s0 = get(state["rho_s0"]); const BoutReal Pnorm = SI::qe * Tnorm * Nnorm; state[std::string("N") + name].setAttributes({{"time_dimension", "t"}, @@ -468,6 +536,7 @@ void NeutralMixed::outputVars(Options& state) { {"standard_name", "temperature"}, {"long_name", name + " temperature"}, {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("Dnn") + name], Dnn, {{"time_dimension", "t"}, {"units", "m^2/s"}, @@ -475,6 +544,21 @@ void NeutralMixed::outputVars(Options& state) { {"standard_name", "diffusion coefficient"}, {"long_name", name + " diffusion coefficient"}, {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("eta_") + name], eta_n, + {{"time_dimension", "t"}, + {"units", "Pa s"}, + {"conversion", SQ(rho_s0) * Omega_ci * SI::Mp * Nnorm}, + {"standard_name", "viscosity"}, + {"long_name", name + " viscosity"}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("kappa_") + name], kappa_n, + {{"time_dimension", "t"}, + {"units", "W / m / eV"}, + {"conversion", SI::qe * Nnorm * rho_s0 * Cs0}, + {"standard_name", "heat conductivity"}, + {"long_name", name + " heat conductivity"}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("SN") + name], Sn, {{"time_dimension", "t"}, {"units", "m^-3 s^-1"}, @@ -513,6 +597,30 @@ void NeutralMixed::outputVars(Options& state) { {"species", name}, {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("particle_flux_factor_") + name], particle_flux_factor, + {{"time_dimension", "t"}, + {"units", ""}, + {"conversion", 1.0}, + {"standard_name", "flux factor"}, + {"long_name", name + " particle flux factor"}, + {"species", name}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("momentum_flux_factor_") + name], momentum_flux_factor, + {{"time_dimension", "t"}, + {"units", ""}, + {"conversion", 1.0}, + {"standard_name", "flux factor"}, + {"long_name", name + " momentum flux factor"}, + {"species", name}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("heat_flux_factor_") + name], heat_flux_factor, + {{"time_dimension", "t"}, + {"units", ""}, + {"conversion", 1.0}, + {"standard_name", "flux factor"}, + {"long_name", name + " heat flux factor"}, + {"species", name}, + {"source", "neutral_mixed"}}); } } @@ -524,7 +632,7 @@ void NeutralMixed::precon(const Options& state, BoutReal gamma) { // Neutral gas diffusion // Solve (1 - gamma*Dnn*Delp2)^{-1} - Field3D coef = -gamma * Dnn; + Field3D coef = -gamma * Dnn * particle_flux_factor; if (state.isSet("scale_timederivs")) { coef *= get(state["scale_timederivs"]); From 769ff8bd94646a56a10488324a155a3da4e9d4e8 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 11 Aug 2023 14:38:13 -0700 Subject: [PATCH 04/14] neutral_mixed: Fix and enable flux limit by default Previous commit accidentally reversed the sign of the cross-field particle flow, resulting in very poor convergence. Flux limits are now enabled by default. --- src/neutral_mixed.cxx | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 13d19599..a674add9 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -52,7 +52,7 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* flux_limit = options["flux_limit"] .doc("Use isotropic flux limiters?") - .withDefault(false); + .withDefault(true); flux_limit_alpha = options["flux_limit_alpha"] .doc("Scale flux limits") @@ -257,13 +257,13 @@ void NeutralMixed::finally(const Options& state) { BoutReal neutral_lmax = 0.1 / get(state["units"]["meters"]); // Normalised length - Field3D Rnn = sqrt(Tnlim / AA) / neutral_lmax; // Neutral-neutral collisions [normalised frequency] + Field3D Rnn = sqrt(Tn / AA) / neutral_lmax; // Neutral-neutral collisions [normalised frequency] if (localstate.isSet("collision_frequency")) { // Dnn = Vth^2 / sigma - Dnn = (Tnlim / AA) / (get(localstate["collision_frequency"]) + Rnn); + Dnn = (Tn / AA) / (get(localstate["collision_frequency"]) + Rnn); } else { - Dnn = (Tnlim / AA) / Rnn; + Dnn = (Tn / AA) / Rnn; } mesh->communicate(Dnn); @@ -393,8 +393,8 @@ void NeutralMixed::finally(const Options& state) { // Note: Parallel and perpendicular flux scaled by limiter ddt(Nn) = -FV::Div_par_mod(Nn, Vn * particle_flux_factor, sound_speed) // Advection - - FV::Div_a_Grad_perp(DnnNn * particle_flux_factor, logPnlim) // Perpendicular diffusion - ; + + FV::Div_a_Grad_perp(DnnNn * particle_flux_factor, logPnlim) // Perpendicular diffusion + ; Sn = density_source; // Save for possible output if (localstate.isSet("density_source")) { From 98c85a024046bf977cc287f7406da491851208e7 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Wed, 16 Aug 2023 13:55:23 +0100 Subject: [PATCH 05/14] Add flags for individual neutral flux limiters --- include/neutral_mixed.hxx | 1 + src/neutral_mixed.cxx | 23 +++++++++++++++++++---- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index f4e83514..d1fed1fc 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -51,6 +51,7 @@ private: BoutReal nn_floor; ///< Minimum Nn used when dividing NVn by Nn to get Vn. bool flux_limit; ///< Impose flux limiter? + bool particle_flux_limiter, heat_flux_limiter, momentum_flux_limiter; ///< Which limiters to impose BoutReal flux_limit_alpha; BoutReal flux_limit_gamma; Field3D particle_flux_factor; ///< Particle flux scaling factor diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index a674add9..d3645bf6 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -54,6 +54,18 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Use isotropic flux limiters?") .withDefault(true); + particle_flux_limiter = options["particle_flux_limiter"] + .doc("Enable particle flux limiter?") + .withDefault(true); + + heat_flux_limiter = options["heat_flux_limiter"] + .doc("Enable heat flux limiter?") + .withDefault(true); + + momentum_flux_limiter = options["momentum_flux_limiter"] + .doc("Enable momentum flux limiter?") + .withDefault(true); + flux_limit_alpha = options["flux_limit_alpha"] .doc("Scale flux limits") .withDefault(1.0); @@ -338,10 +350,13 @@ void NeutralMixed::finally(const Options& state) { Field3D particle_limit = Nnlim * 0.25 * sqrt(8 * Tnlim / (PI * AA)); // Particle flux reduction factor - particle_flux_factor = pow(1. + pow(particle_flux_abs / (flux_limit_alpha * particle_limit), - flux_limit_gamma), - -1./flux_limit_gamma); - + if (particle_flux_limiter){ + particle_flux_factor = pow(1. + pow(particle_flux_abs / (flux_limit_alpha * particle_limit), + flux_limit_gamma), + -1./flux_limit_gamma); + } else { + particle_flux_factor = 1.0; + } // BOUT_FOR(i, particle_flux_factor.getRegion("RGN_NOBNDRY")) { // output.write("particle_flux_factor {}: {} {} {}\n", i, particle_flux_factor[i], particle_flux_abs[i], particle_limit[i]); // } From 321a0a6c7d56091806022741f5c062d05ee13955 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Thu, 17 Aug 2023 19:54:54 +0100 Subject: [PATCH 06/14] Allow suppression of momentum and heat flux factor --- src/neutral_mixed.cxx | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index d3645bf6..254f9c7e 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -374,10 +374,13 @@ void NeutralMixed::finally(const Options& state) { Field3D momentum_flux_abs = sqrt(momentum_flux * momentum_flux); Field3D momentum_limit = Pnlim; - momentum_flux_factor = pow(1. + pow(momentum_flux_abs / (flux_limit_alpha * momentum_limit), - flux_limit_gamma), - -1./flux_limit_gamma); - + if (momentum_flux_limiter) { + momentum_flux_factor = pow(1. + pow(momentum_flux_abs / (flux_limit_alpha * momentum_limit), + flux_limit_gamma), + -1./flux_limit_gamma); + } else { + momentum_flux_factor = 1.0; + } // Flux of heat // Note: // - Limiting the heat flux, not energy flux @@ -390,10 +393,13 @@ void NeutralMixed::finally(const Options& state) { Field3D heat_flux_abs = sqrt(heat_flux * heat_flux); Field3D heat_limit = Pnlim * sqrt(2. * Tnlim / (PI * AA)); - heat_flux_factor = pow(1. + pow(heat_flux_abs / (flux_limit_alpha * heat_limit), - flux_limit_gamma), - -1./flux_limit_gamma); - + if (heat_flux_limiter) { + heat_flux_factor = pow(1. + pow(heat_flux_abs / (flux_limit_alpha * heat_limit), + flux_limit_gamma), + -1./flux_limit_gamma); + } else { + heat_flux_factor = 1.0; + } // Communicate guard cells and apply boundary conditions // because the flux factors will be differentiated mesh->communicate(particle_flux_factor, momentum_flux_factor, heat_flux_factor); From 609b95c9dd30e61d87b0fd5fb5b51d1cb3346ec3 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Tue, 29 Aug 2023 19:23:39 +0100 Subject: [PATCH 07/14] Add option disabling neutral MFP=0.1m - This is for testing and the Rnn term may be removed permanently later - The AFN limiters don't have this. --- include/neutral_mixed.hxx | 1 + src/neutral_mixed.cxx | 10 +++++++++- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index d1fed1fc..372a2b22 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -52,6 +52,7 @@ private: bool flux_limit; ///< Impose flux limiter? bool particle_flux_limiter, heat_flux_limiter, momentum_flux_limiter; ///< Which limiters to impose + bool include_rnn; ///< Include hardcoded max MFP of 0.1m in diffusion calculation? BoutReal flux_limit_alpha; BoutReal flux_limit_gamma; Field3D particle_flux_factor; ///< Particle flux scaling factor diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 254f9c7e..efae05e7 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -78,6 +78,10 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Include neutral gas viscosity?") .withDefault(true); + include_rnn = options["include_rnn"] + .doc("Include hardcoded MFP of 0.1m in diffusion calculation?") + .withDefault(true); + if (precondition) { inv = std::unique_ptr(Laplacian::create(&options["precon_laplace"])); @@ -273,7 +277,11 @@ void NeutralMixed::finally(const Options& state) { if (localstate.isSet("collision_frequency")) { // Dnn = Vth^2 / sigma - Dnn = (Tn / AA) / (get(localstate["collision_frequency"]) + Rnn); + if (include_rnn) { + Dnn = (Tn / AA) / (get(localstate["collision_frequency"]) + Rnn); + } else { + Dnn = (Tn / AA) / (get(localstate["collision_frequency"])); + } } else { Dnn = (Tn / AA) / Rnn; } From ffb536842062f7bbab68feeab2f4c84cd3d28fa8 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Wed, 30 Aug 2023 21:58:42 +0100 Subject: [PATCH 08/14] Make maximum diffusion MFP a user input - For testing --- include/neutral_mixed.hxx | 2 +- src/neutral_mixed.cxx | 23 ++++++++++++++--------- 2 files changed, 15 insertions(+), 10 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 372a2b22..9a1eaa78 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -52,7 +52,7 @@ private: bool flux_limit; ///< Impose flux limiter? bool particle_flux_limiter, heat_flux_limiter, momentum_flux_limiter; ///< Which limiters to impose - bool include_rnn; ///< Include hardcoded max MFP of 0.1m in diffusion calculation? + BoutReal maximum_mfp; ///< Maximum mean free path for diffusion. 0.1 by default, -1 is off. BoutReal flux_limit_alpha; BoutReal flux_limit_gamma; Field3D particle_flux_factor; ///< Particle flux scaling factor diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index efae05e7..cce91999 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -78,9 +78,9 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Include neutral gas viscosity?") .withDefault(true); - include_rnn = options["include_rnn"] - .doc("Include hardcoded MFP of 0.1m in diffusion calculation?") - .withDefault(true); + maximum_mfp = options["maximum_mfp"] + .doc("Optional maximum mean free path in [m] for diffusive processes. -1 is off") + .withDefault(0.1); if (precondition) { inv = std::unique_ptr(Laplacian::create(&options["precon_laplace"])); @@ -277,12 +277,17 @@ void NeutralMixed::finally(const Options& state) { if (localstate.isSet("collision_frequency")) { // Dnn = Vth^2 / sigma - if (include_rnn) { - Dnn = (Tn / AA) / (get(localstate["collision_frequency"]) + Rnn); - } else { - Dnn = (Tn / AA) / (get(localstate["collision_frequency"])); - } - } else { + + if (maximum_mfp != -1) { // MFP limit enabled + Field3D Rnn = sqrt(Tn / AA) / (maximum_mfp / get(state["units"]["meters"])); + Dnn = (Tn / AA) / (get(localstate["collision_frequency"]) + Rnn); + } else { // MFP limit disabled + Dnn = (Tn / AA) / (get(localstate["collision_frequency"])); + } + + } else { // If no collisions, hardcode max MFP to 0.1m + output_warn.write("No collisions set for the neutrals, limiting mean free path to 0.1m"); + Field3D Rnn = sqrt(Tn / AA) / (0.1 / get(state["units"]["meters"])); Dnn = (Tn / AA) / Rnn; } From 2655106522c7f81f0162ca4fdf41873071ab65a8 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 1 Sep 2023 21:28:06 -0700 Subject: [PATCH 09/14] neutral_mixed: Minor tidying of maximum_mfp Check if float option is < 0 rather than equal to -1. Remove some unused and commented-out code. --- src/neutral_mixed.cxx | 14 +++----------- 1 file changed, 3 insertions(+), 11 deletions(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index cce91999..1f546ebf 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -79,7 +79,7 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .withDefault(true); maximum_mfp = options["maximum_mfp"] - .doc("Optional maximum mean free path in [m] for diffusive processes. -1 is off") + .doc("Optional maximum mean free path in [m] for diffusive processes. < 0 is off") .withDefault(0.1); if (precondition) { @@ -270,15 +270,10 @@ void NeutralMixed::finally(const Options& state) { // Calculate cross-field diffusion from collision frequency // // - BoutReal neutral_lmax = - 0.1 / get(state["units"]["meters"]); // Normalised length - - Field3D Rnn = sqrt(Tn / AA) / neutral_lmax; // Neutral-neutral collisions [normalised frequency] - if (localstate.isSet("collision_frequency")) { // Dnn = Vth^2 / sigma - if (maximum_mfp != -1) { // MFP limit enabled + if (maximum_mfp > 0) { // MFP limit enabled Field3D Rnn = sqrt(Tn / AA) / (maximum_mfp / get(state["units"]["meters"])); Dnn = (Tn / AA) / (get(localstate["collision_frequency"]) + Rnn); } else { // MFP limit disabled @@ -370,10 +365,7 @@ void NeutralMixed::finally(const Options& state) { } else { particle_flux_factor = 1.0; } - // BOUT_FOR(i, particle_flux_factor.getRegion("RGN_NOBNDRY")) { - // output.write("particle_flux_factor {}: {} {} {}\n", i, particle_flux_factor[i], particle_flux_abs[i], particle_limit[i]); - // } - + // Flux of parallel momentum // Question: Should flux-limited particle flux be used here, or original flux? // Note: momentum flux here doesn't include pressure gradient term From b456558bfe5c63565f765ee0e93942b051c46071 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 1 Sep 2023 21:30:35 -0700 Subject: [PATCH 10/14] neutral_mixed: Modify flux limiters Flux limited particle flux is used in momentum flux before the momentum flux is limited. Particle flux limiter applied to advection of momentum and pressure. --- src/neutral_mixed.cxx | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 1f546ebf..3eed7f2e 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -367,12 +367,12 @@ void NeutralMixed::finally(const Options& state) { } // Flux of parallel momentum - // Question: Should flux-limited particle flux be used here, or original flux? - // Note: momentum flux here doesn't include pressure gradient term + // Note: Flux-limited particle flux is used in calculating the momentum flux before the + // momentum limiter is applied. // // The resulting momentum_flux_factor is applied to the momentum advection and viscosity terms, // and so also scales the advection of kinetic energy - Vector3D momentum_flux = NVn * v_total; + Vector3D momentum_flux = particle_flux_factor * NVn * v_total; if (neutral_viscosity) { momentum_flux -= eta_n * Grad_perp(Vn); } @@ -393,7 +393,7 @@ void NeutralMixed::finally(const Options& state) { // could be in opposite directions. // - Flux doesn't include compression term, or kinetic energy transport // that is in the momentum equation - Vector3D heat_flux = (3./2) * Pn * v_total + Pn * v_perp - kappa_n * Grad(Tn); + Vector3D heat_flux = (3./2) * Pn * particle_flux_factor * v_total + Pn * particle_flux_factor * v_perp - kappa_n * Grad(Tn); Field3D heat_flux_abs = sqrt(heat_flux * heat_flux); Field3D heat_limit = Pnlim * sqrt(2. * Tnlim / (PI * AA)); @@ -418,7 +418,7 @@ void NeutralMixed::finally(const Options& state) { TRACE("Neutral density"); // Note: Parallel and perpendicular flux scaled by limiter - ddt(Nn) = -FV::Div_par_mod(Nn, Vn * particle_flux_factor, sound_speed) // Advection + ddt(Nn) = -FV::Div_par_mod(Nn * particle_flux_factor, Vn, sound_speed) // Advection + FV::Div_a_Grad_perp(DnnNn * particle_flux_factor, logPnlim) // Perpendicular diffusion ; @@ -433,10 +433,11 @@ void NeutralMixed::finally(const Options& state) { TRACE("Neutral momentum"); ddt(NVn) = - -AA * FV::Div_par_fvv(Nnlim, Vn * momentum_flux_factor, sound_speed) // Momentum flow - - Grad_par(Pn) // Pressure gradient. Not included in flux limit. - + FV::Div_a_Grad_perp(DnnNVn * momentum_flux_factor, logPnlim) // Perpendicular diffusion - ; + // Note: The second argument (Vn) is squared in this operator, so the limiters are applied to the first argument + -AA * FV::Div_par_fvv(Nnlim * particle_flux_factor * momentum_flux_factor, Vn, sound_speed) // Momentum flow + - Grad_par(Pn) // Pressure gradient. Not included in flux limit. + + FV::Div_a_Grad_perp(DnnNVn * particle_flux_factor * momentum_flux_factor, logPnlim) // Perpendicular diffusion + ; if (localstate.isSet("momentum_source")) { Snv = get(localstate["momentum_source"]); @@ -447,9 +448,9 @@ void NeutralMixed::finally(const Options& state) { // Neutral pressure TRACE("Neutral pressure"); - ddt(Pn) = -FV::Div_par_mod(Pn, Vn * heat_flux_factor, sound_speed) // Advection + ddt(Pn) = -FV::Div_par_mod(Pn * particle_flux_factor * heat_flux_factor, Vn, sound_speed) // Advection - (2. / 3) * Pn * Div_par(Vn) // Compression - + FV::Div_a_Grad_perp((5. / 3) * DnnPn * heat_flux_factor, logPnlim) // Perpendicular advection: q = 5/2 p u_perp + + FV::Div_a_Grad_perp((5. / 3) * DnnPn * particle_flux_factor * heat_flux_factor, logPnlim) // Perpendicular advection: q = 5/2 p u_perp + (2. / 3) * (FV::Div_a_Grad_perp(kappa_n * heat_flux_factor, Tn) // Perpendicular Conduction + FV::Div_par_K_Grad_par(kappa_n * heat_flux_factor, Tn)) // Parallel conduction ; From eff93f14b65e91a1a9f11fb98cf24424425be517 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 1 Sep 2023 21:32:24 -0700 Subject: [PATCH 11/14] neutral_mixed: Add flux_factor_timescale option Applies a time-averaging to the flux limit factors, with a given timescale. Doesn't appear to have much benefit in quick tests, but left in and turned off by default. May be removed if it really doesn't do anything useful. --- include/neutral_mixed.hxx | 6 ++++++ src/neutral_mixed.cxx | 29 +++++++++++++++++++++++++++++ 2 files changed, 35 insertions(+) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 9a1eaa78..f260a094 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -59,6 +59,12 @@ private: Field3D momentum_flux_factor; Field3D heat_flux_factor; + BoutReal flux_factor_timescale; ///< Timescale over which flux factors vary + BoutReal flux_factor_time {-1.0}; ///< Last time the flux factors were updated + Field3D particle_flux_avg; ///< Moving average flux factor + Field3D momentum_flux_avg; ///< Moving average flux factor + Field3D heat_flux_avg; ///< Moving average flux factor + bool neutral_viscosity; ///< include viscosity? bool precondition {true}; ///< Enable preconditioner? diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 3eed7f2e..5a2a3821 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -54,6 +54,10 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Use isotropic flux limiters?") .withDefault(true); + flux_factor_timescale = options["flux_factor_timescale"] + .doc("Time constant in flux limitation factor variation [seconds]. < 0 is off.") + .withDefault(-1.0) / seconds; + particle_flux_limiter = options["particle_flux_limiter"] .doc("Enable particle flux limiter?") .withDefault(true); @@ -411,6 +415,31 @@ void NeutralMixed::finally(const Options& state) { particle_flux_factor.applyBoundary("neumann"); momentum_flux_factor.applyBoundary("neumann"); heat_flux_factor.applyBoundary("neumann"); + + if (flux_factor_timescale > 0) { + // Smooth flux limitation factors over time + // Suppress rapid changes in flux limitation factors, while retaining + // the same steady-state solution. + BoutReal time = get(state["time"]); + if (flux_factor_time < 0.0) { + // First time, so just copy values + particle_flux_avg = particle_flux_factor; + momentum_flux_avg = momentum_flux_factor; + heat_flux_avg = heat_flux_factor; + } else { + // Weight by a factor that depends on the elapsed time. + BoutReal weight = exp(-(time - flux_factor_time) / flux_factor_timescale); + + particle_flux_avg = weight * particle_flux_avg + (1. - weight) * particle_flux_factor; + momentum_flux_avg = weight * momentum_flux_avg + (1. - weight) * momentum_flux_factor; + heat_flux_avg = weight * heat_flux_avg + (1. - weight) * heat_flux_factor; + + particle_flux_factor = particle_flux_avg; + momentum_flux_factor = momentum_flux_avg; + heat_flux_factor = heat_flux_avg; + } + flux_factor_time = time; + } } ///////////////////////////////////////////////////// From 40243b6fcb4679a2ecf49b208d28d08c7c97b9b8 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 26 Sep 2023 20:57:14 -0700 Subject: [PATCH 12/14] Add solver linear flag to state The linear flag is set when the time integrator is in a linear solve. Components can disable nonlinear updates to (sometimes) improve convergence. --- hermes-3.cxx | 7 ++++--- hermes-3.hxx | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/hermes-3.cxx b/hermes-3.cxx index 4bf08419..3de85697 100644 --- a/hermes-3.cxx +++ b/hermes-3.cxx @@ -244,12 +244,13 @@ int Hermes::init(bool restarting) { return 0; } -int Hermes::rhs(BoutReal time) { +int Hermes::rhs(BoutReal time, bool linear) { // Need to reset the state, since fields may be modified in transform steps state = Options(); - + set(state["time"], time); - state["units"] = units; + state["units"] = units; + set(state["solver"]["linear"], linear); // Call all the components scheduler->transform(state); diff --git a/hermes-3.hxx b/hermes-3.hxx index 3d273c57..b34a44ee 100644 --- a/hermes-3.hxx +++ b/hermes-3.hxx @@ -33,7 +33,7 @@ public: virtual ~Hermes() {} protected: int init(bool restarting) override; - int rhs(BoutReal t) override; + int rhs(BoutReal t, bool linear) override; int precon(BoutReal t, BoutReal gamma, BoutReal delta); /// Add variables to be written to the output file From 35e85b39870ed3a6d0e5f089e98e8f5d912f32a0 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 26 Sep 2023 21:00:47 -0700 Subject: [PATCH 13/14] neutral_mixed: Improvements to AFN flux limiters - Add a relaxation timescale to the coefficients. Damps oscillations and improves convergence. - Calculate fluxes at cell edges, and use maximum local flux to calculate flux limit factors - Near boundaries use the minimum of the flux limit factor in the last two cells. This is because fluxes in the last cell may not be calculated correctly in the flux limit calculation. - Don't update coefficients in a linear solve (this seems to be a small effect) --- src/neutral_mixed.cxx | 502 +++++++++++++++++++++++++++++------------- 1 file changed, 355 insertions(+), 147 deletions(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 5a2a3821..20d4e66e 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -270,178 +270,386 @@ void NeutralMixed::finally(const Options& state) { logPnlim = log(Pnlim); logPnlim.applyBoundary(); - /////////////////////////////////////////////////////// - // Calculate cross-field diffusion from collision frequency - // - // - if (localstate.isSet("collision_frequency")) { - // Dnn = Vth^2 / sigma + // Note: If linear solve then don't recalculate coefficients + if (!get(state["solver"]["linear"])) { + /////////////////////////////////////////////////////// + // Calculate cross-field diffusion from collision frequency + // + // + if (localstate.isSet("collision_frequency")) { + // Dnn = Vth^2 / sigma - if (maximum_mfp > 0) { // MFP limit enabled + if (maximum_mfp > 0) { // MFP limit enabled Field3D Rnn = sqrt(Tn / AA) / (maximum_mfp / get(state["units"]["meters"])); Dnn = (Tn / AA) / (get(localstate["collision_frequency"]) + Rnn); } else { // MFP limit disabled Dnn = (Tn / AA) / (get(localstate["collision_frequency"])); } - - } else { // If no collisions, hardcode max MFP to 0.1m - output_warn.write("No collisions set for the neutrals, limiting mean free path to 0.1m"); - Field3D Rnn = sqrt(Tn / AA) / (0.1 / get(state["units"]["meters"])); - Dnn = (Tn / AA) / Rnn; - } - mesh->communicate(Dnn); - Dnn.clearParallelSlices(); - Dnn.applyBoundary(); - - // Neutral diffusion parameters have the same boundary condition as Dnn - DnnPn = Dnn * Pn; - DnnPn.applyBoundary(); - DnnNn = Dnn * Nn; - DnnNn.applyBoundary(); - Field3D DnnNVn = Dnn * NVn; - DnnNVn.applyBoundary(); + } else { // If no collisions, hardcode max MFP to 0.1m + output_warn.write("No collisions set for the neutrals, limiting mean free path to 0.1m"); + Field3D Rnn = sqrt(Tn / AA) / (0.1 / get(state["units"]["meters"])); + Dnn = (Tn / AA) / Rnn; + } - if (sheath_ydown) { - for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - Dnn(r.ind, mesh->ystart - 1, jz) = -Dnn(r.ind, mesh->ystart, jz); - DnnPn(r.ind, mesh->ystart - 1, jz) = -DnnPn(r.ind, mesh->ystart, jz); - DnnNn(r.ind, mesh->ystart - 1, jz) = -DnnNn(r.ind, mesh->ystart, jz); - DnnNVn(r.ind, mesh->ystart - 1, jz) = -DnnNVn(r.ind, mesh->ystart, jz); + mesh->communicate(Dnn); + Dnn.clearParallelSlices(); + Dnn.applyBoundary(); + + // Neutral diffusion parameters have the same boundary condition as Dnn + DnnPn = Dnn * Pn; + DnnPn.applyBoundary(); + DnnNn = Dnn * Nn; + DnnNn.applyBoundary(); + DnnNVn = Dnn * NVn; + DnnNVn.applyBoundary(); + + if (sheath_ydown) { + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + Dnn(r.ind, mesh->ystart - 1, jz) = -Dnn(r.ind, mesh->ystart, jz); + DnnPn(r.ind, mesh->ystart - 1, jz) = -DnnPn(r.ind, mesh->ystart, jz); + DnnNn(r.ind, mesh->ystart - 1, jz) = -DnnNn(r.ind, mesh->ystart, jz); + DnnNVn(r.ind, mesh->ystart - 1, jz) = -DnnNVn(r.ind, mesh->ystart, jz); + } } } - } - if (sheath_yup) { - for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - Dnn(r.ind, mesh->yend + 1, jz) = -Dnn(r.ind, mesh->yend, jz); - DnnPn(r.ind, mesh->yend + 1, jz) = -DnnPn(r.ind, mesh->yend, jz); - DnnNn(r.ind, mesh->yend + 1, jz) = -DnnNn(r.ind, mesh->yend, jz); - DnnNVn(r.ind, mesh->yend + 1, jz) = -DnnNVn(r.ind, mesh->yend, jz); + if (sheath_yup) { + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + Dnn(r.ind, mesh->yend + 1, jz) = -Dnn(r.ind, mesh->yend, jz); + DnnPn(r.ind, mesh->yend + 1, jz) = -DnnPn(r.ind, mesh->yend, jz); + DnnNn(r.ind, mesh->yend + 1, jz) = -DnnNn(r.ind, mesh->yend, jz); + DnnNVn(r.ind, mesh->yend + 1, jz) = -DnnNVn(r.ind, mesh->yend, jz); + } } } - } - // 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; + // 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 + eta_n = AA * (2. / 5) * kappa_n; + + // Set factors that multiply the fluxes + particle_flux_factor = 1.0; + momentum_flux_factor = 1.0; + heat_flux_factor = 1.0; + if (flux_limit) { + // Apply flux limiters + // Note: Fluxes calculated here are cell centre, rather than cell edge + + // Cross-field velocity + Vector3D v_perp = -Dnn * Grad_perp(logPnlim); + + // Total velocity, perpendicular + parallel + Vector3D v_total = v_perp; + ASSERT2(v_total.covariant == true); + auto* coord = mesh->getCoordinates(); + v_total.y = Vn * (coord->J * coord->Bxy); // Set parallel component + + Field3D v_sq = v_total * v_total; // v dot v + Field3D v_abs = sqrt(v_sq); // |v| + + // Magnitude of the particle flux + Field3D particle_flux_abs = Nnlim * v_abs; + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + // Calculate flux from i to i+1 + + BoutReal fout = 0.5 * (DnnNn(i, j, k) + DnnNn(i + 1, j, k)) + * (sqrt(coord->g11(i, j, k)) + + sqrt(coord->g11(i + 1, j, k))) + * (logPnlim(i + 1, j, k) - logPnlim(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); + + // Calculate flux from i to i-1 + BoutReal fin = 0.5 * (DnnNn(i, j, k) + DnnNn(i - 1, j, k)) + * (sqrt(coord->g11(i, j, k)) + + sqrt(coord->g11(i - 1, j, k))) + * (logPnlim(i - 1, j, k) - logPnlim(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i - 1, j, k)); + + // Flux from j to j+1 + BoutReal fup = 0.5 * (DnnNn(i, j, k) + DnnNn(i, j + 1, k)) + * (sqrt(coord->g22(i, j, k)) + + sqrt(coord->g22(i, j + 1, k))) + * (logPnlim(i, j + 1, k) - logPnlim(i, j, k)) + / (coord->dy(i, j, k) + coord->dy(i, j + 1, k)); + + // Flux from j to j-1 + BoutReal fdown = 0.5 * (DnnNn(i, j, k) + DnnNn(i, j - 1, k)) + * (sqrt(coord->g22(i, j, k)) + + sqrt(coord->g22(i, j - 1, k))) + * (logPnlim(i, j - 1, k) - logPnlim(i, j, k)) + / (coord->dy(i, j, k) + coord->dy(i, j - 1, k)); + + // Choose the maximum local particle flux, + // using cell center and edge values + particle_flux_abs(i, j, k) = BOUTMAX(particle_flux_abs(i, j, k), + fabs(fin), + fabs(fout), + fabs(fup), + fabs(fdown)); + } + } + } - // Viscosity - eta_n = AA * (2. / 5) * kappa_n; + // Normalised particle flux limit + Field3D particle_limit = Nnlim * 0.25 * sqrt(8 * Tnlim / (PI * AA)); - // Sound speed appearing in Lax flux for advection terms - Field3D sound_speed = sqrt(Tn * (5. / 3) / AA); + // Particle flux reduction factor + if (particle_flux_limiter) { + particle_flux_factor = pow(1. + pow(particle_flux_abs / (flux_limit_alpha * particle_limit), + flux_limit_gamma), + -1./flux_limit_gamma); + } else { + particle_flux_factor = 1.0; + } - // Set factors that multiply the fluxes - particle_flux_factor = 1.0; - momentum_flux_factor = 1.0; - heat_flux_factor = 1.0; - if (flux_limit) { - // Apply flux limiters - // Note: Fluxes calculated here are cell centre, rather than cell edge - - // Cross-field velocity - Vector3D v_perp = -Dnn * Grad_perp(logPnlim); - - // Total velocity, perpendicular + parallel - Vector3D v_total = v_perp; - ASSERT2(v_total.covariant == true); - auto* coord = mesh->getCoordinates(); - v_total.y = Vn * (coord->J * coord->Bxy); // Set parallel component - - Field3D v_sq = v_perp * v_perp + SQ(Vn); // v dot v - Field3D v_abs = sqrt(v_sq); // |v| - - // Magnitude of the particle flux - Field3D particle_flux_abs = Nnlim * v_abs; - - // Normalised particle flux limit - Field3D particle_limit = Nnlim * 0.25 * sqrt(8 * Tnlim / (PI * AA)); - - // Particle flux reduction factor - if (particle_flux_limiter){ - particle_flux_factor = pow(1. + pow(particle_flux_abs / (flux_limit_alpha * particle_limit), - flux_limit_gamma), - -1./flux_limit_gamma); - } else { - particle_flux_factor = 1.0; - } + // Flux of parallel momentum + // Note: Flux-limited particle flux is used in calculating the momentum flux before the + // momentum limiter is applied. + // + // The resulting momentum_flux_factor is applied to the momentum advection and viscosity terms, + // and so also scales the advection of kinetic energy + Vector3D momentum_flux = particle_flux_factor * NVn * v_total; + if (neutral_viscosity) { + momentum_flux -= eta_n * Grad_perp(Vn); + } + Field3D momentum_flux_abs = sqrt(momentum_flux * momentum_flux); + Field3D momentum_limit = Pnlim; - // Flux of parallel momentum - // Note: Flux-limited particle flux is used in calculating the momentum flux before the - // momentum limiter is applied. - // - // The resulting momentum_flux_factor is applied to the momentum advection and viscosity terms, - // and so also scales the advection of kinetic energy - Vector3D momentum_flux = particle_flux_factor * NVn * v_total; - if (neutral_viscosity) { - momentum_flux -= eta_n * Grad_perp(Vn); - } - Field3D momentum_flux_abs = sqrt(momentum_flux * momentum_flux); - Field3D momentum_limit = Pnlim; - - if (momentum_flux_limiter) { - momentum_flux_factor = pow(1. + pow(momentum_flux_abs / (flux_limit_alpha * momentum_limit), - flux_limit_gamma), - -1./flux_limit_gamma); - } else { - momentum_flux_factor = 1.0; - } - // Flux of heat - // Note: - // - Limiting the heat flux, not energy flux - // - Advection and heat conduction are both vectors, and e.g - // could be in opposite directions. - // - Flux doesn't include compression term, or kinetic energy transport - // that is in the momentum equation - Vector3D heat_flux = (3./2) * Pn * particle_flux_factor * v_total + Pn * particle_flux_factor * v_perp - kappa_n * Grad(Tn); - - Field3D heat_flux_abs = sqrt(heat_flux * heat_flux); - Field3D heat_limit = Pnlim * sqrt(2. * Tnlim / (PI * AA)); - - if (heat_flux_limiter) { - heat_flux_factor = pow(1. + pow(heat_flux_abs / (flux_limit_alpha * heat_limit), - flux_limit_gamma), - -1./flux_limit_gamma); - } else { - heat_flux_factor = 1.0; - } - // Communicate guard cells and apply boundary conditions - // because the flux factors will be differentiated - mesh->communicate(particle_flux_factor, momentum_flux_factor, heat_flux_factor); - particle_flux_factor.applyBoundary("neumann"); - momentum_flux_factor.applyBoundary("neumann"); - heat_flux_factor.applyBoundary("neumann"); - - if (flux_factor_timescale > 0) { - // Smooth flux limitation factors over time - // Suppress rapid changes in flux limitation factors, while retaining - // the same steady-state solution. - BoutReal time = get(state["time"]); - if (flux_factor_time < 0.0) { - // First time, so just copy values - particle_flux_avg = particle_flux_factor; - momentum_flux_avg = momentum_flux_factor; - heat_flux_avg = heat_flux_factor; + if (momentum_flux_limiter) { + momentum_flux_factor = pow(1. + pow(momentum_flux_abs / (flux_limit_alpha * momentum_limit), + flux_limit_gamma), + -1./flux_limit_gamma); } else { - // Weight by a factor that depends on the elapsed time. - BoutReal weight = exp(-(time - flux_factor_time) / flux_factor_timescale); + momentum_flux_factor = 1.0; + } + // Flux of heat + // Note: + // - Limiting the heat flux, not energy flux + // - Advection and heat conduction are both vectors, and e.g + // could be in opposite directions. + // - Flux doesn't include compression term, or kinetic energy transport + // that is in the momentum equation + Vector3D heat_flux = (3./2) * Pn * particle_flux_factor * v_total + + Pn * particle_flux_factor * v_perp + - kappa_n * Grad(Tn); + + Field3D heat_flux_abs = sqrt(heat_flux * heat_flux); + + mesh->communicate(particle_flux_factor); + particle_flux_factor.applyBoundary("neumann"); + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + // Calculate energy flux from i to i+1 + + BoutReal fout = + // Convection with particle flux limiter + (5./2) * 0.5 * (particle_flux_factor(i, j, k) * DnnPn(i, j, k) + + particle_flux_factor(i + 1, j, k) * DnnPn(i + 1, j, k)) + * (sqrt(coord->g11(i, j, k)) + + sqrt(coord->g11(i + 1, j, k))) + * (logPnlim(i + 1, j, k) - logPnlim(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)) + + + // Heat conduction + 0.5 * (kappa_n(i, j, k) + kappa_n(i + 1, j, k)) + * (sqrt(coord->g11(i, j, k)) + + sqrt(coord->g11(i + 1, j, k))) + * (Tn(i + 1, j, k) - Tn(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)) + ; + + // Calculate flux from i to i-1 + BoutReal fin = + // Convection + (5./2) * 0.5 * (particle_flux_factor(i, j, k) * DnnPn(i, j, k) + + particle_flux_factor(i - 1, j, k) * DnnPn(i - 1, j, k)) + * (sqrt(coord->g11(i, j, k)) + + sqrt(coord->g11(i - 1, j, k))) + * (logPnlim(i - 1, j, k) - logPnlim(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i - 1, j, k)) + + + // Heat conduction + 0.5 * (kappa_n(i, j, k) + kappa_n(i - 1, j, k)) + * (sqrt(coord->g11(i, j, k)) + + sqrt(coord->g11(i - 1, j, k))) + * (Tn(i - 1, j, k) - Tn(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)) + ; + + // Flux from j to j+1 + BoutReal fup = + // Convection + (5./2) * 0.5 * (particle_flux_factor(i, j, k) * DnnPn(i, j, k) + + particle_flux_factor(i, j + 1, k) * DnnPn(i, j + 1, k)) + * (sqrt(coord->g22(i, j, k)) + + sqrt(coord->g22(i, j + 1, k))) + * (logPnlim(i, j + 1, k) - logPnlim(i, j, k)) + / (coord->dy(i, j, k) + coord->dy(i, j + 1, k)) + + + // Heat conduction + 0.5 * (kappa_n(i, j, k) + kappa_n(i, j + 1, k)) + * (sqrt(coord->g22(i, j, k)) + + sqrt(coord->g22(i, j + 1, k))) + * (Tn(i, j + 1, k) - Tn(i, j, k)) + / (coord->dy(i, j, k) + coord->dy(i, j + 1, k)) + ; + + // Flux from j to j-1 + BoutReal fdown = + // Convection + (5./2) * 0.5 * (particle_flux_factor(i, j, k) * DnnPn(i, j, k) + + particle_flux_factor(i, j - 1, k) * DnnPn(i, j - 1, k)) + * (sqrt(coord->g22(i, j, k)) + + sqrt(coord->g22(i, j - 1, k))) + * (logPnlim(i, j - 1, k) - logPnlim(i, j, k)) + / (coord->dy(i, j, k) + coord->dy(i, j - 1, k)) + + + // Heat conduction + 0.5 * (kappa_n(i, j, k) + kappa_n(i, j - 1, k)) + * (sqrt(coord->g22(i, j, k)) + + sqrt(coord->g22(i, j - 1, k))) + * (Tn(i, j - 1, k) - Tn(i, j, k)) + / (coord->dy(i, j, k) + coord->dy(i, j - 1, k)) + ; + + heat_flux_abs(i, j, k) = BOUTMAX(heat_flux_abs(i, j, k), + fabs(fin), + fabs(fout), + fabs(fup), + fabs(fdown)); + } + } + } - particle_flux_avg = weight * particle_flux_avg + (1. - weight) * particle_flux_factor; - momentum_flux_avg = weight * momentum_flux_avg + (1. - weight) * momentum_flux_factor; - heat_flux_avg = weight * heat_flux_avg + (1. - weight) * heat_flux_factor; + Field3D heat_limit = Pnlim * sqrt(2. * Tnlim / (PI * AA)); - particle_flux_factor = particle_flux_avg; - momentum_flux_factor = momentum_flux_avg; - heat_flux_factor = heat_flux_avg; + if (heat_flux_limiter) { + heat_flux_factor = pow(1. + pow(heat_flux_abs / (flux_limit_alpha * heat_limit), + flux_limit_gamma), + -1./flux_limit_gamma); + } else { + heat_flux_factor = 1.0; + } + // Communicate guard cells and apply boundary conditions + // because the flux factors will be differentiated + mesh->communicate(particle_flux_factor, momentum_flux_factor, heat_flux_factor); + + // Fluxes at the boundary are not calculated correctly, so + // use limiting factors from one cell inwards + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + particle_flux_factor(r.ind, mesh->ystart, jz) = + BOUTMIN(particle_flux_factor(r.ind, mesh->ystart, jz), + particle_flux_factor(r.ind, mesh->ystart + 1, jz)); + momentum_flux_factor(r.ind, mesh->ystart, jz) = + BOUTMIN(momentum_flux_factor(r.ind, mesh->ystart, jz), + momentum_flux_factor(r.ind, mesh->ystart + 1, jz)); + heat_flux_factor(r.ind, mesh->ystart, jz) = + BOUTMIN(heat_flux_factor(r.ind, mesh->ystart, jz), + heat_flux_factor(r.ind, mesh->ystart + 1, jz)); + + // Neumann boundary conditions + particle_flux_factor(r.ind, mesh->ystart - 1, jz) = particle_flux_factor(r.ind, mesh->ystart, jz); + momentum_flux_factor(r.ind, mesh->ystart - 1, jz) = momentum_flux_factor(r.ind, mesh->ystart, jz); + heat_flux_factor(r.ind, mesh->ystart - 1, jz) = heat_flux_factor(r.ind, mesh->ystart, jz); + } + } + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + particle_flux_factor(r.ind, mesh->yend, jz) = + BOUTMIN(particle_flux_factor(r.ind, mesh->yend, jz), + particle_flux_factor(r.ind, mesh->yend - 1, jz)); + momentum_flux_factor(r.ind, mesh->yend, jz) = + BOUTMIN(momentum_flux_factor(r.ind, mesh->yend, jz), + momentum_flux_factor(r.ind, mesh->yend - 1, jz)); + heat_flux_factor(r.ind, mesh->yend, jz) = + BOUTMIN(heat_flux_factor(r.ind, mesh->yend, jz), + heat_flux_factor(r.ind, mesh->yend - 1, jz)); + + // Neumann boundary conditions + particle_flux_factor(r.ind, mesh->yend + 1, jz) = particle_flux_factor(r.ind, mesh->yend, jz); + momentum_flux_factor(r.ind, mesh->yend + 1, jz) = momentum_flux_factor(r.ind, mesh->yend, jz); + heat_flux_factor(r.ind, mesh->yend + 1, jz) = heat_flux_factor(r.ind, mesh->yend, jz); + } + } + + if (mesh->firstX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + particle_flux_factor(mesh->xstart, j, k) = + BOUTMIN(particle_flux_factor(mesh->xstart, j, k), + particle_flux_factor(mesh->xstart + 1, j, k)); + momentum_flux_factor(mesh->xstart, j, k) = + BOUTMIN(momentum_flux_factor(mesh->xstart, j, k), + momentum_flux_factor(mesh->xstart + 1, j, k)); + heat_flux_factor(mesh->xstart, j, k) = + BOUTMIN(heat_flux_factor(mesh->xstart, j, k), + heat_flux_factor(mesh->xstart + 1, j, k)); + } + } + } + if (mesh->lastX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + particle_flux_factor(mesh->xend, j, k) = + BOUTMIN(particle_flux_factor(mesh->xend, j, k), + particle_flux_factor(mesh->xend - 1, j, k)); + momentum_flux_factor(mesh->xend, j, k) = + BOUTMIN(momentum_flux_factor(mesh->xend, j, k), + momentum_flux_factor(mesh->xend - 1, j, k)); + heat_flux_factor(mesh->xend, j, k) = + BOUTMIN(heat_flux_factor(mesh->xend, j, k), + heat_flux_factor(mesh->xend - 1, j, k)); + } + } + } + + particle_flux_factor.applyBoundary("neumann"); + momentum_flux_factor.applyBoundary("neumann"); + heat_flux_factor.applyBoundary("neumann"); + + if (flux_factor_timescale > 0) { + // Smooth flux limitation factors over time + // Suppress rapid changes in flux limitation factors, while retaining + // the same steady-state solution. + BoutReal time = get(state["time"]); + if (flux_factor_time < 0.0) { + // First time, so just copy values + particle_flux_avg = particle_flux_factor; + momentum_flux_avg = momentum_flux_factor; + heat_flux_avg = heat_flux_factor; + } else { + // Weight by a factor that depends on the elapsed time. + BoutReal weight = exp(-(time - flux_factor_time) / flux_factor_timescale); + + particle_flux_avg = weight * particle_flux_avg + (1. - weight) * particle_flux_factor; + momentum_flux_avg = weight * momentum_flux_avg + (1. - weight) * momentum_flux_factor; + heat_flux_avg = weight * heat_flux_avg + (1. - weight) * heat_flux_factor; + + particle_flux_factor = particle_flux_avg; + momentum_flux_factor = momentum_flux_avg; + heat_flux_factor = heat_flux_avg; + } + flux_factor_time = time; } - flux_factor_time = time; } } + // Sound speed appearing in Lax flux for advection terms + Field3D sound_speed = sqrt(Tn * (5. / 3) / AA); + ///////////////////////////////////////////////////// // Neutral density TRACE("Neutral density"); From 3d005400b070b2c3ae7c47350178340d83ce3940 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 26 Sep 2023 21:05:31 -0700 Subject: [PATCH 14/14] neutral_mixed: Change defaults in flux limiter - Disable maximum_mfp by default - Use a timescale of 1e-6 for the flux limit factor relaxation (`flux_factor_timescale` option). --- src/neutral_mixed.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 20d4e66e..db64b818 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -56,7 +56,7 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* flux_factor_timescale = options["flux_factor_timescale"] .doc("Time constant in flux limitation factor variation [seconds]. < 0 is off.") - .withDefault(-1.0) / seconds; + .withDefault(1e-6) / seconds; particle_flux_limiter = options["particle_flux_limiter"] .doc("Enable particle flux limiter?") @@ -84,7 +84,7 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* maximum_mfp = options["maximum_mfp"] .doc("Optional maximum mean free path in [m] for diffusive processes. < 0 is off") - .withDefault(0.1); + .withDefault(-1.0); if (precondition) { inv = std::unique_ptr(Laplacian::create(&options["precon_laplace"]));