From a1d1e79fbee373e31e520fdc3a98313d1197620d Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 1 Mar 2024 19:34:48 -0800 Subject: [PATCH 01/19] neutral_mixed: Add FV::Div_a_Grad_perp_limit operator Apply limiter/upwinding to both X and Y advection --- include/div_ops.hxx | 248 +++++++++++++++++++++++++++++++++++++++--- src/neutral_mixed.cxx | 99 +++++++++-------- 2 files changed, 286 insertions(+), 61 deletions(-) diff --git a/include/div_ops.hxx b/include/div_ops.hxx index 2b359190..3f74677b 100644 --- a/include/div_ops.hxx +++ b/include/div_ops.hxx @@ -23,12 +23,12 @@ */ -#ifndef __DIV_OPS_H__ -#define __DIV_OPS_H__ +#ifndef HERMES_DIV_OPS_H +#define HERMES_DIV_OPS_H #include -#include #include +#include /*! * Diffusion in index space @@ -58,7 +58,8 @@ const Field2D Laplace_FV(const Field2D& k, const Field2D& f); /// Perpendicular diffusion including X and Y directions const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f); /// Version of function that returns flows -const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, Field3D& flux_xlow, Field3D& flux_ylow); +const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, + Field3D& flux_xlow, Field3D& flux_ylow); namespace FV { @@ -100,8 +101,7 @@ struct Superbee { BoutReal sign = SIGN(gL); gL = fabs(gL); gR = fabs(gR); - BoutReal half_slope = sign * BOUTMAX(BOUTMIN(gL, 0.5*gR), - BOUTMIN(gR, 0.5*gL)); + 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; } @@ -210,10 +210,9 @@ const Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, flux = bndryval * vpar * vpar; } else { // Add flux due to difference in boundary values - flux = - s.R * vpar * sv.R - + BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j + 1, k))) - * (s.R * sv.R - bndryval * vpar); + flux = s.R * vpar * sv.R + + BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j + 1, k))) + * (s.R * sv.R - bndryval * vpar); } } else { // Maximum wave speed in the two cells @@ -239,10 +238,9 @@ const Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, flux = bndryval * vpar * vpar; } else { // Add flux due to difference in boundary values - flux = - s.L * vpar * sv.L - - BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j - 1, k))) - * (s.L * sv.L - bndryval * vpar); + flux = s.L * vpar * sv.L + - BOUTMAX(wave_speed(i, j, k), fabs(v(i, j, k)), fabs(v(i, j - 1, k))) + * (s.L * sv.L - bndryval * vpar); } } else { // Maximum wave speed in the two cells @@ -454,6 +452,226 @@ const Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; } +/// Div ( a Grad_perp(f) ) -- diffusion +/// +/// This version uses a slope limiter to calculate cell edge values in X, +/// the advects the upwind cell edge. +/// +/// 1st order upwinding is used in Y. +template +const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& f) { + ASSERT2(a.getLocation() == f.getLocation()); + + Mesh* mesh = a.getMesh(); + + // Requires at least 2 communication guard cells in X, 1 in Y + ASSERT1(mesh->xstart >= 2); + ASSERT1(mesh->ystart >= 1); + + CellEdges cellboundary; + + Field3D result{zeroFrom(f)}; + + Coordinates* coord = f.getCoordinates(); + + // 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 = f(i + 1, j, k) - f(i, j, k); + + BoutReal aedge; // 'a' at the cell edge + if (((i == xs) and mesh->firstX()) or ((i == xe) and mesh->lastX())) { + // Mid-point average boundary value + aedge = 0.5 * (a(i + 1, j, k) + a(i, j, k)); + } else if (gradient > 0) { + // Flux is from (i+1) to (i) + // Reconstruct `a` at left of (i+1, j, k) + + Stencil1D sa; + sa.m = a(i, j, k); + sa.c = a(i + 1, j, k); + sa.p = a(i + 2, j, k); + cellboundary(sa); // Calculate sa.R and sa.L + + aedge = sa.L; + } else { + // Flux is from (i) to (i+1) + // Reconstruct `a` at right of (i, j, k) + + Stencil1D sa; + sa.m = a(i - 1, j, k); + sa.c = a(i, j, k); + sa.p = a(i + 1, j, k); + cellboundary(sa); // Calculate sa.R and sa.L + + aedge = sa.R; + } + + // Flux across cell edge + const BoutReal fout = gradient * aedge + * (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)); + + 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)); + } + } + } + + // 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); + } + + // 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; + + // 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 aedge; + if ((j == mesh->yend) and mesh->lastY(i)) { + // Midpoint boundary value + aedge = 0.5 * (ac(i, j, k) + aup(i, j + 1, k)); + } else if (dfdy > 0) { + // Flux from (j+1) to (j) + aedge = aup(i, j + 1, k); + } else { + // Flux from (j) to (j+1) + aedge = ac(i, j, k); + } + + BoutReal fout = aedge * 0.5 + * (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)); + + if ((j == mesh->ystart) and mesh->firstY(i)) { + aedge = 0.5 * (ac(i, j, k) + adown(i, j - 1, k)); + } else if (dfdy > 0) { + aedge = ac(i, j, k); + } else { + aedge = adown(i, j - 1, k); + } + + fout = aedge * 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)); + } + } + } + + // 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; + + 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)); + + 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); + } + } + } + // Check if we need to transform back + if (f.hasParallelSlices() && a.hasParallelSlices()) { + result += yzresult; + } else { + result += fromFieldAligned(yzresult); + } + + return result; +} + } // namespace FV -#endif // __DIV_OPS_H__ +#endif // HERMES_DIV_OPS_H diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index a2f848a9..c7d1081d 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -1,13 +1,13 @@ #include -#include #include #include +#include #include #include "../include/div_ops.hxx" -#include "../include/neutral_mixed.hxx" #include "../include/hermes_build_config.hxx" +#include "../include/neutral_mixed.hxx" using bout::globals::mesh; @@ -31,16 +31,16 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* // Evolving variables e.g name is "h" or "h+" solver->add(Nn, std::string("N") + name); solver->add(Pn, std::string("P") + name); - evolve_momentum = options["evolve_momentum"] - .doc("Evolve parallel neutral momentum?") - .withDefault(true); + .doc("Evolve parallel neutral momentum?") + .withDefault(true); if (evolve_momentum) { solver->add(NVn, std::string("NV") + name); } else { - output_warn.write("WARNING: Not evolving neutral parallel momentum. NVn and Vn set to zero\n"); + output_warn.write( + "WARNING: Not evolving neutral parallel momentum. NVn and Vn set to zero\n"); NVn = 0.0; Vn = 0.0; } @@ -62,18 +62,19 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Enable preconditioning in neutral model?") .withDefault(true); - flux_limit = options["flux_limit"] - .doc("Limit diffusive fluxes to fraction of thermal speed. <0 means off.") - .withDefault(0.2); + flux_limit = + options["flux_limit"] + .doc("Limit diffusive fluxes to fraction of thermal speed. <0 means off.") + .withDefault(0.2); diffusion_limit = options["diffusion_limit"] - .doc("Upper limit on diffusion coefficient [m^2/s]. <0 means off") - .withDefault(-1.0) - / (meters * meters / seconds); // Normalise + .doc("Upper limit on diffusion coefficient [m^2/s]. <0 means off") + .withDefault(-1.0) + / (meters * meters / seconds); // Normalise neutral_viscosity = options["neutral_viscosity"] - .doc("Include neutral gas viscosity?") - .withDefault(true); + .doc("Include neutral gas viscosity?") + .withDefault(true); if (precondition) { inv = std::unique_ptr(Laplacian::create(&options["precon_laplace"])); @@ -98,10 +99,11 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* density_source = 0.0; mesh->get(density_source, std::string("N") + name + "_src"); // Allow the user to override the source - density_source = alloptions[std::string("N") + name]["source"] - .doc("Source term in ddt(N" + name + std::string("). Units [m^-3/s]")) - .withDefault(density_source) - / (Nnorm * Omega_ci); + density_source = + alloptions[std::string("N") + name]["source"] + .doc("Source term in ddt(N" + name + std::string("). Units [m^-3/s]")) + .withDefault(density_source) + / (Nnorm * Omega_ci); // Try to read the pressure source from the mesh // Units of Pascals per second @@ -109,18 +111,22 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* mesh->get(pressure_source, std::string("P") + name + "_src"); // Allow the user to override the source pressure_source = alloptions[std::string("P") + name]["source"] - .doc(std::string("Source term in ddt(P") + name - + std::string("). Units [N/m^2/s]")) - .withDefault(pressure_source) - / (SI::qe * Nnorm * Tnorm * Omega_ci); + .doc(std::string("Source term in ddt(P") + name + + std::string("). Units [N/m^2/s]")) + .withDefault(pressure_source) + / (SI::qe * Nnorm * Tnorm * Omega_ci); // Set boundary condition defaults: Neumann for all but the diffusivity. // The dirichlet on diffusivity ensures no radial flux. // NV and V are ignored as they are hardcoded in the parallel BC code. - alloptions[std::string("Dnn") + name]["bndry_all"] = alloptions[std::string("Dnn") + name]["bndry_all"].withDefault("dirichlet"); - alloptions[std::string("T") + name]["bndry_all"] = alloptions[std::string("T") + name]["bndry_all"].withDefault("neumann"); - alloptions[std::string("P") + name]["bndry_all"] = alloptions[std::string("P") + name]["bndry_all"].withDefault("neumann"); - alloptions[std::string("N") + name]["bndry_all"] = alloptions[std::string("N") + name]["bndry_all"].withDefault("neumann"); + alloptions[std::string("Dnn") + name]["bndry_all"] = + alloptions[std::string("Dnn") + name]["bndry_all"].withDefault("dirichlet"); + alloptions[std::string("T") + name]["bndry_all"] = + alloptions[std::string("T") + name]["bndry_all"].withDefault("neumann"); + alloptions[std::string("P") + name]["bndry_all"] = + alloptions[std::string("P") + name]["bndry_all"].withDefault("neumann"); + alloptions[std::string("N") + name]["bndry_all"] = + alloptions[std::string("N") + name]["bndry_all"].withDefault("neumann"); // Pick up BCs from input file Dnn.setBoundary(std::string("Dnn") + name); @@ -134,12 +140,12 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* logPnlim.setBoundary(std::string("P") + name); Nnlim.setBoundary(std::string("N") + name); - // Product of Dnn and another parameter has same BC as Dnn - see eqns to see why this is necessary + // Product of Dnn and another parameter has same BC as Dnn - see eqns to see why this is + // necessary DnnNn.setBoundary(std::string("Dnn") + name); DnnPn.setBoundary(std::string("Dnn") + name); DnnTn.setBoundary(std::string("Dnn") + name); DnnNVn.setBoundary(std::string("Dnn") + name); - } void NeutralMixed::transform(Options& state) { @@ -264,7 +270,8 @@ void NeutralMixed::finally(const Options& state) { BoutReal neutral_lmax = 0.1 / get(state["units"]["meters"]); // Normalised length - Field3D Rnn = sqrt(Tn / AA) / neutral_lmax; // Neutral-neutral collisions [normalised frequency] + Field3D Rnn = + sqrt(Tn / AA) / neutral_lmax; // Neutral-neutral collisions [normalised frequency] if (localstate.isSet("collision_frequency")) { // Dnn = Vth^2 / sigma @@ -276,11 +283,8 @@ void NeutralMixed::finally(const Options& state) { if (flux_limit > 0.0) { // Apply flux limit to diffusion, // using the local thermal speed and pressure gradient magnitude - Field3D Dmax = flux_limit * sqrt(Tn / AA) / - (abs(Grad(logPnlim)) + 1. / neutral_lmax); - BOUT_FOR(i, Dmax.getRegion("RGN_NOBNDRY")) { - Dnn[i] = BOUTMIN(Dnn[i], Dmax[i]); - } + Field3D Dmax = flux_limit * sqrt(Tn / AA) / (abs(Grad(logPnlim)) + 1. / neutral_lmax); + BOUT_FOR(i, Dmax.getRegion("RGN_NOBNDRY")) { Dnn[i] = BOUTMIN(Dnn[i], Dmax[i]); } } if (diffusion_limit > 0.0) { @@ -330,8 +334,9 @@ void NeutralMixed::finally(const Options& state) { ///////////////////////////////////////////////////// // Neutral density TRACE("Neutral density"); - ddt(Nn) = -FV::Div_par_mod(Nn, Vn, sound_speed) // Advection - + FV::Div_a_Grad_perp(DnnNn, logPnlim) // Perpendicular diffusion + ddt(Nn) = -FV::Div_par_mod(Nn, Vn, sound_speed) // Parallel advection + + FV::Div_a_Grad_perp_limit( + DnnNn, logPnlim) // Perpendicular advection ; Sn = density_source; // Save for possible output @@ -348,8 +353,9 @@ void NeutralMixed::finally(const Options& state) { ddt(NVn) = -AA * FV::Div_par_fvv(Nnlim, Vn, sound_speed) // Momentum flow - - Grad_par(Pn) // Pressure gradient - + FV::Div_a_Grad_perp(DnnNVn, logPnlim) // Perpendicular diffusion + - Grad_par(Pn) // Pressure gradient + + FV::Div_a_Grad_perp_limit(DnnNVn, + logPnlim) // Perpendicular advection ; if (neutral_viscosity) { @@ -363,9 +369,10 @@ void NeutralMixed::finally(const Options& state) { // eta_n = (2. / 5) * kappa_n; // - ddt(NVn) += AA * FV::Div_a_Grad_perp((2. / 5) * DnnNn, Vn) // Perpendicular viscosity - + AA * FV::Div_par_K_Grad_par((2. / 5) * DnnNn, Vn) // Parallel viscosity - ; + ddt(NVn) += + AA * FV::Div_a_Grad_perp((2. / 5) * DnnNn, Vn) // Perpendicular viscosity + + AA * FV::Div_par_K_Grad_par((2. / 5) * DnnNn, Vn) // Parallel viscosity + ; } if (localstate.isSet("momentum_source")) { @@ -382,11 +389,12 @@ void NeutralMixed::finally(const Options& state) { // Neutral pressure TRACE("Neutral pressure"); - ddt(Pn) = -FV::Div_par_mod(Pn, Vn, sound_speed) // Advection + ddt(Pn) = -FV::Div_par_mod(Pn, Vn, sound_speed) // Parallel advection - (2. / 3) * Pn * Div_par(Vn) // Compression - + FV::Div_a_Grad_perp(DnnPn, logPnlim) // Perpendicular diffusion - + FV::Div_a_Grad_perp(DnnNn, Tn) // Conduction - + FV::Div_par_K_Grad_par(DnnNn, Tn) // Parallel conduction + + FV::Div_a_Grad_perp_limit( + DnnPn, logPnlim) // Perpendicular advection + + FV::Div_a_Grad_perp(DnnNn, Tn) // Perpendicular conduction + + FV::Div_par_K_Grad_par(DnnNn, Tn) // Parallel conduction ; Sp = pressure_source; @@ -531,7 +539,6 @@ void NeutralMixed::outputVars(Options& state) { {"long_name", name + " pressure source"}, {"species", name}, {"source", "neutral_mixed"}}); - } } From 9c2e8719495318a8b4bb39f696f52bc2ac74831e Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 5 Mar 2024 23:01:30 -0800 Subject: [PATCH 02/19] neutral_mixed: Switch perpendicular limiter to Upwind Seems to work better for neutral perpendicular diffusion than MC. --- src/neutral_mixed.cxx | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index c7d1081d..1b4a172d 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -11,6 +11,10 @@ using bout::globals::mesh; +/// The limiter method in the radial pressure-diffusion. +/// Upwind is consistent with the Y (poloidal) advection. +using PerpLimiter = FV::Upwind; + NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* solver) : name(name) { AUTO_TRACE(); @@ -335,7 +339,7 @@ void NeutralMixed::finally(const Options& state) { // Neutral density TRACE("Neutral density"); ddt(Nn) = -FV::Div_par_mod(Nn, Vn, sound_speed) // Parallel advection - + FV::Div_a_Grad_perp_limit( + + FV::Div_a_Grad_perp_limit( DnnNn, logPnlim) // Perpendicular advection ; @@ -354,8 +358,8 @@ void NeutralMixed::finally(const Options& state) { ddt(NVn) = -AA * FV::Div_par_fvv(Nnlim, Vn, sound_speed) // Momentum flow - Grad_par(Pn) // Pressure gradient - + FV::Div_a_Grad_perp_limit(DnnNVn, - logPnlim) // Perpendicular advection + + FV::Div_a_Grad_perp_limit(DnnNVn, + logPnlim) // Perpendicular advection ; if (neutral_viscosity) { @@ -391,7 +395,7 @@ void NeutralMixed::finally(const Options& state) { ddt(Pn) = -FV::Div_par_mod(Pn, Vn, sound_speed) // Parallel advection - (2. / 3) * Pn * Div_par(Vn) // Compression - + FV::Div_a_Grad_perp_limit( + + FV::Div_a_Grad_perp_limit( DnnPn, logPnlim) // Perpendicular advection + FV::Div_a_Grad_perp(DnnNn, Tn) // Perpendicular conduction + FV::Div_par_K_Grad_par(DnnNn, Tn) // Parallel conduction From bcfb1640587494e2b3df6d7833f98c2d87fafc7b Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 26 Mar 2024 22:09:17 -0700 Subject: [PATCH 03/19] neutral_mixed: Modify Div_a_Grad_perp_limit operator Take field (Nn, Pn, NVn) and Dnn separately; only upwind the field, and use cell face average values for Dnn. Also add neutral_conduction switch. --- include/div_ops.hxx | 83 ++++++++++++++++++++++----------------- include/neutral_mixed.hxx | 3 +- src/neutral_mixed.cxx | 47 ++++++++++------------ 3 files changed, 69 insertions(+), 64 deletions(-) diff --git a/include/div_ops.hxx b/include/div_ops.hxx index 3f74677b..fb372a54 100644 --- a/include/div_ops.hxx +++ b/include/div_ops.hxx @@ -452,14 +452,14 @@ const Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; } -/// Div ( a Grad_perp(f) ) -- diffusion +/// Div ( a g Grad_perp(f) ) -- Perpendicular gradient-driven advection /// -/// This version uses a slope limiter to calculate cell edge values in X, +/// This version uses a slope limiter to calculate cell edge values of g in X, /// the advects the upwind cell edge. /// /// 1st order upwinding is used in Y. template -const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& f) { +const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D& f) { ASSERT2(a.getLocation() == f.getLocation()); Mesh* mesh = a.getMesh(); @@ -476,46 +476,46 @@ const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& f) { // Flux in x - int xs = mesh->xstart - 1; - int xe = mesh->xend; - - for (int i = xs; i <= xe; i++) { + 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 const BoutReal gradient = f(i + 1, j, k) - f(i, j, k); - BoutReal aedge; // 'a' at the cell edge - if (((i == xs) and mesh->firstX()) or ((i == xe) and mesh->lastX())) { - // Mid-point average boundary value - aedge = 0.5 * (a(i + 1, j, k) + a(i, j, k)); + // Mid-point average boundary value of 'a' + const BoutReal aedge = 0.5 * (a(i + 1, j, k) + a(i, j, k)); + BoutReal gedge; + if (((i == mesh->xstart - 1) and mesh->firstX()) or + ((i == mesh->xend) and mesh->lastX())) { + // Mid-point average boundary value of 'g' + gedge = 0.5 * (g(i + 1, j, k) + g(i, j, k)); } else if (gradient > 0) { // Flux is from (i+1) to (i) - // Reconstruct `a` at left of (i+1, j, k) + // Reconstruct `g` at left of (i+1, j, k) - Stencil1D sa; - sa.m = a(i, j, k); - sa.c = a(i + 1, j, k); - sa.p = a(i + 2, j, k); - cellboundary(sa); // Calculate sa.R and sa.L + Stencil1D sg; + sg.m = g(i, j, k); + sg.c = g(i + 1, j, k); + sg.p = g(i + 2, j, k); + cellboundary(sg); // Calculate sg.R and sg.L - aedge = sa.L; + gedge = sg.L; } else { // Flux is from (i) to (i+1) - // Reconstruct `a` at right of (i, j, k) + // Reconstruct `g` at right of (i, j, k) - Stencil1D sa; - sa.m = a(i - 1, j, k); - sa.c = a(i, j, k); - sa.p = a(i + 1, j, k); - cellboundary(sa); // Calculate sa.R and sa.L + Stencil1D sg; + sg.m = g(i - 1, j, k); + sg.c = g(i, j, k); + sg.p = g(i + 1, j, k); + cellboundary(sg); // Calculate sg.R and sg.L - aedge = sa.R; + gedge = sg.R; } // Flux across cell edge - const BoutReal fout = gradient * aedge + 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)); @@ -531,11 +531,13 @@ const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& f) { // Fields containing values along the magnetic field Field3D fup(mesh), fdown(mesh); Field3D aup(mesh), adown(mesh); + Field3D gup(mesh), gdown(mesh); // Values on this y slice (centre). // This is needed because toFieldAligned may modify the field - Field3D fc = f; Field3D ac = a; + Field3D gc = g; + Field3D fc = f; // Result of the Y and Z fluxes Field3D yzresult(mesh); @@ -549,12 +551,16 @@ const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& f) { aup = a.yup(); adown = a.ydown(); + + gup = g.yup(); + gdown = g.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); + gup = gdown = gc = toFieldAligned(g); yzresult.setDirectionY(YDirectionType::Aligned); } @@ -587,19 +593,20 @@ const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& f) { BoutReal dfdy = 2. * (fup(i, j + 1, k) - fc(i, j, k)) / (coord->dy(i, j + 1) + coord->dy(i, j)); - BoutReal aedge; + 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 - aedge = 0.5 * (ac(i, j, k) + aup(i, j + 1, k)); + gedge = 0.5 * (gc(i, j, k) + gup(i, j + 1, k)); } else if (dfdy > 0) { // Flux from (j+1) to (j) - aedge = aup(i, j + 1, k); + gedge = gup(i, j + 1, k); } else { // Flux from (j) to (j+1) - aedge = ac(i, j, k); + gedge = gc(i, j, k); } - BoutReal fout = aedge * 0.5 + 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); @@ -614,15 +621,16 @@ const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& f) { dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k)) / (coord->dy(i, j) + coord->dy(i, j - 1)); + aedge = 0.5 * (ac(i, j, k) + adown(i, j - 1, k)); if ((j == mesh->ystart) and mesh->firstY(i)) { - aedge = 0.5 * (ac(i, j, k) + adown(i, j - 1, k)); + gedge = 0.5 * (gc(i, j, k) + gdown(i, j - 1, k)); } else if (dfdy > 0) { - aedge = ac(i, j, k); + gedge = gc(i, j, k); } else { - aedge = adown(i, j - 1, k); + gedge = gdown(i, j - 1, k); } - fout = aedge * 0.5 + 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); @@ -655,7 +663,8 @@ const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& f) { * (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 * 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); diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index f7a6d5a4..e15d3a9a 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -42,7 +42,7 @@ private: BoutReal AA; ///< Atomic mass (proton = 1) Field3D Dnn; ///< Diffusion coefficient - Field3D DnnNn, DnnPn, DnnTn, DnnNVn; ///< Used for operators + Field3D DnnNn; bool sheath_ydown, sheath_yup; @@ -52,6 +52,7 @@ private: BoutReal diffusion_limit; ///< Maximum diffusion coefficient bool neutral_viscosity; ///< include viscosity? + bool neutral_conduction; ///< Include heat conduction? bool evolve_momentum; ///< Evolve parallel momentum? bool precondition {true}; ///< Enable preconditioner? diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 1b4a172d..16a0fd12 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -80,6 +80,10 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Include neutral gas viscosity?") .withDefault(true); + neutral_conduction = options["neutral_conduction"] + .doc("Include neutral gas heat conduction?") + .withDefault(true); + if (precondition) { inv = std::unique_ptr(Laplacian::create(&options["precon_laplace"])); @@ -147,9 +151,6 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* // Product of Dnn and another parameter has same BC as Dnn - see eqns to see why this is // necessary DnnNn.setBoundary(std::string("Dnn") + name); - DnnPn.setBoundary(std::string("Dnn") + name); - DnnTn.setBoundary(std::string("Dnn") + name); - DnnNVn.setBoundary(std::string("Dnn") + name); } void NeutralMixed::transform(Options& state) { @@ -175,7 +176,7 @@ void NeutralMixed::transform(Options& state) { Vn.applyBoundary("neumann"); Vnlim.applyBoundary("neumann"); - Pnlim = floor(Nnlim * Tn, 1e-8); + Pnlim = floor(Pn, 1e-8); Pnlim.applyBoundary(); ///////////////////////////////////////////////////// @@ -275,7 +276,7 @@ void NeutralMixed::finally(const Options& state) { 0.1 / get(state["units"]["meters"]); // Normalised length Field3D Rnn = - sqrt(Tn / AA) / neutral_lmax; // Neutral-neutral collisions [normalised frequency] + sqrt(floor(Tn, 1e-5) / AA) / neutral_lmax; // Neutral-neutral collisions [normalised frequency] if (localstate.isSet("collision_frequency")) { // Dnn = Vth^2 / sigma @@ -303,20 +304,14 @@ void NeutralMixed::finally(const Options& state) { Dnn.applyBoundary(); // Neutral diffusion parameters have the same boundary condition as Dnn - DnnPn = Dnn * Pn; - DnnPn.applyBoundary(); DnnNn = Dnn * Nn; DnnNn.applyBoundary(); - Field3D DnnNVn = Dnn * NVn; - DnnNVn.applyBoundary(); if (sheath_ydown) { for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { Dnn(r.ind, mesh->ystart - 1, jz) = -Dnn(r.ind, mesh->ystart, jz); - DnnPn(r.ind, mesh->ystart - 1, jz) = -DnnPn(r.ind, mesh->ystart, jz); DnnNn(r.ind, mesh->ystart - 1, jz) = -DnnNn(r.ind, mesh->ystart, jz); - DnnNVn(r.ind, mesh->ystart - 1, jz) = -DnnNVn(r.ind, mesh->ystart, jz); } } } @@ -325,9 +320,7 @@ void NeutralMixed::finally(const Options& state) { for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { Dnn(r.ind, mesh->yend + 1, jz) = -Dnn(r.ind, mesh->yend, jz); - DnnPn(r.ind, mesh->yend + 1, jz) = -DnnPn(r.ind, mesh->yend, jz); DnnNn(r.ind, mesh->yend + 1, jz) = -DnnNn(r.ind, mesh->yend, jz); - DnnNVn(r.ind, mesh->yend + 1, jz) = -DnnNVn(r.ind, mesh->yend, jz); } } } @@ -338,10 +331,10 @@ void NeutralMixed::finally(const Options& state) { ///////////////////////////////////////////////////// // Neutral density TRACE("Neutral density"); - ddt(Nn) = -FV::Div_par_mod(Nn, Vn, sound_speed) // Parallel advection - + FV::Div_a_Grad_perp_limit( - DnnNn, logPnlim) // Perpendicular advection - ; + ddt(Nn) = + -FV::Div_par_mod(Nn, Vn, sound_speed) // Parallel advection + + FV::Div_a_Grad_perp_limit(Dnn, Nn, logPnlim) // Perpendicular advection + ; Sn = density_source; // Save for possible output if (localstate.isSet("density_source")) { @@ -358,12 +351,11 @@ void NeutralMixed::finally(const Options& state) { ddt(NVn) = -AA * FV::Div_par_fvv(Nnlim, Vn, sound_speed) // Momentum flow - Grad_par(Pn) // Pressure gradient - + FV::Div_a_Grad_perp_limit(DnnNVn, - logPnlim) // Perpendicular advection - ; + + FV::Div_a_Grad_perp_limit(Dnn, NVn, logPnlim) // Perpendicular advection + ; if (neutral_viscosity) { - // NOTE: The following viscosity terms are are not (yet) balanced + // NOTE: The following viscosity terms are not (yet) balanced // by a viscous heating term // Relationship between heat conduction and viscosity for neutral @@ -395,12 +387,15 @@ void NeutralMixed::finally(const Options& state) { ddt(Pn) = -FV::Div_par_mod(Pn, Vn, sound_speed) // Parallel advection - (2. / 3) * Pn * Div_par(Vn) // Compression - + FV::Div_a_Grad_perp_limit( - DnnPn, logPnlim) // Perpendicular advection - + FV::Div_a_Grad_perp(DnnNn, Tn) // Perpendicular conduction - + FV::Div_par_K_Grad_par(DnnNn, Tn) // Parallel conduction - ; + + FV::Div_a_Grad_perp_limit(Dnn, Pn, logPnlim) // Perpendicular advection + ; + if (neutral_conduction) { + ddt(Pn) += FV::Div_a_Grad_perp(DnnNn, Tn) // Perpendicular conduction + + FV::Div_par_K_Grad_par(DnnNn, Tn) // Parallel conduction + ; + } + Sp = pressure_source; if (localstate.isSet("energy_source")) { Sp += (2. / 3) * get(localstate["energy_source"]); From 1ec6be0a207d040536a06a4cd89cb520b2638514 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 28 Mar 2024 15:01:12 -0700 Subject: [PATCH 04/19] neutral_boundary: Fix typos in fast reflection coefficients target settings used in SOL and PFR regions, rather than SOL and PFR settings. --- src/neutral_boundary.cxx | 37 +++++++++++++++---------------------- 1 file changed, 15 insertions(+), 22 deletions(-) diff --git a/src/neutral_boundary.cxx b/src/neutral_boundary.cxx index 72e625b2..548863d5 100644 --- a/src/neutral_boundary.cxx +++ b/src/neutral_boundary.cxx @@ -175,8 +175,7 @@ void NeutralBoundary::transform(Options& state) { // Outgoing neutral heat flux [W/m^2] // This is rearranged from Power for clarity - note definition of v_th. - BoutReal q = - 2 * nnsheath * tnsheath * v_th // Incident energy + BoutReal q = 2 * nnsheath * tnsheath * v_th // Incident energy - (target_energy_refl_factor * target_fast_refl_fraction ) * 2 * nnsheath * tnsheath * v_th // Fast reflected energy - (1 - target_fast_refl_fraction) * T_FC * nnsheath * v_th; // Thermal reflected energy @@ -198,11 +197,11 @@ void NeutralBoundary::transform(Options& state) { } } - // SOL edge + // SOL edge if (sol) { - if(mesh->lastX()){ // Only do this for the processor which has the edge region - for(int iy=0; iy < mesh->LocalNy ; iy++){ - for(int iz=0; iz < mesh->LocalNz; iz++){ + if (mesh->lastX()) { // Only do this for the processor which has the edge region + for (int iy = 0; iy < mesh->LocalNy; iy++) { + for (int iz = 0; iz < mesh->LocalNz; iz++) { auto i = indexAt(Nn, mesh->xend, iy, iz); // Final domain cell auto ig = indexAt(Nn, mesh->xend+1, iy, iz); // Guard cell @@ -219,11 +218,9 @@ void NeutralBoundary::transform(Options& state) { // Outgoing neutral heat flux [W/m^2] // This is rearranged from Power for clarity - note definition of v_th. - BoutReal q = - 2 * nnsheath * tnsheath * v_th // Incident energy - - (target_energy_refl_factor * target_fast_refl_fraction ) * 2 * nnsheath * tnsheath * v_th // Fast reflected energy - - (1 - target_fast_refl_fraction) * T_FC * nnsheath * v_th; // Thermal reflected energy - + BoutReal q = 2 * nnsheath * tnsheath * v_th // Incident energy + - (sol_energy_refl_factor * sol_fast_refl_fraction ) * 2 * nnsheath * tnsheath * v_th // Fast reflected energy + - (1 - sol_fast_refl_fraction) * T_FC * nnsheath * v_th; // Thermal reflected energy // Multiply by radial cell area to get power // Expanded form of the calculation for clarity @@ -245,7 +242,6 @@ void NeutralBoundary::transform(Options& state) { // Subtract from cell next to boundary energy_source[i] -= cooling_source; wall_energy_source[i] -= cooling_source; - } } } @@ -254,12 +250,12 @@ void NeutralBoundary::transform(Options& state) { // PFR edge if (pfr) { if ((mesh->firstX()) and (!mesh->periodicY(mesh->xstart))) { // do loop if inner edge and not periodic (i.e. PFR) - for(int iy=0; iy < mesh->LocalNy ; iy++){ - for(int iz=0; iz < mesh->LocalNz; iz++){ + for (int iy = 0; iy < mesh->LocalNy ; iy++) { + for (int iz = 0; iz < mesh->LocalNz; iz++) { auto i = indexAt(Nn, mesh->xstart, iy, iz); // Final domain cell auto ig = indexAt(Nn, mesh->xstart-1, iy, iz); // Guard cell - + // Calculate midpoint values at wall const BoutReal nnsheath = 0.5 * (Nn[ig] + Nn[i]); const BoutReal tnsheath = 0.5 * (Tn[ig] + Tn[i]); @@ -267,16 +263,14 @@ void NeutralBoundary::transform(Options& state) { // Thermal speed of static Maxwellian in one direction const BoutReal v_th = 0.25 * sqrt( 8*tnsheath / (PI*AA) ); // Stangeby p.69 eqns. 2.21, 2.24 - // Approach adapted from D. Power thesis 2023 + // Approach adapted from D. Power thesis 2023 BoutReal T_FC = 3 / Tnorm; // Franck-Condon temp (hardcoded for now) // Outgoing neutral heat flux [W/m^2] // This is rearranged from Power for clarity - note definition of v_th. - BoutReal q = - 2 * nnsheath * tnsheath * v_th // Incident energy - - (target_energy_refl_factor * target_fast_refl_fraction ) * 2 * nnsheath * tnsheath * v_th // Fast reflected energy - - (1 - target_fast_refl_fraction) * T_FC * nnsheath * v_th; // Thermal reflected energy - + BoutReal q = 2 * nnsheath * tnsheath * v_th // Incident energy + - (pfr_energy_refl_factor * pfr_fast_refl_fraction ) * 2 * nnsheath * tnsheath * v_th // Fast reflected energy + - (1 - pfr_fast_refl_fraction) * T_FC * nnsheath * v_th; // Thermal reflected energy // Multiply by radial cell area to get power // Expanded form of the calculation for clarity @@ -298,7 +292,6 @@ void NeutralBoundary::transform(Options& state) { // Subtract from cell next to boundary energy_source[i] -= cooling_source; wall_energy_source[i] -= cooling_source; - } } } From fec8232e7d1dcb66b15882139214f6743466b0bc Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 28 Mar 2024 15:02:11 -0700 Subject: [PATCH 05/19] neutral_mixed: Change parallel to Upwind Seems to help with smoothness of profiles --- include/div_ops.hxx | 6 +++--- src/neutral_mixed.cxx | 9 +++++---- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/include/div_ops.hxx b/include/div_ops.hxx index fb372a54..09bcccee 100644 --- a/include/div_ops.hxx +++ b/include/div_ops.hxx @@ -525,7 +525,7 @@ const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Fi } } } - + // Y and Z fluxes require Y derivatives // Fields containing values along the magnetic field @@ -543,8 +543,8 @@ const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Fi Field3D yzresult(mesh); yzresult.allocate(); - if (f.hasParallelSlices() && a.hasParallelSlices()) { - // Both inputs have yup and ydown + if (f.hasParallelSlices() && a.hasParallelSlices() && g.hasParallelSlices()) { + // All inputs have yup and ydown fup = f.yup(); fdown = f.ydown(); diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 16a0fd12..305d0343 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -14,6 +14,7 @@ using bout::globals::mesh; /// The limiter method in the radial pressure-diffusion. /// Upwind is consistent with the Y (poloidal) advection. using PerpLimiter = FV::Upwind; +using ParLimiter = FV::Upwind; NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* solver) : name(name) { @@ -332,7 +333,7 @@ void NeutralMixed::finally(const Options& state) { // Neutral density TRACE("Neutral density"); ddt(Nn) = - -FV::Div_par_mod(Nn, Vn, sound_speed) // Parallel advection + - FV::Div_par_mod(Nn, Vn, sound_speed) // Parallel advection + FV::Div_a_Grad_perp_limit(Dnn, Nn, logPnlim) // Perpendicular advection ; @@ -349,7 +350,7 @@ void NeutralMixed::finally(const Options& state) { TRACE("Neutral momentum"); ddt(NVn) = - -AA * FV::Div_par_fvv(Nnlim, Vn, sound_speed) // Momentum flow + -AA * FV::Div_par_fvv(Nnlim, Vn, sound_speed) // Momentum flow - Grad_par(Pn) // Pressure gradient + FV::Div_a_Grad_perp_limit(Dnn, NVn, logPnlim) // Perpendicular advection ; @@ -385,8 +386,8 @@ void NeutralMixed::finally(const Options& state) { // Neutral pressure TRACE("Neutral pressure"); - ddt(Pn) = -FV::Div_par_mod(Pn, Vn, sound_speed) // Parallel advection - - (2. / 3) * Pn * Div_par(Vn) // Compression + ddt(Pn) = - FV::Div_par_mod(Pn, Vn, sound_speed) // Parallel advection + - (2. / 3) * Pn * Div_par(Vn) // Compression + FV::Div_a_Grad_perp_limit(Dnn, Pn, logPnlim) // Perpendicular advection ; From 48a8db9cc4c694804baacc3ba5d2c16d88cf21b2 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 9 Apr 2024 21:56:00 -0700 Subject: [PATCH 06/19] neutral_mixed: Use Div_a_Grad_perp_upwind_flows Same as anomalous_diffusion, obtain flows in X and Y for output when diagnose = true. Note: Only saves perpendicular flux, not parallel flux. --- include/neutral_mixed.hxx | 7 +++- src/neutral_mixed.cxx | 88 ++++++++++++++++++++++++++++++++++++--- 2 files changed, 88 insertions(+), 7 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index e15d3a9a..442c8010 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -42,7 +42,7 @@ private: BoutReal AA; ///< Atomic mass (proton = 1) Field3D Dnn; ///< Diffusion coefficient - Field3D DnnNn; + Field3D DnnNn, DnnPn, DnnNVn; bool sheath_ydown, sheath_yup; @@ -63,6 +63,11 @@ private: bool output_ddt; ///< Save time derivatives? bool diagnose; ///< Save additional diagnostics? + + // Flow diagnostics + Field3D particle_flow_xlow, particle_flow_ylow; + Field3D momentum_flow_xlow, momentum_flow_ylow; + Field3D energy_flow_xlow, energy_flow_ylow; }; namespace { diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 305d0343..f2508db7 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -11,9 +11,6 @@ using bout::globals::mesh; -/// The limiter method in the radial pressure-diffusion. -/// Upwind is consistent with the Y (poloidal) advection. -using PerpLimiter = FV::Upwind; using ParLimiter = FV::Upwind; NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* solver) @@ -152,6 +149,8 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* // Product of Dnn and another parameter has same BC as Dnn - see eqns to see why this is // necessary DnnNn.setBoundary(std::string("Dnn") + name); + DnnPn.setBoundary(std::string("Dnn") + name); + DnnNVn.setBoundary(std::string("Dnn") + name); } void NeutralMixed::transform(Options& state) { @@ -305,14 +304,20 @@ void NeutralMixed::finally(const Options& state) { Dnn.applyBoundary(); // Neutral diffusion parameters have the same boundary condition as Dnn + DnnPn = Dnn * Pn; + DnnPn.applyBoundary(); DnnNn = Dnn * Nn; DnnNn.applyBoundary(); + DnnNVn = Dnn * NVn; + DnnNVn.applyBoundary(); if (sheath_ydown) { for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { Dnn(r.ind, mesh->ystart - 1, jz) = -Dnn(r.ind, mesh->ystart, jz); DnnNn(r.ind, mesh->ystart - 1, jz) = -DnnNn(r.ind, mesh->ystart, jz); + DnnPn(r.ind, mesh->ystart - 1, jz) = -DnnPn(r.ind, mesh->ystart, jz); + DnnNVn(r.ind, mesh->ystart - 1, jz) = -DnnNVn(r.ind, mesh->ystart, jz); } } } @@ -322,6 +327,8 @@ void NeutralMixed::finally(const Options& state) { for (int jz = 0; jz < mesh->LocalNz; jz++) { Dnn(r.ind, mesh->yend + 1, jz) = -Dnn(r.ind, mesh->yend, jz); DnnNn(r.ind, mesh->yend + 1, jz) = -DnnNn(r.ind, mesh->yend, jz); + DnnPn(r.ind, mesh->yend + 1, jz) = -DnnPn(r.ind, mesh->yend, jz); + DnnNVn(r.ind, mesh->yend + 1, jz) = -DnnNVn(r.ind, mesh->yend, jz); } } } @@ -334,7 +341,9 @@ void NeutralMixed::finally(const Options& state) { TRACE("Neutral density"); ddt(Nn) = - FV::Div_par_mod(Nn, Vn, sound_speed) // Parallel advection - + FV::Div_a_Grad_perp_limit(Dnn, Nn, logPnlim) // Perpendicular advection + + Div_a_Grad_perp_upwind_flows(DnnNn, logPnlim, + particle_flow_xlow, + particle_flow_ylow) // Perpendicular advection ; Sn = density_source; // Save for possible output @@ -352,7 +361,9 @@ void NeutralMixed::finally(const Options& state) { ddt(NVn) = -AA * FV::Div_par_fvv(Nnlim, Vn, sound_speed) // Momentum flow - Grad_par(Pn) // Pressure gradient - + FV::Div_a_Grad_perp_limit(Dnn, NVn, logPnlim) // Perpendicular advection + + Div_a_Grad_perp_upwind_flows(DnnNVn, logPnlim, + momentum_flow_xlow, + momentum_flow_ylow) // Perpendicular advection ; if (neutral_viscosity) { @@ -388,8 +399,11 @@ void NeutralMixed::finally(const Options& state) { ddt(Pn) = - FV::Div_par_mod(Pn, Vn, sound_speed) // Parallel advection - (2. / 3) * Pn * Div_par(Vn) // Compression - + FV::Div_a_Grad_perp_limit(Dnn, Pn, logPnlim) // Perpendicular advection + + Div_a_Grad_perp_upwind_flows(DnnPn, logPnlim, + energy_flow_xlow, energy_flow_ylow) // Perpendicular advection ; + energy_flow_xlow *= 3/2; // Note: Should this be 5/2? + energy_flow_ylow *= 3/2; if (neutral_conduction) { ddt(Pn) += FV::Div_a_Grad_perp(DnnNn, Tn) // Perpendicular conduction @@ -441,6 +455,7 @@ void NeutralMixed::outputVars(Options& state) { auto Tnorm = get(state["Tnorm"]); auto Omega_ci = get(state["Omega_ci"]); auto Cs0 = get(state["Cs0"]); + auto rho_s0 = get(state["rho_s0"]); const BoutReal Pnorm = SI::qe * Tnorm * Nnorm; state[std::string("N") + name].setAttributes({{"time_dimension", "t"}, @@ -539,6 +554,67 @@ void NeutralMixed::outputVars(Options& state) { {"long_name", name + " pressure source"}, {"species", name}, {"source", "neutral_mixed"}}); + + if (particle_flow_xlow.isAllocated()) { + set_with_attrs(state[std::string("ParticleFlow_") + name + std::string("_xlow")], particle_flow_xlow, + {{"time_dimension", "t"}, + {"units", "s^-1"}, + {"conversion", rho_s0 * SQ(rho_s0) * Nnorm * Omega_ci}, + {"standard_name", "particle flow"}, + {"long_name", name + " particle flow in X. Note: May be incomplete."}, + {"species", name}, + {"source", "neutral_mixed"}}); + } + if (particle_flow_ylow.isAllocated()) { + set_with_attrs(state[std::string("ParticleFlow_") + name + std::string("_ylow")], particle_flow_ylow, + {{"time_dimension", "t"}, + {"units", "s^-1"}, + {"conversion", rho_s0 * SQ(rho_s0) * Nnorm * Omega_ci}, + {"standard_name", "particle flow"}, + {"long_name", name + " particle flow in Y. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_density"}}); + } + if (momentum_flow_xlow.isAllocated()) { + set_with_attrs(state[std::string("MomentumFlow_") + name + std::string("_xlow")], momentum_flow_xlow, + {{"time_dimension", "t"}, + {"units", "N"}, + {"conversion", rho_s0 * SQ(rho_s0) * SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"standard_name", "momentum flow"}, + {"long_name", name + " momentum flow in X. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_momentum"}}); + } + if (momentum_flow_ylow.isAllocated()) { + set_with_attrs(state[std::string("MomentumFlow_") + name + std::string("_ylow")], momentum_flow_ylow, + {{"time_dimension", "t"}, + {"units", "N"}, + {"conversion", rho_s0 * SQ(rho_s0) * SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"standard_name", "momentum flow"}, + {"long_name", name + " momentum flow in Y. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_momentum"}}); + } + if (energy_flow_xlow.isAllocated()) { + set_with_attrs(state[std::string("EnergyFlow_") + name + std::string("_xlow")], energy_flow_xlow, + {{"time_dimension", "t"}, + {"units", "W"}, + {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci}, + {"standard_name", "power"}, + {"long_name", name + " power through X cell face. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_pressure"}}); + } + if (energy_flow_ylow.isAllocated()) { + set_with_attrs(state[std::string("EnergyFlow_") + name + std::string("_ylow")], energy_flow_ylow, + {{"time_dimension", "t"}, + {"units", "W"}, + {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci}, + {"standard_name", "power"}, + {"long_name", name + " power through Y cell face. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_pressure"}}); + } } } From 7eb92e7f736baf609b11e1d99509fe8f945713fb Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Fri, 12 Apr 2024 12:45:54 +0100 Subject: [PATCH 07/19] Add flag to suppress lax flux - For debugging purposes --- include/neutral_mixed.hxx | 2 ++ src/neutral_mixed.cxx | 9 ++++++++- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 442c8010..0a17efc1 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -56,10 +56,12 @@ private: bool evolve_momentum; ///< Evolve parallel momentum? bool precondition {true}; ///< Enable preconditioner? + bool lax_flux; ///< Use Lax flux for advection terms std::unique_ptr inv; ///< Laplacian inversion used for preconditioning Field3D density_source, pressure_source; ///< External input source Field3D Sn, Sp, Snv; ///< Particle, pressure and momentum source + Field3D sound_speed; ///< Sound speed for use with Lax flux bool output_ddt; ///< Save time derivatives? bool diagnose; ///< Save additional diagnostics? diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index f2508db7..120353f4 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -64,6 +64,10 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Enable preconditioning in neutral model?") .withDefault(true); + lax_flux = options["lax_flux"] + .doc("Enable stabilising lax flux?") + .withDefault(true); + flux_limit = options["flux_limit"] .doc("Limit diffusive fluxes to fraction of thermal speed. <0 means off.") @@ -334,7 +338,10 @@ void NeutralMixed::finally(const Options& state) { } // Sound speed appearing in Lax flux for advection terms - Field3D sound_speed = sqrt(Tn * (5. / 3) / AA); + sound_speed = 0; + if (lax_flux) { + sound_speed = sqrt(Tn * (5. / 3) / AA); + } ///////////////////////////////////////////////////// // Neutral density From b9d8fb94e33d457673a46644e70768540755c157 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Sat, 13 Apr 2024 15:08:23 +0100 Subject: [PATCH 08/19] neutral_mixed: improvements/diags for flooring - added setting for pressure floor. Added temporary debug flags for using floored Nn and Pn in DnnNn and DnnPn calculation (otherwise they could unbalance and lead to density loss - need more testing). There is a secondary floor which kicks in at zero gradient, this is now linked to the primary floor (but 2 orders of magnitude below). --- include/neutral_mixed.hxx | 4 +++ src/neutral_mixed.cxx | 63 +++++++++++++++++++++++++++++++++------ 2 files changed, 58 insertions(+), 9 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 0a17efc1..8a33a12a 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -47,6 +47,7 @@ private: bool sheath_ydown, sheath_yup; BoutReal nn_floor; ///< Minimum Nn used when dividing NVn by Nn to get Vn. + BoutReal pn_floor; ///< Minimum Pn used when dividing Pn by Nn to get Tn. BoutReal flux_limit; ///< Diffusive flux limit BoutReal diffusion_limit; ///< Maximum diffusion coefficient @@ -62,9 +63,12 @@ private: Field3D density_source, pressure_source; ///< External input source Field3D Sn, Sp, Snv; ///< Particle, pressure and momentum source Field3D sound_speed; ///< Sound speed for use with Lax flux + Field3D perp_nn_adv_src; ///< Source due to perpendicular advection operator + Field3D par_nn_adv_src; ///< Source due to parallel advection operator bool output_ddt; ///< Save time derivatives? bool diagnose; ///< Save additional diagnostics? + bool dnnnnfix, dnnpnfix; // Flow diagnostics Field3D particle_flow_xlow, particle_flow_ylow; diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 120353f4..c76aec36 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -60,6 +60,11 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* "Normalised units.") .withDefault(1e-5); + pn_floor = options["pn_floor"] + .doc("A minimum pressure used when dividing Pn by Nn. " + "Normalised units.") + .withDefault(1e-8); + precondition = options["precondition"] .doc("Enable preconditioning in neutral model?") .withDefault(true); @@ -68,6 +73,14 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Enable stabilising lax flux?") .withDefault(true); + dnnpnfix = options["dnnpnfix"] + .doc("Use DnnPn with Pnlim") + .withDefault(false); + + dnnnnfix = options["dnnnnfix"] + .doc("Use DnnNn with Nnlim") + .withDefault(false); + flux_limit = options["flux_limit"] .doc("Limit diffusive fluxes to fraction of thermal speed. <0 means off.") @@ -180,7 +193,7 @@ void NeutralMixed::transform(Options& state) { Vn.applyBoundary("neumann"); Vnlim.applyBoundary("neumann"); - Pnlim = floor(Pn, 1e-8); + Pnlim = floor(Pn, pn_floor); Pnlim.applyBoundary(); ///////////////////////////////////////////////////// @@ -308,9 +321,20 @@ void NeutralMixed::finally(const Options& state) { Dnn.applyBoundary(); // Neutral diffusion parameters have the same boundary condition as Dnn - DnnPn = Dnn * Pn; + if (dnnnnfix) { + DnnNn = Dnn * Nnlim; + } else { + DnnNn = Dnn * Nn; + } + + if (dnnpnfix) { + DnnPn = Dnn * Pnlim; + } else { + DnnPn = Dnn * Pn; + } + DnnPn.applyBoundary(); - DnnNn = Dnn * Nn; + DnnNn.applyBoundary(); DnnNVn = Dnn * NVn; DnnNVn.applyBoundary(); @@ -346,11 +370,16 @@ void NeutralMixed::finally(const Options& state) { ///////////////////////////////////////////////////// // Neutral density TRACE("Neutral density"); - ddt(Nn) = - - FV::Div_par_mod(Nn, Vn, sound_speed) // Parallel advection - + Div_a_Grad_perp_upwind_flows(DnnNn, logPnlim, + + perp_nn_adv_src = Div_a_Grad_perp_upwind_flows(DnnNn, logPnlim, particle_flow_xlow, - particle_flow_ylow) // Perpendicular advection + particle_flow_ylow); // Perpendicular advection + + par_nn_adv_src = FV::Div_par_mod(Nn, Vn, sound_speed); // Parallel advection + + ddt(Nn) = + - par_nn_adv_src + + perp_nn_adv_src ; Sn = density_source; // Save for possible output @@ -425,10 +454,10 @@ void NeutralMixed::finally(const Options& state) { ddt(Pn) += Sp; BOUT_FOR(i, Pn.getRegion("RGN_ALL")) { - if ((Pn[i] < 1e-9) && (ddt(Pn)[i] < 0.0)) { + if ((Pn[i] < pn_floor * 1e-2) && (ddt(Pn)[i] < 0.0)) { ddt(Pn)[i] = 0.0; } - if ((Nn[i] < 1e-7) && (ddt(Nn)[i] < 0.0)) { + if ((Nn[i] < nn_floor * 1e-2) && (ddt(Nn)[i] < 0.0)) { ddt(Nn)[i] = 0.0; } } @@ -561,6 +590,22 @@ void NeutralMixed::outputVars(Options& state) { {"long_name", name + " pressure source"}, {"species", name}, {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("S") + name + std::string("_perp_adv")], perp_nn_adv_src, + {{"time_dimension", "t"}, + {"units", "m^-3 s^-1"}, + {"conversion", Nnorm * Omega_ci}, + {"standard_name", "density source due to perp advection"}, + {"long_name", name + " number density source due to perp advection"}, + {"species", name}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("S") + name + std::string("_par_adv")], par_nn_adv_src, + {{"time_dimension", "t"}, + {"units", "m^-3 s^-1"}, + {"conversion", Nnorm * Omega_ci}, + {"standard_name", "density source due to par advection"}, + {"long_name", name + " number density source due to par advection"}, + {"species", name}, + {"source", "neutral_mixed"}}); if (particle_flow_xlow.isAllocated()) { set_with_attrs(state[std::string("ParticleFlow_") + name + std::string("_xlow")], particle_flow_xlow, From b50e4fa30bfc812a0fbff9d7d5f4103c0e1f95ac Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Wed, 22 May 2024 11:49:25 +0100 Subject: [PATCH 09/19] Hardcode correct use of Nnlim and Pnlim in perp neutral advection --- include/neutral_mixed.hxx | 1 - src/neutral_mixed.cxx | 24 +++--------------------- 2 files changed, 3 insertions(+), 22 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 8a33a12a..baaff4de 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -68,7 +68,6 @@ private: bool output_ddt; ///< Save time derivatives? bool diagnose; ///< Save additional diagnostics? - bool dnnnnfix, dnnpnfix; // Flow diagnostics Field3D particle_flow_xlow, particle_flow_ylow; diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index c76aec36..7edb45be 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -73,14 +73,6 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Enable stabilising lax flux?") .withDefault(true); - dnnpnfix = options["dnnpnfix"] - .doc("Use DnnPn with Pnlim") - .withDefault(false); - - dnnnnfix = options["dnnnnfix"] - .doc("Use DnnNn with Nnlim") - .withDefault(false); - flux_limit = options["flux_limit"] .doc("Limit diffusive fluxes to fraction of thermal speed. <0 means off.") @@ -321,22 +313,12 @@ void NeutralMixed::finally(const Options& state) { Dnn.applyBoundary(); // Neutral diffusion parameters have the same boundary condition as Dnn - if (dnnnnfix) { - DnnNn = Dnn * Nnlim; - } else { - DnnNn = Dnn * Nn; - } - - if (dnnpnfix) { - DnnPn = Dnn * Pnlim; - } else { - DnnPn = Dnn * Pn; - } + DnnNn = Dnn * Nnlim; + DnnPn = Dnn * Pnlim; + DnnNVn = Dnn * NVn; DnnPn.applyBoundary(); - DnnNn.applyBoundary(); - DnnNVn = Dnn * NVn; DnnNVn.applyBoundary(); if (sheath_ydown) { From b31b18e1e6eb41363abc9bf8f86799b0531e0954 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Wed, 22 May 2024 12:39:38 +0100 Subject: [PATCH 10/19] Calculate neutral Pnfloor from Nnfloor like evolve_density --- src/neutral_mixed.cxx | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 7edb45be..4b6f0471 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -58,13 +58,10 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* nn_floor = options["nn_floor"] .doc("A minimum density used when dividing NVn by Nn. " "Normalised units.") - .withDefault(1e-5); - - pn_floor = options["pn_floor"] - .doc("A minimum pressure used when dividing Pn by Nn. " - "Normalised units.") .withDefault(1e-8); + pn_floor = pressure_floor = nn_floor * (1./get(alloptions["units"]["eV"])); + precondition = options["precondition"] .doc("Enable preconditioning in neutral model?") .withDefault(true); From 967ef3e4f76d4ccdd0e1e830b45d72f9d5460704 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Thu, 23 May 2024 12:57:05 +0100 Subject: [PATCH 11/19] Fix typo --- src/neutral_mixed.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 4b6f0471..0d3958cd 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -60,7 +60,7 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* "Normalised units.") .withDefault(1e-8); - pn_floor = pressure_floor = nn_floor * (1./get(alloptions["units"]["eV"])); + pn_floor = nn_floor * (1./get(alloptions["units"]["eV"])); precondition = options["precondition"] .doc("Enable preconditioning in neutral model?") From b8accfe148813d0539adf9b517d0b015723249ff Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Tue, 28 May 2024 15:04:51 +0100 Subject: [PATCH 12/19] Change energy flow diagnostic factor from 3/2 to 5/2 - Seems to make energy imbalance less imbalanced, but I haven't gone through the maths yet --- src/neutral_mixed.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 0d3958cd..7d9811eb 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -417,8 +417,8 @@ void NeutralMixed::finally(const Options& state) { + Div_a_Grad_perp_upwind_flows(DnnPn, logPnlim, energy_flow_xlow, energy_flow_ylow) // Perpendicular advection ; - energy_flow_xlow *= 3/2; // Note: Should this be 5/2? - energy_flow_ylow *= 3/2; + energy_flow_xlow *= 5/2; // Note: Should this be 5/2? + energy_flow_ylow *= 5/2; if (neutral_conduction) { ddt(Pn) += FV::Div_a_Grad_perp(DnnNn, Tn) // Perpendicular conduction From 0aa087b0c2b21e6445a0f84400b11c813fff573a Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Tue, 28 May 2024 15:21:21 +0100 Subject: [PATCH 13/19] Change conduction to upwind and expose fluxes --- include/neutral_mixed.hxx | 1 + src/neutral_mixed.cxx | 23 ++++++++++++++++++++++- 2 files changed, 23 insertions(+), 1 deletion(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index baaff4de..3747878f 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -73,6 +73,7 @@ private: Field3D particle_flow_xlow, particle_flow_ylow; Field3D momentum_flow_xlow, momentum_flow_ylow; Field3D energy_flow_xlow, energy_flow_ylow; + Field3D conduction_flow_xlow, conduction_flow_ylow; }; namespace { diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 7d9811eb..3b71815b 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -421,7 +421,8 @@ void NeutralMixed::finally(const Options& state) { energy_flow_ylow *= 5/2; if (neutral_conduction) { - ddt(Pn) += FV::Div_a_Grad_perp(DnnNn, Tn) // Perpendicular conduction + ddt(Pn) += Div_a_Grad_perp_upwind_flows(DnnNn, Tn, + conduction_flow_xlow, conduction_flow_ylow) // Perpendicular conduction + FV::Div_par_K_Grad_par(DnnNn, Tn) // Parallel conduction ; } @@ -646,6 +647,26 @@ void NeutralMixed::outputVars(Options& state) { {"species", name}, {"source", "evolve_pressure"}}); } + if (conduction_flow_xlow.isAllocated()) { + set_with_attrs(state[std::string("ConductionFlow_") + name + std::string("_xlow")],conduction_flow_xlow, + {{"time_dimension", "t"}, + {"units", "W"}, + {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci}, + {"standard_name", "power"}, + {"long_name", name + " conducted power through X cell face. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_pressure"}}); + } + if (conduction_flow_ylow.isAllocated()) { + set_with_attrs(state[std::string("ConductionFlow_") + name + std::string("_ylow")], conduction_flow_ylow, + {{"time_dimension", "t"}, + {"units", "W"}, + {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci}, + {"standard_name", "power"}, + {"long_name", name + " conducted power through Y cell face. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_pressure"}}); + } } } From 20271186f853a3ab27db7358ceff365f31722ea3 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Tue, 28 May 2024 15:26:42 +0100 Subject: [PATCH 14/19] Add comments on flow diags and fix conduction one --- src/neutral_mixed.cxx | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 3b71815b..b1d7065c 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -417,7 +417,10 @@ void NeutralMixed::finally(const Options& state) { + Div_a_Grad_perp_upwind_flows(DnnPn, logPnlim, energy_flow_xlow, energy_flow_ylow) // Perpendicular advection ; - energy_flow_xlow *= 5/2; // Note: Should this be 5/2? + + // The factor here is likely 5/2 as we're advecting internal energy and pressure. + // Doing this still leaves a heat imbalance factor of 0.11 in the cells, but better than 0.33 with 3/2. + energy_flow_xlow *= 5/2; energy_flow_ylow *= 5/2; if (neutral_conduction) { @@ -426,6 +429,10 @@ void NeutralMixed::finally(const Options& state) { + FV::Div_par_K_Grad_par(DnnNn, Tn) // Parallel conduction ; } + + // The factor here is likely 3/2 as this is pure energy flow, but needs checking. + energy_flow_xlow *= 3/2; + energy_flow_ylow *= 3/2; Sp = pressure_source; if (localstate.isSet("energy_source")) { From c095e434c5dcae3572d03db8d6e26dfa2917d22f Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Tue, 28 May 2024 15:30:25 +0100 Subject: [PATCH 15/19] Fix typo --- src/neutral_mixed.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index b1d7065c..abf75420 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -431,8 +431,8 @@ void NeutralMixed::finally(const Options& state) { } // The factor here is likely 3/2 as this is pure energy flow, but needs checking. - energy_flow_xlow *= 3/2; - energy_flow_ylow *= 3/2; + conduction_flow_xlow *= 3/2; + conduction_flow_ylow *= 3/2; Sp = pressure_source; if (localstate.isSet("energy_source")) { From 620c9144392238fbc29c25f26245c2a4147a6fa6 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 19 Aug 2024 15:54:06 +0100 Subject: [PATCH 16/19] Revert to Div_a_Grad_perp, add flows Upwinding and limiting with MC/MinMod are causing checkerboarding and/or crashes in the neutral_mixed example and provide seemingly no benefit in production runs. Reverting to original operator but with added flow diagnostics. Left warnings on the other operators, maybe we'll get back to them at some point. --- include/div_ops.hxx | 8 +- src/div_ops.cxx | 231 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 238 insertions(+), 1 deletion(-) diff --git a/include/div_ops.hxx b/include/div_ops.hxx index a04fff93..096044c7 100644 --- a/include/div_ops.hxx +++ b/include/div_ops.hxx @@ -61,8 +61,14 @@ const Field3D D4DZ4_Index(const Field3D& f); const Field2D Laplace_FV(const Field2D& k, const Field2D& f); /// Perpendicular diffusion including X and Y directions +/// Takes Div_a_Grad_perp from BOUT++ and adds flows +const Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, + Field3D& flux_xlow, Field3D& flux_ylow); +/// Same but with upwinding +/// WARNING: Causes checkerboarding in neutral_mixed integrated test const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f); -/// Version of function that returns flows +/// Same but with upwinding and flows +/// WARNING: Causes checkerboarding in neutral_mixed integrated test const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, Field3D& flux_xlow, Field3D& flux_ylow); diff --git a/src/div_ops.cxx b/src/div_ops.cxx index 64ff8143..19e9bb33 100644 --- a/src/div_ops.cxx +++ b/src/div_ops.cxx @@ -792,7 +792,237 @@ const Field2D Laplace_FV(const Field2D &k, const Field2D &f) { return result; } +const Field3D Div_a_Grad_perp_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; + + /* + if(mesh->firstX()) + xs += 1; + */ + /* + if(mesh->lastX()) + xe -= 1; + */ + + 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 + + BoutReal fout = 0.5 * (a(i, j, k) + a(i + 1, j, k)) + * (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)); + + 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); + } + } + } + + // 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 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 Jc = coord->J; + Field3D dyc = coord->dy; + Field3D dzc = coord->dz; + Field3D Bxyc = coord->Bxy; + + // 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); + } + + 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: Maybe communicate in init?"); + } + + g23up = coord->g23.yup(); + g23down = coord->g23.ydown(); + + g_23up = coord->g_23.yup(); + g_23down = coord->g_23.ydown(); + + Jup = coord->J.yup(); + Jdown = coord->J.ydown(); + + dyup = coord->dy.yup(); + dydown = coord->dy.ydown(); + + dzup = coord->dz.yup(); + dzdown = coord->dz.ydown(); + + Bxyup = coord->Bxy.yup(); + Bxydown = coord->Bxy.ydown(); + + } 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); + Jup = Jdown = Jc = toFieldAligned(coord->J); + dyup = dydown = dyc = toFieldAligned(coord->dy); + dzup = dzdown = dzc = toFieldAligned(coord->dz); + Bxyup = Bxydown = Bxyc = toFieldAligned(coord->Bxy); + } + + // Y 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 j and j+1 + int kp = (k + 1) % mesh->LocalNz; + int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; + + BoutReal coef = + 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))); + + // 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)); + + BoutReal fout = + 0.25 * (ac(i, j, k) + aup(i, j + 1, k)) + * (Jc(i, j, k) * g23c(i, j, k) + Jup(i, j + 1, k) * g23up(i, j + 1, k)) + * (dfdz - coef * dfdy); + + yzresult(i, j, k) = fout / (dyc(i, j, k) * Jc(i, j, k)); + + // Calculate flux between j and j-1 + coef = + 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))); + + 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)); + + fout = 0.25 * (ac(i, j, k) + adown(i, j - 1, k)) + * (Jc(i, j, k) * g23c(i, j, k) + Jdown(i, j - 1, k) * g23down(i, j - 1, k)) + * (dfdz - coef * 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) * coord->dz(i, j); + } + } + } + + // 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; + + // Coefficient in front of df/dy term + BoutReal coef = g_23c(i, j, k) + / (dyup(i, j + 1, k) + 2. * dyc(i, j, k) + dydown(i, j - 1, k)) + / SQ(Jc(i, j, k) * Bxyc(i, j, k)); + + BoutReal fout = + 0.25 * (ac(i, j, k) + ac(i, j, kp)) + * (Jc(i, j, k) * coord->g33(i, j, k) + Jc(i, j, kp) * coord->g33(i, j, kp)) + * ( // df/dz + (fc(i, j, kp) - fc(i, j, k)) / dzc(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))); + + 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 + if (f.hasParallelSlices() && a.hasParallelSlices()) { + result += yzresult; + } else { + result += fromFieldAligned(yzresult); + flow_ylow = fromFieldAligned(flow_ylow); + } + + return result; +} + // Div ( a Grad_perp(f) ) -- diffusion +/// WARNING: Causes checkerboarding in neutral_mixed integrated test const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { ASSERT2(a.getLocation() == f.getLocation()); @@ -969,6 +1199,7 @@ const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { /// 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) { From 6eb0a04da3e8cd868a00361084757a19dd5b0439 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 19 Aug 2024 15:54:06 +0100 Subject: [PATCH 17/19] Revert to Div_a_Grad_perp, add flows Upwinding and limiting with MC/MinMod are causing checkerboarding and/or crashes in the neutral_mixed example and provide seemingly no benefit in production runs. Reverting to original operator but with added flow diagnostics. Left warnings on the other operators, maybe we'll get back to them at some point. Squashed with commit to update operator choice. --- include/div_ops.hxx | 8 +- src/div_ops.cxx | 231 ++++++++++++++++++++++++++++++++++++++++++ src/neutral_mixed.cxx | 11 +- 3 files changed, 243 insertions(+), 7 deletions(-) diff --git a/include/div_ops.hxx b/include/div_ops.hxx index a04fff93..096044c7 100644 --- a/include/div_ops.hxx +++ b/include/div_ops.hxx @@ -61,8 +61,14 @@ const Field3D D4DZ4_Index(const Field3D& f); const Field2D Laplace_FV(const Field2D& k, const Field2D& f); /// Perpendicular diffusion including X and Y directions +/// Takes Div_a_Grad_perp from BOUT++ and adds flows +const Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, + Field3D& flux_xlow, Field3D& flux_ylow); +/// Same but with upwinding +/// WARNING: Causes checkerboarding in neutral_mixed integrated test const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f); -/// Version of function that returns flows +/// Same but with upwinding and flows +/// WARNING: Causes checkerboarding in neutral_mixed integrated test const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, Field3D& flux_xlow, Field3D& flux_ylow); diff --git a/src/div_ops.cxx b/src/div_ops.cxx index 64ff8143..19e9bb33 100644 --- a/src/div_ops.cxx +++ b/src/div_ops.cxx @@ -792,7 +792,237 @@ const Field2D Laplace_FV(const Field2D &k, const Field2D &f) { return result; } +const Field3D Div_a_Grad_perp_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; + + /* + if(mesh->firstX()) + xs += 1; + */ + /* + if(mesh->lastX()) + xe -= 1; + */ + + 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 + + BoutReal fout = 0.5 * (a(i, j, k) + a(i + 1, j, k)) + * (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)); + + 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); + } + } + } + + // 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 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 Jc = coord->J; + Field3D dyc = coord->dy; + Field3D dzc = coord->dz; + Field3D Bxyc = coord->Bxy; + + // 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); + } + + 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: Maybe communicate in init?"); + } + + g23up = coord->g23.yup(); + g23down = coord->g23.ydown(); + + g_23up = coord->g_23.yup(); + g_23down = coord->g_23.ydown(); + + Jup = coord->J.yup(); + Jdown = coord->J.ydown(); + + dyup = coord->dy.yup(); + dydown = coord->dy.ydown(); + + dzup = coord->dz.yup(); + dzdown = coord->dz.ydown(); + + Bxyup = coord->Bxy.yup(); + Bxydown = coord->Bxy.ydown(); + + } 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); + Jup = Jdown = Jc = toFieldAligned(coord->J); + dyup = dydown = dyc = toFieldAligned(coord->dy); + dzup = dzdown = dzc = toFieldAligned(coord->dz); + Bxyup = Bxydown = Bxyc = toFieldAligned(coord->Bxy); + } + + // Y 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 j and j+1 + int kp = (k + 1) % mesh->LocalNz; + int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; + + BoutReal coef = + 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))); + + // 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)); + + BoutReal fout = + 0.25 * (ac(i, j, k) + aup(i, j + 1, k)) + * (Jc(i, j, k) * g23c(i, j, k) + Jup(i, j + 1, k) * g23up(i, j + 1, k)) + * (dfdz - coef * dfdy); + + yzresult(i, j, k) = fout / (dyc(i, j, k) * Jc(i, j, k)); + + // Calculate flux between j and j-1 + coef = + 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))); + + 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)); + + fout = 0.25 * (ac(i, j, k) + adown(i, j - 1, k)) + * (Jc(i, j, k) * g23c(i, j, k) + Jdown(i, j - 1, k) * g23down(i, j - 1, k)) + * (dfdz - coef * 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) * coord->dz(i, j); + } + } + } + + // 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; + + // Coefficient in front of df/dy term + BoutReal coef = g_23c(i, j, k) + / (dyup(i, j + 1, k) + 2. * dyc(i, j, k) + dydown(i, j - 1, k)) + / SQ(Jc(i, j, k) * Bxyc(i, j, k)); + + BoutReal fout = + 0.25 * (ac(i, j, k) + ac(i, j, kp)) + * (Jc(i, j, k) * coord->g33(i, j, k) + Jc(i, j, kp) * coord->g33(i, j, kp)) + * ( // df/dz + (fc(i, j, kp) - fc(i, j, k)) / dzc(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))); + + 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 + if (f.hasParallelSlices() && a.hasParallelSlices()) { + result += yzresult; + } else { + result += fromFieldAligned(yzresult); + flow_ylow = fromFieldAligned(flow_ylow); + } + + return result; +} + // Div ( a Grad_perp(f) ) -- diffusion +/// WARNING: Causes checkerboarding in neutral_mixed integrated test const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { ASSERT2(a.getLocation() == f.getLocation()); @@ -969,6 +1199,7 @@ const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { /// 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) { diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index abf75420..4ce16fd8 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -350,7 +350,7 @@ void NeutralMixed::finally(const Options& state) { // Neutral density TRACE("Neutral density"); - perp_nn_adv_src = Div_a_Grad_perp_upwind_flows(DnnNn, logPnlim, + perp_nn_adv_src = Div_a_Grad_perp_flows(DnnNn, logPnlim, particle_flow_xlow, particle_flow_ylow); // Perpendicular advection @@ -376,7 +376,7 @@ void NeutralMixed::finally(const Options& state) { ddt(NVn) = -AA * FV::Div_par_fvv(Nnlim, Vn, sound_speed) // Momentum flow - Grad_par(Pn) // Pressure gradient - + Div_a_Grad_perp_upwind_flows(DnnNVn, logPnlim, + + Div_a_Grad_perp_flows(DnnNVn, logPnlim, momentum_flow_xlow, momentum_flow_ylow) // Perpendicular advection ; @@ -414,17 +414,16 @@ void NeutralMixed::finally(const Options& state) { ddt(Pn) = - FV::Div_par_mod(Pn, Vn, sound_speed) // Parallel advection - (2. / 3) * Pn * Div_par(Vn) // Compression - + Div_a_Grad_perp_upwind_flows(DnnPn, logPnlim, + + Div_a_Grad_perp_flows(DnnPn, logPnlim, energy_flow_xlow, energy_flow_ylow) // Perpendicular advection ; - // The factor here is likely 5/2 as we're advecting internal energy and pressure. - // Doing this still leaves a heat imbalance factor of 0.11 in the cells, but better than 0.33 with 3/2. + // The factor here is 5/2 as we're advecting internal energy and pressure. energy_flow_xlow *= 5/2; energy_flow_ylow *= 5/2; if (neutral_conduction) { - ddt(Pn) += Div_a_Grad_perp_upwind_flows(DnnNn, Tn, + ddt(Pn) += Div_a_Grad_perp_flows(DnnNn, Tn, conduction_flow_xlow, conduction_flow_ylow) // Perpendicular conduction + FV::Div_par_K_Grad_par(DnnNn, Tn) // Parallel conduction ; From e2e19f3367fd244a5b7e4668b0deeb11e0918ddc Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Fri, 23 Aug 2024 13:01:01 +0100 Subject: [PATCH 18/19] Add flags to suppress par or perp neutral transport --- include/neutral_mixed.hxx | 2 ++ src/neutral_mixed.cxx | 69 +++++++++++++++++++++++++-------------- 2 files changed, 47 insertions(+), 24 deletions(-) diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 3747878f..745feb22 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -56,6 +56,8 @@ private: bool neutral_conduction; ///< Include heat conduction? bool evolve_momentum; ///< Evolve parallel momentum? + bool parallel_transport, perpendicular_transport; ///< Include transport in a particular direction? + bool precondition {true}; ///< Enable preconditioner? bool lax_flux; ///< Use Lax flux for advection terms std::unique_ptr inv; ///< Laplacian inversion used for preconditioning diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 4ce16fd8..24fcc2e4 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -88,6 +88,14 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Include neutral gas heat conduction?") .withDefault(true); + perpendicular_transport = options["perpendicular_transport"] + .doc("Solve perpendicular transport?") + .withDefault(true); + + parallel_transport = options["parallel_transport"] + .doc("Solve parallel transport?") + .withDefault(true); + if (precondition) { inv = std::unique_ptr(Laplacian::create(&options["precon_laplace"])); @@ -356,10 +364,12 @@ void NeutralMixed::finally(const Options& state) { par_nn_adv_src = FV::Div_par_mod(Nn, Vn, sound_speed); // Parallel advection - ddt(Nn) = - - par_nn_adv_src - + perp_nn_adv_src - ; + if (parallel_transport) { + ddt(Nn) -= par_nn_adv_src; + } + if (perpendicular_transport) { + ddt(Nn) += perp_nn_adv_src; + } Sn = density_source; // Save for possible output if (localstate.isSet("density_source")) { @@ -373,13 +383,16 @@ void NeutralMixed::finally(const Options& state) { // Neutral momentum TRACE("Neutral momentum"); - ddt(NVn) = - -AA * FV::Div_par_fvv(Nnlim, Vn, sound_speed) // Momentum flow - - Grad_par(Pn) // Pressure gradient - + Div_a_Grad_perp_flows(DnnNVn, logPnlim, + if (parallel_transport) { + ddt(NVn) -= AA * FV::Div_par_fvv(Nnlim, Vn, sound_speed); // Momentum flow + ddt(NVn) -= Grad_par(Pn) // Pressure gradient + } + if (perpendicular_transport) { + ddt(NVn) += Div_a_Grad_perp_flows(DnnNVn, logPnlim, momentum_flow_xlow, - momentum_flow_ylow) // Perpendicular advection - ; + momentum_flow_ylow); // Perpendicular advection + } + if (neutral_viscosity) { // NOTE: The following viscosity terms are not (yet) balanced @@ -392,11 +405,12 @@ void NeutralMixed::finally(const Options& state) { // eta_n = (2. / 5) * kappa_n; // - ddt(NVn) += - AA * FV::Div_a_Grad_perp((2. / 5) * DnnNn, Vn) // Perpendicular viscosity - + AA * FV::Div_par_K_Grad_par((2. / 5) * DnnNn, Vn) // Parallel viscosity - ; - } + if (parallel_transport) { + ddt(NVn) += AA * FV::Div_par_K_Grad_par((2. / 5) * DnnNn, Vn); // Parallel viscosity + } + if (perpendicular_transport) { + ddt(NVn) += AA * FV::Div_a_Grad_perp((2. / 5) * DnnNn, Vn); // Perpendicular viscosity + } if (localstate.isSet("momentum_source")) { Snv = get(localstate["momentum_source"]); @@ -412,21 +426,28 @@ void NeutralMixed::finally(const Options& state) { // Neutral pressure TRACE("Neutral pressure"); - ddt(Pn) = - FV::Div_par_mod(Pn, Vn, sound_speed) // Parallel advection - - (2. / 3) * Pn * Div_par(Vn) // Compression - + Div_a_Grad_perp_flows(DnnPn, logPnlim, - energy_flow_xlow, energy_flow_ylow) // Perpendicular advection - ; + if (parallel_transport) { + ddt(NVn) -= FV::Div_par_mod(Pn, Vn, sound_speed); // Parallel advection + ddt(NVn) -= (2. / 3) * Pn * Div_par(Vn) ; // Parallel advection + } + if (perpendicular_transport) { + ddt(NVn) += Div_a_Grad_perp_flows(DnnPn, logPnlim, + energy_flow_xlow, energy_flow_ylow); // Perpendicular viscosity + } + // The factor here is 5/2 as we're advecting internal energy and pressure. energy_flow_xlow *= 5/2; energy_flow_ylow *= 5/2; if (neutral_conduction) { - ddt(Pn) += Div_a_Grad_perp_flows(DnnNn, Tn, - conduction_flow_xlow, conduction_flow_ylow) // Perpendicular conduction - + FV::Div_par_K_Grad_par(DnnNn, Tn) // Parallel conduction - ; + if (parallel_transport) { + ddt(Pn) += FV::Div_par_K_Grad_par(DnnNn, Tn); // Parallel conduction + } + if (perpendicular_transport) { + ddt(Pn) += Div_a_Grad_perp_flows(DnnNn, Tn, + conduction_flow_xlow, conduction_flow_ylow); // Perpendicular conduction + } } // The factor here is likely 3/2 as this is pure energy flow, but needs checking. From d27b5f50f0eca1281cf32942616a48867c992c31 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Fri, 23 Aug 2024 13:11:35 +0100 Subject: [PATCH 19/19] Fix silly mistakes --- src/neutral_mixed.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index 24fcc2e4..226cdf0a 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -385,7 +385,7 @@ void NeutralMixed::finally(const Options& state) { if (parallel_transport) { ddt(NVn) -= AA * FV::Div_par_fvv(Nnlim, Vn, sound_speed); // Momentum flow - ddt(NVn) -= Grad_par(Pn) // Pressure gradient + ddt(NVn) -= Grad_par(Pn); // Pressure gradient } if (perpendicular_transport) { ddt(NVn) += Div_a_Grad_perp_flows(DnnNVn, logPnlim, @@ -411,6 +411,7 @@ void NeutralMixed::finally(const Options& state) { if (perpendicular_transport) { ddt(NVn) += AA * FV::Div_a_Grad_perp((2. / 5) * DnnNn, Vn); // Perpendicular viscosity } + } if (localstate.isSet("momentum_source")) { Snv = get(localstate["momentum_source"]);