From 5d0b9ca0670382ceb7352bf95a9bc251ba686de0 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 20 May 2026 15:47:49 -0700 Subject: [PATCH 1/2] relax_potential: Add evolve_vorticity switch Rather than evolving both potential and vorticity equations, setting `evolve_vorticity` to `false` evolves only potential. --- include/relax_potential.hxx | 42 +++++++---- src/relax_potential.cxx | 146 ++++++++++++++++++++---------------- 2 files changed, 108 insertions(+), 80 deletions(-) diff --git a/include/relax_potential.hxx b/include/relax_potential.hxx index 227a78942..0fff4598a 100644 --- a/include/relax_potential.hxx +++ b/include/relax_potential.hxx @@ -2,10 +2,16 @@ #ifndef RELAX_POTENTIAL_H #define RELAX_POTENTIAL_H +#include "../include/guarded_options.hxx" + +#include +#include #include #include "component.hxx" +#include + /// Evolve vorticity and potential in time. /// /// Uses a relaxation method for the potential, which is valid for @@ -80,8 +86,14 @@ struct RelaxPotential : public Component { // {{"long_name", "plasma potential"}, // {"source", "vorticity"}}); // } + + /// Calculate vorticity from potential and species + Field3D vorticity(const Field3D& phi, GuardedOptions& allspecies); + private: - Field3D Vort; // Evolving vorticity + bool evolve_vorticity; ///< Evolve vorticity? + Field3D Vort; // Evolving vorticity + Field3D Vort_from_phi; ///< Vort calculated from phi Field3D phi1; // Scaled electrostatic potential, evolving in time ϕ_1 = λ_2 ϕ Field3D phi; // Electrostatic potential @@ -102,25 +114,25 @@ private: bool sheath_boundary; ///< Set outer boundary to j=0? - bool vort_dissipation; ///< Parallel dissipation of vorticity - bool phi_dissipation; ///< Parallel dissipation of potential + bool vort_dissipation; ///< Parallel dissipation of vorticity + bool phi_dissipation; ///< Parallel dissipation of potential bool phi_sheath_dissipation; ///< Dissipation at the sheath if phi < 0 - bool damp_core_vorticity; ///< Damp axisymmetric component of vorticity + bool damp_core_vorticity; ///< Damp axisymmetric component of vorticity - bool phi_boundary_relax; ///< Relax boundary to zero-gradient - BoutReal phi_boundary_timescale; ///< Relaxation timescale [normalised] + bool phi_boundary_relax; ///< Relax boundary to zero-gradient + BoutReal phi_boundary_timescale; ///< Relaxation timescale [normalised] BoutReal phi_boundary_last_update; ///< Time when last updated - bool phi_core_averagey; ///< Average phi core boundary in Y? + bool phi_core_averagey; ///< Average phi core boundary in Y? - Field2D Bsq; ///< SQ(coord->Bxy) - Vector2D Curlb_B; ///< Curvature vector Curl(b/B) - BoutReal hyper_z; ///< Hyper-viscosity in Z - Field2D viscosity; ///< Perpendicular Kinematic viscosity - Field2D viscosity_par; ///< Parallel Kinematic viscosity + Field2D Bsq; ///< SQ(coord->Bxy) + Vector2D Curlb_B; ///< Curvature vector Curl(b/B) + BoutReal hyper_z; ///< Hyper-viscosity in Z + Field2D viscosity; ///< Perpendicular Kinematic viscosity + Field2D viscosity_par; ///< Parallel Kinematic viscosity - // Relax-potential related variables - BoutReal lambda_1; ///< Relaxation parameters. NOTE: lambda_1 has dimensions! - BoutReal lambda_2; ///< Relaxation parameters + // Relax-potential related variables + BoutReal lambda_1; ///< Relaxation parameters. NOTE: lambda_1 has dimensions! + BoutReal lambda_2; ///< Relaxation parameters // Diagnostic outputs Field3D DivJdia, DivJcol; // Divergence of diamagnetic and collisional current diff --git a/src/relax_potential.cxx b/src/relax_potential.cxx index 3c4c53faf..31c3a1840 100644 --- a/src/relax_potential.cxx +++ b/src/relax_potential.cxx @@ -20,9 +20,11 @@ #include #include #include +#include #include #include +#include #include using bout::globals::mesh; @@ -58,10 +60,15 @@ BoutReal limitFree(BoutReal fm, BoutReal fc) { RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* solver) : Component({readWrite("fields:vorticity"), readWrite("fields:phi")}) { - solver->add(Vort, "Vort"); // Vorticity evolving + auto& options = alloptions[name]; + + evolve_vorticity = options["evolve_vorticity"].doc("").withDefault(true); + + if (evolve_vorticity) { + solver->add(Vort, "Vort"); // Vorticity evolving + } solver->add(phi1, "phi1"); // Evolving scaled potential ϕ_1 = λ_2 ϕ - auto& options = alloptions[name]; // Normalisations const Options& units = alloptions["units"]; const BoutReal Omega_ci = 1. / units["seconds"].as(); @@ -206,10 +213,9 @@ RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* so if (diamagnetic) { // Need curvature throw; - } else { - output_warn.write("No curvature vector in input grid"); - Curlb_B = 0.0; } + output_warn.write("No curvature vector in input grid"); + Curlb_B = 0.0; } } @@ -245,9 +251,15 @@ void RelaxPotential::transform_impl(GuardedOptions& state) { // Scale potential phi = phi1 / lambda_2; - // Set the boundary of Vort. . Needed only if dissipation terms are included. + // Calculate vorticity from phi. + this->Vort_from_phi = vorticity(phi, allspecies); + if (!evolve_vorticity) { + // Not evolving vorticity so calculate from phi + Vort = this->Vort_from_phi; + } + + // Set the boundary of Vort. Needed only if dissipation terms are included. Vort.applyBoundary("neumann"); - // Vort.applyBoundary("dirichlet"); if (Vort.hasParallelSlices()) { Field3D& Vort_ydown = Vort.ydown(); @@ -287,9 +299,6 @@ void RelaxPotential::transform_impl(GuardedOptions& state) { Vort = fromFieldAligned(Vort_fa); } - // Set the boundary of phi. - // phi.applyBoundary("neumann"); - // Note: For now the boundary values are all at the midpoint, // and only phi is considered, not phi + Pi which is handled in Boussinesq solves Pi_hat = 0.0; // Contribution from ion pressure, weighted by atomic mass / charge @@ -642,13 +651,9 @@ void RelaxPotential::transform_impl(GuardedOptions& state) { } void RelaxPotential::finally(const Options& state) { - - const Options& allspecies = state["species"]; - const auto* coord = mesh->getCoordinates(); phi = get(state["fields"]["phi"]); - // Vort = get(state["fields"]["vorticity"]); // Solve vorticity equation @@ -787,60 +792,16 @@ void RelaxPotential::finally(const Options& state) { } } - // Solve diffusion equation for potential - - if (boussinesq) { - ddt(phi1) = lambda_1 * (FV::Div_a_Grad_perp(average_atomic_mass / Bsq, phi) - Vort); + if (evolve_vorticity) { + // Solve diffusion equation for potential, relaxing towards + // the evolving vorticity. Uses Vort_from_phi calculated in transform() + ddt(phi1) = lambda_1 * (this->Vort_from_phi - Vort); - if (diamagnetic_polarisation) { - for (const auto& kv : allspecies.getChildren()) { - // Note: includes electrons (should it?) - - const Options& species = kv.second; - if (!species.isSet("charge")) { - continue; // Not charged - } - const BoutReal Z = get(species["charge"]); - if (fabs(Z) < 1e-5) { - continue; // Not charged - } - if (!species.isSet("pressure")) { - continue; // No pressure - } - const BoutReal A = get(species["AA"]); - const Field3D P = get(species["pressure"]); - ddt(phi1) += lambda_1 * FV::Div_a_Grad_perp(A / Bsq, P); - } - } } else { - // Non-Boussinesq. Calculate mass density by summing over species - - // Calculate vorticity from potential phi - Field3D phi_vort = 0.0; - for (const auto& kv : allspecies.getChildren()) { - const Options& species = kv.second; - - if (!species.isSet("charge")) { - continue; // Not charged - } - const BoutReal Zi = get(species["charge"]); - if (fabs(Zi) < 1e-5) { - continue; // Not charged - } - - const BoutReal Ai = get(species["AA"]); - const Field3D Ni = get(species["density"]); - - phi_vort += FV::Div_a_Grad_perp((Ai / Bsq) * Ni, phi); + // Evolve potential not vorticity. + // This is only valid in steady state (ddt -> 0) - if (diamagnetic_polarisation and species.isSet("pressure")) { - // Calculate the diamagnetic flow contribution - const Field3D Pi = get(species["pressure"]); - phi_vort += FV::Div_a_Grad_perp(Ai / Bsq / Zi, Pi); - } - } - - ddt(phi1) = lambda_1 * (phi_vort - Vort); + ddt(phi1) = -lambda_1 * ddt(Vort); } } @@ -895,3 +856,58 @@ void RelaxPotential::outputVars(Options& state) { } } } + +Field3D RelaxPotential::vorticity(const Field3D& phi, GuardedOptions& allspecies) { + if (boussinesq) { + Field3D phi_vort = FV::Div_a_Grad_perp(average_atomic_mass / Bsq, phi); + + if (diamagnetic_polarisation) { + for (const auto& kv : allspecies.getChildren()) { + // Note: includes electrons (should it?) + + const GuardedOptions& species = kv.second; + if (!species.isSet("charge")) { + continue; // Not charged + } + const BoutReal Z = get(species["charge"]); + if (fabs(Z) < 1e-5) { + continue; // Not charged + } + if (!species.isSet("pressure")) { + continue; // No pressure + } + const BoutReal A = get(species["AA"]); + const Field3D P = get(species["pressure"]); + phi_vort += FV::Div_a_Grad_perp(A / Bsq, P); + } + } + return phi_vort; + } + // Non-Boussinesq. Calculate mass density by summing over species + + // Calculate vorticity from potential phi + Field3D phi_vort = 0.0; + for (const auto& kv : allspecies.getChildren()) { + const GuardedOptions& species = kv.second; + + if (!species.isSet("charge")) { + continue; // Not charged + } + const BoutReal Zi = get(species["charge"]); + if (fabs(Zi) < 1e-5) { + continue; // Not charged + } + + const BoutReal Ai = get(species["AA"]); + const Field3D Ni = get(species["density"]); + + phi_vort += FV::Div_a_Grad_perp((Ai / Bsq) * Ni, phi); + + if (diamagnetic_polarisation and species.isSet("pressure")) { + // Calculate the diamagnetic flow contribution + const Field3D Pi = get(species["pressure"]); + phi_vort += FV::Div_a_Grad_perp(Ai / Bsq / Zi, Pi); + } + } + return phi_vort; +} From 71d2c008d6f2a7207450d4416f45944cace17f4e Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 22 May 2026 09:02:29 -0700 Subject: [PATCH 2/2] relax_potential: Update permissions Vorticity is now calculated from phi in `transform` rather than `finally`, so additional permissions are needed. --- src/relax_potential.cxx | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/relax_potential.cxx b/src/relax_potential.cxx index 31c3a1840..0024d20a8 100644 --- a/src/relax_potential.cxx +++ b/src/relax_potential.cxx @@ -58,7 +58,12 @@ BoutReal limitFree(BoutReal fm, BoutReal fc) { } // namespace RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* solver) - : Component({readWrite("fields:vorticity"), readWrite("fields:phi")}) { + : Component({ + readWrite("fields:vorticity"), + readWrite("fields:phi"), + readIfSet("species:{all_species}:charge"), + readOnly("species:{charged}:AA"), + }) { auto& options = alloptions[name]; @@ -189,6 +194,8 @@ RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* so .doc("Timescale for phi boundary relaxation [seconds]") .withDefault(1e-4) / get(alloptions["units"]["seconds"]); + } else { + setPermissions(readIfSet("species:{charged}:temperature", Regions::Interior)); } if (diamagnetic) { @@ -877,7 +884,7 @@ Field3D RelaxPotential::vorticity(const Field3D& phi, GuardedOptions& allspecies continue; // No pressure } const BoutReal A = get(species["AA"]); - const Field3D P = get(species["pressure"]); + const Field3D P = GET_NOBOUNDARY(Field3D, species["pressure"]); phi_vort += FV::Div_a_Grad_perp(A / Bsq, P); } }