From 5ef65bf50a2bbceffb3ada0b7475c25556675ba9 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 18 Jun 2026 14:39:23 +0200 Subject: [PATCH 01/24] Update BOUT++ to include more fci bits --- external/BOUT-dev | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/BOUT-dev b/external/BOUT-dev index 8218f71a2..3b0117ea1 160000 --- a/external/BOUT-dev +++ b/external/BOUT-dev @@ -1 +1 @@ -Subproject commit 8218f71a2fe495bc138c7ddf7c33bb1be9837525 +Subproject commit 3b0117ea11652b576fbcf8e2dcb2be9163227ce3 From f99fc41e2b0e0ab1e817ab3b547f52b9b76c9b73 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 22 Jun 2026 10:04:45 +0200 Subject: [PATCH 02/24] Fix hermes-3.cxx for 3D compilation --- hermes-3.cxx | 220 ++++++++++++++++++++++++--------------------------- 1 file changed, 105 insertions(+), 115 deletions(-) diff --git a/hermes-3.cxx b/hermes-3.cxx index 5a06834b9..b1582091d 100644 --- a/hermes-3.cxx +++ b/hermes-3.cxx @@ -96,8 +96,8 @@ BOUT_OVERRIDE_DEFAULT_OPTION("mesh:calcParallelSlices_on_communicate", false); class DecayLengthBoundary : public BoundaryOp { public: DecayLengthBoundary() : gen(nullptr) {} - DecayLengthBoundary(BoundaryRegion* region, std::shared_ptr g) - : BoundaryOp(region), gen(std::move(g)) {} + DecayLengthBoundary(BoundaryRegion* region, std::shared_ptr g) + : BoundaryOp(region), gen(std::move(g)) {} using BoundaryOp::clone; /// Create a copy of this boundary condition @@ -118,10 +118,11 @@ class DecayLengthBoundary : public BoundaryOp { ASSERT1(mesh == f.getMesh()); // Get cell radial length - Coordinates *coord = mesh->getCoordinates(); - Field2D dx = coord->dx; - Field2D g11 = coord->g11; - Field2D dr = dx / sqrt(g11); // cell radial length. dr = dx/(Bpol * R) and g11 = (Bpol*R)**2 + Coordinates* coord = mesh->getCoordinates(); + auto dx = coord->dx; + auto g11 = coord->g11; + auto dr = + dx / sqrt(g11); // cell radial length. dr = dx/(Bpol * R) and g11 = (Bpol*R)**2 // Only implemented for cell centre quantities ASSERT1(f.getLocation() == CELL_CENTRE); @@ -129,11 +130,13 @@ class DecayLengthBoundary : public BoundaryOp { // This loop goes over the first row of boundary cells (in X and Y) for (bndry->first(); !bndry->isDone(); bndry->next1d()) { for (int zk = 0; zk < mesh->LocalNz; zk++) { // Loop over Z points - BoutReal decay_length = 3; // Default decay length is 3 normalised units (usually ~3mm) + BoutReal decay_length = + 3; // Default decay length is 3 normalised units (usually ~3mm) if (gen) { // Pick up the boundary condition setting from the input file // Must be specified in normalised units like the other BCs inputs - decay_length = gen->generate(bout::generator::Context(bndry, zk, CELL_CENTRE, 0.0, mesh)); + decay_length = + gen->generate(bout::generator::Context(bndry, zk, CELL_CENTRE, 0.0, mesh)); } // Set value in inner guard cell f(bndry->x, bndry->y, zk) // using the final domain cell value f(bndry->x - bndry->bx, bndry->y - bndry->by, zk) @@ -142,19 +145,21 @@ class DecayLengthBoundary : public BoundaryOp { // (-1, 0) X inner boundary (Core or PF) // (0, 1) Y upper boundary (outer lower target) // (0, -1) Y lower boundary (inner lower target) - + // Distance between final cell centre and inner guard cell centre in normalised units - BoutReal distance = 0.5 * (dr(bndry->x, bndry->y) + - dr(bndry->x - bndry->bx, bndry->y - bndry->by)); + BoutReal distance = 0.5 + * (dr(bndry->x, bndry->y, zk) + + dr(bndry->x - bndry->bx, bndry->y - bndry->by, zk)); - // Exponential decay - f(bndry->x, bndry->y, zk) = - f(bndry->x - bndry->bx, bndry->y - bndry->by, zk) * exp(-1 * distance / decay_length); + // Exponential decay + f(bndry->x, bndry->y, zk) = f(bndry->x - bndry->bx, bndry->y - bndry->by, zk) + * exp(-1 * distance / decay_length); // Set any remaining guard cells (i.e. the outer guards) to the same value // Should the outer guards have the decay continue, or just copy what the inners have? for (int i = 1; i < bndry->width; i++) { - f(bndry->x + i * bndry->bx, bndry->y + i * bndry->by, zk) = f(bndry->x, bndry->y, zk); + f(bndry->x + i * bndry->bx, bndry->y + i * bndry->by, zk) = + f(bndry->x, bndry->y, zk); } } } @@ -170,8 +175,8 @@ class DecayLengthBoundary : public BoundaryOp { int Hermes::init(bool restarting) { - auto &options = Options::root()["hermes"]; - + auto& options = Options::root()["hermes"]; + output.write("\nGit Version of Hermes: {:s}\n", hermes::version::revision); options["revision"] = hermes::version::revision; options["revision"].setConditionallyUsed(); @@ -193,7 +198,7 @@ int Hermes::init(bool restarting) { units["inv_meters_cubed"] = Nnorm; units["eV"] = Tnorm; units["Tesla"] = Bnorm; - units["seconds"] = 1./Omega_ci; + units["seconds"] = 1. / Omega_ci; units["meters"] = rho_s0; // Put into the options tree, so quantities can be normalised @@ -212,11 +217,13 @@ int Hermes::init(bool restarting) { TRACE("Loading metric tensor"); if (options.isSet("loadmetric")) { - throw BoutException("Error: The loadmetric option has been renamed to recalculate_metric.\n" - "Note: The default (true/false) is inverted for the new option.\n" - "Setting recalculate_metric=false (the default) uses the metric in the grid file.\n" - "Setting recalculate_metric=true loads Rxy, Bpxy etc from the grid file.\n" - " This assumes an orthogonal coordinate system. See manual for details."); + throw BoutException( + "Error: The loadmetric option has been renamed to recalculate_metric.\n" + "Note: The default (true/false) is inverted for the new option.\n" + "Setting recalculate_metric=false (the default) uses the metric in the grid " + "file.\n" + "Setting recalculate_metric=true loads Rxy, Bpxy etc from the grid file.\n" + " This assumes an orthogonal coordinate system. See manual for details."); } if (options["recalculate_metric"] @@ -227,28 +234,23 @@ int Hermes::init(bool restarting) { } else { // Check that the grid file contains at least one metric tensor component // Note: Older grid files did not, so would silently default to the identity metric - if (!(mesh->sourceHasVar("J") or - mesh->sourceHasVar("g11") or - mesh->sourceHasVar("g22") or - mesh->sourceHasVar("g33") or - mesh->sourceHasVar("g12") or - mesh->sourceHasVar("g23") or - mesh->sourceHasVar("g13") or - mesh->sourceHasVar("g_11") or - mesh->sourceHasVar("g_22") or - mesh->sourceHasVar("g_33") or - mesh->sourceHasVar("g_12") or - mesh->sourceHasVar("g_23") or - mesh->sourceHasVar("g_13"))) { - throw BoutException("Grid input does not contain any metric components (J, g11, g_22 etc).\n" - "Set hermes:recalculate_metric=true to calculate from Rxy, Bpxy etc.\n" - "If the default identity metric is intended then set e.g. mesh:J=1\n"); + if (!(mesh->sourceHasVar("J") or mesh->sourceHasVar("g11") + or mesh->sourceHasVar("g22") or mesh->sourceHasVar("g33") + or mesh->sourceHasVar("g12") or mesh->sourceHasVar("g23") + or mesh->sourceHasVar("g13") or mesh->sourceHasVar("g_11") + or mesh->sourceHasVar("g_22") or mesh->sourceHasVar("g_33") + or mesh->sourceHasVar("g_12") or mesh->sourceHasVar("g_23") + or mesh->sourceHasVar("g_13"))) { + throw BoutException( + "Grid input does not contain any metric components (J, g11, g_22 etc).\n" + "Set hermes:recalculate_metric=true to calculate from Rxy, Bpxy etc.\n" + "If the default identity metric is intended then set e.g. mesh:J=1\n"); } if (options["normalise_metric"] - .doc("Normalise input metric tensor? (assumes input is in SI units)") - .withDefault(true)) { - Coordinates *coord = mesh->getCoordinates(); + .doc("Normalise input metric tensor? (assumes input is in SI units)") + .withDefault(true)) { + Coordinates* coord = mesh->getCoordinates(); // To use non-orthogonal metric // Normalise coord->dx /= rho_s0 * rho_s0 * Bnorm; @@ -296,7 +298,7 @@ int Hermes::init(bool restarting) { 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.copy(); set(state["linear"], linear); @@ -330,82 +332,70 @@ void Hermes::outputVars(Options& options) { // Save normalisation quantities. These may be used by components // to calculate conversion factors to SI units - set_with_attrs(options["Tnorm"], Tnorm, { - {"units", "eV"}, - {"conversion", 1}, // Already in SI units - {"standard_name", "temperature normalisation"}, - {"long_name", "temperature normalisation"} - }); - set_with_attrs(options["Nnorm"], Nnorm, { - {"units", "m^-3"}, - {"conversion", 1}, - {"standard_name", "density normalisation"}, - {"long_name", "Number density normalisation"} - }); - set_with_attrs(options["Bnorm"], Bnorm, { - {"units", "T"}, - {"conversion", 1}, - {"standard_name", "magnetic field normalisation"}, - {"long_name", "Magnetic field normalisation"} - }); - set_with_attrs(options["Cs0"], Cs0, { - {"units", "m/s"}, - {"conversion", 1}, - {"standard_name", "velocity normalisation"}, - {"long_name", "Sound speed normalisation"} - }); - set_with_attrs(options["Omega_ci"], Omega_ci, { - {"units", "s^-1"}, - {"conversion", 1}, - {"standard_name", "frequency normalisation"}, - {"long_name", "Cyclotron frequency normalisation"} - }); - set_with_attrs(options["rho_s0"], rho_s0, { - {"units", "m"}, - {"conversion", 1}, - {"standard_name", "length normalisation"}, - {"long_name", "Gyro-radius length normalisation"} - }); + set_with_attrs(options["Tnorm"], Tnorm, + {{"units", "eV"}, + {"conversion", 1}, // Already in SI units + {"standard_name", "temperature normalisation"}, + {"long_name", "temperature normalisation"}}); + set_with_attrs(options["Nnorm"], Nnorm, + {{"units", "m^-3"}, + {"conversion", 1}, + {"standard_name", "density normalisation"}, + {"long_name", "Number density normalisation"}}); + set_with_attrs(options["Bnorm"], Bnorm, + {{"units", "T"}, + {"conversion", 1}, + {"standard_name", "magnetic field normalisation"}, + {"long_name", "Magnetic field normalisation"}}); + set_with_attrs(options["Cs0"], Cs0, + {{"units", "m/s"}, + {"conversion", 1}, + {"standard_name", "velocity normalisation"}, + {"long_name", "Sound speed normalisation"}}); + set_with_attrs(options["Omega_ci"], Omega_ci, + {{"units", "s^-1"}, + {"conversion", 1}, + {"standard_name", "frequency normalisation"}, + {"long_name", "Cyclotron frequency normalisation"}}); + set_with_attrs(options["rho_s0"], rho_s0, + {{"units", "m"}, + {"conversion", 1}, + {"standard_name", "length normalisation"}, + {"long_name", "Gyro-radius length normalisation"}}); scheduler->outputVars(options); } void Hermes::restartVars(Options& options) { - set_with_attrs(options["Tnorm"], Tnorm, { - {"units", "eV"}, - {"conversion", 1}, // Already in SI units - {"standard_name", "temperature normalisation"}, - {"long_name", "temperature normalisation"} - }); - set_with_attrs(options["Nnorm"], Nnorm, { - {"units", "m^-3"}, - {"conversion", 1}, - {"standard_name", "density normalisation"}, - {"long_name", "Number density normalisation"} - }); - set_with_attrs(options["Bnorm"], Bnorm, { - {"units", "T"}, - {"conversion", 1}, - {"standard_name", "magnetic field normalisation"}, - {"long_name", "Magnetic field normalisation"} - }); - set_with_attrs(options["Cs0"], Cs0, { - {"units", "m/s"}, - {"conversion", 1}, - {"standard_name", "velocity normalisation"}, - {"long_name", "Sound speed normalisation"} - }); - set_with_attrs(options["Omega_ci"], Omega_ci, { - {"units", "s^-1"}, - {"conversion", 1}, - {"standard_name", "frequency normalisation"}, - {"long_name", "Cyclotron frequency normalisation"} - }); - set_with_attrs(options["rho_s0"], rho_s0, { - {"units", "m"}, - {"conversion", 1}, - {"standard_name", "length normalisation"}, - {"long_name", "Gyro-radius length normalisation"} - }); + set_with_attrs(options["Tnorm"], Tnorm, + {{"units", "eV"}, + {"conversion", 1}, // Already in SI units + {"standard_name", "temperature normalisation"}, + {"long_name", "temperature normalisation"}}); + set_with_attrs(options["Nnorm"], Nnorm, + {{"units", "m^-3"}, + {"conversion", 1}, + {"standard_name", "density normalisation"}, + {"long_name", "Number density normalisation"}}); + set_with_attrs(options["Bnorm"], Bnorm, + {{"units", "T"}, + {"conversion", 1}, + {"standard_name", "magnetic field normalisation"}, + {"long_name", "Magnetic field normalisation"}}); + set_with_attrs(options["Cs0"], Cs0, + {{"units", "m/s"}, + {"conversion", 1}, + {"standard_name", "velocity normalisation"}, + {"long_name", "Sound speed normalisation"}}); + set_with_attrs(options["Omega_ci"], Omega_ci, + {{"units", "s^-1"}, + {"conversion", 1}, + {"standard_name", "frequency normalisation"}, + {"long_name", "Cyclotron frequency normalisation"}}); + set_with_attrs(options["rho_s0"], rho_s0, + {{"units", "m"}, + {"conversion", 1}, + {"standard_name", "length normalisation"}, + {"long_name", "Gyro-radius length normalisation"}}); scheduler->restartVars(options); } From 224caf1004fe2e9849e9525514d1563329658ce7 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 22 Jun 2026 10:05:22 +0200 Subject: [PATCH 03/24] Fix BraginskiiIonViscosity for 3D compilation --- include/braginskii_ion_viscosity.hxx | 4 +- src/braginskii_ion_viscosity.cxx | 61 ++++++++++++++-------------- 2 files changed, 32 insertions(+), 33 deletions(-) diff --git a/include/braginskii_ion_viscosity.hxx b/include/braginskii_ion_viscosity.hxx index 242759959..109e3b99c 100644 --- a/include/braginskii_ion_viscosity.hxx +++ b/include/braginskii_ion_viscosity.hxx @@ -8,7 +8,7 @@ #include #include -#include +#include #include "component.hxx" @@ -54,7 +54,7 @@ private: std::string viscosity_collisions_mode; ///< Collision selection, either multispecies or ///< braginskii Field3D nu; ///< Collision frequency for conduction - Vector2D Curlb_B; ///< Curvature vector Curl(b/B) + VectorMetric Curlb_B; ///< Curvature vector Curl(b/B) bool bounce_frequency; ///< Modify the collision time with the bounce frequency? BoutReal bounce_frequency_q95; ///< Input q95 for when including bounce frequency change BoutReal bounce_frequency_epsilon; ///< Input inverse aspect ratio for including bounce diff --git a/src/braginskii_ion_viscosity.cxx b/src/braginskii_ion_viscosity.cxx index 9578d1af5..9f90e26a9 100644 --- a/src/braginskii_ion_viscosity.cxx +++ b/src/braginskii_ion_viscosity.cxx @@ -21,6 +21,7 @@ #include #include #include +#include #include "../include/braginskii_ion_viscosity.hxx" #include "../include/component.hxx" @@ -32,16 +33,16 @@ using bout::globals::mesh; BraginskiiIonViscosity::BraginskiiIonViscosity(const std::string& name, Options& alloptions, Solver*) : Component({ - readIfSet("species:{non_electrons}:pressure"), - readIfSet("species:{non_electrons}:temperature"), - readIfSet("species:{non_electrons}:density"), - readIfSet("species:{non_electrons}:velocity"), - readIfSet("species:{non_electrons}:charge"), - readIfSet("species:{non_electrons}:collision_frequencies:{coll_type}"), - readWrite("species:{non_electrons}:momentum_source"), - readWrite("species:{non_electrons}:energy_source"), - readWrite("fields:DivJextra"), - }) { + readIfSet("species:{non_electrons}:pressure"), + readIfSet("species:{non_electrons}:temperature"), + readIfSet("species:{non_electrons}:density"), + readIfSet("species:{non_electrons}:velocity"), + readIfSet("species:{non_electrons}:charge"), + readIfSet("species:{non_electrons}:collision_frequencies:{coll_type}"), + readWrite("species:{non_electrons}:momentum_source"), + readWrite("species:{non_electrons}:energy_source"), + readWrite("fields:DivJextra"), + }) { auto& options = alloptions[name]; eta_limit_alpha = options["eta_limit_alpha"] @@ -52,8 +53,8 @@ BraginskiiIonViscosity::BraginskiiIonViscosity(const std::string& name, options["diagnose"].doc("Output additional diagnostics?").withDefault(false); parallel = options["parallel"] - .doc("Include parallel viscosity (requires parallel velocity)") - .withDefault(true); + .doc("Include parallel viscosity (requires parallel velocity)") + .withDefault(true); perpendicular = options["perpendicular"] .doc("Include perpendicular flow? (Requires phi)") @@ -93,7 +94,7 @@ BraginskiiIonViscosity::BraginskiiIonViscosity(const std::string& name, mesh->get(Curlb_B, "bxcv"); } catch (BoutException& e) { // May be 2D, reading as 3D - Vector2D curv2d; + VectorMetric curv2d; curv2d.covariant = false; mesh->get(curv2d, "bxcv"); Curlb_B = curv2d; @@ -143,9 +144,9 @@ void BraginskiiIonViscosity::transform_impl(GuardedOptions& state) { GuardedOptions allspecies = state["species"]; auto coord = mesh->getCoordinates(); - const Field2D Bxy = coord->Bxy; - const Field2D sqrtB = sqrt(Bxy); - const Field2D Grad_par_logB = Grad_par(log(Bxy)); + const auto Bxy = coord->Bxy; + const auto sqrtB = sqrt(Bxy); + const auto Grad_par_logB = Grad_par(log(Bxy)); // Loop through all species for (auto& kv : allspecies.getChildren()) { @@ -180,8 +181,8 @@ void BraginskiiIonViscosity::transform_impl(GuardedOptions& state) { const std::string collision_name = collision.first; if ( // Self-collisions - (collisionSpeciesMatch(collision_name, species_name, species_name, - "coll", "exact"))) { + (collisionSpeciesMatch(collision_name, species_name, species_name, "coll", + "exact"))) { collision_names[species_name].push_back(collision_name); } } @@ -237,10 +238,10 @@ void BraginskiiIonViscosity::transform_impl(GuardedOptions& state) { // Parallel ion viscosity (4/3 * 0.96 coefficient) Field3D eta = 1.28 * P * tau; - const Field2D tau_av = DC(tau); + const auto tau_av = DC(tau); Field2D bounce_factor = - 1.0; // if bounce_frequency = false, this factor does nothing to anything + 1.0; // if bounce_frequency = false, this factor does nothing to anything Field2D nu_star = 1.0; if (bounce_frequency) { @@ -255,11 +256,11 @@ void BraginskiiIonViscosity::transform_impl(GuardedOptions& state) { const Field2D v_thermal = sqrt(2.0 * T_av / mass); nu_star = (bounce_frequency_R * bounce_frequency_q95) - / (tau_av * pow(bounce_frequency_epsilon, 1.5) * v_thermal); + / (tau_av * pow(bounce_frequency_epsilon, 1.5) * v_thermal); bounce_factor *= - (1 / (1 + (1. / nu_star))) - * (1 / (1 + (1. / pow(bounce_frequency_epsilon, 1.5)) * (1. / nu_star))); + (1 / (1 + (1. / nu_star))) + * (1 / (1 + (1. / pow(bounce_frequency_epsilon, 1.5)) * (1. / nu_star))); eta *= bounce_factor; } @@ -289,8 +290,7 @@ void BraginskiiIonViscosity::transform_impl(GuardedOptions& state) { subtract(species["energy_source"], V * div_Pi_cipar); // Internal energy // Parallel ion stress tensor component - Pi_cipar = -0.96 * P * tau * bounce_factor - * (2. * Grad_par(V) + V * Grad_par_logB); + Pi_cipar = -0.96 * P * tau * bounce_factor * (2. * Grad_par(V) + V * Grad_par_logB); // Could also be written as: // Pi_cipar = -0.96*Pi*tau*2.*Grad_par(sqrt(Bxy)*Vi)/sqrt(Bxy); @@ -298,7 +298,7 @@ void BraginskiiIonViscosity::transform_impl(GuardedOptions& state) { const Field3D Pi_ciperp = 0.0; const Field3D DivJ = 0.0; diagnostics[species_name] = - Diagnostics{Pi_ciperp, Pi_cipar, DivJ, bounce_factor, nu_star}; + Diagnostics{Pi_ciperp, Pi_cipar, DivJ, bounce_factor, nu_star}; } } @@ -316,9 +316,8 @@ void BraginskiiIonViscosity::transform_impl(GuardedOptions& state) { // Perpendicular ion stress tensor // 0.96 P tau kappa * (V_E + V_di + 1.61 b x Grad(T)/B ) // Note: Heat flux terms are neglected for now - Field3D Pi_ciperp = - -0.5 * 0.96 * P * tau * bounce_factor - * (Curlb_B * Grad(phi + 1.61 * T) - Curlb_B * Grad(P) / Nlim); + Field3D Pi_ciperp = -0.5 * 0.96 * P * tau * bounce_factor + * (Curlb_B * Grad(phi + 1.61 * T) - Curlb_B * Grad(P) / Nlim); // Limit size of stress tensor components // If the off-diagonal components of the pressure tensor are large compared @@ -343,7 +342,7 @@ void BraginskiiIonViscosity::transform_impl(GuardedOptions& state) { const Field3D V = get(species["velocity"]); const Field3D div_Pi_ciperp = - -(2. / 3) * Grad_par(Pi_ciperp) + Pi_ciperp * Grad_par_logB; + -(2. / 3) * Grad_par(Pi_ciperp) + Pi_ciperp * Grad_par_logB; // const Field3D div_Pi_ciperp = - (2. / 3) * B32 * Grad_par(Pi_ciperp / B32); add(species["momentum_source"], div_Pi_ciperp); @@ -374,7 +373,7 @@ void BraginskiiIonViscosity::transform_impl(GuardedOptions& state) { // Transfer of energy between ion internal energy and ExB flow subtract(species["energy_source"], 0.5 * Pi_ci * Curlb_B * Grad(phi + P) - - (1. / 3.) * bracket(Pi_ci, phi + P, BRACKET_STD)); + - (1. / 3.) * bracket(Pi_ci, phi + P, BRACKET_STD)); if (diagnose) { diagnostics[species_name] = From 5b76db75dbbd65f1285418f64f7729392945a6d7 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 22 Jun 2026 10:05:53 +0200 Subject: [PATCH 04/24] Fix ClassicalDiffusion for 3D compilation --- include/classical_diffusion.hxx | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/include/classical_diffusion.hxx b/include/classical_diffusion.hxx index a87817a94..89e5e699b 100644 --- a/include/classical_diffusion.hxx +++ b/include/classical_diffusion.hxx @@ -2,17 +2,20 @@ #ifndef CLASSICAL_DIFFUSION_H #define CLASSICAL_DIFFUSION_H +#include + #include "component.hxx" struct ClassicalDiffusion : public Component { ClassicalDiffusion(std::string name, Options& alloptions, Solver*); - void outputVars(Options &state) override; + void outputVars(Options& state) override; + private: - Field2D Bsq; // Magnetic field squared + Coordinates::FieldMetric Bsq; // Magnetic field squared - bool diagnose; ///< Output additional diagnostics? - Field3D Dn; ///< Particle diffusion coefficient + bool diagnose; ///< Output additional diagnostics? + Field3D Dn; ///< Particle diffusion coefficient BoutReal custom_D; ///< User-set particle diffusion coefficient override // Flow diagnostics @@ -24,7 +27,8 @@ private: }; namespace { -RegisterComponent registercomponentclassicaldiffusion("classical_diffusion"); +RegisterComponent + registercomponentclassicaldiffusion("classical_diffusion"); } #endif // CLASSICAL_DIFFUSION_H From fef76a02877c1a236c959020ce3fdd20c188d1b5 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 22 Jun 2026 10:06:15 +0200 Subject: [PATCH 05/24] Fix DiamagneticDrift for 3D compilation --- include/diamagnetic_drift.hxx | 10 +++++----- src/diamagnetic_drift.cxx | 28 ++++++++++++++++++---------- 2 files changed, 23 insertions(+), 15 deletions(-) diff --git a/include/diamagnetic_drift.hxx b/include/diamagnetic_drift.hxx index 0ca0d8481..a5525d4e2 100644 --- a/include/diamagnetic_drift.hxx +++ b/include/diamagnetic_drift.hxx @@ -2,17 +2,19 @@ #ifndef DIAMAGNETIC_DRIFT_H #define DIAMAGNETIC_DRIFT_H +#include + #include "component.hxx" /// Calculate diamagnetic flows struct DiamagneticDrift : public Component { - DiamagneticDrift(std::string name, Options &options, Solver *UNUSED(solver)); + DiamagneticDrift(std::string name, Options& options, Solver* UNUSED(solver)); private: - Vector2D Curlb_B; + VectorMetric Curlb_B; bool bndry_flux; - Field2D diamag_form; + Coordinates::FieldMetric diamag_form; /// For every species, if it has: /// - temperature @@ -30,5 +32,3 @@ RegisterComponent registercomponentdiamagnetic("diamagnetic_dr } #endif // DIAMAGNETIC_DRIFT_H - - diff --git a/src/diamagnetic_drift.cxx b/src/diamagnetic_drift.cxx index 79f3c192c..bc2b0fb2d 100644 --- a/src/diamagnetic_drift.cxx +++ b/src/diamagnetic_drift.cxx @@ -17,8 +17,8 @@ DiamagneticDrift::DiamagneticDrift(std::string name, Options& alloptions, options["bndry_flux"].doc("Allow fluxes through boundary?").withDefault(true); diamag_form = options["diamag_form"] - .doc("Form of diamagnetic drift: 0 = gradient; 1 = divergence") - .withDefault(Field2D(1.0)); + .doc("Form of diamagnetic drift: 0 = gradient; 1 = divergence") + .withDefault(Field2D(1.0)); // Read curvature vector Curlb_B.covariant = false; // Contravariant @@ -27,8 +27,8 @@ DiamagneticDrift::DiamagneticDrift(std::string name, Options& alloptions, } Options& paralleltransform = Options::root()["mesh"]["paralleltransform"]; - if (paralleltransform.isSet("type") and - paralleltransform["type"].as() == "shifted") { + if (paralleltransform.isSet("type") + and paralleltransform["type"].as() == "shifted") { Field2D I; if (mesh->get(I, "sinty")) { I = 0.0; @@ -52,10 +52,14 @@ DiamagneticDrift::DiamagneticDrift(std::string name, Options& alloptions, // Set drift to zero through sheath boundaries. // Flux through those cell faces should be set by sheath. for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { - Curlb_B.y(r.ind, mesh->ystart - 1) = -Curlb_B.y(r.ind, mesh->ystart); + for (int k = 0; k < -Curlb_B.y.getNz(); ++k) { + Curlb_B.y(r.ind, mesh->ystart - 1, k) = -Curlb_B.y(r.ind, mesh->ystart, k); + } } for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { - Curlb_B.y(r.ind, mesh->yend + 1) = -Curlb_B.y(r.ind, mesh->yend); + for (int k = 0; k < -Curlb_B.y.getNz(); ++k) { + Curlb_B.y(r.ind, mesh->yend + 1, k) = -Curlb_B.y(r.ind, mesh->yend, k); + } } // FIXME: density, pressure, and momentum will not be read even if @@ -76,8 +80,9 @@ void DiamagneticDrift::transform_impl(GuardedOptions& state) { for (auto& kv : allspecies.getChildren()) { GuardedOptions species = allspecies[kv.first]; // Note: Need non-const - if (!(species.isSet("charge") and species.isSet("temperature"))) + if (!(species.isSet("charge") and species.isSet("temperature"))) { continue; // Skip, go to next species + } // Calculate diamagnetic drift velocity for this species auto q = get(species["charge"]); @@ -97,7 +102,8 @@ void DiamagneticDrift::transform_impl(GuardedOptions& state) { // Gradient form: Curlb_B dot Grad(N T / q) Field3D grad_form = Curlb_B * Grad(N * T / q); - subtract(species["density_source"], diamag_form * div_form + (1. - diamag_form) * grad_form); + subtract(species["density_source"], + diamag_form * div_form + (1. - diamag_form) * grad_form); } if (IS_SET(species["pressure"])) { @@ -105,14 +111,16 @@ void DiamagneticDrift::transform_impl(GuardedOptions& state) { Field3D div_form = FV::Div_f_v(P, vD, bndry_flux); Field3D grad_form = Curlb_B * Grad(P * T / q); - subtract(species["energy_source"], (5. / 2) * (diamag_form * div_form + (1. - diamag_form) * grad_form)); + subtract(species["energy_source"], + (5. / 2) * (diamag_form * div_form + (1. - diamag_form) * grad_form)); } if (IS_SET(species["momentum"])) { auto NV = get(species["momentum"]); Field3D div_form = FV::Div_f_v(NV, vD, bndry_flux); Field3D grad_form = Curlb_B * Grad(NV * T / q); - subtract(species["momentum_source"], diamag_form * div_form + (1. - diamag_form) * grad_form); + subtract(species["momentum_source"], + diamag_form * div_form + (1. - diamag_form) * grad_form); } } } From cd0afa1a2261308a144cfd874957b44a3502956c Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 22 Jun 2026 10:44:59 +0200 Subject: [PATCH 06/24] Div_par_mod has been moved to BOUT++ --- include/div_ops.hxx | 216 -------------------------------------------- 1 file changed, 216 deletions(-) diff --git a/include/div_ops.hxx b/include/div_ops.hxx index 6eaf81207..8e8f98519 100644 --- a/include/div_ops.hxx +++ b/include/div_ops.hxx @@ -612,222 +612,6 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in, return fromFieldAligned(result, "RGN_NOBNDRY"); } -/// Finite volume parallel divergence -/// -/// NOTE: Modified version, applies limiter to velocity and field -/// Performs better (smaller overshoots) than Div_par -/// -/// Preserves the sum of f*J*dx*dy*dz over the domain -/// -/// @param[in] f_in The field being advected. -/// This will be reconstructed at cell faces -/// using the given CellEdges method -/// @param[in] v_in The advection velocity. -/// This will be interpolated to cell boundaries -/// using linear interpolation -/// @param[in] wave_speed_in Local maximum speed of all waves in the system at each -// point in space -/// @param[in] fixflux Fix the flux at the boundary to be the value at the -/// midpoint (for boundary conditions) -/// -/// @param[out] flow_ylow Flow at the lower Y cell boundary -/// Already includes area factor * flux -/// -/// NB: Uses to/from FieldAligned coordinates -template -Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, - const Field3D& wave_speed_in, Field3D& flow_ylow, - bool fixflux = true) { - - ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); - ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in); - - Mesh* mesh = f_in.getMesh(); - - CellEdges cellboundary; - - ASSERT2(f_in.getDirectionY() == v_in.getDirectionY()); - ASSERT2(f_in.getDirectionY() == wave_speed_in.getDirectionY()); - const bool are_unaligned = - ((f_in.getDirectionY() == YDirectionType::Standard) - and (v_in.getDirectionY() == YDirectionType::Standard) - and (wave_speed_in.getDirectionY() == YDirectionType::Standard)); - - Field3D f = are_unaligned ? toFieldAligned(f_in, "RGN_NOX") : f_in; - Field3D v = are_unaligned ? toFieldAligned(v_in, "RGN_NOX") : v_in; - Field3D wave_speed = - are_unaligned ? toFieldAligned(wave_speed_in, "RGN_NOX") : wave_speed_in; - - Coordinates* coord = f_in.getCoordinates(); - - Field3D result{zeroFrom(f)}; - flow_ylow = zeroFrom(f); - - // Only need one guard cell, so no need to communicate fluxes - // Instead calculate in guard cells to preserve fluxes - int ys = mesh->ystart - 1; - int ye = mesh->yend + 1; - - for (int i = mesh->xstart; i <= mesh->xend; i++) { - - if (!mesh->firstY(i) || mesh->periodicY(i)) { - // Calculate in guard cell to get fluxes consistent between processors - ys = mesh->ystart - 1; - } else { - // Don't include the boundary cell. Note that this implies special - // handling of boundaries later - ys = mesh->ystart; - } - - if (!mesh->lastY(i) || mesh->periodicY(i)) { - // Calculate in guard cells - ye = mesh->yend + 1; - } else { - // Not in boundary cells - ye = mesh->yend; - } - - for (int j = ys; j <= ye; j++) { - // Pre-calculate factors which multiply fluxes -#if not(BOUT_USE_METRIC_3D) - // For right cell boundaries - BoutReal common_factor = (coord->J(i, j) + coord->J(i, j + 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); - - BoutReal flux_factor_rc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - BoutReal flux_factor_rp = - common_factor / (coord->dy(i, j + 1) * coord->J(i, j + 1)); - - BoutReal area_rp = common_factor * coord->dx(i, j + 1) * coord->dz(i, j + 1); - - // For left cell boundaries - common_factor = (coord->J(i, j) + coord->J(i, j - 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); - - BoutReal flux_factor_lc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - BoutReal flux_factor_lm = - common_factor / (coord->dy(i, j - 1) * coord->J(i, j - 1)); - - BoutReal area_lc = common_factor * coord->dx(i, j) * coord->dz(i, j); -#endif - for (int k = 0; k < mesh->LocalNz; k++) { -#if BOUT_USE_METRIC_3D - // For right cell boundaries - BoutReal common_factor = - (coord->J(i, j, k) + coord->J(i, j + 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); - - BoutReal flux_factor_rc = - common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); - BoutReal flux_factor_rp = - common_factor / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); - - BoutReal area_rp = - common_factor * coord->dx(i, j + 1, k) * coord->dz(i, j + 1, k); - - // For left cell boundaries - common_factor = (coord->J(i, j, k) + coord->J(i, j - 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); - - BoutReal flux_factor_lc = - common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); - BoutReal flux_factor_lm = - common_factor / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); - - BoutReal area_lc = common_factor * coord->dx(i, j, k) * coord->dz(i, j, k); -#endif - - //////////////////////////////////////////// - // Reconstruct f at the cell faces - // This calculates s.R and s.L for the Right and Left - // face values on this cell - - // Reconstruct f at the cell faces - Stencil1D s; - s.c = f(i, j, k); - s.m = f(i, j - 1, k); - s.p = f(i, j + 1, k); - - cellboundary(s); // Calculate s.R and s.L - - //////////////////////////////////////////// - // Reconstruct v at the cell faces - Stencil1D sv; - sv.c = v(i, j, k); - sv.m = v(i, j - 1, k); - sv.p = v(i, j + 1, k); - - cellboundary(sv); // Calculate sv.R and sv.L - - //////////////////////////////////////////// - // Right boundary - - BoutReal flux; - - if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) { - // Last point in domain - - // Calculate velocity at right boundary (y+1/2) - BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j + 1, k)); - - BoutReal bndryval = 0.5 * (s.c + s.p); - if (fixflux) { - // Use mid-point to be consistent with boundary conditions - flux = bndryval * vpar; - } else { - // Add flux due to difference in boundary values - flux = s.R * vpar + wave_speed(i, j, k) * (s.R - bndryval); - } - - } else { - // Maximum wave speed in the two cells - BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k), - fabs(v(i, j, k)), fabs(v(i, j + 1, k))); - - flux = s.R * 0.5 * (sv.R + amax); - } - - result(i, j, k) += flux * flux_factor_rc; - result(i, j + 1, k) -= flux * flux_factor_rp; - - flow_ylow(i, j + 1, k) += flux * area_rp; - - //////////////////////////////////////////// - // Calculate at left boundary - - if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) { - // First point in domain - BoutReal bndryval = 0.5 * (s.c + s.m); - BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j - 1, k)); - if (fixflux) { - // Use mid-point to be consistent with boundary conditions - flux = bndryval * vpar; - } else { - // Add flux due to difference in boundary values - flux = s.L * vpar - wave_speed(i, j, k) * (s.L - bndryval); - } - } else { - - // Maximum wave speed in the two cells - BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k), - fabs(v(i, j, k)), fabs(v(i, j - 1, k))); - - flux = s.L * 0.5 * (sv.L - amax); - } - - result(i, j, k) -= flux * flux_factor_lc; - result(i, j - 1, k) += flux * flux_factor_lm; - - flow_ylow(i, j, k) += flux * area_lc; - } - } - } - if (are_unaligned) { - flow_ylow = fromFieldAligned(flow_ylow, "RGN_NOBNDRY"); - } - return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; -} - /// Div ( a g Grad_perp(f) ) -- Perpendicular gradient-driven advection /// /// This version uses a slope limiter to calculate cell edge values of g in X, From 3c93c2df91db66e6e5e2fbf9ea2cddfd8e735824 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 22 Jun 2026 10:46:11 +0200 Subject: [PATCH 07/24] Div_par_fvv has been moved to BOUT++ --- include/div_ops.hxx | 207 -------------------------------------------- 1 file changed, 207 deletions(-) diff --git a/include/div_ops.hxx b/include/div_ops.hxx index 8e8f98519..783e42fdc 100644 --- a/include/div_ops.hxx +++ b/include/div_ops.hxx @@ -144,213 +144,6 @@ struct Superbee { } }; -/// This operator calculates Div_par(f v v) -/// It is used primarily (only?) in the parallel momentum equation. -/// -/// This operator is used rather than Div(f fv) so that the values of -/// f and v are consistent with other advection equations: The product -/// fv is not interpolated to cell boundaries. -template -Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, - const Field3D& wave_speed_in, bool fixflux = true) { - ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); - Mesh* mesh = f_in.getMesh(); - Coordinates* coord = f_in.getCoordinates(); - CellEdges cellboundary; - - if (f_in.isFci()) { - // FCI version, using yup/down fields - ASSERT1(f_in.hasParallelSlices()); - ASSERT1(v_in.hasParallelSlices()); - - const auto B = coord->Bxy; - const auto B_up = coord->Bxy.yup(); - const auto B_down = coord->Bxy.ydown(); - - const auto& f_up = f_in.yup(); - const auto& f_down = f_in.ydown(); - - const auto& v_up = v_in.yup(); - const auto& v_down = v_in.ydown(); - - const auto g_22 = coord->g_22; - const auto dy = coord->dy; - - Field3D result{emptyFrom(f_in)}; - BOUT_FOR(i, f_in.getRegion("RGN_NOBNDRY")) { - const auto iyp = i.yp(); - const auto iym = i.ym(); - - //Maximum local wave speed - const BoutReal amax = - BOUTMAX(wave_speed_in[i], fabs(v_in[i]), fabs(v_up[iyp]), fabs(v_down[iym])); - - result[i] = - B[i] - * ((f_up[iyp] * v_up[iyp] * v_up[iyp] / B_up[iyp]) - - (f_down[iym] * v_down[iym] * v_down[iym] / B_down[iym]) - // Penalty terms. This implementation is very dissipative. - // Note: This version adds a viscosity that damps gradients of velocity - + amax * (f_in[i] + f_up[iyp]) * (v_in[i] - v_up[iyp]) / (B[i] + B_up[iyp]) - + amax * (f_in[i] + f_down[iym]) * (v_in[i] - v_down[iym]) - / (B[i] + B_down[iym])) - / (2 * dy[i] * sqrt(g_22[i])); - -#if CHECK > 0 - if (!std::isfinite(result[i])) { - throw BoutException("Non-finite value in Div_par_fvv at {}\n" - "fup {} vup {} fdown {} vdown {} amax {}\n", - "B {} Bup {} Bdown {} dy {} sqrt(g_22} {}", i, f_up[i], - v_up[i], f_down[i], v_down[i], amax, B[i], B_up[i], B_down[i], - dy[i], sqrt(g_22[i])); - } -#endif - } - return result; - } - - ASSERT1(areFieldsCompatible(f_in, wave_speed_in)); - - /// Ensure that f, v and wave_speed are field aligned - Field3D f = toFieldAligned(f_in, "RGN_NOX"); - Field3D v = toFieldAligned(v_in, "RGN_NOX"); - Field3D wave_speed = toFieldAligned(wave_speed_in, "RGN_NOX"); - - Field3D result{zeroFrom(f)}; - - // Only need one guard cell, so no need to communicate fluxes - // Instead calculate in guard cells to preserve fluxes - int ys = mesh->ystart - 1; - int ye = mesh->yend + 1; - - for (int i = mesh->xstart; i <= mesh->xend; i++) { - - if (!mesh->firstY(i) || mesh->periodicY(i)) { - // Calculate in guard cell to get fluxes consistent between processors - ys = mesh->ystart - 1; - } else { - // Don't include the boundary cell. Note that this implies special - // handling of boundaries later - ys = mesh->ystart; - } - - if (!mesh->lastY(i) || mesh->periodicY(i)) { - // Calculate in guard cells - ye = mesh->yend + 1; - } else { - // Not in boundary cells - ye = mesh->yend; - } - - for (int j = ys; j <= ye; j++) { - // Pre-calculate factors which multiply fluxes - - for (int k = 0; k < mesh->LocalNz; k++) { - // For right cell boundaries - BoutReal common_factor = - (coord->J(i, j, k) + coord->J(i, j + 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); - - const BoutReal flux_factor_rc = - common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); - const BoutReal flux_factor_rp = - common_factor / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); - - // For left cell boundaries - common_factor = (coord->J(i, j, k) + coord->J(i, j - 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); - - const BoutReal flux_factor_lc = - common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); - const BoutReal flux_factor_lm = - common_factor / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); - - //////////////////////////////////////////// - // Reconstruct f at the cell faces - // This calculates s.R and s.L for the Right and Left - // face values on this cell - - // Reconstruct f at the cell faces - Stencil1D s; - s.c = f(i, j, k); - s.m = f(i, j - 1, k); - s.p = f(i, j + 1, k); - - cellboundary(s); // Calculate s.R and s.L - - // Reconstruct v at the cell faces - Stencil1D sv; - sv.c = v(i, j, k); - sv.m = v(i, j - 1, k); - sv.p = v(i, j + 1, k); - - cellboundary(sv); - - //////////////////////////////////////////// - // Right boundary - - // Calculate velocity at right boundary (y+1/2) - BoutReal v_mid = 0.5 * (sv.c + sv.p); - // And mid-point density at right boundary - BoutReal n_mid = 0.5 * (s.c + s.p); - BoutReal flux; - - if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) { - // Last point in domain - - if (fixflux) { - // Use mid-point to be consistent with boundary conditions - flux = n_mid * v_mid * v_mid; - } else { - // Add flux due to difference in boundary values - flux = s.R * sv.R * sv.R // Use right cell edge values - + BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.p)) * n_mid - * (sv.R - v_mid); // Damp differences in velocity, not flux - } - } else { - // Maximum wave speed in the two cells - const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k), - fabs(sv.c), fabs(sv.p)); - - flux = s.R * 0.5 * (sv.R + amax) * sv.R; - } - - result(i, j, k) += flux * flux_factor_rc; - result(i, j + 1, k) -= flux * flux_factor_rp; - - //////////////////////////////////////////// - // Calculate at left boundary - - v_mid = 0.5 * (sv.c + sv.m); - n_mid = 0.5 * (s.c + s.m); - - if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) { - // First point in domain - if (fixflux) { - // Use mid-point to be consistent with boundary conditions - flux = n_mid * v_mid * v_mid; - } else { - // Add flux due to difference in boundary values - flux = s.L * sv.L * sv.L - - BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.m)) * n_mid - * (sv.L - v_mid); - } - } else { - // Maximum wave speed in the two cells - const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k), - fabs(sv.c), fabs(sv.m)); - - flux = s.L * 0.5 * (sv.L - amax) * sv.L; - } - - result(i, j, k) -= flux * flux_factor_lc; - result(i, j - 1, k) += flux * flux_factor_lm; - } - } - } - return fromFieldAligned(result, "RGN_NOBNDRY"); -} - // Calculates viscous heating due to numerical momentum fluxes // and flow of kinetic energy (in flow_ylow) template From b544b6d5ad6087e03a6718154fbee43429d9623a Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 22 Jun 2026 10:47:31 +0200 Subject: [PATCH 08/24] Div_par_K_Grad_par_mod has been moved to BOUT++ --- include/div_ops.hxx | 4 - src/div_ops.cxx | 193 +++++++++++++++----------------------------- 2 files changed, 66 insertions(+), 131 deletions(-) diff --git a/include/div_ops.hxx b/include/div_ops.hxx index 783e42fdc..b24ff28bc 100644 --- a/include/div_ops.hxx +++ b/include/div_ops.hxx @@ -82,10 +82,6 @@ const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f); const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, Field3D& flux_xlow, Field3D& flux_ylow); -/// Version with energy flow diagnostic -const Field3D Div_par_K_Grad_par_mod(const Field3D& k, const Field3D& f, - Field3D& flow_ylow, bool bndry_flux = true); - /*! * Div ( a Grad_perp(f) ) -- ∇⊥ ( a ⋅ ∇⊥ f) -- Vorticity * diff --git a/src/div_ops.cxx b/src/div_ops.cxx index adcaf8572..d83c023d7 100644 --- a/src/div_ops.cxx +++ b/src/div_ops.cxx @@ -42,17 +42,19 @@ const Field3D Div_par_diffusion_index(const Field3D& f, bool bndry_flux) { Coordinates* coord = mesh->getCoordinates(); - for (int i = mesh->xstart; i <= mesh->xend; i++) - for (int j = mesh->ystart - 1; j <= mesh->yend; j++) + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart - 1; j <= mesh->yend; j++) { for (int k = 0; k < mesh->LocalNz; k++) { // Calculate flux at upper surface if (!bndry_flux && !mesh->periodicY(i)) { - if ((j == mesh->yend) && mesh->lastY(i)) + if ((j == mesh->yend) && mesh->lastY(i)) { continue; + } - if ((j == mesh->ystart - 1) && mesh->firstY(i)) + if ((j == mesh->ystart - 1) && mesh->firstY(i)) { continue; + } } BoutReal J = 0.5 * (coord->J(i, j) + coord->J(i, j + 1)); // Jacobian at boundary @@ -63,6 +65,8 @@ const Field3D Div_par_diffusion_index(const Field3D& f, bool bndry_flux) { result(i, j, k) += flux / coord->J(i, j); result(i, j + 1, k) -= flux / coord->J(i, j + 1); } + } + } return result; } @@ -91,8 +95,7 @@ BoutReal minmod(BoutReal a, BoutReal b, BoutReal c) { // Monotonized Central limiter (Van-Leer) void MC(Stencil1D& n) { - BoutReal slope = - minmod(2. * (n.p - n.c), 0.5 * (n.p - n.m), 2. * (n.c - n.m)); + BoutReal slope = minmod(2. * (n.p - n.c), 0.5 * (n.p - n.m), 2. * (n.c - n.m)); n.L = n.c - 0.5 * slope; n.R = n.c + 0.5 * slope; } @@ -126,8 +129,8 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd // int nz = mesh->LocalNz; - for (int i = mesh->xstart; i <= mesh->xend; i++) - for (int j = mesh->ystart; j <= mesh->yend; j++) + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { for (int k = 0; k < nz; k++) { int kp = (k + 1) % nz; int kpp = (kp + 1) % nz; @@ -249,6 +252,8 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd result(i, j, km) += flux; } } + } + } FV::communicateFluxes(result); ////////////////////////////////////////// @@ -288,8 +293,8 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd } } - for (int i = xs; i <= xe; i++) - for (int j = mesh->ystart - 1; j <= mesh->yend; j++) + for (int i = xs; i <= xe; i++) { + for (int j = mesh->ystart - 1; j <= mesh->yend; j++) { for (int k = 0; k < mesh->LocalNz; k++) { // Average dfdy to right X boundary @@ -325,6 +330,8 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); result(i + 1, j, k) -= flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); } + } + } } if (poloidal) { @@ -431,8 +438,8 @@ const Field3D Div_Perp_Lap_FV_Index(const Field3D& as, const Field3D& fs) { Coordinates* coord = mesh->getCoordinates(); - for (int i = mesh->xstart; i <= mesh->xend; i++) - for (int j = mesh->ystart; j <= mesh->yend; j++) + 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++) { int kp = (k + 1) % mesh->LocalNz; int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; @@ -469,17 +476,19 @@ const Field3D Div_Perp_Lap_FV_Index(const Field3D& as, const Field3D& fs) { flux = gD * 0.5 * (as(i, j, k) + as(i, j, km)); result(i, j, k) -= flux; } + } + } return result; } /// Z diffusion in index space -const Field3D Div_Z_FV_Index(const Field3D &as, const Field3D &fs) { +const Field3D Div_Z_FV_Index(const Field3D& as, const Field3D& fs) { Field3D result = 0.0; - for (int i = mesh->xstart; i <= mesh->xend; i++) - for (int j = mesh->ystart; j <= mesh->yend; j++) + 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++) { int kp = (k + 1) % mesh->LocalNz; int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; @@ -494,7 +503,9 @@ const Field3D Div_Z_FV_Index(const Field3D &as, const Field3D &fs) { result(i, j, k) -= gD * 0.5 * (as(i, j, k) + as(i, j, km)); } - + } + } + return result; } @@ -504,7 +515,7 @@ const Field3D D4DX4_FV_Index(const Field3D& f, bool bndry_flux) { Coordinates* coord = mesh->getCoordinates(); - for (int i = mesh->xstart; i <= mesh->xend; i++) + 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++) { @@ -576,6 +587,7 @@ const Field3D D4DX4_FV_Index(const Field3D& f, bool bndry_flux) { } } } + } return result; } @@ -584,7 +596,7 @@ const Field3D D4DZ4_Index(const Field3D& f) { Field3D result; result.allocate(); BOUT_FOR(i, f.getRegion("RGN_NOBNDRY")) { - result[i] = f[i.zp(2)] - 4.*f[i.zp()] + 6 * f[i] - 4 * f[i.zm()] + f[i.zm(2)]; + result[i] = f[i.zp(2)] - 4. * f[i.zp()] + 6 * f[i] - 4 * f[i.zm()] + f[i.zm(2)]; } return result; } @@ -602,7 +614,7 @@ const Field2D Laplace_FV(const Field2D& k, const Field2D& f) { Coordinates* coord = mesh->getCoordinates(); - for (int i = mesh->xstart; i <= mesh->xend; i++) + for (int i = mesh->xstart; i <= mesh->xend; i++) { for (int j = mesh->ystart; j <= mesh->yend; j++) { // Calculate gradients on cell faces @@ -641,12 +653,12 @@ const Field2D Laplace_FV(const Field2D& k, const Field2D& f) { flux = gD * 0.25 * (coord->J(i, j - 1) + coord->J(i, j)) * (k(i, j - 1) + k(i, j)); result(i, j) -= flux / (coord->dy(i, j) * coord->J(i, j)); } + } return result; } const Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, - Field3D &flow_xlow, - Field3D &flow_ylow) { + Field3D& flow_xlow, Field3D& flow_ylow) { ASSERT2(a.getLocation() == f.getLocation()); Mesh* mesh = a.getMesh(); @@ -890,7 +902,7 @@ const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { int xs = mesh->xstart - 1; int xe = mesh->xend; - for (int i = xs; i <= xe; i++) + for (int i = xs; i <= xe; 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 @@ -907,6 +919,7 @@ const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { result(i + 1, j, k) -= fout / (coord->dx(i + 1, j) * coord->J(i + 1, j)); } } + } // Y and Z fluxes require Y derivatives @@ -970,9 +983,9 @@ const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { / (coord->dy(i, j + 1) + coord->dy(i, j)); BoutReal fout = 0.25 * (ac(i, j, k) + aup(i, j + 1, k)) - * (coord->J(i, j) * coord->g23(i, j) - + coord->J(i, j + 1) * coord->g23(i, j + 1)) - * (dfdz - coef_u * dfdy); + * (coord->J(i, j) * coord->g23(i, j) + + coord->J(i, j + 1) * coord->g23(i, j + 1)) + * (dfdz - coef_u * dfdy); yzresult(i, j, k) = fout / (coord->dy(i, j) * coord->J(i, j)); @@ -1398,7 +1411,7 @@ Field3D Div_a_Grad_perp_nonorthog(const Field3D& a, const Field3D& f, Field3D& f /// Flows are always in the positive {x,y} direction /// i.e xlow(i,j) is the flow into cell (i,j) from the left, /// and the flow out of cell (i-1,j) to the right -/// +/// /// ylow(i,j+1) /// ^ /// +---|---+ @@ -1411,8 +1424,7 @@ Field3D Div_a_Grad_perp_nonorthog(const Field3D& a, const Field3D& f, Field3D& f /// /// WARNING: Causes checkerboarding in neutral_mixed integrated test const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, - Field3D &flow_xlow, - Field3D &flow_ylow) { + Field3D& flow_xlow, Field3D& flow_ylow) { ASSERT2(a.getLocation() == f.getLocation()); Mesh* mesh = a.getMesh(); @@ -1430,13 +1442,13 @@ const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, int xs = mesh->xstart - 1; int xe = mesh->xend; - for (int i = xs; i <= xe; i++) + for (int i = xs; i <= xe; 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 const BoutReal gradient = (coord->J(i, j) * coord->g11(i, j) - + coord->J(i + 1, j) * coord->g11(i + 1, j)) + + coord->J(i + 1, j) * coord->g11(i + 1, j)) * (f(i + 1, j, k) - f(i, j, k)) / (coord->dx(i, j) + coord->dx(i + 1, j)); @@ -1450,6 +1462,7 @@ const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, flow_xlow(i + 1, j, k) = -1.0 * fout * coord->dy(i, j) * coord->dz(i, j); } } + } // Y and Z fluxes require Y derivatives @@ -1514,9 +1527,9 @@ const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, / (coord->dy(i, j + 1) + coord->dy(i, j)); BoutReal fout = 0.25 * (ac(i, j, k) + aup(i, j + 1, k)) - * (coord->J(i, j) * coord->g23(i, j) - + coord->J(i, j + 1) * coord->g23(i, j + 1)) - * (dfdz - coef_u * dfdy); + * (coord->J(i, j) * coord->g23(i, j) + + coord->J(i, j + 1) * coord->g23(i, j + 1)) + * (dfdz - coef_u * dfdy); yzresult(i, j, k) = fout / (coord->dy(i, j) * coord->J(i, j)); @@ -1585,8 +1598,8 @@ const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, const Fi bool bndry_flux) { Field3D result{0.0}; - Coordinates *coord = mesh->getCoordinates(); - + Coordinates* coord = mesh->getCoordinates(); + ////////////////////////////////////////// // X-Z advection. // @@ -1613,24 +1626,25 @@ const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, const Fi // 1) Interpolate stream function f onto corners fmp, fpp, fpm - BoutReal fmm = 0.25 * (f(i, j, k) + f(i - 1, j, k) + f(i, j, km) + - f(i - 1, j, km)); - BoutReal fmp = 0.25 * (f(i, j, k) + f(i, j, kp) + f(i - 1, j, k) + - f(i - 1, j, kp)); // 2nd order accurate - BoutReal fpp = 0.25 * (f(i, j, k) + f(i, j, kp) + f(i + 1, j, k) + - f(i + 1, j, kp)); - BoutReal fpm = 0.25 * (f(i, j, k) + f(i + 1, j, k) + f(i, j, km) + - f(i + 1, j, km)); + BoutReal fmm = + 0.25 * (f(i, j, k) + f(i - 1, j, k) + f(i, j, km) + f(i - 1, j, km)); + BoutReal fmp = 0.25 + * (f(i, j, k) + f(i, j, kp) + f(i - 1, j, k) + + f(i - 1, j, kp)); // 2nd order accurate + BoutReal fpp = + 0.25 * (f(i, j, k) + f(i, j, kp) + f(i + 1, j, k) + f(i + 1, j, kp)); + BoutReal fpm = + 0.25 * (f(i, j, k) + f(i + 1, j, k) + f(i, j, km) + f(i + 1, j, km)); // 2) Calculate velocities on cell faces BoutReal vU = coord->J(i, j) * (fmp - fpp) / coord->dx(i, j); // -J*df/dx BoutReal vD = coord->J(i, j) * (fmm - fpm) / coord->dx(i, j); // -J*df/dx - BoutReal vR = 0.5 * (coord->J(i, j) + coord->J(i + 1, j)) * (fpp - fpm) / - coord->dz(i, j); // J*df/dz - BoutReal vL = 0.5 * (coord->J(i, j) + coord->J(i - 1, j)) * (fmp - fmm) / - coord->dz(i, j); // J*df/dz + BoutReal vR = 0.5 * (coord->J(i, j) + coord->J(i + 1, j)) * (fpp - fpm) + / coord->dz(i, j); // J*df/dz + BoutReal vL = 0.5 * (coord->J(i, j) + coord->J(i - 1, j)) * (fmp - fmm) + / coord->dz(i, j); // J*df/dz // 3) Calculate g on cell faces @@ -1667,8 +1681,7 @@ const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, const Fi } result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= - flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + result(i + 1, j, k) -= flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); } } else { // Not at a boundary @@ -1676,8 +1689,7 @@ const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, const Fi // Flux out into next cell BoutReal flux = vR * s.R * gR; result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= - flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + result(i + 1, j, k) -= flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); } } @@ -1698,8 +1710,7 @@ const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, const Fi flux = vL * 0.5 * (n(i - 1, j, k) + n(i, j, k)) * gL; } result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); - result(i - 1, j, k) += - flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); + result(i - 1, j, k) += flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); } } else { // Not at a boundary @@ -1707,8 +1718,7 @@ const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, const Fi if (vL < 0.0) { BoutReal flux = vL * s.L * gL; result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); - result(i - 1, j, k) += - flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); + result(i - 1, j, k) += flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); } } @@ -1723,7 +1733,7 @@ const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, const Fi MC(s); if (vU > 0.0) { - BoutReal flux = vU * s.R * gU/ (coord->J(i, j) * coord->dz(i, j)); + BoutReal flux = vU * s.R * gU / (coord->J(i, j) * coord->dz(i, j)); result(i, j, k) += flux; result(i, j, kp) -= flux; } @@ -1738,74 +1748,3 @@ const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, const Fi FV::communicateFluxes(result); return result; } - -const Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, - Field3D& flow_ylow, - bool bndry_flux) { - TRACE("FV::Div_par_K_Grad_par_mod"); - - ASSERT2(Kin.getLocation() == fin.getLocation()); - - Mesh* mesh = Kin.getMesh(); - - bool use_parallel_slices = (Kin.hasParallelSlices() && fin.hasParallelSlices()); - - const auto& K = use_parallel_slices ? Kin : toFieldAligned(Kin, "RGN_NOX"); - const auto& f = use_parallel_slices ? fin : toFieldAligned(fin, "RGN_NOX"); - - Field3D result{zeroFrom(f)}; - flow_ylow = zeroFrom(f); - - // K and f fields in yup and ydown directions - const auto& Kup = use_parallel_slices ? Kin.yup() : K; - const auto& Kdown = use_parallel_slices ? Kin.ydown() : K; - const auto& fup = use_parallel_slices ? fin.yup() : f; - const auto& fdown = use_parallel_slices ? fin.ydown() : f; - - Coordinates* coord = fin.getCoordinates(); - - BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { - // Calculate flux at upper surface - - const auto iyp = i.yp(); - const auto iym = i.ym(); - - if (bndry_flux || mesh->periodicY(i.x()) || !mesh->lastY(i.x()) - || (i.y() != mesh->yend)) { - - BoutReal c = 0.5 * (K[i] + Kup[iyp]); // K at the upper boundary - BoutReal J = 0.5 * (coord->J[i] + coord->J[iyp]); // Jacobian at boundary - BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); - - BoutReal gradient = 2. * (fup[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); - - BoutReal flux = c * J * gradient / g_22; - - result[i] += flux / (coord->dy[i] * coord->J[i]); - } - - // Calculate flux at lower surface - if (bndry_flux || mesh->periodicY(i.x()) || !mesh->firstY(i.x()) - || (i.y() != mesh->ystart)) { - BoutReal c = 0.5 * (K[i] + Kdown[iym]); // K at the lower boundary - BoutReal J = 0.5 * (coord->J[i] + coord->J[iym]); // Jacobian at boundary - - BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); - - BoutReal gradient = 2. * (f[i] - fdown[iym]) / (coord->dy[i] + coord->dy[iym]); - - BoutReal flux = c * J * gradient / g_22; - - result[i] -= flux / (coord->dy[i] * coord->J[i]); - flow_ylow[i] = - flux * coord->dx[i] * coord->dz[i]; - } - } - - if (!use_parallel_slices) { - // Shifted to field aligned coordinates, so need to shift back - result = fromFieldAligned(result, "RGN_NOBNDRY"); - flow_ylow = fromFieldAligned(flow_ylow); - } - - return result; -} From d98655966870b5a6f67f3aacc9bd9d45a0a6fd9c Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 22 Jun 2026 10:48:13 +0200 Subject: [PATCH 09/24] FV::Superbee has been moved to BOUT++ --- include/div_ops.hxx | 47 --------------------------------------------- 1 file changed, 47 deletions(-) diff --git a/include/div_ops.hxx b/include/div_ops.hxx index b24ff28bc..4e00c1536 100644 --- a/include/div_ops.hxx +++ b/include/div_ops.hxx @@ -93,53 +93,6 @@ Field3D Div_a_Grad_perp_nonorthog(const Field3D& a, const Field3D& x, Field3D& f Field3D& flux_ylow); namespace FV { - -/// Superbee limiter -/// -/// This corresponds to the limiter function -/// φ(r) = max(0, min(2r, 1), min(r,2) -/// -/// The value at cell right (i.e. i + 1/2) is: -/// -/// n.R = n.c - φ(r) (n.c - (n.p + n.c)/2) -/// = n.c + φ(r) (n.p - n.c)/2 -/// -/// Four regimes: -/// a) r < 1/2 -> φ(r) = 2r -/// n.R = n.c + gL -/// b) 1/2 < r < 1 -> φ(r) = 1 -/// n.R = n.c + gR/2 -/// c) 1 < r < 2 -> φ(r) = r -/// n.R = n.c + gL/2 -/// d) 2 < r -> φ(r) = 2 -/// n.R = n.c + gR -/// -/// where the left and right gradients are: -/// gL = n.c - n.m -/// gR = n.p - n.c -/// -struct Superbee { - void operator()(Stencil1D& n) { - BoutReal gL = n.c - n.L; - BoutReal gR = n.R - n.c; - - // r = gL / gR - // Limiter is φ(r) - if (gL * gR < 0) { - // Different signs => Zero gradient - n.L = n.R = n.c; - } else { - const BoutReal sign = SIGN(gL); - gL = fabs(gL); - gR = fabs(gR); - const BoutReal half_slope = - sign * BOUTMAX(BOUTMIN(gL, 0.5 * gR), BOUTMIN(gR, 0.5 * gL)); - n.L = n.c - half_slope; - n.R = n.c + half_slope; - } - } -}; - // Calculates viscous heating due to numerical momentum fluxes // and flow of kinetic energy (in flow_ylow) template From 52ec710da8563459ff4704cbab06cfff28eaf51f Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 23 Jun 2026 13:51:23 +0200 Subject: [PATCH 10/24] CI: Add 3D metric build --- .github/workflows/tests.yml | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index c181f0756..a9bc176f5 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -31,14 +31,19 @@ jobs: config: - name: "CMake default options" os: ubuntu-22.04 - configure_options: "-DPACKAGE_TESTS=OFF -DHERMES_ERROR_ON_WARNINGS=ON" build_type: - cmake_build_type: Debug test_type: unit coverage: ON + configure_options: "-DPACKAGE_TESTS=OFF -DHERMES_ERROR_ON_WARNINGS=ON" - cmake_build_type: Release test_type: integration coverage: OFF + configure_options: "-DPACKAGE_TESTS=OFF -DHERMES_ERROR_ON_WARNINGS=ON" + - cmake_build_type: Debug + test_type: None + configure_options: "-PACKAGE_TESTS=OFF -DHERMES_ERROR_ON_WARNINGS=ON -DBOUT_ENABLE_METRIC_3D=ON" + coverage: OFF steps: - name: Job information From 52908237920be66b789651dfb7d15263dc18d16a Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 22 Jun 2026 11:00:34 +0200 Subject: [PATCH 11/24] Fix div_ops.hxx for 3D compilation --- include/div_ops.hxx | 221 ++++++++++++++++++++++++-------------------- 1 file changed, 122 insertions(+), 99 deletions(-) diff --git a/include/div_ops.hxx b/include/div_ops.hxx index 4e00c1536..0b648d7b8 100644 --- a/include/div_ops.hxx +++ b/include/div_ops.hxx @@ -419,12 +419,12 @@ const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, // Flux across cell edge const BoutReal fout = gradient * aedge * gedge - * (coord->J(i, j) * coord->g11(i, j) - + coord->J(i + 1, j) * coord->g11(i + 1, j)) - / (coord->dx(i, j) + coord->dx(i + 1, j)); + * (coord->J(i, j, k) * coord->g11(i, j, k) + + coord->J(i + 1, j, k) * coord->g11(i + 1, j, k)) + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); - result(i, j, k) += fout / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= fout / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + result(i, j, k) += fout / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= fout / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); } } } @@ -471,119 +471,142 @@ const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, for (int i = mesh->xstart; i <= mesh->xend; i++) { for (int j = mesh->ystart; j <= mesh->yend; j++) { - - BoutReal coef_u = - 0.5 - * (coord->g_23(i, j) / SQ(coord->J(i, j) * coord->Bxy(i, j)) - + coord->g_23(i, j + 1) / SQ(coord->J(i, j + 1) * coord->Bxy(i, j + 1))); - - BoutReal coef_d = - 0.5 - * (coord->g_23(i, j) / SQ(coord->J(i, j) * coord->Bxy(i, j)) - + coord->g_23(i, j - 1) / SQ(coord->J(i, j - 1) * coord->Bxy(i, j - 1))); - +#if BOUT_USE_METRIC_3D for (int k = 0; k < mesh->LocalNz; k++) { - // Calculate flux between j and j+1 - int kp = (k + 1) % mesh->LocalNz; - int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; - - // Calculate Z derivative at y boundary - BoutReal dfdz = - 0.25 * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km)) - / coord->dz(i, j); +#else + int k = 0; +#endif + BoutReal coef_u = + 0.5 + * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)) + + coord->g_23(i, j + 1, k) + / SQ(coord->J(i, j + 1, k) * coord->Bxy(i, j + 1, k))); + + BoutReal coef_d = + 0.5 + * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)) + + coord->g_23(i, j - 1, k) + / SQ(coord->J(i, j - 1, k) * coord->Bxy(i, j - 1, k))); + +#if not BOUT_USE_METRIC_3D + for (int k = 0; k < mesh->LocalNz; k++) { +#endif + // Calculate flux between j and j+1 + int kp = (k + 1) % mesh->LocalNz; + int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; + + // Calculate Z derivative at y boundary + BoutReal dfdz = + 0.25 * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km)) + / coord->dz(i, j, k); + + // Y derivative + BoutReal dfdy = 2. * (fup(i, j + 1, k) - fc(i, j, k)) + / (coord->dy(i, j + 1, k) + coord->dy(i, j, k)); + + BoutReal aedge = 0.5 * (ac(i, j, k) + aup(i, j + 1, k)); + BoutReal gedge; + if ((j == mesh->yend) and mesh->lastY(i)) { + // Midpoint boundary value + gedge = 0.5 * (gc(i, j, k) + gup(i, j + 1, k)); + } else if (dfdy > 0) { + // Flux from (j+1) to (j) + gedge = gup(i, j + 1, k); + } else { + // Flux from (j) to (j+1) + gedge = gc(i, j, k); + } - // Y derivative - BoutReal dfdy = 2. * (fup(i, j + 1, k) - fc(i, j, k)) - / (coord->dy(i, j + 1) + coord->dy(i, j)); + BoutReal fout = aedge * gedge * 0.5 + * (coord->J(i, j, k) * coord->g23(i, j, k) + + coord->J(i, j + 1, k) * coord->g23(i, j + 1, k)) + * (dfdz - coef_u * dfdy); - BoutReal aedge = 0.5 * (ac(i, j, k) + aup(i, j + 1, k)); - BoutReal gedge; - if ((j == mesh->yend) and mesh->lastY(i)) { - // Midpoint boundary value - gedge = 0.5 * (gc(i, j, k) + gup(i, j + 1, k)); - } else if (dfdy > 0) { - // Flux from (j+1) to (j) - gedge = gup(i, j + 1, k); - } else { - // Flux from (j) to (j+1) - gedge = gc(i, j, k); - } + yzresult(i, j, k) = fout / (coord->dy(i, j, k) * coord->J(i, j, k)); - BoutReal fout = aedge * gedge * 0.5 - * (coord->J(i, j) * coord->g23(i, j) - + coord->J(i, j + 1) * coord->g23(i, j + 1)) - * (dfdz - coef_u * dfdy); + // Calculate flux between j and j-1 + dfdz = + 0.25 + * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) - fdown(i, j - 1, km)) + / coord->dz(i, j, k); - yzresult(i, j, k) = fout / (coord->dy(i, j) * coord->J(i, j)); + dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k)) + / (coord->dy(i, j, k) + coord->dy(i, j - 1, k)); - // Calculate flux between j and j-1 - dfdz = 0.25 - * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) - fdown(i, j - 1, km)) - / coord->dz(i, j); + aedge = 0.5 * (ac(i, j, k) + adown(i, j - 1, k)); + if ((j == mesh->ystart) and mesh->firstY(i)) { + gedge = 0.5 * (gc(i, j, k) + gdown(i, j - 1, k)); + } else if (dfdy > 0) { + gedge = gc(i, j, k); + } else { + gedge = gdown(i, j - 1, k); + } - dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k)) - / (coord->dy(i, j) + coord->dy(i, j - 1)); + fout = aedge * gedge * 0.5 + * (coord->J(i, j, k) * coord->g23(i, j, k) + + coord->J(i, j - 1, k) * coord->g23(i, j - 1, k)) + * (dfdz - coef_d * dfdy); - aedge = 0.5 * (ac(i, j, k) + adown(i, j - 1, k)); - if ((j == mesh->ystart) and mesh->firstY(i)) { - gedge = 0.5 * (gc(i, j, k) + gdown(i, j - 1, k)); - } else if (dfdy > 0) { - gedge = gc(i, j, k); - } else { - gedge = gdown(i, j - 1, k); + yzresult(i, j, k) -= fout / (coord->dy(i, j, k) * coord->J(i, j, k)); +#if BOUT_USE_METRIC_3D } - - fout = aedge * gedge * 0.5 - * (coord->J(i, j) * coord->g23(i, j) - + coord->J(i, j - 1) * coord->g23(i, j - 1)) - * (dfdz - coef_d * dfdy); - - yzresult(i, j, k) -= fout / (coord->dy(i, j) * coord->J(i, j)); +#else + } +#endif } } - } - // Z flux - // Easier since all metrics constant in Z + // Z flux + // Easier since all metrics constant in Z - for (int i = mesh->xstart; i <= mesh->xend; i++) { - for (int j = mesh->ystart; j <= mesh->yend; j++) { - // Coefficient in front of df/dy term - BoutReal coef = coord->g_23(i, j) - / (coord->dy(i, j + 1) + 2. * coord->dy(i, j) + coord->dy(i, j - 1)) - / SQ(coord->J(i, j) * coord->Bxy(i, j)); - - for (int k = 0; k < mesh->LocalNz; k++) { - // Calculate flux between k and k+1 - int kp = (k + 1) % mesh->LocalNz; + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { +#if BOUT_USE_METRIC_3D + for (int k = 0; k < mesh->LocalNz; k++) { +#else + int k = 0; +#endif + // Coefficient in front of df/dy term + BoutReal coef = coord->g_23(i, j, k) + / (coord->dy(i, j + 1, k) + 2. * coord->dy(i, j, k) + + coord->dy(i, j - 1, k)) + / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)); +#if not BOUT_USE_METRIC_3D + for (int k = 0; k < mesh->LocalNz; k++) { +#endif + // Calculate flux between k and k+1 + int kp = (k + 1) % mesh->LocalNz; - BoutReal gradient = - // df/dz - (fc(i, j, kp) - fc(i, j, k)) / coord->dz(i, j) + BoutReal gradient = + // df/dz + (fc(i, j, kp) - fc(i, j, k)) / coord->dz(i, j, k) - // - g_yz * df/dy / SQ(J*B) - - coef - * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k) - - fdown(i, j - 1, kp)); + // - g_yz * df/dy / SQ(J*B) + - coef + * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k) + - fdown(i, j - 1, kp)); - BoutReal fout = gradient * 0.5 * (ac(i, j, kp) + ac(i, j, k)) - * ((gradient > 0) ? gc(i, j, kp) : gc(i, j, k)); + BoutReal fout = gradient * 0.5 * (ac(i, j, kp) + ac(i, j, k)) + * ((gradient > 0) ? gc(i, j, kp) : gc(i, j, k)); - yzresult(i, j, k) += fout / coord->dz(i, j); - yzresult(i, j, kp) -= fout / coord->dz(i, j); + yzresult(i, j, k) += fout / coord->dz(i, j, k); + yzresult(i, j, kp) -= fout / coord->dz(i, j, kp); +#if BOUT_USE_METRIC_3D + } +#else + } +#endif + } } + // Check if we need to transform back + if (f.hasParallelSlices() && a.hasParallelSlices()) { + result += yzresult; + } else { + result += fromFieldAligned(yzresult); + } + return result; } - } - // Check if we need to transform back - if (f.hasParallelSlices() && a.hasParallelSlices()) { - result += yzresult; - } else { - result += fromFieldAligned(yzresult); - } - - return result; -} -} // namespace FV + } // namespace FV #endif // DIV_OPS_H From ebc484e69b95d1431520129a15c509f15a518677 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 23 Jun 2026 13:53:31 +0200 Subject: [PATCH 12/24] Fix div_ops for 3D compilation --- include/div_ops.hxx | 87 ++- src/div_ops.cxx | 1680 +++++++++++++++++++++++-------------------- 2 files changed, 959 insertions(+), 808 deletions(-) diff --git a/include/div_ops.hxx b/include/div_ops.hxx index 0b648d7b8..eec7e6523 100644 --- a/include/div_ops.hxx +++ b/include/div_ops.hxx @@ -93,6 +93,7 @@ Field3D Div_a_Grad_perp_nonorthog(const Field3D& a, const Field3D& x, Field3D& f Field3D& flux_ylow); namespace FV { + // Calculates viscous heating due to numerical momentum fluxes // and flow of kinetic energy (in flow_ylow) template @@ -472,10 +473,11 @@ const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, for (int i = mesh->xstart; i <= mesh->xend; i++) { for (int j = mesh->ystart; j <= mesh->yend; j++) { #if BOUT_USE_METRIC_3D - for (int k = 0; k < mesh->LocalNz; k++) { + for (int k = 0; k < mesh->LocalNz; k++) #else int k = 0; #endif + { BoutReal coef_u = 0.5 * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)) @@ -489,8 +491,9 @@ const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, / SQ(coord->J(i, j - 1, k) * coord->Bxy(i, j - 1, k))); #if not BOUT_USE_METRIC_3D - for (int k = 0; k < mesh->LocalNz; k++) { + for (int k = 0; k < mesh->LocalNz; k++) #endif + { // Calculate flux between j and j+1 int kp = (k + 1) % mesh->LocalNz; int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; @@ -548,65 +551,61 @@ const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, * (dfdz - coef_d * dfdy); yzresult(i, j, k) -= fout / (coord->dy(i, j, k) * coord->J(i, j, k)); -#if BOUT_USE_METRIC_3D } -#else - } -#endif } } + } - // Z flux - // Easier since all metrics constant in Z + // Z flux + // Easier since all metrics constant in Z - for (int i = mesh->xstart; i <= mesh->xend; i++) { - for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { #if BOUT_USE_METRIC_3D - for (int k = 0; k < mesh->LocalNz; k++) { + for (int k = 0; k < mesh->LocalNz; k++) #else - int k = 0; + int k = 0; #endif - // Coefficient in front of df/dy term - BoutReal coef = coord->g_23(i, j, k) - / (coord->dy(i, j + 1, k) + 2. * coord->dy(i, j, k) - + coord->dy(i, j - 1, k)) - / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)); + { + // Coefficient in front of df/dy term + BoutReal coef = + coord->g_23(i, j, k) + / (coord->dy(i, j + 1, k) + 2. * coord->dy(i, j, k) + coord->dy(i, j - 1, k)) + / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)); #if not BOUT_USE_METRIC_3D - for (int k = 0; k < mesh->LocalNz; k++) { + for (int k = 0; k < mesh->LocalNz; k++) #endif - // Calculate flux between k and k+1 - int kp = (k + 1) % mesh->LocalNz; + { + // Calculate flux between k and k+1 + int kp = (k + 1) % mesh->LocalNz; - BoutReal gradient = - // df/dz - (fc(i, j, kp) - fc(i, j, k)) / coord->dz(i, j, k) + BoutReal gradient = + // df/dz + (fc(i, j, kp) - fc(i, j, k)) / coord->dz(i, j, k) - // - g_yz * df/dy / SQ(J*B) - - coef - * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k) - - fdown(i, j - 1, kp)); + // - g_yz * df/dy / SQ(J*B) + - coef + * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k) + - fdown(i, j - 1, kp)); - BoutReal fout = gradient * 0.5 * (ac(i, j, kp) + ac(i, j, k)) - * ((gradient > 0) ? gc(i, j, kp) : gc(i, j, k)); + BoutReal fout = gradient * 0.5 * (ac(i, j, kp) + ac(i, j, k)) + * ((gradient > 0) ? gc(i, j, kp) : gc(i, j, k)); - yzresult(i, j, k) += fout / coord->dz(i, j, k); - yzresult(i, j, kp) -= fout / coord->dz(i, j, kp); -#if BOUT_USE_METRIC_3D - } -#else - } -#endif + yzresult(i, j, k) += fout / coord->dz(i, j, k); + yzresult(i, j, kp) -= fout / coord->dz(i, j, kp); } } - // Check if we need to transform back - if (f.hasParallelSlices() && a.hasParallelSlices()) { - result += yzresult; - } else { - result += fromFieldAligned(yzresult); - } - return result; } + } + // Check if we need to transform back + if (f.hasParallelSlices() && a.hasParallelSlices()) { + result += yzresult; + } else { + result += fromFieldAligned(yzresult); + } + return result; +} - } // namespace FV +} // namespace FV #endif // DIV_OPS_H diff --git a/src/div_ops.cxx b/src/div_ops.cxx index d83c023d7..f28f99a5f 100644 --- a/src/div_ops.cxx +++ b/src/div_ops.cxx @@ -56,14 +56,15 @@ const Field3D Div_par_diffusion_index(const Field3D& f, bool bndry_flux) { continue; } } - BoutReal J = 0.5 * (coord->J(i, j) + coord->J(i, j + 1)); // Jacobian at boundary + BoutReal J = + 0.5 * (coord->J(i, j, k) + coord->J(i, j + 1, k)); // Jacobian at boundary BoutReal gradient = f(i, j + 1, k) - f(i, j, k); BoutReal flux = J * gradient; - result(i, j, k) += flux / coord->J(i, j); - result(i, j + 1, k) -= flux / coord->J(i, j + 1); + result(i, j, k) += flux / coord->J(i, j, k); + result(i, j + 1, k) -= flux / coord->J(i, j + 1, k); } } } @@ -151,13 +152,13 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd // 2) Calculate velocities on cell faces - BoutReal vU = coord->J(i, j) * (fmp - fpp) / coord->dx(i, j); // -J*df/dx - BoutReal vD = coord->J(i, j) * (fmm - fpm) / coord->dx(i, j); // -J*df/dx + BoutReal vU = coord->J(i, j, k) * (fmp - fpp) / coord->dx(i, j, k); // -J*df/dx + BoutReal vD = coord->J(i, j, k) * (fmm - fpm) / coord->dx(i, j, k); // -J*df/dx - BoutReal vR = 0.5 * (coord->J(i, j) + coord->J(i + 1, j)) * (fpp - fpm) - / coord->dz(i, j); // J*df/dz - BoutReal vL = 0.5 * (coord->J(i, j) + coord->J(i - 1, j)) * (fmp - fmm) - / coord->dz(i, j); // J*df/dz + BoutReal vR = 0.5 * (coord->J(i, j, k) + coord->J(i + 1, j, k)) * (fpp - fpm) + / coord->dz(i, j, k); // J*df/dz + BoutReal vL = 0.5 * (coord->J(i, j, k) + coord->J(i - 1, j, k)) * (fmp - fmm) + / coord->dz(i, j, k); // J*df/dz // 3) Calculate n on the cell faces. The sign of the // velocity determines which side is used. @@ -185,16 +186,18 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd // Flux in from boundary flux = vR * 0.5 * (n(i + 1, j, k) + n(i, j, k)); } - result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + result(i, j, k) += flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= + flux / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); } } else { // Not at a boundary if (vR > 0.0) { // Flux out into next cell BoutReal flux = vR * s.R; - result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + result(i, j, k) += flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= + flux / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); // if(i==mesh->xend) // output.write("Setting flux (%d,%d) : %e\n", @@ -218,16 +221,18 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd // Flux in from boundary flux = vL * 0.5 * (n(i - 1, j, k) + n(i, j, k)); } - result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); - result(i - 1, j, k) += flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); + result(i, j, k) -= flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i - 1, j, k) += + flux / (coord->dx(i - 1, j, k) * coord->J(i - 1, j, k)); } } else { // Not at a boundary if (vL < 0.0) { BoutReal flux = vL * s.L; - result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); - result(i - 1, j, k) += flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); + result(i, j, k) -= flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i - 1, j, k) += + flux / (coord->dx(i - 1, j, k) * coord->J(i - 1, j, k)); } } @@ -242,12 +247,12 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd MC(s); if (vU > 0.0) { - BoutReal flux = vU * s.R / (coord->J(i, j) * coord->dz(i, j)); + BoutReal flux = vU * s.R / (coord->J(i, j, k) * coord->dz(i, j, k)); result(i, j, k) += flux; result(i, j, kp) -= flux; } if (vD < 0.0) { - BoutReal flux = vD * s.L / (coord->J(i, j) * coord->dz(i, j)); + BoutReal flux = vD * s.L / (coord->J(i, j, k) * coord->dz(i, j, k)); result(i, j, k) -= flux; result(i, j, km) += flux; } @@ -300,13 +305,14 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd // Average dfdy to right X boundary BoutReal f_R = 0.5 - * ((coord->g11(i + 1, j) * coord->g23(i + 1, j) / SQ(coord->Bxy(i + 1, j))) + * ((coord->g11(i + 1, j, k) * coord->g23(i + 1, j, k) + / SQ(coord->Bxy(i + 1, j, k))) * dfdy(i + 1, j, k) - + (coord->g11(i, j) * coord->g23(i, j) / SQ(coord->Bxy(i, j))) + + (coord->g11(i, j, k) * coord->g23(i, j, k) / SQ(coord->Bxy(i, j, k))) * dfdy(i, j, k)); // Advection velocity across cell face - BoutReal Vx = 0.5 * (coord->J(i + 1, j) + coord->J(i, j)) * f_R; + BoutReal Vx = 0.5 * (coord->J(i + 1, j, k) + coord->J(i, j, k)) * f_R; // Fromm method BoutReal flux = Vx; @@ -327,8 +333,8 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd flux *= nval; } - result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + result(i, j, k) += flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= flux / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); } } } @@ -370,12 +376,13 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd // Average dfdx to upper Y boundary BoutReal f_U = 0.5 - * ((coord->g11(i, j + 1) * coord->g23(i, j + 1) / SQ(coord->Bxy(i, j + 1))) + * ((coord->g11(i, j + 1, k) * coord->g23(i, j + 1, k) + / SQ(coord->Bxy(i, j + 1, k))) * dfdx(i, j + 1, k) - + (coord->g11(i, j) * coord->g23(i, j) / SQ(coord->Bxy(i, j))) + + (coord->g11(i, j, k) * coord->g23(i, j, k) / SQ(coord->Bxy(i, j, k))) * dfdx(i, j, k)); - BoutReal Vy = -0.5 * (coord->J(i, j + 1) + coord->J(i, j)) * f_U; + BoutReal Vy = -0.5 * (coord->J(i, j + 1, k) + coord->J(i, j, k)) * f_U; if (mesh->firstY(i) && !mesh->periodicY(i) && (j == mesh->ystart - 1)) { // Lower y boundary. Allow flows out of the domain only @@ -406,8 +413,8 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd flux *= nval; } - yresult(i, j, k) += flux / (coord->dy(i, j) * coord->J(i, j)); - yresult(i, j + 1, k) -= flux / (coord->dy(i, j + 1) * coord->J(i, j + 1)); + yresult(i, j, k) += flux / (coord->dy(i, j, k) * coord->J(i, j, k)); + yresult(i, j + 1, k) -= flux / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); } } } @@ -455,18 +462,18 @@ const Field3D Div_Perp_Lap_FV_Index(const Field3D& as, const Field3D& fs) { BoutReal gU = fs(i, j, kp) - fs(i, j, k); // Flow right - BoutReal flux = gR * 0.25 * (coord->J(i + 1, j) + coord->J(i, j)) - * (coord->dx(i + 1, j) + coord->dx(i, j)) + BoutReal flux = gR * 0.25 * (coord->J(i + 1, j, k) + coord->J(i, j, k)) + * (coord->dx(i + 1, j, k) + coord->dx(i, j, k)) * (as(i + 1, j, k) + as(i, j, k)); - result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); + result(i, j, k) += flux / (coord->dx(i, j, k) * coord->J(i, j, k)); // Flow left - flux = gL * 0.25 * (coord->J(i - 1, j) + coord->J(i, j)) - * (coord->dx(i - 1, j) + coord->dx(i, j)) + flux = gL * 0.25 * (coord->J(i - 1, j, k) + coord->J(i, j, k)) + * (coord->dx(i - 1, j, k) + coord->dx(i, j, k)) * (as(i - 1, j, k) + as(i, j, k)); - result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); + result(i, j, k) -= flux / (coord->dx(i, j, k) * coord->J(i, j, k)); // Flow up @@ -524,8 +531,8 @@ const Field3D D4DX4_FV_Index(const Field3D& f, bool bndry_flux) { BoutReal d3fdx3 = (f(i + 2, j, k) - 3. * f(i + 1, j, k) + 3. * f(i, j, k) - f(i - 1, j, k)); - BoutReal flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) - * (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; + BoutReal flux = 0.25 * (coord->dx(i, j, k) + coord->dx(i + 1, j, k)) + * (coord->J(i, j, k) + coord->J(i + 1, j, k)) * d3fdx3; if (mesh->lastX() && (i == mesh->xend)) { // Boundary @@ -540,8 +547,8 @@ const Field3D D4DX4_FV_Index(const Field3D& f, bool bndry_flux) { - (6. / 5) * f(i - 2, j, k) // f_2 ); - flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) - * (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; + flux = 0.25 * (coord->dx(i, j, k) + coord->dx(i + 1, j, k)) + * (coord->J(i, j, k) + coord->J(i + 1, j, k)) * d3fdx3; } else { // No fluxes through boundary @@ -549,8 +556,8 @@ const Field3D D4DX4_FV_Index(const Field3D& f, bool bndry_flux) { } } - result(i, j, k) += flux / (coord->J(i, j) * coord->dx(i, j)); - result(i + 1, j, k) -= flux / (coord->J(i + 1, j) * coord->dx(i + 1, j)); + result(i, j, k) += flux / (coord->J(i, j, k) * coord->dx(i, j, k)); + result(i + 1, j, k) -= flux / (coord->J(i + 1, j, k) * coord->dx(i + 1, j, k)); if (j == mesh->xstart) { // Left cell boundary, no flux through boundaries @@ -566,11 +573,12 @@ const Field3D D4DX4_FV_Index(const Field3D& f, bool bndry_flux) { + (6. / 5) * f(i + 2, j, k) // f_2 ); - flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) - * (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; + flux = 0.25 * (coord->dx(i, j, k) + coord->dx(i + 1, j, k)) + * (coord->J(i, j, k) + coord->J(i + 1, j, k)) * d3fdx3; - result(i, j, k) -= flux / (coord->J(i, j) * coord->dx(i, j)); - result(i - 1, j, k) += flux / (coord->J(i - 1, j) * coord->dx(i - 1, j)); + result(i, j, k) -= flux / (coord->J(i, j, k) * coord->dx(i, j, k)); + result(i - 1, j, k) += + flux / (coord->J(i - 1, j, k) * coord->dx(i - 1, j, k)); } } else { @@ -578,11 +586,12 @@ const Field3D D4DX4_FV_Index(const Field3D& f, bool bndry_flux) { d3fdx3 = (f(i + 1, j, k) - 3. * f(i, j, k) + 3. * f(i - 1, j, k) - f(i - 2, j, k)); - flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) - * (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; + flux = 0.25 * (coord->dx(i, j, k) + coord->dx(i + 1, j, k)) + * (coord->J(i, j, k) + coord->J(i + 1, j, k)) * d3fdx3; - result(i, j, k) -= flux / (coord->J(i, j) * coord->dx(i, j)); - result(i - 1, j, k) += flux / (coord->J(i - 1, j) * coord->dx(i - 1, j)); + result(i, j, k) -= flux / (coord->J(i, j, k) * coord->dx(i, j, k)); + result(i - 1, j, k) += + flux / (coord->J(i - 1, j, k) * coord->dx(i - 1, j, k)); } } } @@ -608,6 +617,12 @@ const Field3D D4DZ4_Index(const Field3D& f) { * we would need the corner cell values to take Y derivatives along X edges * */ +#if BOUT_USE_METRIC_3D +const Field2D Laplace_FV([[maybe_unused]] const Field2D& k, + [[maybe_unused]] const Field2D& f) { + throw BoutException("Not Implemented!"); +} +#else const Field2D Laplace_FV(const Field2D& k, const Field2D& f) { Field2D result; result.allocate(); @@ -619,43 +634,51 @@ const Field2D Laplace_FV(const Field2D& k, const Field2D& f) { // Calculate gradients on cell faces - BoutReal gR = (coord->g11(i, j) + coord->g11(i + 1, j)) * (f(i + 1, j) - f(i, j)) - / (coord->dx(i + 1, j) + coord->dx(i, j)); + BoutReal gR = (coord->g11(i, j, k) + coord->g11(i + 1, j, k)) + * (f(i + 1, j, k) - f(i, j, k)) + / (coord->dx(i + 1, j, k) + coord->dx(i, j, k)); - BoutReal gL = (coord->g11(i - 1, j) + coord->g11(i, j)) * (f(i, j) - f(i - 1, j)) - / (coord->dx(i - 1, j) + coord->dx(i, j)); + BoutReal gL = (coord->g11(i - 1, j, k) + coord->g11(i, j, k)) + * (f(i, j, k) - f(i - 1, j, k)) + / (coord->dx(i - 1, j, k) + coord->dx(i, j, k)); - BoutReal gU = (coord->g22(i, j) + coord->g22(i, j + 1)) * (f(i, j + 1) - f(i, j)) - / (coord->dy(i, j + 1) + coord->dy(i, j)); + BoutReal gU = (coord->g22(i, j, k) + coord->g22(i, j + 1)) + * (f(i, j + 1) - f(i, j, k)) + / (coord->dy(i, j + 1) + coord->dy(i, j, k)); - BoutReal gD = (coord->g22(i, j - 1) + coord->g22(i, j)) * (f(i, j) - f(i, j - 1)) - / (coord->dy(i, j) + coord->dy(i, j - 1)); + BoutReal gD = (coord->g22(i, j - 1) + coord->g22(i, j, k)) + * (f(i, j, k) - f(i, j - 1)) + / (coord->dy(i, j, k) + coord->dy(i, j - 1)); // Flow right - BoutReal flux = - gR * 0.25 * (coord->J(i + 1, j) + coord->J(i, j)) * (k(i + 1, j) + k(i, j)); + BoutReal flux = gR * 0.25 * (coord->J(i + 1, j, k) + coord->J(i, j, k)) + * (k(i + 1, j, k) + k(i, j, k)); - result(i, j) = flux / (coord->dx(i, j) * coord->J(i, j)); + result(i, j, k) = flux / (coord->dx(i, j, k) * coord->J(i, j, k)); // Flow left - flux = gL * 0.25 * (coord->J(i - 1, j) + coord->J(i, j)) * (k(i - 1, j) + k(i, j)); - result(i, j) -= flux / (coord->dx(i, j) * coord->J(i, j)); + flux = gL * 0.25 * (coord->J(i - 1, j, k) + coord->J(i, j, k)) + * (k(i - 1, j, k) + k(i, j, k)); + result(i, j, k) -= flux / (coord->dx(i, j, k) * coord->J(i, j, k)); // Flow up - flux = gU * 0.25 * (coord->J(i, j + 1) + coord->J(i, j)) * (k(i, j + 1) + k(i, j)); - result(i, j) += flux / (coord->dy(i, j) * coord->J(i, j)); + flux = gU * 0.25 * (coord->J(i, j + 1) + coord->J(i, j, k)) + * (k(i, j + 1) + k(i, j, k)); + result(i, j, k) += flux / (coord->dy(i, j, k) * coord->J(i, j, k)); // Flow down - flux = gD * 0.25 * (coord->J(i, j - 1) + coord->J(i, j)) * (k(i, j - 1) + k(i, j)); - result(i, j) -= flux / (coord->dy(i, j) * coord->J(i, j)); + flux = gD * 0.25 * (coord->J(i, j - 1) + coord->J(i, j, k)) + * (k(i, j - 1) + k(i, j, k)); + result(i, j, k) -= flux / (coord->dy(i, j, k) * coord->J(i, j, k)); } } return result; } +#endif const Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, Field3D& flow_xlow, Field3D& flow_ylow) { @@ -700,7 +723,7 @@ const Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, result(i + 1, j, k) -= fout / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); // Flow will be positive in the positive coordinate direction - flow_xlow(i + 1, j, k) = -1.0 * fout * coord->dy(i, j) * coord->dz(i, j); + flow_xlow(i + 1, j, k) = -1.0 * fout * coord->dy(i, j, k) * coord->dz(i, j, k); } } } @@ -842,7 +865,7 @@ const Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, yzresult(i, j, k) -= fout / (dyc(i, j, k) * Jc(i, j, k)); // Flow will be positive in the positive coordinate direction - flow_ylow(i, j, k) = -1.0 * fout * coord->dx(i, j) * coord->dz(i, j); + flow_ylow(i, j, k) = -1.0 * fout * coord->dx(i, j, k) * coord->dz(i, j, k); } } } @@ -907,16 +930,16 @@ const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { for (int k = 0; k < mesh->LocalNz; k++) { // Calculate flux from i to i+1 - const BoutReal gradient = (coord->J(i, j) * coord->g11(i, j) - + coord->J(i + 1, j) * coord->g11(i + 1, j)) + const BoutReal gradient = (coord->J(i, j, k) * coord->g11(i, j, k) + + coord->J(i + 1, j, k) * coord->g11(i + 1, j, k)) * (f(i + 1, j, k) - f(i, j, k)) - / (coord->dx(i, j) + coord->dx(i + 1, j)); + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); // Use the upwind coefficient const BoutReal fout = gradient * ((gradient > 0) ? a(i + 1, j, k) : a(i, j, k)); - result(i, j, k) += fout / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= fout / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + result(i, j, k) += fout / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= fout / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); } } } @@ -958,793 +981,922 @@ const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { for (int i = mesh->xstart; i <= mesh->xend; i++) { for (int j = mesh->ystart; j <= mesh->yend; j++) { - BoutReal coef_u = - 0.5 - * (coord->g_23(i, j) / SQ(coord->J(i, j) * coord->Bxy(i, j)) - + coord->g_23(i, j + 1) / SQ(coord->J(i, j + 1) * coord->Bxy(i, j + 1))); - - BoutReal coef_d = - 0.5 - * (coord->g_23(i, j) / SQ(coord->J(i, j) * coord->Bxy(i, j)) - + coord->g_23(i, j - 1) / SQ(coord->J(i, j - 1) * coord->Bxy(i, j - 1))); - +#if BOUT_USE_METRIC_3D for (int k = 0; k < mesh->LocalNz; k++) { - // Calculate flux between j and j+1 - int kp = (k + 1) % mesh->LocalNz; - int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; - - // Calculate Z derivative at y boundary - BoutReal dfdz = - 0.25 * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km)) - / coord->dz(i, j); - - // Y derivative - BoutReal dfdy = 2. * (fup(i, j + 1, k) - fc(i, j, k)) - / (coord->dy(i, j + 1) + coord->dy(i, j)); - - BoutReal fout = 0.25 * (ac(i, j, k) + aup(i, j + 1, k)) - * (coord->J(i, j) * coord->g23(i, j) - + coord->J(i, j + 1) * coord->g23(i, j + 1)) - * (dfdz - coef_u * dfdy); - - yzresult(i, j, k) = fout / (coord->dy(i, j) * coord->J(i, j)); - - // Calculate flux between j and j-1 - dfdz = 0.25 - * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) - fdown(i, j - 1, km)) - / coord->dz(i, j); - - dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k)) - / (coord->dy(i, j) + coord->dy(i, j - 1)); +#else + int k = 0; +#endif + BoutReal coef_u = + 0.5 + * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)) + + coord->g_23(i, j + 1, k) + / SQ(coord->J(i, j + 1, k) * coord->Bxy(i, j + 1, k))); - fout = 0.25 * (ac(i, j, k) + adown(i, j - 1, k)) - * (coord->J(i, j) * coord->g23(i, j) - + coord->J(i, j - 1) * coord->g23(i, j - 1)) - * (dfdz - coef_d * dfdy); + BoutReal coef_d = + 0.5 + * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)) + + coord->g_23(i, j - 1, k) + / SQ(coord->J(i, j - 1, k) * coord->Bxy(i, j - 1, k))); - yzresult(i, j, k) -= fout / (coord->dy(i, j) * coord->J(i, j)); +#if not BOUT_USE_METRIC_3D + for (int k = 0; k < mesh->LocalNz; k++) { +#endif + // Calculate flux between j and j+1 + int kp = (k + 1) % mesh->LocalNz; + int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; + + // Calculate Z derivative at y boundary + BoutReal dfdz = + 0.25 * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km)) + / coord->dz(i, j, k); + + // Y derivative + BoutReal dfdy = 2. * (fup(i, j + 1, k) - fc(i, j, k)) + / (coord->dy(i, j + 1, k) + coord->dy(i, j, k)); + + BoutReal fout = 0.25 * (ac(i, j, k) + aup(i, j + 1, k)) + * (coord->J(i, j, k) * coord->g23(i, j, k) + + coord->J(i, j + 1, k) * coord->g23(i, j + 1, k)) + * (dfdz - coef_u * dfdy); + + yzresult(i, j, k) = fout / (coord->dy(i, j, k) * coord->J(i, j, k)); + + // Calculate flux between j and j-1 + dfdz = + 0.25 + * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) - fdown(i, j - 1, km)) + / coord->dz(i, j, k); + + dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k)) + / (coord->dy(i, j, k) + coord->dy(i, j - 1, k)); + + fout = 0.25 * (ac(i, j, k) + adown(i, j - 1, k)) + * (coord->J(i, j, k) * coord->g23(i, j, k) + + coord->J(i, j - 1, k) * coord->g23(i, j - 1, k)) + * (dfdz - coef_d * dfdy); + + yzresult(i, j, k) -= fout / (coord->dy(i, j, k) * coord->J(i, j, k)); + } } } - } - - // Z flux - // Easier since all metrics constant in Z - - for (int i = mesh->xstart; i <= mesh->xend; i++) { - for (int j = mesh->ystart; j <= mesh->yend; j++) { - // Coefficient in front of df/dy term - BoutReal coef = coord->g_23(i, j) - / (coord->dy(i, j + 1) + 2. * coord->dy(i, j) + coord->dy(i, j - 1)) - / SQ(coord->J(i, j) * coord->Bxy(i, j)); - for (int k = 0; k < mesh->LocalNz; k++) { - // Calculate flux between k and k+1 - int kp = (k + 1) % mesh->LocalNz; + // Z flux + // Easier since all metrics constant in Z - BoutReal gradient = - // df/dz - (fc(i, j, kp) - fc(i, j, k)) / coord->dz(i, j) + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { +#if BOUT_USE_METRIC_3D + for (int k = 0; k < mesh->LocalNz; k++) { +#else + int k = 0; +#endif + // Coefficient in front of df/dy term + BoutReal coef = coord->g_23(i, j, k) + / (coord->dy(i, j + 1, k) + 2. * coord->dy(i, j, k) + + coord->dy(i, j - 1, k)) + / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)); + +#if not BOUT_USE_METRIC_3D + for (int k = 0; k < mesh->LocalNz; k++) { +#endif + // Calculate flux between k and k+1 + int kp = (k + 1) % mesh->LocalNz; + + BoutReal gradient = + // df/dz + (fc(i, j, kp) - fc(i, j, k)) / coord->dz(i, j, k) - // - g_yz * df/dy / SQ(J*B) - - coef - * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k) - - fdown(i, j - 1, kp)); + // - g_yz * df/dy / SQ(J*B) + - coef + * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k) + - fdown(i, j - 1, kp)); - BoutReal fout = gradient * ((gradient > 0) ? ac(i, j, kp) : ac(i, j, k)); + BoutReal fout = gradient * ((gradient > 0) ? ac(i, j, kp) : ac(i, j, k)); - yzresult(i, j, k) += fout / coord->dz(i, j); - yzresult(i, j, kp) -= fout / coord->dz(i, j); + yzresult(i, j, k) += fout / coord->dz(i, j, k); + yzresult(i, j, kp) -= fout / coord->dz(i, j, k); + } + } } - } - } - // Check if we need to transform back - if (f.hasParallelSlices() && a.hasParallelSlices()) { - result += yzresult; - } else { - result += fromFieldAligned(yzresult); - } - - return result; -} - -// Div ( a Grad_perp(f) ) -- ∇⊥ ( a ⋅ ∇⊥ f) -- Vorticity -// Includes nonorthogonal X-Y and X-Z metric corrections -// -Field3D Div_a_Grad_perp_nonorthog(const Field3D& a, const Field3D& f, Field3D& flow_xlow, - Field3D& flow_ylow) { - ASSERT2(a.getLocation() == f.getLocation()); - - Mesh* mesh = a.getMesh(); - - Coordinates* coord = f.getCoordinates(); - - // Zero all flows - flow_xlow = 0.0; - flow_ylow = 0.0; - - // Y and Z fluxes require Y derivatives - - // Fields containing values along the magnetic field - Field3D fup(mesh), fdown(mesh); - Field3D aup(mesh), adown(mesh); - - Field3D g23up(mesh), g23down(mesh); - Field3D g_23up(mesh), g_23down(mesh); - Field3D g12up(mesh), g12down(mesh); - Field3D g_12up(mesh), g_12down(mesh); - Field3D Jup(mesh), Jdown(mesh); - Field3D dyup(mesh), dydown(mesh); - Field3D dzup(mesh), dzdown(mesh); - Field3D Bxyup(mesh), Bxydown(mesh); - - // Values on this y slice (centre). - // This is needed because toFieldAligned may modify the field - Field3D fc = f; - Field3D ac = a; - - Field3D g23c = coord->g23; - Field3D g_23c = coord->g_23; - Field3D g12c = coord->g12; - Field3D g_12c = coord->g_12; - Field3D Jc = coord->J; - Field3D dxc = coord->dx; - Field3D dyc = coord->dy; - Field3D dzc = coord->dz; - Field3D Bxyc = coord->Bxy; - - // Calculate the X derivative at cell edge (X + 1/2), including in Y guard cells - // This is used to calculate Y flux contribution from g21 * d/dx - Field3D fddx_xhigh(mesh); - fddx_xhigh.allocate(); - for (int i = mesh->xstart - 1; i <= mesh->xend; i++) { - for (int j = mesh->ystart - 1; j <= mesh->yend + 1; - j++) { // Note: Including one guard cell - for (int k = 0; k < mesh->LocalNz; k++) { - fddx_xhigh(i, j, k) = 2. * (f(i + 1, j, k) - f(i, j, k)) - / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); + // Check if we need to transform back + if (f.hasParallelSlices() && a.hasParallelSlices()) { + result += yzresult; + } else { + result += fromFieldAligned(yzresult); } - } - } - Field3D fddx_xhigh_up(mesh), fddx_xhigh_down(mesh); - - // Result of the Y and Z fluxes - Field3D yzresult(mesh); - yzresult.allocate(); - - if (f.hasParallelSlices() && a.hasParallelSlices()) { - // Both inputs have yup and ydown - - fup = f.yup(); - fdown = f.ydown(); - - aup = a.yup(); - adown = a.ydown(); - - mesh->communicate(fddx_xhigh); - fddx_xhigh_up = fddx_xhigh.yup(); - fddx_xhigh_down = fddx_xhigh.ydown(); - } else { - // At least one input doesn't have yup/ydown fields. - // Need to shift to/from field aligned coordinates - - fup = fdown = fc = toFieldAligned(f); - aup = adown = ac = toFieldAligned(a); - - fddx_xhigh_up = fddx_xhigh_down = toFieldAligned(fddx_xhigh); - - yzresult.setDirectionY(YDirectionType::Aligned); - flow_ylow.setDirectionY(YDirectionType::Aligned); - } - if (bout::build::use_metric_3d) { - // 3D Metric, need yup/ydown fields. - // Requires previous communication of metrics - // -- should insert communication here? - if (!coord->g23.hasParallelSlices() || !coord->g_23.hasParallelSlices() - || !coord->dy.hasParallelSlices() || !coord->dz.hasParallelSlices() - || !coord->Bxy.hasParallelSlices() || !coord->J.hasParallelSlices()) { - throw BoutException("metrics have no yup/down"); + return result; } - g23up = coord->g23.yup(); - g23down = coord->g23.ydown(); + // Div ( a Grad_perp(f) ) -- ∇⊥ ( a ⋅ ∇⊥ f) -- Vorticity + // Includes nonorthogonal X-Y and X-Z metric corrections + // + Field3D Div_a_Grad_perp_nonorthog(const Field3D& a, const Field3D& f, + Field3D& flow_xlow, Field3D& flow_ylow) { + ASSERT2(a.getLocation() == f.getLocation()); + + Mesh* mesh = a.getMesh(); + + Coordinates* coord = f.getCoordinates(); + + // Zero all flows + flow_xlow = 0.0; + flow_ylow = 0.0; + + // Y and Z fluxes require Y derivatives + + // Fields containing values along the magnetic field + Field3D fup(mesh), fdown(mesh); + Field3D aup(mesh), adown(mesh); + + Field3D g23up(mesh), g23down(mesh); + Field3D g_23up(mesh), g_23down(mesh); + Field3D g12up(mesh), g12down(mesh); + Field3D g_12up(mesh), g_12down(mesh); + Field3D Jup(mesh), Jdown(mesh); + Field3D dyup(mesh), dydown(mesh); + Field3D dzup(mesh), dzdown(mesh); + Field3D Bxyup(mesh), Bxydown(mesh); + + // Values on this y slice (centre). + // This is needed because toFieldAligned may modify the field + Field3D fc = f; + Field3D ac = a; + + Field3D g23c = coord->g23; + Field3D g_23c = coord->g_23; + Field3D g12c = coord->g12; + Field3D g_12c = coord->g_12; + Field3D Jc = coord->J; + Field3D dxc = coord->dx; + Field3D dyc = coord->dy; + Field3D dzc = coord->dz; + Field3D Bxyc = coord->Bxy; + + // Calculate the X derivative at cell edge (X + 1/2), including in Y guard cells + // This is used to calculate Y flux contribution from g21 * d/dx + Field3D fddx_xhigh(mesh); + fddx_xhigh.allocate(); + for (int i = mesh->xstart - 1; i <= mesh->xend; i++) { + for (int j = mesh->ystart - 1; j <= mesh->yend + 1; + j++) { // Note: Including one guard cell + for (int k = 0; k < mesh->LocalNz; k++) { + fddx_xhigh(i, j, k) = 2. * (f(i + 1, j, k) - f(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); + } + } + } + Field3D fddx_xhigh_up(mesh), fddx_xhigh_down(mesh); - g_23up = coord->g_23.yup(); - g_23down = coord->g_23.ydown(); + // Result of the Y and Z fluxes + Field3D yzresult(mesh); + yzresult.allocate(); - g12up = coord->g12.yup(); - g12down = coord->g12.ydown(); + if (f.hasParallelSlices() && a.hasParallelSlices()) { + // Both inputs have yup and ydown - g_12up = coord->g_12.yup(); - g_12down = coord->g_12.ydown(); + fup = f.yup(); + fdown = f.ydown(); - Jup = coord->J.yup(); - Jdown = coord->J.ydown(); + aup = a.yup(); + adown = a.ydown(); - dyup = coord->dy.yup(); - dydown = coord->dy.ydown(); + mesh->communicate(fddx_xhigh); + fddx_xhigh_up = fddx_xhigh.yup(); + fddx_xhigh_down = fddx_xhigh.ydown(); + } else { + // At least one input doesn't have yup/ydown fields. + // Need to shift to/from field aligned coordinates - dzup = coord->dz.yup(); - dzdown = coord->dz.ydown(); + fup = fdown = fc = toFieldAligned(f); + aup = adown = ac = toFieldAligned(a); - Bxyup = coord->Bxy.yup(); - Bxydown = coord->Bxy.ydown(); + fddx_xhigh_up = fddx_xhigh_down = toFieldAligned(fddx_xhigh); - } else { - // No 3D metrics - // Need to shift to/from field aligned coordinates - g23up = g23down = g23c = toFieldAligned(coord->g23); - g_23up = g_23down = g_23c = toFieldAligned(coord->g_23); - g12up = g12down = g12c = toFieldAligned(coord->g12); - g_12up = g_12down = g_12c = toFieldAligned(coord->g_12); - Jup = Jdown = Jc = toFieldAligned(coord->J); - dxc = toFieldAligned(coord->dx); - dyup = dydown = dyc = toFieldAligned(coord->dy); - dzup = dzdown = dzc = toFieldAligned(coord->dz); - Bxyup = Bxydown = Bxyc = toFieldAligned(coord->Bxy); - } + yzresult.setDirectionY(YDirectionType::Aligned); + flow_ylow.setDirectionY(YDirectionType::Aligned); + } - // Y flux - // Includes fluxes due to Z derivatives (non-zero g23 metric) - // and due to X derivatives if grid is nonorthogonal (non-zero g12 metric) + if (bout::build::use_metric_3d) { + // 3D Metric, need yup/ydown fields. + // Requires previous communication of metrics + // -- should insert communication here? + if (!coord->g23.hasParallelSlices() || !coord->g_23.hasParallelSlices() + || !coord->dy.hasParallelSlices() || !coord->dz.hasParallelSlices() + || !coord->Bxy.hasParallelSlices() || !coord->J.hasParallelSlices()) { + throw BoutException("metrics have no yup/down"); + } - 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 between j and j+1 - int kp = (k + 1) % mesh->LocalNz; - int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; + g23up = coord->g23.yup(); + g23down = coord->g23.ydown(); - BoutReal coef_yz = - 0.5 - * (g_23c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) - + g_23up(i, j + 1, k) / SQ(Jup(i, j + 1, k) * Bxyup(i, j + 1, k))); + g_23up = coord->g_23.yup(); + g_23down = coord->g_23.ydown(); - BoutReal coef_xy = - 0.5 - * (g_12c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) - + g_12up(i, j + 1, k) / SQ(Jup(i, j + 1, k) * Bxyup(i, j + 1, k))); + g12up = coord->g12.yup(); + g12down = coord->g12.ydown(); - // Calculate Z derivative at y boundary - BoutReal dfdz = - 0.5 * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km)) - / (dzc(i, j, k) + dzup(i, j + 1, k)); + g_12up = coord->g_12.yup(); + g_12down = coord->g_12.ydown(); - // Y derivative - BoutReal dfdy = - 2. * (fup(i, j + 1, k) - fc(i, j, k)) / (dyup(i, j + 1, k) + dyc(i, j, k)); + Jup = coord->J.yup(); + Jdown = coord->J.ydown(); - // X derivative at Y boundary - BoutReal dfdx = 0.25 - * (fddx_xhigh(i - 1, j, k) + fddx_xhigh(i, j, k) - + fddx_xhigh_up(i - 1, j + 1, k) + fddx_xhigh_up(i, j + 1, k)); + dyup = coord->dy.yup(); + dydown = coord->dy.ydown(); - BoutReal fout = - 0.25 * (ac(i, j, k) + aup(i, j + 1, k)) - * ( - // Component of flux from d/dz and d/dy - (Jc(i, j, k) * g23c(i, j, k) + Jup(i, j + 1, k) * g23up(i, j + 1, k)) - * (dfdz - coef_yz * dfdy) - // Non-orthogonal mesh correction with g12 metric - + (Jc(i, j, k) * g12c(i, j, k) + Jup(i, j + 1, k) * g12up(i, j + 1, k)) - * (dfdx - coef_xy * dfdy)); + dzup = coord->dz.yup(); + dzdown = coord->dz.ydown(); - yzresult(i, j, k) = fout / (dyc(i, j, k) * Jc(i, j, k)); + Bxyup = coord->Bxy.yup(); + Bxydown = coord->Bxy.ydown(); - // Calculate flux between j and j-1 - coef_yz = - 0.5 - * (g_23c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) - + g_23down(i, j - 1, k) / SQ(Jdown(i, j - 1, k) * Bxydown(i, j - 1, k))); + } else { + // No 3D metrics + // Need to shift to/from field aligned coordinates + g23up = g23down = g23c = toFieldAligned(coord->g23); + g_23up = g_23down = g_23c = toFieldAligned(coord->g_23); + g12up = g12down = g12c = toFieldAligned(coord->g12); + g_12up = g_12down = g_12c = toFieldAligned(coord->g_12); + Jup = Jdown = Jc = toFieldAligned(coord->J); + dxc = toFieldAligned(coord->dx); + dyup = dydown = dyc = toFieldAligned(coord->dy); + dzup = dzdown = dzc = toFieldAligned(coord->dz); + Bxyup = Bxydown = Bxyc = toFieldAligned(coord->Bxy); + } - coef_xy = - 0.5 - * (g_12c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) - + g_12down(i, j - 1, k) / SQ(Jdown(i, j - 1, k) * Bxydown(i, j - 1, k))); + // Y flux + // Includes fluxes due to Z derivatives (non-zero g23 metric) + // and due to X derivatives if grid is nonorthogonal (non-zero g12 metric) - dfdz = 0.5 - * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) - fdown(i, j - 1, km)) - / (dzc(i, j, k) + dzdown(i, j - 1, k)); + 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 between j and j+1 + int kp = (k + 1) % mesh->LocalNz; + int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; - dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k)) - / (dyc(i, j, k) + dydown(i, j - 1, k)); - - dfdx = 0.25 - * (fddx_xhigh(i - 1, j, k) + fddx_xhigh(i, j, k) - + fddx_xhigh_down(i - 1, j - 1, k) + fddx_xhigh_down(i, j - 1, k)); - - fout = - 0.25 * (ac(i, j, k) + adown(i, j - 1, k)) - * ( - // Component of flux from d/dz and d/dy - (Jc(i, j, k) * g23c(i, j, k) + Jdown(i, j - 1, k) * g23down(i, j - 1, k)) - * (dfdz - coef_yz * dfdy) - // Non-orthogonal mesh correction with g12 metric - + (Jc(i, j, k) * g12c(i, j, k) - + Jdown(i, j - 1, k) * g12down(i, j - 1, k)) - * (dfdx - coef_xy * dfdy)); + BoutReal coef_yz = + 0.5 + * (g_23c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) + + g_23up(i, j + 1, k) / SQ(Jup(i, j + 1, k) * Bxyup(i, j + 1, k))); - yzresult(i, j, k) -= fout / (dyc(i, j, k) * Jc(i, j, k)); + BoutReal coef_xy = + 0.5 + * (g_12c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) + + g_12up(i, j + 1, k) / SQ(Jup(i, j + 1, k) * Bxyup(i, j + 1, k))); - // Flow will be positive in the positive coordinate direction - flow_ylow(i, j, k) = -1.0 * fout * coord->dx(i, j) * coord->dz(i, j); + // Calculate Z derivative at y boundary + BoutReal dfdz = + 0.5 + * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km)) + / (dzc(i, j, k) + dzup(i, j + 1, k)); + + // Y derivative + BoutReal dfdy = 2. * (fup(i, j + 1, k) - fc(i, j, k)) + / (dyup(i, j + 1, k) + dyc(i, j, k)); + + // X derivative at Y boundary + BoutReal dfdx = + 0.25 + * (fddx_xhigh(i - 1, j, k) + fddx_xhigh(i, j, k) + + fddx_xhigh_up(i - 1, j + 1, k) + fddx_xhigh_up(i, j + 1, k)); + + BoutReal fout = + 0.25 * (ac(i, j, k) + aup(i, j + 1, k)) + * ( + // Component of flux from d/dz and d/dy + (Jc(i, j, k) * g23c(i, j, k) + Jup(i, j + 1, k) * g23up(i, j + 1, k)) + * (dfdz - coef_yz * dfdy) + // Non-orthogonal mesh correction with g12 metric + + (Jc(i, j, k) * g12c(i, j, k) + + Jup(i, j + 1, k) * g12up(i, j + 1, k)) + * (dfdx - coef_xy * dfdy)); + + yzresult(i, j, k) = fout / (dyc(i, j, k) * Jc(i, j, k)); + + // Calculate flux between j and j-1 + coef_yz = 0.5 + * (g_23c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) + + g_23down(i, j - 1, k) + / SQ(Jdown(i, j - 1, k) * Bxydown(i, j - 1, k))); + + coef_xy = 0.5 + * (g_12c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) + + g_12down(i, j - 1, k) + / SQ(Jdown(i, j - 1, k) * Bxydown(i, j - 1, k))); + + dfdz = 0.5 + * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) + - fdown(i, j - 1, km)) + / (dzc(i, j, k) + dzdown(i, j - 1, k)); + + dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k)) + / (dyc(i, j, k) + dydown(i, j - 1, k)); + + dfdx = 0.25 + * (fddx_xhigh(i - 1, j, k) + fddx_xhigh(i, j, k) + + fddx_xhigh_down(i - 1, j - 1, k) + fddx_xhigh_down(i, j - 1, k)); + + fout = 0.25 * (ac(i, j, k) + adown(i, j - 1, k)) + * ( + // Component of flux from d/dz and d/dy + (Jc(i, j, k) * g23c(i, j, k) + + Jdown(i, j - 1, k) * g23down(i, j - 1, k)) + * (dfdz - coef_yz * dfdy) + // Non-orthogonal mesh correction with g12 metric + + (Jc(i, j, k) * g12c(i, j, k) + + Jdown(i, j - 1, k) * g12down(i, j - 1, k)) + * (dfdx - coef_xy * dfdy)); + + yzresult(i, j, k) -= fout / (dyc(i, j, k) * Jc(i, j, k)); + + // Flow will be positive in the positive coordinate direction + flow_ylow(i, j, k) = -1.0 * fout * coord->dx(i, j, k) * coord->dz(i, j, k); + } + } } - } - } - // Calculate the Y derivative, including in X guard cells - // This is used to calculate X flux contribution from g12 * d/dy - Field3D fddy(mesh); - fddy.allocate(); - fddy.setDirectionY(yzresult.getDirectionY()); - for (int i = mesh->xstart - 1; i <= mesh->xend + 1; i++) { - for (int j = mesh->ystart; j <= mesh->yend; j++) { - for (int k = 0; k < mesh->LocalNz; k++) { - fddy(i, j, k) = - (fup(i, j + 1, k) - fdown(i, j - 1, k)) - / (0.5 * dyup(i, j + 1, k) + dyc(i, j, k) + 0.5 * dydown(i, j - 1, k)); + // Calculate the Y derivative, including in X guard cells + // This is used to calculate X flux contribution from g12 * d/dy + Field3D fddy(mesh); + fddy.allocate(); + fddy.setDirectionY(yzresult.getDirectionY()); + for (int i = mesh->xstart - 1; i <= mesh->xend + 1; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + fddy(i, j, k) = + (fup(i, j + 1, k) - fdown(i, j - 1, k)) + / (0.5 * dyup(i, j + 1, k) + dyc(i, j, k) + 0.5 * dydown(i, j - 1, k)); + } + } } - } - } - // Z flux + // Z flux - 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 between k and k+1 - int kp = (k + 1) % mesh->LocalNz; - - BoutReal ddx = - 0.5 - * ((fc(i + 1, j, k) - fc(i - 1, j, k)) - / (0.5 * dxc(i + 1, j, k) + dxc(i, j, k) + 0.5 * dxc(i - 1, j, k)) - + (fc(i + 1, j, kp) - fc(i - 1, j, kp)) - / (0.5 * dxc(i + 1, j, kp) + dxc(i, j, kp) - + 0.5 * dxc(i - 1, j, kp))); + 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 between k and k+1 + int kp = (k + 1) % mesh->LocalNz; - BoutReal ddy = - 0.5 - * ((fup(i, j + 1, k) - fdown(i, j - 1, k)) - / (0.5 * dyup(i, j + 1, k) + dyc(i, j, k) + 0.5 * dydown(i, j - 1, k)) - + (fup(i, j + 1, kp) - fdown(i, j - 1, kp)) - / (0.5 * dyup(i, j + 1, kp) + dyc(i, j, kp) - + 0.5 * dydown(i, j - 1, kp))); - - BoutReal ddz = 2. * (fc(i, j, kp) - fc(i, j, k)) / (dzc(i, j, k) + dzc(i, j, kp)); - - BoutReal fout = - 0.5 * (ac(i, j, k) + ac(i, j, kp)) - * ( - // Jg^xz (d/dx - g_xy / g_yy d/dy) + BoutReal ddx = 0.5 - * (Jc(i, j, k) * coord->g13(i, j, k) - + Jc(i, j, kp) * coord->g13(i, j, kp)) - * (ddx - - ddy * 0.5 - * (g_12c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) - + g_12c(i, j, kp) / SQ(Jc(i, j, kp) * Bxyc(i, j, kp)))) - // Jg^zz (d/dz - g_yz / g_yy d/dy) - + 0.5 - * (Jc(i, j, k) * coord->g33(i, j, k) - + Jc(i, j, kp) * coord->g33(i, j, kp)) - * (ddz - - ddy * 0.5 - * (g_23c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) - + g_23c(i, j, kp) - / SQ(Jc(i, j, kp) * Bxyc(i, j, kp))))); - - yzresult(i, j, k) += fout / (Jc(i, j, k) * dzc(i, j, k)); - yzresult(i, j, kp) -= fout / (Jc(i, j, kp) * dzc(i, j, kp)); + * ((fc(i + 1, j, k) - fc(i - 1, j, k)) + / (0.5 * dxc(i + 1, j, k) + dxc(i, j, k) + 0.5 * dxc(i - 1, j, k)) + + (fc(i + 1, j, kp) - fc(i - 1, j, kp)) + / (0.5 * dxc(i + 1, j, kp) + dxc(i, j, kp) + + 0.5 * dxc(i - 1, j, kp))); + + BoutReal ddy = 0.5 + * ((fup(i, j + 1, k) - fdown(i, j - 1, k)) + / (0.5 * dyup(i, j + 1, k) + dyc(i, j, k) + + 0.5 * dydown(i, j - 1, k)) + + (fup(i, j + 1, kp) - fdown(i, j - 1, kp)) + / (0.5 * dyup(i, j + 1, kp) + dyc(i, j, kp) + + 0.5 * dydown(i, j - 1, kp))); + + BoutReal ddz = + 2. * (fc(i, j, kp) - fc(i, j, k)) / (dzc(i, j, k) + dzc(i, j, kp)); + + BoutReal fout = + 0.5 * (ac(i, j, k) + ac(i, j, kp)) + * ( + // Jg^xz (d/dx - g_xy / g_yy d/dy) + 0.5 + * (Jc(i, j, k) * coord->g13(i, j, k) + + Jc(i, j, kp) * coord->g13(i, j, kp)) + * (ddx + - ddy * 0.5 + * (g_12c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) + + g_12c(i, j, kp) + / SQ(Jc(i, j, kp) * Bxyc(i, j, kp)))) + // Jg^zz (d/dz - g_yz / g_yy d/dy) + + 0.5 + * (Jc(i, j, k) * coord->g33(i, j, k) + + Jc(i, j, kp) * coord->g33(i, j, kp)) + * (ddz + - ddy * 0.5 + * (g_23c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) + + g_23c(i, j, kp) + / SQ(Jc(i, j, kp) * Bxyc(i, j, kp))))); + + yzresult(i, j, k) += fout / (Jc(i, j, k) * dzc(i, j, k)); + yzresult(i, j, kp) -= fout / (Jc(i, j, kp) * dzc(i, j, kp)); + } + } } - } - } - // Check if we need to transform back - Field3D result; - if (f.hasParallelSlices() && a.hasParallelSlices()) { - result = yzresult; - } else { - result = fromFieldAligned(yzresult); - fddy = fromFieldAligned(fddy); - } - - // Flux in X + // Check if we need to transform back + Field3D result; + if (f.hasParallelSlices() && a.hasParallelSlices()) { + result = yzresult; + } else { + result = fromFieldAligned(yzresult); + fddy = fromFieldAligned(fddy); + } - for (int i = mesh->xstart - 1; 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 + // Flux in X - const int kp = (k + 1) % mesh->LocalNz; - const int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; + for (int i = mesh->xstart - 1; 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 ddx = 2 * (f(i + 1, j, k) - f(i, j, k)) - / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); - BoutReal ddy = 0.5 * (fddy(i, j, k) + fddy(i + 1, j, k)); + const int kp = (k + 1) % mesh->LocalNz; + const int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; - BoutReal ddz = 0.5 - * ((f(i, j, kp) - f(i, j, km)) - / (0.5 * coord->dz(i, j, kp) + coord->dz(i, j, k) - + 0.5 * coord->dz(i, j, km)) - + (f(i + 1, j, kp) - f(i + 1, j, km)) - / (0.5 * coord->dz(i + 1, j, kp) + coord->dz(i + 1, j, k) - + 0.5 * coord->dz(i + 1, j, km))); + BoutReal ddx = 2 * (f(i + 1, j, k) - f(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); + BoutReal ddy = 0.5 * (fddy(i, j, k) + fddy(i + 1, j, k)); - BoutReal fout = - 0.5 * (a(i, j, k) + a(i + 1, j, k)) - * ( - // Jg^xx (d/dx - g_xy / g_yy d/dy) + BoutReal ddz = 0.5 - * (coord->J(i, j, k) * coord->g11(i, j, k) - + coord->J(i + 1, j, k) * coord->g11(i + 1, j, k)) - * (ddx - - ddy * 0.5 - * (coord->g_12(i, j, k) / coord->g_22(i, j, k) - + coord->g_12(i + 1, j, k) / coord->g_22(i + 1, j, k))) - // Jg^xz (d/dz - g_yz / g_yy d/dy) - + 0.5 - * ((coord->J(i, j, k) * coord->g13(i, j, k) - + coord->J(i + 1, j, k) * coord->g13(i + 1, j, k)) - * (ddz - - ddy * 0.5 - * (coord->g_23(i, j, k) / coord->g_22(i, j, k) - + coord->g_23(i + 1, j, k) - / coord->g_22(i + 1, j, k))))); - - result(i, j, k) += fout / (coord->dx(i, j, k) * coord->J(i, j, k)); - result(i + 1, j, k) -= fout / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); - - // Flow will be positive in the positive coordinate direction - flow_xlow(i + 1, j, k) = -1.0 * fout * coord->dy(i, j) * coord->dz(i, j); + * ((f(i, j, kp) - f(i, j, km)) + / (0.5 * coord->dz(i, j, kp) + coord->dz(i, j, k) + + 0.5 * coord->dz(i, j, km)) + + (f(i + 1, j, kp) - f(i + 1, j, km)) + / (0.5 * coord->dz(i + 1, j, kp) + coord->dz(i + 1, j, k) + + 0.5 * coord->dz(i + 1, j, km))); + + BoutReal fout = + 0.5 * (a(i, j, k) + a(i + 1, j, k)) + * ( + // Jg^xx (d/dx - g_xy / g_yy d/dy) + 0.5 + * (coord->J(i, j, k) * coord->g11(i, j, k) + + coord->J(i + 1, j, k) * coord->g11(i + 1, j, k)) + * (ddx + - ddy * 0.5 + * (coord->g_12(i, j, k) / coord->g_22(i, j, k) + + coord->g_12(i + 1, j, k) + / coord->g_22(i + 1, j, k))) + // Jg^xz (d/dz - g_yz / g_yy d/dy) + + 0.5 + * ((coord->J(i, j, k) * coord->g13(i, j, k) + + coord->J(i + 1, j, k) * coord->g13(i + 1, j, k)) + * (ddz + - ddy * 0.5 + * (coord->g_23(i, j, k) / coord->g_22(i, j, k) + + coord->g_23(i + 1, j, k) + / coord->g_22(i + 1, j, k))))); + + result(i, j, k) += fout / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= + fout / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); + + // Flow will be positive in the positive coordinate direction + flow_xlow(i + 1, j, k) = + -1.0 * fout * coord->dy(i, j, k) * coord->dz(i, j, k); + } + } } - } - } - - return result; -} - -/// Div ( a Grad_perp(f) ) -- diffusion -/// -/// Returns the flows in the final arguments -/// -/// Flows are always in the positive {x,y} direction -/// i.e xlow(i,j) is the flow into cell (i,j) from the left, -/// and the flow out of cell (i-1,j) to the right -/// -/// ylow(i,j+1) -/// ^ -/// +---|---+ -/// | | -/// xlow(i,j) -> (i,j) -> xlow(i+1,j) -/// | ^ | -/// +---|---+ -/// ylow(i,j) -/// -/// -/// WARNING: Causes checkerboarding in neutral_mixed integrated test -const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, - Field3D& flow_xlow, Field3D& flow_ylow) { - ASSERT2(a.getLocation() == f.getLocation()); - - Mesh* mesh = a.getMesh(); - - Field3D result{zeroFrom(f)}; - - Coordinates* coord = f.getCoordinates(); - - // Zero all flows - flow_xlow = 0.0; - flow_ylow = 0.0; - - // Flux in x - - int xs = mesh->xstart - 1; - int xe = mesh->xend; - for (int i = xs; i <= xe; 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 - - const BoutReal gradient = (coord->J(i, j) * coord->g11(i, j) - + coord->J(i + 1, j) * coord->g11(i + 1, j)) - * (f(i + 1, j, k) - f(i, j, k)) - / (coord->dx(i, j) + coord->dx(i + 1, j)); - - // Use the upwind coefficient - const BoutReal fout = gradient * ((gradient > 0) ? a(i + 1, j, k) : a(i, j, k)); - - result(i, j, k) += fout / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= fout / (coord->dx(i + 1, j) * coord->J(i + 1, j)); - - // Flow will be positive in the positive coordinate direction - flow_xlow(i + 1, j, k) = -1.0 * fout * coord->dy(i, j) * coord->dz(i, j); - } + return result; } - } - - // Y and Z fluxes require Y derivatives - - // Fields containing values along the magnetic field - Field3D fup(mesh), fdown(mesh); - Field3D aup(mesh), adown(mesh); - - // Values on this y slice (centre). - // This is needed because toFieldAligned may modify the field - Field3D fc = f; - Field3D ac = a; - - // Result of the Y and Z fluxes - Field3D yzresult(mesh); - yzresult.allocate(); - if (f.hasParallelSlices() && a.hasParallelSlices()) { - // Both inputs have yup and ydown - - fup = f.yup(); - fdown = f.ydown(); - - aup = a.yup(); - adown = a.ydown(); - } else { - // At least one input doesn't have yup/ydown fields. - // Need to shift to/from field aligned coordinates - - fup = fdown = fc = toFieldAligned(f); - aup = adown = ac = toFieldAligned(a); - yzresult.setDirectionY(YDirectionType::Aligned); - flow_ylow.setDirectionY(YDirectionType::Aligned); - } - - // Y flux - - for (int i = mesh->xstart; i <= mesh->xend; i++) { - for (int j = mesh->ystart; j <= mesh->yend; j++) { - - BoutReal coef_u = - 0.5 - * (coord->g_23(i, j) / SQ(coord->J(i, j) * coord->Bxy(i, j)) - + coord->g_23(i, j + 1) / SQ(coord->J(i, j + 1) * coord->Bxy(i, j + 1))); - - BoutReal coef_d = - 0.5 - * (coord->g_23(i, j) / SQ(coord->J(i, j) * coord->Bxy(i, j)) - + coord->g_23(i, j - 1) / SQ(coord->J(i, j - 1) * coord->Bxy(i, j - 1))); - - for (int k = 0; k < mesh->LocalNz; k++) { - // Calculate flux between j and j+1 - int kp = (k + 1) % mesh->LocalNz; - int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; + /// Div ( a Grad_perp(f) ) -- diffusion + /// + /// Returns the flows in the final arguments + /// + /// Flows are always in the positive {x,y} direction + /// i.e xlow(i,j, k) is the flow into cell (i,j, k) from the left, + /// and the flow out of cell (i-1,j, k) to the right + /// + /// ylow(i,j+1) + /// ^ + /// +---|---+ + /// | | + /// xlow(i,j, k) -> (i,j, k) -> xlow(i+1,j, k) + /// | ^ | + /// +---|---+ + /// ylow(i,j, k) + /// + /// + /// WARNING: Causes checkerboarding in neutral_mixed integrated test + const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, + Field3D& flow_xlow, Field3D& flow_ylow) { + ASSERT2(a.getLocation() == f.getLocation()); + + Mesh* mesh = a.getMesh(); + + Field3D result{zeroFrom(f)}; + + Coordinates* coord = f.getCoordinates(); + + // Zero all flows + flow_xlow = 0.0; + flow_ylow = 0.0; + + // Flux in x + + int xs = mesh->xstart - 1; + int xe = mesh->xend; + + for (int i = xs; i <= xe; 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 + + const BoutReal gradient = (coord->J(i, j, k) * coord->g11(i, j, k) + + coord->J(i + 1, j, k) * coord->g11(i + 1, j, k)) + * (f(i + 1, j, k) - f(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); + + // Use the upwind coefficient + const BoutReal fout = + gradient * ((gradient > 0) ? a(i + 1, j, k) : a(i, j, k)); + + result(i, j, k) += fout / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= + fout / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); + + // Flow will be positive in the positive coordinate direction + flow_xlow(i + 1, j, k) = + -1.0 * fout * coord->dy(i, j, k) * coord->dz(i, j, k); + } + } + } - // Calculate Z derivative at y boundary - BoutReal dfdz = - 0.25 * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km)) - / coord->dz(i, j); + // Y and Z fluxes require Y derivatives - // Y derivative - BoutReal dfdy = 2. * (fup(i, j + 1, k) - fc(i, j, k)) - / (coord->dy(i, j + 1) + coord->dy(i, j)); + // Fields containing values along the magnetic field + Field3D fup(mesh), fdown(mesh); + Field3D aup(mesh), adown(mesh); - BoutReal fout = 0.25 * (ac(i, j, k) + aup(i, j + 1, k)) - * (coord->J(i, j) * coord->g23(i, j) - + coord->J(i, j + 1) * coord->g23(i, j + 1)) - * (dfdz - coef_u * dfdy); + // Values on this y slice (centre). + // This is needed because toFieldAligned may modify the field + Field3D fc = f; + Field3D ac = a; - yzresult(i, j, k) = fout / (coord->dy(i, j) * coord->J(i, j)); + // Result of the Y and Z fluxes + Field3D yzresult(mesh); + yzresult.allocate(); - // Calculate flux between j and j-1 - dfdz = 0.25 - * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) - fdown(i, j - 1, km)) - / coord->dz(i, j); + if (f.hasParallelSlices() && a.hasParallelSlices()) { + // Both inputs have yup and ydown - dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k)) - / (coord->dy(i, j) + coord->dy(i, j - 1)); - - fout = 0.25 * (ac(i, j, k) + adown(i, j - 1, k)) - * (coord->J(i, j) * coord->g23(i, j) - + coord->J(i, j - 1) * coord->g23(i, j - 1)) - * (dfdz - coef_d * dfdy); + fup = f.yup(); + fdown = f.ydown(); - yzresult(i, j, k) -= fout / (coord->dy(i, j) * coord->J(i, j)); + aup = a.yup(); + adown = a.ydown(); + } else { + // At least one input doesn't have yup/ydown fields. + // Need to shift to/from field aligned coordinates - // Flow will be positive in the positive coordinate direction - flow_ylow(i, j, k) = -1.0 * fout * coord->dx(i, j) * coord->dz(i, j); + fup = fdown = fc = toFieldAligned(f); + aup = adown = ac = toFieldAligned(a); + yzresult.setDirectionY(YDirectionType::Aligned); + flow_ylow.setDirectionY(YDirectionType::Aligned); } - } - } - // Z flux - // Easier since all metrics constant in Z + // Y flux - for (int i = mesh->xstart; i <= mesh->xend; i++) { - for (int j = mesh->ystart; j <= mesh->yend; j++) { - // Coefficient in front of df/dy term - BoutReal coef = coord->g_23(i, j) - / (coord->dy(i, j + 1) + 2. * coord->dy(i, j) + coord->dy(i, j - 1)) - / SQ(coord->J(i, j) * coord->Bxy(i, j)); - for (int k = 0; k < mesh->LocalNz; k++) { - // Calculate flux between k and k+1 - int kp = (k + 1) % mesh->LocalNz; - - BoutReal gradient = - // df/dz - (fc(i, j, kp) - fc(i, j, k)) / coord->dz(i, j) - - // - g_yz * df/dy / SQ(J*B) - - coef - * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k) - - fdown(i, j - 1, kp)); + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { +#if BOUT_USE_METRIC_3D + for (int k = 0; k < mesh->LocalNz; k++) { +#else + int k = 0; +#endif + BoutReal coef_u = + 0.5 + * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)) + + coord->g_23(i, j + 1, k) + / SQ(coord->J(i, j + 1, k) * coord->Bxy(i, j + 1, k))); - BoutReal fout = gradient * ((gradient > 0) ? ac(i, j, kp) : ac(i, j, k)); + BoutReal coef_d = + 0.5 + * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)) + + coord->g_23(i, j - 1, k) + / SQ(coord->J(i, j - 1, k) * coord->Bxy(i, j - 1, k))); + +#if not BOUT_USE_METRIC_3D + for (int k = 0; k < mesh->LocalNz; k++) { +#endif + // Calculate flux between j and j+1 + int kp = (k + 1) % mesh->LocalNz; + int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; + + // Calculate Z derivative at y boundary + BoutReal dfdz = + 0.25 + * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km)) + / coord->dz(i, j, k); + + // Y derivative + BoutReal dfdy = 2. * (fup(i, j + 1, k) - fc(i, j, k)) + / (coord->dy(i, j + 1, k) + coord->dy(i, j, k)); + + BoutReal fout = 0.25 * (ac(i, j, k) + aup(i, j + 1, k)) + * (coord->J(i, j, k) * coord->g23(i, j, k) + + coord->J(i, j + 1, k) * coord->g23(i, j + 1, k)) + * (dfdz - coef_u * dfdy); + + yzresult(i, j, k) = fout / (coord->dy(i, j, k) * coord->J(i, j, k)); + + // Calculate flux between j and j-1 + dfdz = 0.25 + * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) + - fdown(i, j - 1, km)) + / coord->dz(i, j, k); + + dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k)) + / (coord->dy(i, j, k) + coord->dy(i, j - 1, k)); + + fout = 0.25 * (ac(i, j, k) + adown(i, j - 1, k)) + * (coord->J(i, j, k) * coord->g23(i, j, k) + + coord->J(i, j - 1, k) * coord->g23(i, j - 1, k)) + * (dfdz - coef_d * dfdy); + + yzresult(i, j, k) -= fout / (coord->dy(i, j, k) * coord->J(i, j, k)); + + // Flow will be positive in the positive coordinate direction + flow_ylow(i, j, k) = -1.0 * fout * coord->dx(i, j, k) * coord->dz(i, j, k); + } + } + } - yzresult(i, j, k) += fout / coord->dz(i, j); - yzresult(i, j, kp) -= fout / coord->dz(i, j); - } - } - } - // Check if we need to transform back - if (f.hasParallelSlices() && a.hasParallelSlices()) { - result += yzresult; - } else { - result += fromFieldAligned(yzresult); - flow_ylow = fromFieldAligned(flow_ylow); - } + // Z flux + // Easier since all metrics constant in Z + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { +#if BOUT_USE_METRIC_3D + for (int k = 0; k < mesh->LocalNz; k++) { +#else + int k = 0; +#endif + // Coefficient in front of df/dy term + BoutReal coef = coord->g_23(i, j, k) + / (coord->dy(i, j + 1, k) + 2. * coord->dy(i, j, k) + + coord->dy(i, j - 1, k)) + / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)); +#if not BOUT_USE_METRIC_3D + for (int k = 0; k < mesh->LocalNz; k++) { +#endif + // Calculate flux between k and k+1 + int kp = (k + 1) % mesh->LocalNz; + + BoutReal gradient = + // df/dz + (fc(i, j, kp) - fc(i, j, k)) / coord->dz(i, j, k) + + // - g_yz * df/dy / SQ(J*B) + - coef + * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k) + - fdown(i, j - 1, kp)); + + BoutReal fout = gradient * ((gradient > 0) ? ac(i, j, kp) : ac(i, j, k)); + + yzresult(i, j, k) += fout / coord->dz(i, j, k); + yzresult(i, j, kp) -= fout / coord->dz(i, j, k); + } + } + } + // Check if we need to transform back + if (f.hasParallelSlices() && a.hasParallelSlices()) { + result += yzresult; + } else { + result += fromFieldAligned(yzresult); + flow_ylow = fromFieldAligned(flow_ylow); + } - return result; -} + return result; + } -const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, const Field3D& f, - bool bndry_flux) { - Field3D result{0.0}; + const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, + const Field3D& f, bool bndry_flux) { + Field3D result{0.0}; + + Coordinates* coord = mesh->getCoordinates(); + + ////////////////////////////////////////// + // X-Z advection. + // + // Z + // | + // + // fmp --- vU --- fpp + // | nU | + // | | + // vL nL nR vR -> X + // | | + // | nD | + // fmm --- vD --- fpm + // + + int nz = mesh->LocalNz; + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < nz; k++) { + int kp = (k + 1) % nz; + int kpp = (kp + 1) % nz; + int km = (k - 1 + nz) % nz; + int kmm = (km - 1 + nz) % nz; + + // 1) Interpolate stream function f onto corners fmp, fpp, fpm + + BoutReal fmm = + 0.25 * (f(i, j, k) + f(i - 1, j, k) + f(i, j, km) + f(i - 1, j, km)); + BoutReal fmp = 0.25 + * (f(i, j, k) + f(i, j, kp) + f(i - 1, j, k) + + f(i - 1, j, kp)); // 2nd order accurate + BoutReal fpp = + 0.25 * (f(i, j, k) + f(i, j, kp) + f(i + 1, j, k) + f(i + 1, j, kp)); + BoutReal fpm = + 0.25 * (f(i, j, k) + f(i + 1, j, k) + f(i, j, km) + f(i + 1, j, km)); + + // 2) Calculate velocities on cell faces + + BoutReal vU = + coord->J(i, j, k) * (fmp - fpp) / coord->dx(i, j, k); // -J*df/dx + BoutReal vD = + coord->J(i, j, k) * (fmm - fpm) / coord->dx(i, j, k); // -J*df/dx + + BoutReal vR = 0.5 * (coord->J(i, j, k) + coord->J(i + 1, j, k)) + * (fpp - fpm) / coord->dz(i, j, k); // J*df/dz + BoutReal vL = 0.5 * (coord->J(i, j, k) + coord->J(i - 1, j, k)) + * (fmp - fmm) / coord->dz(i, j, k); // J*df/dz + + // 3) Calculate g on cell faces + + BoutReal gU = 0.5 * (g(i, j, kp) + g(i, j, k)); + BoutReal gD = 0.5 * (g(i, j, km) + g(i, j, k)); + BoutReal gR = 0.5 * (g(i + 1, j, k) + g(i, j, k)); + BoutReal gL = 0.5 * (g(i - 1, j, k) + g(i, j, k)); + + // 4) Calculate n on the cell faces. The sign of the + // velocity determines which side is used. + + // X direction + Stencil1D s; + s.c = n(i, j, k); + s.m = n(i - 1, j, k); + s.mm = n(i - 2, j, k); + s.p = n(i + 1, j, k); + s.pp = n(i + 2, j, k); + + MC(s); + + // Right side + if ((i == mesh->xend) && (mesh->lastX())) { + // At right boundary in X + + if (bndry_flux) { + BoutReal flux; + if (vR > 0.0) { + // Flux to boundary + flux = vR * s.R * gR; + } else { + // Flux in from boundary + flux = vR * 0.5 * (n(i + 1, j, k) + n(i, j, k)) * gR; + } + + result(i, j, k) += flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= + flux / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); + } + } else { + // Not at a boundary + if (vR > 0.0) { + // Flux out into next cell + BoutReal flux = vR * s.R * gR; + result(i, j, k) += flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= + flux / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); + } + } + + // Left side + + if ((i == mesh->xstart) && (mesh->firstX())) { + // At left boundary in X + + if (bndry_flux) { + BoutReal flux; + + if (vL < 0.0) { + // Flux to boundary + flux = vL * s.L * gL; + + } else { + // Flux in from boundary + flux = vL * 0.5 * (n(i - 1, j, k) + n(i, j, k)) * gL; + } + result(i, j, k) -= flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i - 1, j, k) += + flux / (coord->dx(i - 1, j, k) * coord->J(i - 1, j, k)); + } + } else { + // Not at a boundary + + if (vL < 0.0) { + BoutReal flux = vL * s.L * gL; + result(i, j, k) -= flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i - 1, j, k) += + flux / (coord->dx(i - 1, j, k) * coord->J(i - 1, j, k)); + } + } + + /// NOTE: Need to communicate fluxes + + // Z direction + s.m = n(i, j, km); + s.mm = n(i, j, kmm); + s.p = n(i, j, kp); + s.pp = n(i, j, kpp); + + MC(s); + + if (vU > 0.0) { + BoutReal flux = + vU * s.R * gU / (coord->J(i, j, k) * coord->dz(i, j, k)); + result(i, j, k) += flux; + result(i, j, kp) -= flux; + } + if (vD < 0.0) { + BoutReal flux = + vD * s.L * gD / (coord->J(i, j, k) * coord->dz(i, j, k)); + result(i, j, k) -= flux; + result(i, j, km) += flux; + } + } + } + } + FV::communicateFluxes(result); + return result; + } - Coordinates* coord = mesh->getCoordinates(); + const Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, + Field3D& flow_ylow, bool bndry_flux) { + TRACE("FV::Div_par_K_Grad_par_mod"); - ////////////////////////////////////////// - // X-Z advection. - // - // Z - // | - // - // fmp --- vU --- fpp - // | nU | - // | | - // vL nL nR vR -> X - // | | - // | nD | - // fmm --- vD --- fpm - // + ASSERT2(Kin.getLocation() == fin.getLocation()); - int nz = mesh->LocalNz; - for (int i = mesh->xstart; i <= mesh->xend; i++) { - for (int j = mesh->ystart; j <= mesh->yend; j++) { - for (int k = 0; k < nz; k++) { - int kp = (k + 1) % nz; - int kpp = (kp + 1) % nz; - int km = (k - 1 + nz) % nz; - int kmm = (km - 1 + nz) % nz; + Mesh* mesh = Kin.getMesh(); - // 1) Interpolate stream function f onto corners fmp, fpp, fpm + bool use_parallel_slices = (Kin.hasParallelSlices() && fin.hasParallelSlices()); - BoutReal fmm = - 0.25 * (f(i, j, k) + f(i - 1, j, k) + f(i, j, km) + f(i - 1, j, km)); - BoutReal fmp = 0.25 - * (f(i, j, k) + f(i, j, kp) + f(i - 1, j, k) - + f(i - 1, j, kp)); // 2nd order accurate - BoutReal fpp = - 0.25 * (f(i, j, k) + f(i, j, kp) + f(i + 1, j, k) + f(i + 1, j, kp)); - BoutReal fpm = - 0.25 * (f(i, j, k) + f(i + 1, j, k) + f(i, j, km) + f(i + 1, j, km)); + const auto& K = use_parallel_slices ? Kin : toFieldAligned(Kin, "RGN_NOX"); + const auto& f = use_parallel_slices ? fin : toFieldAligned(fin, "RGN_NOX"); - // 2) Calculate velocities on cell faces + Field3D result{zeroFrom(f)}; + flow_ylow = zeroFrom(f); - BoutReal vU = coord->J(i, j) * (fmp - fpp) / coord->dx(i, j); // -J*df/dx - BoutReal vD = coord->J(i, j) * (fmm - fpm) / coord->dx(i, j); // -J*df/dx + // K and f fields in yup and ydown directions + const auto& Kup = use_parallel_slices ? Kin.yup() : K; + const auto& Kdown = use_parallel_slices ? Kin.ydown() : K; + const auto& fup = use_parallel_slices ? fin.yup() : f; + const auto& fdown = use_parallel_slices ? fin.ydown() : f; - BoutReal vR = 0.5 * (coord->J(i, j) + coord->J(i + 1, j)) * (fpp - fpm) - / coord->dz(i, j); // J*df/dz - BoutReal vL = 0.5 * (coord->J(i, j) + coord->J(i - 1, j)) * (fmp - fmm) - / coord->dz(i, j); // J*df/dz + Coordinates* coord = fin.getCoordinates(); - // 3) Calculate g on cell faces + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + // Calculate flux at upper surface - BoutReal gU = 0.5 * (g(i, j, kp) + g(i, j, k)); - BoutReal gD = 0.5 * (g(i, j, km) + g(i, j, k)); - BoutReal gR = 0.5 * (g(i + 1, j, k) + g(i, j, k)); - BoutReal gL = 0.5 * (g(i - 1, j, k) + g(i, j, k)); + const auto iyp = i.yp(); + const auto iym = i.ym(); - // 4) Calculate n on the cell faces. The sign of the - // velocity determines which side is used. + if (bndry_flux || mesh->periodicY(i.x()) || !mesh->lastY(i.x()) + || (i.y() != mesh->yend)) { - // X direction - Stencil1D s; - s.c = n(i, j, k); - s.m = n(i - 1, j, k); - s.mm = n(i - 2, j, k); - s.p = n(i + 1, j, k); - s.pp = n(i + 2, j, k); + BoutReal c = 0.5 * (K[i] + Kup[iyp]); // K at the upper boundary + BoutReal J = 0.5 * (coord->J[i] + coord->J[iyp]); // Jacobian at boundary + BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); - MC(s); + BoutReal gradient = + 2. * (fup[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); - // Right side - if ((i == mesh->xend) && (mesh->lastX())) { - // At right boundary in X + BoutReal flux = c * J * gradient / g_22; - if (bndry_flux) { - BoutReal flux; - if (vR > 0.0) { - // Flux to boundary - flux = vR * s.R * gR; - } else { - // Flux in from boundary - flux = vR * 0.5 * (n(i + 1, j, k) + n(i, j, k)) * gR; + result[i] += flux / (coord->dy[i] * coord->J[i]); } - result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); - } - } else { - // Not at a boundary - if (vR > 0.0) { - // Flux out into next cell - BoutReal flux = vR * s.R * gR; - result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); - } - } - - // Left side + // Calculate flux at lower surface + if (bndry_flux || mesh->periodicY(i.x()) || !mesh->firstY(i.x()) + || (i.y() != mesh->ystart)) { + BoutReal c = 0.5 * (K[i] + Kdown[iym]); // K at the lower boundary + BoutReal J = 0.5 * (coord->J[i] + coord->J[iym]); // Jacobian at boundary - if ((i == mesh->xstart) && (mesh->firstX())) { - // At left boundary in X + BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); - if (bndry_flux) { - BoutReal flux; + BoutReal gradient = + 2. * (f[i] - fdown[iym]) / (coord->dy[i] + coord->dy[iym]); - if (vL < 0.0) { - // Flux to boundary - flux = vL * s.L * gL; + BoutReal flux = c * J * gradient / g_22; - } else { - // Flux in from boundary - flux = vL * 0.5 * (n(i - 1, j, k) + n(i, j, k)) * gL; + result[i] -= flux / (coord->dy[i] * coord->J[i]); + flow_ylow[i] = -flux * coord->dx[i] * coord->dz[i]; } - result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); - result(i - 1, j, k) += flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); } - } else { - // Not at a boundary - if (vL < 0.0) { - BoutReal flux = vL * s.L * gL; - result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); - result(i - 1, j, k) += flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); + if (!use_parallel_slices) { + // Shifted to field aligned coordinates, so need to shift back + result = fromFieldAligned(result, "RGN_NOBNDRY"); + flow_ylow = fromFieldAligned(flow_ylow); } - } - - /// NOTE: Need to communicate fluxes - - // Z direction - s.m = n(i, j, km); - s.mm = n(i, j, kmm); - s.p = n(i, j, kp); - s.pp = n(i, j, kpp); - - MC(s); - if (vU > 0.0) { - BoutReal flux = vU * s.R * gU / (coord->J(i, j) * coord->dz(i, j)); - result(i, j, k) += flux; - result(i, j, kp) -= flux; - } - if (vD < 0.0) { - BoutReal flux = vD * s.L * gD / (coord->J(i, j) * coord->dz(i, j)); - result(i, j, k) -= flux; - result(i, j, km) += flux; + return result; } - } - } - } - FV::communicateFluxes(result); - return result; -} From 00907b4a4363f7f55ab108247da16424502d134e Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 23 Jun 2026 13:59:08 +0200 Subject: [PATCH 13/24] Disable full neutral model for 3D metrics --- include/neutral_full_velocity.hxx | 4 +++- src/neutral_full_velocity.cxx | 4 ++++ tests/unit/test_neutral_full_velocity.cxx | 7 ++++++- 3 files changed, 13 insertions(+), 2 deletions(-) diff --git a/include/neutral_full_velocity.hxx b/include/neutral_full_velocity.hxx index ccca38537..175a96f82 100644 --- a/include/neutral_full_velocity.hxx +++ b/include/neutral_full_velocity.hxx @@ -9,6 +9,8 @@ #include +#if not BOUT_USE_METRIC_3D + /// Neutral gas model, evolving three components of velocity as axisymmetric fields /// /// Evolves neutral density, pressure and velocity as Field2D quantities @@ -100,5 +102,5 @@ namespace { RegisterComponent registersolverneutralfullvelocity("neutral_full_velocity"); } - +#endif // fu #endif // NEUTRAL_FULL_VELOCITY_H diff --git a/src/neutral_full_velocity.cxx b/src/neutral_full_velocity.cxx index ff690f9a5..c5a630944 100644 --- a/src/neutral_full_velocity.cxx +++ b/src/neutral_full_velocity.cxx @@ -10,6 +10,8 @@ #include +#if not BOUT_USE_METRIC_3D + using bout::globals::mesh; NeutralFullVelocity::NeutralFullVelocity(const std::string& name, Options& alloptions, @@ -855,3 +857,5 @@ void NeutralFullVelocity::outputVars(Options& state) { {"source", "neutral_full_velocity"}}); } } + +#endif diff --git a/tests/unit/test_neutral_full_velocity.cxx b/tests/unit/test_neutral_full_velocity.cxx index 1c10e8211..043d1273d 100644 --- a/tests/unit/test_neutral_full_velocity.cxx +++ b/tests/unit/test_neutral_full_velocity.cxx @@ -8,6 +8,7 @@ #include // For generating functions +#if not BOUT_USE_METRIC_3D /// Global mesh namespace bout { namespace globals { @@ -89,7 +90,9 @@ TEST_F(NeutralFullVelocityTest, CreateComponentRequiresBpxy) { static_cast(bout::globals::mesh) ->setGridDataSource(new FakeGridDataSource{{ - {"Rxy", 0.0}, {"Zxy", 1.0}, {"hthe", 1.0}, + {"Rxy", 0.0}, + {"Zxy", 1.0}, + {"hthe", 1.0}, // Missing Bpxy }}); @@ -316,3 +319,5 @@ TEST_F(NeutralFullVelocityTest, OutputVars) { ASSERT_TRUE(outputs.isSet("Urx")); ASSERT_TRUE(outputs.isSet("Tyr")); } + +#endif From 8de567b0f048d84a8e1cc5cb59af9c0ede64236b Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 23 Jun 2026 13:59:56 +0200 Subject: [PATCH 14/24] Fix PolarisationDrift for 3D compilation --- include/polarisation_drift.hxx | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/include/polarisation_drift.hxx b/include/polarisation_drift.hxx index 8d0a5a3c6..293f2d3d6 100644 --- a/include/polarisation_drift.hxx +++ b/include/polarisation_drift.hxx @@ -2,6 +2,8 @@ #ifndef POLARISATION_DRIFT_H #define POLARISATION_DRIFT_H +#include + #include "component.hxx" class Laplacian; @@ -30,9 +32,9 @@ class Laplacian; /// charged species (ions and electrons) struct PolarisationDrift : public Component { // - PolarisationDrift(std::string name, Options &options, Solver *UNUSED(solver)); + PolarisationDrift(std::string name, Options& options, Solver* UNUSED(solver)); - void outputVars(Options &state) override; + void outputVars(Options& state) override; // The following functions are public for unit testing @@ -61,21 +63,22 @@ struct PolarisationDrift : public Component { /// Sets density_source, energy_source and momentum_source /// for all species with mass and charge. void polarisationAdvection(GuardedOptions& state, Field3D phi_pol); + private: std::unique_ptr phiSolver; // Laplacian solver in X-Z - Field2D Bsq; // Cached SQ(coord->Bxy) - + Coordinates::FieldMetric Bsq; // Cached SQ(coord->Bxy) + // Diagnostic outputs - bool diagnose; ///< Save diagnostic outputs? - Field3D DivJ; ///< Divergence of all other currents - Field3D phi_pol; ///< Polarisation drift potential + bool diagnose; ///< Save diagnostic outputs? + Field3D DivJ; ///< Divergence of all other currents + Field3D phi_pol; ///< Polarisation drift potential Options diagnostics; ///< Other diagnostic outputs - bool boussinesq; // If true, assume a constant mass density in Jpol - BoutReal average_atomic_mass; // If boussinesq=true, mass density to use - BoutReal density_floor; // Minimum mass density if boussinesq=false - bool advection; // Advect fluids by an approximate polarisation velocity? + bool boussinesq; // If true, assume a constant mass density in Jpol + BoutReal average_atomic_mass; // If boussinesq=true, mass density to use + BoutReal density_floor; // Minimum mass density if boussinesq=false + bool advection; // Advect fluids by an approximate polarisation velocity? bool diamagnetic_polarisation; // Calculate compression terms? /// Inputs @@ -104,7 +107,8 @@ private: }; namespace { -RegisterComponent registercomponentpolarisationdrift("polarisation_drift"); +RegisterComponent + registercomponentpolarisationdrift("polarisation_drift"); } #endif // POLARISATION_DRIFT_H From e5a530e0d6b9ba0f55e9be2eecaea0eccf6744cf Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 23 Jun 2026 14:00:41 +0200 Subject: [PATCH 15/24] Fix RelaxPotential for 3D compilation --- include/relax_potential.hxx | 31 ++++++++++++++++--------------- src/relax_potential.cxx | 2 +- 2 files changed, 17 insertions(+), 16 deletions(-) diff --git a/include/relax_potential.hxx b/include/relax_potential.hxx index 227a78942..36912cf45 100644 --- a/include/relax_potential.hxx +++ b/include/relax_potential.hxx @@ -2,7 +2,8 @@ #ifndef RELAX_POTENTIAL_H #define RELAX_POTENTIAL_H -#include +#include +#include #include "component.hxx" @@ -102,25 +103,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 + Coordinates::FieldMetric Bsq; ///< SQ(coord->Bxy) + VectorMetric Curlb_B; ///< Curvature vector Curl(b/B) + BoutReal hyper_z; ///< Hyper-viscosity in Z + Coordinates::FieldMetric viscosity; ///< Perpendicular Kinematic viscosity + Coordinates::FieldMetric 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..571f74e8a 100644 --- a/src/relax_potential.cxx +++ b/src/relax_potential.cxx @@ -198,7 +198,7 @@ RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* so // Read curvature vector try { // May be 2D, reading as 3D - Vector2D curv2d; + VectorMetric curv2d; curv2d.covariant = false; mesh->get(curv2d, "bxcv"); Curlb_B = curv2d; From 3e11b80d3b855253dfa4bf04f0a5f8b2d1654098 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 23 Jun 2026 14:01:13 +0200 Subject: [PATCH 16/24] Fix ScaleTimederivs for 3D compilation --- include/scale_timederivs.hxx | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/include/scale_timederivs.hxx b/include/scale_timederivs.hxx index 8476e3654..25cc02498 100644 --- a/include/scale_timederivs.hxx +++ b/include/scale_timederivs.hxx @@ -15,12 +15,12 @@ struct ScaleTimeDerivs : public Component { : Component({readOnly("species:e:temperature"), writeFinal("scale_timederivs")}) {} void outputVars(Options& state) override { - set_with_attrs( - state["scale_timederivs"], scaling, - {{"time_dimension", "t"}, - {"long_name", "Scaling factor applied to all time derivatives"}, - {"source", "scale_timederivs"}}); + set_with_attrs(state["scale_timederivs"], scaling, + {{"time_dimension", "t"}, + {"long_name", "Scaling factor applied to all time derivatives"}, + {"source", "scale_timederivs"}}); } + private: Field3D scaling; // The scaling factor applied to each cell @@ -37,7 +37,7 @@ private: void transform_impl(GuardedOptions& state) override { auto* coord = bout::globals::mesh->getCoordinates(); - Field2D dl2 = coord->g_22 * SQ(coord->dy); + auto dl2 = coord->g_22 * SQ(coord->dy); // Scale by parallel heat conduction CFL timescale auto Te = get(state["species"]["e"]["temperature"]); @@ -49,8 +49,7 @@ private: }; namespace { - RegisterComponent registercomponentscaletimederivs("scale_timederivs"); +RegisterComponent registercomponentscaletimederivs("scale_timederivs"); } #endif // SCALE_TIMEDERIVS_H - From 0318c5ca7256edac0559406996f7b248078d149e Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 23 Jun 2026 14:01:36 +0200 Subject: [PATCH 17/24] Fix Vorticity for 3D compilation --- include/vorticity.hxx | 65 ++++++++++++++++++++++--------------------- src/vorticity.cxx | 2 +- 2 files changed, 34 insertions(+), 33 deletions(-) diff --git a/include/vorticity.hxx b/include/vorticity.hxx index 62d5d7899..8021a6b70 100644 --- a/include/vorticity.hxx +++ b/include/vorticity.hxx @@ -2,16 +2,17 @@ #ifndef VORTICITY_H #define VORTICITY_H -#include +#include + +#include #include -#include +#include +#include #include "component.hxx" -#include - /// Evolve electron density in time -/// +/// struct Vorticity : public Component { /// Options /// @@ -48,23 +49,23 @@ struct Vorticity : public Component { /// If phi_boundary_relax is false, set the radial boundary to the sheath potential? /// - split_n0: bool, default false /// Split phi into n=0 and n!=0 components? - /// - viscosity: Field2D, default 0.0 + /// - viscosity: FieldMetric, default 0.0 /// Kinematic viscosity [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 /// - Vorticity(std::string name, Options &options, Solver *solver); + Vorticity(std::string name, Options& options, Solver* solver); /// Optional inputs /// - fields /// - DivJextra Divergence of current, including parallel current /// Not including diamagnetic or polarisation currents - /// - void finally(const Options &state) override; + /// + void finally(const Options& state) override; - void outputVars(Options &state) override; + void outputVars(Options& state) override; // Save and restore potential phi void restartVars(Options& state) override { @@ -80,8 +81,7 @@ struct Vorticity : public Component { // Save the potential set_with_attrs(state["phi"], phi, - {{"long_name", "plasma potential"}, - {"source", "vorticity"}}); + {{"long_name", "plasma potential"}, {"source", "vorticity"}}); } // The following are public functions for unit testing @@ -97,42 +97,43 @@ struct Vorticity : public Component { private: Field3D Vort; // Evolving vorticity - Field3D phi; // Electrostatic potential + Field3D phi; // Electrostatic potential std::unique_ptr phiSolver; // Laplacian solver in X-Z Field3D Pi_hat; ///< Contribution from ion pressure, weighted by atomic mass / charge - bool exb_advection; // Include nonlinear ExB advection? + bool exb_advection; // Include nonlinear ExB advection? bool exb_advection_simplified; // Simplify nonlinear ExB advection form? - bool diamagnetic; // Include diamagnetic current? + 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) - bool poloidal_flows; ///< Include poloidal ExB flow? - bool bndry_flux; ///< Allow flows through radial boundaries? + 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? bool collisional_friction; ///< Damping of vorticity due to collisional friction 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 split_n0; // Split phi into n=0 and n!=0 components + bool phi_core_averagey; ///< Average phi core boundary in Y? + + bool split_n0; // Split phi into n=0 and n!=0 components std::unique_ptr laplacexy; // Laplacian solver in X-Y (n=0) - Field2D Bsq; // SQ(coord->Bxy) - Vector2D Curlb_B; // Curvature vector Curl(b/B) - BoutReal hyper_z; ///< Hyper-viscosity in Z - Field2D viscosity; ///< Kinematic viscosity - Field3D viscous_heating; ///< Heating due to kinematic viscosity - bool include_viscosity; ///< Is viscosity > 0? + Coordinates::FieldMetric Bsq; // SQ(coord->Bxy) + VectorMetric Curlb_B; // Curvature vector Curl(b/B) + BoutReal hyper_z; ///< Hyper-viscosity in Z + Coordinates::FieldMetric viscosity; ///< Kinematic viscosity + Field3D viscous_heating; ///< Heating due to kinematic viscosity + bool include_viscosity; ///< Is viscosity > 0? // Diagnostic outputs Field3D DivJdia, DivJcol; // Divergence of diamagnetic and collisional current diff --git a/src/vorticity.cxx b/src/vorticity.cxx index d2e9bae70..0bf29e148 100644 --- a/src/vorticity.cxx +++ b/src/vorticity.cxx @@ -190,7 +190,7 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) // Create an XY solver for n=0 component laplacexy = LaplaceXY::create(mesh); // Set coefficients for Boussinesq solve - laplacexy->setCoefs(average_atomic_mass / SQ(coord->Bxy), 0.0); + laplacexy->setCoefs(average_atomic_mass / SQ(DC(coord->Bxy)), 0.0); } phiSolver = Laplacian::create(&options["laplacian"]); // Set coefficients for Boussinesq solve From 9e6a10cfe7e9ef18d6b0a5c8cc32766c5c12fe3a Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 23 Jun 2026 14:02:01 +0200 Subject: [PATCH 18/24] Fix Electromagnetic for 3D compilation --- src/electromagnetic.cxx | 96 ++++++++++++++++++++--------------------- 1 file changed, 48 insertions(+), 48 deletions(-) diff --git a/src/electromagnetic.cxx b/src/electromagnetic.cxx index c51a9d6ab..2854ee5b0 100644 --- a/src/electromagnetic.cxx +++ b/src/electromagnetic.cxx @@ -1,8 +1,8 @@ #include "../include/electromagnetic.hxx" #include -#include #include +#include // Set the default acceptance tolerances for the Naulin solver. // These are used if the maximum iterations is reached. @@ -43,8 +43,8 @@ Electromagnetic::Electromagnetic(std::string name, Options& alloptions, Solver* } const_gradient = options["const_gradient"] - .doc("Extrapolate gradient of Apar into all radial boundaries?") - .withDefault(false); + .doc("Extrapolate gradient of Apar into all radial boundaries?") + .withDefault(false); // Give Apar an initial value because we solve Apar by iteration // starting from the previous solution @@ -58,35 +58,35 @@ Electromagnetic::Electromagnetic(std::string name, Options& alloptions, Solver* last_time = 0.0; apar_boundary_timescale = options["apar_boundary_timescale"] - .doc("Timescale for Apar boundary relaxation [seconds]") - .withDefault(1e-8) - / get(alloptions["units"]["seconds"]); + .doc("Timescale for Apar boundary relaxation [seconds]") + .withDefault(1e-8) + / get(alloptions["units"]["seconds"]); } else if (options["apar_boundary_neumann"] - .doc("Neumann on all radial boundaries?") - .withDefault(false)) { + .doc("Neumann on all radial boundaries?") + .withDefault(false)) { // Set zero-gradient (neumann) boundary condition DC on the core aparSolver->setInnerBoundaryFlags(INVERT_DC_GRAD + INVERT_AC_GRAD); aparSolver->setOuterBoundaryFlags(INVERT_DC_GRAD + INVERT_AC_GRAD); } else if (options["apar_core_neumann"] - .doc("Neumann radial boundary in the core? False => Dirichlet") - .withDefault(true) + .doc("Neumann radial boundary in the core? False => Dirichlet") + .withDefault(true) and bout::globals::mesh->periodicY(bout::globals::mesh->xstart)) { // Set zero-gradient (neumann) boundary condition DC on the core aparSolver->setInnerBoundaryFlags(INVERT_DC_GRAD); } - diagnose = options["diagnose"] - .doc("Output additional diagnostics?") - .withDefault(false); + diagnose = + options["diagnose"].doc("Output additional diagnostics?").withDefault(false); magnetic_flutter = options["magnetic_flutter"] - .doc("Set magnetic flutter terms (Apar_flutter)?") - .withDefault(false); + .doc("Set magnetic flutter terms (Apar_flutter)?") + .withDefault(false); - if (magnetic_flutter) + if (magnetic_flutter) { setPermissions(readWrite("fields:Apar_flutter")); + } } void Electromagnetic::restartVars(Options& state) { @@ -98,7 +98,8 @@ void Electromagnetic::restartVars(Options& state) { if (first and state.isSet("Apar")) { first = false; Apar = state["Apar"].as(); - output.write("\nElectromagnetic: Read Apar from restart file (min {}, max {})\n", min(Apar), max(Apar)); + output.write("\nElectromagnetic: Read Apar from restart file (min {}, max {})\n", + min(Apar), max(Apar)); } // Save the Apar field. It is solved using an iterative method, @@ -163,9 +164,9 @@ void Electromagnetic::transform_impl(GuardedOptions& state) { const int x = mesh->xstart - 1; for (int y = mesh->ystart; y <= mesh->yend; y++) { for (int z = mesh->zstart; z <= mesh->zend; z++) { - rhs(x, y, z) = (weight * (Apar(x + 1, y, z) - Apar(x, y, z)) + - (1 - weight) * (Apar(x + 2, y, z) - Apar(x + 1, y, z))) / - (sqrt(coords->g_11(x, y)) * coords->dx(x, y)); + rhs(x, y, z) = (weight * (Apar(x + 1, y, z) - Apar(x, y, z)) + + (1 - weight) * (Apar(x + 2, y, z) - Apar(x + 1, y, z))) + / (sqrt(coords->g_11(x, y, z)) * coords->dx(x, y, z)); } } } @@ -173,9 +174,9 @@ void Electromagnetic::transform_impl(GuardedOptions& state) { const int x = mesh->xend + 1; for (int y = mesh->ystart; y <= mesh->yend; y++) { for (int z = mesh->zstart; z <= mesh->zend; z++) { - rhs(x, y, z) = (weight * (Apar(x, y, z) - Apar(x - 1, y, z)) + - (1 - weight) * (Apar(x - 1, y, z) - Apar(x - 2, y, z))) / - sqrt(coords->g_11(x, y)) / coords->dx(x, y); + rhs(x, y, z) = (weight * (Apar(x, y, z) - Apar(x - 1, y, z)) + + (1 - weight) * (Apar(x - 1, y, z) - Apar(x - 2, y, z))) + / sqrt(coords->g_11(x, y, z)) / coords->dx(x, y, z); } } } @@ -247,40 +248,39 @@ void Electromagnetic::transform_impl(GuardedOptions& state) { } } -void Electromagnetic::outputVars(Options &state) { +void Electromagnetic::outputVars(Options& state) { // Normalisations auto Bnorm = get(state["Bnorm"]); auto rho_s0 = get(state["rho_s0"]); - set_with_attrs(state["beta_em"], beta_em, { - {"long_name", "Helmholtz equation parameter"} - }); + set_with_attrs(state["beta_em"], beta_em, + {{"long_name", "Helmholtz equation parameter"}}); - set_with_attrs(state["Apar"], Apar, { - {"time_dimension", "t"}, - {"units", "T m"}, - {"conversion", Bnorm * rho_s0}, - {"standard_name", "b dot A"}, - {"long_name", "Parallel component of vector potential A"} - }); + set_with_attrs(state["Apar"], Apar, + {{"time_dimension", "t"}, + {"units", "T m"}, + {"conversion", Bnorm * rho_s0}, + {"standard_name", "b dot A"}, + {"long_name", "Parallel component of vector potential A"}}); if (magnetic_flutter) { - set_with_attrs(state["Apar_flutter"], Apar_flutter, { - {"time_dimension", "t"}, - {"units", "T m"}, - {"conversion", Bnorm * rho_s0}, - {"standard_name", "b dot A"}, - {"long_name", "Vector potential A|| used in flutter terms"} - }); + set_with_attrs(state["Apar_flutter"], Apar_flutter, + {{"time_dimension", "t"}, + {"units", "T m"}, + {"conversion", Bnorm * rho_s0}, + {"standard_name", "b dot A"}, + {"long_name", "Vector potential A|| used in flutter terms"}}); } if (diagnose) { - set_with_attrs(state["Ajpar"], Ajpar, { - {"time_dimension", "t"}, - }); - - set_with_attrs(state["alpha_em"], alpha_em, { - {"time_dimension", "t"}, - }); + set_with_attrs(state["Ajpar"], Ajpar, + { + {"time_dimension", "t"}, + }); + + set_with_attrs(state["alpha_em"], alpha_em, + { + {"time_dimension", "t"}, + }); } } From 0a94913a39373da189fec1d22957498ce25c5207 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 23 Jun 2026 14:02:20 +0200 Subject: [PATCH 19/24] Fix Recycling for 3D compilation --- src/recycling.cxx | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/src/recycling.cxx b/src/recycling.cxx index 8cf6cf5f6..f66ff808e 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -54,8 +54,9 @@ Recycling::Recycling(std::string name, Options& alloptions, Solver*) for (const auto& species : species_list) { std::string from = trim(species, " \t\r()"); // The species name in the list - if (from.empty()) + if (from.empty()) { continue; // Missing + } // Get the options for this species Options& from_options = alloptions[from]; @@ -195,11 +196,11 @@ void Recycling::transform_impl(GuardedOptions& state) { // Get metric tensor components Coordinates* coord = mesh->getCoordinates(); - const Field2D& J = coord->J; - const Field2D& dy = coord->dy; - const Field2D& dx = coord->dx; - const Field2D& dz = coord->dz; - const Field2D& g_22 = coord->g_22; + const auto& J = coord->J; + const auto& dy = coord->dy; + const auto& dx = coord->dx; + const auto& dz = coord->dz; + const auto& g_22 = coord->g_22; for (auto& channel : channels) { const GuardedOptions species_from = state["species"][channel.from]; @@ -464,8 +465,8 @@ void Recycling::transform_impl(GuardedOptions& state) { for (int iz = 0; iz < mesh->LocalNz; iz++) { // Volume of cell adjacent to wall which will receive source - BoutReal volume = J(mesh->xend, iy) * dx(mesh->xend, iy) * dy(mesh->xend, iy) - * dz(mesh->xend, iy); + BoutReal volume = J(mesh->xend, iy, iz) * dx(mesh->xend, iy, iz) + * dy(mesh->xend, iy, iz) * dz(mesh->xend, iy, iz); // If cell is a pump, overwrite multiplier with pump multiplier BoutReal multiplier = channel.sol_multiplier; @@ -587,8 +588,8 @@ void Recycling::transform_impl(GuardedOptions& state) { for (int iz = 0; iz < mesh->LocalNz; iz++) { // Volume of cell adjacent to wall which will receive source - BoutReal volume = J(mesh->xstart, iy) * dx(mesh->xstart, iy) - * dy(mesh->xstart, iy) * dz(mesh->xstart, iy); + BoutReal volume = J(mesh->xstart, iy, iz) * dx(mesh->xstart, iy, iz) + * dy(mesh->xstart, iy, iz) * dz(mesh->xstart, iy, iz); // If cell is a pump, overwrite multiplier with pump multiplier BoutReal multiplier = channel.pfr_multiplier; From 566fce0fd0b2998e0d70d93b869225e527f3bd2c Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 23 Jun 2026 14:02:35 +0200 Subject: [PATCH 20/24] Fix SNBConduction for 3D compilation --- src/snb_conduction.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/snb_conduction.cxx b/src/snb_conduction.cxx index 562037caa..0bce5f876 100644 --- a/src/snb_conduction.cxx +++ b/src/snb_conduction.cxx @@ -1,5 +1,5 @@ -#include #include "../include/snb_conduction.hxx" +#include #include using bout::globals::mesh; @@ -13,7 +13,7 @@ void SNBConduction::transform_impl(GuardedOptions& state) { // SNB non-local heat flux. Also returns the Spitzer-Harm value for comparison // Note: Te in eV, Ne in Nnorm - Field2D dy_orig = mesh->getCoordinates()->dy; + auto dy_orig = mesh->getCoordinates()->dy; mesh->getCoordinates()->dy *= rho_s0; // Convert distances to m // Inputs in eV and m^-3 From 715964b3de3df9e63b4f226d94bab41f816c253b Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 23 Jun 2026 14:02:52 +0200 Subject: [PATCH 21/24] Fix tests for 3D compilation --- tests/unit/fake_mesh.hxx | 24 +++++++++-------- tests/unit/test_anomalous_diffusion.cxx | 34 +++++++++++++------------ 2 files changed, 32 insertions(+), 26 deletions(-) diff --git a/tests/unit/fake_mesh.hxx b/tests/unit/fake_mesh.hxx index d3fb0ee61..a43d681d5 100644 --- a/tests/unit/fake_mesh.hxx +++ b/tests/unit/fake_mesh.hxx @@ -1,6 +1,7 @@ #pragma once #include +#include #include #include #include @@ -19,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -95,10 +97,11 @@ public: // Use this if the FakeMesh needs x- and y-boundaries void createBoundaries() { - addBoundary(new BoundaryRegionXIn("core", ystart, yend, this)); - addBoundary(new BoundaryRegionXOut("sol", ystart, yend, this)); - addBoundary(new BoundaryRegionYUp("upper_target", xstart, xend, this)); - addBoundary(new BoundaryRegionYDown("lower_target", xstart, xend, this)); + addBoundary(bout::boundary::NewBoundaryRegionXIn("core", ystart, yend, this)); + addBoundary(bout::boundary::NewBoundaryRegionXOut("sol", ystart, yend, this)); + addBoundary(bout::boundary::NewBoundaryRegionYUp("upper_target", xstart, xend, this)); + addBoundary( + bout::boundary::NewBoundaryRegionYDown("lower_target", xstart, xend, this)); } comm_handle send(FieldGroup& UNUSED(g)) override { return nullptr; } @@ -121,6 +124,7 @@ public: [[maybe_unused]] int Z) const override { return 0; } + std::set getPossibleBoundaries() const override { return {}; } bool firstX() const override { return true; } bool lastX() const override { return true; } int sendXOut(BoutReal* UNUSED(buffer), int UNUSED(size), int UNUSED(tag)) override { @@ -171,11 +175,11 @@ public: RangeIterator iterateBndryUpperInnerY() const override { return RangeIterator(); } bool hasBndryLowerY() const override { return false; } bool hasBndryUpperY() const override { return false; } - void addBoundary(BoundaryRegion* region) override { boundaries.push_back(region); } - std::vector getBoundaries() override { return boundaries; } - std::vector> - getBoundariesPar(BoundaryParType UNUSED(type)) override { - return std::vector>(); + void addBoundary(BoundaryRegionBase* region) override { boundaries.push_back(region); } + std::vector getBoundaries() override { return boundaries; } + std::vector> + getBoundariesPar(BoundaryParType UNUSED(type)) const override { + return std::vector>(); } BoutReal GlobalX(int jx) const override { return jx; } BoutReal GlobalY(int jy) const override { return jy; } @@ -264,7 +268,7 @@ public: using Mesh::msg_len; private: - std::vector boundaries; + std::vector boundaries; }; /// FakeGridDataSource provides a non-null GridDataSource* source to use with FakeMesh, to diff --git a/tests/unit/test_anomalous_diffusion.cxx b/tests/unit/test_anomalous_diffusion.cxx index 197a9c22b..c642175ab 100644 --- a/tests/unit/test_anomalous_diffusion.cxx +++ b/tests/unit/test_anomalous_diffusion.cxx @@ -1,22 +1,22 @@ #include "gtest/gtest.h" -#include "test_extras.hxx" // FakeMesh #include "fake_mesh_fixture.hxx" +#include "test_extras.hxx" // FakeMesh #include "../../include/anomalous_diffusion.hxx" /// Global mesh -namespace bout{ -namespace globals{ -extern Mesh *mesh; +namespace bout { +namespace globals { +extern Mesh* mesh; } // namespace globals } // namespace bout // The unit tests use the global mesh using namespace bout::globals; -#include // For generating functions +#include // For generating functions // Reuse the "standard" fixture for FakeMesh using AnomalousDiffusionTest = FakeMeshFixture; @@ -38,10 +38,10 @@ TEST_F(AnomalousDiffusionTest, NoDiffusion) { Field3D N = FieldFactory::get()->create3D("1 + y * (x - 0.5)", &options, mesh); mesh->communicate(N); - + Options state; state["species"]["h"]["density"] = N; - + // If D is not set, then the diffusion should not be calculated component.transform(state); @@ -51,23 +51,23 @@ TEST_F(AnomalousDiffusionTest, NoDiffusion) { } TEST_F(AnomalousDiffusionTest, ParticleDiffusion) { - - Coordinates *coords = mesh->getCoordinates(); + + Coordinates* coords = mesh->getCoordinates(); coords->Bxy = 1.0; // Note: This is non-finite or zero? - + Options options; options["units"]["meters"] = 1.0; options["units"]["seconds"] = 1.0; - options["h"]["anomalous_D"] = 1.0; // Set particle diffusion for "h" species + options["h"]["anomalous_D"] = 1.0; // Set particle diffusion for "h" species AnomalousDiffusion component("h", options, nullptr); Options state; state["species"]["h"]["density"] = - FieldFactory::get()->create3D("1 + y * (x - 0.5)", &options, mesh); + FieldFactory::get()->create3D("1 + y * (x - 0.5)", &options, mesh); state["species"]["h"]["AA"] = 1.0; // Atomic mass number - + component.transform(state); // Expect all sources to be set @@ -82,13 +82,15 @@ TEST_F(AnomalousDiffusionTest, ParticleDiffusion) { "RGN_NOBNDRY")); // Expect the sum over all cells of density source to be zero - Field2D dV = coords->J * coords->dx * coords->dy * coords->dz; // Cell volume + auto dV = coords->J * coords->dx * coords->dy * coords->dz; // Cell volume Field3D source = get(state["species"]["h"]["density_source"]); BoutReal integral = 0.0; BOUT_FOR_SERIAL(i, source.getRegion("RGN_NOBNDRY")) { - ASSERT_TRUE(std::isfinite(dV[i])) << "Volume element not finite at " << i.x() << ", " << i.y() << ", " << i.z(); - ASSERT_TRUE(std::isfinite(source[i])) << "Density source not finite at " << i.x() << ", " << i.y() << ", " << i.z(); + ASSERT_TRUE(std::isfinite(dV[i])) + << "Volume element not finite at " << i.x() << ", " << i.y() << ", " << i.z(); + ASSERT_TRUE(std::isfinite(source[i])) + << "Density source not finite at " << i.x() << ", " << i.y() << ", " << i.z(); integral += source[i] * dV[i]; } ASSERT_LT(abs(integral), 1e-3) << "Integral of density source should be close to zero"; From 181c93717a5d78cb9c756c4b11fc3ad953e220a4 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 23 Jun 2026 14:34:00 +0200 Subject: [PATCH 22/24] Fixup for 2D metrics --- src/div_ops.cxx | 41 +++++++++++++++++------------------------ 1 file changed, 17 insertions(+), 24 deletions(-) diff --git a/src/div_ops.cxx b/src/div_ops.cxx index f28f99a5f..c901cbc76 100644 --- a/src/div_ops.cxx +++ b/src/div_ops.cxx @@ -634,46 +634,39 @@ const Field2D Laplace_FV(const Field2D& k, const Field2D& f) { // Calculate gradients on cell faces - BoutReal gR = (coord->g11(i, j, k) + coord->g11(i + 1, j, k)) - * (f(i + 1, j, k) - f(i, j, k)) - / (coord->dx(i + 1, j, k) + coord->dx(i, j, k)); + BoutReal gR = (coord->g11(i, j) + coord->g11(i + 1, j)) * (f(i + 1, j) - f(i, j)) + / (coord->dx(i + 1, j) + coord->dx(i, j)); - BoutReal gL = (coord->g11(i - 1, j, k) + coord->g11(i, j, k)) - * (f(i, j, k) - f(i - 1, j, k)) - / (coord->dx(i - 1, j, k) + coord->dx(i, j, k)); + BoutReal gL = (coord->g11(i - 1, j) + coord->g11(i, j)) * (f(i, j) - f(i - 1, j)) + / (coord->dx(i - 1, j) + coord->dx(i, j)); - BoutReal gU = (coord->g22(i, j, k) + coord->g22(i, j + 1)) - * (f(i, j + 1) - f(i, j, k)) - / (coord->dy(i, j + 1) + coord->dy(i, j, k)); + BoutReal gU = (coord->g22(i, j) + coord->g22(i, j + 1)) * (f(i, j + 1) - f(i, j)) + / (coord->dy(i, j + 1) + coord->dy(i, j)); - BoutReal gD = (coord->g22(i, j - 1) + coord->g22(i, j, k)) - * (f(i, j, k) - f(i, j - 1)) - / (coord->dy(i, j, k) + coord->dy(i, j - 1)); + BoutReal gD = (coord->g22(i, j - 1) + coord->g22(i, j)) * (f(i, j) - f(i, j - 1)) + / (coord->dy(i, j) + coord->dy(i, j - 1)); // Flow right - BoutReal flux = gR * 0.25 * (coord->J(i + 1, j, k) + coord->J(i, j, k)) - * (k(i + 1, j, k) + k(i, j, k)); + BoutReal flux = + gR * 0.25 * (coord->J(i + 1, j) + coord->J(i, j)) * (k(i + 1, j) + k(i, j)); - result(i, j, k) = flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i, j) = flux / (coord->dx(i, j) * coord->J(i, j)); // Flow left - flux = gL * 0.25 * (coord->J(i - 1, j, k) + coord->J(i, j, k)) - * (k(i - 1, j, k) + k(i, j, k)); - result(i, j, k) -= flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + flux = gL * 0.25 * (coord->J(i - 1, j) + coord->J(i, j)) * (k(i - 1, j) + k(i, j)); + result(i, j) -= flux / (coord->dx(i, j) * coord->J(i, j)); // Flow up - flux = gU * 0.25 * (coord->J(i, j + 1) + coord->J(i, j, k)) - * (k(i, j + 1) + k(i, j, k)); - result(i, j, k) += flux / (coord->dy(i, j, k) * coord->J(i, j, k)); + flux = gU * 0.25 * (coord->J(i, j + 1) + coord->J(i, j)) * (k(i, j + 1) + k(i, j)); + result(i, j) += flux / (coord->dy(i, j) * coord->J(i, j)); // Flow down - flux = gD * 0.25 * (coord->J(i, j - 1) + coord->J(i, j, k)) - * (k(i, j - 1) + k(i, j, k)); - result(i, j, k) -= flux / (coord->dy(i, j, k) * coord->J(i, j, k)); + flux = gD * 0.25 * (coord->J(i, j - 1) + coord->J(i, j)) * (k(i, j - 1) + k(i, j)); + result(i, j) -= flux / (coord->dy(i, j) * coord->J(i, j)); } } return result; From 9784bd1a5870d59085b2b33d92a49d2c97fd72be Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 23 Jun 2026 14:36:32 +0200 Subject: [PATCH 23/24] CI: use configure_options from build_type section --- .github/workflows/tests.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index a9bc176f5..a63bfe862 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -22,7 +22,7 @@ jobs: OMPI_MCA_rmaps_base_oversubscribe: yes MPIRUN: mpiexec -np LD_LIBRARY_PATH: /home/runner/local/lib:$LD_LIBRARY_PATH - BOUT_CONFIGURE_OPTIONS: ${{ matrix.config.configure_options }} -DHERMES_COVERAGE=${{ matrix.build_type.coverage }} + BOUT_CONFIGURE_OPTIONS: ${{ matrix.build_type.configure_options }} -DHERMES_COVERAGE=${{ matrix.build_type.coverage }} BUILD_TYPE: ${{ matrix.build_type.cmake_build_type }} TEST_TYPE: ${{ matrix.build_type.test_type }} strategy: @@ -49,7 +49,7 @@ jobs: - name: Job information run: | echo Build: ${{ matrix.config.name }}, ${{ matrix.config.os }} - echo Configure options: ${{ matrix.config.configure_options }} + echo Configure options: ${{ matrix.build_type.configure_options }} echo Build type: ${{ matrix.build_type.cmake_build_type }} echo Coverage: ${{ matrix.build_type.coverage }} echo Tests to run: ${{ matrix.build_type.test_type }} From 8ddace8d5a22c9ea5bcff5ca4bc7e9be9cf2022d Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 23 Jun 2026 14:38:57 +0200 Subject: [PATCH 24/24] CI: Fix typo --- .github/workflows/tests.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index a63bfe862..a2e03161e 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -42,7 +42,7 @@ jobs: configure_options: "-DPACKAGE_TESTS=OFF -DHERMES_ERROR_ON_WARNINGS=ON" - cmake_build_type: Debug test_type: None - configure_options: "-PACKAGE_TESTS=OFF -DHERMES_ERROR_ON_WARNINGS=ON -DBOUT_ENABLE_METRIC_3D=ON" + configure_options: "-DPACKAGE_TESTS=OFF -DHERMES_ERROR_ON_WARNINGS=ON -DBOUT_ENABLE_METRIC_3D=ON" coverage: OFF steps: