diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 94c599c58c..2ae0d3d547 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -5,17 +5,11 @@ #ifndef BOUT_FV_OPS_H #define BOUT_FV_OPS_H -#include "bout/assert.hxx" #include "bout/bout_types.hxx" -#include "bout/build_defines.hxx" -#include "bout/coordinates.hxx" #include "bout/field.hxx" #include "bout/field3d.hxx" -#include "bout/globals.hxx" #include "bout/mesh.hxx" #include "bout/output_bout_types.hxx" // NOLINT(unused-includes, misc-include-cleaner) -#include "bout/region.hxx" -#include "bout/utils.hxx" #include "bout/vector2d.hxx" namespace FV { @@ -166,5 +160,22 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, bool fixflux = true); + +/// Calculates viscous heating due to numerical momentum fluxes +/// and flow of kinetic energy (in flow_ylow) +template +Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, Field3D& flow_ylow, + bool fixflux = true); + +/// 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, +/// the advects the upwind cell edge. +/// +/// 1st order upwinding is used in Y. +template +Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D& f); + } // namespace FV #endif // BOUT_FV_OPS_H diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx index 7da980e584..1b7a36f874 100644 --- a/include/bout/fv_ops_impl.hxx +++ b/include/bout/fv_ops_impl.hxx @@ -91,7 +91,7 @@ private: return 0.0; } - if (fabs(a) < fabs(b)) { + if (std::abs(a) < std::abs(b)) { return a; } return b; @@ -126,7 +126,7 @@ private: } // Return the minimum absolute value - return SIGN(a) * BOUTMIN(fabs(a), fabs(b), fabs(c)); + return SIGN(a) * BOUTMIN(std::abs(a), std::abs(b), std::abs(c)); } }; @@ -167,8 +167,8 @@ struct Superbee { n.L = n.R = n.c; } else { const BoutReal sign = SIGN(gL); - const BoutReal abs_gL = fabs(gL); - const BoutReal abs_gR = fabs(gR); + const BoutReal abs_gL = std::abs(gL); + const BoutReal abs_gR = std::abs(gR); const BoutReal half_slope = sign * BOUTMAX(BOUTMIN(abs_gL, 0.5 * abs_gR), BOUTMIN(abs_gR, 0.5 * abs_gL)); n.L = n.c - half_slope; @@ -270,11 +270,21 @@ struct WENO3 { } }; -/*! - * Communicate fluxes between processors - * Takes values in guard cells, and adds them to cells - */ -void communicateFluxes(Field3D& f); +template +struct Slices { + T c; + T up; + T down; + + Slices(bool use_slices, const T& field) + : c(use_slices ? field : toFieldAligned(field)), up(use_slices ? field.yup() : c), + down(use_slices ? field.ydown() : c) {} +}; + +template +Slices makeslices(bool use_slices, const T& field) { + return Slices(use_slices, field); +} /// Finite volume parallel divergence /// @@ -373,7 +383,7 @@ Field3D Div_par(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_sp // Calculate velocity at right boundary (y+1/2) BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j + 1, k)); - BoutReal flux = NAN; + BoutReal flux = BoutNaN; if (is_last_y && (j == mesh->yend) && !is_periodic_y) { // Last point in domain @@ -499,7 +509,7 @@ Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { if ((i.x() == mesh->xend) && (mesh->lastX())) { // At right boundary in X if (bndry_flux) { - BoutReal flux = NAN; + BoutReal flux = BoutNaN; if (vR > 0.0) { // Flux to boundary flux = vR * s.R; @@ -526,7 +536,7 @@ Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { // At left boundary in X if (bndry_flux) { - BoutReal flux = NAN; + BoutReal flux = BoutNaN; if (vL < 0.0) { // Flux to boundary flux = vL * s.L; @@ -748,7 +758,7 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, } else { // Maximum wave speed in the two cells const 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))); + std::abs(v(i, j, k)), std::abs(v(i, j + 1, k))); flux = s.R * 0.5 * (sv.R + amax); } @@ -776,7 +786,7 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, // Maximum wave speed in the two cells const 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))); + std::abs(v(i, j, k)), std::abs(v(i, j - 1, k))); flux = s.L * 0.5 * (sv.L - amax); } @@ -832,8 +842,8 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, 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])); + const BoutReal amax = BOUTMAX(wave_speed_in[i], std::abs(v_in[i]), + std::abs(v_up[iyp]), std::abs(v_down[iym])); const BoutReal term = (f_up[iyp] * v_up[iyp] * v_up[iyp] / B_up[iyp]) - (f_down[iym] * v_down[iym] * v_down[iym] / B_down[iym]); @@ -844,7 +854,7 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, + (amax * (f_in[i] * v_in[i] - f_down[iym] * v_down[iym]) / (B[i] + B_down[iym])); - if (fabs(penalty) > fabs(term) and penalty * v_in[i] > 0) { + if (std::abs(penalty) > std::abs(term) and penalty * v_in[i] > 0) { if (term * penalty > 0) { penalty = term; } else { @@ -925,7 +935,7 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, const BoutReal v_mid_r = 0.5 * (sv.c + sv.p); // And mid-point density at right boundary const BoutReal n_mid_r = 0.5 * (s.c + s.p); - BoutReal flux = NAN; + BoutReal flux = BoutNaN; if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) { // Last point in domain @@ -935,14 +945,15 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, flux = n_mid_r * v_mid_r * v_mid_r; } 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_r - * (sv.R - v_mid_r)); // Damp differences in velocity, not flux + flux = + (s.R * sv.R * sv.R) // Use right cell edge values + + (BOUTMAX(wave_speed(i, j, k), std::abs(sv.c), std::abs(sv.p)) * n_mid_r + * (sv.R - v_mid_r)); // 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)); + std::abs(sv.c), std::abs(sv.p)); flux = s.R * 0.5 * (sv.R + amax) * sv.R; } @@ -964,13 +975,13 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, } 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_l - * (sv.L - v_mid_l)); + - (BOUTMAX(wave_speed(i, j, k), std::abs(sv.c), std::abs(sv.m)) + * n_mid_l * (sv.L - v_mid_l)); } } 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)); + std::abs(sv.c), std::abs(sv.m)); flux = s.L * 0.5 * (sv.L - amax) * sv.L; } @@ -982,5 +993,521 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, } return fromFieldAligned(result, "RGN_NOBNDRY"); } + +// Calculates viscous heating due to numerical momentum fluxes +// and flow of kinetic energy (in flow_ylow) +template +Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, Field3D& flow_ylow, + bool fixflux) { + + ASSERT1(areFieldsCompatible(f_in, v_in)); + ASSERT1(areFieldsCompatible(f_in, wave_speed_in)); + + const 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)}; + flow_ylow = zeroFrom(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], std::abs(v_in[i]), + std::abs(v_up[iyp]), std::abs(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; + } + + /// 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"); + + // result and flow_ylow are field-aligned. + // Will be converted to non-aligned before return. + 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++) { + for (int k = 0; k < mesh->LocalNz; k++) { + // Pre-calculate factors which multiply fluxes + // Note: In 3D metric geometries these quantities can depend on (i,j,k), + // so calculate inside the k loop. + + // 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 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))); + + const BoutReal flux_factor_lc = + common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); + const BoutReal area_lc = common_factor * coord->dx(i, j, k) * coord->dz(i, j, 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); + + if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) { + // Last point in domain + + // Expected loss of kinetic energy into boundary + // This is used in the sheath boundary condition to calculate + // energy losses. + const BoutReal expected_ke = 0.5 * n_mid * v_mid * v_mid * v_mid; + + BoutReal flux_mom = BoutNaN; + if (fixflux) { + // Mid-point consistent with boundary conditions + // but kinetic energy loss will not match expected + // -> Adjust energy balance in pressure equation + flux_mom = n_mid * v_mid * v_mid; + } else { + flux_mom = (s.R * sv.R * sv.R) + + (BOUTMAX(wave_speed(i, j, k), std::abs(sv.c), std::abs(sv.p)) + * (s.R * sv.R - n_mid * v_mid)); + } + + // Assume that particle flux is fixed to boundary value + const BoutReal flux_part = n_mid * v_mid; + + // d/dt(1/2 m n v^2) = v * d/dt(mnv) - 1/2 m v^2 * dn/dt + const BoutReal actual_ke = (sv.c * flux_mom) - (0.5 * sv.c * sv.c * flux_part); + + // Note: If the actual loss was higher than expected, then + // plasma heating is needed to compensate + result(i, j, k) += (actual_ke - expected_ke) * flux_factor_rc; + + // Final flow through boundary is the expected value + flow_ylow(i, j + 1, k) += expected_ke * area_rp; //expected_ke * area_rp; + + } else { + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k), + std::abs(sv.c), std::abs(sv.p)); + + // Viscous heating due to relaxation of velocity towards midpoint + result(i, j, k) += + (amax + 0.5 * sv.R) * s.R * (sv.c - sv.p) * (sv.R - v_mid) * flux_factor_rc; + + // Kinetic energy flow into next cell. + // Note: Different from flow out of this cell; the difference + // is in the viscous heating. + const BoutReal flux_part = s.R * 0.5 * (sv.R + amax); + const BoutReal flux_mom = flux_part * sv.R; + + flow_ylow(i, j + 1, k) += + (sv.p * flux_mom - 0.5 * SQ(sv.p) * flux_part) * area_rp; + } + + //////////////////////////////////////////// + // Calculate at left boundary + + v_mid = 0.5 * (sv.c + sv.m); + n_mid = 0.5 * (s.c + s.m); + + // Expected KE loss. Note minus sign because negative v into boundary + const BoutReal expected_ke = -0.5 * n_mid * v_mid * v_mid * v_mid; + + if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) { + // First point in domain + BoutReal flux_mom = BoutNaN; + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux_mom = n_mid * v_mid * v_mid; + } else { + // Add flux due to difference in boundary values + flux_mom = (s.L * sv.L * sv.L) + - (BOUTMAX(wave_speed(i, j, k), std::abs(sv.c), std::abs(sv.m)) + * (s.L * sv.L - n_mid * v_mid)); + } + + // Assume that density flux is fixed to boundary value + const BoutReal flux_part = n_mid * v_mid; + + // d/dt(1/2 m n v^2) = v * d/dt(mnv) - 1/2 m v^2 * dn/dt + const BoutReal actual_ke = (-sv.c * flux_mom) + (0.5 * sv.c * sv.c * flux_part); + + result(i, j, k) += (actual_ke - expected_ke) * flux_factor_lc; + + flow_ylow(i, j, k) -= expected_ke * area_lc; + } else { + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k), + std::abs(sv.c), std::abs(sv.m)); + + // Viscous heating due to relaxation + result(i, j, k) += + (amax - 0.5 * sv.L) * s.L * (sv.c - sv.m) * (sv.L - v_mid) * flux_factor_lc; + + // Kinetic energy flow into this cell. + // Note: Different from flow out of left cell; the difference + // is in the viscous heating. + const BoutReal flux_part = s.L * 0.5 * (sv.L - amax); + const BoutReal flux_mom = flux_part * sv.L; + + flow_ylow(i, j, k) += (sv.c * flux_mom - 0.5 * SQ(sv.c) * flux_part) * area_lc; + } + } + } + } + flow_ylow = fromFieldAligned(flow_ylow, "RGN_NOBNDRY"); + return fromFieldAligned(result, "RGN_NOBNDRY"); +} + +/// 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, +/// the advects the upwind cell edge. +/// +/// 1st order upwinding is used in Y. +template +Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, 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 + + 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); + + // Mid-point average boundary value of 'a' + const BoutReal aedge = 0.5 * (a(i + 1, j, k) + a(i, j, k)); + BoutReal gedge = BoutNaN; + 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 `g` at left of (i+1, j, k) + + 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 + + gedge = sg.L; + } else { + // Flux is from (i) to (i+1) + // Reconstruct `g` at right of (i, j, k) + + 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 + + gedge = sg.R; + } + + // Flux across cell edge + const BoutReal fout = gradient * aedge * gedge + * (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, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= fout / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); + } + } + } + + const bool fci = + f.hasParallelSlices() && a.hasParallelSlices() && g.hasParallelSlices(); + +#if BOUT_USE_METRIC_3D + if (fci) { + // 3D Metric, need yup/ydown fields. + // Requires previous communication of metrics. + 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?"); + } + } +#endif + + // Y and Z fluxes require Y derivatives + + // Values on this y slice (centre). + // This is needed because toFieldAligned may modify the field + const auto f_slice = makeslices(fci, f); + const auto a_slice = makeslices(fci, a); + const auto g_slice = makeslices(fci, g); + +#if BOUT_USE_METRIC_3D + const bool metric_fci = fci; +#else + constexpr bool metric_fci = false; +#endif + const auto g23 = makeslices(metric_fci, coord->g23); + const auto g_23 = makeslices(metric_fci, coord->g_23); + const auto J = makeslices(metric_fci, coord->J); + const auto dy = makeslices(metric_fci, coord->dy); + const auto dz = makeslices(metric_fci, coord->dz); + const auto Bxy = makeslices(metric_fci, coord->Bxy); + + // Result of the Y and Z fluxes + Field3D yzresult(0.0, mesh); + if (!fci) { + yzresult.setDirectionY(YDirectionType::Aligned); + } + + // Y flux + + 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 + const int k = 0; +#endif + { + const BoutReal coef_u = + 0.5 + * (g_23.c(i, j, k) / SQ(J.c(i, j, k) * Bxy.c(i, j, k)) + + g_23.up(i, j + 1, k) / SQ(J.up(i, j + 1, k) * Bxy.up(i, j + 1, k))); + + const BoutReal coef_d = + 0.5 + * (g_23.c(i, j, k) / SQ(J.c(i, j, k) * Bxy.c(i, j, k)) + + g_23.down(i, j - 1, k) + / SQ(J.down(i, j - 1, k) * Bxy.down(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 + const int kp = (k + 1) % mesh->LocalNz; + const int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; + + // Calculate Z derivative at y boundary + BoutReal dfdz = 0.5 + * (f_slice.c(i, j, kp) - f_slice.c(i, j, km) + + f_slice.up(i, j + 1, kp) - f_slice.up(i, j + 1, km)) + / (dz.c(i, j, k) + dz.up(i, j + 1, k)); + + // Y derivative + BoutReal dfdy = 2. * (f_slice.up(i, j + 1, k) - f_slice.c(i, j, k)) + / (dy.up(i, j + 1, k) + dy.c(i, j, k)); + + BoutReal aedge = 0.5 * (a_slice.c(i, j, k) + a_slice.up(i, j + 1, k)); + BoutReal gedge = BoutNaN; + if ((j == mesh->yend) and mesh->lastY(i)) { + // Midpoint boundary value + gedge = 0.5 * (g_slice.c(i, j, k) + g_slice.up(i, j + 1, k)); + } else if (dfdy > 0) { + // Flux from (j+1) to (j) + gedge = g_slice.up(i, j + 1, k); + } else { + // Flux from (j) to (j+1) + gedge = g_slice.c(i, j, k); + } + + BoutReal fout = + aedge * gedge * 0.5 + * (J.c(i, j, k) * g23.c(i, j, k) + J.up(i, j + 1, k) * g23.up(i, j + 1, k)) + * (dfdz - coef_u * dfdy); + + yzresult(i, j, k) += fout / (dy.c(i, j, k) * J.c(i, j, k)); + + // Calculate flux between j and j-1 + dfdz = 0.5 + * (f_slice.c(i, j, kp) - f_slice.c(i, j, km) + f_slice.down(i, j - 1, kp) + - f_slice.down(i, j - 1, km)) + / (dz.c(i, j, k) + dz.down(i, j - 1, k)); + + dfdy = 2. * (f_slice.c(i, j, k) - f_slice.down(i, j - 1, k)) + / (dy.c(i, j, k) + dy.down(i, j - 1, k)); + + aedge = 0.5 * (a_slice.c(i, j, k) + a_slice.down(i, j - 1, k)); + if ((j == mesh->ystart) and mesh->firstY(i)) { + gedge = 0.5 * (g_slice.c(i, j, k) + g_slice.down(i, j - 1, k)); + } else if (dfdy > 0) { + gedge = g_slice.c(i, j, k); + } else { + gedge = g_slice.down(i, j - 1, k); + } + + fout = aedge * gedge * 0.5 + * (J.c(i, j, k) * g23.c(i, j, k) + + J.down(i, j - 1, k) * g23.down(i, j - 1, k)) + * (dfdz - coef_d * dfdy); + + yzresult(i, j, k) -= fout / (dy.c(i, j, k) * J.c(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++) { +#if BOUT_USE_METRIC_3D + for (int k = 0; k < mesh->LocalNz; k++) +#else + const int k = 0; +#endif + { + // Coefficient in front of df/dy term + const BoutReal coef = + g_23.c(i, j, k) + / (dy.up(i, j + 1, k) + 2. * dy.c(i, j, k) + dy.down(i, j - 1, k)) + / SQ(J.c(i, j, k) * Bxy.c(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 + const int kp = (k + 1) % mesh->LocalNz; + + const BoutReal gradient = + // df/dz + ((f_slice.c(i, j, kp) - f_slice.c(i, j, k)) / dz.c(i, j, k)) + + // - g_yz * df/dy / SQ(J*B) + - (coef + * (f_slice.up(i, j + 1, k) + f_slice.up(i, j + 1, kp) + - f_slice.down(i, j - 1, k) - f_slice.down(i, j - 1, kp))); + + const BoutReal fout = + gradient * 0.5 * (a_slice.c(i, j, kp) + a_slice.c(i, j, k)) + * ((gradient > 0) ? g_slice.c(i, j, kp) : g_slice.c(i, j, k)); + + yzresult(i, j, k) += fout / dz.c(i, j, k); + yzresult(i, j, kp) -= fout / dz.c(i, j, kp); + } + } + } + } + // Check if we need to transform back + if (fci) { + result += yzresult; + } else { + result += fromFieldAligned(yzresult); + } + return result; +} + } // namespace FV #endif // BOUT_FV_OPS_H diff --git a/manual/sphinx/user_docs/coordinates.rst b/manual/sphinx/user_docs/coordinates.rst index 2ac086e05e..d270dda468 100644 --- a/manual/sphinx/user_docs/coordinates.rst +++ b/manual/sphinx/user_docs/coordinates.rst @@ -77,7 +77,7 @@ The cross products are: .. math:: :label: eq:psithetazetacrossproducts - \boldsymbol{e}_\psi\times\boldsymbol{e}_\theta = J_{\psi\theta\zeta} \nabla\zeta \qquad + \boldsymbol{e}_\psi\times\boldsymbol{e}_\theta = J_{\psi\theta\zeta} \nabla\zeta \qquad \boldsymbol{e}_\psi\times\boldsymbol{e}_\zeta = -J_{\psi\theta\zeta} \nabla\theta \qquad \boldsymbol{e}_\theta\times\boldsymbol{e}_\zeta = J_{\psi\theta\zeta} \nabla\psi @@ -122,7 +122,7 @@ To align to the magnetic field we define a local field line pitch `\nu`: \frac{{B_{\text{tor}}}{h_\theta}}{{B_{\text{pol}}}R} \end{aligned} -The sign of the poloidal field `{B_{\text{pol}}}` and toroidal field +The sign of the poloidal field `{B_{\text{pol}}}` and toroidal field `{B_{\text{tor}}}` can be either + or -. The field-aligned coordinates `\left(x,y,z\right)` are defined by: @@ -158,7 +158,7 @@ The reciprocal basis vectors are \nabla z = \nabla\zeta - \sigma_{B\text{pol}}\left[\int_{\theta_0}^\theta{\frac{\partial \nu\left(\psi,\theta\right)}{\partial \psi}} d\theta\right] \nabla\psi - \sigma_{B\text{pol}}\nu\left(\psi, \theta\right)\nabla\theta \end{aligned} - + The term in square brackets is the integrated local shear: .. math:: @@ -167,18 +167,18 @@ The term in square brackets is the integrated local shear: \begin{aligned} I = \int_{y_0}^y\frac{\partial\nu\left(x, y\right)}{\partial\psi}dy\end{aligned} - + The basis vectors are: .. math:: :label: eq:basisvectors - + \begin{aligned} \boldsymbol{e}_x =& J_{xyz}\left(\nabla y \times \nabla z\right) = {\sigma_{B\text{pol}}} {\boldsymbol{e}}_\psi + I{\boldsymbol{e}}_\zeta \\ \boldsymbol{e}_y =& J_{xyz}\left(\nabla z \times \nabla x\right) = {\boldsymbol{e}}_\theta + \nu{\boldsymbol{e}}_\phi \\ \boldsymbol{e}_z =& J_{xyz}\left(\nabla x \times \nabla y\right) = {\boldsymbol{e}}_\zeta \end{aligned} - + where `{\boldsymbol{e}}_\phi = {\sigma_{B\text{pol}}}{\boldsymbol{e}}_\zeta` is always anticlockwise when seen from above the tokamak looking down. The direction of @@ -458,7 +458,7 @@ The reciprocal basis vectors are .. math:: :label: eq:reciprocalbasisvectorsrighthanded - + \begin{aligned} \nabla x =& {\sigma_{B\text{pol}}} \nabla \psi \\ \nabla \eta =& {\sigma_{B\text{pol}}} \nabla \theta \\ @@ -469,7 +469,7 @@ and basis vectors .. math:: :label: eq:basisvectorsrighthanded - + \begin{aligned} \boldsymbol{e}_x =& J_{x\eta z}\left(\nabla y \times \nabla z\right) = {\sigma_{B\text{pol}}} {\boldsymbol{e}}_\psi + I{\boldsymbol{e}}_\zeta \\ \boldsymbol{e}_\eta =& J_{x\eta z}\left(\nabla z \times \nabla x\right) = {\sigma_{B\text{pol}}} {\boldsymbol{e}}_\theta + \nu{\boldsymbol{e}}_\zeta \\ @@ -853,7 +853,7 @@ Components of `\nabla\times{\boldsymbol{b}}` are [#curvature]_: \begin{aligned} \left(\nabla\times{\boldsymbol{b}}\right)^x =& {\sigma_y}\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}} - R}{B}\right) \\ + R}{B}\right) \\ \left(\nabla\times{\boldsymbol{b}}\right)^y =& -{\sigma_y}\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}R}{B}\right) \\ \left(\nabla\times{\boldsymbol{b}}\right)^z =& \frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) - \frac{{B_{\text{pol}}}{B_{\text{tor}}} R}{{h_\theta}B}{\frac{\partial \nu}{\partial x}} - {\sigma_y}I\frac{{|B_{\text{pol}}|}}{{h_\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}} R}{B}\right) \\ \end{aligned} @@ -919,7 +919,7 @@ we can re-write the above components as: flip, the `x`-coordinate stays the same and the `z`-coordinate flips sign. The `y`-coordinate stays the same, or the `\eta`-coordinate flips sign. - Therefore the x-component of `\nabla\times\boldsymbol{b}` + Therefore the x-component of `\nabla\times\boldsymbol{b}` should flip sign, the `z`-component should not flip sign (product of two sign flips), and the `y`-component should flip sign if '`y`' is `y` and not @@ -1383,7 +1383,7 @@ and \frac{{B_{\text{pol}}}}{B{h_\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}R}{B}\right)\end{aligned} The second and third terms partly cancel, and using -`\sigma_y\sigma_{B\text{pol}}{\frac{\partial I}{\partial y}} = +`\sigma_y\sigma_{B\text{pol}}{\frac{\partial I}{\partial y}} = {\frac{\partial \nu}{\partial x}}` .. math:: @@ -1750,10 +1750,10 @@ Differential geometry ===================== .. warning:: The following are notes from [haeseleer]_. If you are new to this - topic it is strongly suggested to read [haeseleer]_ chapter 2 - instead, as here not all terms are defined, and the discussion of - co- and contra variant components is incomplete. Similiarly, the - notation is based on [haeseleer]_ and not explained in detail. + topic it is strongly suggested to read [haeseleer]_ chapter 2 + instead, as here not all terms are defined, and the discussion of + co- and contra variant components is incomplete. Similiarly, the + notation is based on [haeseleer]_ and not explained in detail. Sets of vectors `\left\{\mathbf{A, B, C}\right\}` and `\left\{\mathbf{a, b, c}\right\}` are reciprocal if @@ -2229,8 +2229,9 @@ The perpendicular Laplacian can therefore be written in divergence form as: +& \frac{1}{J}\frac{\partial}{\partial z}\left[Jg^{xz}\left(\frac{\partial}{\partial x} - \frac{g_{xy}}{g_{yy}}\frac{\partial}{\partial y}\right) + Jg^{zz}\left(\frac{\partial}{\partial z} - \frac{g_{yz}}{g_{yy}}\frac{\partial}{\partial y}\right)\right] \end{aligned} -This form is currently implemented in ``FV::Div_a_Grad_perp`` -(``bout/fv_ops.hxx``) but that operator currently assumes that the +This form is currently implemented in ``FV::Div_a_Grad_perp`` and the +limited variant ``FV::Div_a_Grad_perp_limit`` (both in +``bout/fv_ops.hxx``), but these operators currently assume that the off-diagonal terms `g^{xz}` and `g^{xy}` are zero, which is the case for orthogonal grids with shifted metrics but not in general. diff --git a/manual/sphinx/user_docs/differential_operators.rst b/manual/sphinx/user_docs/differential_operators.rst index bcc04dcc29..e3a2d03d79 100644 --- a/manual/sphinx/user_docs/differential_operators.rst +++ b/manual/sphinx/user_docs/differential_operators.rst @@ -613,6 +613,65 @@ reconstructed values of :math:`f` and :math:`v` remain consistent with the other finite-volume advection operators. +Parallel momentum-flux heating ``Div_par_fvv_heating`` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This operator is a companion to ``FV::Div_par_fvv``. It calculates the heating +associated with the numerical momentum-flux discretisation, and returns a +diagnostic for the kinetic-energy flow through the lower :math:`y` cell face: + +:: + + template + Field3D Div_par_fvv_heating(const Field3D &f_in, const Field3D &v_in, + const Field3D &a, Field3D &flow_ylow, + bool fixflux=true); + + +The returned field is intended for use in parallel momentum and energy systems +where ``FV::Div_par_fvv`` is used in the momentum equation. The extra output +argument ``flow_ylow`` stores the lower-face kinetic-energy flow, including the +area factor. For FCI fields this diagnostic is currently returned as zero. + +As with the other finite-volume parallel operators, ``fixflux=true`` adjusts +boundary fluxes to remain consistent with the imposed boundary values. + + +Perpendicular finite-volume operators +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The finite-volume operators in ``bout/fv_ops.hxx`` also include conservative +perpendicular operators. The baseline operator +``FV::Div_a_Grad_perp(a, f)`` calculates +:math:`\nabla \cdot \left( a \nabla_\perp f \right)` in flux-conservative form: + +:: + + Field3D Div_a_Grad_perp(const Field3D &a, const Field3D &f); + + +There is also a limited variant: + +:: + + template + Field3D Div_a_Grad_perp_limit(const Field3D &a, const Field3D &g, + const Field3D &f); + + +This calculates :math:`\nabla \cdot \left( a g \nabla_\perp f \right)` using a +limited reconstruction of :math:`g` on cell faces. The arguments are: + +* ``a`` - the transport coefficient +* ``g`` - the factor reconstructed and upwinded at cell faces +* ``f`` - the field whose perpendicular gradient is used + +In the :math:`x` direction, the chosen ``CellEdges`` limiter is used to +reconstruct cell-edge values of :math:`g`. In the :math:`y` direction, +first-order upwinding is used. As with ``FV::Div_par``, the template parameter +selects the limiter family described in :ref:`sec-slope-limiters`. + + Example and convergence test ++++++++++++++++++++++++++++ diff --git a/src/mesh/fv_ops.cxx b/src/mesh/fv_ops.cxx index 1bd00d34a9..01c7ca8352 100644 --- a/src/mesh/fv_ops.cxx +++ b/src/mesh/fv_ops.cxx @@ -573,6 +573,11 @@ template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bn template Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, Field3D& flow_ylow, bool fixflux = true); +template Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, + Field3D& flow_ylow, bool fixflux = true); +template Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, + const Field3D& f); template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, bool fixflux = true); @@ -582,6 +587,11 @@ template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bnd template Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, Field3D& flow_ylow, bool fixflux = true); +template Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, + Field3D& flow_ylow, bool fixflux = true); +template Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, + const Field3D& f); template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, bool fixflux = true); @@ -591,6 +601,11 @@ template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bn template Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, Field3D& flow_ylow, bool fixflux = true); +template Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, + Field3D& flow_ylow, bool fixflux = true); +template Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, + const Field3D& f); template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, bool fixflux = true); @@ -600,6 +615,11 @@ template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_ template Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, Field3D& flow_ylow, bool fixflux = true); +template Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, Field3D& flow_ylow, + bool fixflux = true); +template Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, + const Field3D& f); template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, bool fixflux = true); @@ -610,6 +630,11 @@ template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, template Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, Field3D& flow_ylow, bool fixflux = true); +template Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, + Field3D& flow_ylow, bool fixflux = true); +template Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, + const Field3D& f); template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, @@ -621,6 +646,11 @@ template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, template Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, Field3D& flow_ylow, bool fixflux = true); +template Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, + Field3D& flow_ylow, bool fixflux = true); +template Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, + const Field3D& f); template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, bool fixflux = true); @@ -630,5 +660,9 @@ template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bnd template Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, Field3D& flow_ylow, bool fixflux = true); - +template Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, + Field3D& flow_ylow, bool fixflux = true); +template Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, + const Field3D& f); } // Namespace FV