diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 789fdb2a1..b0b14d33d 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -14,7 +14,7 @@ defaults: jobs: cmake_tests: - name: ${{ matrix.build_type.test_type }} tests, ${{ matrix.config.name }} + name: ${{ matrix.build_type.test_type }} tests, ${{ matrix.build_type.name }} runs-on: ${{ matrix.config.os }} env: OMP_NUM_THREADS: ${{ matrix.config.omp_num_threads }} @@ -22,29 +22,36 @@ jobs: OMPI_MCA_rmaps_base_oversubscribe: yes MPIRUN: mpiexec -np LD_LIBRARY_PATH: /home/runner/local/lib:$LD_LIBRARY_PATH - BOUT_CONFIGURE_OPTIONS: ${{ matrix.config.configure_options }} -DHERMES_COVERAGE=${{ matrix.build_type.coverage }} + BOUT_CONFIGURE_OPTIONS: ${{ matrix.build_type.configure_options }} -DHERMES_COVERAGE=${{ matrix.build_type.coverage }} BUILD_TYPE: ${{ matrix.build_type.cmake_build_type }} TEST_TYPE: ${{ matrix.build_type.test_type }} strategy: fail-fast: false matrix: config: - - name: "CMake default options" - os: ubuntu-latest - configure_options: "-DPACKAGE_TESTS=OFF -DHERMES_ERROR_ON_WARNINGS=ON" + - os: ubuntu-latest build_type: - cmake_build_type: Debug test_type: unit coverage: ON + name: "CMake default options" + configure_options: "-DPACKAGE_TESTS=OFF -DHERMES_ERROR_ON_WARNINGS=ON" - cmake_build_type: Release test_type: integration coverage: OFF + name: "CMake default options" + configure_options: "-DPACKAGE_TESTS=OFF -DHERMES_ERROR_ON_WARNINGS=ON" + - cmake_build_type: Debug + test_type: None + name: "3D metrics (compile only)" + configure_options: "-DPACKAGE_TESTS=OFF -DHERMES_ERROR_ON_WARNINGS=ON -DBOUT_ENABLE_METRIC_3D=ON" + coverage: OFF steps: - name: Job information run: | echo Build: ${{ matrix.config.name }}, ${{ matrix.config.os }} - echo Configure options: ${{ matrix.config.configure_options }} + echo Configure options: ${{ matrix.build_type.configure_options }} echo Build type: ${{ matrix.build_type.cmake_build_type }} echo Coverage: ${{ matrix.build_type.coverage }} echo Tests to run: ${{ matrix.build_type.test_type }} diff --git a/CMakeLists.txt b/CMakeLists.txt index 8ddb7466c..8a9f8e12b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -430,14 +430,6 @@ set(HERMES_SLOPE_LIMITER set_property(CACHE HERMES_SLOPE_LIMITER PROPERTY STRINGS ${SLOPE_LIMITERS}) message(STATUS "Slope limiter: ${HERMES_SLOPE_LIMITER}") -set(CONDUCTION_METHODS Original ProductJK Harmonic) -set(HERMES_CONDUCTION_METHOD - Original - CACHE STRING "Parallel conduction method") -set_property(CACHE HERMES_CONDUCTION_METHOD PROPERTY STRINGS - ${CONDUCTION_METHODS}) -message(STATUS "Conduction method: ${HERMES_CONDUCTION_METHOD}") - # Generate the build config header if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/include/hermes_build_config.hxx") diff --git a/docs/sphinx/solver_numerics.rst b/docs/sphinx/solver_numerics.rst index aba550570..cd75df4cd 100644 --- a/docs/sphinx/solver_numerics.rst +++ b/docs/sphinx/solver_numerics.rst @@ -230,21 +230,22 @@ The divergence of the flux, and so the rate of change of :math:`f` in cell Selecting slope limiter and conduction method ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -Hermes-3 selects both the parallel advection slope limiter and the -``Div_par_K_Grad_par_mod`` conduction discretisation at compile time +Hermes-3 selects the parallel advection slope limiter at compile time through CMake cache options: .. code-block:: bash cmake -S . -B build \ - -DHERMES_SLOPE_LIMITER=MC \ - -DHERMES_CONDUCTION_METHOD=Original + -DHERMES_SLOPE_LIMITER=MC ``HERMES_SLOPE_LIMITER`` Available values are ``MC``, ``VanAlbada``, ``WENO3``, ``MinMod``, ``Upwind``, and ``Superbee``. -``HERMES_CONDUCTION_METHOD`` + +The +``Div_par_K_Grad_par_mod`` conduction discretisation can be choosen at run time with: +``conduction_method`` Selects the parallel heat conduction discretisation used in :cpp:func:`Div_par_K_Grad_par_mod`. Available values are: diff --git a/external/BOUT-dev b/external/BOUT-dev index 2384555b6..a511a6632 160000 --- a/external/BOUT-dev +++ b/external/BOUT-dev @@ -1 +1 @@ -Subproject commit 2384555b6ad9763cddc3e99a99ed4c4c1037faa2 +Subproject commit a511a6632640ebfe2df2f1fd483dcd72c6702319 diff --git a/hermes-3.cxx b/hermes-3.cxx index 74fd0343b..eb0dd72ff 100644 --- a/hermes-3.cxx +++ b/hermes-3.cxx @@ -121,9 +121,9 @@ class DecayLengthBoundary : public BoundaryOp { // Get cell radial length Coordinates* coord = mesh->getCoordinates(); - Field2D dx = coord->dx; - Field2D g11 = coord->g11; - Field2D dr = + auto dx = coord->dx; + auto g11 = coord->g11; + Coordinates::FieldMetric dr = dx / sqrt(g11); // cell radial length. dr = dx/(Bpol * R) and g11 = (Bpol*R)**2 // Only implemented for cell centre quantities @@ -149,9 +149,10 @@ class DecayLengthBoundary : public BoundaryOp { // (0, -1) Y lower boundary (inner lower target) // Distance between final cell centre and inner guard cell centre in normalised units - BoutReal distance = - 0.5 - * (dr(bndry->x, bndry->y) + dr(bndry->x - bndry->bx, bndry->y - bndry->by)); + + BoutReal distance = 0.5 + * (dr(bndry->x, bndry->y, zk) + + dr(bndry->x - bndry->bx, bndry->y - bndry->by, zk)); // Exponential decay f(bndry->x, bndry->y, zk) = f(bndry->x - bndry->bx, bndry->y - bndry->by, zk) @@ -186,9 +187,6 @@ int Hermes::init(bool restarting) { output.write("Slope limiter: {}\n", hermes::limiter_typename); options["slope_limiter"] = hermes::limiter_typename; options["slope_limiter"].setConditionallyUsed(); - output.write("Conduction method: {}\n", hermes::conduction_typename); - options["conduction_method"] = hermes::conduction_typename; - options["conduction_method"].setConditionallyUsed(); // Choose normalisations Tnorm = options["Tnorm"].doc("Reference temperature [eV]").withDefault(100.); diff --git a/include/braginskii_conduction.hxx b/include/braginskii_conduction.hxx index 3a2efbcd7..6356a9c8b 100644 --- a/include/braginskii_conduction.hxx +++ b/include/braginskii_conduction.hxx @@ -78,6 +78,8 @@ private: /// - energy_flow_ylow Energy flow diagnostics. /// void transform_impl(GuardedOptions& state) override; + + std::string conduction_method{"orginal"}; }; namespace { diff --git a/include/braginskii_ion_viscosity.hxx b/include/braginskii_ion_viscosity.hxx index 0aacf49bb..9f5b5fee7 100644 --- a/include/braginskii_ion_viscosity.hxx +++ b/include/braginskii_ion_viscosity.hxx @@ -8,7 +8,7 @@ #include #include -#include +#include #include "component.hxx" @@ -56,7 +56,7 @@ private: std::string viscosity_collisions_mode; ///< Collision selection, either multispecies or ///< braginskii Field3D nu; ///< Collision frequency for conduction - Vector2D Curlb_B; ///< Curvature vector Curl(b/B) + VectorMetric Curlb_B; ///< Curvature vector Curl(b/B) bool bounce_frequency; ///< Modify the collision time with the bounce frequency? BoutReal bounce_frequency_q95; ///< Input q95 for when including bounce frequency change BoutReal bounce_frequency_epsilon; ///< Input inverse aspect ratio for including bounce diff --git a/include/classical_diffusion.hxx b/include/classical_diffusion.hxx index f73ce8731..4a4ced497 100644 --- a/include/classical_diffusion.hxx +++ b/include/classical_diffusion.hxx @@ -2,6 +2,8 @@ #ifndef CLASSICAL_DIFFUSION_H #define CLASSICAL_DIFFUSION_H +#include + #include "component.hxx" struct ClassicalDiffusion : public NamedComponent { @@ -12,7 +14,7 @@ struct ClassicalDiffusion : public NamedComponent { static constexpr auto type = "classical_diffusion"; private: - Field2D Bsq; // Magnetic field squared + Coordinates::FieldMetric Bsq; // Magnetic field squared bool diagnose; ///< Output additional diagnostics? Field3D Dn; ///< Particle diffusion coefficient diff --git a/include/component.hxx b/include/component.hxx index fa4c32ed7..fefa4c0e4 100644 --- a/include/component.hxx +++ b/include/component.hxx @@ -236,6 +236,7 @@ struct RegisterComponent : public ComponentFactory::RegisterInFactory(std::string(DerivedType::type)) { } }; +using RegisterUnavailableComponent = ComponentFactory::RegisterUnavailableInFactory; /// Faster non-printing getter for Options /// If this fails, it will throw BoutException diff --git a/include/diamagnetic_drift.hxx b/include/diamagnetic_drift.hxx index d91f8fca5..8487cdf96 100644 --- a/include/diamagnetic_drift.hxx +++ b/include/diamagnetic_drift.hxx @@ -2,6 +2,8 @@ #ifndef DIAMAGNETIC_DRIFT_H #define DIAMAGNETIC_DRIFT_H +#include + #include "component.hxx" /// Calculate diamagnetic flows @@ -12,9 +14,9 @@ struct DiamagneticDrift : public NamedComponent { static constexpr auto type = "diamagnetic_drift"; private: - Vector2D Curlb_B; + VectorMetric Curlb_B; bool bndry_flux; - Field2D diamag_form; + Coordinates::FieldMetric diamag_form; /// For every species, if it has: /// - temperature diff --git a/include/div_ops.hxx b/include/div_ops.hxx index 4c929c312..d02f054c3 100644 --- a/include/div_ops.hxx +++ b/include/div_ops.hxx @@ -26,10 +26,20 @@ #ifndef DIV_OPS_H #define DIV_OPS_H +#include +#include +#include +#include +#include #include -#include +#include +#include +#include +#include #include +#include + /*! * Diffusion in index space * @@ -38,43 +48,38 @@ * @param[in] The field to be differentiated * @param[in] bndry_flux Are fluxes through the boundary calculated? */ -const Field3D Div_par_diffusion_index(const Field3D& f, bool bndry_flux = true); +Field3D Div_par_diffusion_index(const Field3D& f, bool bndry_flux = true); -const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, - bool bndry_flux = true, bool poloidal = false, - bool positive = false); +Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bndry_flux = true, + bool poloidal = false, bool positive = false); /// This version has an extra coefficient 'g' that is linearly interpolated /// onto cell faces -const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, const Field3D& f, - bool bndry_flux = true); +Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, const Field3D& f, + bool bndry_flux = true); -const Field3D Div_Perp_Lap_FV_Index(const Field3D& a, const Field3D& f); +Field3D Div_Perp_Lap_FV_Index(const Field3D& a, const Field3D& f); -const Field3D Div_Z_FV_Index(const Field3D& a, const Field3D& f); +Field3D Div_Z_FV_Index(const Field3D& a, const Field3D& f); // 4th-order flux conserving term, in index space -const Field3D D4DX4_FV_Index(const Field3D& f, bool bndry_flux = false); -const Field3D D4DZ4_Index(const Field3D& f); +Field3D D4DX4_FV_Index(const Field3D& f, bool bndry_flux = false); +Field3D D4DZ4_Index(const Field3D& f); // Div ( k * Grad(f) ) -const Field2D Laplace_FV(const Field2D& k, const Field2D& f); +Field2D Laplace_FV(const Field2D& k, const Field2D& f); /// Perpendicular diffusion including X and Y directions /// Takes Div_a_Grad_perp from BOUT++ and adds flows -const Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, - Field3D& flux_xlow, Field3D& flux_ylow); +Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, Field3D& flux_xlow, + Field3D& flux_ylow); /// Same but with upwinding /// WARNING: Causes checkerboarding in neutral_mixed integrated test -const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f); +Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f); /// Same but with upwinding and flows /// WARNING: Causes checkerboarding in neutral_mixed integrated test -const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, - Field3D& flux_xlow, Field3D& flux_ylow); - -/// Version with energy flow diagnostic -Field3D Div_par_K_Grad_par_mod(const Field3D& k, const Field3D& f, Field3D& flow_ylow, - bool bndry_flux = true); +Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, + Field3D& flux_xlow, Field3D& flux_ylow); /*! * Div ( a Grad_perp(f) ) -- ∇⊥ ( a ⋅ ∇⊥ f) -- Vorticity @@ -86,853 +91,4 @@ Field3D Div_par_K_Grad_par_mod(const Field3D& k, const Field3D& f, Field3D& flow Field3D Div_a_Grad_perp_nonorthog(const Field3D& a, const Field3D& x, Field3D& flux_xlow, Field3D& flux_ylow); -namespace FV { - -/// Superbee limiter -/// -/// This corresponds to the limiter function -/// φ(r) = max(0, min(2r, 1), min(r,2) -/// -/// The value at cell right (i.e. i + 1/2) is: -/// -/// n.R = n.c - φ(r) (n.c - (n.p + n.c)/2) -/// = n.c + φ(r) (n.p - n.c)/2 -/// -/// Four regimes: -/// a) r < 1/2 -> φ(r) = 2r -/// n.R = n.c + gL -/// b) 1/2 < r < 1 -> φ(r) = 1 -/// n.R = n.c + gR/2 -/// c) 1 < r < 2 -> φ(r) = r -/// n.R = n.c + gL/2 -/// d) 2 < r -> φ(r) = 2 -/// n.R = n.c + gR -/// -/// where the left and right gradients are: -/// gL = n.c - n.m -/// gR = n.p - n.c -/// -struct Superbee { - void operator()(Stencil1D& n) { - BoutReal gL = n.c - n.L; - BoutReal gR = n.R - n.c; - - // r = gL / gR - // Limiter is φ(r) - if (gL * gR < 0) { - // Different signs => Zero gradient - n.L = n.R = n.c; - } else { - BoutReal sign = SIGN(gL); - gL = fabs(gL); - gR = fabs(gR); - BoutReal half_slope = sign * BOUTMAX(BOUTMIN(gL, 0.5 * gR), BOUTMIN(gR, 0.5 * gL)); - n.L = n.c - half_slope; - n.R = n.c + half_slope; - } - } -}; - -template -const Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, - const Field3D& wave_speed_in, bool fixflux = true) { - - ASSERT1(areFieldsCompatible(f_in, v_in)); - ASSERT1(areFieldsCompatible(f_in, wave_speed_in)); - - Mesh* mesh = f_in.getMesh(); - - CellEdges cellboundary; - - /// 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"); - - 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 - - // For right cell boundaries - BoutReal common_factor = (coord->J(i, j) + coord->J(i, j + 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); - - BoutReal flux_factor_rc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - BoutReal flux_factor_rp = - common_factor / (coord->dy(i, j + 1) * coord->J(i, j + 1)); - - // For left cell boundaries - common_factor = (coord->J(i, j) + coord->J(i, j - 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); - - BoutReal flux_factor_lc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - BoutReal flux_factor_lm = - common_factor / (coord->dy(i, j - 1) * coord->J(i, j - 1)); - - for (int k = 0; k < mesh->LocalNz; k++) { - - //////////////////////////////////////////// - // Reconstruct f at the cell faces - // This calculates s.R and s.L for the Right and Left - // face values on this cell - - // Reconstruct f at the cell faces - Stencil1D s; - s.c = f(i, j, k); - s.m = f(i, j - 1, k); - s.p = f(i, j + 1, k); - - cellboundary(s); // Calculate s.R and s.L - - // Reconstruct v at the cell faces - Stencil1D sv; - sv.c = v(i, j, k); - sv.m = v(i, j - 1, k); - sv.p = v(i, j + 1, k); - - cellboundary(sv); - - //////////////////////////////////////////// - // Right boundary - - // Calculate velocity at right boundary (y+1/2) - BoutReal v_mid = 0.5 * (sv.c + sv.p); - // And mid-point density at right boundary - BoutReal n_mid = 0.5 * (s.c + s.p); - BoutReal flux; - - if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) { - // Last point in domain - - if (fixflux) { - // Use mid-point to be consistent with boundary conditions - flux = n_mid * v_mid * v_mid; - } else { - // Add flux due to difference in boundary values - flux = s.R * sv.R * sv.R // Use right cell edge values - + BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.p)) * n_mid - * (sv.R - v_mid); // Damp differences in velocity, not flux - } - } else { - // Maximum wave speed in the two cells - BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k), - fabs(sv.c), fabs(sv.p)); - - flux = s.R * 0.5 * (sv.R + amax) * sv.R; - } - - result(i, j, k) += flux * flux_factor_rc; - result(i, j + 1, k) -= flux * flux_factor_rp; - - //////////////////////////////////////////// - // Calculate at left boundary - - v_mid = 0.5 * (sv.c + sv.m); - n_mid = 0.5 * (s.c + s.m); - - if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) { - // First point in domain - if (fixflux) { - // Use mid-point to be consistent with boundary conditions - flux = n_mid * v_mid * v_mid; - } else { - // Add flux due to difference in boundary values - flux = s.L * sv.L * sv.L - - BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.m)) * n_mid - * (sv.L - v_mid); - } - } else { - // Maximum wave speed in the two cells - BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k), - fabs(sv.c), fabs(sv.m)); - - flux = s.L * 0.5 * (sv.L - amax) * sv.L; - } - - result(i, j, k) -= flux * flux_factor_lc; - result(i, j - 1, k) += flux * flux_factor_lm; - } - } - } - return fromFieldAligned(result, "RGN_NOBNDRY"); -} - -// Calculates viscous heating due to numerical momentum fluxes -// and flow of kinetic energy (in flow_ylow) -template -const Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in, - const Field3D& wave_speed_in, Field3D& flow_ylow, - bool fixflux = true) { - - ASSERT1(areFieldsCompatible(f_in, v_in)); - ASSERT1(areFieldsCompatible(f_in, wave_speed_in)); - - Mesh* mesh = f_in.getMesh(); - - CellEdges cellboundary; - - /// 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"); - - Coordinates* coord = f_in.getCoordinates(); - - Field3D result{zeroFrom(f)}; - flow_ylow = zeroFrom(f); - - // Only need one guard cell, so no need to communicate fluxes - // Instead calculate in guard cells to preserve fluxes - int ys = mesh->ystart - 1; - int ye = mesh->yend + 1; - - for (int i = mesh->xstart; i <= mesh->xend; i++) { - - if (!mesh->firstY(i) || mesh->periodicY(i)) { - // Calculate in guard cell to get fluxes consistent between processors - ys = mesh->ystart - 1; - } else { - // Don't include the boundary cell. Note that this implies special - // handling of boundaries later - ys = mesh->ystart; - } - - if (!mesh->lastY(i) || mesh->periodicY(i)) { - // Calculate in guard cells - ye = mesh->yend + 1; - } else { - // Not in boundary cells - ye = mesh->yend; - } - - for (int j = ys; j <= ye; j++) { - // Pre-calculate factors which multiply fluxes - - // For right cell boundaries - BoutReal common_factor = (coord->J(i, j) + coord->J(i, j + 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); - - BoutReal flux_factor_rc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - BoutReal area_rp = common_factor * coord->dx(i, j + 1) * coord->dz(i, j + 1); - - // For left cell boundaries - common_factor = (coord->J(i, j) + coord->J(i, j - 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); - - BoutReal flux_factor_lc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - BoutReal area_lc = common_factor * coord->dx(i, j) * coord->dz(i, j); - - for (int k = 0; k < mesh->LocalNz; k++) { - - //////////////////////////////////////////// - // Reconstruct f at the cell faces - // This calculates s.R and s.L for the Right and Left - // face values on this cell - - // Reconstruct f at the cell faces - Stencil1D s; - s.c = f(i, j, k); - s.m = f(i, j - 1, k); - s.p = f(i, j + 1, k); - - cellboundary(s); // Calculate s.R and s.L - - // Reconstruct v at the cell faces - Stencil1D sv; - sv.c = v(i, j, k); - sv.m = v(i, j - 1, k); - sv.p = v(i, j + 1, k); - - cellboundary(sv); - - //////////////////////////////////////////// - // Right boundary - - // Calculate velocity at right boundary (y+1/2) - BoutReal v_mid = 0.5 * (sv.c + sv.p); - // And mid-point density at right boundary - BoutReal n_mid = 0.5 * (s.c + s.p); - - if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) { - // Last point in domain - - // Expected loss of kinetic energy into boundary - // This is used in the sheath boundary condition to calculate - // energy losses. - BoutReal expected_ke = 0.5 * n_mid * v_mid * v_mid * v_mid; - - BoutReal flux_mom; - if (fixflux) { - // Mid-point consistent with boundary conditions - // but kinetic energy loss will not match expected - // -> Adjust energy balance in pressure equation - flux_mom = n_mid * v_mid * v_mid; - } else { - flux_mom = s.R * sv.R * sv.R - + BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.p)) - * (s.R * sv.R - n_mid * v_mid); - } - - // Assume that particle flux is fixed to boundary value - const BoutReal flux_part = n_mid * v_mid; - - // d/dt(1/2 m n v^2) = v * d/dt(mnv) - 1/2 m v^2 * dn/dt - BoutReal actual_ke = sv.c * flux_mom - 0.5 * sv.c * sv.c * flux_part; - - // Note: If the actual loss was higher than expected, then - // plasma heating is needed to compensate - result(i, j, k) += (actual_ke - expected_ke) * flux_factor_rc; - - // Final flow through boundary is the expected value - flow_ylow(i, j + 1, k) += expected_ke * area_rp; //expected_ke * area_rp; - - } else { - // Maximum wave speed in the two cells - BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k), - fabs(sv.c), fabs(sv.p)); - - // Viscous heating due to relaxation of velocity towards midpoint - result(i, j, k) += - (amax + 0.5 * sv.R) * s.R * (sv.c - sv.p) * (sv.R - v_mid) * flux_factor_rc; - - // Kinetic energy flow into next cell. - // Note: Different from flow out of this cell; the difference - // is in the viscous heating. - BoutReal flux_part = s.R * 0.5 * (sv.R + amax); - BoutReal flux_mom = flux_part * sv.R; - - flow_ylow(i, j + 1, k) += - (sv.p * flux_mom - 0.5 * SQ(sv.p) * flux_part) * area_rp; - } - - //////////////////////////////////////////// - // Calculate at left boundary - - v_mid = 0.5 * (sv.c + sv.m); - n_mid = 0.5 * (s.c + s.m); - - // Expected KE loss. Note minus sign because negative v into boundary - BoutReal expected_ke = -0.5 * n_mid * v_mid * v_mid * v_mid; - - if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) { - // First point in domain - BoutReal flux_mom; - if (fixflux) { - // Use mid-point to be consistent with boundary conditions - flux_mom = n_mid * v_mid * v_mid; - } else { - // Add flux due to difference in boundary values - flux_mom = s.L * sv.L * sv.L - - BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.m)) - * (s.L * sv.L - n_mid * v_mid); - } - - // Assume that density flux is fixed to boundary value - const BoutReal flux_part = n_mid * v_mid; - - // d/dt(1/2 m n v^2) = v * d/dt(mnv) - 1/2 m v^2 * dn/dt - BoutReal actual_ke = -sv.c * flux_mom + 0.5 * sv.c * sv.c * flux_part; - - result(i, j, k) += (actual_ke - expected_ke) * flux_factor_lc; - - flow_ylow(i, j, k) -= expected_ke * area_lc; - } else { - // Maximum wave speed in the two cells - BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k), - fabs(sv.c), fabs(sv.m)); - - // Viscous heating due to relaxation - result(i, j, k) += - (amax - 0.5 * sv.L) * s.L * (sv.c - sv.m) * (sv.L - v_mid) * flux_factor_lc; - - // Kinetic energy flow into this cell. - // Note: Different from flow out of left cell; the difference - // is in the viscous heating. - BoutReal flux_part = s.L * 0.5 * (sv.L - amax); - BoutReal flux_mom = flux_part * sv.L; - - flow_ylow(i, j, k) += (sv.c * flux_mom - 0.5 * SQ(sv.c) * flux_part) * area_lc; - } - } - } - } - flow_ylow = fromFieldAligned(flow_ylow, "RGN_NOBNDRY"); - return fromFieldAligned(result, "RGN_NOBNDRY"); -} - -/// Finite volume parallel divergence -/// -/// NOTE: Modified version, applies limiter to velocity and field -/// Performs better (smaller overshoots) than Div_par -/// -/// Preserves the sum of f*J*dx*dy*dz over the domain -/// -/// @param[in] f_in The field being advected. -/// This will be reconstructed at cell faces -/// using the given CellEdges method -/// @param[in] v_in The advection velocity. -/// This will be interpolated to cell boundaries -/// using linear interpolation -/// @param[in] wave_speed_in Local maximum speed of all waves in the system at each -// point in space -/// @param[in] fixflux Fix the flux at the boundary to be the value at the -/// midpoint (for boundary conditions) -/// -/// @param[out] flow_ylow Flow at the lower Y cell boundary -/// Already includes area factor * flux -/// -/// NB: Uses to/from FieldAligned coordinates -template -const Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, - const Field3D& wave_speed_in, Field3D& flow_ylow, - bool fixflux = true) { - - ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); - ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in); - - Mesh* mesh = f_in.getMesh(); - - CellEdges cellboundary; - - ASSERT2(f_in.getDirectionY() == v_in.getDirectionY()); - ASSERT2(f_in.getDirectionY() == wave_speed_in.getDirectionY()); - const bool are_unaligned = - ((f_in.getDirectionY() == YDirectionType::Standard) - and (v_in.getDirectionY() == YDirectionType::Standard) - and (wave_speed_in.getDirectionY() == YDirectionType::Standard)); - - Field3D f = are_unaligned ? toFieldAligned(f_in, "RGN_NOX") : f_in; - Field3D v = are_unaligned ? toFieldAligned(v_in, "RGN_NOX") : v_in; - Field3D wave_speed = - are_unaligned ? toFieldAligned(wave_speed_in, "RGN_NOX") : wave_speed_in; - - Coordinates* coord = f_in.getCoordinates(); - - Field3D result{zeroFrom(f)}; - flow_ylow = zeroFrom(f); - - // Only need one guard cell, so no need to communicate fluxes - // Instead calculate in guard cells to preserve fluxes - int ys = mesh->ystart - 1; - int ye = mesh->yend + 1; - - for (int i = mesh->xstart; i <= mesh->xend; i++) { - - if (!mesh->firstY(i) || mesh->periodicY(i)) { - // Calculate in guard cell to get fluxes consistent between processors - ys = mesh->ystart - 1; - } else { - // Don't include the boundary cell. Note that this implies special - // handling of boundaries later - ys = mesh->ystart; - } - - if (!mesh->lastY(i) || mesh->periodicY(i)) { - // Calculate in guard cells - ye = mesh->yend + 1; - } else { - // Not in boundary cells - ye = mesh->yend; - } - - for (int j = ys; j <= ye; j++) { - // Pre-calculate factors which multiply fluxes -#if not(BOUT_USE_METRIC_3D) - // For right cell boundaries - BoutReal common_factor = (coord->J(i, j) + coord->J(i, j + 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); - - BoutReal flux_factor_rc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - BoutReal flux_factor_rp = - common_factor / (coord->dy(i, j + 1) * coord->J(i, j + 1)); - - BoutReal area_rp = common_factor * coord->dx(i, j + 1) * coord->dz(i, j + 1); - - // For left cell boundaries - common_factor = (coord->J(i, j) + coord->J(i, j - 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); - - BoutReal flux_factor_lc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - BoutReal flux_factor_lm = - common_factor / (coord->dy(i, j - 1) * coord->J(i, j - 1)); - - BoutReal area_lc = common_factor * coord->dx(i, j) * coord->dz(i, j); -#endif - for (int k = 0; k < mesh->LocalNz; k++) { -#if BOUT_USE_METRIC_3D - // For right cell boundaries - BoutReal common_factor = - (coord->J(i, j, k) + coord->J(i, j + 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); - - BoutReal flux_factor_rc = - common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); - BoutReal flux_factor_rp = - common_factor / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); - - BoutReal area_rp = - common_factor * coord->dx(i, j + 1, k) * coord->dz(i, j + 1, k); - - // For left cell boundaries - common_factor = (coord->J(i, j, k) + coord->J(i, j - 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); - - BoutReal flux_factor_lc = - common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); - BoutReal flux_factor_lm = - common_factor / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); - - BoutReal area_lc = common_factor * coord->dx(i, j, k) * coord->dz(i, j, k); -#endif - - //////////////////////////////////////////// - // Reconstruct f at the cell faces - // This calculates s.R and s.L for the Right and Left - // face values on this cell - - // Reconstruct f at the cell faces - Stencil1D s; - s.c = f(i, j, k); - s.m = f(i, j - 1, k); - s.p = f(i, j + 1, k); - - cellboundary(s); // Calculate s.R and s.L - - //////////////////////////////////////////// - // Reconstruct v at the cell faces - Stencil1D sv; - sv.c = v(i, j, k); - sv.m = v(i, j - 1, k); - sv.p = v(i, j + 1, k); - - cellboundary(sv); // Calculate sv.R and sv.L - - //////////////////////////////////////////// - // Right boundary - - BoutReal flux; - - if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) { - // Last point in domain - - // Calculate velocity at right boundary (y+1/2) - BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j + 1, k)); - - BoutReal bndryval = 0.5 * (s.c + s.p); - if (fixflux) { - // Use mid-point to be consistent with boundary conditions - flux = bndryval * vpar; - } else { - // Add flux due to difference in boundary values - flux = s.R * vpar + wave_speed(i, j, k) * (s.R - bndryval); - } - - } else { - // Maximum wave speed in the two cells - BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k), - fabs(v(i, j, k)), fabs(v(i, j + 1, k))); - - flux = s.R * 0.5 * (sv.R + amax); - } - - result(i, j, k) += flux * flux_factor_rc; - result(i, j + 1, k) -= flux * flux_factor_rp; - - flow_ylow(i, j + 1, k) += flux * area_rp; - - //////////////////////////////////////////// - // Calculate at left boundary - - if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) { - // First point in domain - BoutReal bndryval = 0.5 * (s.c + s.m); - BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j - 1, k)); - if (fixflux) { - // Use mid-point to be consistent with boundary conditions - flux = bndryval * vpar; - } else { - // Add flux due to difference in boundary values - flux = s.L * vpar - wave_speed(i, j, k) * (s.L - bndryval); - } - } else { - - // Maximum wave speed in the two cells - BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k), - fabs(v(i, j, k)), fabs(v(i, j - 1, k))); - - flux = s.L * 0.5 * (sv.L - amax); - } - - result(i, j, k) -= flux * flux_factor_lc; - result(i, j - 1, k) += flux * flux_factor_lm; - - flow_ylow(i, j, k) += flux * area_lc; - } - } - } - if (are_unaligned) { - flow_ylow = fromFieldAligned(flow_ylow, "RGN_NOBNDRY"); - } - return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; -} - -/// Div ( a g Grad_perp(f) ) -- Perpendicular gradient-driven advection -/// -/// This version uses a slope limiter to calculate cell edge values of g in X, -/// the advects the upwind cell edge. -/// -/// 1st order upwinding is used in Y. -template -const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, - const Field3D& f) { - ASSERT2(a.getLocation() == f.getLocation()); - - Mesh* mesh = a.getMesh(); - - // Requires at least 2 communication guard cells in X, 1 in Y - ASSERT1(mesh->xstart >= 2); - ASSERT1(mesh->ystart >= 1); - - CellEdges cellboundary; - - Field3D result{zeroFrom(f)}; - - Coordinates* coord = f.getCoordinates(); - - // Flux in x - - for (int i = mesh->xstart - 1; i <= mesh->xend; i++) { - for (int j = mesh->ystart; j <= mesh->yend; j++) { - for (int k = 0; k < mesh->LocalNz; k++) { - // Calculate flux from i to i+1 - - const BoutReal gradient = f(i + 1, j, k) - f(i, j, k); - - // Mid-point average boundary value of 'a' - const BoutReal aedge = 0.5 * (a(i + 1, j, k) + a(i, j, k)); - BoutReal gedge; - if (((i == mesh->xstart - 1) and mesh->firstX()) - or ((i == mesh->xend) and mesh->lastX())) { - // Mid-point average boundary value of 'g' - gedge = 0.5 * (g(i + 1, j, k) + g(i, j, k)); - } else if (gradient > 0) { - // Flux is from (i+1) to (i) - // Reconstruct `g` at left of (i+1, j, k) - - Stencil1D sg; - sg.m = g(i, j, k); - sg.c = g(i + 1, j, k); - sg.p = g(i + 2, j, k); - cellboundary(sg); // Calculate sg.R and sg.L - - gedge = sg.L; - } else { - // Flux is from (i) to (i+1) - // Reconstruct `g` at right of (i, j, k) - - Stencil1D sg; - sg.m = g(i - 1, j, k); - sg.c = g(i, j, k); - sg.p = g(i + 1, j, k); - cellboundary(sg); // Calculate sg.R and sg.L - - gedge = sg.R; - } - - // Flux across cell edge - const BoutReal fout = gradient * aedge * gedge - * (coord->J(i, j) * coord->g11(i, j) - + coord->J(i + 1, j) * coord->g11(i + 1, j)) - / (coord->dx(i, j) + coord->dx(i + 1, j)); - - result(i, j, k) += fout / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= fout / (coord->dx(i + 1, j) * coord->J(i + 1, j)); - } - } - } - - // Y and Z fluxes require Y derivatives - - // Fields containing values along the magnetic field - Field3D fup(mesh), fdown(mesh); - Field3D aup(mesh), adown(mesh); - Field3D gup(mesh), gdown(mesh); - - // Values on this y slice (centre). - // This is needed because toFieldAligned may modify the field - Field3D ac = a; - Field3D gc = g; - Field3D fc = f; - - // Result of the Y and Z fluxes - Field3D yzresult(mesh); - yzresult.allocate(); - - if (f.hasParallelSlices() && a.hasParallelSlices() && g.hasParallelSlices()) { - // All inputs have yup and ydown - - fup = f.yup(); - fdown = f.ydown(); - - aup = a.yup(); - adown = a.ydown(); - - gup = g.yup(); - gdown = g.ydown(); - } else { - // At least one input doesn't have yup/ydown fields. - // Need to shift to/from field aligned coordinates - - fup = fdown = fc = toFieldAligned(f); - aup = adown = ac = toFieldAligned(a); - gup = gdown = gc = toFieldAligned(g); - yzresult.setDirectionY(YDirectionType::Aligned); - } - - // Y flux - - for (int i = mesh->xstart; i <= mesh->xend; i++) { - for (int j = mesh->ystart; j <= mesh->yend; j++) { - - BoutReal coef_u = - 0.5 - * (coord->g_23(i, j) / SQ(coord->J(i, j) * coord->Bxy(i, j)) - + coord->g_23(i, j + 1) / SQ(coord->J(i, j + 1) * coord->Bxy(i, j + 1))); - - BoutReal coef_d = - 0.5 - * (coord->g_23(i, j) / SQ(coord->J(i, j) * coord->Bxy(i, j)) - + coord->g_23(i, j - 1) / SQ(coord->J(i, j - 1) * coord->Bxy(i, j - 1))); - - for (int k = 0; k < mesh->LocalNz; k++) { - // Calculate flux between j and j+1 - int kp = (k + 1) % mesh->LocalNz; - int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; - - // Calculate Z derivative at y boundary - BoutReal dfdz = - 0.25 * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km)) - / coord->dz(i, j); - - // Y derivative - BoutReal dfdy = 2. * (fup(i, j + 1, k) - fc(i, j, k)) - / (coord->dy(i, j + 1) + coord->dy(i, j)); - - BoutReal aedge = 0.5 * (ac(i, j, k) + aup(i, j + 1, k)); - BoutReal gedge; - if ((j == mesh->yend) and mesh->lastY(i)) { - // Midpoint boundary value - gedge = 0.5 * (gc(i, j, k) + gup(i, j + 1, k)); - } else if (dfdy > 0) { - // Flux from (j+1) to (j) - gedge = gup(i, j + 1, k); - } else { - // Flux from (j) to (j+1) - gedge = gc(i, j, k); - } - - BoutReal fout = aedge * gedge * 0.5 - * (coord->J(i, j) * coord->g23(i, j) - + coord->J(i, j + 1) * coord->g23(i, j + 1)) - * (dfdz - coef_u * dfdy); - - yzresult(i, j, k) = fout / (coord->dy(i, j) * coord->J(i, j)); - - // Calculate flux between j and j-1 - dfdz = 0.25 - * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) - fdown(i, j - 1, km)) - / coord->dz(i, j); - - dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k)) - / (coord->dy(i, j) + coord->dy(i, j - 1)); - - aedge = 0.5 * (ac(i, j, k) + adown(i, j - 1, k)); - if ((j == mesh->ystart) and mesh->firstY(i)) { - gedge = 0.5 * (gc(i, j, k) + gdown(i, j - 1, k)); - } else if (dfdy > 0) { - gedge = gc(i, j, k); - } else { - gedge = gdown(i, j - 1, k); - } - - fout = aedge * gedge * 0.5 - * (coord->J(i, j) * coord->g23(i, j) - + coord->J(i, j - 1) * coord->g23(i, j - 1)) - * (dfdz - coef_d * dfdy); - - yzresult(i, j, k) -= fout / (coord->dy(i, j) * coord->J(i, j)); - } - } - } - - // Z flux - // Easier since all metrics constant in Z - - for (int i = mesh->xstart; i <= mesh->xend; i++) { - for (int j = mesh->ystart; j <= mesh->yend; j++) { - // Coefficient in front of df/dy term - BoutReal coef = coord->g_23(i, j) - / (coord->dy(i, j + 1) + 2. * coord->dy(i, j) + coord->dy(i, j - 1)) - / SQ(coord->J(i, j) * coord->Bxy(i, j)); - - for (int k = 0; k < mesh->LocalNz; k++) { - // Calculate flux between k and k+1 - int kp = (k + 1) % mesh->LocalNz; - - BoutReal gradient = - // df/dz - (fc(i, j, kp) - fc(i, j, k)) / coord->dz(i, j) - - // - g_yz * df/dy / SQ(J*B) - - coef - * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k) - - fdown(i, j - 1, kp)); - - BoutReal fout = gradient * 0.5 * (ac(i, j, kp) + ac(i, j, k)) - * ((gradient > 0) ? gc(i, j, kp) : gc(i, j, k)); - - yzresult(i, j, k) += fout / coord->dz(i, j); - yzresult(i, j, kp) -= fout / coord->dz(i, j); - } - } - } - // Check if we need to transform back - if (f.hasParallelSlices() && a.hasParallelSlices()) { - result += yzresult; - } else { - result += fromFieldAligned(yzresult); - } - - return result; -} - -} // namespace FV - #endif // DIV_OPS_H diff --git a/include/hermes_build_config.hxx.in b/include/hermes_build_config.hxx.in index 750100b00..e084258ef 100644 --- a/include/hermes_build_config.hxx.in +++ b/include/hermes_build_config.hxx.in @@ -10,12 +10,6 @@ namespace hermes { /// Slope limiter to use in advection operators using Limiter=FV::@HERMES_SLOPE_LIMITER@; const char* const limiter_typename = "@HERMES_SLOPE_LIMITER@"; - - // Method used in Div_par_K_Grad_par_mod - enum class ConductionMethod {Original, ProductJK, Harmonic}; - - constexpr ConductionMethod conduction_method = ConductionMethod::@HERMES_CONDUCTION_METHOD@; - const char* const conduction_typename = "@HERMES_CONDUCTION_METHOD@"; } #endif // HERMES_BUILD_CONFIG_HXX diff --git a/include/integrate.hxx b/include/integrate.hxx index ac18520c7..09a10114d 100644 --- a/include/integrate.hxx +++ b/include/integrate.hxx @@ -5,13 +5,13 @@ #include #include #include -#include +#include #include "../include/hermes_build_config.hxx" /// Get the first argument from a parameter pack template -auto firstArg(const Head &head, Tail... ) { +auto firstArg(const Head& head, Tail...) { return head; } @@ -57,7 +57,7 @@ BoutReal cellRight(BoutReal c, BoutReal m, BoutReal p) { /// /// Example /// Field3D Ne = ..., Te = ...; -/// +/// /// Field3D result = cellAverage( /// [](BoutReal Ne, BoutReal Te) {return Ne*Te;} // The function to evaluate /// Ne.getRegion("RGN_NOBNDRY") // The region to iterate over @@ -67,9 +67,9 @@ BoutReal cellRight(BoutReal c, BoutReal m, BoutReal p) { /// is the same as the input fields. /// template -auto cellAverage(Function func, const RegionType ®ion) { +auto cellAverage(Function func, const RegionType& region) { // Note: Capture by value or func and region go out of scope - return [=](const auto &... args) { + return [=](const auto&... args) { // Use the first argument to set the result mesh etc. Field3D result{emptyFrom(firstArg(args...))}; result.allocate(); @@ -84,10 +84,11 @@ auto cellAverage(Function func, const RegionType ®ion) { // Integrate in Y using Simpson's rule // Using limiter to calculate cell edge values - result[i] = - 4. / 6 * func((args[i])...) + - (Ji + J[ym]) / (12. * Ji) * func(cellLeft(args[i], args[ym], args[yp])...) + - (Ji + J[yp]) / (12. * Ji) * func(cellRight(args[i], args[ym], args[yp])...); + result[i] = 4. / 6 * func((args[i])...) + + (Ji + J[ym]) / (12. * Ji) + * func(cellLeft(args[i], args[ym], args[yp])...) + + (Ji + J[yp]) / (12. * Ji) + * func(cellRight(args[i], args[ym], args[yp])...); } return result; }; diff --git a/include/neutral_full_velocity.hxx b/include/neutral_full_velocity.hxx index 263dad807..ae3a03cb3 100644 --- a/include/neutral_full_velocity.hxx +++ b/include/neutral_full_velocity.hxx @@ -9,6 +9,8 @@ #include +#if not BOUT_USE_METRIC_3D + /// Neutral gas model, evolving three components of velocity as axisymmetric fields /// /// Evolves neutral density, pressure and velocity as Field2D quantities @@ -100,5 +102,11 @@ private: namespace { RegisterComponent registersolverneutralfullvelocity; } - +#else // not BOUT_USE_METRIC_3D +namespace { +RegisterUnavailableComponent + registersolverneutralfullvelocity("neutral_full_velocity", + "This not (yet) compatible with 3D metrics"); +} +#endif // not BOUT_USE_METRIC_3D #endif // NEUTRAL_FULL_VELOCITY_H diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index 34446056d..6b164bc0a 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -82,6 +82,8 @@ private: bool output_ddt; ///< Save time derivatives? bool diagnose; ///< Save additional diagnostics? + std::string conduction_method{"orginal"}; + // Flow diagnostics Field3D pf_adv_perp_xlow, pf_adv_perp_ylow, pf_adv_par_ylow; Field3D mf_adv_perp_xlow, mf_adv_perp_ylow, mf_adv_par_ylow; diff --git a/include/polarisation_drift.hxx b/include/polarisation_drift.hxx index 54e53c5a6..a302ab22a 100644 --- a/include/polarisation_drift.hxx +++ b/include/polarisation_drift.hxx @@ -2,6 +2,8 @@ #ifndef POLARISATION_DRIFT_H #define POLARISATION_DRIFT_H +#include + #include "component.hxx" class Laplacian; @@ -67,7 +69,7 @@ struct PolarisationDrift : public NamedComponent { private: std::unique_ptr phiSolver; // Laplacian solver in X-Z - Field2D Bsq; // Cached SQ(coord->Bxy) + Coordinates::FieldMetric Bsq; // Cached SQ(coord->Bxy) // Diagnostic outputs bool diagnose; ///< Save diagnostic outputs? diff --git a/include/relax_potential.hxx b/include/relax_potential.hxx index 5f0b75c68..ed7717497 100644 --- a/include/relax_potential.hxx +++ b/include/relax_potential.hxx @@ -2,7 +2,8 @@ #ifndef RELAX_POTENTIAL_H #define RELAX_POTENTIAL_H -#include +#include +#include #include "component.hxx" @@ -115,11 +116,11 @@ private: BoutReal phi_boundary_last_update; ///< Time when last updated bool phi_core_averagey; ///< Average phi core boundary in Y? - Field2D Bsq; ///< SQ(coord->Bxy) - Vector2D Curlb_B; ///< Curvature vector Curl(b/B) - BoutReal hyper_z; ///< Hyper-viscosity in Z - Field2D viscosity; ///< Perpendicular Kinematic viscosity - Field2D viscosity_par; ///< Parallel Kinematic viscosity + Coordinates::FieldMetric Bsq; ///< SQ(coord->Bxy) + VectorMetric Curlb_B; ///< Curvature vector Curl(b/B) + BoutReal hyper_z; ///< Hyper-viscosity in Z + Coordinates::FieldMetric viscosity; ///< Perpendicular Kinematic viscosity + Coordinates::FieldMetric viscosity_par; ///< Parallel Kinematic viscosity // Relax-potential related variables BoutReal lambda_1; ///< Relaxation parameters. NOTE: lambda_1 has dimensions! diff --git a/include/scale_timederivs.hxx b/include/scale_timederivs.hxx index 8eae9ae9c..5889d4ace 100644 --- a/include/scale_timederivs.hxx +++ b/include/scale_timederivs.hxx @@ -40,7 +40,7 @@ private: void transform_impl(GuardedOptions& state) override { auto* coord = bout::globals::mesh->getCoordinates(); - Field2D dl2 = coord->g_22 * SQ(coord->dy); + auto dl2 = coord->g_22 * SQ(coord->dy); // Scale by parallel heat conduction CFL timescale auto Te = get(state["species"]["e"]["temperature"]); diff --git a/include/vorticity.hxx b/include/vorticity.hxx index b1638ef09..d461bae58 100644 --- a/include/vorticity.hxx +++ b/include/vorticity.hxx @@ -4,12 +4,10 @@ #include #include -#include +#include #include "component.hxx" -#include - /// Evolve electron density in time /// struct Vorticity : public NamedComponent { @@ -48,7 +46,7 @@ struct Vorticity : public NamedComponent { /// If phi_boundary_relax is false, set the radial boundary to the sheath potential? /// - split_n0: bool, default false /// Split phi into n=0 and n!=0 components? - /// - viscosity: Field2D, default 0.0 + /// - viscosity: FieldMetric, default 0.0 /// Kinematic viscosity [m^2/s] /// - vort_dissipation: bool, default false /// Parallel dissipation of vorticity? @@ -129,12 +127,12 @@ private: bool split_n0; // Split phi into n=0 and n!=0 components std::unique_ptr laplacexy; // Laplacian solver in X-Y (n=0) - Field2D Bsq; // SQ(coord->Bxy) - Vector2D Curlb_B; // Curvature vector Curl(b/B) - BoutReal hyper_z; ///< Hyper-viscosity in Z - Field2D viscosity; ///< Kinematic viscosity - Field3D viscous_heating; ///< Heating due to kinematic viscosity - bool include_viscosity; ///< Is viscosity > 0? + Coordinates::FieldMetric Bsq; // SQ(coord->Bxy) + VectorMetric Curlb_B; // Curvature vector Curl(b/B) + BoutReal hyper_z; ///< Hyper-viscosity in Z + Coordinates::FieldMetric viscosity; ///< Kinematic viscosity + Field3D viscous_heating; ///< Heating due to kinematic viscosity + bool include_viscosity; ///< Is viscosity > 0? // Diagnostic outputs Field3D DivJdia, DivJcol; // Divergence of diamagnetic and collisional current diff --git a/src/braginskii_conduction.cxx b/src/braginskii_conduction.cxx index 45a2005d3..c29bd45b6 100644 --- a/src/braginskii_conduction.cxx +++ b/src/braginskii_conduction.cxx @@ -113,6 +113,9 @@ BraginskiiConduction::BraginskiiConduction(const std::string& name, Options& all } } substitutePermissions("sp", species); + + conduction_method = + alloptions["conduction_method"].withDefault(conduction_method); } void BraginskiiConduction::transform_impl(GuardedOptions& state) { @@ -292,7 +295,8 @@ void BraginskiiConduction::transform_impl(GuardedOptions& state) { // is calculated and removed separately set(species["kappa_par"], kappa_par); add(species["energy_source"], - Div_par_K_Grad_par_mod(kappa_par, T, flow_ylow_conduction, false)); + Div_par_K_Grad_par_mod(kappa_par, T, flow_ylow_conduction, false, + conduction_method)); add(species["energy_flow_ylow"], flow_ylow_conduction); if (state.isSection("fields") and state["fields"].isSet("Apar_flutter")) { diff --git a/src/braginskii_ion_viscosity.cxx b/src/braginskii_ion_viscosity.cxx index cb7a33748..b104a3184 100644 --- a/src/braginskii_ion_viscosity.cxx +++ b/src/braginskii_ion_viscosity.cxx @@ -21,6 +21,7 @@ #include #include #include +#include #include "../include/braginskii_ion_viscosity.hxx" #include "../include/component.hxx" @@ -95,7 +96,7 @@ BraginskiiIonViscosity::BraginskiiIonViscosity(const std::string& name, mesh->get(Curlb_B, "bxcv"); } catch (BoutException& e) { // May be 2D, reading as 3D - Vector2D curv2d; + VectorMetric curv2d; curv2d.covariant = false; mesh->get(curv2d, "bxcv"); Curlb_B = curv2d; @@ -145,9 +146,9 @@ void BraginskiiIonViscosity::transform_impl(GuardedOptions& state) { GuardedOptions allspecies = state["species"]; auto coord = mesh->getCoordinates(); - const Field2D Bxy = coord->Bxy; - const Field2D sqrtB = sqrt(Bxy); - const Field2D Grad_par_logB = Grad_par(log(Bxy)); + const auto Bxy = coord->Bxy; + const auto sqrtB = sqrt(Bxy); + const auto Grad_par_logB = Grad_par(log(Bxy)); // Loop through all species for (auto& kv : allspecies.getChildren()) { @@ -239,7 +240,7 @@ void BraginskiiIonViscosity::transform_impl(GuardedOptions& state) { // Parallel ion viscosity (4/3 * 0.96 coefficient) Field3D eta = 1.28 * P * tau; - const Field2D tau_av = DC(tau); + const auto tau_av = DC(tau); Field2D bounce_factor = 1.0; // if bounce_frequency = false, this factor does nothing to anything diff --git a/src/diamagnetic_drift.cxx b/src/diamagnetic_drift.cxx index b95fc3df9..041e4a6f1 100644 --- a/src/diamagnetic_drift.cxx +++ b/src/diamagnetic_drift.cxx @@ -52,10 +52,14 @@ DiamagneticDrift::DiamagneticDrift(std::string name, Options& alloptions, // Set drift to zero through sheath boundaries. // Flux through those cell faces should be set by sheath. for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { - Curlb_B.y(r.ind, mesh->ystart - 1) = -Curlb_B.y(r.ind, mesh->ystart); + for (int k = 0; k < Curlb_B.y.getNz(); ++k) { + Curlb_B.y(r.ind, mesh->ystart - 1, k) = -Curlb_B.y(r.ind, mesh->ystart, k); + } } for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { - Curlb_B.y(r.ind, mesh->yend + 1) = -Curlb_B.y(r.ind, mesh->yend); + for (int k = 0; k < Curlb_B.y.getNz(); ++k) { + Curlb_B.y(r.ind, mesh->yend + 1, k) = -Curlb_B.y(r.ind, mesh->yend, k); + } } // FIXME: density, pressure, and momentum will not be read even if diff --git a/src/div_ops.cxx b/src/div_ops.cxx index 56fed7315..2160e6606 100644 --- a/src/div_ops.cxx +++ b/src/div_ops.cxx @@ -30,7 +30,7 @@ #include #include #include -#include +#include #include #include #include @@ -43,7 +43,7 @@ using bout::globals::mesh; -const Field3D Div_par_diffusion_index(const Field3D& f, bool bndry_flux) { +Field3D Div_par_diffusion_index(const Field3D& f, bool bndry_flux) { Field3D result; result = 0.0; @@ -63,14 +63,15 @@ const Field3D Div_par_diffusion_index(const Field3D& f, bool bndry_flux) { continue; } } - BoutReal J = 0.5 * (coord->J(i, j) + coord->J(i, j + 1)); // Jacobian at boundary + BoutReal J = + 0.5 * (coord->J(i, j, k) + coord->J(i, j + 1, k)); // Jacobian at boundary BoutReal gradient = f(i, j + 1, k) - f(i, j, k); BoutReal flux = J * gradient; - result(i, j, k) += flux / coord->J(i, j); - result(i, j + 1, k) -= flux / coord->J(i, j + 1); + result(i, j, k) += flux / coord->J(i, j, k); + result(i, j + 1, k) -= flux / coord->J(i, j + 1, k); } } } @@ -114,8 +115,8 @@ void MC(Stencil1D& n) { * poloidal - If true, includes X-Y flows * positive - If true, limit advected quantity (n_in) to be positive */ -const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bndry_flux, - bool poloidal, bool positive) { +Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bndry_flux, + bool poloidal, bool positive) { Field3D result{0.0}; Coordinates* coord = mesh->getCoordinates(); @@ -158,13 +159,13 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd // 2) Calculate velocities on cell faces - BoutReal vU = coord->J(i, j) * (fmp - fpp) / coord->dx(i, j); // -J*df/dx - BoutReal vD = coord->J(i, j) * (fmm - fpm) / coord->dx(i, j); // -J*df/dx + BoutReal vU = coord->J(i, j, k) * (fmp - fpp) / coord->dx(i, j, k); // -J*df/dx + BoutReal vD = coord->J(i, j, k) * (fmm - fpm) / coord->dx(i, j, k); // -J*df/dx - BoutReal vR = 0.5 * (coord->J(i, j) + coord->J(i + 1, j)) * (fpp - fpm) - / coord->dz(i, j); // J*df/dz - BoutReal vL = 0.5 * (coord->J(i, j) + coord->J(i - 1, j)) * (fmp - fmm) - / coord->dz(i, j); // J*df/dz + BoutReal vR = 0.5 * (coord->J(i, j, k) + coord->J(i + 1, j, k)) * (fpp - fpm) + / coord->dz(i, j, k); // J*df/dz + BoutReal vL = 0.5 * (coord->J(i, j, k) + coord->J(i - 1, j, k)) * (fmp - fmm) + / coord->dz(i, j, k); // J*df/dz // 3) Calculate n on the cell faces. The sign of the // velocity determines which side is used. @@ -192,16 +193,18 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd // Flux in from boundary flux = vR * 0.5 * (n(i + 1, j, k) + n(i, j, k)); } - result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + result(i, j, k) += flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= + flux / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); } } else { // Not at a boundary if (vR > 0.0) { // Flux out into next cell BoutReal flux = vR * s.R; - result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + result(i, j, k) += flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= + flux / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); // if(i==mesh->xend) // output.write("Setting flux (%d,%d) : %e\n", @@ -225,16 +228,18 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd // Flux in from boundary flux = vL * 0.5 * (n(i - 1, j, k) + n(i, j, k)); } - result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); - result(i - 1, j, k) += flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); + result(i, j, k) -= flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i - 1, j, k) += + flux / (coord->dx(i - 1, j, k) * coord->J(i - 1, j, k)); } } else { // Not at a boundary if (vL < 0.0) { BoutReal flux = vL * s.L; - result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); - result(i - 1, j, k) += flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); + result(i, j, k) -= flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i - 1, j, k) += + flux / (coord->dx(i - 1, j, k) * coord->J(i - 1, j, k)); } } @@ -249,12 +254,12 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd MC(s); if (vU > 0.0) { - BoutReal flux = vU * s.R / (coord->J(i, j) * coord->dz(i, j)); + BoutReal flux = vU * s.R / (coord->J(i, j, k) * coord->dz(i, j, k)); result(i, j, k) += flux; result(i, j, kp) -= flux; } if (vD < 0.0) { - BoutReal flux = vD * s.L / (coord->J(i, j) * coord->dz(i, j)); + BoutReal flux = vD * s.L / (coord->J(i, j, k) * coord->dz(i, j, k)); result(i, j, k) -= flux; result(i, j, km) += flux; } @@ -307,13 +312,14 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd // Average dfdy to right X boundary BoutReal f_R = 0.5 - * ((coord->g11(i + 1, j) * coord->g23(i + 1, j) / SQ(coord->Bxy(i + 1, j))) + * ((coord->g11(i + 1, j, k) * coord->g23(i + 1, j, k) + / SQ(coord->Bxy(i + 1, j, k))) * dfdy(i + 1, j, k) - + (coord->g11(i, j) * coord->g23(i, j) / SQ(coord->Bxy(i, j))) + + (coord->g11(i, j, k) * coord->g23(i, j, k) / SQ(coord->Bxy(i, j, k))) * dfdy(i, j, k)); // Advection velocity across cell face - BoutReal Vx = 0.5 * (coord->J(i + 1, j) + coord->J(i, j)) * f_R; + BoutReal Vx = 0.5 * (coord->J(i + 1, j, k) + coord->J(i, j, k)) * f_R; // Fromm method BoutReal flux = Vx; @@ -334,8 +340,8 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd flux *= nval; } - result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + result(i, j, k) += flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= flux / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); } } } @@ -377,12 +383,13 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd // Average dfdx to upper Y boundary BoutReal f_U = 0.5 - * ((coord->g11(i, j + 1) * coord->g23(i, j + 1) / SQ(coord->Bxy(i, j + 1))) + * ((coord->g11(i, j + 1, k) * coord->g23(i, j + 1, k) + / SQ(coord->Bxy(i, j + 1, k))) * dfdx(i, j + 1, k) - + (coord->g11(i, j) * coord->g23(i, j) / SQ(coord->Bxy(i, j))) + + (coord->g11(i, j, k) * coord->g23(i, j, k) / SQ(coord->Bxy(i, j, k))) * dfdx(i, j, k)); - BoutReal Vy = -0.5 * (coord->J(i, j + 1) + coord->J(i, j)) * f_U; + BoutReal Vy = -0.5 * (coord->J(i, j + 1, k) + coord->J(i, j, k)) * f_U; if (mesh->firstY(i) && !mesh->periodicY(i) && (j == mesh->ystart - 1)) { // Lower y boundary. Allow flows out of the domain only @@ -413,8 +420,8 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd flux *= nval; } - yresult(i, j, k) += flux / (coord->dy(i, j) * coord->J(i, j)); - yresult(i, j + 1, k) -= flux / (coord->dy(i, j + 1) * coord->J(i, j + 1)); + yresult(i, j, k) += flux / (coord->dy(i, j, k) * coord->J(i, j, k)); + yresult(i, j + 1, k) -= flux / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); } } } @@ -424,7 +431,7 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D& n, const Field3D& f, bool bnd return result; } -const Field3D Div_Perp_Lap_FV_Index(const Field3D& as, const Field3D& fs) { +Field3D Div_Perp_Lap_FV_Index(const Field3D& as, const Field3D& fs) { Field3D result = 0.0; @@ -462,18 +469,18 @@ const Field3D Div_Perp_Lap_FV_Index(const Field3D& as, const Field3D& fs) { BoutReal gU = fs(i, j, kp) - fs(i, j, k); // Flow right - BoutReal flux = gR * 0.25 * (coord->J(i + 1, j) + coord->J(i, j)) - * (coord->dx(i + 1, j) + coord->dx(i, j)) + BoutReal flux = gR * 0.25 * (coord->J(i + 1, j, k) + coord->J(i, j, k)) + * (coord->dx(i + 1, j, k) + coord->dx(i, j, k)) * (as(i + 1, j, k) + as(i, j, k)); - result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); + result(i, j, k) += flux / (coord->dx(i, j, k) * coord->J(i, j, k)); // Flow left - flux = gL * 0.25 * (coord->J(i - 1, j) + coord->J(i, j)) - * (coord->dx(i - 1, j) + coord->dx(i, j)) + flux = gL * 0.25 * (coord->J(i - 1, j, k) + coord->J(i, j, k)) + * (coord->dx(i - 1, j, k) + coord->dx(i, j, k)) * (as(i - 1, j, k) + as(i, j, k)); - result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); + result(i, j, k) -= flux / (coord->dx(i, j, k) * coord->J(i, j, k)); // Flow up @@ -490,7 +497,7 @@ const Field3D Div_Perp_Lap_FV_Index(const Field3D& as, const Field3D& fs) { } /// Z diffusion in index space -const Field3D Div_Z_FV_Index(const Field3D& as, const Field3D& fs) { +Field3D Div_Z_FV_Index(const Field3D& as, const Field3D& fs) { Field3D result = 0.0; @@ -517,7 +524,7 @@ const Field3D Div_Z_FV_Index(const Field3D& as, const Field3D& fs) { } // *** USED *** -const Field3D D4DX4_FV_Index(const Field3D& f, bool bndry_flux) { +Field3D D4DX4_FV_Index(const Field3D& f, bool bndry_flux) { Field3D result = 0.0; Coordinates* coord = mesh->getCoordinates(); @@ -531,8 +538,8 @@ const Field3D D4DX4_FV_Index(const Field3D& f, bool bndry_flux) { BoutReal d3fdx3 = (f(i + 2, j, k) - 3. * f(i + 1, j, k) + 3. * f(i, j, k) - f(i - 1, j, k)); - BoutReal flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) - * (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; + BoutReal flux = 0.25 * (coord->dx(i, j, k) + coord->dx(i + 1, j, k)) + * (coord->J(i, j, k) + coord->J(i + 1, j, k)) * d3fdx3; if (mesh->lastX() && (i == mesh->xend)) { // Boundary @@ -547,8 +554,8 @@ const Field3D D4DX4_FV_Index(const Field3D& f, bool bndry_flux) { - (6. / 5) * f(i - 2, j, k) // f_2 ); - flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) - * (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; + flux = 0.25 * (coord->dx(i, j, k) + coord->dx(i + 1, j, k)) + * (coord->J(i, j, k) + coord->J(i + 1, j, k)) * d3fdx3; } else { // No fluxes through boundary @@ -556,8 +563,8 @@ const Field3D D4DX4_FV_Index(const Field3D& f, bool bndry_flux) { } } - result(i, j, k) += flux / (coord->J(i, j) * coord->dx(i, j)); - result(i + 1, j, k) -= flux / (coord->J(i + 1, j) * coord->dx(i + 1, j)); + result(i, j, k) += flux / (coord->J(i, j, k) * coord->dx(i, j, k)); + result(i + 1, j, k) -= flux / (coord->J(i + 1, j, k) * coord->dx(i + 1, j, k)); if (j == mesh->xstart) { // Left cell boundary, no flux through boundaries @@ -573,11 +580,12 @@ const Field3D D4DX4_FV_Index(const Field3D& f, bool bndry_flux) { + (6. / 5) * f(i + 2, j, k) // f_2 ); - flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) - * (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; + flux = 0.25 * (coord->dx(i, j, k) + coord->dx(i + 1, j, k)) + * (coord->J(i, j, k) + coord->J(i + 1, j, k)) * d3fdx3; - result(i, j, k) -= flux / (coord->J(i, j) * coord->dx(i, j)); - result(i - 1, j, k) += flux / (coord->J(i - 1, j) * coord->dx(i - 1, j)); + result(i, j, k) -= flux / (coord->J(i, j, k) * coord->dx(i, j, k)); + result(i - 1, j, k) += + flux / (coord->J(i - 1, j, k) * coord->dx(i - 1, j, k)); } } else { @@ -585,11 +593,12 @@ const Field3D D4DX4_FV_Index(const Field3D& f, bool bndry_flux) { d3fdx3 = (f(i + 1, j, k) - 3. * f(i, j, k) + 3. * f(i - 1, j, k) - f(i - 2, j, k)); - flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) - * (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; + flux = 0.25 * (coord->dx(i, j, k) + coord->dx(i + 1, j, k)) + * (coord->J(i, j, k) + coord->J(i + 1, j, k)) * d3fdx3; - result(i, j, k) -= flux / (coord->J(i, j) * coord->dx(i, j)); - result(i - 1, j, k) += flux / (coord->J(i - 1, j) * coord->dx(i - 1, j)); + result(i, j, k) -= flux / (coord->J(i, j, k) * coord->dx(i, j, k)); + result(i - 1, j, k) += + flux / (coord->J(i - 1, j, k) * coord->dx(i - 1, j, k)); } } } @@ -599,7 +608,7 @@ const Field3D D4DX4_FV_Index(const Field3D& f, bool bndry_flux) { return result; } -const Field3D D4DZ4_Index(const Field3D& f) { +Field3D D4DZ4_Index(const Field3D& f) { Field3D result; result.allocate(); BOUT_FOR(i, f.getRegion("RGN_NOBNDRY")) { @@ -615,7 +624,12 @@ const Field3D D4DZ4_Index(const Field3D& f) { * we would need the corner cell values to take Y derivatives along X edges * */ -const Field2D Laplace_FV(const Field2D& k, const Field2D& f) { +#if BOUT_USE_METRIC_3D +Field2D Laplace_FV([[maybe_unused]] const Field2D& k, [[maybe_unused]] const Field2D& f) { + throw BoutException("Not Implemented!"); +} +#else +Field2D Laplace_FV(const Field2D& k, const Field2D& f) { Field2D result; result.allocate(); @@ -663,9 +677,10 @@ const Field2D Laplace_FV(const Field2D& k, const Field2D& f) { } return result; } +#endif -const Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, - Field3D& flow_xlow, Field3D& flow_ylow) { +Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, Field3D& flow_xlow, + Field3D& flow_ylow) { ASSERT2(a.getLocation() == f.getLocation()); Mesh* mesh = a.getMesh(); @@ -707,7 +722,7 @@ const Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, result(i + 1, j, k) -= fout / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); // Flow will be positive in the positive coordinate direction - flow_xlow(i + 1, j, k) = -1.0 * fout * coord->dy(i, j) * coord->dz(i, j); + flow_xlow(i + 1, j, k) = -1.0 * fout * coord->dy(i, j, k) * coord->dz(i, j, k); } } } @@ -760,7 +775,7 @@ const Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, flow_ylow.setDirectionY(YDirectionType::Aligned); } - if (bout::build::use_metric_3d) { + if constexpr (bout::build::use_metric_3d) { // 3D Metric, need yup/ydown fields. // Requires previous communication of metrics // -- should insert communication here? @@ -849,7 +864,7 @@ const Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, yzresult(i, j, k) -= fout / (dyc(i, j, k) * Jc(i, j, k)); // Flow will be positive in the positive coordinate direction - flow_ylow(i, j, k) = -1.0 * fout * coord->dx(i, j) * coord->dz(i, j); + flow_ylow(i, j, k) = -1.0 * fout * coord->dx(i, j, k) * coord->dz(i, j, k); } } } @@ -895,7 +910,7 @@ const Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, // Div ( a Grad_perp(f) ) -- diffusion /// WARNING: Causes checkerboarding in neutral_mixed integrated test -const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { +Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { ASSERT2(a.getLocation() == f.getLocation()); Mesh* mesh = a.getMesh(); @@ -914,16 +929,16 @@ const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { for (int k = 0; k < mesh->LocalNz; k++) { // Calculate flux from i to i+1 - const BoutReal gradient = (coord->J(i, j) * coord->g11(i, j) - + coord->J(i + 1, j) * coord->g11(i + 1, j)) + const BoutReal gradient = (coord->J(i, j, k) * coord->g11(i, j, k) + + coord->J(i + 1, j, k) * coord->g11(i + 1, j, k)) * (f(i + 1, j, k) - f(i, j, k)) - / (coord->dx(i, j) + coord->dx(i + 1, j)); + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); // Use the upwind coefficient const BoutReal fout = gradient * ((gradient > 0) ? a(i + 1, j, k) : a(i, j, k)); - result(i, j, k) += fout / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= fout / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + result(i, j, k) += fout / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= fout / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); } } } @@ -961,55 +976,63 @@ const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { } // Y flux + const auto zstart = bout::build::use_metric_3d ? mesh->zstart : 0; + const auto zend = bout::build::use_metric_3d ? mesh->zend : 0; for (int i = mesh->xstart; i <= mesh->xend; i++) { for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = zstart; k <= zend; k++) { + BoutReal coef_u = + 0.5 + * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)) + + coord->g_23(i, j + 1, k) + / SQ(coord->J(i, j + 1, k) * coord->Bxy(i, j + 1, k))); - BoutReal coef_u = - 0.5 - * (coord->g_23(i, j) / SQ(coord->J(i, j) * coord->Bxy(i, j)) - + coord->g_23(i, j + 1) / SQ(coord->J(i, j + 1) * coord->Bxy(i, j + 1))); - - BoutReal coef_d = - 0.5 - * (coord->g_23(i, j) / SQ(coord->J(i, j) * coord->Bxy(i, j)) - + coord->g_23(i, j - 1) / SQ(coord->J(i, j - 1) * coord->Bxy(i, j - 1))); - - for (int k = 0; k < mesh->LocalNz; k++) { - // Calculate flux between j and j+1 - int kp = (k + 1) % mesh->LocalNz; - int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; - - // Calculate Z derivative at y boundary - BoutReal dfdz = - 0.25 * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km)) - / coord->dz(i, j); - - // Y derivative - BoutReal dfdy = 2. * (fup(i, j + 1, k) - fc(i, j, k)) - / (coord->dy(i, j + 1) + coord->dy(i, j)); - - BoutReal fout = 0.25 * (ac(i, j, k) + aup(i, j + 1, k)) - * (coord->J(i, j) * coord->g23(i, j) - + coord->J(i, j + 1) * coord->g23(i, j + 1)) - * (dfdz - coef_u * dfdy); - - yzresult(i, j, k) = fout / (coord->dy(i, j) * coord->J(i, j)); - - // Calculate flux between j and j-1 - dfdz = 0.25 - * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) - fdown(i, j - 1, km)) - / coord->dz(i, j); - - dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k)) - / (coord->dy(i, j) + coord->dy(i, j - 1)); - - fout = 0.25 * (ac(i, j, k) + adown(i, j - 1, k)) - * (coord->J(i, j) * coord->g23(i, j) - + coord->J(i, j - 1) * coord->g23(i, j - 1)) - * (dfdz - coef_d * dfdy); - - yzresult(i, j, k) -= fout / (coord->dy(i, j) * coord->J(i, j)); + BoutReal coef_d = + 0.5 + * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)) + + coord->g_23(i, j - 1, k) + / SQ(coord->J(i, j - 1, k) * coord->Bxy(i, j - 1, k))); + + const auto zstart = bout::build::use_metric_3d ? k : mesh->zstart; + const auto zend = bout::build::use_metric_3d ? k : mesh->zend; + for (int k = zstart; k <= zend; k++) { + // Calculate flux between j and j+1 + int kp = (k + 1) % mesh->LocalNz; + int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; + + // Calculate Z derivative at y boundary + BoutReal dfdz = + 0.25 * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km)) + / coord->dz(i, j, k); + + // Y derivative + BoutReal dfdy = 2. * (fup(i, j + 1, k) - fc(i, j, k)) + / (coord->dy(i, j + 1, k) + coord->dy(i, j, k)); + + BoutReal fout = 0.25 * (ac(i, j, k) + aup(i, j + 1, k)) + * (coord->J(i, j, k) * coord->g23(i, j, k) + + coord->J(i, j + 1, k) * coord->g23(i, j + 1, k)) + * (dfdz - coef_u * dfdy); + + yzresult(i, j, k) = fout / (coord->dy(i, j, k) * coord->J(i, j, k)); + + // Calculate flux between j and j-1 + dfdz = + 0.25 + * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) - fdown(i, j - 1, km)) + / coord->dz(i, j, k); + + dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k)) + / (coord->dy(i, j, k) + coord->dy(i, j - 1, k)); + + fout = 0.25 * (ac(i, j, k) + adown(i, j - 1, k)) + * (coord->J(i, j, k) * coord->g23(i, j, k) + + coord->J(i, j - 1, k) * coord->g23(i, j - 1, k)) + * (dfdz - coef_d * dfdy); + + yzresult(i, j, k) -= fout / (coord->dy(i, j, k) * coord->J(i, j, k)); + } } } } @@ -1019,28 +1042,33 @@ const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { for (int i = mesh->xstart; i <= mesh->xend; i++) { for (int j = mesh->ystart; j <= mesh->yend; j++) { - // Coefficient in front of df/dy term - BoutReal coef = coord->g_23(i, j) - / (coord->dy(i, j + 1) + 2. * coord->dy(i, j) + coord->dy(i, j - 1)) - / SQ(coord->J(i, j) * coord->Bxy(i, j)); + for (int k = zstart; k <= zend; k++) { + // Coefficient in front of df/dy term + BoutReal coef = + coord->g_23(i, j, k) + / (coord->dy(i, j + 1, k) + 2. * coord->dy(i, j, k) + coord->dy(i, j - 1, k)) + / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)); - for (int k = 0; k < mesh->LocalNz; k++) { - // Calculate flux between k and k+1 - int kp = (k + 1) % mesh->LocalNz; + const auto zstart = bout::build::use_metric_3d ? k : mesh->zstart; + const auto zend = bout::build::use_metric_3d ? k : mesh->zend; + for (int k = zstart; k <= zend; k++) { + // Calculate flux between k and k+1 + int kp = (k + 1) % mesh->LocalNz; - BoutReal gradient = - // df/dz - (fc(i, j, kp) - fc(i, j, k)) / coord->dz(i, j) + BoutReal gradient = + // df/dz + (fc(i, j, kp) - fc(i, j, k)) / coord->dz(i, j, k) - // - g_yz * df/dy / SQ(J*B) - - coef - * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k) - - fdown(i, j - 1, kp)); + // - g_yz * df/dy / SQ(J*B) + - coef + * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k) + - fdown(i, j - 1, kp)); - BoutReal fout = gradient * ((gradient > 0) ? ac(i, j, kp) : ac(i, j, k)); + BoutReal fout = gradient * ((gradient > 0) ? ac(i, j, kp) : ac(i, j, k)); - yzresult(i, j, k) += fout / coord->dz(i, j); - yzresult(i, j, kp) -= fout / coord->dz(i, j); + yzresult(i, j, k) += fout / coord->dz(i, j, k); + yzresult(i, j, kp) -= fout / coord->dz(i, j, k); + } } } } @@ -1104,9 +1132,10 @@ Field3D Div_a_Grad_perp_nonorthog(const Field3D& a, const Field3D& f, Field3D& f Field3D fddx_xhigh(mesh); fddx_xhigh.allocate(); for (int i = mesh->xstart - 1; i <= mesh->xend; i++) { - for (int j = mesh->ystart - 1; j <= mesh->yend + 1; - j++) { // Note: Including one guard cell - for (int k = 0; k < mesh->LocalNz; k++) { + // Note: Including one guard cell + for (int j = mesh->ystart - 1; j <= mesh->yend + 1; j++) { + for (int k = mesh->zstart; k <= mesh->zend; k++) { + fddx_xhigh(i, j, k) = 2. * (f(i + 1, j, k) - f(i, j, k)) / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); } @@ -1143,7 +1172,7 @@ Field3D Div_a_Grad_perp_nonorthog(const Field3D& a, const Field3D& f, Field3D& f flow_ylow.setDirectionY(YDirectionType::Aligned); } - if (bout::build::use_metric_3d) { + if constexpr (bout::build::use_metric_3d) { // 3D Metric, need yup/ydown fields. // Requires previous communication of metrics // -- should insert communication here? @@ -1274,7 +1303,7 @@ Field3D Div_a_Grad_perp_nonorthog(const Field3D& a, const Field3D& f, Field3D& f yzresult(i, j, k) -= fout / (dyc(i, j, k) * Jc(i, j, k)); // Flow will be positive in the positive coordinate direction - flow_ylow(i, j, k) = -1.0 * fout * coord->dx(i, j) * coord->dz(i, j); + flow_ylow(i, j, k) = -1.0 * fout * coord->dx(i, j, k) * coord->dz(i, j, k); } } } @@ -1403,7 +1432,7 @@ Field3D Div_a_Grad_perp_nonorthog(const Field3D& a, const Field3D& f, Field3D& f result(i + 1, j, k) -= fout / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); // Flow will be positive in the positive coordinate direction - flow_xlow(i + 1, j, k) = -1.0 * fout * coord->dy(i, j) * coord->dz(i, j); + flow_xlow(i + 1, j, k) = -1.0 * fout * coord->dy(i, j, k) * coord->dz(i, j, k); } } } @@ -1416,22 +1445,22 @@ Field3D Div_a_Grad_perp_nonorthog(const Field3D& a, const Field3D& f, Field3D& f /// Returns the flows in the final arguments /// /// Flows are always in the positive {x,y} direction -/// i.e xlow(i,j) is the flow into cell (i,j) from the left, -/// and the flow out of cell (i-1,j) to the right +/// i.e xlow(i,j, k) is the flow into cell (i,j, k) from the left, +/// and the flow out of cell (i-1,j, k) to the right /// /// ylow(i,j+1) /// ^ /// +---|---+ /// | | -/// xlow(i,j) -> (i,j) -> xlow(i+1,j) +/// xlow(i,j, k) -> (i,j, k) -> xlow(i+1,j, k) /// | ^ | /// +---|---+ -/// ylow(i,j) +/// ylow(i,j, k) /// /// /// WARNING: Causes checkerboarding in neutral_mixed integrated test -const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, - Field3D& flow_xlow, Field3D& flow_ylow) { +Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, + Field3D& flow_xlow, Field3D& flow_ylow) { ASSERT2(a.getLocation() == f.getLocation()); Mesh* mesh = a.getMesh(); @@ -1454,19 +1483,19 @@ const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, for (int k = 0; k < mesh->LocalNz; k++) { // Calculate flux from i to i+1 - const BoutReal gradient = (coord->J(i, j) * coord->g11(i, j) - + coord->J(i + 1, j) * coord->g11(i + 1, j)) + const BoutReal gradient = (coord->J(i, j, k) * coord->g11(i, j, k) + + coord->J(i + 1, j, k) * coord->g11(i + 1, j, k)) * (f(i + 1, j, k) - f(i, j, k)) - / (coord->dx(i, j) + coord->dx(i + 1, j)); + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); // Use the upwind coefficient const BoutReal fout = gradient * ((gradient > 0) ? a(i + 1, j, k) : a(i, j, k)); - result(i, j, k) += fout / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= fout / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + result(i, j, k) += fout / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= fout / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); // Flow will be positive in the positive coordinate direction - flow_xlow(i + 1, j, k) = -1.0 * fout * coord->dy(i, j) * coord->dz(i, j); + flow_xlow(i + 1, j, k) = -1.0 * fout * coord->dy(i, j, k) * coord->dz(i, j, k); } } } @@ -1505,88 +1534,101 @@ const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, } // Y flux + const auto zstart = bout::build::use_metric_3d ? mesh->zstart : 0; + const auto zend = bout::build::use_metric_3d ? mesh->zend : 0; for (int i = mesh->xstart; i <= mesh->xend; i++) { for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = zstart; k <= zend; k++) { + BoutReal coef_u = + 0.5 + * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)) + + coord->g_23(i, j + 1, k) + / SQ(coord->J(i, j + 1, k) * coord->Bxy(i, j + 1, k))); - BoutReal coef_u = - 0.5 - * (coord->g_23(i, j) / SQ(coord->J(i, j) * coord->Bxy(i, j)) - + coord->g_23(i, j + 1) / SQ(coord->J(i, j + 1) * coord->Bxy(i, j + 1))); - - BoutReal coef_d = - 0.5 - * (coord->g_23(i, j) / SQ(coord->J(i, j) * coord->Bxy(i, j)) - + coord->g_23(i, j - 1) / SQ(coord->J(i, j - 1) * coord->Bxy(i, j - 1))); - - for (int k = 0; k < mesh->LocalNz; k++) { - // Calculate flux between j and j+1 - int kp = (k + 1) % mesh->LocalNz; - int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; - - // Calculate Z derivative at y boundary - BoutReal dfdz = - 0.25 * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km)) - / coord->dz(i, j); - - // Y derivative - BoutReal dfdy = 2. * (fup(i, j + 1, k) - fc(i, j, k)) - / (coord->dy(i, j + 1) + coord->dy(i, j)); - - BoutReal fout = 0.25 * (ac(i, j, k) + aup(i, j + 1, k)) - * (coord->J(i, j) * coord->g23(i, j) - + coord->J(i, j + 1) * coord->g23(i, j + 1)) - * (dfdz - coef_u * dfdy); - - yzresult(i, j, k) = fout / (coord->dy(i, j) * coord->J(i, j)); - - // Calculate flux between j and j-1 - dfdz = 0.25 - * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) - fdown(i, j - 1, km)) - / coord->dz(i, j); - - dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k)) - / (coord->dy(i, j) + coord->dy(i, j - 1)); - - fout = 0.25 * (ac(i, j, k) + adown(i, j - 1, k)) - * (coord->J(i, j) * coord->g23(i, j) - + coord->J(i, j - 1) * coord->g23(i, j - 1)) - * (dfdz - coef_d * dfdy); - - yzresult(i, j, k) -= fout / (coord->dy(i, j) * coord->J(i, j)); - - // Flow will be positive in the positive coordinate direction - flow_ylow(i, j, k) = -1.0 * fout * coord->dx(i, j) * coord->dz(i, j); + BoutReal coef_d = + 0.5 + * (coord->g_23(i, j, k) / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)) + + coord->g_23(i, j - 1, k) + / SQ(coord->J(i, j - 1, k) * coord->Bxy(i, j - 1, k))); + + const auto zstart = bout::build::use_metric_3d ? k : mesh->zstart; + const auto zend = bout::build::use_metric_3d ? k : mesh->zend; + for (int k = zstart; k <= zend; k++) { + // Calculate flux between j and j+1 + int kp = (k + 1) % mesh->LocalNz; + int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; + + // Calculate Z derivative at y boundary + BoutReal dfdz = + 0.25 * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km)) + / coord->dz(i, j, k); + + // Y derivative + BoutReal dfdy = 2. * (fup(i, j + 1, k) - fc(i, j, k)) + / (coord->dy(i, j + 1, k) + coord->dy(i, j, k)); + + BoutReal fout = 0.25 * (ac(i, j, k) + aup(i, j + 1, k)) + * (coord->J(i, j, k) * coord->g23(i, j, k) + + coord->J(i, j + 1, k) * coord->g23(i, j + 1, k)) + * (dfdz - coef_u * dfdy); + + yzresult(i, j, k) = fout / (coord->dy(i, j, k) * coord->J(i, j, k)); + + // Calculate flux between j and j-1 + dfdz = + 0.25 + * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) - fdown(i, j - 1, km)) + / coord->dz(i, j, k); + + dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k)) + / (coord->dy(i, j, k) + coord->dy(i, j - 1, k)); + + fout = 0.25 * (ac(i, j, k) + adown(i, j - 1, k)) + * (coord->J(i, j, k) * coord->g23(i, j, k) + + coord->J(i, j - 1, k) * coord->g23(i, j - 1, k)) + * (dfdz - coef_d * dfdy); + + yzresult(i, j, k) -= fout / (coord->dy(i, j, k) * coord->J(i, j, k)); + + // Flow will be positive in the positive coordinate direction + flow_ylow(i, j, k) = -1.0 * fout * coord->dx(i, j, k) * coord->dz(i, j, k); + } } } } // Z flux // Easier since all metrics constant in Z - for (int i = mesh->xstart; i <= mesh->xend; i++) { for (int j = mesh->ystart; j <= mesh->yend; j++) { - // Coefficient in front of df/dy term - BoutReal coef = coord->g_23(i, j) - / (coord->dy(i, j + 1) + 2. * coord->dy(i, j) + coord->dy(i, j - 1)) - / SQ(coord->J(i, j) * coord->Bxy(i, j)); - for (int k = 0; k < mesh->LocalNz; k++) { - // Calculate flux between k and k+1 - int kp = (k + 1) % mesh->LocalNz; + for (int k = zstart; k <= zend; k++) { + // Coefficient in front of df/dy term + BoutReal coef = + coord->g_23(i, j, k) + / (coord->dy(i, j + 1, k) + 2. * coord->dy(i, j, k) + coord->dy(i, j - 1, k)) + / SQ(coord->J(i, j, k) * coord->Bxy(i, j, k)); - BoutReal gradient = - // df/dz - (fc(i, j, kp) - fc(i, j, k)) / coord->dz(i, j) + const auto zstart = bout::build::use_metric_3d ? k : mesh->zstart; + const auto zend = bout::build::use_metric_3d ? k : mesh->zend; + for (int k = zstart; k <= zend; k++) { + // Calculate flux between k and k+1 + int kp = (k + 1) % mesh->LocalNz; - // - g_yz * df/dy / SQ(J*B) - - coef - * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k) - - fdown(i, j - 1, kp)); + BoutReal gradient = + // df/dz + (fc(i, j, kp) - fc(i, j, k)) / coord->dz(i, j, k) - BoutReal fout = gradient * ((gradient > 0) ? ac(i, j, kp) : ac(i, j, k)); + // - g_yz * df/dy / SQ(J*B) + - coef + * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k) + - fdown(i, j - 1, kp)); - yzresult(i, j, k) += fout / coord->dz(i, j); - yzresult(i, j, kp) -= fout / coord->dz(i, j); + BoutReal fout = gradient * ((gradient > 0) ? ac(i, j, kp) : ac(i, j, k)); + + yzresult(i, j, k) += fout / coord->dz(i, j, k); + yzresult(i, j, kp) -= fout / coord->dz(i, j, k); + } } } } @@ -1601,8 +1643,8 @@ const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, return result; } -const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, const Field3D& f, - bool bndry_flux) { +Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, const Field3D& f, + bool bndry_flux) { Field3D result{0.0}; Coordinates* coord = mesh->getCoordinates(); @@ -1645,13 +1687,13 @@ const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, const Fi // 2) Calculate velocities on cell faces - BoutReal vU = coord->J(i, j) * (fmp - fpp) / coord->dx(i, j); // -J*df/dx - BoutReal vD = coord->J(i, j) * (fmm - fpm) / coord->dx(i, j); // -J*df/dx + BoutReal vU = coord->J(i, j, k) * (fmp - fpp) / coord->dx(i, j, k); // -J*df/dx + BoutReal vD = coord->J(i, j, k) * (fmm - fpm) / coord->dx(i, j, k); // -J*df/dx - BoutReal vR = 0.5 * (coord->J(i, j) + coord->J(i + 1, j)) * (fpp - fpm) - / coord->dz(i, j); // J*df/dz - BoutReal vL = 0.5 * (coord->J(i, j) + coord->J(i - 1, j)) * (fmp - fmm) - / coord->dz(i, j); // J*df/dz + BoutReal vR = 0.5 * (coord->J(i, j, k) + coord->J(i + 1, j, k)) * (fpp - fpm) + / coord->dz(i, j, k); // J*df/dz + BoutReal vL = 0.5 * (coord->J(i, j, k) + coord->J(i - 1, j, k)) * (fmp - fmm) + / coord->dz(i, j, k); // J*df/dz // 3) Calculate g on cell faces @@ -1687,16 +1729,18 @@ const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, const Fi flux = vR * 0.5 * (n(i + 1, j, k) + n(i, j, k)) * gR; } - result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + result(i, j, k) += flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= + flux / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); } } else { // Not at a boundary if (vR > 0.0) { // Flux out into next cell - const BoutReal flux = vR * s.R * gR; - result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); - result(i + 1, j, k) -= flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + BoutReal flux = vR * s.R * gR; + result(i, j, k) += flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= + flux / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); } } @@ -1716,16 +1760,18 @@ const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, const Fi // Flux in from boundary flux = vL * 0.5 * (n(i - 1, j, k) + n(i, j, k)) * gL; } - result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); - result(i - 1, j, k) += flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); + result(i, j, k) -= flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i - 1, j, k) += + flux / (coord->dx(i - 1, j, k) * coord->J(i - 1, j, k)); } } else { // Not at a boundary if (vL < 0.0) { - const BoutReal flux = vL * s.L * gL; - result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); - result(i - 1, j, k) += flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); + BoutReal flux = vL * s.L * gL; + result(i, j, k) -= flux / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i - 1, j, k) += + flux / (coord->dx(i - 1, j, k) * coord->J(i - 1, j, k)); } } @@ -1740,12 +1786,12 @@ const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, const Fi MC(s); if (vU > 0.0) { - const BoutReal flux = vU * s.R * gU / (coord->J(i, j) * coord->dz(i, j)); + BoutReal flux = vU * s.R * gU / (coord->J(i, j, k) * coord->dz(i, j, k)); result(i, j, k) += flux; result(i, j, kp) -= flux; } if (vD < 0.0) { - const BoutReal flux = vD * s.L * gD / (coord->J(i, j) * coord->dz(i, j)); + BoutReal flux = vD * s.L * gD / (coord->J(i, j, k) * coord->dz(i, j, k)); result(i, j, k) -= flux; result(i, j, km) += flux; } @@ -1755,129 +1801,3 @@ const Field3D Div_n_g_bxGrad_f_B_XZ(const Field3D& n, const Field3D& g, const Fi FV::communicateFluxes(result); return result; } - -Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, Field3D& flow_ylow, - bool bndry_flux) { - TRACE("FV::Div_par_K_Grad_par_mod"); - - ASSERT2(Kin.getLocation() == fin.getLocation()); - - Mesh* mesh = Kin.getMesh(); - - 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"); - - Field3D result{zeroFrom(f)}; - flow_ylow = zeroFrom(f); - - // K and f fields in yup and ydown directions - const auto& Kup = use_parallel_slices ? Kin.yup() : K; - const auto& Kdown = use_parallel_slices ? Kin.ydown() : K; - const auto& fup = use_parallel_slices ? fin.yup() : f; - const auto& fdown = use_parallel_slices ? fin.ydown() : f; - - Coordinates* coord = fin.getCoordinates(); - - BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { - // Calculate flux at upper surface - - const auto iyp = i.yp(); - const auto iym = i.ym(); - - if (bndry_flux || mesh->periodicY(i.x()) || !mesh->lastY(i.x()) - || (i.y() != mesh->yend)) { - BoutReal flux = 0.0; - - if constexpr (hermes::conduction_method == hermes::ConductionMethod::Original) { - 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]); - - const BoutReal gradient = - 2. * (fup[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); - - flux = c * J * gradient / g_22; - } else if constexpr (hermes::conduction_method - == hermes::ConductionMethod::ProductJK) { - // Intended to reduce sensitivity of result to K in small cells - const BoutReal cJ = - 0.5 * (K[i] * coord->J[i] + Kup[iyp] * coord->J[iyp]); // K * J at boundary - const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); - - const BoutReal gradient = - 2. * (fup[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); - - flux = cJ * gradient / g_22; - } else if constexpr (hermes::conduction_method - == hermes::ConductionMethod::Harmonic) { - // Harmonic average (serial resistance) - const BoutReal cond_i = K[i] * coord->J[i] / (coord->g_22[i] * coord->dy[i]); - const BoutReal cond_iyp = - Kup[iyp] * coord->J[iyp] / (coord->g_22[iyp] * coord->dy[iyp]); - const BoutReal denom = cond_i + cond_iyp; - - // Harmonic mean: series resistance of two half-cells - const BoutReal C_edge = - (std::abs(denom) > std::numeric_limits::epsilon()) - ? 2.0 * cond_i * cond_iyp / denom - : 0.0; - - flux = C_edge * (fup[iyp] - f[i]); - } - result[i] += flux / (coord->dy[i] * coord->J[i]); - } - - // Calculate flux at lower surface - if (bndry_flux || mesh->periodicY(i.x()) || !mesh->firstY(i.x()) - || (i.y() != mesh->ystart)) { - BoutReal flux = 0.0; - - if constexpr (hermes::conduction_method == hermes::ConductionMethod::Original) { - 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 - const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); - - const BoutReal gradient = - 2. * (f[i] - fdown[iym]) / (coord->dy[i] + coord->dy[iym]); - - flux = c * J * gradient / g_22; - } else if constexpr (hermes::conduction_method - == hermes::ConductionMethod::ProductJK) { - const BoutReal cJ = - 0.5 * (K[i] * coord->J[i] + Kdown[iym] * coord->J[iym]); // K * J at boundary - const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); - - const BoutReal gradient = - 2. * (f[i] - fdown[iym]) / (coord->dy[i] + coord->dy[iym]); - - flux = cJ * gradient / g_22; - } else if constexpr (hermes::conduction_method - == hermes::ConductionMethod::Harmonic) { - const BoutReal cond_i = K[i] * coord->J[i] / (coord->g_22[i] * coord->dy[i]); - const BoutReal cond_iym = - Kdown[iym] * coord->J[iym] / (coord->g_22[iym] * coord->dy[iym]); - const BoutReal denom = cond_i + cond_iym; - - const BoutReal C_edge = - (std::abs(denom) > std::numeric_limits::epsilon()) - ? 2.0 * cond_i * cond_iym / denom - : 0.0; - - flux = C_edge * (f[i] - fdown[iym]); - } - - result[i] -= flux / (coord->dy[i] * coord->J[i]); - flow_ylow[i] = -flux * coord->dx[i] * coord->dz[i]; - } - } - - if (!use_parallel_slices) { - // Shifted to field aligned coordinates, so need to shift back - result = fromFieldAligned(result, "RGN_NOBNDRY"); - flow_ylow = fromFieldAligned(flow_ylow); - } - - return result; -} diff --git a/src/electromagnetic.cxx b/src/electromagnetic.cxx index da8478559..6fce27005 100644 --- a/src/electromagnetic.cxx +++ b/src/electromagnetic.cxx @@ -1,8 +1,18 @@ #include "../include/electromagnetic.hxx" +#include "../include/component.hxx" +#include "../include/guarded_options.hxx" +#include "../include/permissions.hxx" + +#include #include +#include #include #include +#include +#include + +#include // Set the default acceptance tolerances for the Naulin solver. // These are used if the maximum iterations is reached. @@ -166,7 +176,7 @@ void Electromagnetic::transform_impl(GuardedOptions& state) { for (int z = mesh->zstart; z <= mesh->zend; z++) { rhs(x, y, z) = (weight * (Apar(x + 1, y, z) - Apar(x, y, z)) + (1 - weight) * (Apar(x + 2, y, z) - Apar(x + 1, y, z))) - / (sqrt(coords->g_11(x, y)) * coords->dx(x, y)); + / (sqrt(coords->g_11(x, y, z)) * coords->dx(x, y, z)); } } } @@ -176,7 +186,7 @@ void Electromagnetic::transform_impl(GuardedOptions& state) { for (int z = mesh->zstart; z <= mesh->zend; z++) { rhs(x, y, z) = (weight * (Apar(x, y, z) - Apar(x - 1, y, z)) + (1 - weight) * (Apar(x - 1, y, z) - Apar(x - 2, y, z))) - / sqrt(coords->g_11(x, y)) / coords->dx(x, y); + / sqrt(coords->g_11(x, y, z)) / coords->dx(x, y, z); } } } diff --git a/src/fieldline_geometry.cxx b/src/fieldline_geometry.cxx index bdef1cba8..399bce426 100644 --- a/src/fieldline_geometry.cxx +++ b/src/fieldline_geometry.cxx @@ -225,13 +225,14 @@ FieldlineGeometry::FieldlineGeometry(std::string name, Options& options, Solver* // inter-processor boundaries degrades the solution (and roughly doubles the // spread between different y-decompositions). effective_magnetic_field_strength // is defined over RGN_ALL, so the guard cells get a consistent value here. - for (int j = 0; j < mesh->LocalNy; ++j) { + BOUT_FOR(i, effective_magnetic_field_strength.getRegion("RGN_ALL")) { // N.b. // The Jacobian has units of [m / radian T], which is why we need an extra factor of Lnorm. - coord->J(0, j) = 1 / effective_magnetic_field_strength(0, j, 0) / Lnorm; + coord->J[i] = 1 / effective_magnetic_field_strength[i] / Lnorm; } // Fill the inter-processor guard cells with the neighbour's Jacobian so the // metric is continuous across rank boundaries in a parallel y-decomposition. + ASSERT2(not coord->J.isFci()); mesh->communicate(coord->J); // Parallel length of cell diff --git a/src/neutral_full_velocity.cxx b/src/neutral_full_velocity.cxx index 60118bb86..d5d285a74 100644 --- a/src/neutral_full_velocity.cxx +++ b/src/neutral_full_velocity.cxx @@ -10,6 +10,8 @@ #include +#if not BOUT_USE_METRIC_3D + using bout::globals::mesh; NeutralFullVelocity::NeutralFullVelocity(const std::string& name, Options& alloptions, @@ -857,3 +859,5 @@ void NeutralFullVelocity::outputVars(Options& state) { {"source", "neutral_full_velocity"}}); } } + +#endif diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index c12249ca7..61fd6a0ab 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -118,6 +118,9 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* .doc("Include neutral gas heat conduction?") .withDefault(true); + conduction_method = + options["conduction_method"].withDefault(conduction_method); + collisionality_override = options["collisionality_override"] .doc( @@ -541,7 +544,8 @@ void NeutralMixed::finally(const Options& state) { ddt(Pn) += (2. / 3) * Div_par_K_Grad_par_mod(kappa_n, Tn, // Parallel conduction ef_cond_par_ylow, - false); // No conduction through target boundary + false, // No conduction through target boundary + conduction_method); // Perpendicular conduction if (nonorthogonal_operators) { @@ -598,8 +602,8 @@ void NeutralMixed::finally(const Options& state) { Field3D viscosity_source = Div_par_K_Grad_par_mod( // Parallel viscosity eta_n, Vn, mf_visc_par_ylow, - false) // No viscosity through target boundary - ; + false, // No viscosity through target boundary + conduction_method); // Perpendicular viscosity if (nonorthogonal_operators) { diff --git a/src/recycling.cxx b/src/recycling.cxx index e7d737c16..81a2a053f 100644 --- a/src/recycling.cxx +++ b/src/recycling.cxx @@ -197,11 +197,11 @@ void Recycling::transform_impl(GuardedOptions& state) { // Get metric tensor components Coordinates* coord = mesh->getCoordinates(); - const Field2D& J = coord->J; - const Field2D& dy = coord->dy; - const Field2D& dx = coord->dx; - const Field2D& dz = coord->dz; - const Field2D& g_22 = coord->g_22; + const auto& J = coord->J; + const auto& dy = coord->dy; + const auto& dx = coord->dx; + const auto& dz = coord->dz; + const auto& g_22 = coord->g_22; for (auto& channel : channels) { const GuardedOptions species_from = state["species"][channel.from]; @@ -466,8 +466,8 @@ void Recycling::transform_impl(GuardedOptions& state) { for (int iz = 0; iz < mesh->LocalNz; iz++) { // Volume of cell adjacent to wall which will receive source - BoutReal volume = J(mesh->xend, iy) * dx(mesh->xend, iy) * dy(mesh->xend, iy) - * dz(mesh->xend, iy); + BoutReal volume = J(mesh->xend, iy, iz) * dx(mesh->xend, iy, iz) + * dy(mesh->xend, iy, iz) * dz(mesh->xend, iy, iz); // If cell is a pump, overwrite multiplier with pump multiplier BoutReal multiplier = channel.sol_multiplier; @@ -589,8 +589,8 @@ void Recycling::transform_impl(GuardedOptions& state) { for (int iz = 0; iz < mesh->LocalNz; iz++) { // Volume of cell adjacent to wall which will receive source - BoutReal volume = J(mesh->xstart, iy) * dx(mesh->xstart, iy) - * dy(mesh->xstart, iy) * dz(mesh->xstart, iy); + BoutReal volume = J(mesh->xstart, iy, iz) * dx(mesh->xstart, iy, iz) + * dy(mesh->xstart, iy, iz) * dz(mesh->xstart, iy, iz); // If cell is a pump, overwrite multiplier with pump multiplier BoutReal multiplier = channel.pfr_multiplier; diff --git a/src/relax_potential.cxx b/src/relax_potential.cxx index ace162277..0fc905d5d 100644 --- a/src/relax_potential.cxx +++ b/src/relax_potential.cxx @@ -198,7 +198,7 @@ RelaxPotential::RelaxPotential(std::string name, Options& alloptions, Solver* so // Read curvature vector try { // May be 2D, reading as 3D - Vector2D curv2d; + VectorMetric curv2d; curv2d.covariant = false; mesh->get(curv2d, "bxcv"); Curlb_B = curv2d; diff --git a/src/sheath_boundary.cxx b/src/sheath_boundary.cxx index 94fa0d558..fd79625fd 100644 --- a/src/sheath_boundary.cxx +++ b/src/sheath_boundary.cxx @@ -169,7 +169,8 @@ void SheathBoundary::transform_impl(GuardedOptions& state) { // Need electron properties // Not const because boundary conditions will be set - Field3D Ne = toFieldAligned(floor(GET_NOBOUNDARY(Field3D, electrons["density"]), 0.0)); + Field3D Ne = + toFieldAligned(Field3D(floor(GET_NOBOUNDARY(Field3D, electrons["density"]), 0.0))); Field3D Te = toFieldAligned(GET_NOBOUNDARY(Field3D, electrons["temperature"])); Field3D Pe = IS_SET_NOBOUNDARY(electrons["pressure"]) ? toFieldAligned(getNoBoundary(electrons["pressure"])) @@ -231,8 +232,8 @@ void SheathBoundary::transform_impl(GuardedOptions& state) { continue; // Skip electrons and non-charged ions } - const Field3D Ni = - toFieldAligned(floor(GET_NOBOUNDARY(Field3D, species["density"]), 0.0)); + const Field3D Ni = toFieldAligned( + Field3D(floor(GET_NOBOUNDARY(Field3D, species["density"]), 0.0))); const Field3D Ti = toFieldAligned(GET_NOBOUNDARY(Field3D, species["temperature"])); const BoutReal Mi = GET_NOBOUNDARY(BoutReal, species["AA"]); const BoutReal Zi = GET_NOBOUNDARY(BoutReal, species["charge"]); @@ -683,7 +684,8 @@ void SheathBoundary::transform_impl(GuardedOptions& state) { : 5. / 3; // Ratio of specific heats (ideal gas) // Density and temperature boundary conditions will be imposed (free) - Field3D Ni = toFieldAligned(floor(getNoBoundary(species["density"]), 0.0)); + Field3D Ni = + toFieldAligned(Field3D(floor(getNoBoundary(species["density"]), 0.0))); Field3D Ti = toFieldAligned(getNoBoundary(species["temperature"])); Field3D Pi = species.isSet("pressure") ? toFieldAligned(getNoBoundary(species["pressure"])) diff --git a/src/sheath_boundary_insulating.cxx b/src/sheath_boundary_insulating.cxx index 63d7d76ad..7d08d90e5 100644 --- a/src/sheath_boundary_insulating.cxx +++ b/src/sheath_boundary_insulating.cxx @@ -106,7 +106,8 @@ void SheathBoundaryInsulating::transform_impl(GuardedOptions& state) { // Need electron properties // Not const because boundary conditions will be set - Field3D Ne = toFieldAligned(floor(GET_NOBOUNDARY(Field3D, electrons["density"]), 0.0)); + Field3D Ne = + toFieldAligned(Field3D(floor(GET_NOBOUNDARY(Field3D, electrons["density"]), 0.0))); Field3D Te = toFieldAligned(GET_NOBOUNDARY(Field3D, electrons["temperature"])); Field3D Pe = IS_SET_NOBOUNDARY(electrons["pressure"]) ? toFieldAligned(getNoBoundary(electrons["pressure"])) @@ -240,7 +241,8 @@ void SheathBoundaryInsulating::transform_impl(GuardedOptions& state) { : 5. / 3; // Ratio of specific heats (ideal gas) // Density and temperature boundary conditions will be imposed (free) - Field3D Ni = toFieldAligned(floor(getNoBoundary(species["density"]), 0.0)); + Field3D Ni = + toFieldAligned(Field3D(floor(getNoBoundary(species["density"]), 0.0))); Field3D Ti = toFieldAligned(getNoBoundary(species["temperature"])); Field3D Pi = species.isSet("pressure") ? toFieldAligned(getNoBoundary(species["pressure"])) diff --git a/src/sheath_boundary_simple.cxx b/src/sheath_boundary_simple.cxx index f51d07436..fb0828a1b 100644 --- a/src/sheath_boundary_simple.cxx +++ b/src/sheath_boundary_simple.cxx @@ -161,7 +161,8 @@ void SheathBoundarySimple::transform_impl(GuardedOptions& state) { // Need electron properties // Not const because boundary conditions will be set - Field3D Ne = toFieldAligned(floor(GET_NOBOUNDARY(Field3D, electrons["density"]), 0.0)); + Field3D Ne = + toFieldAligned(Field3D(floor(GET_NOBOUNDARY(Field3D, electrons["density"]), 0.0))); Field3D Te = toFieldAligned(GET_NOBOUNDARY(Field3D, electrons["temperature"])); Field3D Pe = IS_SET_NOBOUNDARY(electrons["pressure"]) ? toFieldAligned(getNoBoundary(electrons["pressure"])) @@ -532,7 +533,8 @@ void SheathBoundarySimple::transform_impl(GuardedOptions& state) { const BoutReal Mi = get(species["AA"]); // Density and temperature boundary conditions will be imposed (free) - Field3D Ni = toFieldAligned(floor(getNoBoundary(species["density"]), 0.0)); + Field3D Ni = + toFieldAligned(Field3D(floor(getNoBoundary(species["density"]), 0.0))); Field3D Ti = toFieldAligned(getNoBoundary(species["temperature"])); Field3D Pi = species.isSet("pressure") ? toFieldAligned(getNoBoundary(species["pressure"])) diff --git a/src/snb_conduction.cxx b/src/snb_conduction.cxx index 562037caa..0bce5f876 100644 --- a/src/snb_conduction.cxx +++ b/src/snb_conduction.cxx @@ -1,5 +1,5 @@ -#include #include "../include/snb_conduction.hxx" +#include #include using bout::globals::mesh; @@ -13,7 +13,7 @@ void SNBConduction::transform_impl(GuardedOptions& state) { // SNB non-local heat flux. Also returns the Spitzer-Harm value for comparison // Note: Te in eV, Ne in Nnorm - Field2D dy_orig = mesh->getCoordinates()->dy; + auto dy_orig = mesh->getCoordinates()->dy; mesh->getCoordinates()->dy *= rho_s0; // Convert distances to m // Inputs in eV and m^-3 diff --git a/src/vorticity.cxx b/src/vorticity.cxx index 6fa351f0c..c723cc4af 100644 --- a/src/vorticity.cxx +++ b/src/vorticity.cxx @@ -190,7 +190,7 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) // Create an XY solver for n=0 component laplacexy = LaplaceXY::create(mesh); // Set coefficients for Boussinesq solve - laplacexy->setCoefs(average_atomic_mass / SQ(coord->Bxy), 0.0); + laplacexy->setCoefs(average_atomic_mass / SQ(DC(coord->Bxy)), 0.0); } phiSolver = Laplacian::create(&options["laplacian"]); // Set coefficients for Boussinesq solve diff --git a/tests/integrated/2D-production/runtest b/tests/integrated/2D-production/runtest index f0ffa1007..7b383704f 100755 --- a/tests/integrated/2D-production/runtest +++ b/tests/integrated/2D-production/runtest @@ -162,39 +162,39 @@ regions["sol_ring"] = ds.hermes.select_region( expected = {} expected["inner_lower_target"] = {} -expected["inner_lower_target"]["ddt(Pe)"] = [-8849462223.88488, -8831051488.549446, -8753984989.447958, -8691651051.890438, -8648458002.247812, -8626430667.731733, -8632920959.625614, -8680321583.759829, -8768088551.261625, -8910898239.90236, -9155532468.54299, -9577380530.673615, -10306375902.070868, -11556904030.756124, -13506109862.898067, -15915303326.324772] -expected["inner_lower_target"]["ddt(NVd+)"] = [-98778.13652086735, -98741.6946791423, -97872.80930592745, -97167.63081640242, -96679.2573085041, -96430.41335521628, -96504.07955420233, -97040.17129483374, -98032.75589236667, -99648.3845251484, -102417.9100220749, -107200.12878450433, -115484.96486891803, -129758.4815239607, -152179.87954933156, -180158.5386572483] -expected["inner_lower_target"]["ddt(Pd)"] = [9146362.327422312, 8953437.748600213, 8874293.65374496, 8810078.63725109, 8765619.543826174, 8742981.126763903, 8749713.236555804, 8798549.864154402, 8887870.312733471, 9029521.738910213, 9274148.969614554, 9701892.519143388, 10447952.419195276, 11736975.010078121, 13764025.118308637, 16316056.430827195] +expected["inner_lower_target"]["ddt(Pe)"] = [-8845832605.979357, -8827676518.631382, -8750535250.14947, -8688129612.055775, -8644878854.200537, -8622808808.260796, -8629285592.37284, -8676715550.043682, -8764551128.53605, -8907476976.73953, -9152310480.514866, -9574492059.159468, -10304063533.546902, -11555248968.147617, -13505015937.799826, -15914485257.738344] +expected["inner_lower_target"]["ddt(NVd+)"] = [-98705.23708360767, -98673.9007546511, -97803.51694181342, -97096.90190919275, -96607.37188077759, -96357.67143168558, -96431.06598997173, -96967.74415410355, -97961.7018230309, -99579.65574271265, -102353.17132880865, -107142.07033965031, -115438.45485979918, -129725.1504577338, -152157.8002683891, -180141.97641620558] +expected["inner_lower_target"]["ddt(Pd)"] = [9146362.356310116, 8953437.774993332, 8874293.680351352, 8810078.664210789, 8765619.57106727, 8742981.154214574, 8749713.264078587, 8798549.891571185, 8887870.339878503, 9029521.765566679, 9274148.99540396, 9701892.543372998, 10447952.44027508, 11736975.027389083, 13764025.132318461, 16316056.443996936] expected["inner_upper_target"] = {} -expected["inner_upper_target"]["ddt(Pe)"] = [-8823943202.471989, -8789962558.824818, -8701462180.468233, -8631501884.342407, -8584096490.494612, -8560632815.489453, -8568230476.558625, -8618987809.323042, -8711967060.982983, -8862219398.87059, -9117470243.350624, -9555162105.26555, -10310499744.83234, -11595894771.356598, -13589521540.56696, -16050099401.450678] -expected["inner_upper_target"]["ddt(NVd+)"] = [98491.93446997949, 98278.29513137124, 97280.76266255787, 96489.8861081698, 95954.22466234666, 95689.2920339197, 95775.43919427939, 96349.24642880238, 97400.37873759655, 99099.64444049557, 101988.63531253753, 106949.66414838922, 115533.3311221996, 130206.48717574387, 153147.28547105496, 181744.18707505654] -expected["inner_upper_target"]["ddt(Pd)"] = [9118155.256683575, 8911160.373716276, 8820314.935151355, 8748306.937132085, 8699548.710075507, 8675448.084100712, 8683314.865534816, 8735578.88173725, 8830237.124684855, 8979565.178662885, 9235121.321032396, 9679133.439105643, 10452293.114633715, 11777512.500410842, 13851591.682986833, 16459518.70842988] +expected["inner_upper_target"]["ddt(Pe)"] = [-8820439221.188353, -8786698782.75518, -8698119649.235191, -8628087028.10581, -8580623645.7536745, -8557115685.751221, -8564697641.114288, -8615482212.728687, -8708528794.11472, -8858896969.916811, -9114345190.07144, -9552366221.9814, -10308278335.600462, -11594322653.359823, -13588504141.79072, -16049346561.855524] +expected["inner_upper_target"]["ddt(NVd+)"] = [98421.55921469684, 98212.73630774661, 97213.62572461924, 96421.30029781346, 95884.47681712266, 95618.65613444758, 95704.48749041579, 96278.83901176142, 97331.31855863592, 99032.9029636494, 101925.84564307997, 106893.46729790306, 115488.65005047154, 130174.82486401015, 153126.74822205945, 181728.94183818184] +expected["inner_upper_target"]["ddt(Pd)"] = [9118155.284563677, 8911160.399155032, 8820314.96078519, 8748306.963089367, 8699548.736295765, 8675448.110532202, 8683314.89205486, 8735578.908176228, 8830237.15087853, 8979565.204394676, 9235121.34593747, 9679133.46250711, 10452293.134917185, 11777512.516958566, 13851591.696197322, 16459518.720820397] expected["outer_upper_target"] = {} -expected["outer_upper_target"]["ddt(Pe)"] = [-23404631798.069298, -23137193767.44033, -22923179006.227013, -22727204584.171326, -22546433346.453415, -22378877181.75347, -22223252346.68647, -22078008047.831448, -21932747513.197884, -21755225656.335003, -21503129658.7458, -21162306428.565918, -20783966457.00139, -20524064179.522243, -20647913060.552975, -21300873226.238827] -expected["outer_upper_target"]["ddt(NVd+)"] = [-270587.26283424057, -267654.541184525, -264982.8748315516, -262541.6545805871, -260293.2730227858, -258212.00558490504, -256281.63738355093, -254482.6700133599, -252686.72227618776, -250497.22032824415, -247397.45100650765, -243223.44476478483, -238615.5202110081, -235485.60084327328, -237063.90635116014, -245128.16869089202] -expected["outer_upper_target"]["ddt(Pd)"] = [25069506.380063523, 24303219.10077882, 24059044.550968695, 23835946.86831967, 23630485.45543071, 23440307.295676913, 23263928.078312278, 23099564.63759126, 22927387.415258802, 22697017.872459196, 22389646.36438664, 21996572.71869353, 21571316.86335293, 21283445.529979847, 21423114.968473796, 22175127.401863053] +expected["outer_upper_target"]["ddt(Pe)"] = [-23405141654.867226, -23137705857.450264, -22923691483.52866, -22727715961.565666, -22546942700.238037, -22379383762.18214, -22223755105.486065, -22078505082.629696, -21933235457.569817, -21755698424.56308, -21503574538.06294, -21162702232.697865, -20784286373.426586, -20524283070.39638, -20648026129.763287, -21300907861.139652] +expected["outer_upper_target"]["ddt(NVd+)"] = [-270597.71500810067, -267665.0342186318, -264993.3713478301, -262552.1245426391, -260303.69787063624, -258222.3702753069, -256291.9207626431, -254492.83344874342, -252696.6970482766, -250506.8816432458, -247406.53819903306, -243231.52465214147, -238622.04676400786, -235490.06468358944, -237066.21293258137, -245128.87611144106] +expected["outer_upper_target"]["ddt(Pd)"] = [25069506.369957812, 24303219.090774454, 24059044.541057177, 23835946.858519178, 23630485.445751753, 23440307.28612849, 23263928.068909798, 23099564.6283654, 22927387.406274293, 22697017.863850918, 22389646.356406078, 21996572.71173118, 21571316.857857622, 21283445.526306726, 21423114.966608323, 22175127.40130783] expected["outer_lower_target"] = {} -expected["outer_lower_target"]["ddt(Pe)"] = [-23023636859.765823, -22762414988.433784, -22552922729.9815, -22360407781.89368, -22182715047.34058, -22018546055.155735, -21865280266.065907, -21717135679.161137, -21568352970.922516, -21389971539.477703, -21141440412.13172, -20813185082.50465, -20463748473.246456, -20265071996.760475, -20492303272.269432, -21270149563.758247] -expected["outer_lower_target"]["ddt(NVd+)"] = [265848.7376228907, 262993.30363926815, 260388.13857312407, 257997.5625475745, 255794.2166940457, 253761.58908917048, 251867.11848241117, 250040.054005325, 248208.5951704684, 246017.33029848168, 242973.17640620537, 238967.84788367205, 234728.27653112903, 232354.6134366621, 235192.0088846579, 244782.71419614818] -expected["outer_lower_target"]["ddt(Pd)"] = [24623794.360936105, 23878010.145980082, 23639959.190310214, 23421532.86883631, 23220226.501549855, 23034527.89361795, 22861461.792449873, 22694571.06565057, 22519464.03739986, 22290014.549381863, 21988599.595274065, 21611326.215515602, 21219680.516759127, 21000344.294578698, 21253905.216931745, 22143468.831178702] +expected["outer_lower_target"]["ddt(Pe)"] = [-23024211769.782314, -22762988907.912186, -22553493728.780678, -22360973892.765034, -22183274498.013798, -22019097045.133156, -21865820839.130302, -21717663870.77216, -21568866580.914787, -21390465974.51861, -21141904016.868923, -20813596972.0476, -20464079753.56235, -20265295777.762238, -20492415301.17124, -21270181655.155045] +expected["outer_lower_target"]["ddt(NVd+)"] = [265860.5145628087, 263005.0549888104, 260399.8253579506, 258009.1449599475, 255805.6589073013, 253772.85468968746, 251878.16785758038, 250050.84729107004, 248219.0875902036, 246027.42770163465, 242982.6399843201, 238976.2510132908, 234735.0312343737, 232359.17504283693, 235194.29365029107, 244783.3695939909] +expected["outer_lower_target"]["ddt(Pd)"] = [24623794.349818617, 23878010.135046236, 23639959.17954063, 23421532.85825515, 23220226.49118077, 23034527.883486442, 22861461.78258554, 22694571.056085724, 22519464.028175615, 22290014.540602535, 21988599.587166544, 21611326.208455052, 21219680.511215176, 21000344.29092123, 21253905.215135977, 22143468.830689732] expected["outer_midplane_a"] = {} -expected["outer_midplane_a"]["ddt(Pe)"] = [2597007.0953292064, -116.85701867093881, -116.7954446943411, -116.76650568244983, -116.74901391282451, -116.73961111640538, -116.73516959560204, -116.73281146047266, -116.73123583333035, -116.72938895439181, -116.72644474856148, -116.72174378294572, -116.71324882279168, -116.69715025275985, -116.73011219751741, -3921049.789799078] -expected["outer_midplane_a"]["ddt(NVd+)"] = [2.6992534364870027e-18, -3.610414591271963e-12, -3.1075429659090278e-12, -3.795448520946963e-12, -3.18942306493015e-12, -3.1988862106514207e-12, -3.2015546508800496e-12, -3.2020293072139276e-12, -3.841930295861274e-12, -3.840529329431264e-12, -3.837461369408508e-12, -3.83039633854506e-12, -6.357833046136906e-12, -6.303257607375951e-12, -1.057203645504073e-11, 1.386011930097718e-06] -expected["outer_midplane_a"]["ddt(Pd)"] = [453.46419892649334, 453.265467039755, 453.0143813193228, 452.90058491335833, 452.83363677954753, 452.7971853131224, 452.7799642835784, 452.77081962798945, 452.7642718360932, 452.75667119841034, 452.7452496196667, 452.72701094106804, 452.6957748604165, 452.63327970182706, 452.6243535084463, 59504.39222505075] +expected["outer_midplane_a"]["ddt(Pe)"] = [2597007.0953292064, -116.85702230583456, -116.79544832923685, -116.76650568244983, -116.74901391282451, -116.73961111640538, -116.73516959560204, -116.73281146047266, -116.73123583333035, -116.72938895439181, -116.72644474856148, -116.72174378294572, -116.71324882279168, -116.69715025275985, -116.73011222437968, -3921049.789282036] +expected["outer_midplane_a"]["ddt(NVd+)"] = [5.704099965413376e-13, -3.6104145908564327e-12, -3.107542965908425e-12, -3.795448521186892e-12, -3.1894230652114743e-12, -3.198886210916926e-12, -3.201554650935935e-12, -3.2020293070870943e-12, -3.841930295734914e-12, -3.84052932943429e-12, -3.837461369497229e-12, -3.830396338769282e-12, -6.357833046350431e-12, -6.3032576075306815e-12, -1.0572036455487379e-11, 1.386012543092309e-06] +expected["outer_midplane_a"]["ddt(Pd)"] = [453.46419892649334, 453.2654670397951, 453.0143813193492, 452.90058491335833, 452.83363677954753, 452.7971853131224, 452.7799642835784, 452.77081962798945, 452.7642718360932, 452.75667119841034, 452.7452496196667, 452.72701094106804, 452.6957748604165, 452.6332797018316, 452.6243535084692, 59504.39222506831] expected["sol_boundary"] = {} -expected["sol_boundary"]["ddt(Pe)"] = [-15915303326.324772, -235004386.52630776, -2573359.417000814, -1080769.5432358312, -1160274.5074505645, -1169357.893068723, -1605960.9965739148, -2418089.977143572, -2828866.790314278, -2917258.6331448057, -2917258.633134451, -2828866.7903316407, -2418089.9771426367, -1605961.000760319, -1169357.870510763, -1159503.5861332845, -1079816.5590735574, -2652835.4499015827, -242161513.28153047, -16050099401.450678, -21300873226.238827, -727748830.3001844, -19488176.829850353, -1686865.3067635575, -1355938.0073427209, -1293342.7138832766, -1195775.7262589014, -1044754.9939097713, -984816.8616117283, -1051878.9814284493, -1047056.7567408985, -1436305.9698601882, -2592335.3721857346, -3526547.029351511, -3921049.789799078, -3921049.7897891193, -3526547.0293668155, -2592335.3721706825, -1436305.9697091824, -1047056.6968764504, -1051033.1578108056, -982797.697528408, -1040701.940399521, -1187733.2773297138, -1280968.5147299608, -1339991.8048691915, -1686713.8485614809, -20070876.73981682, -733201604.2318095, -21270149563.758247] -expected["sol_boundary"]["ddt(NVd+)"] = [-180158.5386572483, -0.25678617942429444, -0.00099210409754064, -4.730562770331288e-06, 8.201740934830418e-07, 9.071343171134917e-07, 1.5673919541158623e-06, 1.8075385890103469e-06, 8.912983712899969e-07, 1.6781590250093625e-07, -1.6781590256127795e-07, -8.912983716074666e-07, -1.807538589541688e-06, -1.5673917349744658e-06, -9.087030581108094e-07, -8.409599630718918e-07, 5.0119526728306216e-06, 0.0010377121754149231, 0.26274704018145206, 181744.18707505654, -245128.16869089202, -0.5579039030205757, -0.006801708347879607, -0.00011294478540997627, -4.426478006706961e-06, -3.51149723291996e-06, -4.1834992420089365e-06, -2.789998289367938e-06, 9.286804680494488e-08, 5.243480446088509e-07, 1.293417682170279e-06, 3.089026944996127e-06, 5.325450464967266e-06, 4.310220087640809e-06, 1.386011930097718e-06, -1.3860119307688146e-06, -4.310220089556725e-06, -5.325450466626116e-06, -3.0890272941982084e-06, -1.2962608521074998e-06, -5.482091823302218e-07, -1.37832989304325e-07, 2.780265264444143e-06, 4.141854889773126e-06, 3.404530834308758e-06, 4.382294851948851e-06, 0.00011876924073506516, 0.006958348608064335, 0.5609368596638125, 244782.71419614818] -expected["sol_boundary"]["ddt(Pd)"] = [16316056.430827195, 21167.94216123154, 18681.83487842866, 16683.89164464725, 18039.02195128105, 18053.13776185012, 24611.808523026742, 36869.25092966416, 43061.35535722714, 44388.440296405955, 44388.44029624992, 43061.35535746125, 36869.25092967708, 24611.808523346168, 18053.051358399396, 18029.121045960066, 16663.34235631647, 18637.633071300974, 21093.002229460708, 16459518.70842988, 22175127.401863053, 23919.878677851637, 22832.292622289282, 21795.015479328762, 20840.938423331074, 19931.234684259114, 18458.65729775797, 16175.701158480531, 15299.952784999443, 16333.711180658913, 16222.741743923862, 29991.14782618698, 47424.91238805747, 61500.55355367606, 59504.39222505075, 59504.39222495491, 61500.553553907375, 47424.91238794103, 29991.14782577931, 16222.554830613079, 16321.883642221796, 15270.150245949651, 16114.756381368254, 18337.532127800674, 19744.767699159085, 20596.91620828044, 21495.159752840038, 22479.439154042615, 23523.0810865654, 22143468.831178702] +expected["sol_boundary"]["ddt(Pe)"] = [-15914485257.738344, -234996802.73753864, -2573315.227434359, -1080769.396462712, -1160274.5162140653, -1169357.8931050405, -1605960.991131004, -2418089.9828573, -2828866.79127515, -2917258.6329875784, -2917258.6329772244, -2828866.7912925114, -2418089.9828563654, -1605960.9953075105, -1169357.8704219938, -1159503.6022381529, -1079816.3841319336, -2652792.608114189, -242154363.29138586, -16049346561.855524, -21300907861.139652, -727749518.8840947, -19488192.881806802, -1686865.7389404273, -1355937.9510951568, -1293342.7653578392, -1195775.718318139, -1044754.9799597275, -984816.8790375756, -1051878.9748483081, -1047056.7570313151, -1436305.9760177017, -2592335.366126241, -3526547.0281921206, -3921049.789282036, -3921049.7892757105, -3526547.028211061, -2592335.366114824, -1436305.9759502434, -1047056.6961715296, -1051033.1577193856, -982797.7071436741, -1040701.9326985655, -1187733.279879586, -1280968.4707431267, -1339991.8925005512, -1686714.1835470526, -20070892.057414804, -733202243.8484284, -21270181655.155045] +expected["sol_boundary"]["ddt(NVd+)"] = [-180141.97641620558, -0.2567663657539139, -0.0009920711309736195, -4.730446252244835e-06, 8.201802743722614e-07, 9.071343170783025e-07, 1.5673912971104926e-06, 1.8075380733161987e-06, 8.912989942713912e-07, 1.6781689626505285e-07, -1.6781689632539508e-07, -8.912989945888611e-07, -1.8075380738476062e-06, -1.5673912969662867e-06, -9.087034132484913e-07, -8.409599632286418e-07, 5.011837140237842e-06, 0.001037680632844035, 0.2627285876565175, 181728.94183818184, -245128.87611144106, -0.5579052555285287, -0.006801714785069861, -0.00011294493956370598, -4.426478006711512e-06, -3.5115011054837166e-06, -4.183490433228065e-06, -2.7900005968722983e-06, 9.286577921222597e-08, 5.243480446318717e-07, 1.2934170950549884e-06, 3.0890272939460754e-06, 5.325451353821953e-06, 4.310221785574232e-06, 1.386012543092309e-06, -1.3860125437634054e-06, -4.3102217874901485e-06, -5.325451355480662e-06, -3.089027294371886e-06, -1.2962608520583168e-06, -5.482091823544743e-07, -1.378306665525741e-07, 2.7802676314924144e-06, 4.14185789741942e-06, 3.404522918140009e-06, 4.382289906067185e-06, 0.00011876939800720628, 0.006958354664927835, 0.560938121559757, 244783.3695939909] +expected["sol_boundary"]["ddt(Pd)"] = [16316056.443996936, 21167.93953075991, 18681.834878808244, 16683.891644035477, 18039.021951616105, 18053.137761893468, 24611.80852280045, 36869.2509299062, 43061.35535725219, 44388.440296405955, 44388.44029624992, 43061.355357486296, 36869.25092991913, 24611.80852312003, 18053.05135843915, 18029.12104631104, 16663.342355703175, 18637.63307167044, 21092.999642608156, 16459518.720820397, 22175127.40130783, 23919.87865164244, 22832.292622285826, 21795.015479329148, 20840.938423336334, 19931.234684257288, 18458.657297737205, 16175.701158495007, 15299.952785020578, 16333.711180651033, 16222.741742983584, 29991.147826846307, 47424.91238816511, 61500.5535533456, 59504.39222506831, 59504.39222497247, 61500.55355357692, 47424.91238804867, 29991.147826438646, 16222.554829672934, 16321.88364221309, 15270.15024597168, 16114.75638138239, 18337.532127780454, 19744.767699157575, 20596.91620828544, 21495.159752840398, 22479.439154039752, 23523.081036799325, 22143468.830689732] expected["sol_ring"] = {} -expected["sol_ring"]["ddt(Pe)"] = [-8768088551.261625, -33197517.771432474, -34270.18480201952, -145.42679598603908, -112.78370438917457, -106.29256396898377, -135.79529603103316, -120.70842578547348, -117.5461850393553, -117.31954732388641, -117.3195473246545, -117.5461850395169, -120.70842578416652, -135.79531943747222, -106.26697363481617, -112.84132894906931, -145.3306157431527, -33966.474474843926, -32819414.138547227, -8711967060.982983, -21932747513.197884, -738543641.2932615, -16021286.994439086, -177929.79112491212, -965.1632234118106, -119.97909820884219, -119.91148533793915, -126.7233360862852, -131.600846346435, -111.23505036790579, -108.28315514484545, -130.41884347444665, -119.67717269647282, -117.20648507595878, -116.73123583333035, -116.7312358336583, -117.20648507481415, -119.67717269739042, -130.41884529165083, -108.23762204582864, -111.30025478061458, -131.60183969856658, -126.70037238478437, -119.8753386378501, -119.80520778772954, -886.4583145534157, -163512.84535809705, -14973432.53666446, -705914601.3896351, -21568352970.922516] -expected["sol_ring"]["ddt(NVd+)"] = [-98032.75589236667, -0.06937109113426526, -5.7621060006721315e-05, -1.5416103301579583e-08, -5.800811204853931e-11, 2.662459414558923e-11, 4.860318164906013e-11, -8.766558496349684e-11, -1.9647936340500678e-11, -1.2645971729581726e-12, 1.2645971729941652e-12, 1.96479363406296e-11, 8.766558496122037e-11, -4.8368383634201935e-11, -2.654190906792441e-11, 5.812606233673052e-11, 1.5326944688842124e-08, 5.7061367935128635e-05, 0.06916818230231037, 97400.37873759655, -252686.72227618776, -0.5638888269552669, -0.006351576349484462, -8.142885231450911e-05, -4.922229247682577e-07, -1.0550940068772505e-09, 1.848476187145342e-10, 1.4909329445488978e-10, -7.276699197769393e-11, -3.17216721134545e-11, 2.2562199253482067e-11, 4.70379605630196e-11, -8.95091107599747e-11, -2.3647983200904405e-11, -3.841930295861274e-12, 3.8419302958427794e-12, 2.364798319988606e-11, 8.950911076145584e-11, -4.703796055409992e-11, -2.2646386571154443e-11, 3.1895228416937104e-11, 7.34011163613448e-11, -1.4934926274033144e-10, -1.8453534248570585e-10, 9.45965969566738e-10, 4.493438197213855e-07, 7.546030398845405e-05, 0.006002655490103674, 0.5460008969160347, 248208.5951704684] -expected["sol_ring"]["ddt(Pd)"] = [8887870.312733471, 472.20514261453576, 507.04336100390327, 551.0798383315962, 437.45608173504183, 412.2775048346568, 526.6828154678757, 468.20426453097656, 455.9287460035725, 455.0441495392622, 455.0441495422412, 455.92874600419924, 468.20426452590783, 526.6828154883442, 412.17833580226494, 437.6794613945601, 550.7391898636769, 506.8512228811758, 472.086393313632, 8830237.124684855, 22927387.415258802, 468.0371805773448, 454.02111354714, 455.1378079742271, 457.02345381985754, 459.89133973573183, 465.2014846146613, 491.44871128065256, 510.4186700624353, 431.44671945599487, 419.99751237069904, 505.8302295134958, 464.2012908035447, 454.61399787695984, 452.7642718360932, 452.76427183736547, 454.61399787252054, 464.2012908071038, 505.8302295055703, 419.820913680721, 431.6996314410605, 510.4223924249515, 491.3595099916829, 465.05802958550476, 459.7804349367682, 456.9945327577362, 455.1349952772563, 454.0161354495289, 467.66861473231893, 22519464.03739986] +expected["sol_ring"]["ddt(Pe)"] = [-8764551128.53605, -33189045.712066155, -34261.4735287216, -145.42582352877366, -112.78370529789852, -106.29256124281197, -135.79530693572042, -120.70842578547348, -117.5461850393553, -117.31954732388641, -117.3195473246545, -117.5461850395169, -120.70842578416652, -135.79531655951223, -106.26697986457468, -112.84132258822429, -145.3297674846184, -33958.029479271056, -32811221.083995894, -8708528794.11472, -21933235457.569817, -738553316.4868612, -16021495.819874318, -177932.21578270002, -965.1762376045174, -119.97909820884199, -119.91149260773062, -126.72334337045871, -131.6008408940914, -111.23505173099169, -108.28315832537922, -130.41883620465515, -119.67717269627522, -117.20648507595878, -116.73123583333035, -116.7312358336583, -117.20648870970989, -119.67717269719277, -130.41883620441146, -108.23762477200044, -111.3002556893385, -131.60183606367085, -126.70037965457584, -119.8753386378501, -119.80520778772939, -886.4711695079761, -163515.20187899197, -14973641.535221634, -705924501.9418006, -21568866580.914787] +expected["sol_ring"]["ddt(NVd+)"] = [-97961.7018230309, -0.06932889049195282, -5.760596074158181e-05, -1.5411746630051635e-08, -5.800811209296945e-11, 2.6624594145130448e-11, 4.860318169889781e-11, -8.766558505430314e-11, -1.9647936332050324e-11, -1.2645971731316317e-12, 1.2645971731676244e-12, 1.9647936332179244e-11, 8.766558505202667e-11, -4.8368383684039606e-11, -2.6541909067464324e-11, 5.812606238102515e-11, 1.532301145909176e-08, 5.704674781596614e-05, 0.06912723932410765, 97331.31855863592, -252696.6970482766, -0.5639075183097414, -0.006351662039486548, -8.142996522013759e-05, -4.922302565585815e-07, -1.05509400680854e-09, 1.8484761871453394e-10, 1.490932942458208e-10, -7.276699187293421e-11, -3.172167211938507e-11, 2.2562199256073062e-11, 4.674579921196455e-11, -8.950911067460031e-11, -2.305678370471812e-11, -3.841930295734914e-12, 3.8419302957164195e-12, 2.3056783703725228e-11, 8.950911067608148e-11, -4.6745799203118146e-11, -2.26463865737254e-11, 3.189522842284073e-11, 7.340111625726334e-11, -1.4934926252416104e-10, -1.8453534248570554e-10, 9.459659694531288e-10, 4.493511083644922e-07, 7.546139396995623e-05, 0.006002742175163397, 0.5460202864500181, 248219.0875902036] +expected["sol_ring"]["ddt(Pd)"] = [8887870.339878503, 472.20244442311423, 507.0433609748652, 551.0798384810611, 437.456081703489, 412.27750480309055, 526.6828155876699, 468.2042645123386, 455.9287460036652, 455.0441495392622, 455.0441495422412, 455.92874600429195, 468.2042645072699, 526.6828156082327, 412.17833577028745, 437.67946136343176, 550.7391900104835, 506.8512228527787, 472.0837885828307, 8830237.15087853, 22927387.406274293, 468.0376577843274, 454.02111354746745, 455.137807974162, 457.02345381976073, 459.89133973502845, 465.2014846199154, 491.4487112836086, 510.41867001494194, 431.4467194700936, 419.9975124058364, 505.83022940655655, 464.2012908152233, 454.6139978768262, 452.7642718360932, 452.76427183736547, 454.61399787238696, 464.2012908187824, 505.8302293984267, 419.82091371593066, 431.6996314548898, 510.42239237834195, 491.3595099946218, 465.0580295906958, 459.78043493608965, 456.99453275763125, 455.1349952772932, 454.0161354498959, 467.66905229922514, 22519464.028175615] # fmt: on # Generate test data or run test diff --git a/tests/unit/fake_mesh.hxx b/tests/unit/fake_mesh.hxx index cd44d4785..1be1cdd43 100644 --- a/tests/unit/fake_mesh.hxx +++ b/tests/unit/fake_mesh.hxx @@ -175,8 +175,12 @@ public: RangeIterator iterateBndryUpperInnerY() const override { return RangeIterator(); } bool hasBndryLowerY() const override { return false; } bool hasBndryUpperY() const override { return false; } - void addBoundary(BoundaryRegionBase* region) override { boundaries.push_back(region); } - std::vector getBoundaries() const override { return boundaries; } + void addBoundary(std::shared_ptr region) override { + boundaries.push_back(region); + } + std::vector> getBoundaries() const override { + return boundaries; + } std::vector> getBoundariesPar(BoundaryParType UNUSED(type)) const override { return std::vector>(); @@ -268,7 +272,7 @@ public: using Mesh::msg_len; private: - std::vector boundaries; + std::vector> boundaries; }; /// FakeGridDataSource provides a non-null GridDataSource* source to use with FakeMesh, to diff --git a/tests/unit/test_anomalous_diffusion.cxx b/tests/unit/test_anomalous_diffusion.cxx index 8147cd1d8..561757953 100644 --- a/tests/unit/test_anomalous_diffusion.cxx +++ b/tests/unit/test_anomalous_diffusion.cxx @@ -74,7 +74,7 @@ TEST_F(AnomalousDiffusionTest, ParticleDiffusion) { "RGN_NOBNDRY")); // Expect the sum over all cells of density source to be zero - Field2D dV = coords->J * coords->dx * coords->dy * coords->dz; // Cell volume + auto dV = coords->J * coords->dx * coords->dy * coords->dz; // Cell volume Field3D source = get(state["species"]["h"]["density_source"]); BoutReal integral = 0.0; diff --git a/tests/unit/test_fieldline_geometry.cxx b/tests/unit/test_fieldline_geometry.cxx index 3f4c1dfbe..c7d5c42ba 100644 --- a/tests/unit/test_fieldline_geometry.cxx +++ b/tests/unit/test_fieldline_geometry.cxx @@ -140,14 +140,14 @@ TEST_F(FieldlineGeometryTest, SetsCoordinatesJacobianAndBxy) { // and Lnorm = 1, so J = 1/Btotal. const BoutReal Btotal = sqrt(0.5 * 0.5 + 1.0 * 1.0); for (int j = mesh->ystart; j <= mesh->yend; ++j) { - ASSERT_NEAR(coord->J(0, j), 1.0 / Btotal, 1e-12); + ASSERT_NEAR(coord->J(0, j, 0), 1.0 / Btotal, 1e-12); } // J must also be set in the guard cells (not left at the default mesh value), // otherwise the metric is discontinuous at domain / inter-processor // boundaries. Check the full local range including guard cells. for (int j = 0; j < mesh->LocalNy; ++j) { - ASSERT_NEAR(coord->J(0, j), 1.0 / Btotal, 1e-12); + ASSERT_NEAR(coord->J(0, j, 0), 1.0 / Btotal, 1e-12); } } @@ -205,7 +205,7 @@ TEST_F(FieldlineGeometryTest, JacobianNotOneWhenBnormDoesNotMatchUpstreamField) Coordinates* coord = mesh->getCoordinates(); const BoutReal expected_J = 1.0 / sqrt(10.0); for (int j = mesh->ystart; j <= mesh->yend; ++j) { - ASSERT_NEAR(coord->J(0, j), expected_J, 1e-12); + ASSERT_NEAR(coord->J(0, j, 0), expected_J, 1e-12); } } @@ -219,7 +219,7 @@ TEST_F(FieldlineGeometryTest, JacobianIsOneWhenBnormMatchesUpstreamField) { Coordinates* coord = mesh->getCoordinates(); for (int j = mesh->ystart; j <= mesh->yend; ++j) { - ASSERT_NEAR(coord->J(0, j), 1.0, 1e-12); + ASSERT_NEAR(coord->J(0, j, 0), 1.0, 1e-12); } } diff --git a/tests/unit/test_neutral_full_velocity.cxx b/tests/unit/test_neutral_full_velocity.cxx index 1c10e8211..043d1273d 100644 --- a/tests/unit/test_neutral_full_velocity.cxx +++ b/tests/unit/test_neutral_full_velocity.cxx @@ -8,6 +8,7 @@ #include // For generating functions +#if not BOUT_USE_METRIC_3D /// Global mesh namespace bout { namespace globals { @@ -89,7 +90,9 @@ TEST_F(NeutralFullVelocityTest, CreateComponentRequiresBpxy) { static_cast(bout::globals::mesh) ->setGridDataSource(new FakeGridDataSource{{ - {"Rxy", 0.0}, {"Zxy", 1.0}, {"hthe", 1.0}, + {"Rxy", 0.0}, + {"Zxy", 1.0}, + {"hthe", 1.0}, // Missing Bpxy }}); @@ -316,3 +319,5 @@ TEST_F(NeutralFullVelocityTest, OutputVars) { ASSERT_TRUE(outputs.isSet("Urx")); ASSERT_TRUE(outputs.isSet("Tyr")); } + +#endif