diff --git a/CMakeLists.txt b/CMakeLists.txt index 194f97e845..ff310834d5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -140,6 +140,7 @@ set(BOUT_SOURCES ./include/bout/fieldgroup.hxx ./include/bout/fieldperp.hxx ./include/bout/fv_ops.hxx + ./include/bout/fv_ops_impl.hxx ./include/bout/generic_factory.hxx ./include/bout/globalfield.hxx ./include/bout/globalindexer.hxx diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index 6d803ab28a..14019728d0 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -1,15 +1,8 @@ /************************************************************************** * Describes coordinate systems * - * ChangeLog - * ========= - * - * 2014-11-10 Ben Dudson - * * Created by separating metric from Mesh - * - * ************************************************************************** - * Copyright 2014-2025 BOUT++ contributors + * Copyright 2014-2026 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov * @@ -116,101 +109,115 @@ public: FieldMetric& g_22_yhigh(); // Cell Areas const FieldMetric& cell_area_xlow() const { - if (!_cell_area_xlow.has_value()) { - _compute_cell_area_x(); + if (_cell_area_xlow.has_value()) { + return *_cell_area_xlow; } + _compute_cell_area_x(); ASSERT2(_cell_area_xlow.has_value()); return *_cell_area_xlow; } const FieldMetric& cell_area_xhigh() const { - if (!_cell_area_xhigh.has_value()) { - _compute_cell_area_x(); + if (_cell_area_xhigh.has_value()) { + return *_cell_area_xhigh; } + _compute_cell_area_x(); ASSERT2(_cell_area_xhigh.has_value()); return *_cell_area_xhigh; } const FieldMetric& cell_area_ylow() const { - if (!_cell_area_ylow.has_value()) { - _compute_cell_area_y(); + if (_cell_area_ylow.has_value()) { + return *_cell_area_ylow; } + _compute_cell_area_y(); ASSERT2(_cell_area_ylow.has_value()); return *_cell_area_ylow; } const FieldMetric& cell_area_yhigh() const { - if (!_cell_area_yhigh.has_value()) { - _compute_cell_area_y(); + if (_cell_area_yhigh.has_value()) { + return *_cell_area_yhigh; } + _compute_cell_area_y(); ASSERT2(_cell_area_yhigh.has_value()); return *_cell_area_yhigh; } const FieldMetric& cell_area_zlow() const { - if (!_cell_area_zlow.has_value()) { - _compute_cell_area_z(); + if (_cell_area_zlow.has_value()) { + return *_cell_area_zlow; } + _compute_cell_area_z(); ASSERT2(_cell_area_zlow.has_value()); return *_cell_area_zlow; } const FieldMetric& cell_area_zhigh() const { - if (!_cell_area_zhigh.has_value()) { - _compute_cell_area_z(); + if (_cell_area_zhigh.has_value()) { + return *_cell_area_zhigh; } + _compute_cell_area_z(); ASSERT2(_cell_area_zhigh.has_value()); return *_cell_area_zhigh; } FieldMetric& cell_area_xlow() { - if (!_cell_area_xlow.has_value()) { - _compute_cell_area_x(); + if (_cell_area_xlow.has_value()) { + return *_cell_area_xlow; } + _compute_cell_area_x(); ASSERT2(_cell_area_xlow.has_value()); return *_cell_area_xlow; } FieldMetric& cell_area_xhigh() { - if (!_cell_area_xhigh.has_value()) { - _compute_cell_area_x(); + if (_cell_area_xhigh.has_value()) { + return *_cell_area_xhigh; } + _compute_cell_area_x(); ASSERT2(_cell_area_xhigh.has_value()); return *_cell_area_xhigh; } FieldMetric& cell_area_ylow() { - if (!_cell_area_ylow.has_value()) { - _compute_cell_area_y(); + if (_cell_area_ylow.has_value()) { + return *_cell_area_ylow; } + _compute_cell_area_y(); ASSERT2(_cell_area_ylow.has_value()); return *_cell_area_ylow; } FieldMetric& cell_area_yhigh() { - if (!_cell_area_yhigh.has_value()) { - _compute_cell_area_y(); + if (_cell_area_yhigh.has_value()) { + return *_cell_area_yhigh; } + _compute_cell_area_y(); ASSERT2(_cell_area_yhigh.has_value()); return *_cell_area_yhigh; } FieldMetric& cell_area_zlow() { - if (!_cell_area_zlow.has_value()) { - _compute_cell_area_z(); + if (_cell_area_zlow.has_value()) { + return *_cell_area_zlow; } + _compute_cell_area_z(); ASSERT2(_cell_area_zlow.has_value()); return *_cell_area_zlow; } FieldMetric& cell_area_zhigh() { - if (!_cell_area_zhigh.has_value()) { - _compute_cell_area_z(); + if (_cell_area_zhigh.has_value()) { + return *_cell_area_zhigh; } + _compute_cell_area_z(); ASSERT2(_cell_area_zhigh.has_value()); return *_cell_area_zhigh; } // Cell Volume const FieldMetric& cell_volume() const { - if (!_cell_volume.has_value()) { - _compute_cell_volume(); + if (_cell_volume.has_value()) { + return *_cell_volume; } + _compute_cell_volume(); ASSERT2(_cell_volume.has_value()); return *_cell_volume; } FieldMetric& cell_volume() { - if (!_cell_volume.has_value()) { - _compute_cell_volume(); + if (_cell_volume.has_value()) { + return *_cell_volume; } + _compute_cell_volume(); ASSERT2(_cell_volume.has_value()); return *_cell_volume; } diff --git a/include/bout/derivs.hxx b/include/bout/derivs.hxx index a8d9279378..14bc5c2824 100644 --- a/include/bout/derivs.hxx +++ b/include/bout/derivs.hxx @@ -35,6 +35,7 @@ #include "bout/vector3d.hxx" #include "bout/bout_types.hxx" +#include ////////// FIRST DERIVATIVES ////////// @@ -86,6 +87,13 @@ Field3D DDY(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", const std::string& region = "RGN_NOBNDRY"); +template +std::enable_if_t && !bout::utils::is_Field3D_v, const Field3D> +DDY(const E& expr, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", + const std::string& region = "RGN_NOBNDRY") { + return DDY(Field3D{expr}, outloc, method, region); +} + /// Calculate first partial derivative in Y /// /// \f$\partial / \partial y\f$ diff --git a/include/bout/difops.hxx b/include/bout/difops.hxx index 18220b63ad..0754c1feb6 100644 --- a/include/bout/difops.hxx +++ b/include/bout/difops.hxx @@ -1,11 +1,11 @@ /*!****************************************************************************** * \file difops.hxx - * + * * Differential operators * * Changelog: * - * 2009-01 Ben Dudson + * 2009-01 Ben Dudson * * Added two optional parameters which can be put in any order * These determine the method to use (DIFF_METHOD) * and CELL_LOC location of the result. @@ -15,7 +15,7 @@ * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu * * Contact: Ben Dudson, bd512@york.ac.uk - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -30,7 +30,7 @@ * * You should have received a copy of the GNU Lesser General Public License * along with BOUT++. If not, see . - * + * *******************************************************************************/ #ifndef BOUT_DIFOPS_H @@ -81,7 +81,7 @@ Field3D Grad_parP(const Field3D& apar, const Field3D& f); * \f[ * v\mathbf{b}_0 \cdot \nabla f * \f] - * + * * * @param[in] v The velocity in y direction * @param[in] f The scalar field to be differentiated @@ -175,7 +175,7 @@ inline Field3D Grad2_par2(const Field3D& f, CELL_LOC outloc, DIFF_METHOD method) /*! * Parallel divergence of diffusive flux, K*Grad_par - * + * * \f[ * \nabla \cdot ( \mathbf{b}_0 kY (\mathbf{b}_0 \cdot \nabla) f ) * \f] @@ -195,12 +195,17 @@ Field3D Div_par_K_Grad_par(const Field3D& kY, const Field2D& f, Field3D Div_par_K_Grad_par(const Field3D& kY, const Field3D& f, CELL_LOC outloc = CELL_DEFAULT); +/// Version with energy flow diagnostic +/// For FCI fields, `flow_ylow` is currently returned as zero. +Field3D Div_par_K_Grad_par_mod(const Field3D& k, const Field3D& f, Field3D& flow_ylow, + bool bndry_flux = true); + /*! * Perpendicular Laplacian operator * * This version only includes terms in X and Z, dropping - * derivatives in Y. This is the inverse operation to - * the Laplacian inversion class. + * derivatives in Y. This is the inverse operation to + * the Laplacian inversion class. * * For the full perpendicular Laplacian, use Laplace_perp */ @@ -212,7 +217,7 @@ FieldPerp Delp2(const FieldPerp& f, CELL_LOC outloc = CELL_DEFAULT, bool useFFT /*! * Perpendicular Laplacian, keeping y derivatives * - * + * */ Coordinates::FieldMetric Laplace_perp(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT, @@ -246,18 +251,18 @@ Field2D Laplace_perpXY(const Field2D& A, const Field2D& f); /*! * Terms of form b0 x Grad(phi) dot Grad(A) - * + * */ Coordinates::FieldMetric b0xGrad_dot_Grad(const Field2D& phi, const Field2D& A, CELL_LOC outloc = CELL_DEFAULT); /*! - * Terms of form + * Terms of form * * \f[ * \mathbf{b}_0 \times \nabla \phi \cdot \nabla A * \f] - * + * * @param[in] phi The scalar potential * @param[in] A The field being advected * @param[in] outloc The cell location where the result is defined. By default the same as A. @@ -291,13 +296,13 @@ constexpr BRACKET_METHOD BRACKET_CTU = BRACKET_METHOD::ctu; * \f[ * [f, g] = (1/B) \mathbf{b}_0 \times \nabla f \cdot \nabla g * \f] - * + * * @param[in] f The potential * @param[in] g The field being advected * @param[in] method The method to use * @param[in] outloc The cell location where the result is defined. Default is the same as g * @param[in] solver Pointer to the time integration solver - * + * */ Coordinates::FieldMetric bracket(const Field2D& f, const Field2D& g, BRACKET_METHOD method = BRACKET_STD, diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 323e73f58e..94c599c58c 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -5,7 +5,6 @@ #ifndef BOUT_FV_OPS_H #define BOUT_FV_OPS_H -#include "boutexception.hxx" #include "bout/assert.hxx" #include "bout/bout_types.hxx" #include "bout/build_defines.hxx" @@ -13,27 +12,27 @@ #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" -#include namespace FV { /*! * Div ( a Grad_perp(f) ) -- ∇⊥ ( a ⋅ ∇⊥ f) -- Vorticity */ -Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& x); +Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f); [[deprecated("Please use Div_a_Grad_perp instead")]] inline Field3D -Div_a_Laplace_perp(const Field3D& a, const Field3D& x) { - return Div_a_Grad_perp(a, x); +Div_a_Laplace_perp(const Field3D& a, const Field3D& f) { + return Div_a_Grad_perp(a, f); } /*! * Divergence of a parallel diffusion Div( k * Grad_par(f) ) */ -const Field3D Div_par_K_Grad_par(const Field3D& k, const Field3D& f, - bool bndry_flux = true); +Field3D Div_par_K_Grad_par(const Field3D& k, const Field3D& f, bool bndry_flux = true); /*! * 4th-order derivative in Y, using derivatives @@ -55,7 +54,7 @@ const Field3D Div_par_K_Grad_par(const Field3D& k, const Field3D& f, * * No fluxes through domain boundaries */ -const Field3D D4DY4(const Field3D& d, const Field3D& f); +Field3D D4DY4(const Field3D& d, const Field3D& f); /*! * 4th-order dissipation term @@ -73,195 +72,18 @@ const Field3D D4DY4(const Field3D& d, const Field3D& f); * f_2 | f_1 | f_0 | * f_b */ -const Field3D D4DY4_Index(const Field3D& f, bool bndry_flux = true); - -/*! - * Stencil used for Finite Volume calculations - * which includes cell face values L and R - */ -struct Stencil1D { - // Cell centre values - BoutReal c, m, p, mm, pp; - - // Left and right cell face values - BoutReal L, R; -}; - -/*! - * First order upwind for testing - */ -struct Upwind { - void operator()(Stencil1D& n) { n.L = n.R = n.c; } -}; - -/*! - * Fromm method - */ -struct Fromm { - void operator()(Stencil1D& n) { - n.L = n.c - 0.25 * (n.p - n.m); - n.R = n.c + 0.25 * (n.p - n.m); - } -}; - -/*! - * Second order slope limiter method - * - * Limits slope to minimum absolute value - * of left and right gradients. If at a maximum - * or minimum slope set to zero, i.e. reverts - * to first order upwinding - */ -struct MinMod { - void operator()(Stencil1D& n) { - // Choose the gradient within the cell - // as the minimum (smoothest) solution - const BoutReal slope = _minmod(n.p - n.c, n.c - n.m); - n.L = n.c - 0.5 * slope; - n.R = n.c + 0.5 * slope; - } - -private: - /*! - * Internal helper function for minmod slope limiter - * - * If the inputs have different signs then - * returns zero, otherwise chooses the value - * with the minimum magnitude. - */ - BoutReal _minmod(BoutReal a, BoutReal b) { - if (a * b <= 0.0) { - return 0.0; - } - - if (fabs(a) < fabs(b)) { - return a; - } - return b; - } -}; - -/*! - * Monotonised Central (MC) second order slope limiter (Van Leer) - * - * Limits the slope based on taking the slope with - * the minimum absolute value from central, 2*left and - * 2*right. If any of these slopes have different signs - * then the slope reverts to zero (i.e. 1st-order upwinding). - */ -struct MC { - void operator()(Stencil1D& n) { - const BoutReal slope = minmod(2. * (n.p - n.c), // 2*right difference - 0.5 * (n.p - n.m), // Central difference - 2. * (n.c - n.m)); // 2*left difference - n.L = n.c - 0.5 * slope; - n.R = n.c + 0.5 * slope; - } - -private: - // Return zero if any signs are different - // otherwise return the value with the minimum magnitude - BoutReal minmod(BoutReal a, BoutReal b, BoutReal c) { - // if any of the signs are different, return zero gradient - if ((a * b <= 0.0) || (a * c <= 0.0)) { - return 0.0; - } - - // Return the minimum absolute value - return SIGN(a) * BOUTMIN(fabs(a), fabs(b), fabs(c)); - } -}; - -/*! - * Symmetric Van Albada second order slope limiter - * - * Uses a smooth (differentiable) approximation to `max(a*b, 0)` to avoid - * introducing a kink at extrema, which can be helpful for nonlinear solvers - * and finite-difference Jacobian calculations. - * - * The limited slope is calculated from the left and right differences - * `dl = c - m` and `dr = p - c` as - * - * slope = (pos(dl*dr) * (dl + dr)) / (dl^2 + dr^2) - * - * where `pos(x)` is a smooth approximation to `max(x, 0)`. - */ -struct VanAlbada { - void operator()(Stencil1D& n) { - const BoutReal dl = n.c - n.m; - const BoutReal dr = n.p - n.c; - - const BoutReal denom = dl * dl + dr * dr; - - // Smoothness parameters: - // - keep division well-defined when dl=dr=0 - // - provide a differentiable approximation to max(dl*dr, 0) - const BoutReal eps = 1e-12 * denom + 1e-30; - - const BoutReal ab = dl * dr; - const BoutReal ab_pos = 0.5 * (ab + sqrt(ab * ab + eps * eps)); - - const BoutReal slope = (ab_pos * (dl + dr)) / (denom + eps); - - n.L = n.c - 0.5 * slope; - n.R = n.c + 0.5 * slope; - } -}; - -/*! - * WENO3-JS (Jiang-Shu) reconstruction to cell faces - * - * This is a third-order essentially non-oscillatory reconstruction using two - * candidate second-order polynomials and smoothness-weighted blending. - * - * Unlike TVD slope limiters (e.g. ``MC``), WENO reconstruction is generally - * smooth (differentiable) for all inputs, but it does not enforce strict - * monotonicity. - * - * Uses only the three-point stencil (`m`, `c`, `p`), so it is a drop-in - * replacement anywhere `Stencil1D` is populated with those values. - */ -struct WENO3 { - void operator()(Stencil1D& n) { - // Right face (between c and p): value from cell c (left state at i+1/2) - const BoutReal p0_r = 0.5 * (-n.m + 3.0 * n.c); - const BoutReal p1_r = 0.5 * (n.c + n.p); - - const BoutReal beta0_r = SQ(n.c - n.m); - const BoutReal beta1_r = SQ(n.p - n.c); - - // Left face (between m and c): value from cell c (right state at i-1/2) - const BoutReal p0_l = 0.5 * (-n.p + 3.0 * n.c); - const BoutReal p1_l = 0.5 * (n.m + n.c); - - const BoutReal beta0_l = beta1_r; - const BoutReal beta1_l = beta0_r; - - // Smoothness parameter (scaled to local variation) - const BoutReal eps = 1e-12 * (beta0_r + beta1_r) + 1e-30; - - // Linear weights for WENO3-JS - constexpr BoutReal d0 = 1.0 / 3.0; - constexpr BoutReal d1 = 2.0 / 3.0; - - // Right face weights - const BoutReal a0_r = d0 / SQ(eps + beta0_r); - const BoutReal a1_r = d1 / SQ(eps + beta1_r); - const BoutReal wsum_r = a0_r + a1_r; - const BoutReal w0_r = a0_r / wsum_r; - const BoutReal w1_r = a1_r / wsum_r; - - // Left face weights (mirrored) - const BoutReal a0_l = d0 / SQ(eps + beta0_l); - const BoutReal a1_l = d1 / SQ(eps + beta1_l); - const BoutReal wsum_l = a0_l + a1_l; - const BoutReal w0_l = a0_l / wsum_l; - const BoutReal w1_l = a1_l / wsum_l; - - n.R = w0_r * p0_r + w1_r * p1_r; - n.L = w0_l * p0_l + w1_l * p1_l; - } -}; +Field3D D4DY4_Index(const Field3D& f, bool bndry_flux = true); + +// Forward declarations of flux limiters +// If you want to use your own flux limiter, you need to +// #include to instantiate the templates. +class Upwind; +class Fromm; +class MinMod; +class MC; +class Superbee; +class VanAlbada; +class WENO3; /*! * Communicate fluxes between processors @@ -287,185 +109,7 @@ void communicateFluxes(Field3D& f); /// NB: Uses to/from FieldAligned coordinates template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, - bool fixflux = true) { - - ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); - ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in); - - Mesh* mesh = f_in.getMesh(); - - CellEdges cellboundary; - - ASSERT2(f_in.getDirectionY() == v_in.getDirectionY()); - ASSERT2(f_in.getDirectionY() == wave_speed_in.getDirectionY()); - const bool are_unaligned = - ((f_in.getDirectionY() == YDirectionType::Standard) - and (v_in.getDirectionY() == YDirectionType::Standard) - and (wave_speed_in.getDirectionY() == YDirectionType::Standard)); - - Field3D f = are_unaligned ? toFieldAligned(f_in, "RGN_NOX") : f_in; - Field3D v = are_unaligned ? toFieldAligned(v_in, "RGN_NOX") : v_in; - Field3D wave_speed = - are_unaligned ? toFieldAligned(wave_speed_in, "RGN_NOX") : wave_speed_in; - - Coordinates* coord = f_in.getCoordinates(); - - Field3D result{zeroFrom(f)}; - - // Only need one guard cell, so no need to communicate fluxes - // Instead calculate in guard cells to preserve fluxes - int ys = mesh->ystart - 1; - int ye = mesh->yend + 1; - - for (int i = mesh->xstart; i <= mesh->xend; i++) { - - if (!mesh->firstY(i) || mesh->periodicY(i)) { - // Calculate in guard cell to get fluxes consistent between processors - ys = mesh->ystart - 1; - } else { - // Don't include the boundary cell. Note that this implies special - // handling of boundaries later - ys = mesh->ystart; - } - - if (!mesh->lastY(i) || mesh->periodicY(i)) { - // Calculate in guard cells - ye = mesh->yend + 1; - } else { - // Not in boundary cells - ye = mesh->yend; - } - - for (int j = ys; j <= ye; j++) { - // Pre-calculate factors which multiply fluxes -#if not(BOUT_USE_METRIC_3D) - // For right cell boundaries - BoutReal common_factor = (coord->J(i, j) + coord->J(i, j + 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); - - const BoutReal flux_factor_rc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - const BoutReal flux_factor_rp = - common_factor / (coord->dy(i, j + 1) * coord->J(i, j + 1)); - - // For left cell boundaries - common_factor = (coord->J(i, j) + coord->J(i, j - 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); - - const BoutReal flux_factor_lc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - const BoutReal flux_factor_lm = - common_factor / (coord->dy(i, j - 1) * coord->J(i, j - 1)); -#endif - for (int k = mesh->zstart; k <= mesh->zend; k++) { -#if BOUT_USE_METRIC_3D - // For right cell boundaries - BoutReal common_factor = - (coord->J(i, j, k) + coord->J(i, j + 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); - - BoutReal flux_factor_rc = - common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); - BoutReal flux_factor_rp = - common_factor / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); - - // For left cell boundaries - common_factor = (coord->J(i, j, k) + coord->J(i, j - 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); - - BoutReal flux_factor_lc = - common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); - BoutReal flux_factor_lm = - common_factor / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); -#endif - - //////////////////////////////////////////// - // Reconstruct f at the cell faces - // This calculates s.R and s.L for the Right and Left - // face values on this cell - - // Reconstruct f at the cell faces - Stencil1D s; - s.c = f(i, j, k); - s.m = f(i, j - 1, k); - s.p = f(i, j + 1, k); - - cellboundary(s); // Calculate s.R and s.L - - //////////////////////////////////////////// - // Right boundary - - // Calculate velocity at right boundary (y+1/2) - BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j + 1, k)); - BoutReal flux; - - if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) { - // Last point in domain - - const BoutReal bndryval = 0.5 * (s.c + s.p); - if (fixflux) { - // Use mid-point to be consistent with boundary conditions - flux = bndryval * vpar; - } else { - // Add flux due to difference in boundary values - flux = s.R * vpar + wave_speed(i, j, k) * (s.R - bndryval); - } - } else { - - // Maximum wave speed in the two cells - const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k)); - - if (vpar > amax) { - // Supersonic flow out of this cell - flux = s.R * vpar; - } else if (vpar < -amax) { - // Supersonic flow into this cell - flux = 0.0; - } else { - // Subsonic flow, so a mix of right and left fluxes - flux = s.R * 0.5 * (vpar + amax); - } - } - - result(i, j, k) += flux * flux_factor_rc; - result(i, j + 1, k) -= flux * flux_factor_rp; - - //////////////////////////////////////////// - // Calculate at left boundary - - vpar = 0.5 * (v(i, j, k) + v(i, j - 1, k)); - - if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) { - // First point in domain - const BoutReal bndryval = 0.5 * (s.c + s.m); - if (fixflux) { - // Use mid-point to be consistent with boundary conditions - flux = bndryval * vpar; - } else { - // Add flux due to difference in boundary values - flux = s.L * vpar - wave_speed(i, j, k) * (s.L - bndryval); - } - } else { - - // Maximum wave speed in the two cells - const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k)); - - if (vpar < -amax) { - // Supersonic out of this cell - flux = s.L * vpar; - } else if (vpar > amax) { - // Supersonic into this cell - flux = 0.0; - } else { - flux = s.L * 0.5 * (vpar - amax); - } - } - - result(i, j, k) -= flux * flux_factor_lc; - result(i, j - 1, k) += flux * flux_factor_lm; - } - } - } - return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; -} + bool fixflux = true); /*! * Div ( n * v ) -- Magnetic drifts @@ -480,146 +124,47 @@ Field3D Div_par(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_sp * */ template -Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { - ASSERT1(n_in.getLocation() == v.getLocation()); - ASSERT1_FIELDS_COMPATIBLE(n_in, v.x); - - Mesh* mesh = n_in.getMesh(); - - CellEdges cellboundary; - - Coordinates* coord = n_in.getCoordinates(); - - if (v.covariant) { - // Got a covariant vector instead - throw BoutException("Div_f_v passed a covariant v"); - } - - Field3D result{zeroFrom(n_in)}; - - Field3D vx = v.x; - Field3D vz = v.z; - Field3D n = n_in; - - BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { - // Calculate velocities - const BoutReal vU = 0.25 * (vz[i.zp()] + vz[i]) * (coord->J[i.zp()] + coord->J[i]); - const BoutReal vD = 0.25 * (vz[i.zm()] + vz[i]) * (coord->J[i.zm()] + coord->J[i]); - const BoutReal vL = 0.25 * (vx[i.xm()] + vx[i]) * (coord->J[i.xm()] + coord->J[i]); - const BoutReal vR = 0.25 * (vx[i.xp()] + vx[i]) * (coord->J[i.xp()] + coord->J[i]); - - // X direction - Stencil1D s; - s.c = n[i]; - s.m = n[i.xm()]; - s.mm = n[i.xmm()]; - s.p = n[i.xp()]; - s.pp = n[i.xpp()]; - - cellboundary(s); - - if ((i.x() == mesh->xend) && (mesh->lastX())) { - // At right boundary in X - if (bndry_flux) { - BoutReal flux; - if (vR > 0.0) { - // Flux to boundary - flux = vR * s.R; - } else { - // Flux in from boundary - flux = vR * 0.5 * (n[i.xp()] + n[i]); - } - result[i] += flux / (coord->dx[i] * coord->J[i]); - result[i.xp()] -= flux / (coord->dx[i.xp()] * coord->J[i.xp()]); - } - } else { - // Not at a boundary - if (vR > 0.0) { - // Flux out into next cell - const BoutReal flux = vR * s.R; - result[i] += flux / (coord->dx[i] * coord->J[i]); - result[i.xp()] -= flux / (coord->dx[i.xp()] * coord->J[i.xp()]); - } - } - - // Left side - - if ((i.x() == mesh->xstart) && (mesh->firstX())) { - // At left boundary in X - - if (bndry_flux) { - BoutReal flux; - if (vL < 0.0) { - // Flux to boundary - flux = vL * s.L; - } else { - // Flux in from boundary - flux = vL * 0.5 * (n[i.xm()] + n[i]); - } - result[i] -= flux / (coord->dx[i] * coord->J[i]); - result[i.xm()] += flux / (coord->dx[i.xm()] * coord->J[i.xm()]); - } - } else { - // Not at a boundary - if (vL < 0.0) { - const BoutReal flux = vL * s.L; - result[i] -= flux / (coord->dx[i] * coord->J[i]); - result[i.xm()] += flux / (coord->dx[i.xm()] * coord->J[i.xm()]); - } - } - - /// NOTE: Need to communicate fluxes - - // Z direction - s.m = n[i.zm()]; - s.mm = n[i.zmm()]; - s.p = n[i.zp()]; - s.pp = n[i.zpp()]; - - cellboundary(s); - - if (vU > 0.0) { - const BoutReal flux = vU * s.R; - result[i] += flux / (coord->J[i] * coord->dz[i]); - result[i.zp()] -= flux / (coord->J[i.zp()] * coord->dz[i.zp()]); - } - if (vD < 0.0) { - const BoutReal flux = vD * s.L; - result[i] -= flux / (coord->J[i] * coord->dz[i]); - result[i.zm()] += flux / (coord->J[i.zm()] * coord->dz[i.zm()]); - } - } - - communicateFluxes(result); - - // Y advection - // Currently just using simple centered differences - // so no fluxes need to be exchanged - - n = toFieldAligned(n_in, "RGN_NOX"); - Field3D vy = toFieldAligned(v.y, "RGN_NOX"); - - Field3D yresult = 0.0; - yresult.setDirectionY(YDirectionType::Aligned); - - BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { - // Y velocities on y boundaries - const BoutReal vU = 0.25 * (vy[i] + vy[i.yp()]) * (coord->J[i] + coord->J[i.yp()]); - const BoutReal vD = 0.25 * (vy[i] + vy[i.ym()]) * (coord->J[i] + coord->J[i.ym()]); - - // n (advected quantity) on y boundaries - // Note: Use unshifted n_in variable - const BoutReal nU = 0.5 * (n[i] + n[i.yp()]); - const BoutReal nD = 0.5 * (n[i] + n[i.ym()]); - - yresult[i] = (nU * vU - nD * vD) / (coord->J[i] * coord->dy[i]); - } - return result + fromFieldAligned(yresult, "RGN_NOBNDRY"); -} +Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux); /*! * X-Z Finite Volume diffusion operator */ Field3D Div_Perp_Lap(const Field3D& a, const Field3D& f, CELL_LOC outloc = CELL_DEFAULT); + +/// Finite volume parallel divergence +/// +/// NOTE: Modified version, applies limiter to velocity and field +/// Performs better (smaller overshoots) than Div_par +/// +/// Preserves the sum of f*J*dx*dy*dz over the domain +/// +/// @param[in] f_in The field being advected. +/// This will be reconstructed at cell faces +/// using the given CellEdges method +/// @param[in] v_in The advection velocity. +/// This will be interpolated to cell boundaries +/// using linear interpolation +/// @param[in] wave_speed_in Local maximum speed of all waves in the system at each +// point in space +/// @param[in] fixflux Fix the flux at the boundary to be the value at the +/// midpoint (for boundary conditions) +/// +/// @param[out] flow_ylow Flow at the lower Y cell boundary +/// Already includes area factor * flux +/// For FCI fields this diagnostic is currently set to zero. +template +Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, Field3D& flow_ylow, + bool fixflux = true); + +/// This operator calculates Div_par(f v v) +/// It is used primarily (only?) in the parallel momentum equation. +/// +/// This operator is used rather than Div(f fv) so that the values of +/// f and v are consistent with other advection equations: The product +/// fv is not interpolated to cell boundaries. +template +Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); } // namespace FV #endif // BOUT_FV_OPS_H diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx new file mode 100644 index 0000000000..7da980e584 --- /dev/null +++ b/include/bout/fv_ops_impl.hxx @@ -0,0 +1,986 @@ +/* + Finite-volume discretisation methods. Flux-conservative form + */ + +#ifndef BOUT_FV_OPS_IMPL_H +#define BOUT_FV_OPS_IMPL_H + +#include "bout/assert.hxx" +#include "bout/bout_types.hxx" +#include "bout/boutexception.hxx" +#include "bout/build_defines.hxx" +#include "bout/coordinates.hxx" +#include "bout/field.hxx" +#include "bout/field3d.hxx" +#include "bout/fv_ops.hxx" // NOLINT(unused-includes, misc-include-cleaner) +#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" + +#include + +namespace FV { +/*! + * Stencil used for Finite Volume calculations + * which includes cell face values L and R + */ +struct Stencil1D { + /// Cell centre values + BoutReal c{}; + BoutReal m{}; + BoutReal p{}; + BoutReal mm = BoutNaN; + BoutReal pp = BoutNaN; + + /// Left cell face value + BoutReal L = BoutNaN; + /// Right cell face value + BoutReal R = BoutNaN; +}; + +/*! + * First order upwind for testing + */ +// NB: The templates need to be explicitly instantiated in fv_ops.cxx +struct Upwind { + void operator()(Stencil1D& n) { n.L = n.R = n.c; } +}; + +/*! + * Fromm method + */ +// NB: The templates need to be explicitly instantiated in fv_ops.cxx +struct Fromm { + void operator()(Stencil1D& n) { + n.L = n.c - (0.25 * (n.p - n.m)); + n.R = n.c + (0.25 * (n.p - n.m)); + } +}; + +/*! + * Second order slope limiter method + * + * Limits slope to minimum absolute value + * of left and right gradients. If at a maximum + * or minimum slope set to zero, i.e. reverts + * to first order upwinding + */ +// NB: The templates need to be explicitly instantiated in fv_ops.cxx +struct MinMod { + void operator()(Stencil1D& n) { + // Choose the gradient within the cell + // as the minimum (smoothest) solution + const BoutReal slope = _minmod(n.p - n.c, n.c - n.m); + n.L = n.c - (0.5 * slope); + n.R = n.c + (0.5 * slope); + } + +private: + /*! + * Internal helper function for minmod slope limiter + * + * If the inputs have different signs then + * returns zero, otherwise chooses the value + * with the minimum magnitude. + */ + static BoutReal _minmod(BoutReal a, BoutReal b) { + if (a * b <= 0.0) { + return 0.0; + } + + if (fabs(a) < fabs(b)) { + return a; + } + return b; + } +}; + +/*! + * Monotonised Central (MC) second order slope limiter (Van Leer) + * + * Limits the slope based on taking the slope with + * the minimum absolute value from central, 2*left and + * 2*right. If any of these slopes have different signs + * then the slope reverts to zero (i.e. 1st-order upwinding). + */ +// NB: The templates need to be explicitly instantiated in fv_ops.cxx +struct MC { + void operator()(Stencil1D& n) { + const BoutReal slope = minmod(2. * (n.p - n.c), // 2*right difference + 0.5 * (n.p - n.m), // Central difference + 2. * (n.c - n.m)); // 2*left difference + n.L = n.c - (0.5 * slope); + n.R = n.c + (0.5 * slope); + } + +private: + // Return zero if any signs are different + // otherwise return the value with the minimum magnitude + static BoutReal minmod(BoutReal a, BoutReal b, BoutReal c) { + // if any of the signs are different, return zero gradient + if ((a * b <= 0.0) || (a * c <= 0.0)) { + return 0.0; + } + + // Return the minimum absolute value + return SIGN(a) * BOUTMIN(fabs(a), fabs(b), fabs(c)); + } +}; + +/// Superbee limiter +/// +/// This corresponds to the limiter function +/// φ(r) = max(0, min(2r, 1), min(r,2) +/// +/// The value at cell right (i.e. i + 1/2) is: +/// +/// n.R = n.c - φ(r) (n.c - (n.p + n.c)/2) +/// = n.c + φ(r) (n.p - n.c)/2 +/// +/// Four regimes: +/// a) r < 1/2 -> φ(r) = 2r +/// n.R = n.c + gL +/// b) 1/2 < r < 1 -> φ(r) = 1 +/// n.R = n.c + gR/2 +/// c) 1 < r < 2 -> φ(r) = r +/// n.R = n.c + gL/2 +/// d) 2 < r -> φ(r) = 2 +/// n.R = n.c + gR +/// +/// where the left and right gradients are: +/// gL = n.c - n.m +/// gR = n.p - n.c +/// +// NB: The templates need to be explicitly instantiated in fv_ops.cxx +struct Superbee { + void operator()(Stencil1D& n) { + const BoutReal gL = n.c - n.m; + const BoutReal gR = n.p - n.c; + + // r = gL / gR + // Limiter is φ(r) + if (gL * gR < 0) { + // Different signs => Zero gradient + n.L = n.R = n.c; + } else { + const BoutReal sign = SIGN(gL); + const BoutReal abs_gL = fabs(gL); + const BoutReal abs_gR = fabs(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; + n.R = n.c + half_slope; + } + } +}; + +/*! + * Symmetric Van Albada second order slope limiter + * + * Uses a smooth (differentiable) approximation to `max(a*b, 0)` to avoid + * introducing a kink at extrema, which can be helpful for nonlinear solvers + * and finite-difference Jacobian calculations. + * + * The limited slope is calculated from the left and right differences + * `dl = c - m` and `dr = p - c` as + * + * slope = (pos(dl*dr) * (dl + dr)) / (dl^2 + dr^2) + * + * where `pos(x)` is a smooth approximation to `max(x, 0)`. + */ +// NB: The templates need to be explicitly instantiated in fv_ops.cxx +struct VanAlbada { + void operator()(Stencil1D& n) { + const BoutReal dl = n.c - n.m; + const BoutReal dr = n.p - n.c; + + const BoutReal denom = (dl * dl) + (dr * dr); + + // Smoothness parameters: + // - keep division well-defined when dl=dr=0 + // - provide a differentiable approximation to max(dl*dr, 0) + const BoutReal eps = (1e-12 * denom) + 1e-30; + + const BoutReal ab = dl * dr; + const BoutReal ab_pos = 0.5 * (ab + sqrt((ab * ab) + (eps * eps))); + + const BoutReal slope = (ab_pos * (dl + dr)) / (denom + eps); + + n.L = n.c - (0.5 * slope); + n.R = n.c + (0.5 * slope); + } +}; + +/*! + * WENO3-JS (Jiang-Shu) reconstruction to cell faces + * + * This is a third-order essentially non-oscillatory reconstruction using two + * candidate second-order polynomials and smoothness-weighted blending. + * + * Unlike TVD slope limiters (e.g. ``MC``), WENO reconstruction is generally + * smooth (differentiable) for all inputs, but it does not enforce strict + * monotonicity. + * + * Uses only the three-point stencil (`m`, `c`, `p`), so it is a drop-in + * replacement anywhere `Stencil1D` is populated with those values. + */ +// NB: The templates need to be explicitly instantiated in fv_ops.cxx +struct WENO3 { + void operator()(Stencil1D& n) { + // Right face (between c and p): value from cell c (left state at i+1/2) + const BoutReal p0_r = 0.5 * (-n.m + 3.0 * n.c); + const BoutReal p1_r = 0.5 * (n.c + n.p); + + const BoutReal beta0_r = SQ(n.c - n.m); + const BoutReal beta1_r = SQ(n.p - n.c); + + // Left face (between m and c): value from cell c (right state at i-1/2) + const BoutReal p0_l = 0.5 * (-n.p + 3.0 * n.c); + const BoutReal p1_l = 0.5 * (n.m + n.c); + + const BoutReal beta0_l = beta1_r; + const BoutReal beta1_l = beta0_r; + + // Smoothness parameter (scaled to local variation) + const BoutReal eps = (1e-12 * (beta0_r + beta1_r)) + 1e-30; + + // Linear weights for WENO3-JS + constexpr BoutReal d0 = 1.0 / 3.0; + constexpr BoutReal d1 = 2.0 / 3.0; + + // Right face weights + const BoutReal a0_r = d0 / SQ(eps + beta0_r); + const BoutReal a1_r = d1 / SQ(eps + beta1_r); + const BoutReal wsum_r = a0_r + a1_r; + const BoutReal w0_r = a0_r / wsum_r; + const BoutReal w1_r = a1_r / wsum_r; + + // Left face weights (mirrored) + const BoutReal a0_l = d0 / SQ(eps + beta0_l); + const BoutReal a1_l = d1 / SQ(eps + beta1_l); + const BoutReal wsum_l = a0_l + a1_l; + const BoutReal w0_l = a0_l / wsum_l; + const BoutReal w1_l = a1_l / wsum_l; + + n.R = (w0_r * p0_r) + (w1_r * p1_r); + n.L = (w0_l * p0_l) + (w1_l * p1_l); + } +}; + +/*! + * Communicate fluxes between processors + * Takes values in guard cells, and adds them to cells + */ +void communicateFluxes(Field3D& f); + +/// Finite volume parallel divergence +/// +/// Preserves the sum of f*J*dx*dy*dz over the domain +/// +/// @param[in] f_in The field being advected. +/// This will be reconstructed at cell faces +/// using the given CellEdges method +/// @param[in] v_in The advection velocity. +/// This will be interpolated to cell boundaries +/// using linear interpolation +/// @param[in] wave_speed_in Local maximum speed of all waves in the system at each +// point in space +/// @param[in] fixflux Fix the flux at the boundary to be the value at the +/// midpoint (for boundary conditions) +/// +/// NB: Uses to/from FieldAligned coordinates +template +Field3D Div_par(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, + bool fixflux) { + + ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); + ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in); + + const Mesh* mesh = f_in.getMesh(); + + CellEdges cellboundary; + + ASSERT2(f_in.getDirectionY() == v_in.getDirectionY()); + ASSERT2(f_in.getDirectionY() == wave_speed_in.getDirectionY()); + const bool are_unaligned = + ((f_in.getDirectionY() == YDirectionType::Standard) + and (v_in.getDirectionY() == YDirectionType::Standard) + and (wave_speed_in.getDirectionY() == YDirectionType::Standard)); + + Field3D f = are_unaligned ? toFieldAligned(f_in, "RGN_NOX") : f_in; + Field3D v = are_unaligned ? toFieldAligned(v_in, "RGN_NOX") : v_in; + Field3D wave_speed = + are_unaligned ? toFieldAligned(wave_speed_in, "RGN_NOX") : wave_speed_in; + + Coordinates* coord = f_in.getCoordinates(); + + Field3D result{zeroFrom(f)}; + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + const bool is_periodic_y = mesh->periodicY(i); + const bool is_first_y = mesh->firstY(i); + const bool is_last_y = mesh->lastY(i); + + // Only need one guard cell, so no need to communicate fluxes Instead + // calculate in guard cells to get fluxes consistent between processors, but + // don't include the boundary cell. Note that this implies special handling + // of boundaries later + const int ys = (!is_first_y || is_periodic_y) ? mesh->ystart - 1 : mesh->ystart; + const int ye = (!is_last_y || is_periodic_y) ? mesh->yend + 1 : mesh->yend; + + for (int j = ys; j <= ye; j++) { + // Pre-calculate factors which multiply fluxes +#if not(BOUT_USE_METRIC_3D) + // For right cell boundaries + const BoutReal area_r = coord->cell_area_yhigh()(i, j); + const BoutReal flux_factor_rc = area_r / coord->cell_volume()(i, j); + const BoutReal flux_factor_rp = area_r / coord->cell_volume()(i, j + 1); + // For left cell boundaries + const BoutReal area_l = coord->cell_area_ylow()(i, j); + const BoutReal flux_factor_lc = area_l / coord->cell_volume()(i, j); + const BoutReal flux_factor_lm = area_l / coord->cell_volume()(i, j - 1); +#endif + for (int k = mesh->zstart; k <= mesh->zend; k++) { +#if BOUT_USE_METRIC_3D + // For right cell boundaries + const BoutReal area_r = coord->cell_area_yhigh()(i, j, k); + const BoutReal flux_factor_rc = area_r / coord->cell_volume()(i, j, k); + const BoutReal flux_factor_rp = area_r / coord->cell_volume()(i, j + 1, k); + // For left cell boundaries + const BoutReal area_l = coord->cell_area_ylow()(i, j, k); + const BoutReal flux_factor_lc = area_l / coord->cell_volume()(i, j, k); + const BoutReal flux_factor_lm = area_l / coord->cell_volume()(i, j - 1, k); +#endif + + //////////////////////////////////////////// + // Reconstruct f at the cell faces + // This calculates s.R and s.L for the Right and Left + // face values on this cell + + // Reconstruct f at the cell faces + Stencil1D s; + s.c = f(i, j, k); + s.m = f(i, j - 1, k); + s.p = f(i, j + 1, k); + + cellboundary(s); // Calculate s.R and s.L + + //////////////////////////////////////////// + // Right boundary + + // Calculate velocity at right boundary (y+1/2) + BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j + 1, k)); + BoutReal flux = NAN; + + if (is_last_y && (j == mesh->yend) && !is_periodic_y) { + // Last point in domain + + const BoutReal bndryval = 0.5 * (s.c + s.p); + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = bndryval * vpar; + } else { + // Add flux due to difference in boundary values + flux = (s.R * vpar) + (wave_speed(i, j, k) * (s.R - bndryval)); + } + } else { + + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k)); + + if (vpar > amax) { + // Supersonic flow out of this cell + flux = s.R * vpar; + } else if (vpar < -amax) { + // Supersonic flow into this cell + flux = 0.0; + } else { + // Subsonic flow, so a mix of right and left fluxes + flux = s.R * 0.5 * (vpar + amax); + } + } + + result(i, j, k) += flux * flux_factor_rc; + result(i, j + 1, k) -= flux * flux_factor_rp; + + //////////////////////////////////////////// + // Calculate at left boundary + + vpar = 0.5 * (v(i, j, k) + v(i, j - 1, k)); + + if (is_first_y && (j == mesh->ystart) && !is_periodic_y) { + // First point in domain + const BoutReal bndryval = 0.5 * (s.c + s.m); + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = bndryval * vpar; + } else { + // Add flux due to difference in boundary values + flux = (s.L * vpar) - (wave_speed(i, j, k) * (s.L - bndryval)); + } + } else { + + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k)); + + if (vpar < -amax) { + // Supersonic out of this cell + flux = s.L * vpar; + } else if (vpar > amax) { + // Supersonic into this cell + flux = 0.0; + } else { + flux = s.L * 0.5 * (vpar - amax); + } + } + + result(i, j, k) -= flux * flux_factor_lc; + result(i, j - 1, k) += flux * flux_factor_lm; + } + } + } + return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; +} + +/*! + * Div ( n * v ) -- Magnetic drifts + * + * This uses the expression + * + * Div( A ) = 1/J * d/di ( J * A^i ) + * + * Hence the input vector should be contravariant + * + * Note: Uses to/from FieldAligned + * + */ +template +Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { + ASSERT1(n_in.getLocation() == v.getLocation()); + ASSERT1_FIELDS_COMPATIBLE(n_in, v.x); + + const Mesh* mesh = n_in.getMesh(); + + CellEdges cellboundary; + + Coordinates* coord = n_in.getCoordinates(); + + if (v.covariant) { + // Got a covariant vector instead + throw BoutException("Div_f_v passed a covariant v"); + } + + Field3D result{zeroFrom(n_in)}; + + Field3D vx = v.x; + Field3D vz = v.z; + Field3D n = n_in; + + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + // Calculate velocities + const BoutReal vU = 0.25 * (vz[i.zp()] + vz[i]) * (coord->J[i.zp()] + coord->J[i]); + const BoutReal vD = 0.25 * (vz[i.zm()] + vz[i]) * (coord->J[i.zm()] + coord->J[i]); + const BoutReal vL = 0.25 * (vx[i.xm()] + vx[i]) * (coord->J[i.xm()] + coord->J[i]); + const BoutReal vR = 0.25 * (vx[i.xp()] + vx[i]) * (coord->J[i.xp()] + coord->J[i]); + + // X direction + Stencil1D s; + s.c = n[i]; + s.m = n[i.xm()]; + s.mm = n[i.xmm()]; + s.p = n[i.xp()]; + s.pp = n[i.xpp()]; + + cellboundary(s); + + if ((i.x() == mesh->xend) && (mesh->lastX())) { + // At right boundary in X + if (bndry_flux) { + BoutReal flux = NAN; + if (vR > 0.0) { + // Flux to boundary + flux = vR * s.R; + } else { + // Flux in from boundary + flux = vR * 0.5 * (n[i.xp()] + n[i]); + } + result[i] += flux / (coord->dx[i] * coord->J[i]); + result[i.xp()] -= flux / (coord->dx[i.xp()] * coord->J[i.xp()]); + } + } else { + // Not at a boundary + if (vR > 0.0) { + // Flux out into next cell + const BoutReal flux = vR * s.R; + result[i] += flux / (coord->dx[i] * coord->J[i]); + result[i.xp()] -= flux / (coord->dx[i.xp()] * coord->J[i.xp()]); + } + } + + // Left side + + if ((i.x() == mesh->xstart) && (mesh->firstX())) { + // At left boundary in X + + if (bndry_flux) { + BoutReal flux = NAN; + if (vL < 0.0) { + // Flux to boundary + flux = vL * s.L; + } else { + // Flux in from boundary + flux = vL * 0.5 * (n[i.xm()] + n[i]); + } + result[i] -= flux / (coord->dx[i] * coord->J[i]); + result[i.xm()] += flux / (coord->dx[i.xm()] * coord->J[i.xm()]); + } + } else { + // Not at a boundary + if (vL < 0.0) { + const BoutReal flux = vL * s.L; + result[i] -= flux / (coord->dx[i] * coord->J[i]); + result[i.xm()] += flux / (coord->dx[i.xm()] * coord->J[i.xm()]); + } + } + + /// NOTE: Need to communicate fluxes + + // Z direction + s.m = n[i.zm()]; + s.mm = n[i.zmm()]; + s.p = n[i.zp()]; + s.pp = n[i.zpp()]; + + cellboundary(s); + + if (vU > 0.0) { + const BoutReal flux = vU * s.R; + result[i] += flux / (coord->J[i] * coord->dz[i]); + result[i.zp()] -= flux / (coord->J[i.zp()] * coord->dz[i.zp()]); + } + if (vD < 0.0) { + const BoutReal flux = vD * s.L; + result[i] -= flux / (coord->J[i] * coord->dz[i]); + result[i.zm()] += flux / (coord->J[i.zm()] * coord->dz[i.zm()]); + } + } + + communicateFluxes(result); + + // Y advection + // Currently just using simple centered differences + // so no fluxes need to be exchanged + + n = toFieldAligned(n_in, "RGN_NOX"); + Field3D vy = toFieldAligned(v.y, "RGN_NOX"); + + Field3D yresult = 0.0; + yresult.setDirectionY(YDirectionType::Aligned); + + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + // Y velocities on y boundaries + const BoutReal vU = 0.25 * (vy[i] + vy[i.yp()]) * (coord->J[i] + coord->J[i.yp()]); + const BoutReal vD = 0.25 * (vy[i] + vy[i.ym()]) * (coord->J[i] + coord->J[i.ym()]); + + // n (advected quantity) on y boundaries + // Note: Use unshifted n_in variable + const BoutReal nU = 0.5 * (n[i] + n[i.yp()]); + const BoutReal nD = 0.5 * (n[i] + n[i.ym()]); + + yresult[i] = (nU * vU - nD * vD) / (coord->J[i] * coord->dy[i]); + } + return result + fromFieldAligned(yresult, "RGN_NOBNDRY"); +} + +/// Finite volume parallel divergence +/// +/// NOTE: Modified version, applies limiter to velocity and field +/// Performs better (smaller overshoots) than Div_par +/// +/// Preserves the sum of f*J*dx*dy*dz over the domain +/// +/// @param[in] f_in The field being advected. +/// This will be reconstructed at cell faces +/// using the given CellEdges method +/// @param[in] v_in The advection velocity. +/// This will be interpolated to cell boundaries +/// using linear interpolation +/// @param[in] wave_speed_in Local maximum speed of all waves in the system at each +// point in space +/// @param[in] fixflux Fix the flux at the boundary to be the value at the +/// midpoint (for boundary conditions) +/// +/// @param[out] flow_ylow Flow at the lower Y cell boundary +/// Already includes area factor * flux +template +Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, Field3D& flow_ylow, bool fixflux) { + + Coordinates* coord = f_in.getCoordinates(); + ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); + + if (f_in.isFci()) { + // Use mid-point (cell boundary) averages + + ASSERT1(f_in.hasParallelSlices()); + ASSERT1(v_in.hasParallelSlices()); + + 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(); + + 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(); + + result[i] = (0.25 * (f_in[i] + f_up[iyp]) * (v_in[i] + v_up[iyp]) + * coord->cell_area_yhigh()[i] + - 0.25 * (f_in[i] + f_down[iym]) * (v_in[i] + v_down[iym]) + * coord->cell_area_ylow()[i]) + / coord->cell_volume()[i]; + } + return result; + } + ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in); + + const Mesh* mesh = f_in.getMesh(); + + CellEdges cellboundary; + + ASSERT2(f_in.getDirectionY() == v_in.getDirectionY()); + ASSERT2(f_in.getDirectionY() == wave_speed_in.getDirectionY()); + const bool are_unaligned = + ((f_in.getDirectionY() == YDirectionType::Standard) + and (v_in.getDirectionY() == YDirectionType::Standard) + and (wave_speed_in.getDirectionY() == YDirectionType::Standard)); + + const Field3D f = are_unaligned ? toFieldAligned(f_in, "RGN_NOX") : f_in; + const Field3D v = are_unaligned ? toFieldAligned(v_in, "RGN_NOX") : v_in; + const Field3D wave_speed = + are_unaligned ? toFieldAligned(wave_speed_in, "RGN_NOX") : wave_speed_in; + + Field3D result{zeroFrom(f)}; + flow_ylow = zeroFrom(f); + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + const bool is_periodic_y = mesh->periodicY(i); + const bool is_first_y = mesh->firstY(i); + const bool is_last_y = mesh->lastY(i); + + // Only need one guard cell, so no need to communicate fluxes Instead + // calculate in guard cells to get fluxes consistent between processors, but + // don't include the boundary cell. Note that this implies special handling + // of boundaries later + const int ys = (!is_first_y || is_periodic_y) ? mesh->ystart - 1 : mesh->ystart; + const int ye = (!is_last_y || is_periodic_y) ? mesh->yend + 1 : mesh->yend; + + for (int j = ys; j <= ye; j++) { + // Pre-calculate factors which multiply fluxes +#if not(BOUT_USE_METRIC_3D) + // For right cell boundaries + const BoutReal area_rp = coord->cell_area_yhigh()(i, j); + + const BoutReal flux_factor_rc = area_rp / coord->cell_volume()(i, j); + const BoutReal flux_factor_rp = area_rp / coord->cell_volume()(i, j + 1); + + // For left cell boundaries + const BoutReal area_lc = coord->cell_area_ylow()(i, j); + + const BoutReal flux_factor_lc = area_lc / coord->cell_volume()(i, j); + const BoutReal flux_factor_lm = area_lc / coord->cell_volume()(i, j - 1); +#endif + for (int k = 0; k < mesh->LocalNz; k++) { +#if BOUT_USE_METRIC_3D + // For right cell boundaries + const BoutReal area_rp = coord->cell_area_yhigh()(i, j, k); + + const BoutReal flux_factor_rc = area_rp / coord->cell_volume()(i, j, k); + const BoutReal flux_factor_rp = area_rp / coord->cell_volume()(i, j + 1, k); + + // For left cell boundaries + const BoutReal area_lc = coord->cell_area_ylow()(i, j, k); + + const BoutReal flux_factor_lc = area_lc / coord->cell_volume()(i, j, k); + const BoutReal flux_factor_lm = area_lc / coord->cell_volume()(i, j - 1, k); +#endif + + //////////////////////////////////////////// + // Reconstruct f at the cell faces + // This calculates s.R and s.L for the Right and Left + // face values on this cell + + // Reconstruct f at the cell faces + Stencil1D s{.c = f(i, j, k), .m = f(i, j - 1, k), .p = f(i, j + 1, k)}; + cellboundary(s); // Calculate s.R and s.L + + //////////////////////////////////////////// + // Reconstruct v at the cell faces + Stencil1D sv{.c = v(i, j, k), .m = v(i, j - 1, k), .p = v(i, j + 1, k)}; + cellboundary(sv); // Calculate sv.R and sv.L + + //////////////////////////////////////////// + // Right boundary + + BoutReal flux = BoutNaN; + + if (is_last_y && (j == mesh->yend) && !is_periodic_y) { + // Last point in domain + + // Calculate velocity at right boundary (y+1/2) + const BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j + 1, k)); + + const BoutReal bndryval = 0.5 * (s.c + s.p); + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = bndryval * vpar; + } else { + // Add flux due to difference in boundary values + flux = (s.R * vpar) + (wave_speed(i, j, k) * (s.R - bndryval)); + } + + } else { + // Maximum wave speed in the two cells + 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))); + + flux = s.R * 0.5 * (sv.R + amax); + } + + result(i, j, k) += flux * flux_factor_rc; + result(i, j + 1, k) -= flux * flux_factor_rp; + + flow_ylow(i, j + 1, k) += flux * area_rp; + + //////////////////////////////////////////// + // Calculate at left boundary + + if (is_first_y && (j == mesh->ystart) && !is_periodic_y) { + // First point in domain + const BoutReal bndryval = 0.5 * (s.c + s.m); + const BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j - 1, k)); + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = bndryval * vpar; + } else { + // Add flux due to difference in boundary values + flux = (s.L * vpar) - (wave_speed(i, j, k) * (s.L - bndryval)); + } + } else { + + // Maximum wave speed in the two cells + 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))); + + flux = s.L * 0.5 * (sv.L - amax); + } + + result(i, j, k) -= flux * flux_factor_lc; + result(i, j - 1, k) += flux * flux_factor_lm; + + flow_ylow(i, j, k) += flux * area_lc; + } + } + } + if (are_unaligned) { + flow_ylow = fromFieldAligned(flow_ylow, "RGN_NOBNDRY"); + } + return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; +} + +/// This operator calculates Div_par(f v v) +/// It is used primarily (only?) in the parallel momentum equation. +/// +/// This operator is used rather than Div(f fv) so that the values of +/// f and v are consistent with other advection equations: The product +/// fv is not interpolated to cell boundaries. +template +Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux) { + ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); + const Mesh* mesh = f_in.getMesh(); + const Coordinates* coord = f_in.getCoordinates(); + CellEdges cellboundary; + + if (f_in.isFci()) { + // FCI version, using yup/down fields + ASSERT1(f_in.hasParallelSlices()); + ASSERT1(v_in.hasParallelSlices()); + + const auto& B = coord->Bxy; + const auto& B_up = coord->Bxy.yup(); + const auto& B_down = coord->Bxy.ydown(); + + const auto& f_up = f_in.yup(); + const auto& f_down = f_in.ydown(); + + const auto& v_up = v_in.yup(); + const auto& v_down = v_in.ydown(); + + const auto& g_22 = coord->g_22; + const auto& dy = coord->dy; + + Field3D result{emptyFrom(f_in)}; + BOUT_FOR(i, f_in.getRegion("RGN_NOBNDRY")) { + const auto iyp = i.yp(); + const auto iym = i.ym(); + + // Maximum local wave speed + const BoutReal amax = + BOUTMAX(wave_speed_in[i], fabs(v_in[i]), fabs(v_up[iyp]), fabs(v_down[iym])); + + 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]); + + // Penalty terms. This implementation is very dissipative. + BoutReal penalty = + (amax * (f_in[i] * v_in[i] - f_up[iyp] * v_up[iyp]) / (B[i] + B_up[iyp])) + + (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 (term * penalty > 0) { + penalty = term; + } else { + penalty = -term; + } + } + + result[i] = B[i] * (term + penalty) / (2 * dy[i] * sqrt(g_22[i])); + +#if CHECK > 0 + if (!std::isfinite(result[i])) { + throw BoutException("Non-finite value in Div_par_fvv at {}\n" + "fup {} vup {} fdown {} vdown {} amax {}\n", + "B {} Bup {} Bdown {} dy {} sqrt(g_22} {}", i, f_up[i], + v_up[i], f_down[i], v_down[i], amax, B[i], B_up[i], B_down[i], + dy[i], sqrt(g_22[i])); + } +#endif + } + return result; + } + + ASSERT1(areFieldsCompatible(f_in, wave_speed_in)); + + /// Ensure that f, v and wave_speed are field aligned + Field3D f = toFieldAligned(f_in, "RGN_NOX"); + Field3D v = toFieldAligned(v_in, "RGN_NOX"); + Field3D wave_speed = toFieldAligned(wave_speed_in, "RGN_NOX"); + + Field3D result{zeroFrom(f)}; + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + const bool is_periodic_y = mesh->periodicY(i); + const bool is_first_y = mesh->firstY(i); + const bool is_last_y = mesh->lastY(i); + + // Only need one guard cell, so no need to communicate fluxes Instead + // calculate in guard cells to get fluxes consistent between processors, but + // don't include the boundary cell. Note that this implies special handling + // of boundaries later + const int ys = (!is_first_y || is_periodic_y) ? mesh->ystart - 1 : mesh->ystart; + const int ye = (!is_last_y || is_periodic_y) ? mesh->yend + 1 : mesh->yend; + + for (int j = ys; j <= ye; j++) { + // Pre-calculate factors which multiply fluxes + + for (int k = 0; k < mesh->LocalNz; k++) { + // For right cell boundaries + const BoutReal area_r = coord->cell_area_yhigh()(i, j, k); + + const BoutReal flux_factor_rc = area_r / coord->cell_volume()(i, j, k); + const BoutReal flux_factor_rp = area_r / coord->cell_volume()(i, j + 1, k); + + // For left cell boundaries + const BoutReal area_l = coord->cell_area_ylow()(i, j, k); + + const BoutReal flux_factor_lc = area_l / coord->cell_volume()(i, j, k); + const BoutReal flux_factor_lm = area_l / coord->cell_volume()(i, j - 1, k); + + //////////////////////////////////////////// + // Reconstruct f at the cell faces + // This calculates s.R and s.L for the Right and Left + // face values on this cell + + // Reconstruct f at the cell faces + Stencil1D s{.c = f(i, j, k), .m = f(i, j - 1, k), .p = f(i, j + 1, k)}; + cellboundary(s); // Calculate s.R and s.L + + //////////////////////////////////////////// + // Reconstruct v at the cell faces + Stencil1D sv{.c = v(i, j, k), .m = v(i, j - 1, k), .p = v(i, j + 1, k)}; + cellboundary(sv); + + //////////////////////////////////////////// + // Right boundary + + // Calculate velocity at right boundary (y+1/2) + 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; + + if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) { + // Last point in domain + + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = n_mid_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 + } + } else { + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k), + fabs(sv.c), fabs(sv.p)); + + flux = s.R * 0.5 * (sv.R + amax) * sv.R; + } + + result(i, j, k) += flux * flux_factor_rc; + result(i, j + 1, k) -= flux * flux_factor_rp; + + //////////////////////////////////////////// + // Calculate at left boundary + + const BoutReal v_mid_l = 0.5 * (sv.c + sv.m); + const BoutReal n_mid_l = 0.5 * (s.c + s.m); + + if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) { + // First point in domain + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = n_mid_l * v_mid_l * v_mid_l; + } 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)); + } + } else { + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k), + fabs(sv.c), fabs(sv.m)); + + flux = s.L * 0.5 * (sv.L - amax) * sv.L; + } + + result(i, j, k) -= flux * flux_factor_lc; + result(i, j - 1, k) += flux * flux_factor_lm; + } + } + } + return fromFieldAligned(result, "RGN_NOBNDRY"); +} +} // namespace FV +#endif // BOUT_FV_OPS_H diff --git a/manual/sphinx/developer_docs/mesh.rst b/manual/sphinx/developer_docs/mesh.rst index ff3e40d4b0..b1861fc416 100644 --- a/manual/sphinx/developer_docs/mesh.rst +++ b/manual/sphinx/developer_docs/mesh.rst @@ -307,3 +307,26 @@ because they are needed in a lot of the code. They shouldn’t change after initialisation, unless the physics model starts doing fancy things with deforming meshes. In that case it is up to the user to ensure they are updated. + +.. _sec-derived-geometric-quantities: + +Derived geometric quantities +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +`Coordinates` also provides several quantities derived from the metric tensor +and Jacobian, rather than requiring each operator to reconstruct them locally. +Examples include: + +* `g_22` evaluated at the lower and upper `y` cell faces +* cell-face areas in the `x`, `y`, and `z` directions +* cell volumes + +These are used by conservative operators, especially the finite-volume +operators documented in :ref:`sec-finite-volume-operators`, where fluxes are +naturally expressed as a face area multiplied by a face flux and divided by a +cell volume. + +In the current implementation these quantities are computed lazily from the +current metric data and then cached for reuse. This keeps the operator code +closer to the discrete flux expressions while centralising the geometry +construction in `Coordinates`. diff --git a/manual/sphinx/user_docs/differential_operators.rst b/manual/sphinx/user_docs/differential_operators.rst index 71dbeb4f9e..3de13001ab 100644 --- a/manual/sphinx/user_docs/differential_operators.rst +++ b/manual/sphinx/user_docs/differential_operators.rst @@ -503,6 +503,8 @@ neglected if :math:`g_{xy}` and :math:`g_{yz}` are non-zero. An example of usage of the brackets can be found in for example ``examples/MMS/advection`` or ``examples/blob2d``. +.. _sec-finite-volume-operators: + Finite volume, conservative finite difference methods ----------------------------------------------------- @@ -537,6 +539,8 @@ The methods can be used by including the Some methods (those with templates) are defined in the header, but others are defined in :doc:`src/mesh/fv_ops.cxx<../_breathe_autogen/file/fv__ops_8cxx>`. +These operators use derived geometric quantities such as cell-face areas and +cell volumes; see :ref:`sec-derived-geometric-quantities`. Parallel divergence ``Div_par`` @@ -570,6 +574,45 @@ be changed at compile time e.g:: A list of available limiters is given in section :ref:`sec-slope-limiters` below. +Modified parallel divergence ``Div_par_mod`` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This is a modified version of ``FV::Div_par`` which applies the limiter to both +the advected field and the velocity. This typically gives smaller overshoots +than ``FV::Div_par`` for the same limiter choice. + +:: + + template + Field3D Div_par_mod(const Field3D &f_in, const Field3D &v_in, + const Field3D &a, Field3D &flow_ylow, + bool fixflux=true); + + +The extra output argument ``flow_ylow`` stores the flow through the lower +:math:`y` cell boundary, including the area factor. This can be useful as a +diagnostic in energy or flux budgets. For FCI fields this diagnostic is +currently returned as zero. + + +Parallel momentum flux ``Div_par_fvv`` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This operator calculates the divergence of :math:`f v v`, and is mainly useful +for the parallel momentum equation: + +:: + + template + Field3D Div_par_fvv(const Field3D &f_in, const Field3D &v_in, + const Field3D &a, bool fixflux=true); + + +This is provided separately rather than forming :math:`fv` first, so that the +reconstructed values of :math:`f` and :math:`v` remain consistent with the +other finite-volume advection operators. + + Example and convergence test ++++++++++++++++++++++++++++ @@ -620,6 +663,24 @@ This is done by calculating the flux :math:`k\partial_{||}\left(f\right)` on cel using central differencing. +Parallel diffusion with flow diagnostic ``Div_par_K_Grad_par_mod`` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This variant of the parallel diffusion operator also returns a lower-face flow +diagnostic: + +:: + + Field3D Div_par_K_Grad_par_mod(const Field3D &k, const Field3D &f, + Field3D &flow_ylow, + bool bndry_flux=true); + + +As for ``Div_par_mod``, ``flow_ylow`` stores the lower :math:`y` boundary flow +including the area factor. For FCI fields this diagnostic is currently returned +as zero. + + Advection in 3D ~~~~~~~~~~~~~~~ @@ -667,6 +728,10 @@ values. Several slope limiters are defined in ``fv_ops.hxx``: to ``MinMod``. It has smaller dissipation than ``MinMod`` so is the default. +* ``Superbee`` - A more compressive TVD limiter than ``MC`` or ``MinMod``. + It tends to sharpen steep gradients and contacts more aggressively, at the + cost of being less smooth. + * ``VanAlbada`` - A smooth (differentiable) symmetric slope limiter which avoids piecewise branches at extrema. This can be useful for nonlinear solvers and finite-difference Jacobian calculations. diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 057ffa65b5..e7a7579d45 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -2015,17 +2015,23 @@ const Coordinates::FieldMetric& Coordinates::g_22_ylow() const { if (_g_22_ylow.has_value()) { return *_g_22_ylow; } - _g_22_ylow.emplace(emptyFrom(g_22)); - //_g_22_ylow->setLocation(CELL_YLOW); - auto* mesh = Bxy.getMesh(); - if (Bxy.isFci()) { - if (mesh->get(_g_22_ylow.value(), "g_22_cell_ylow", 0.0, false) != 0) { - throw BoutException("The grid file does not contain `g_22_cell_ylow`."); - } - } else { - ASSERT0(mesh->ystart > 0); - BOUT_FOR(i, g_22.getRegion("RGN_NOY")) { - _g_22_ylow.value()[i] = SQ(0.5 * (std::sqrt(g_22[i]) + std::sqrt(g_22[i.ym()]))); + BOUT_OMP_SAFE(critical) + { + if (!_g_22_ylow.has_value()) { + _g_22_ylow.emplace(emptyFrom(g_22)); + //_g_22_ylow->setLocation(CELL_YLOW); + auto* mesh = Bxy.getMesh(); + if (Bxy.isFci()) { + if (mesh->get(_g_22_ylow.value(), "g_22_cell_ylow", 0.0, false) != 0) { + throw BoutException("The grid file does not contain `g_22_cell_ylow`."); + } + } else { + ASSERT0(mesh->ystart > 0); + BOUT_FOR(i, g_22.getRegion("RGN_NOY")) { + _g_22_ylow.value()[i] = + SQ(0.5 * (std::sqrt(g_22[i]) + std::sqrt(g_22[i.ym()]))); + } + } } } return g_22_ylow(); @@ -2035,99 +2041,127 @@ const Coordinates::FieldMetric& Coordinates::g_22_yhigh() const { if (_g_22_yhigh.has_value()) { return *_g_22_yhigh; } - _g_22_yhigh.emplace(emptyFrom(g_22)); - auto* mesh = Bxy.getMesh(); - if (Bxy.isFci()) { - if (mesh->get(_g_22_yhigh.value(), "g_22_cell_yhigh", 0.0, false) != 0) { - throw BoutException("The grid file does not contain `g_22_cell_yhigh`."); - } - } else { - ASSERT0(mesh->ystart > 0); - BOUT_FOR(i, g_22.getRegion("RGN_NOY")) { - _g_22_yhigh.value()[i] = SQ(0.5 * (std::sqrt(g_22[i]) + std::sqrt(g_22[i.yp()]))); + BOUT_OMP_SAFE(critical) + { + if (!_g_22_yhigh.has_value()) { + _g_22_yhigh.emplace(emptyFrom(g_22)); + auto* mesh = Bxy.getMesh(); + if (Bxy.isFci()) { + if (mesh->get(_g_22_yhigh.value(), "g_22_cell_yhigh", 0.0, false) != 0) { + throw BoutException("The grid file does not contain `g_22_cell_yhigh`."); + } + } else { + ASSERT0(mesh->ystart > 0); + BOUT_FOR(i, g_22.getRegion("RGN_NOY")) { + _g_22_yhigh.value()[i] = + SQ(0.5 * (std::sqrt(g_22[i]) + std::sqrt(g_22[i.yp()]))); + } + } } } return g_22_yhigh(); } void Coordinates::_compute_cell_area_x() const { - const FieldMetric area_centre = sqrt(g_22 * g_33 - SQ(g_23)) * dy * dz; - _cell_area_xlow.emplace(emptyFrom(area_centre)); - _cell_area_xhigh.emplace(emptyFrom(area_centre)); - // We cannot setLocation, as that would trigger the computation of staggered - // metrics. - auto* mesh = Bxy.getMesh(); - ASSERT0(mesh->xstart > 0); - BOUT_FOR(i, area_centre.getRegion("RGN_NOX")) { - (*_cell_area_xlow)[i] = 0.5 * (area_centre[i] + area_centre[i.xm()]); - (*_cell_area_xhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.xp()]); + BOUT_OMP_SAFE(critical) + { + if (!_cell_area_xlow.has_value()) { + const FieldMetric area_centre = J / sqrt(g_11) * dy * dz; + _cell_area_xlow.emplace(emptyFrom(area_centre)); + _cell_area_xhigh.emplace(emptyFrom(area_centre)); + // We cannot setLocation, as that would trigger the computation of staggered + // metrics. + auto* mesh = Bxy.getMesh(); + ASSERT0(mesh->xstart > 0); + BOUT_FOR(i, area_centre.getRegion("RGN_NOX")) { + (*_cell_area_xlow)[i] = 0.5 * (area_centre[i] + area_centre[i.xm()]); + (*_cell_area_xhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.xp()]); + } + } } } void Coordinates::_compute_cell_area_y() const { - auto* mesh = Bxy.getMesh(); - if (g_11.isFci()) { - const FieldMetric jxz_centre = sqrt(g_11 * g_33 - SQ(g_13)); - auto jxz_ylow = emptyFrom(jxz_centre); - auto jxz_yhigh = emptyFrom(jxz_centre); - - auto By_c = emptyFrom(jxz_centre); - auto By_h = emptyFrom(jxz_yhigh); - auto By_l = emptyFrom(jxz_ylow); - if (mesh->get(By_c, "By", 0.0, false, CELL_CENTRE) != 0) { - throw BoutException("The grid file does not contain `By`."); - } - if (mesh->get(By_l, "By_cell_ylow", 0.0, false) != 0) { - throw BoutException("The grid file does not contain `By_cell_ylow`."); - } - if (mesh->get(By_h, "By_cell_yhigh", 0.0, false) != 0) { - throw BoutException("The grid file does not contain `By_cell_yhigh`."); - } - BOUT_FOR(i, By_c.getRegion("RGN_NOY")) { - jxz_ylow[i] = By_c[i] / By_l[i] * jxz_centre[i]; - jxz_yhigh[i] = By_c[i] / By_h[i] * jxz_centre[i]; - } - ASSERT3(isUniform(dx, true, "RGN_ALL")); - ASSERT2(isUniform(dx, false, "RGN_ALL")); - ASSERT3(isUniform(dz, true, "RGN_ALL")); - ASSERT2(isUniform(dz, false, "RGN_ALL")); - _cell_area_ylow.emplace(jxz_ylow * dx * dz); - _cell_area_yhigh.emplace(jxz_yhigh * dx * dz); - } else { - // Field aligned - const FieldMetric area_centre = sqrt(g_11 * g_33 - SQ(g_13)) * dx * dz; - _cell_area_ylow.emplace(emptyFrom(area_centre)); - _cell_area_yhigh.emplace(emptyFrom(area_centre)); - // We cannot setLocation, as that would trigger the computation of staggered - // metrics. - BOUT_FOR(i, mesh->getRegion("RGN_ALL")) { - if (i.y() > 0) { - (*_cell_area_ylow)[i] = 0.5 * (area_centre[i] + area_centre[i.ym()]); - } else { - (*_cell_area_ylow)[i] = BoutNaN; - } - if (i.y() < mesh->LocalNy - 1) { - (*_cell_area_yhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.yp()]); + BOUT_OMP_SAFE(critical) + { + if (!_cell_area_ylow.has_value()) { + auto* mesh = Bxy.getMesh(); + if (g_11.isFci()) { + const FieldMetric jxz_centre = J / sqrt(g_22); + auto jxz_ylow = emptyFrom(jxz_centre); + auto jxz_yhigh = emptyFrom(jxz_centre); + + auto By_c = emptyFrom(jxz_centre); + auto By_h = emptyFrom(jxz_yhigh); + auto By_l = emptyFrom(jxz_ylow); + if (mesh->get(By_c, "By", 0.0, false, CELL_CENTRE) != 0) { + throw BoutException("The grid file does not contain `By`."); + } + if (mesh->get(By_l, "By_cell_ylow", 0.0, false) != 0) { + throw BoutException("The grid file does not contain `By_cell_ylow`."); + } + if (mesh->get(By_h, "By_cell_yhigh", 0.0, false) != 0) { + throw BoutException("The grid file does not contain `By_cell_yhigh`."); + } + BOUT_FOR(i, By_c.getRegion("RGN_NOY")) { + jxz_ylow[i] = By_c[i] / By_l[i] * jxz_centre[i]; + jxz_yhigh[i] = By_c[i] / By_h[i] * jxz_centre[i]; + } + ASSERT3(isUniform(dx, true, "RGN_ALL")); + ASSERT2(isUniform(dx, false, "RGN_ALL")); + ASSERT3(isUniform(dz, true, "RGN_ALL")); + ASSERT2(isUniform(dz, false, "RGN_ALL")); + _cell_area_ylow.emplace(jxz_ylow * dx * dz); + _cell_area_yhigh.emplace(jxz_yhigh * dx * dz); } else { - (*_cell_area_yhigh)[i] = BoutNaN; + // Field aligned + const FieldMetric area_centre = J / sqrt(g_22) * dx * dz; + _cell_area_ylow.emplace(emptyFrom(area_centre)); + _cell_area_yhigh.emplace(emptyFrom(area_centre)); + // We cannot setLocation, as that would trigger the computation of staggered + // metrics. + BOUT_FOR(i, mesh->getRegion("RGN_ALL")) { + if (i.y() > 0) { + (*_cell_area_ylow)[i] = 0.5 * (area_centre[i] + area_centre[i.ym()]); + } else { + (*_cell_area_ylow)[i] = BoutNaN; + } + if (i.y() < mesh->LocalNy - 1) { + (*_cell_area_yhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.yp()]); + } else { + (*_cell_area_yhigh)[i] = BoutNaN; + } + } } } } } void Coordinates::_compute_cell_area_z() const { - const FieldMetric area_centre = sqrt(g_11 * g_22 - SQ(g_12)) * dx * dy; - _cell_area_zlow.emplace(emptyFrom(area_centre)); - _cell_area_zhigh.emplace(emptyFrom(area_centre)); - // We cannot setLocation, as that would trigger the computation of staggered - // metrics. - BOUT_FOR(i, area_centre.getRegion("RGN_NOZ")) { - (*_cell_area_zlow)[i] = 0.5 * (area_centre[i] + area_centre[i.zm()]); - (*_cell_area_zhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.zp()]); + BOUT_OMP_SAFE(critical) + { + if (!_cell_area_zlow.has_value()) { + const FieldMetric area_centre = J / sqrt(g_33) * dx * dy; + _cell_area_zlow.emplace(emptyFrom(area_centre)); + _cell_area_zhigh.emplace(emptyFrom(area_centre)); + // We cannot setLocation, as that would trigger the computation of staggered + // metrics. + BOUT_FOR(i, area_centre.getRegion("RGN_NOZ")) { + (*_cell_area_zlow)[i] = 0.5 * (area_centre[i] + area_centre[i.zm()]); + (*_cell_area_zhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.zp()]); + } + } } } -void Coordinates::_compute_cell_volume() const { _cell_volume.emplace(J * dx * dy * dz); } +void Coordinates::_compute_cell_volume() const { + BOUT_OMP_SAFE(critical) + { + if (!_cell_volume.has_value()) { + _cell_volume.emplace(J * dx * dy * dz); + } + } +} std::shared_ptr Coordinates::makeYBoundary(YBndryType type) const { return std::make_shared(type, localoptions, *localmesh); diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 3f7aec08bf..8ecd9d64ff 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -36,6 +36,7 @@ #include #include #include +#include #include #include #include @@ -373,6 +374,108 @@ Field3D Div_par_K_Grad_par(const Field3D& kY, const Field3D& f, CELL_LOC outloc) + Div_par(kY, outloc) * Grad_par(f, outloc); } +Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, Field3D& flow_ylow, + bool bndry_flux) { + ASSERT2(Kin.getLocation() == fin.getLocation()); + + const Mesh* mesh = Kin.getMesh(); + const Coordinates* coord = fin.getCoordinates(); + + if (Kin.isFci()) { + ASSERT1(Kin.hasParallelSlices()); + ASSERT1(fin.hasParallelSlices()); + // Using parallel slices. + // Note: Y slices may use different coordinate systems + // -> Only B, dy and g_22 can be used in yup/ydown + // Others (e.g J) may not be averaged between y planes. + + const auto& K_up = Kin.yup(); + const auto& K_down = Kin.ydown(); + + const auto& f_up = fin.yup(); + const auto& f_down = fin.ydown(); + + Field3D result{zeroFrom(fin)}; + flow_ylow = zeroFrom(fin); + + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + const auto iyp = i.yp(); + const auto iym = i.ym(); + + // Upper cell edge + const BoutReal c_up = 0.5 * (Kin[i] + K_up[iyp]); // K at the upper boundary + const BoutReal J_up = + 0.5 * (coord->J[i] + coord->J.yup()[iyp]); // Jacobian at boundary + const BoutReal g_22_up = 0.5 * (coord->g_22[i] + coord->g_22.yup()[iyp]); + const BoutReal gradient_up = + 2. * (f_up[iyp] - fin[i]) / (coord->dy[i] + coord->dy.yup()[iyp]); + + const BoutReal flux_up = c_up * J_up * gradient_up / g_22_up; + + // Lower cell edge + const BoutReal c_down = 0.5 * (Kin[i] + K_down[iym]); // K at the lower boundary + const BoutReal J_down = + 0.5 * (coord->J[i] + coord->J.ydown()[iym]); // Jacobian at boundary + const BoutReal g_22_down = 0.5 * (coord->g_22[i] + coord->g_22.ydown()[iym]); + const BoutReal gradient_down = + 2. * (fin[i] - f_down[iym]) / (coord->dy[i] + coord->dy.ydown()[iym]); + + const BoutReal flux_down = c_down * J_down * gradient_down / g_22_down; + + result[i] = (flux_up - flux_down) / (coord->dy[i] * coord->J[i]); + } + + return result; + } + + // Calculate in field-aligned coordinates + const auto& K = toFieldAligned(Kin, "RGN_NOX"); + const auto& f = toFieldAligned(fin, "RGN_NOX"); + + Field3D result{zeroFrom(f)}; + flow_ylow = zeroFrom(f); + + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + // Calculate flux at upper surface + const auto ix = i.x(); + const auto iy = i.y(); + const auto iyp = i.yp(); + const auto iym = i.ym(); + + const bool is_periodic_y = mesh->periodicY(ix); + + if (bndry_flux || is_periodic_y || !mesh->lastY(ix) || (iy != mesh->yend)) { + const BoutReal c = 0.5 * (K[i] + K[iyp]); // K at the upper boundary + const BoutReal J = 0.5 * (coord->J[i] + coord->J[iyp]); // Jacobian at boundary + const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); + const BoutReal gradient = 2. * (f[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); + + const BoutReal flux = c * J * gradient / g_22; + + result[i] += flux / (coord->dy[i] * coord->J[i]); + } + + // Calculate flux at lower surface + if (bndry_flux || is_periodic_y || !mesh->firstY(ix) || (iy != mesh->ystart)) { + const BoutReal c = 0.5 * (K[i] + K[iym]); // K at the lower boundary + const BoutReal J = 0.5 * (coord->J[i] + coord->J[iym]); // Jacobian at boundary + const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); + const BoutReal gradient = 2. * (f[i] - f[iym]) / (coord->dy[i] + coord->dy[iym]); + + const BoutReal flux = c * J * gradient / g_22; + + result[i] -= flux / (coord->dy[i] * coord->J[i]); + flow_ylow[i] = -flux * coord->dx[i] * coord->dz[i]; + } + } + + // Shifted to field aligned coordinates, so need to shift back + result = fromFieldAligned(result, "RGN_NOBNDRY"); + flow_ylow = fromFieldAligned(flow_ylow); + + return result; +} + /******************************************************************************* * Delp2 * perpendicular Laplacian operator diff --git a/src/mesh/fv_ops.cxx b/src/mesh/fv_ops.cxx index fab8beb794..1bd00d34a9 100644 --- a/src/mesh/fv_ops.cxx +++ b/src/mesh/fv_ops.cxx @@ -1,7 +1,17 @@ -#include -#include -#include -#include +#include "bout/fv_ops.hxx" +#include "bout/fv_ops_impl.hxx" + +#include "bout/assert.hxx" +#include "bout/bout_types.hxx" +#include "bout/boutexception.hxx" +#include "bout/build_config.hxx" +#include "bout/coordinates.hxx" +#include "bout/field2d.hxx" +#include "bout/field3d.hxx" +#include "bout/globals.hxx" +#include "bout/msg_stack.hxx" +#include "bout/region.hxx" +#include "bout/utils.hxx" namespace { template @@ -33,28 +43,19 @@ Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f) { // Flux in x - int xs = mesh->xstart - 1; - int xe = mesh->xend; - - /* - if(mesh->firstX()) - xs += 1; - */ - /* - if(mesh->lastX()) - xe -= 1; - */ + const int xs = mesh->xstart - 1; + const int xe = mesh->xend; for (int i = xs; i <= xe; i++) { for (int j = mesh->ystart; j <= mesh->yend; j++) { for (int k = mesh->zstart; k <= mesh->zend; 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)); + const 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)); @@ -178,14 +179,13 @@ Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f) { return result; } -const Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, - bool bndry_flux) { +Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, bool bndry_flux) { ASSERT2(Kin.getLocation() == fin.getLocation()); - Mesh* mesh = Kin.getMesh(); + const Mesh* mesh = Kin.getMesh(); - bool use_parallel_slices = (Kin.hasParallelSlices() && fin.hasParallelSlices()); + const bool use_parallel_slices = (Kin.hasParallelSlices() && fin.hasParallelSlices()); const auto& K = use_parallel_slices ? Kin : toFieldAligned(Kin, "RGN_NOX"); const auto& f = use_parallel_slices ? fin : toFieldAligned(fin, "RGN_NOX"); @@ -209,13 +209,13 @@ const Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, if (bndry_flux || mesh->periodicY(i.x()) || !mesh->lastY(i.x()) || (i.y() != mesh->yend)) { - BoutReal c = 0.5 * (K[i] + Kup[iyp]); // K at the upper boundary - BoutReal J = 0.5 * (coord->J[i] + coord->J[iyp]); // Jacobian at boundary - BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); + const BoutReal c = 0.5 * (K[i] + Kup[iyp]); // K at the upper boundary + const BoutReal J = 0.5 * (coord->J[i] + coord->J[iyp]); // Jacobian at boundary + const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); - BoutReal gradient = 2. * (fup[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); + const BoutReal gradient = 2. * (fup[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); - BoutReal flux = c * J * gradient / g_22; + const BoutReal flux = c * J * gradient / g_22; result[i] += flux / (coord->dy[i] * coord->J[i]); } @@ -223,14 +223,15 @@ const Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, // Calculate flux at lower surface if (bndry_flux || mesh->periodicY(i.x()) || !mesh->firstY(i.x()) || (i.y() != mesh->ystart)) { - BoutReal c = 0.5 * (K[i] + Kdown[iym]); // K at the lower boundary - BoutReal J = 0.5 * (coord->J[i] + coord->J[iym]); // Jacobian at boundary + const BoutReal c = 0.5 * (K[i] + Kdown[iym]); // K at the lower boundary + const BoutReal J = 0.5 * (coord->J[i] + coord->J[iym]); // Jacobian at boundary - BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); + const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); - BoutReal gradient = 2. * (f[i] - fdown[iym]) / (coord->dy[i] + coord->dy[iym]); + const BoutReal gradient = + 2. * (f[i] - fdown[iym]) / (coord->dy[i] + coord->dy[iym]); - BoutReal flux = c * J * gradient / g_22; + const BoutReal flux = c * J * gradient / g_22; result[i] -= flux / (coord->dy[i] * coord->J[i]); } @@ -244,10 +245,10 @@ const Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, return result; } -const Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { +Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { ASSERT1_FIELDS_COMPATIBLE(d_in, f_in); - Mesh* mesh = d_in.getMesh(); + const Mesh* mesh = d_in.getMesh(); Coordinates* coord = f_in.getCoordinates(); @@ -263,9 +264,9 @@ const Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { for (int i = mesh->xstart; i <= mesh->xend; i++) { // Check for boundaries - bool yperiodic = mesh->periodicY(i); - bool has_upper_boundary = !yperiodic && mesh->lastY(i); - bool has_lower_boundary = !yperiodic && mesh->firstY(i); + const bool yperiodic = mesh->periodicY(i); + const bool has_upper_boundary = !yperiodic && mesh->lastY(i); + const bool has_lower_boundary = !yperiodic && mesh->firstY(i); // Always calculate fluxes at upper Y cell boundary const int ystart = @@ -281,15 +282,15 @@ const Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { for (int j = ystart; j <= yend; j++) { for (int k = mesh->zstart; k <= mesh->zend; k++) { - BoutReal dy3 = SQ(coord->dy(i, j, k)) * coord->dy(i, j, k); + const BoutReal dy3 = SQ(coord->dy(i, j, k)) * coord->dy(i, j, k); // 3rd derivative at upper boundary - BoutReal d3fdy3 = + const BoutReal d3fdy3 = (f(i, j + 2, k) - 3. * f(i, j + 1, k) + 3. * f(i, j, k) - f(i, j - 1, k)) / dy3; - BoutReal flux = 0.5 * (d(i, j, k) + d(i, j + 1, k)) - * (coord->J(i, j, k) + coord->J(i, j + 1, k)) * d3fdy3; + const BoutReal flux = 0.5 * (d(i, j, k) + d(i, j + 1, k)) + * (coord->J(i, j, k) + coord->J(i, j + 1, k)) * d3fdy3; result(i, j, k) += flux / (coord->J(i, j, k) * coord->dy(i, j, k)); result(i, j + 1, k) -= flux / (coord->J(i, j + 1, k) * coord->dy(i, j + 1, k)); @@ -301,8 +302,8 @@ const Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; } -const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { - Mesh* mesh = f_in.getMesh(); +Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { + const Mesh* mesh = f_in.getMesh(); // Convert to field aligned coordinates const bool is_unaligned = (f_in.getDirectionY() == YDirectionType::Standard); @@ -313,10 +314,10 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { Coordinates* coord = f_in.getCoordinates(); for (int i = mesh->xstart; i <= mesh->xend; i++) { - bool yperiodic = mesh->periodicY(i); + const bool yperiodic = mesh->periodicY(i); - bool has_upper_boundary = !yperiodic && mesh->lastY(i); - bool has_lower_boundary = !yperiodic && mesh->firstY(i); + const bool has_upper_boundary = !yperiodic && mesh->lastY(i); + const bool has_lower_boundary = !yperiodic && mesh->firstY(i); for (int j = mesh->ystart; j <= mesh->yend; j++) { @@ -341,8 +342,8 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { // Not on domain boundary // 3rd derivative at right cell boundary - const BoutReal d3fdx3 = - (f(i, j + 2, k) - 3. * f(i, j + 1, k) + 3. * f(i, j, k) - f(i, j - 1, k)); + const BoutReal d3fdx3 = (f(i, j + 2, k) - (3. * f(i, j + 1, k)) + + (3. * f(i, j, k)) - f(i, j - 1, k)); result(i, j, k) += d3fdx3 * factor_rc; result(i, j + 1, k) -= d3fdx3 * factor_rp; @@ -363,10 +364,10 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { common_factor / (coord->J(i, j + 1, k) * coord->dy(i, j + 1, k)); const BoutReal d3fdx3 = - -((16. / 5) * 0.5 * (f(i, j + 1, k) + f(i, j, k)) // Boundary value f_b - - 6. * f(i, j, k) // f_0 - + 4. * f(i, j - 1, k) // f_1 - - (6. / 5) * f(i, j - 2, k) // f_2 + -(((16. / 5) * 0.5 * (f(i, j + 1, k) + f(i, j, k))) // Boundary value f_b + - (6. * f(i, j, k)) // f_0 + + (4. * f(i, j - 1, k)) // f_1 + - ((6. / 5) * f(i, j - 2, k)) // f_2 ); result(i, j, k) += d3fdx3 * factor_rc; @@ -392,8 +393,8 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { common_factor / (coord->J(i, j - 1, k) * coord->dy(i, j - 1, k)); // Not on a domain boundary - const BoutReal d3fdx3 = - (f(i, j + 1, k) - 3. * f(i, j, k) + 3. * f(i, j - 1, k) - f(i, j - 2, k)); + const BoutReal d3fdx3 = (f(i, j + 1, k) - (3. * f(i, j, k)) + + (3. * f(i, j - 1, k)) - f(i, j - 2, k)); result(i, j, k) -= d3fdx3 * factor_lc; result(i, j - 1, k) += d3fdx3 * factor_lm; @@ -410,10 +411,10 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { const BoutReal factor_lm = common_factor / (coord->J(i, j - 1, k) * coord->dy(i, j - 1, k)); const BoutReal d3fdx3 = - -(-(16. / 5) * 0.5 * (f(i, j - 1, k) + f(i, j, k)) // Boundary value f_b - + 6. * f(i, j, k) // f_0 - - 4. * f(i, j + 1, k) // f_1 - + (6. / 5) * f(i, j + 2, k) // f_2 + -((-(16. / 5) * 0.5 * (f(i, j - 1, k) + f(i, j, k))) // Boundary value f_b + + (6. * f(i, j, k)) // f_0 + - (4. * f(i, j + 1, k)) // f_1 + + ((6. / 5) * f(i, j + 2, k)) // f_2 ); result(i, j, k) -= d3fdx3 * factor_lc; @@ -436,8 +437,9 @@ void communicateFluxes(Field3D& f) { throw BoutException("communicateFluxes: Sorry!"); } - int size = mesh->LocalNy * mesh->LocalNz; - comm_handle xin, xout; + const int size = mesh->LocalNy * mesh->LocalNz; + comm_handle xin = nullptr; + comm_handle xout = nullptr; // Cache results to silence spurious compiler warning about xin, // xout possibly being uninitialised when used const bool not_first = mesh->periodicX || !mesh->firstX(); @@ -496,45 +498,45 @@ Field3D Div_Perp_Lap(const Field3D& a, const Field3D& f, CELL_LOC outloc) { // o --- gD --- o // Coordinates* coords = a.getCoordinates(outloc); - Mesh* mesh = f.getMesh(); + const Mesh* mesh = f.getMesh(); 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++) { // wrap k-index around as Z is (currently) periodic. - int kp = (k + 1) % (mesh->LocalNz); - int km = (k - 1 + mesh->LocalNz) % (mesh->LocalNz); + const int kp = (k + 1) % (mesh->LocalNz); + const int km = (k - 1 + mesh->LocalNz) % (mesh->LocalNz); // Calculate gradients on cell faces -- assumes constant grid spacing - BoutReal gR = - (coords->g11(i, j, k) + coords->g11(i + 1, j, k)) - * (f(i + 1, j, k) - f(i, j, k)) - / (coords->dx(i + 1, j, k) + coords->dx(i, j, k)) - + 0.5 * (coords->g13(i, j, k) + coords->g13(i + 1, j, k)) - * (f(i + 1, j, kp) - f(i + 1, j, km) + f(i, j, kp) - f(i, j, km)) - / (4. * coords->dz(i, j, k)); - - BoutReal gL = - (coords->g11(i - 1, j, k) + coords->g11(i, j, k)) - * (f(i, j, k) - f(i - 1, j, k)) - / (coords->dx(i - 1, j, k) + coords->dx(i, j, k)) - + 0.5 * (coords->g13(i - 1, j, k) + coords->g13(i, j, k)) - * (f(i - 1, j, kp) - f(i - 1, j, km) + f(i, j, kp) - f(i, j, km)) - / (4 * coords->dz(i, j, k)); - - BoutReal gD = - coords->g13(i, j, k) - * (f(i + 1, j, km) - f(i - 1, j, km) + f(i + 1, j, k) - f(i - 1, j, k)) - / (4. * coords->dx(i, j, k)) - + coords->g33(i, j, k) * (f(i, j, k) - f(i, j, km)) / coords->dz(i, j, k); - - BoutReal gU = - coords->g13(i, j, k) - * (f(i + 1, j, kp) - f(i - 1, j, kp) + f(i + 1, j, k) - f(i - 1, j, k)) - / (4. * coords->dx(i, j, k)) - + coords->g33(i, j, k) * (f(i, j, kp) - f(i, j, k)) / coords->dz(i, j, k); + const BoutReal gR = + ((coords->g11(i, j, k) + coords->g11(i + 1, j, k)) + * (f(i + 1, j, k) - f(i, j, k)) + / (coords->dx(i + 1, j, k) + coords->dx(i, j, k))) + + (0.5 * (coords->g13(i, j, k) + coords->g13(i + 1, j, k)) + * (f(i + 1, j, kp) - f(i + 1, j, km) + f(i, j, kp) - f(i, j, km)) + / (4. * coords->dz(i, j, k))); + + const BoutReal gL = + ((coords->g11(i - 1, j, k) + coords->g11(i, j, k)) + * (f(i, j, k) - f(i - 1, j, k)) + / (coords->dx(i - 1, j, k) + coords->dx(i, j, k))) + + (0.5 * (coords->g13(i - 1, j, k) + coords->g13(i, j, k)) + * (f(i - 1, j, kp) - f(i - 1, j, km) + f(i, j, kp) - f(i, j, km)) + / (4 * coords->dz(i, j, k))); + + const BoutReal gD = + (coords->g13(i, j, k) + * (f(i + 1, j, km) - f(i - 1, j, km) + f(i + 1, j, k) - f(i - 1, j, k)) + / (4. * coords->dx(i, j, k))) + + (coords->g33(i, j, k) * (f(i, j, k) - f(i, j, km)) / coords->dz(i, j, k)); + + const BoutReal gU = + (coords->g13(i, j, k) + * (f(i + 1, j, kp) - f(i - 1, j, kp) + f(i + 1, j, k) - f(i - 1, j, k)) + / (4. * coords->dx(i, j, k))) + + (coords->g33(i, j, k) * (f(i, j, kp) - f(i, j, k)) / coords->dz(i, j, k)); // Flow right BoutReal flux = gR * 0.25 * (coords->J(i + 1, j, k) + coords->J(i, j, k)) @@ -562,4 +564,71 @@ Field3D Div_Perp_Lap(const Field3D& a, const Field3D& f, CELL_LOC outloc) { return result; } +// Explicit instantiations of flux-limited finite volume methods +template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux); +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(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux); +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(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux); +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(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux); +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(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, + bool bndry_flux); +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(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, + bool fixflux = true); +template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, + bool bndry_flux); +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(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux); +template Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, Field3D& flow_ylow, + bool fixflux = true); + } // Namespace FV diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index afbd336bb2..706fea9c7a 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -120,8 +120,12 @@ XZHermiteSplineBase::XZHermiteSplineBase(int y_offset, Mesh if (options == nullptr) { options = &Options::root()["mesh:paralleltransform:xzinterpolation"]; } - abs_fac_monotonic = (*options)["abs_tol"].withDefault(abs_fac_monotonic); - rel_fac_monotonic = (*options)["rel_tol"].withDefault(rel_fac_monotonic); + abs_fac_monotonic = (*options)["atol"] + .doc("Absolute tolerance for clipping overshoot") + .withDefault(abs_fac_monotonic); + rel_fac_monotonic = (*options)["rtol"] + .doc("Relative tolerance for clipping overshoot") + .withDefault(rel_fac_monotonic); } // Index arrays contain guard cells in order to get subscripts right diff --git a/tests/MMS/CMakeLists.txt b/tests/MMS/CMakeLists.txt index a6667bcfa5..510385d54c 100644 --- a/tests/MMS/CMakeLists.txt +++ b/tests/MMS/CMakeLists.txt @@ -8,6 +8,7 @@ add_subdirectory(spatial/d2dx2) add_subdirectory(spatial/d2dz2) add_subdirectory(spatial/diffusion) add_subdirectory(spatial/fci) +add_subdirectory(spatial/finite-volume) add_subdirectory(time) add_subdirectory(time-petsc) add_subdirectory(wave-1d) diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp index 93e2101473..76ac3035c9 100644 --- a/tests/MMS/spatial/fci/data/BOUT.inp +++ b/tests/MMS/spatial/fci/data/BOUT.inp @@ -5,6 +5,9 @@ div_par_solution = (0.01*x + 0.045)*(-12.5663706143592*cos(y - 2*z) - 6.28318530 div_par_K_grad_par_solution = (0.01*x + 0.045)*(6.28318530717959*sin(y - z) - 0.628318530717959*sin(y - z)/(0.01*x + 0.045))*(6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z))/((0.01*x + 0.045)^2 + 1.0) + (6.28318530717959*(0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0))*cos(y - z)/sqrt((0.01*x + 0.045)^2 + 1.0) K = cos(y - z) laplace_par_solution = (0.01*x + 0.045)*(6.28318530717959*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/((0.01*x + 0.045)*sqrt((0.01*x + 0.045)^2 + 1.0)))/sqrt((0.01*x + 0.045)^2 + 1.0) +FV_div_par_mod_solution = (0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*((sin(y - 2*z) + sin(y - z))*sin(y - z)/(0.01*x + 0.045) + (-2*cos(y - 2*z) - cos(y - z))*cos(y - z)/(0.01*x + 0.045)) - 0.628318530717959*(sin(y - 2*z) + sin(y - z))*sin(y - z)/(0.01*x + 0.045) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))*cos(y - z)/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0) +FV_div_par_fvv_solution = (0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(2*(sin(y - 2*z) + sin(y - z))*sin(y - z)*cos(y - z)/(0.01*x + 0.045) + (-2*cos(y - 2*z) - cos(y - z))*cos(y - z)^2/(0.01*x + 0.045)) - 1.25663706143592*(sin(y - 2*z) + sin(y - z))*sin(y - z)*cos(y - z)/(0.01*x + 0.045) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))*cos(y - z)^2/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0) +div_par_K_grad_par_mod_solution = (0.01*x + 0.045)*(6.28318530717959*sin(y - z) - 0.628318530717959*sin(y - z)/(0.01*x + 0.045))*(6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z))/((0.01*x + 0.045)^2 + 1.0) + (6.28318530717959*(0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0))*cos(y - z)/sqrt((0.01*x + 0.045)^2 + 1.0) [mesh] symmetricglobalx = true diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx index 90bf107f3c..e75f857a5f 100644 --- a/tests/MMS/spatial/fci/fci_mms.cxx +++ b/tests/MMS/spatial/fci/fci_mms.cxx @@ -4,6 +4,7 @@ #include "bout/field.hxx" #include "bout/field3d.hxx" #include "bout/field_factory.hxx" +#include "bout/fv_ops.hxx" #include "bout/globals.hxx" #include "bout/options.hxx" #include "bout/options_io.hxx" @@ -52,12 +53,21 @@ int main(int argc, char** argv) { // Add mesh geometry variables mesh->outputVars(dump); + // Dummy variable for *_mod overloads + Field3D flow_ylow; + fci_op_test("grad_par", dump, input, Grad_par(input)); fci_op_test("grad2_par2", dump, input, Grad2_par2(input)); fci_op_test("div_par", dump, input, Div_par(input)); fci_op_test("div_par_K_grad_par", dump, input, Div_par_K_Grad_par(K, input)); + fci_op_test("div_par_K_grad_par_mod", dump, input, + Div_par_K_Grad_par_mod(K, input, flow_ylow)); fci_op_test("laplace_par", dump, input, Laplace_par(input)); + // Finite volume methods + fci_op_test("FV_div_par_mod", dump, input, FV::Div_par_mod(input, K, K, flow_ylow)); + fci_op_test("FV_div_par_fvv", dump, input, FV::Div_par_fvv(input, K, K)); + bout::writeDefaultOutputFile(dump); BoutFinalise(); diff --git a/tests/MMS/spatial/fci/mms.py b/tests/MMS/spatial/fci/mms.py index b28e337ac0..801a8d3f26 100755 --- a/tests/MMS/spatial/fci/mms.py +++ b/tests/MMS/spatial/fci/mms.py @@ -16,6 +16,7 @@ f = sin(y - z) + sin(y - 2 * z) K = cos(z - y) + Lx = 0.1 Ly = 10.0 Lz = 1.0 @@ -60,7 +61,10 @@ def FCI_Laplace_par(f: Expr) -> Expr: ("grad2_par2_solution", FCI_grad2_par2(f)), ("div_par_solution", FCI_div_par(f)), ("div_par_K_grad_par_solution", FCI_div_par_K_grad_par(f, K)), + ("div_par_K_grad_par_mod_solution", FCI_div_par_K_grad_par(f, K)), ("laplace_par_solution", FCI_Laplace_par(f)), + ("FV_div_par_mod_solution", FCI_div_par(f * K)), + ("FV_div_par_fvv_solution", FCI_div_par(f * K * K)), ): expr_str = exprToStr(expr) print(f"{name} = {expr_str}") diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 9dd4d1817f..ed22dd70ce 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -22,7 +22,21 @@ from numpy import array, linspace, log, polyfit DIRECTORY = "data" NPROC = 2 MTHREAD = 2 -OPERATORS = ("grad_par", "grad2_par2", "div_par", "div_par_K_grad_par", "laplace_par") +OPERATORS = ( + "grad_par", + "grad2_par2", + "div_par", + "div_par_K_grad_par", + "div_par_K_grad_par_mod", + "laplace_par", + "FV_div_par_mod", + "FV_div_par_fvv", +) + +# div_par_fvv is also tested in ../finite-volume, where 1.5th order is achieved +operator_order = { + "FV_div_par_fvv": 1, +} # Note that we need at least _2_ interior points for hermite spline # interpolation due to an awkwardness with the boundaries NX = 4 @@ -165,7 +179,10 @@ def check_fci_operators(name: str, case: dict) -> bool: for operator in OPERATORS: test_name = f"{operator} {name}" success = assert_convergence( - final_errors[operator]["l_2"], dx, test_name, order + final_errors[operator]["l_2"], + dx, + test_name, + operator_order.get(operator, order), ) if not success: failures.append(test_name) diff --git a/tests/MMS/spatial/finite-volume/CMakeLists.txt b/tests/MMS/spatial/finite-volume/CMakeLists.txt new file mode 100644 index 0000000000..7180eb7d98 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/CMakeLists.txt @@ -0,0 +1,6 @@ +bout_add_mms_test( + MMS-spatial-finite-volume + SOURCES fv_mms.cxx + USE_RUNTEST USE_DATA_BOUT_INP + PROCESSORS 2 +) diff --git a/tests/MMS/spatial/finite-volume/data/BOUT.inp b/tests/MMS/spatial/finite-volume/data/BOUT.inp new file mode 100644 index 0000000000..029011e437 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/data/BOUT.inp @@ -0,0 +1,23 @@ +input_field = 0.1*sin(2.0*y) + 1 +FV_Div_par_mod_solution = -0.188495559215388*sin(3.0*y) +FV_Div_par_fvv_solution = -0.376991118430775*(0.1*cos(3.0*y) + 1)*sin(3.0*y)/(0.1*sin(2.0*y) + 1) - 0.125663706143592*(0.1*cos(3.0*y) + 1)^2*cos(2.0*y)/(0.1*sin(2.0*y) + 1)^2 +FV_Div_par_solution = -0.188495559215388*sin(3.0*y) +FV_Div_par_K_Grad_par_solution = 0.125663706143592*(-0.188495559215388*sin(3.0*y)/(0.1*sin(2.0*y) + 1) - 0.125663706143592*(0.1*cos(3.0*y) + 1)*cos(2.0*y)/(0.1*sin(2.0*y) + 1)^2)*cos(2.0*y) - 0.15791367041743*(0.1*cos(3.0*y) + 1)*sin(2.0*y)/(0.1*sin(2.0*y) + 1) +FV_Div_par_K_Grad_par_mod_solution = 0.125663706143592*(-0.188495559215388*sin(3.0*y)/(0.1*sin(2.0*y) + 1) - 0.125663706143592*(0.1*cos(3.0*y) + 1)*cos(2.0*y)/(0.1*sin(2.0*y) + 1)^2)*cos(2.0*y) - 0.15791367041743*(0.1*cos(3.0*y) + 1)*sin(2.0*y)/(0.1*sin(2.0*y) + 1) +v = (0.1*cos(3.0*y) + 1)/(0.1*sin(2.0*y) + 1) + +[mesh] +MXG = 0 + +nx = 1 +ny = 128 +nz = 1 + +Ly = 10 + +dy = Ly / ny +J = 1 # Identity metric + +[mesh:ddy] +first = C2 +second = C2 diff --git a/tests/MMS/spatial/finite-volume/fv_mms.cxx b/tests/MMS/spatial/finite-volume/fv_mms.cxx new file mode 100644 index 0000000000..1de6052e89 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/fv_mms.cxx @@ -0,0 +1,91 @@ +#include "bout/bout.hxx" +#include "bout/difops.hxx" +#include "bout/field.hxx" +#include "bout/field3d.hxx" +#include "bout/field_factory.hxx" +#include "bout/fv_ops.hxx" +#include "bout/globals.hxx" +#include "bout/options.hxx" +#include "bout/options_io.hxx" + +#include + +#include +#include + +namespace { +auto fv_op_test(const std::string& name, Options& dump, const Field3D& input, + const Field3D& result, std::string suffix = "") { + auto* mesh = input.getMesh(); + const Field3D solution{FieldFactory::get()->create3D(fmt::format("{}_solution", name), + Options::getRoot(), mesh)}; + const Field3D error{result - solution}; + + dump[fmt::format("{}{}_l_2", name, suffix)] = + sqrt(mean(SQ(error), true, "RGN_NOBNDRY")); + dump[fmt::format("{}{}_l_inf", name, suffix)] = max(abs(error), true, "RGN_NOBNDRY"); + + dump[fmt::format("{}{}_result", name, suffix)] = result; + dump[fmt::format("{}{}_error", name, suffix)] = error; + dump[fmt::format("{}{}_input", name, suffix)] = input; + dump[fmt::format("{}{}_solution", name, suffix)] = solution; +} +} // namespace + +int main(int argc, char** argv) { + BoutInitialise(argc, argv); + + using bout::globals::mesh; + + Field3D input{FieldFactory::get()->create3D("input_field", Options::getRoot(), mesh)}; + Field3D v{FieldFactory::get()->create3D("v", Options::getRoot(), mesh)}; + + // Communicate to calculate parallel transform. + mesh->communicate(input, v); + + Options dump; + // Add mesh geometry variables + mesh->outputVars(dump); + dump["v"] = v; + + // Dummy variable for *_mod overloads + Field3D flow_ylow; + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), "_MC"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_MC"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), "_MC"); + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), "_Upwind"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_Upwind"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), + "_Upwind"); + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), "_Fromm"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_Fromm"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), + "_Fromm"); + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), "_MinMod"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_MinMod"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), + "_MinMod"); + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), + "_Superbee"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_Superbee"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), + "_Superbee"); + + fv_op_test("FV_Div_par_K_Grad_par", dump, input, FV::Div_par_K_Grad_par(v, input)); + fv_op_test("FV_Div_par_K_Grad_par_mod", dump, input, + Div_par_K_Grad_par_mod(v, input, flow_ylow)); + + bout::writeDefaultOutputFile(dump); + + BoutFinalise(); +} diff --git a/tests/MMS/spatial/finite-volume/makefile b/tests/MMS/spatial/finite-volume/makefile new file mode 100644 index 0000000000..7111f5ed8e --- /dev/null +++ b/tests/MMS/spatial/finite-volume/makefile @@ -0,0 +1,6 @@ + +BOUT_TOP = ../../../.. + +SOURCEC = fv_mms.cxx + +include $(BOUT_TOP)/make.config diff --git a/tests/MMS/spatial/finite-volume/mms.py b/tests/MMS/spatial/finite-volume/mms.py new file mode 100755 index 0000000000..a95ecc1328 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/mms.py @@ -0,0 +1,64 @@ +#!/usr/bin/env python3 +# +# Generate manufactured solution and sources for FCI test +# + +from math import pi +import warnings + +from boututils.boutwarnings import AlwaysWarning +from boutdata.data import BoutOptionsFile +from boutdata.mms import exprToStr, y, Grad_par, Div_par, Metric +from sympy import sin, cos, Expr + +warnings.simplefilter("ignore", AlwaysWarning) + +# Length of the y domain +Ly = 10.0 + +# Identity +metric = Metric() + +# Define solution in terms of input x,y,z +f = 1 + 0.1 * sin(2 * y) +fv = 1 + 0.1 * cos(3 * y) + + +# Turn solution into real x and z coordinates +replace = [(y, metric.y * 2 * pi / Ly)] + +f = f.subs(replace) +fv = fv.subs(replace) +v = fv / f + +# Substitute back to get input y coordinates +replace = [(metric.y, y * Ly / (2 * pi))] + + +def Grad2_par2(f: Expr) -> Expr: + return Grad_par(Grad_par(f)) + + +def Div_par_K_Grad_par(f: Expr, K: Expr) -> Expr: + return (K * Grad2_par2(f)) + (Div_par(K) * Grad_par(f)) + + +############################################ +# Equations solved + +options = BoutOptionsFile("data/BOUT.inp") + +for name, expr in ( + ("input_field", f), + ("v", v), + ("FV_Div_par_solution", Div_par(f * v)), + ("FV_Div_par_K_Grad_par_solution", Div_par_K_Grad_par(f, v)), + ("FV_Div_par_K_Grad_par_mod_solution", Div_par_K_Grad_par(f, v)), + ("FV_Div_par_mod_solution", Div_par(f * v)), + ("FV_Div_par_fvv_solution", Div_par(f * v * v)), +): + expr_str = exprToStr(expr.subs(replace)) + print(f"{name} = {expr_str}") + options[name] = expr_str + +options.write("data/BOUT.inp", overwrite=True) diff --git a/tests/MMS/spatial/finite-volume/runtest b/tests/MMS/spatial/finite-volume/runtest new file mode 100755 index 0000000000..bcd4672545 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/runtest @@ -0,0 +1,201 @@ +#!/usr/bin/env python3 +# +# Python script to run and analyse MMS test +# + +import argparse +import json +import pathlib +import sys +from time import time + +from boutdata.collect import collect +from boututils.run_wrapper import build_and_log, launch_safe +from numpy import array, log, polyfit + +# Global parameters +DIRECTORY = "data" +NPROC = 2 +MTHREAD = 2 +OPERATORS = { + # Slope-limiters necessarily reduce the accuracy in places + "FV_Div_par_MC": 1.5, + "FV_Div_par_mod_MC": 1.5, + "FV_Div_par_fvv_MC": 1.5, + "FV_Div_par_Upwind": 1, + "FV_Div_par_mod_Upwind": 1, + "FV_Div_par_fvv_Upwind": 1, + "FV_Div_par_Fromm": 1.5, + "FV_Div_par_mod_Fromm": 1.5, + "FV_Div_par_fvv_Fromm": 1.5, + "FV_Div_par_MinMod": 1.5, + "FV_Div_par_mod_MinMod": 1.5, + "FV_Div_par_fvv_MinMod": 1.5, + "FV_Div_par_Superbee": 1.5, + "FV_Div_par_mod_Superbee": 1.5, + "FV_Div_par_fvv_Superbee": 1.5, + "FV_Div_par_K_Grad_par": 2, + "FV_Div_par_K_Grad_par_mod": 2, +} +# Resolution in y and z +NLIST = [8, 16, 32, 64] +dx = 1.0 / array(NLIST) + + +def quiet_collect(name: str) -> float: + # Index to return a plain (numpy) float rather than `BoutArray` + return collect( + name, + tind=[1, 1], + info=False, + path=DIRECTORY, + xguards=False, + yguards=False, + )[()] + + +def assert_convergence(error, dx, name, expected) -> bool: + fit = polyfit(log(dx), log(error), 1) + order = fit[0] + print(f"{name} convergence order = {order:f} (fit)", end="") + + order = log(error[-2] / error[-1]) / log(dx[-2] / dx[-1]) + print(f", {order:f} (small spacing)", end="") + + # Should be close to the expected order + success = order > expected * 0.95 + print(f"\t............ {'PASS' if success else 'FAIL'}") + + return success + + +def run_fv_operators(nz: int) -> dict[str, float]: + # Command to run + cmd = f"./fv_mms MZ={nz} mesh:ny={nz}" + print(f"Running command: {cmd}", end="") + + # Launch using MPI + start = time() + status, out = launch_safe(cmd, nproc=NPROC, mthread=MTHREAD, pipe=True) + print(f" ... done in {time() - start:.3}s") + + # Save output to log file + pathlib.Path(f"run.log.{nz}").write_text(out) + + if status: + print(f"Run failed!\nOutput was:\n{out}") + sys.exit(status) + + return { + operator: { + "l_2": quiet_collect(f"{operator}_l_2"), + "l_inf": quiet_collect(f"{operator}_l_inf"), + } + for operator in OPERATORS + } + + +def transpose( + errors: list[dict[str, dict[str, float]]], +) -> dict[str, dict[str, list[float]]]: + """Turn a list of {operator: error} into a dict of {operator: [errors]}""" + + kinds = ("l_2", "l_inf") + result = {operator: {kind: [] for kind in kinds} for operator in OPERATORS} + for error in errors: + for k, v in error.items(): + for kind in kinds: + result[k][kind].append(v[kind]) + return result + + +def test_fv_operators() -> bool: + failures = [] + + all_errors = [] + + for n in NLIST: + errors = run_fv_operators(n) + all_errors.append(errors) + + for operator in OPERATORS: + l_2 = errors[operator]["l_2"] + l_inf = errors[operator]["l_inf"] + + print(f"{operator} errors: l-2 {l_2:f} l-inf {l_inf:f}") + + final_errors = transpose(all_errors) + for operator, order in OPERATORS.items(): + success = assert_convergence(final_errors[operator]["l_2"], dx, operator, order) + if not success: + failures.append(operator) + + return final_errors, failures + + +def make_plots(cases: dict[str, dict]): + try: + import matplotlib.pyplot as plt + except ImportError: + print("No matplotlib") + return + + num_operators = len(OPERATORS) + fig, axes = plt.subplots(1, num_operators, figsize=(num_operators * 4, 4)) + + for ax, operator in zip(axes, OPERATORS): + ax.loglog(dx, cases[operator]["l_2"], "-", label="$l_2$") + ax.loglog(dx, cases[operator]["l_inf"], "--", label="$l_\\inf$") + ax.legend(loc="upper left") + ax.grid() + ax.set_title(f"Error scaling for {operator}") + ax.set_xlabel(r"Mesh spacing $\delta x$") + ax.set_ylabel("Error norm") + + fig.tight_layout() + fig.savefig("fv_mms.pdf") + print("Plot saved to fv_mms.pdf") + + if args.show_plots: + plt.show() + plt.close() + + +if __name__ == "__main__": + build_and_log("Finite volume MMS test") + + parser = argparse.ArgumentParser("Error scaling test for finite volume operators") + parser.add_argument( + "--make-plots", action="store_true", help="Create plots of error scaling" + ) + parser.add_argument( + "--show-plots", + action="store_true", + help="Stop and show plots, implies --make-plots", + ) + parser.add_argument( + "--dump-errors", + type=str, + help="Output file to dump errors as JSON", + default="fv_operator_errors.json", + ) + + args = parser.parse_args() + + error2, failures = test_fv_operators() + success = len(failures) == 0 + + if args.dump_errors: + pathlib.Path(args.dump_errors).write_text(json.dumps(error2)) + + if args.make_plots or args.show_plots: + make_plots(error2) + + if success: + print("\nAll tests passed") + else: + print("\nSome tests failed:") + for failure in failures: + print(f"\t{failure}") + + sys.exit(0 if success else 1) diff --git a/tests/unit/include/bout/test_fv_ops.cxx b/tests/unit/include/bout/test_fv_ops.cxx index e78c0cca41..5064b076b9 100644 --- a/tests/unit/include/bout/test_fv_ops.cxx +++ b/tests/unit/include/bout/test_fv_ops.cxx @@ -1,6 +1,6 @@ #include "gtest/gtest.h" -#include +#include TEST(FVOpsLimiterTest, VanAlbadaConstant) { FV::Stencil1D s{}; diff --git a/tests/unit/mesh/test_coordinates.cxx b/tests/unit/mesh/test_coordinates.cxx index ea233cd9ca..8e1aa2e13f 100644 --- a/tests/unit/mesh/test_coordinates.cxx +++ b/tests/unit/mesh/test_coordinates.cxx @@ -374,3 +374,30 @@ TEST_F(CoordinatesTest, CellAreasUpdate) { EXPECT_TRUE(IsFieldEqual(coords.cell_volume(), 8.0)); } + +TEST_F(CoordinatesTest, CellAreaZComputedAfterCellVolume) { + Coordinates coords{mesh, + FieldMetric{1.0}, // dx + FieldMetric{1.0}, // dy + FieldMetric{1.0}, // dz + FieldMetric{6.0}, // J + FieldMetric{1.0}, // Bxy + FieldMetric{1.0}, // g11 + FieldMetric{1.0}, // g22 + FieldMetric{9.0}, // g33 + FieldMetric{0.0}, // g12 + FieldMetric{0.0}, // g13 + FieldMetric{0.0}, // g23 + FieldMetric{4.0}, // g_11 + FieldMetric{1.0}, // g_22 + FieldMetric{9.0}, // g_33 + FieldMetric{0.0}, // g_12 + FieldMetric{0.0}, // g_13 + FieldMetric{0.0}, // g_23 + FieldMetric{0.0}, // ShiftTorsion + FieldMetric{0.0}}; // IntShiftTorsion + + EXPECT_TRUE(IsFieldEqual(coords.cell_volume(), 6.0)); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_zlow(), 2.0, "RGN_NOZ")); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_zhigh(), 2.0, "RGN_NOZ")); +}