From de487f833448dcce814aad0dcef63055c05d42ef Mon Sep 17 00:00:00 2001 From: malamast Date: Tue, 30 Dec 2025 10:52:19 -0800 Subject: [PATCH 1/4] vorticity: precalculate Grad(phi) instead of calculating it for everyspecies. --- src/vorticity.cxx | 65 +++++++++++++++++++++++++---------------------- 1 file changed, 34 insertions(+), 31 deletions(-) diff --git a/src/vorticity.cxx b/src/vorticity.cxx index e663d73e7..b4ed63d23 100644 --- a/src/vorticity.cxx +++ b/src/vorticity.cxx @@ -482,6 +482,39 @@ void Vorticity::transform(Options& state) { Jdia.z = 0.0; Jdia.covariant = Curlb_B.covariant; + // Pre-calculate this rather than calculate for each species + // Note: The below calculation requires phi derivatives at the Y boundaries + // Setting to free boundaries + if (phi.hasParallelSlices()) { + Field3D &phi_ydown = phi.ydown(); + Field3D &phi_yup = phi.yup(); + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + phi_ydown(r.ind, mesh->ystart - 1, jz) = 2 * phi(r.ind, mesh->ystart, jz) - phi_yup(r.ind, mesh->ystart + 1, jz); + } + } + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + phi_yup(r.ind, mesh->yend + 1, jz) = 2 * phi(r.ind, mesh->yend, jz) - phi_ydown(r.ind, mesh->yend - 1, jz); + } + } + } else { + Field3D phi_fa = toFieldAligned(phi); + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + phi_fa(r.ind, mesh->ystart - 1, jz) = 2 * phi_fa(r.ind, mesh->ystart, jz) - phi_fa(r.ind, mesh->ystart + 1, jz); + } + } + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + phi_fa(r.ind, mesh->yend + 1, jz) = 2 * phi_fa(r.ind, mesh->yend, jz) - phi_fa(r.ind, mesh->yend - 1, jz); + } + } + phi = fromFieldAligned(phi_fa); + } + + Vector3D Grad_phi = Grad(phi); + Options& allspecies = state["species"]; for (auto& kv : allspecies.getChildren()) { @@ -532,41 +565,11 @@ void Vorticity::transform(Options& state) { P = fromFieldAligned(P_fa); } - // Note: This calculation requires phi derivatives at the Y boundaries - // Setting to free boundaries - if (phi.hasParallelSlices()) { - Field3D &phi_ydown = phi.ydown(); - Field3D &phi_yup = phi.yup(); - for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - phi_ydown(r.ind, mesh->ystart - 1, jz) = 2 * phi(r.ind, mesh->ystart, jz) - phi_yup(r.ind, mesh->ystart + 1, jz); - } - } - for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - phi_yup(r.ind, mesh->yend + 1, jz) = 2 * phi(r.ind, mesh->yend, jz) - phi_ydown(r.ind, mesh->yend - 1, jz); - } - } - } else { - Field3D phi_fa = toFieldAligned(phi); - for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - phi_fa(r.ind, mesh->ystart - 1, jz) = 2 * phi_fa(r.ind, mesh->ystart, jz) - phi_fa(r.ind, mesh->ystart + 1, jz); - } - } - for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - phi_fa(r.ind, mesh->yend + 1, jz) = 2 * phi_fa(r.ind, mesh->yend, jz) - phi_fa(r.ind, mesh->yend - 1, jz); - } - } - phi = fromFieldAligned(phi_fa); - } - Vector3D Jdia_species = P * Curlb_B; // Diamagnetic current for this species // This term energetically balances diamagnetic term // in the vorticity equation - subtract(species["energy_source"], Jdia_species * Grad(phi)); + subtract(species["energy_source"], Jdia_species * Grad_phi); Jdia += Jdia_species; // Collect total diamagnetic current } From 10159064e721dd6bf381c5c3dda418c32bb13171 Mon Sep 17 00:00:00 2001 From: malamast Date: Thu, 8 Jan 2026 15:29:53 -0800 Subject: [PATCH 2/4] relax_potential: major revision. Coppied many functionalities frm the vorticity component. --- include/relax_potential.hxx | 92 ++++- src/relax_potential.cxx | 726 +++++++++++++++++++++++++++++++++--- 2 files changed, 755 insertions(+), 63 deletions(-) diff --git a/include/relax_potential.hxx b/include/relax_potential.hxx index bb6668d54..f35299ea5 100644 --- a/include/relax_potential.hxx +++ b/include/relax_potential.hxx @@ -16,15 +16,42 @@ struct RelaxPotential : public Component { /// Options /// /// - - /// - diamagnetic - /// - diamagnetic_polarisation - /// - average_atomic_mass - /// - bndry_flux - /// - poloidal_flows - /// - split_n0 - /// - laplacian - /// Options for the Laplacian phi solver + /// - average_atomic_mass: float, default 2.0 + /// Weighted average ion atomic mass for polarisation current + /// - bndry_flux: bool, default true + /// Allow flows through radial (X) boundaries? + /// - collisional_friction: bool, default false + /// Damp vorticity based on mass-weighted collision frequency? + /// - diagnose: bool, false + /// Output additional diagnostics? + /// - diamagnetic: bool, default true + /// Include diamagnetic current, using curvature vector? + /// - diamagnetic_polarisation: bool, default true + /// Include ion diamagnetic drift in polarisation current? + /// - exb_advection: bool, default true + /// Include ExB advection (nonlinear term)? + /// - phi_boundary_relax: bool, default false + /// Relax radial phi boundaries towards zero-gradient? + /// - phi_boundary_timescale: float, 1e-4 + /// Timescale for phi boundary relaxation [seconds] + /// - phi_core_averagey: bool, default false + /// Average phi core boundary in Y? (if phi_boundary_relax) + /// - phi_dissipation: bool, default true + /// Parallel dissipation of potential (Recommended) + /// - poloidal_flows: bool, default true + /// Include poloidal ExB flow? + /// - sheath_boundary: bool, default false + /// If phi_boundary_relax is false, set the radial boundary to the sheath potential? + /// - viscosity_perp: Field2D, default 0.0 + /// Kinematic viscosity in perpendicular diraction [m^2/s] + /// - viscosity_par: Field2D, default 0.0 + /// Kinematic viscosity in parallel diraction [m^2/s] + /// - vort_dissipation: bool, default false + /// Parallel dissipation of vorticity? + /// - damp_core_vorticity: bool, default false + /// Damp axisymmetric component of vorticity in cell next to core boundary /// + RelaxPotential(std::string name, Options& options, Solver* solver); /// Optional inputs @@ -55,13 +82,35 @@ struct RelaxPotential : public Component { void finally(const Options& state) override; void outputVars(Options& state) override; + + // // Save and restore potential phi + // void restartVars(Options& state) override { + // AUTO_TRACE(); + + // // NOTE: This is a hack because we know that the loaded restart file + // // is passed into restartVars in PhysicsModel::postInit + // // The restart value should be used in init() rather than here + // static bool first = true; + // if (first and state.isSet("phi")) { + // first = false; + // phi = state["phi"].as(); + // } + + // // Save the potential + // set_with_attrs(state["phi"], phi, + // {{"long_name", "plasma potential"}, + // {"source", "vorticity"}}); + // } private: Field3D Vort; // Evolving vorticity Field3D phi1; // Scaled electrostatic potential, evolving in time ϕ_1 = λ_2 ϕ Field3D phi; // Electrostatic potential + Field3D Pi_hat; ///< Contribution from ion pressure, weighted by atomic mass / charge + bool exb_advection; //< Include nonlinear ExB advection? + bool exb_advection_simplified; // Simplify nonlinear ExB advection form? bool diamagnetic; //< Include diamagnetic current? bool diamagnetic_polarisation; //< Include diamagnetic drift in polarisation current? bool boussinesq; ///< Use the Boussinesq approximation? @@ -70,12 +119,37 @@ private: bool poloidal_flows; ///< Include poloidal ExB flow? bool bndry_flux; ///< Allow flows through radial boundaries? + bool collisional_friction; ///< Damping of vorticity due to collisional friction + bool sheath_boundary; ///< Set outer boundary to j=0? + BoutReal Ge; // Secondary electron emission coefficient + BoutReal sin_alpha; // sin of angle between magnetic field and wall. + Field3D wall_potential; ///< Voltage at the wall. Normalised units. + + 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 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? + Field2D Bsq; ///< SQ(coord->Bxy) Vector2D Curlb_B; ///< Curvature vector Curl(b/B) + Field2D viscosity; ///< Perpendicular Kinematic viscosity + Field2D viscosity_par; ///< Parallel Kinematic viscosity + + // Relax-potential related variables + BoutReal lambda_1; ///< Relaxation parameters. NOTE: lambda_1 has dimentions! + BoutReal lambda_2; ///< Relaxation parameters + + // Diagnostic outputs + Field3D DivJdia, DivJcol; // Divergence of diamagnetic and collisional current - BoutReal lambda_1, lambda_2; ///< Relaxation parameters + bool diagnose; ///< Output additional diagnostics? }; namespace { diff --git a/src/relax_potential.cxx b/src/relax_potential.cxx index aa7bbb061..266b70af9 100644 --- a/src/relax_potential.cxx +++ b/src/relax_potential.cxx @@ -1,15 +1,59 @@ -#include #include -using bout::globals::mesh; +#include "../include/relax_potential.hxx" #include "../include/div_ops.hxx" -#include "../include/relax_potential.hxx" +#include + +#include +#include "../include/hermes_build_config.hxx" +#include "hermes_utils.hxx" + +using bout::globals::mesh; + + +namespace { +/// Limited free gradient of log of a quantity +/// This ensures that the guard cell values remain positive +/// while also ensuring that the quantity never increases +/// +/// fm fc | fp +/// ^ boundary +/// +/// exp( 2*log(fc) - log(fm) ) +/// +BoutReal limitFree(BoutReal fm, BoutReal fc) { + if (fm < fc) { + return fc; // Neumann rather than increasing into boundary + } + if (fm < 1e-10) { + return fc; // Low / no density condition + } + BoutReal fp = SQ(fc) / fm; +#if CHECKLEVEL >= 2 + if (!std::isfinite(fp)) { + throw BoutException("SheathBoundary limitFree: {}, {} -> {}", fm, fc, fp); + } +#endif + + return fp; +} +} // namespace + RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* solver) { AUTO_TRACE(); - auto* coord = mesh->getCoordinates(); + solver->add(Vort, "Vort"); // Vorticity evolving + solver->add(phi1, "phi1"); // Evolving scaled potential ϕ_1 = λ_2 ϕ + + const Options& units = alloptions["units"]; + // Normalisations + const BoutReal Omega_ci = 1. / units["seconds"].as(); + const BoutReal Bnorm = units["Tesla"]; + const BoutReal Tnorm = units["eV"]; + const BoutReal Nnorm = units["inv_meters_cubed"]; + const BoutReal Lnorm = units["meters"]; auto& options = alloptions[name]; @@ -17,14 +61,57 @@ RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* so .doc("Include nonlinear ExB advection?") .withDefault(true); + exb_advection_simplified = options["exb_advection_simplified"] + .doc("Simplify nonlinear ExB advection form?") + .withDefault(true); + diamagnetic = options["diamagnetic"].doc("Include diamagnetic current?").withDefault(true); + sheath_boundary = options["sheath_boundary"] + .doc("Set potential to j=0 sheath at radial boundaries? (default = 0)") + .withDefault(false); + + Ge = options["secondary_electron_coef"] + .doc("Effective secondary electron emission coefficient") + .withDefault(0.0); + + if ((Ge < 0.0) or (Ge > 1.0)) { + throw BoutException("Secondary electron emission must be between 0 and 1 ({:e})", Ge); + } + + sin_alpha = options["sin_alpha"] + .doc("Sin of the angle between magnetic field line and wall surface. " + "Should be between 0 and 1") + .withDefault(1.0); + + if ((sin_alpha < 0.0) or (sin_alpha > 1.0)) { + throw BoutException("Range of sin_alpha must be between 0 and 1"); + } + + // Read wall voltage, convert to normalised units + wall_potential = options["wall_potential"] + .doc("Voltage of the wall [Volts]") + .withDefault(Field3D(0.0)) + / Tnorm; + // Convert to field aligned coordinates + wall_potential = toFieldAligned(wall_potential); + + // Note: wall potential at the last cell before the boundary is used, + // not the value at the boundary half-way between cells. This is due + // to how twist-shift boundary conditions and non-aligned inputs are + // treated; using the cell boundary gives incorrect results. + diamagnetic_polarisation = options["diamagnetic_polarisation"] .doc("Include diamagnetic drift in polarisation current?") .withDefault(true); + collisional_friction = + options["collisional_friction"] + .doc("Damp vorticity based on mass-weighted collision frequency") + .withDefault(false); + boussinesq = options["boussinesq"] .doc("Use the Boussinesq approximation?") .withDefault(true); @@ -34,59 +121,384 @@ RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* so "(Boussinesq approximation)") .withDefault(2.0); // Deuterium + bndry_flux = options["bndry_flux"] + .doc("Allow flows through radial boundaries") + .withDefault(true); + poloidal_flows = options["poloidal_flows"].doc("Include poloidal ExB flow").withDefault(true); - lambda_1 = options["lambda_1"].doc("λ_1 > λ_2 > 1").withDefault(1.0); - lambda_2 = options["lambda_2"].doc("λ_2 > 1").withDefault(10.0); + viscosity = 0.0; + viscosity = options["viscosity"] + .doc("Perpendicular Kinematic viscosity [m^2/s]") + .withDefault(viscosity) + / (Lnorm * Lnorm * Omega_ci); + viscosity.applyBoundary("dirichlet"); - solver->add(Vort, "Vort"); // Vorticity evolving - solver->add(phi1, "phi1"); // Evolving scaled potential ϕ_1 = λ_2 ϕ + viscosity_par = 0.0; + viscosity_par = options["viscosity_par"] + .doc("Parallel Kinematic viscosity [m^2/s]") + .withDefault(viscosity_par) + / (Lnorm * Lnorm * Omega_ci); + viscosity_par.applyBoundary("dirichlet"); + viscosity_par.applyParallelBoundary("parallel_dirichlet_o2"); - if (diamagnetic) { - // Read curvature vector - try { - Curlb_B.covariant = false; // Contravariant - mesh->get(Curlb_B, "bxcv"); + vort_dissipation = options["vort_dissipation"] + .doc("Parallel dissipation of vorticity") + .withDefault(false); - } catch (BoutException& e) { + phi_dissipation = options["phi_dissipation"] + .doc("Parallel dissipation of potential [Recommended]") + .withDefault(true); + + phi_boundary_relax = options["phi_boundary_relax"] + .doc("Relax x boundaries of phi towards Neumann?") + .withDefault(false); + + phi_sheath_dissipation = options["phi_sheath_dissipation"] + .doc("Add dissipation when phi < 0.0 at the sheath") + .withDefault(false); + + damp_core_vorticity = options["damp_core_vorticity"] + .doc("Damp vorticity at the core boundary?") + .withDefault(false); + + lambda_1 = options["lambda_1"].doc("λ_1 > 1").withDefault(3e8) / (Tnorm * Omega_ci / SI::qe / Nnorm); + lambda_2 = options["lambda_2"].doc("λ_2 > λ_1").withDefault(1.0); + + + // NOTE(malamast): How do we do that? + // Add phi to restart files so that the value in the boundaries + // is restored on restart. This is done even when phi is not evolving, + // so that phi can be saved and re-loaded + + // Set initial value. Will be overwritten if restarting + phi1 = 0.0; + Vort = 0.0; + + auto* coord = mesh->getCoordinates(); + + if (phi_boundary_relax) { + // Set the last update time to -1, so it will reset + // the first time RHS function is called + phi_boundary_last_update = -1.; + + phi_core_averagey = options["phi_core_averagey"] + .doc("Average phi core boundary in Y?") + .withDefault(false) and mesh->periodicY(mesh->xstart); + + phi_boundary_timescale = options["phi_boundary_timescale"] + .doc("Timescale for phi boundary relaxation [seconds]") + .withDefault(1e-4) + / get(alloptions["units"]["seconds"]); + } + + + // Read curvature vector + try { + Curlb_B.covariant = false; // Contravariant + mesh->get(Curlb_B, "bxcv"); + + } catch (BoutException& e) { + try { // May be 2D, reading as 3D Vector2D curv2d; curv2d.covariant = false; mesh->get(curv2d, "bxcv"); Curlb_B = curv2d; + } catch (BoutException& e) { + if (diamagnetic) { + // Need curvature + throw; + } else { + output_warn.write("No curvature vector in input grid"); + Curlb_B = 0.0; + } } + } - if (Options::root()["mesh"]["paralleltransform"]["type"].as() - == "shifted") { - Field2D I; - mesh->get(I, "sinty"); - Curlb_B.z += I * Curlb_B.x; - } - - Options& units = alloptions["units"]; - BoutReal Bnorm = units["Tesla"]; - BoutReal Lnorm = units["meters"]; + if (Options::root()["mesh"]["paralleltransform"]["type"].as() + == "shifted") { + Field2D I; + mesh->get(I, "sinty"); + Curlb_B.z += I * Curlb_B.x; + } - Curlb_B.x /= Bnorm; - Curlb_B.y *= SQ(Lnorm); - Curlb_B.z *= SQ(Lnorm); + Curlb_B.x /= Bnorm; + Curlb_B.y *= SQ(Lnorm); + Curlb_B.z *= SQ(Lnorm); - Curlb_B *= 2. / coord->Bxy; - } + Curlb_B *= 2. / coord->Bxy; Bsq = SQ(coord->Bxy); + + diagnose = options["diagnose"] + .doc("Output additional diagnostics?") + .withDefault(false); } void RelaxPotential::transform(Options& state) { AUTO_TRACE(); + Options& allspecies = state["species"]; + + // phi.name = "phi"; + auto& fields = state["fields"]; + + // Set the boundary of phi1. + phi1.applyBoundary("neumann"); + // Scale potential phi = phi1 / lambda_2; - phi.applyBoundary("neumann"); + + // 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(); + Field3D &Vort_yup = Vort.yup(); + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + // Vort_ydown(r.ind, mesh->ystart - 1, jz) = 2 * Vort(r.ind, mesh->ystart, jz) - Vort_yup(r.ind, mesh->ystart + 1, jz); + Vort_ydown(r.ind, mesh->ystart - 1, jz) = Vort(r.ind, mesh->ystart, jz); // TODO: iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + // Vort_yup(r.ind, mesh->yend + 1, jz) = 2 * Vort(r.ind, mesh->yend, jz) - Vort_ydown(r.ind, mesh->yend - 1, jz); + Vort_yup(r.ind, mesh->yend + 1, jz) = Vort(r.ind, mesh->yend, jz); + } + } + } else { + Field3D Vort_fa = toFieldAligned(Vort); + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + // Vort_fa(r.ind, mesh->ystart - 1, jz) = 2 * Vort_fa(r.ind, mesh->ystart, jz) - Vort_fa(r.ind, mesh->ystart + 1, jz); + Vort_fa(r.ind, mesh->ystart - 1, jz) = Vort_fa(r.ind, mesh->ystart, jz); // Neumann BC + } + } + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + // Vort_fa(r.ind, mesh->yend + 1, jz) = 2 * Vort_fa(r.ind, mesh->yend, jz) - Vort_fa(r.ind, mesh->yend - 1, jz); + Vort_fa(r.ind, mesh->yend + 1, jz) = Vort_fa(r.ind, mesh->yend, jz); // Neumann BC + } + } + Vort = fromFieldAligned(Vort_fa); + } - auto& fields = state["fields"]; + // 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 + if (diamagnetic_polarisation) { + // Diamagnetic term in vorticity. Note this is weighted by the mass + // This includes all species, including electrons + // Options& allspecies = state["species"]; + for (auto& kv : allspecies.getChildren()) { + Options& species = allspecies[kv.first]; // Note: need non-const + + if (!(IS_SET_NOBOUNDARY(species["pressure"]) and species.isSet("charge") + and species.isSet("AA"))) { + continue; // No pressure, charge or mass -> no polarisation current + } + + const auto charge = get(species["charge"]); + if (fabs(charge) < 1e-5) { + // No charge + continue; + } + + // Don't need sheath boundary + const auto P = GET_NOBOUNDARY(Field3D, species["pressure"]); + const auto AA = get(species["AA"]); + + Pi_hat += P * (AA / average_atomic_mass / charge); + } + } + + Pi_hat.applyBoundary("neumann"); + + if (phi_boundary_relax) { + // Update the boundary regions by relaxing towards zero gradient + // on a given timescale. + + BoutReal time = get(state["time"]); + + if (phi_boundary_last_update < 0.0) { + // First time this has been called. + phi_boundary_last_update = time; + + } else if (time > phi_boundary_last_update) { + // Only update if time has advanced + // Uses an exponential decay of the weighting of the value in the boundary + // so that the solution is well behaved for arbitrary steps + BoutReal weight = exp(-(time - phi_boundary_last_update) / phi_boundary_timescale); + phi_boundary_last_update = time; + + if (mesh->firstX()) { + BoutReal phivalue = 0.0; + if (phi_core_averagey) { + BoutReal philocal = 0.0; + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + philocal += phi(mesh->xstart, j, k); + } + } + MPI_Comm comm_inner = mesh->getYcomm(0); + int np; + MPI_Comm_size(comm_inner, &np); + MPI_Allreduce(&philocal, + &phivalue, + 1, MPI_DOUBLE, + MPI_SUM, comm_inner); + phivalue /= (np * mesh->LocalNz * mesh->LocalNy); + } + + for (int j = mesh->ystart; j <= mesh->yend; j++) { + if (!phi_core_averagey) { + phivalue = 0.0; // Calculate phi boundary for each Y index separately + for (int k = 0; k < mesh->LocalNz; k++) { + phivalue += phi(mesh->xstart, j, k); + } + phivalue /= mesh->LocalNz; // Average in Z of point next to boundary + } + + // Old value of phi at boundary + BoutReal oldvalue = + 0.5 * (phi(mesh->xstart - 1, j, 0) + phi(mesh->xstart, j, 0)); + + // New value of phi at boundary, relaxing towards phivalue + BoutReal newvalue = weight * oldvalue + (1. - weight) * phivalue; + + // Set phi at the boundary to this value + for (int k = 0; k < mesh->LocalNz; k++) { + phi(mesh->xstart - 1, j, k) = 2. * newvalue - phi(mesh->xstart, j, k); + + // Note: This seems to make a difference, but don't know why. + // Without this, get convergence failures with no apparent instability + // (all fields apparently smooth, well behaved) + phi(mesh->xstart - 2, j, k) = phi(mesh->xstart - 1, j, k); + } + } + } + + if (mesh->lastX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + BoutReal phivalue = 0.0; + for (int k = 0; k < mesh->LocalNz; k++) { + phivalue += phi(mesh->xend, j, k); + } + phivalue /= mesh->LocalNz; // Average in Z of point next to boundary + + // Old value of phi at boundary + BoutReal oldvalue = 0.5 * (phi(mesh->xend + 1, j, 0) + phi(mesh->xend, j, 0)); + + // New value of phi at boundary, relaxing towards phivalue + BoutReal newvalue = weight * oldvalue + (1. - weight) * phivalue; + + // Set phi at the boundary to this value + for (int k = 0; k < mesh->LocalNz; k++) { + phi(mesh->xend + 1, j, k) = 2. * newvalue - phi(mesh->xend, j, k); + + // Note: This seems to make a difference, but don't know why. + // Without this, get convergence failures with no apparent instability + // (all fields apparently smooth, well behaved) + phi(mesh->xend + 2, j, k) = phi(mesh->xend + 1, j, k); + } + } + } + } + } else { + // phi_boundary_relax = false + // + // Set boundary from temperature, to be consistent with j=0 at sheath + + // Sheath multiplier Te -> phi (2.84522 for Deuterium) + BoutReal sheathmult = 0.0; + if (sheath_boundary) { + BoutReal Me_Mp = get(state["species"]["e"]["AA"]); + sheathmult = log(0.5 * sqrt(1. / (Me_Mp * PI))); + } + + Field3D Te; // Electron temperature, use for outer boundary conditions + if (state["species"]["e"].isSet("temperature")) { + // Electron temperature set + Te = GET_NOBOUNDARY(Field3D, state["species"]["e"]["temperature"]); + } else { + Te = 0.0; + } + + // Sheath multiplier Te -> phi (2.84522 for Deuterium if Ti = 0) + if (mesh->firstX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + BoutReal teavg = 0.0; // Average Te in Z + for (int k = 0; k < mesh->LocalNz; k++) { + teavg += Te(mesh->xstart, j, k); + } + teavg /= mesh->LocalNz; + BoutReal phivalue = sheathmult * teavg; + // Set midpoint (boundary) value + for (int k = 0; k < mesh->LocalNz; k++) { + phi(mesh->xstart - 1, j, k) = 2. * phivalue - phi(mesh->xstart, j, k); + + // Note: This seems to make a difference, but don't know why. + // Without this, get convergence failures with no apparent instability + // (all fields apparently smooth, well behaved) + phi(mesh->xstart - 2, j, k) = phi(mesh->xstart - 1, j, k); + } + } + } + + if (mesh->lastX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + BoutReal teavg = 0.0; // Average Te in Z + for (int k = 0; k < mesh->LocalNz; k++) { + teavg += Te(mesh->xend, j, k); + } + teavg /= mesh->LocalNz; + BoutReal phivalue = sheathmult * teavg; + // Set midpoint (boundary) value + for (int k = 0; k < mesh->LocalNz; k++) { + phi(mesh->xend + 1, j, k) = 2. * phivalue - phi(mesh->xend, j, k); + + // Note: This seems to make a difference, but don't know why. + // Without this, get convergence failures with no apparent instability + // (all fields apparently smooth, well behaved) + phi(mesh->xend + 2, j, k) = phi(mesh->xend + 1, j, k); + } + } + } + } + + + // Outer boundary cells + if (mesh->firstX()) { + for (int i = mesh->xstart - 2; i >= 0; --i) { + for (int j = mesh->ystart; j <= mesh->yend; ++j) { + for (int k = 0; k < mesh->LocalNz; ++k) { + phi(i, j, k) = phi(i + 1, j, k); + } + } + } + } + if (mesh->lastX()) { + for (int i = mesh->xend + 2; i < mesh->LocalNx; ++i) { + for (int j = mesh->ystart; j <= mesh->yend; ++j) { + for (int k = 0; k < mesh->LocalNz; ++k) { + phi(i, j, k) = phi(i - 1, j, k); + } + } + } + } + + + // Ensure that potential is set in the communication guard cells + mesh->communicate(phi, phi1, Vort); //NOTE(malamast): Should we include phi1? + + // Vorticity equation ddt(Vort) = 0.0; @@ -100,18 +512,52 @@ void RelaxPotential::transform(Options& state) { Jdia.z = 0.0; Jdia.covariant = Curlb_B.covariant; - Options& allspecies = state["species"]; - // Pre-calculate this rather than calculate for each species + // Note: The below calculation requires phi derivatives at the Y boundaries + // NOTE: Do we actually need that? + // Setting to free boundaries + + if (phi.hasParallelSlices()) { + Field3D &phi_ydown = phi.ydown(); + Field3D &phi_yup = phi.yup(); + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + phi_ydown(r.ind, mesh->ystart - 1, jz) = 2 * phi(r.ind, mesh->ystart, jz) - phi_yup(r.ind, mesh->ystart + 1, jz); + } + } + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + phi_yup(r.ind, mesh->yend + 1, jz) = 2 * phi(r.ind, mesh->yend, jz) - phi_ydown(r.ind, mesh->yend - 1, jz); + } + } + } else { + Field3D phi_fa = toFieldAligned(phi); + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + phi_fa(r.ind, mesh->ystart - 1, jz) = 2 * phi_fa(r.ind, mesh->ystart, jz) - phi_fa(r.ind, mesh->ystart + 1, jz); + } + } + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + phi_fa(r.ind, mesh->yend + 1, jz) = 2 * phi_fa(r.ind, mesh->yend, jz) - phi_fa(r.ind, mesh->yend - 1, jz); + } + } + phi = fromFieldAligned(phi_fa); + } + Vector3D Grad_phi = Grad(phi); for (auto& kv : allspecies.getChildren()) { Options& species = allspecies[kv.first]; // Note: need non-const - if (!(IS_SET_NOBOUNDARY(species["pressure"]) and IS_SET(species["charge"]) - and (get(species["charge"]) != 0.0))) { + if (!(IS_SET_NOBOUNDARY(species["pressure"]) and IS_SET(species["charge"]))) { continue; // No pressure or charge -> no diamagnetic current } + if (fabs(get(species["charge"])) < 1e-5) { + // No charge + continue; + } + // Note that the species must have a charge, but charge is not used, // because it cancels out in the expression for current @@ -128,27 +574,48 @@ void RelaxPotential::transform(Options& state) { // Note: This term is central differencing so that it balances // the corresponding compression term in the species pressure equations - Field3D DivJdia = Div(Jdia); + DivJdia = Div(Jdia); ddt(Vort) += DivJdia; - if (diamagnetic_polarisation) { - // Calculate energy exchange term nonlinear in pressure - // ddt(Pi) += Pi * Div((Pe + Pi) * Curlb_B); - for (auto& kv : allspecies.getChildren()) { - Options& species = allspecies[kv.first]; // Note: need non-const + set(fields["DivJdia"], DivJdia); + } - if (!(IS_SET_NOBOUNDARY(species["pressure"]) and IS_SET(species["charge"]) - and IS_SET(species["AA"]))) { - continue; // No pressure, charge or mass -> no polarisation current due to - // rate of change of diamagnetic flow - } - auto P = GET_NOBOUNDARY(Field3D, species["pressure"]); + if (collisional_friction) { + // Damping of vorticity due to collisions + + // Calculate a mass-weighted collision frequency + Field3D sum_A_nu_n = + zeroFrom(Vort); // Sum of atomic mass * collision frequency * density + Field3D sum_A_n = zeroFrom(Vort); // Sum of atomic mass * density + + // const Options& allspecies = state["species"]; + for (const auto& kv : allspecies.getChildren()) { + const Options& species = kv.second; + + if (!(species.isSet("charge") and species.isSet("AA"))) { + continue; // No charge or mass -> no current + } + if (fabs(get(species["charge"])) < 1e-5) { + continue; // Zero charge + } - add(species["energy_source"], (3. / 2) * P * DivJdia); + const BoutReal A = get(species["AA"]); + const Field3D N = GET_NOBOUNDARY(Field3D, species["density"]); + const Field3D AN = A * N; + sum_A_n += AN; + if (IS_SET(species["collision_frequency"])) { + sum_A_nu_n += AN * GET_VALUE(Field3D, species["collision_frequency"]); } } - set(fields["DivJdia"], DivJdia); + Field3D weighted_collision_frequency = sum_A_nu_n / sum_A_n; + weighted_collision_frequency.applyBoundary("neumann"); + + DivJcol = -FV::Div_a_Grad_perp( + weighted_collision_frequency * average_atomic_mass / Bsq, phi + Pi_hat); + + ddt(Vort) += DivJcol; + set(fields["DivJcol"], DivJcol); } set(fields["vorticity"], Vort); @@ -160,16 +627,48 @@ void RelaxPotential::finally(const Options& state) { const Options& allspecies = state["species"]; + auto* coord = mesh->getCoordinates(); + phi = get(state["fields"]["phi"]); Vort = get(state["fields"]["vorticity"]); + // Solve vorticity equation + if (exb_advection) { - ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(Vort, phi, bndry_flux, poloidal_flows); + // These terms come from divergence of polarisation current + + if (exb_advection_simplified) { + // By default this is a simplified nonlinear term + ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(Vort, phi, bndry_flux, poloidal_flows); + + } else { + // If diamagnetic_polarisation = false and B is constant, then + // this term reduces to the simplified form above. + // + // Because this is implemented in terms of an operation on the result + // of an operation, we need to communicate and the resulting stencil is + // wider than the simple form. + ddt(Vort) -= + Div_n_bxGrad_f_B_XPPM(0.5 * Vort, phi, bndry_flux, poloidal_flows); + + // V_ExB dot Grad(Pi) + Field3D vEdotGradPi = bracket(phi, Pi_hat, BRACKET_ARAKAWA); + vEdotGradPi.applyBoundary("free_o2"); + + // delp2(phi) term + Field3D DelpPhi_2B2 = 0.5 * average_atomic_mass * Delp2(phi) / Bsq; + DelpPhi_2B2.applyBoundary("free_o2"); + + mesh->communicate(vEdotGradPi, DelpPhi_2B2); + + ddt(Vort) -= FV::Div_a_Grad_perp(0.5 * average_atomic_mass / Bsq, vEdotGradPi); + ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(DelpPhi_2B2, phi + Pi_hat, bndry_flux, + poloidal_flows); + } } if (state.isSection("fields") and state["fields"].isSet("DivJextra")) { auto DivJextra = get(state["fields"]["DivJextra"]); - // Parallel current is handled here, to allow different 2D or 3D closures // to be used ddt(Vort) += DivJextra; @@ -192,9 +691,89 @@ void RelaxPotential::finally(const Options& state) { const BoutReal A = get(species["AA"]); // Note: Using NV rather than N*V so that the cell boundary flux is correct - ddt(Vort) += Div_par((Z / A) * NV); + const Field3D jpar = (Z / A) * NV; + ddt(Vort) += Div_par(jpar); + + if (state["fields"].isSet("Apar_flutter")) { + // Magnetic flutter term + const Field3D Apar_flutter = get(state["fields"]["Apar_flutter"]); + + // Div_par(jpar) = B * Grad_par(jpar / B) + // Using the approximation for small delta-B/B + // b dot Grad(jpar) = Grad_par(jpar) + [jpar, Apar] + ddt(Vort) += coord->Bxy * bracket(jpar / coord->Bxy, Apar_flutter, BRACKET_ARAKAWA); + } + } + + // Viscosity + ddt(Vort) += FV::Div_a_Grad_perp(viscosity, Vort); + + ddt(Vort) += FV::Div_par_K_Grad_par(viscosity_par, Vort); + // NOTE(malamast): Need to check if Div_par_K_Grad_par is equivalent with the Laplace_par operator. + // ddt(Vort) += viscosity_par * Laplace_par(Vort); + + + if (vort_dissipation) { + // Adds dissipation term like in other equations + Field3D sound_speed = get(state["sound_speed"]); + ddt(Vort) -= FV::Div_par(Vort, 0.0, sound_speed); + } + + if (phi_dissipation) { + // Adds dissipation term like in other equations, but depending on gradient of + // potential + Field3D sound_speed = get(state["sound_speed"]); + + Field3D zero {0.0}; + zero.splitParallelSlices(); + zero.yup() = 0.0; + zero.ydown() = 0.0; + + Field3D dummy; + ddt(Vort) -= FV::Div_par_mod(-phi, zero, sound_speed, dummy); + } + + if (phi_sheath_dissipation) { + // Dissipation when phi < 0.0 at the sheath + + auto phi_fa = toFieldAligned(phi); + Field3D dissipation{zeroFrom(phi_fa)}; + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + auto i = indexAt(phi_fa, r.ind, mesh->ystart, jz); + BoutReal phisheath = 0.5*(phi_fa[i] + phi_fa[i.ym()]); + dissipation[i] = -floor(-phisheath, 0.0); + } + } + + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + auto i = indexAt(phi_fa, r.ind, mesh->yend, jz); + BoutReal phisheath = 0.5*(phi_fa[i] + phi_fa[i.yp()]); + dissipation[i] = -floor(-phisheath, 0.0); + } + } + ddt(Vort) += fromFieldAligned(dissipation); + } + + if (damp_core_vorticity) { + // Damp axisymmetric vorticity near core boundary + if (mesh->firstX() and mesh->periodicY(mesh->xstart)) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + BoutReal vort_avg = 0.0; // Average Vort in Z + for (int k = 0; k < mesh->LocalNz; k++) { + vort_avg += Vort(mesh->xstart, j, k); + } + vort_avg /= mesh->LocalNz; + for (int k = 0; k < mesh->LocalNz; k++) { + ddt(Vort)(mesh->xstart, j, k) -= 0.01 * vort_avg; + } + } + } } + + // Solve diffusion equation for potential if (boussinesq) { @@ -245,7 +824,7 @@ void RelaxPotential::finally(const Options& state) { 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, Pi); + phi_vort += FV::Div_a_Grad_perp(Ai / Bsq / Zi, Pi); } } @@ -256,7 +835,15 @@ void RelaxPotential::finally(const Options& state) { void RelaxPotential::outputVars(Options& state) { AUTO_TRACE(); // Normalisations - auto Tnorm = state["Tnorm"].as(); + auto Nnorm = get(state["Nnorm"]); + auto Tnorm = get(state["Tnorm"]); + auto Omega_ci = get(state["Omega_ci"]); + + state["Vort"].setAttributes({{"time_dimension", "t"}, + {"units", "C m^-3"}, + {"conversion", SI::qe * Nnorm}, + {"long_name", "vorticity"}, + {"source", "relax_potential"}}); set_with_attrs(state["phi"], phi, {{"time_dimension", "t"}, @@ -265,4 +852,35 @@ void RelaxPotential::outputVars(Options& state) { {"standard_name", "potential"}, {"long_name", "plasma potential"}, {"source", "relax_potential"}}); + + if (diagnose) { + set_with_attrs(state["ddt(Vort)"], ddt(Vort), + {{"time_dimension", "t"}, + {"units", "A m^-3"}, + {"conversion", SI::qe * Nnorm * Omega_ci}, + {"long_name", "Rate of change of vorticity"}, + {"source", "relax_potential"}}); + set_with_attrs(state["ddt(phi)"], ddt(phi1), + {{"time_dimension", "t"}, + {"units", "V/s"}, + {"conversion", Tnorm * Omega_ci}, + {"long_name", "Rate of change of electrostatic potential"}, + {"source", "relax_potential"}}); + if (diamagnetic) { + set_with_attrs(state["DivJdia"], DivJdia, + {{"time_dimension", "t"}, + {"units", "A m^-3"}, + {"conversion", SI::qe * Nnorm * Omega_ci}, + {"long_name", "Divergence of diamagnetic current"}, + {"source", "relax_potential"}}); + } + if (collisional_friction) { + set_with_attrs(state["DivJcol"], DivJcol, + {{"time_dimension", "t"}, + {"units", "A m^-3"}, + {"conversion", SI::qe * Nnorm * Omega_ci}, + {"long_name", "Divergence of collisional current"}, + {"source", "relax_potential"}}); + } + } } From 14e89bf69691c5b275d13c37771a8761e7c4b735 Mon Sep 17 00:00:00 2001 From: malamast Date: Thu, 8 Jan 2026 16:17:47 -0800 Subject: [PATCH 3/4] relax_potential: tidied up the code. Added parallel BC for presure as it was done in Vorticity component --- include/relax_potential.hxx | 5 +- src/relax_potential.cxx | 124 ++++++++++++++++++------------------ src/vorticity.cxx | 2 +- 3 files changed, 63 insertions(+), 68 deletions(-) diff --git a/include/relax_potential.hxx b/include/relax_potential.hxx index f35299ea5..e1341e384 100644 --- a/include/relax_potential.hxx +++ b/include/relax_potential.hxx @@ -123,10 +123,6 @@ private: bool sheath_boundary; ///< Set outer boundary to j=0? - BoutReal Ge; // Secondary electron emission coefficient - BoutReal sin_alpha; // sin of angle between magnetic field and wall. - Field3D wall_potential; ///< Voltage at the wall. Normalised units. - 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 @@ -139,6 +135,7 @@ private: 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 diff --git a/src/relax_potential.cxx b/src/relax_potential.cxx index 266b70af9..f61988f31 100644 --- a/src/relax_potential.cxx +++ b/src/relax_potential.cxx @@ -1,17 +1,15 @@ -#include #include "../include/relax_potential.hxx" - #include "../include/div_ops.hxx" -#include +#include "../include/hermes_utils.hxx" #include -#include "../include/hermes_build_config.hxx" -#include "hermes_utils.hxx" +#include +#include +#include using bout::globals::mesh; - namespace { /// Limited free gradient of log of a quantity /// This ensures that the guard cell values remain positive @@ -40,22 +38,20 @@ BoutReal limitFree(BoutReal fm, BoutReal fc) { } } // namespace - RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* solver) { AUTO_TRACE(); solver->add(Vort, "Vort"); // Vorticity evolving solver->add(phi1, "phi1"); // Evolving scaled potential ϕ_1 = λ_2 ϕ - const Options& units = alloptions["units"]; + auto& options = alloptions[name]; // Normalisations + const Options& units = alloptions["units"]; const BoutReal Omega_ci = 1. / units["seconds"].as(); const BoutReal Bnorm = units["Tesla"]; + const BoutReal Lnorm = units["meters"]; const BoutReal Tnorm = units["eV"]; const BoutReal Nnorm = units["inv_meters_cubed"]; - const BoutReal Lnorm = units["meters"]; - - auto& options = alloptions[name]; exb_advection = options["exb_advection"] .doc("Include nonlinear ExB advection?") @@ -72,36 +68,6 @@ RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* so .doc("Set potential to j=0 sheath at radial boundaries? (default = 0)") .withDefault(false); - Ge = options["secondary_electron_coef"] - .doc("Effective secondary electron emission coefficient") - .withDefault(0.0); - - if ((Ge < 0.0) or (Ge > 1.0)) { - throw BoutException("Secondary electron emission must be between 0 and 1 ({:e})", Ge); - } - - sin_alpha = options["sin_alpha"] - .doc("Sin of the angle between magnetic field line and wall surface. " - "Should be between 0 and 1") - .withDefault(1.0); - - if ((sin_alpha < 0.0) or (sin_alpha > 1.0)) { - throw BoutException("Range of sin_alpha must be between 0 and 1"); - } - - // Read wall voltage, convert to normalised units - wall_potential = options["wall_potential"] - .doc("Voltage of the wall [Volts]") - .withDefault(Field3D(0.0)) - / Tnorm; - // Convert to field aligned coordinates - wall_potential = toFieldAligned(wall_potential); - - // Note: wall potential at the last cell before the boundary is used, - // not the value at the boundary half-way between cells. This is due - // to how twist-shift boundary conditions and non-aligned inputs are - // treated; using the cell boundary gives incorrect results. - diamagnetic_polarisation = options["diamagnetic_polarisation"] .doc("Include diamagnetic drift in polarisation current?") @@ -117,7 +83,7 @@ RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* so .withDefault(true); average_atomic_mass = options["average_atomic_mass"] - .doc("Weighted average atomic mass, for polarisaion current " + .doc("Weighted average atomic mass, for polarisation current " "(Boussinesq approximation)") .withDefault(2.0); // Deuterium @@ -143,6 +109,13 @@ RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* so viscosity_par.applyBoundary("dirichlet"); viscosity_par.applyParallelBoundary("parallel_dirichlet_o2"); + hyper_z = options["hyper_z"].doc("Hyper-viscosity in Z. < 0 -> off").withDefault(-1.0); + + // Numerical dissipation terms + // These are required to suppress parallel zig-zags in + // cell centred formulations. Essentially adds (hopefully small) + // parallel currents + vort_dissipation = options["vort_dissipation"] .doc("Parallel dissipation of vorticity") .withDefault(false); @@ -167,7 +140,6 @@ RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* so lambda_2 = options["lambda_2"].doc("λ_2 > λ_1").withDefault(1.0); - // NOTE(malamast): How do we do that? // Add phi to restart files so that the value in the boundaries // is restored on restart. This is done even when phi is not evolving, // so that phi can be saved and re-loaded @@ -176,7 +148,7 @@ RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* so phi1 = 0.0; Vort = 0.0; - auto* coord = mesh->getCoordinates(); + auto coord = mesh->getCoordinates(); if (phi_boundary_relax) { // Set the last update time to -1, so it will reset @@ -193,7 +165,6 @@ RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* so / get(alloptions["units"]["seconds"]); } - // Read curvature vector try { Curlb_B.covariant = false; // Contravariant @@ -242,7 +213,7 @@ void RelaxPotential::transform(Options& state) { Options& allspecies = state["species"]; - // phi.name = "phi"; + phi.name = "phi"; auto& fields = state["fields"]; // Set the boundary of phi1. @@ -514,9 +485,7 @@ void RelaxPotential::transform(Options& state) { // Pre-calculate this rather than calculate for each species // Note: The below calculation requires phi derivatives at the Y boundaries - // NOTE: Do we actually need that? // Setting to free boundaries - if (phi.hasParallelSlices()) { Field3D &phi_ydown = phi.ydown(); Field3D &phi_yup = phi.yup(); @@ -563,6 +532,38 @@ void RelaxPotential::transform(Options& state) { auto P = GET_NOBOUNDARY(Field3D, species["pressure"]); + // Note: We need boundary conditions on P, so apply the same + // free boundary condition as sheath_boundary. + if (P.hasParallelSlices()) { + Field3D &P_ydown = P.ydown(); + Field3D &P_yup = P.yup(); + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + P_ydown(r.ind, mesh->ystart - 1, jz) = 2 * P(r.ind, mesh->ystart, jz) - P_yup(r.ind, mesh->ystart + 1, jz); + } + } + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + P_yup(r.ind, mesh->yend + 1, jz) = 2 * P(r.ind, mesh->yend, jz) - P_ydown(r.ind, mesh->yend - 1, jz); + } + } + } else { + Field3D P_fa = toFieldAligned(P); + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + auto i = indexAt(P_fa, r.ind, mesh->ystart, jz); + P_fa[i.ym()] = limitFree(P_fa[i.yp()], P_fa[i]); + } + } + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + auto i = indexAt(P_fa, r.ind, mesh->yend, jz); + P_fa[i.yp()] = limitFree(P_fa[i.ym()], P_fa[i]); + } + } + P = fromFieldAligned(P_fa); + } + Vector3D Jdia_species = P * Curlb_B; // Diamagnetic current for this species // This term energetically balances diamagnetic term @@ -588,7 +589,6 @@ void RelaxPotential::transform(Options& state) { zeroFrom(Vort); // Sum of atomic mass * collision frequency * density Field3D sum_A_n = zeroFrom(Vort); // Sum of atomic mass * density - // const Options& allspecies = state["species"]; for (const auto& kv : allspecies.getChildren()) { const Options& species = kv.second; @@ -627,10 +627,10 @@ void RelaxPotential::finally(const Options& state) { const Options& allspecies = state["species"]; - auto* coord = mesh->getCoordinates(); + auto coord = mesh->getCoordinates(); phi = get(state["fields"]["phi"]); - Vort = get(state["fields"]["vorticity"]); + // Vort = get(state["fields"]["vorticity"]); // Solve vorticity equation @@ -667,15 +667,16 @@ void RelaxPotential::finally(const Options& state) { } } - if (state.isSection("fields") and state["fields"].isSet("DivJextra")) { + if (state["fields"].isSet("DivJextra")) { auto DivJextra = get(state["fields"]["DivJextra"]); + // Parallel current is handled here, to allow different 2D or 3D closures // to be used ddt(Vort) += DivJextra; } // Parallel current due to species parallel flow - for (auto& kv : allspecies.getChildren()) { + for (auto& kv : state["species"].getChildren()) { const Options& species = kv.second; if (!species.isSet("charge") or !species.isSet("momentum")) { @@ -709,7 +710,6 @@ void RelaxPotential::finally(const Options& state) { ddt(Vort) += FV::Div_a_Grad_perp(viscosity, Vort); ddt(Vort) += FV::Div_par_K_Grad_par(viscosity_par, Vort); - // NOTE(malamast): Need to check if Div_par_K_Grad_par is equivalent with the Laplace_par operator. // ddt(Vort) += viscosity_par * Laplace_par(Vort); @@ -723,14 +723,13 @@ void RelaxPotential::finally(const Options& state) { // Adds dissipation term like in other equations, but depending on gradient of // potential Field3D sound_speed = get(state["sound_speed"]); + ddt(Vort) -= FV::Div_par(-phi, 0.0, sound_speed); + } - Field3D zero {0.0}; - zero.splitParallelSlices(); - zero.yup() = 0.0; - zero.ydown() = 0.0; - - Field3D dummy; - ddt(Vort) -= FV::Div_par_mod(-phi, zero, sound_speed, dummy); + if (hyper_z > 0) { + // Form of hyper-viscosity to suppress zig-zags in Z + auto* coord = Vort.getCoordinates(); + ddt(Vort) -= hyper_z * SQ(SQ(coord->dz)) * D4DZ4(Vort); } if (phi_sheath_dissipation) { @@ -772,8 +771,7 @@ void RelaxPotential::finally(const Options& state) { } } - - + // Solve diffusion equation for potential if (boussinesq) { diff --git a/src/vorticity.cxx b/src/vorticity.cxx index b4ed63d23..e95fc95e0 100644 --- a/src/vorticity.cxx +++ b/src/vorticity.cxx @@ -53,7 +53,7 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) { const BoutReal Lnorm = units["meters"]; exb_advection = options["exb_advection"] - .doc("Include ExB advection (nonlinear term)?") + .doc("Include nonlinear ExB advection?") .withDefault(true); exb_advection_simplified = options["exb_advection_simplified"] From 3ddf6e37951b92088dc8c10626f4015f11927d68 Mon Sep 17 00:00:00 2001 From: malamast Date: Wed, 14 Jan 2026 14:43:22 -0800 Subject: [PATCH 4/4] vorticity: Added an option to add phi in the system of equations and solve it as a constraint eq instead of solving it with the laplacian solver. --- include/relax_potential.hxx | 5 +- include/vorticity.hxx | 9 +- src/relax_potential.cxx | 14 +-- src/vorticity.cxx | 228 ++++++++++++++++++++++++++++-------- 4 files changed, 197 insertions(+), 59 deletions(-) diff --git a/include/relax_potential.hxx b/include/relax_potential.hxx index e1341e384..1644b2d38 100644 --- a/include/relax_potential.hxx +++ b/include/relax_potential.hxx @@ -109,11 +109,12 @@ private: Field3D Pi_hat; ///< Contribution from ion pressure, weighted by atomic mass / charge + bool boussinesq; ///< Use the Boussinesq approximation? + bool exb_advection; //< Include nonlinear ExB advection? bool exb_advection_simplified; // Simplify nonlinear ExB advection form? bool diamagnetic; //< Include diamagnetic current? bool diamagnetic_polarisation; //< Include diamagnetic drift in polarisation current? - bool boussinesq; ///< Use the Boussinesq approximation? BoutReal average_atomic_mass; //< Weighted average atomic mass, for polarisaion current // (Boussinesq approximation) bool poloidal_flows; ///< Include poloidal ExB flow? @@ -137,7 +138,7 @@ private: 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 viscosity_par; ///< Parallel Kinematic viscosity // Relax-potential related variables BoutReal lambda_1; ///< Relaxation parameters. NOTE: lambda_1 has dimentions! diff --git a/include/vorticity.hxx b/include/vorticity.hxx index 42218caab..54dba2565 100644 --- a/include/vorticity.hxx +++ b/include/vorticity.hxx @@ -111,11 +111,15 @@ private: Field3D Pi_hat; ///< Contribution from ion pressure, weighted by atomic mass / charge + bool constraint; + bool boussinesq; ///< Use the Boussinesq approximation? + bool exb_advection; // Include nonlinear ExB advection? bool exb_advection_simplified; // Simplify nonlinear ExB advection form? bool diamagnetic; // Include diamagnetic current? bool diamagnetic_polarisation; // Include diamagnetic drift in polarisation current - BoutReal average_atomic_mass; // Weighted average atomic mass, for polarisaion current (Boussinesq approximation) + BoutReal average_atomic_mass; // Weighted average atomic mass, for polarisaion current + // (Boussinesq approximation) bool poloidal_flows; ///< Include poloidal ExB flow? bool bndry_flux; ///< Allow flows through radial boundaries? @@ -139,7 +143,8 @@ private: Field2D Bsq; // SQ(coord->Bxy) Vector2D Curlb_B; // Curvature vector Curl(b/B) BoutReal hyper_z; ///< Hyper-viscosity in Z - Field2D viscosity; ///< Kinematic viscosity + Field2D viscosity; ///< Perpendicular Kinematic viscosity + Field2D viscosity_par; ///< Parallel Kinematic viscosity // Diagnostic outputs Field3D DivJdia, DivJcol; // Divergence of diamagnetic and collisional current diff --git a/src/relax_potential.cxx b/src/relax_potential.cxx index f61988f31..a6f5f0760 100644 --- a/src/relax_potential.cxx +++ b/src/relax_potential.cxx @@ -52,6 +52,10 @@ RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* so const BoutReal Lnorm = units["meters"]; const BoutReal Tnorm = units["eV"]; const BoutReal Nnorm = units["inv_meters_cubed"]; + + boussinesq = options["boussinesq"] + .doc("Use the Boussinesq approximation?") + .withDefault(true); exb_advection = options["exb_advection"] .doc("Include nonlinear ExB advection?") @@ -78,10 +82,6 @@ RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* so .doc("Damp vorticity based on mass-weighted collision frequency") .withDefault(false); - boussinesq = options["boussinesq"] - .doc("Use the Boussinesq approximation?") - .withDefault(true); - average_atomic_mass = options["average_atomic_mass"] .doc("Weighted average atomic mass, for polarisation current " "(Boussinesq approximation)") @@ -139,14 +139,12 @@ RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* so lambda_1 = options["lambda_1"].doc("λ_1 > 1").withDefault(3e8) / (Tnorm * Omega_ci / SI::qe / Nnorm); lambda_2 = options["lambda_2"].doc("λ_2 > λ_1").withDefault(1.0); - // Add phi to restart files so that the value in the boundaries // is restored on restart. This is done even when phi is not evolving, // so that phi can be saved and re-loaded // Set initial value. Will be overwritten if restarting - phi1 = 0.0; - Vort = 0.0; + phi1 = 0.0; phi = 0.0; Vort = 0.0; auto coord = mesh->getCoordinates(); @@ -163,6 +161,7 @@ RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* so .doc("Timescale for phi boundary relaxation [seconds]") .withDefault(1e-4) / get(alloptions["units"]["seconds"]); + // Normalise to internal time units } // Read curvature vector @@ -771,7 +770,6 @@ void RelaxPotential::finally(const Options& state) { } } - // Solve diffusion equation for potential if (boussinesq) { diff --git a/src/vorticity.cxx b/src/vorticity.cxx index e95fc95e0..14ad98c22 100644 --- a/src/vorticity.cxx +++ b/src/vorticity.cxx @@ -52,6 +52,14 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) { const BoutReal Bnorm = units["Tesla"]; const BoutReal Lnorm = units["meters"]; + constraint = options["constraint"] + .doc("Solve DAE with constraint for phi?") + .withDefault(false); + + boussinesq = options["boussinesq"] + .doc("Use the Boussinesq approximation when phi is constrainted?") + .withDefault(true); + exb_advection = options["exb_advection"] .doc("Include nonlinear ExB advection?") .withDefault(true); @@ -93,12 +101,21 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) { .doc("Split phi into n=0 and n!=0 components") .withDefault(false); + viscosity = 0.0; viscosity = options["viscosity"] - .doc("Kinematic viscosity [m^2/s]") - .withDefault(0.0) + .doc("Perpendicular Kinematic viscosity [m^2/s]") + .withDefault(viscosity) / (Lnorm * Lnorm * Omega_ci); viscosity.applyBoundary("dirichlet"); + viscosity_par = 0.0; + viscosity_par = options["viscosity_par"] + .doc("Parallel Kinematic viscosity [m^2/s]") + .withDefault(viscosity_par) + / (Lnorm * Lnorm * Omega_ci); + viscosity_par.applyBoundary("dirichlet"); + viscosity_par.applyParallelBoundary("parallel_dirichlet_o2"); + hyper_z = options["hyper_z"].doc("Hyper-viscosity in Z. < 0 -> off").withDefault(-1.0); // Numerical dissipation terms @@ -131,7 +148,7 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) { // so that phi can be saved and re-loaded // Set initial value. Will be overwritten if restarting - phi = 0.0; + phi = 0.0; Vort = 0.0; auto coord = mesh->getCoordinates(); @@ -164,6 +181,12 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) { phiSolver->setOuterBoundaryFlags(INVERT_SET); } + if (constraint) { + // phi = phiSolver->solve(Vort * (SQ(coord->Bxy) / average_atomic_mass)); //NOTE: get initial value based on P? + // Add phi equation as a constraint + solver->constraint(phi, ddt(phi), "phi"); + } + // Read curvature vector try { Curlb_B.covariant = false; // Contravariant @@ -210,10 +233,51 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) { void Vorticity::transform(Options& state) { AUTO_TRACE(); + Options& allspecies = state["species"]; + phi.name = "phi"; auto& fields = state["fields"]; - // Set the boundary of phi. Both 2D and 3D fields are kept, though the 3D field + // 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(); + Field3D &Vort_yup = Vort.yup(); + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + // Vort_ydown(r.ind, mesh->ystart - 1, jz) = 2 * Vort(r.ind, mesh->ystart, jz) - Vort_yup(r.ind, mesh->ystart + 1, jz); + Vort_ydown(r.ind, mesh->ystart - 1, jz) = Vort(r.ind, mesh->ystart, jz); // TODO: iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + // Vort_yup(r.ind, mesh->yend + 1, jz) = 2 * Vort(r.ind, mesh->yend, jz) - Vort_ydown(r.ind, mesh->yend - 1, jz); + Vort_yup(r.ind, mesh->yend + 1, jz) = Vort(r.ind, mesh->yend, jz); + } + } + } else { + Field3D Vort_fa = toFieldAligned(Vort); + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + // Vort_fa(r.ind, mesh->ystart - 1, jz) = 2 * Vort_fa(r.ind, mesh->ystart, jz) - Vort_fa(r.ind, mesh->ystart + 1, jz); + Vort_fa(r.ind, mesh->ystart - 1, jz) = Vort_fa(r.ind, mesh->ystart, jz); // Neumann BC + } + } + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + // Vort_fa(r.ind, mesh->yend + 1, jz) = 2 * Vort_fa(r.ind, mesh->yend, jz) - Vort_fa(r.ind, mesh->yend - 1, jz); + Vort_fa(r.ind, mesh->yend + 1, jz) = Vort_fa(r.ind, mesh->yend, jz); // Neumann BC + } + } + Vort = fromFieldAligned(Vort_fa); + } + + // Set the boundary of phi. + // phi.applyBoundary("neumann"); + + // Both 2D and 3D fields are kept, though the 3D field // is constant in Z. This is for efficiency, to reduce the number of conversions. // 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 @@ -221,7 +285,6 @@ void Vorticity::transform(Options& state) { if (diamagnetic_polarisation) { // Diamagnetic term in vorticity. Note this is weighted by the mass // This includes all species, including electrons - Options& allspecies = state["species"]; for (auto& kv : allspecies.getChildren()) { Options& species = allspecies[kv.first]; // Note: need non-const @@ -398,57 +461,58 @@ void Vorticity::transform(Options& state) { } } - // Update boundary conditions. Two issues: - // 1) Solving here for phi + Pi, and then subtracting Pi from the result - // The boundary values should therefore include Pi - // 2) The INVERT_SET flag takes the value in the guard (boundary) cell - // and sets the boundary between cells to this value. - // This shift by 1/2 grid cell is important. + if (!constraint) { - Field3D phi_plus_pi = phi + Pi_hat; + // Update boundary conditions. Two issues: + // 1) Solving here for phi + Pi, and then subtracting Pi from the result + // The boundary values should therefore include Pi + // 2) The INVERT_SET flag takes the value in the guard (boundary) cell + // and sets the boundary between cells to this value. + // This shift by 1/2 grid cell is important. - if (mesh->firstX()) { - for (int j = mesh->ystart; j <= mesh->yend; j++) { - for (int k = 0; k < mesh->LocalNz; k++) { - // Average phi + Pi at the boundary, and set the boundary cell - // to this value. The phi solver will then put the value back - // onto the cell mid-point - phi_plus_pi(mesh->xstart - 1, j, k) = - 0.5 * (phi_plus_pi(mesh->xstart - 1, j, k) + phi_plus_pi(mesh->xstart, j, k)); + Field3D phi_plus_pi = phi + Pi_hat; + + if (mesh->firstX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + // Average phi + Pi at the boundary, and set the boundary cell + // to this value. The phi solver will then put the value back + // onto the cell mid-point + phi_plus_pi(mesh->xstart - 1, j, k) = + 0.5 * (phi_plus_pi(mesh->xstart - 1, j, k) + phi_plus_pi(mesh->xstart, j, k)); + } } } - } - if (mesh->lastX()) { - for (int j = mesh->ystart; j <= mesh->yend; j++) { - for (int k = 0; k < mesh->LocalNz; k++) { - phi_plus_pi(mesh->xend + 1, j, k) = - 0.5 * (phi_plus_pi(mesh->xend + 1, j, k) + phi_plus_pi(mesh->xend, j, k)); + if (mesh->lastX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + phi_plus_pi(mesh->xend + 1, j, k) = + 0.5 * (phi_plus_pi(mesh->xend + 1, j, k) + phi_plus_pi(mesh->xend, j, k)); + } } } - } - // Calculate potential - if (split_n0) { - //////////////////////////////////////////// - // Split into axisymmetric and non-axisymmetric components - Field2D Vort2D = DC(Vort); // n=0 component - Field2D phi_plus_pi_2d = DC(phi_plus_pi); - phi_plus_pi -= phi_plus_pi_2d; + // Calculate potential + if (split_n0) { + //////////////////////////////////////////// + // Split into axisymmetric and non-axisymmetric components + Field2D Vort2D = DC(Vort); // n=0 component + Field2D phi_plus_pi_2d = DC(phi_plus_pi); + phi_plus_pi -= phi_plus_pi_2d; - phi_plus_pi_2d = laplacexy->solve(Vort2D, phi_plus_pi_2d); + phi_plus_pi_2d = laplacexy->solve(Vort2D, phi_plus_pi_2d); - // Solve non-axisymmetric part using X-Z solver - phi = phi_plus_pi_2d - + phiSolver->solve((Vort - Vort2D) * (Bsq / average_atomic_mass), phi_plus_pi) - - Pi_hat; + // Solve non-axisymmetric part using X-Z solver + phi = phi_plus_pi_2d + + phiSolver->solve((Vort - Vort2D) * (Bsq / average_atomic_mass), phi_plus_pi) + - Pi_hat; - } else { - phi = phiSolver->solve(Vort * (Bsq / average_atomic_mass), phi_plus_pi) - Pi_hat; - } + } else { + phi = phiSolver->solve(Vort * (Bsq / average_atomic_mass), phi_plus_pi) - Pi_hat; + } - // Ensure that potential is set in the communication guard cells - mesh->communicate(phi); + } // Outer boundary cells if (mesh->firstX()) { @@ -470,6 +534,11 @@ void Vorticity::transform(Options& state) { } } + // Ensure that potential is set in the communication guard cells + // mesh->communicate(phi); + mesh->communicate(phi, Vort); //NOTE: Should we include Vort? + + // Vorticity equation ddt(Vort) = 0.0; if (diamagnetic) { @@ -515,8 +584,6 @@ void Vorticity::transform(Options& state) { Vector3D Grad_phi = Grad(phi); - Options& allspecies = state["species"]; - for (auto& kv : allspecies.getChildren()) { Options& species = allspecies[kv.first]; // Note: need non-const @@ -590,7 +657,6 @@ void Vorticity::transform(Options& state) { zeroFrom(Vort); // Sum of atomic mass * collision frequency * density Field3D sum_A_n = zeroFrom(Vort); // Sum of atomic mass * density - const Options& allspecies = state["species"]; for (const auto& kv : allspecies.getChildren()) { const Options& species = kv.second; @@ -626,9 +692,15 @@ void Vorticity::transform(Options& state) { void Vorticity::finally(const Options& state) { AUTO_TRACE(); + + const Options& allspecies = state["species"]; + auto coord = mesh->getCoordinates(); phi = get(state["fields"]["phi"]); + // Vort = get(state["fields"]["vorticity"]); + + // Solve vorticity equation if (exb_advection) { // These terms come from divergence of polarisation current @@ -705,6 +777,9 @@ void Vorticity::finally(const Options& state) { // Viscosity ddt(Vort) += FV::Div_a_Grad_perp(viscosity, Vort); + ddt(Vort) += FV::Div_par_K_Grad_par(viscosity_par, Vort); + + if (vort_dissipation) { // Adds dissipation term like in other equations Field3D sound_speed = get(state["sound_speed"]); @@ -762,6 +837,65 @@ void Vorticity::finally(const Options& state) { } } } + + if (constraint) { + + // Solve diffusion equation for potential + if (boussinesq) { + ddt(phi) = (FV::Div_a_Grad_perp(average_atomic_mass / Bsq, phi) - Vort); + + if (diamagnetic_polarisation) { + for (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(phi) += 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 (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); + + 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(phi) = phi_vort - Vort; + } + + } } void Vorticity::outputVars(Options& state) {