diff --git a/.ci_fedora.sh b/.ci_fedora.sh index e1dc563c6c..8214a40ccc 100755 --- a/.ci_fedora.sh +++ b/.ci_fedora.sh @@ -47,6 +47,8 @@ test . != ".$1" && mpi="$1" || mpi=openmpi export OMP_NUM_THREADS=1 cd cd BOUT-dev + python3 -m ensurepip + python3 -m pip install git+https://github.com/boutproject/zoidberg@better-metric echo "starting configure" time cmake -S . -B build -DBOUT_USE_PETSC=ON \ -DBOUT_UPDATE_GIT_SUBMODULE=OFF \ diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 02f018b873..8919f82df1 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -166,6 +166,7 @@ jobs: run: | python -m pip install --upgrade pip setuptools python -m pip install -r requirements.txt + python -m pip install git+https://github.com/boutproject/zoidberg@better-metric - name: Cache Zenodo test data uses: actions/cache@v5 @@ -187,7 +188,7 @@ jobs: run: BUILD_PETSC=${{ matrix.config.build_petsc }} ./.build_petsc_for_ci.sh ${{ matrix.config.build_petsc_branch }} - name: Build ADIOS2 - run: BUILD_ADIOS2=${{ matrix.config.build_adios2 }} ./.build_adios2_for_ci.sh + run: BUILD_ADIOS2=${{ matrix.config.build_adios2 }} ./.build_adios2_for_ci.sh - name: Build BOUT++ run: UNIT_ONLY=${{ matrix.config.unit_only }} ./.ci_with_cmake.sh ${{ matrix.config.cmake_options }} diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index e7ead42ee5..6fc97a03dd 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -3,16 +3,16 @@ * * ChangeLog * ========= - * + * * 2014-11-10 Ben Dudson * * Created by separating metric from Mesh * - * + * ************************************************************************** * Copyright 2014-2025 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -160,7 +160,7 @@ public: const std::string& method = "DEFAULT", const std::string& region = "RGN_NOBNDRY"); - Field3D DDY(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, + Field3D DDY(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", const std::string& region = "RGN_NOBNDRY") const; @@ -172,7 +172,7 @@ public: FieldMetric Grad_par(const Field2D& var, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT"); - Field3D Grad_par(const Field3D& var, CELL_LOC outloc = CELL_DEFAULT, + Field3D Grad_par(const Field3DParallel& var, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT"); /// Advection along magnetic field V*b.Grad(f) @@ -180,7 +180,7 @@ public: CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT"); - Field3D Vpar_Grad_par(const Field3D& v, const Field3D& f, + Field3D Vpar_Grad_par(const Field3D& v, const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT"); @@ -188,14 +188,14 @@ public: FieldMetric Div_par(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT"); - Field3D Div_par(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, + Field3D Div_par(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT"); // Second derivative along magnetic field FieldMetric Grad2_par2(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT"); - Field3D Grad2_par2(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, + Field3D Grad2_par2(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT"); // Perpendicular Laplacian operator, using only X-Z derivatives // NOTE: This might be better bundled with the Laplacian inversion code @@ -207,13 +207,13 @@ public: // Full parallel Laplacian operator on scalar field // Laplace_par(f) = Div( b (b dot Grad(f)) ) FieldMetric Laplace_par(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT); - Field3D Laplace_par(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT); + Field3D Laplace_par(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT); // Full Laplacian operator on scalar field FieldMetric Laplace(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& dfdy_boundary_conditions = "free_o3", const std::string& dfdy_dy_region = ""); - Field3D Laplace(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, + Field3D Laplace(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& dfdy_boundary_conditions = "free_o3", const std::string& dfdy_dy_region = ""); @@ -256,10 +256,10 @@ private: class TokamakCoordinates : public Coordinates { public: TokamakCoordinates(Mesh *mesh) : Coordinates(mesh) { - + } private: - + }; */ diff --git a/include/bout/deriv_store.hxx b/include/bout/deriv_store.hxx index b2b7928c2e..83bc34d7c8 100644 --- a/include/bout/deriv_store.hxx +++ b/include/bout/deriv_store.hxx @@ -522,4 +522,13 @@ private: } }; +template +auto& getStore() { + if constexpr (std::is_same::value) { + return DerivativeStore::getInstance(); + } else { + return DerivativeStore::getInstance(); + } +} + #endif diff --git a/include/bout/derivs.hxx b/include/bout/derivs.hxx index 1c360bb9cd..a8d9279378 100644 --- a/include/bout/derivs.hxx +++ b/include/bout/derivs.hxx @@ -82,7 +82,7 @@ Coordinates::FieldMetric DDX(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT, /// If not given, defaults to DIFF_DEFAULT /// @param[in] region What region is expected to be calculated /// If not given, defaults to RGN_NOBNDRY -Field3D DDY(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, +Field3D DDY(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", const std::string& region = "RGN_NOBNDRY"); @@ -410,7 +410,7 @@ Coordinates::FieldMetric VDDX(const Field2D& v, const Field2D& f, /// If not given, defaults to DIFF_DEFAULT /// @param[in] region What region is expected to be calculated /// If not given, defaults to RGN_NOBNDRY -Field3D VDDY(const Field3D& v, const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, +Field3D VDDY(const Field3D& v, const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", const std::string& region = "RGN_NOBNDRY"); @@ -533,7 +533,7 @@ Coordinates::FieldMetric FDDX(const Field2D& v, const Field2D& f, /// If not given, defaults to DIFF_DEFAULT /// @param[in] region What region is expected to be calculated /// If not given, defaults to RGN_NOBNDRY -Field3D FDDY(const Field3D& v, const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, +Field3D FDDY(const Field3D& v, const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", const std::string& region = "RGN_NOBNDRY"); diff --git a/include/bout/difops.hxx b/include/bout/difops.hxx index 18220b63ad..2070cc30d3 100644 --- a/include/bout/difops.hxx +++ b/include/bout/difops.hxx @@ -1,11 +1,11 @@ /*!****************************************************************************** * \file difops.hxx - * + * * Differential operators * * Changelog: * - * 2009-01 Ben Dudson + * 2009-01 Ben Dudson * * Added two optional parameters which can be put in any order * These determine the method to use (DIFF_METHOD) * and CELL_LOC location of the result. @@ -15,7 +15,7 @@ * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu * * Contact: Ben Dudson, bd512@york.ac.uk - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -30,7 +30,7 @@ * * You should have received a copy of the GNU Lesser General Public License * along with BOUT++. If not, see . - * + * *******************************************************************************/ #ifndef BOUT_DIFOPS_H @@ -81,7 +81,7 @@ Field3D Grad_parP(const Field3D& apar, const Field3D& f); * \f[ * v\mathbf{b}_0 \cdot \nabla f * \f] - * + * * * @param[in] v The velocity in y direction * @param[in] f The scalar field to be differentiated @@ -175,7 +175,7 @@ inline Field3D Grad2_par2(const Field3D& f, CELL_LOC outloc, DIFF_METHOD method) /*! * Parallel divergence of diffusive flux, K*Grad_par - * + * * \f[ * \nabla \cdot ( \mathbf{b}_0 kY (\mathbf{b}_0 \cdot \nabla) f ) * \f] @@ -195,12 +195,16 @@ Field3D Div_par_K_Grad_par(const Field3D& kY, const Field2D& f, Field3D Div_par_K_Grad_par(const Field3D& kY, const Field3D& f, CELL_LOC outloc = CELL_DEFAULT); +/// Version with energy flow diagnostic +Field3D Div_par_K_Grad_par_mod(const Field3D& k, const Field3D& f, Field3D& flow_ylow, + bool bndry_flux = true); + /*! * Perpendicular Laplacian operator * * This version only includes terms in X and Z, dropping - * derivatives in Y. This is the inverse operation to - * the Laplacian inversion class. + * derivatives in Y. This is the inverse operation to + * the Laplacian inversion class. * * For the full perpendicular Laplacian, use Laplace_perp */ @@ -212,7 +216,7 @@ FieldPerp Delp2(const FieldPerp& f, CELL_LOC outloc = CELL_DEFAULT, bool useFFT /*! * Perpendicular Laplacian, keeping y derivatives * - * + * */ Coordinates::FieldMetric Laplace_perp(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT, @@ -246,18 +250,18 @@ Field2D Laplace_perpXY(const Field2D& A, const Field2D& f); /*! * Terms of form b0 x Grad(phi) dot Grad(A) - * + * */ Coordinates::FieldMetric b0xGrad_dot_Grad(const Field2D& phi, const Field2D& A, CELL_LOC outloc = CELL_DEFAULT); /*! - * Terms of form + * Terms of form * * \f[ * \mathbf{b}_0 \times \nabla \phi \cdot \nabla A * \f] - * + * * @param[in] phi The scalar potential * @param[in] A The field being advected * @param[in] outloc The cell location where the result is defined. By default the same as A. @@ -291,13 +295,13 @@ constexpr BRACKET_METHOD BRACKET_CTU = BRACKET_METHOD::ctu; * \f[ * [f, g] = (1/B) \mathbf{b}_0 \times \nabla f \cdot \nabla g * \f] - * + * * @param[in] f The potential * @param[in] g The field being advected * @param[in] method The method to use * @param[in] outloc The cell location where the result is defined. Default is the same as g * @param[in] solver Pointer to the time integration solver - * + * */ Coordinates::FieldMetric bracket(const Field2D& f, const Field2D& g, BRACKET_METHOD method = BRACKET_STD, diff --git a/include/bout/field.hxx b/include/bout/field.hxx index 273cc08342..4ad69f5b78 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -137,6 +137,7 @@ public: virtual void calcParallelSlices() {} virtual void splitParallelSlices() {} virtual void clearParallelSlices() {} + virtual size_t numberParallelSlices() const { return 0; } private: /// Labels for the type of coordinate system this field is defined over @@ -189,8 +190,8 @@ inline bool areFieldsCompatible(const Field& field1, const Field& field2) { template inline T emptyFrom(const T& f) { static_assert(bout::utils::is_Field_v, "emptyFrom only works on Fields"); - return T(f.getMesh(), f.getLocation(), {f.getDirectionY(), f.getDirectionZ()}, - f.getRegionID()) + return T(f.getMesh(), f.getLocation(), + DirectionTypes{f.getDirectionY(), f.getDirectionZ()}, f.getRegionID()) .allocate(); } @@ -524,20 +525,27 @@ T pow(BoutReal lhs, const T& rhs, const std::string& rgn = "RGN_ALL") { * result for non-finite numbers * */ +class Field3DParallel; #ifdef FIELD_FUNC #error This macro has already been defined #else -#define FIELD_FUNC(name, func) \ - template > \ - inline T name(const T& f, const std::string& rgn = "RGN_ALL") { \ - \ - /* Check if the input is allocated */ \ - checkData(f); \ - /* Define and allocate the output result */ \ - T result{emptyFrom(f)}; \ - BOUT_FOR(d, result.getRegion(rgn)) { result[d] = func(f[d]); } \ - checkData(result); \ - return result; \ +#define FIELD_FUNC(_name, func) \ + template > \ + inline T _name(const T& f, const std::string& rgn = "RGN_ALL") { \ + /* Check if the input is allocated */ \ + checkData(f); \ + /* Define and allocate the output result */ \ + T result{emptyFrom(f)}; \ + BOUT_FOR(d, result.getRegion(rgn)) { result[d] = func(f[d]); } \ + if constexpr (std::is_base_of_v) { \ + for (int i = 0; i < f.numberParallelSlices(); ++i) { \ + result.yup(i) = func(f.yup(i)); \ + result.ydown(i) = func(f.ydown(i)); \ + } \ + } \ + result.name = std::string(#_name "(") + f.name + std::string(")"); \ + checkData(result); \ + return result; \ } #endif @@ -661,6 +669,8 @@ T copy(const T& f) { return result; } +class Field3DParallel; + /// Apply a floor value \p f to a field \p var. Any value lower than /// the floor is set to the floor. /// @@ -677,7 +687,36 @@ inline T floor(const T& var, BoutReal f, const std::string& rgn = "RGN_ALL") { result[d] = f; } } - + if constexpr (std::is_same_v) { + if (var.hasParallelSlices()) { + for (size_t i = 0; i < result.numberParallelSlices(); ++i) { + if (result.yup(i).isAllocated()) { + BOUT_FOR(d, result.yup(i).getRegion(rgn)) { + if (result.yup(i)[d] < f) { + result.yup(i)[d] = f; + } + } + } else { + if (result.isFci()) { + throw BoutException("Expected parallel slice to be allocated"); + } + } + if (result.ydown(i).isAllocated()) { + BOUT_FOR(d, result.ydown(i).getRegion(rgn)) { + if (result.ydown(i)[d] < f) { + result.ydown(i)[d] = f; + } + } + } else { + if (result.isFci()) { + throw BoutException("Expected parallel slice to be allocated"); + } + } + } + } + } else { + result.clearParallelSlices(); + } return result; } diff --git a/include/bout/field2d.hxx b/include/bout/field2d.hxx index 2af5c519d2..025285b84a 100644 --- a/include/bout/field2d.hxx +++ b/include/bout/field2d.hxx @@ -70,7 +70,7 @@ public: Field2D(Mesh* localmesh = nullptr, CELL_LOC location_in = CELL_CENTRE, DirectionTypes directions_in = {YDirectionType::Standard, ZDirectionType::Average}, - std::optional region = {}); + std::optional regionID = {}); /*! * Copy constructor. After this both fields @@ -274,6 +274,9 @@ public: int size() const override { return nx * ny; } + Field2D& asField3DParallel() { return *this; } + const Field2D& asField3DParallel() const { return *this; } + private: /// Internal data array. Handles allocation/freeing of memory Array data; @@ -286,6 +289,10 @@ private: }; // Non-member overloaded operators +FieldPerp operator+(const Field2D& lhs, const FieldPerp& rhs); +FieldPerp operator-(const Field2D& lhs, const FieldPerp& rhs); +FieldPerp operator*(const Field2D& lhs, const FieldPerp& rhs); +FieldPerp operator/(const Field2D& lhs, const FieldPerp& rhs); Field2D operator+(const Field2D& lhs, const Field2D& rhs); Field2D operator-(const Field2D& lhs, const Field2D& rhs); diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index 7ba6ae9575..68dc84d1a7 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -37,13 +37,16 @@ class Field3D; #include "bout/region.hxx" #include "bout/traits.hxx" +#include #include #include #include +#include #include class Mesh; class Options; +class Field3DParallel; /// Class for 3D X-Y-Z scalar fields /*! @@ -264,6 +267,13 @@ public: #endif } + size_t numberParallelSlices() const override { +#if CHECK > 0 + hasParallelSlices(); +#endif + return yup_fields.size(); + } + /// Check if this field has yup and ydown fields /// Return reference to yup field Field3D& yup(size_t index = 0) { @@ -332,7 +342,7 @@ public: const Region& getValidRegionWithDefault(const std::string& region_name) const; void setRegion(const std::string& region_name) override; void resetRegion() override { regionID.reset(); }; - void resetRegionParallel(); + void resetRegionParallel(bool force = false); void setRegion(size_t id) override { regionID = id; }; void setRegion(std::optional id) override { regionID = id; }; std::optional getRegionID() const override { return regionID; }; @@ -473,6 +483,19 @@ public: Field3D& operator/=(BoutReal rhs); ///@} + Field3D& update_multiplication_inplace(const Field3D& rhs); + Field3D& update_division_inplace(const Field3D& rhs); + Field3D& update_addition_inplace(const Field3D& rhs); + Field3D& update_subtraction_inplace(const Field3D& rhs); + Field3D& update_multiplication_inplace(const Field2D& rhs); + Field3D& update_division_inplace(const Field2D& rhs); + Field3D& update_addition_inplace(const Field2D& rhs); + Field3D& update_subtraction_inplace(const Field2D& rhs); + Field3D& update_multiplication_inplace(BoutReal rhs); + Field3D& update_division_inplace(BoutReal rhs); + Field3D& update_addition_inplace(BoutReal rhs); + Field3D& update_subtraction_inplace(BoutReal rhs); + // FieldData virtual functions FieldType field_type() const override { return FieldType::field3d; } @@ -518,7 +541,10 @@ public: bool areCalcParallelSlicesAllowed() const { return _allowCalcParallelSlices; }; void disallowCalcParallelSlices() { _allowCalcParallelSlices = false; }; -private: + inline Field3DParallel asField3DParallel(); + inline const Field3DParallel asField3DParallel() const; + +protected: /// Array sizes (from fieldmesh). These are valid only if fieldmesh is not null int nx{-1}, ny{-1}, nz{-1}; @@ -581,6 +607,31 @@ Field3D operator-(BoutReal lhs, const Field3D& rhs); Field3D operator*(BoutReal lhs, const Field3D& rhs); Field3D operator/(BoutReal lhs, const Field3D& rhs); +Field3DParallel operator+(const Field3D& lhs, const Field3DParallel& rhs); +Field3DParallel operator-(const Field3D& lhs, const Field3DParallel& rhs); +Field3DParallel operator*(const Field3D& lhs, const Field3DParallel& rhs); +Field3DParallel operator/(const Field3D& lhs, const Field3DParallel& rhs); + +Field3DParallel operator+(const Field3DParallel& lhs, const Field3D& rhs); +Field3DParallel operator-(const Field3DParallel& lhs, const Field3D& rhs); +Field3DParallel operator*(const Field3DParallel& lhs, const Field3D& rhs); +Field3DParallel operator/(const Field3DParallel& lhs, const Field3D& rhs); + +Field3DParallel operator+(const Field3DParallel& lhs, const Field3DParallel& rhs); +Field3DParallel operator-(const Field3DParallel& lhs, const Field3DParallel& rhs); +Field3DParallel operator*(const Field3DParallel& lhs, const Field3DParallel& rhs); +Field3DParallel operator/(const Field3DParallel& lhs, const Field3DParallel& rhs); + +Field3DParallel operator+(BoutReal lhs, const Field3DParallel& rhs); +Field3DParallel operator-(BoutReal lhs, const Field3DParallel& rhs); +Field3DParallel operator*(BoutReal lhs, const Field3DParallel& rhs); +Field3DParallel operator/(BoutReal lhs, const Field3DParallel& rhs); + +Field3DParallel operator+(const Field3DParallel& lhs, BoutReal rhs); +Field3DParallel operator-(const Field3DParallel& lhs, BoutReal rhs); +Field3DParallel operator*(const Field3DParallel& lhs, BoutReal rhs); +Field3DParallel operator/(const Field3DParallel& lhs, BoutReal rhs); + /*! * Unary minus. Returns the negative of given field, * iterates over whole domain including guard/boundary cells. @@ -693,4 +744,150 @@ bool operator==(const Field3D& a, const Field3D& b); /// Output a string describing a Field3D to a stream std::ostream& operator<<(std::ostream& out, const Field3D& value); +inline Field3D copy(const Field3D& f) { + Field3D result{f}; + result.allocate(); + for (size_t i = 0; i < result.numberParallelSlices(); ++i) { + result.yup(i).allocate(); + result.ydown(i).allocate(); + } + return result; +} + +/// Field3DParallel is intended to behave like Field3D, but preserve parallel +/// Fields. +/// Operations on Field3D, like multiplication, exp and floor only work on the +/// "main" field, Field3DParallel will retain the parallel slices. +class Field3DParallel : public Field3D { +public: + template + explicit Field3DParallel(Types... args) : Field3D(std::move(args)...) { + ensureFieldAligned(); + } + Field3DParallel(const Field3D& f) : Field3D(f) { ensureFieldAligned(); } + Field3DParallel(const Field3D& f, bool isRef) : Field3D(f), isRef(isRef) { + ensureFieldAligned(); + } + Field3DParallel(const Field2D& f) : Field3D(f) { ensureFieldAligned(); } + // Explicitly needed, as DirectionTypes is sometimes constructed from a + // brace enclosed list + explicit Field3DParallel(Mesh* localmesh = nullptr, CELL_LOC location_in = CELL_CENTRE, + DirectionTypes directions_in = {YDirectionType::Standard, + ZDirectionType::Standard}, + std::optional regionID = {}) + : Field3D(localmesh, location_in, directions_in, regionID) { + if (isFci()) { + splitParallelSlices(); + } + ensureFieldAligned(); + } + explicit Field3DParallel(Array data, Mesh* localmesh, + CELL_LOC location = CELL_CENTRE, + DirectionTypes directions_in = {YDirectionType::Standard, + ZDirectionType::Standard}) + : Field3D(std::move(data), localmesh, location, directions_in) { + ensureFieldAligned(); + } + explicit Field3DParallel(BoutReal, Mesh* mesh = nullptr); + Field3D& asField3D() { return *this; } + const Field3D& asField3D() const { return *this; } + + Field3DParallel& operator*=(const Field3D&); + Field3DParallel& operator/=(const Field3D&); + Field3DParallel& operator+=(const Field3D&); + Field3DParallel& operator-=(const Field3D&); + Field3DParallel& operator*=(const Field3DParallel&); + Field3DParallel& operator/=(const Field3DParallel&); + Field3DParallel& operator+=(const Field3DParallel&); + Field3DParallel& operator-=(const Field3DParallel&); + Field3DParallel& operator*=(BoutReal); + Field3DParallel& operator/=(BoutReal); + Field3DParallel& operator+=(BoutReal); + Field3DParallel& operator-=(BoutReal); + Field3DParallel& operator=(const Field3D& rhs) { + Field3D::operator=(rhs); + ensureFieldAligned(); + return *this; + } + Field3DParallel& operator=(Field3D&& rhs) { + Field3D::operator=(std::move(rhs)); + ensureFieldAligned(); + return *this; + } + Field3DParallel& operator=(BoutReal); + Field3DParallel& allocate(); + +private: + void ensureFieldAligned(); + bool isRef{false}; +}; + +Field3DParallel Field3D::asField3DParallel() { + if (isAllocated()) { + allocate(); + for (size_t i = 0; i < numberParallelSlices(); ++i) { + if (yup(i).isAllocated()) { + yup(i).allocate(); + } + if (ydown(i).isAllocated()) { + ydown(i).allocate(); + } + } + } + return Field3DParallel(*this, true); +} +const Field3DParallel Field3D::asField3DParallel() const { + return Field3DParallel(*this); +} + +inline Field3D operator+(const Field2D& lhs, const Field3DParallel& rhs) { + return lhs + rhs.asField3D(); +} +inline Field3D operator-(const Field2D& lhs, const Field3DParallel& rhs) { + return lhs + rhs.asField3D(); +} +inline Field3D operator*(const Field2D& lhs, const Field3DParallel& rhs) { + return lhs + rhs.asField3D(); +} +inline Field3D operator/(const Field2D& lhs, const Field3DParallel& rhs) { + return lhs + rhs.asField3D(); +} + +inline Field3D operator+(const Field3DParallel& lhs, const Field2D& rhs) { + return lhs.asField3D() + rhs; +} +inline Field3D operator-(const Field3DParallel& lhs, const Field2D& rhs) { + return lhs.asField3D() - rhs; +} +inline Field3D operator*(const Field3DParallel& lhs, const Field2D& rhs) { + return lhs.asField3D() * rhs; +} +inline Field3D operator/(const Field3DParallel& lhs, const Field2D& rhs) { + return lhs.asField3D() / rhs; +} + +inline Field3DParallel +filledFrom(const Field3DParallel& f, + const std::function& func) { + auto result{emptyFrom(f)}; + if (f.isFci()) { + BOUT_FOR(i, result.getRegion("RGN_NOY")) { result[i] = func(0, i); } + + for (size_t i = 0; i < result.numberParallelSlices(); ++i) { + result.yup(i).allocate(); + BOUT_FOR(d, result.yup(i).getValidRegionWithDefault("RGN_INVALID")) { + result.yup(i)[d] = func(i + 1, d); + } + result.ydown(i).allocate(); + BOUT_FOR(d, result.ydown(i).getValidRegionWithDefault("RGN_INVALID")) { + result.ydown(i)[d] = func(-i - 1, d); + } + } + } else { + BOUT_FOR(i, result.getRegion("RGN_ALL")) { result[i] = func(0, i); } + } + + return result; +} + #endif /* BOUT_FIELD3D_H */ diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 0ec1fbe3ad..027a1ef77b 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -5,29 +5,38 @@ #ifndef BOUT_FV_OPS_H #define BOUT_FV_OPS_H +#include "bout/assert.hxx" +#include "bout/bout_types.hxx" +#include "bout/boutexception.hxx" #include "bout/build_defines.hxx" +#include "bout/coordinates.hxx" +#include "bout/field.hxx" +#include "bout/field2d.hxx" #include "bout/field3d.hxx" #include "bout/globals.hxx" +#include "bout/mesh.hxx" +#include "bout/output_bout_types.hxx" // NOLINT(unused-includes, misc-include-cleaner) +#include "bout/region.hxx" #include "bout/utils.hxx" #include "bout/vector2d.hxx" -#include + +#include namespace FV { /*! * Div ( a Grad_perp(f) ) -- ∇⊥ ( a ⋅ ∇⊥ f) -- Vorticity */ -Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& x); +Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f); [[deprecated("Please use Div_a_Grad_perp instead")]] inline Field3D -Div_a_Laplace_perp(const Field3D& a, const Field3D& x) { - return Div_a_Grad_perp(a, x); +Div_a_Laplace_perp(const Field3D& a, const Field3D& f) { + return Div_a_Grad_perp(a, f); } /*! * Divergence of a parallel diffusion Div( k * Grad_par(f) ) */ -const Field3D Div_par_K_Grad_par(const Field3D& k, const Field3D& f, - bool bndry_flux = true); +Field3D Div_par_K_Grad_par(const Field3D& k, const Field3D& f, bool bndry_flux = true); /*! * 4th-order derivative in Y, using derivatives @@ -49,7 +58,7 @@ const Field3D Div_par_K_Grad_par(const Field3D& k, const Field3D& f, * * No fluxes through domain boundaries */ -const Field3D D4DY4(const Field3D& d, const Field3D& f); +Field3D D4DY4(const Field3D& d, const Field3D& f); /*! * 4th-order dissipation term @@ -67,18 +76,24 @@ const Field3D D4DY4(const Field3D& d, const Field3D& f); * f_2 | f_1 | f_0 | * f_b */ -const Field3D D4DY4_Index(const Field3D& f, bool bndry_flux = true); +Field3D D4DY4_Index(const Field3D& f, bool bndry_flux = true); /*! * Stencil used for Finite Volume calculations * which includes cell face values L and R */ struct Stencil1D { - // Cell centre values - BoutReal c, m, p, mm, pp; - - // Left and right cell face values - BoutReal L, R; + /// Cell centre values + BoutReal c; + BoutReal m; + BoutReal p; + BoutReal mm = BoutNaN; + BoutReal pp = BoutNaN; + + /// Left cell face value + BoutReal L = BoutNaN; + /// Right cell face value + BoutReal R = BoutNaN; }; /*! @@ -93,14 +108,14 @@ struct Upwind { */ struct Fromm { void operator()(Stencil1D& n) { - n.L = n.c - 0.25 * (n.p - n.m); - n.R = n.c + 0.25 * (n.p - n.m); + n.L = n.c - (0.25 * (n.p - n.m)); + n.R = n.c + (0.25 * (n.p - n.m)); } }; /*! * Second order slope limiter method - * + * * Limits slope to minimum absolute value * of left and right gradients. If at a maximum * or minimum slope set to zero, i.e. reverts @@ -110,9 +125,9 @@ struct MinMod { void operator()(Stencil1D& n) { // Choose the gradient within the cell // as the minimum (smoothest) solution - BoutReal slope = _minmod(n.p - n.c, n.c - n.m); - n.L = n.c - 0.5 * slope; - n.R = n.c + 0.5 * slope; + const BoutReal slope = _minmod(n.p - n.c, n.c - n.m); + n.L = n.c - (0.5 * slope); + n.R = n.c + (0.5 * slope); } private: @@ -123,7 +138,7 @@ private: * returns zero, otherwise chooses the value * with the minimum magnitude. */ - BoutReal _minmod(BoutReal a, BoutReal b) { + static BoutReal _minmod(BoutReal a, BoutReal b) { if (a * b <= 0.0) { return 0.0; } @@ -137,25 +152,25 @@ private: /*! * Monotonised Central (MC) second order slope limiter (Van Leer) - * - * Limits the slope based on taking the slope with + * + * Limits the slope based on taking the slope with * the minimum absolute value from central, 2*left and * 2*right. If any of these slopes have different signs * then the slope reverts to zero (i.e. 1st-order upwinding). */ struct MC { void operator()(Stencil1D& n) { - BoutReal slope = minmod(2. * (n.p - n.c), // 2*right difference - 0.5 * (n.p - n.m), // Central difference - 2. * (n.c - n.m)); // 2*left difference - n.L = n.c - 0.5 * slope; - n.R = n.c + 0.5 * slope; + const BoutReal slope = minmod(2. * (n.p - n.c), // 2*right difference + 0.5 * (n.p - n.m), // Central difference + 2. * (n.c - n.m)); // 2*left difference + n.L = n.c - (0.5 * slope); + n.R = n.c + (0.5 * slope); } private: // Return zero if any signs are different // otherwise return the value with the minimum magnitude - BoutReal minmod(BoutReal a, BoutReal b, BoutReal c) { + static BoutReal minmod(BoutReal a, BoutReal b, BoutReal c) { // if any of the signs are different, return zero gradient if ((a * b <= 0.0) || (a * c <= 0.0)) { return 0.0; @@ -166,6 +181,52 @@ private: } }; +/// 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) { + const BoutReal gL = n.c - n.m; + const BoutReal gR = n.p - n.c; + + // r = gL / gR + // Limiter is φ(r) + if (gL * gR < 0) { + // Different signs => Zero gradient + n.L = n.R = n.c; + } else { + const BoutReal sign = SIGN(gL); + const BoutReal abs_gL = fabs(gL); + const BoutReal abs_gR = fabs(gR); + const BoutReal half_slope = + sign * BOUTMAX(BOUTMIN(abs_gL, 0.5 * abs_gR), BOUTMIN(abs_gR, 0.5 * abs_gL)); + n.L = n.c - half_slope; + n.R = n.c + half_slope; + } + } +}; + /*! * Communicate fluxes between processors * Takes values in guard cells, and adds them to cells @@ -189,13 +250,13 @@ void communicateFluxes(Field3D& f); /// /// NB: Uses to/from FieldAligned coordinates template -const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, - const Field3D& wave_speed_in, bool fixflux = true) { +Field3D Div_par(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, + bool fixflux = true) { ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in); - Mesh* mesh = f_in.getMesh(); + const Mesh* mesh = f_in.getMesh(); CellEdges cellboundary; @@ -215,29 +276,17 @@ const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, 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++) { + const bool is_periodic_y = mesh->periodicY(i); + const bool is_first_y = mesh->firstY(i); + const bool is_last_y = mesh->lastY(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; - } + // Only need one guard cell, so no need to communicate fluxes Instead + // calculate in guard cells to get fluxes consistent between processors, but + // don't include the boundary cell. Note that this implies special handling + // of boundaries later + const int ys = (!is_first_y || is_periodic_y) ? mesh->ystart - 1 : mesh->ystart; + const int ye = (!is_last_y || is_periodic_y) ? mesh->yend + 1 : mesh->yend; for (int j = ys; j <= ye; j++) { // Pre-calculate factors which multiply fluxes @@ -246,16 +295,16 @@ const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, 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 = + const BoutReal flux_factor_rc = common_factor / (coord->dy(i, j) * coord->J(i, j)); + const BoutReal flux_factor_rp = common_factor / (coord->dy(i, j + 1) * coord->J(i, j + 1)); // For left cell boundaries common_factor = (coord->J(i, j) + coord->J(i, j - 1)) / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); - BoutReal flux_factor_lc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - BoutReal flux_factor_lm = + const BoutReal flux_factor_lc = common_factor / (coord->dy(i, j) * coord->J(i, j)); + const BoutReal flux_factor_lm = common_factor / (coord->dy(i, j - 1) * coord->J(i, j - 1)); #endif for (int k = mesh->zstart; k <= mesh->zend; k++) { @@ -298,23 +347,23 @@ const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, // Calculate velocity at right boundary (y+1/2) BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j + 1, k)); - BoutReal flux; + BoutReal flux = NAN; - if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) { + if (is_last_y && (j == mesh->yend) && !is_periodic_y) { // Last point in domain - BoutReal bndryval = 0.5 * (s.c + s.p); + const BoutReal bndryval = 0.5 * (s.c + s.p); if (fixflux) { // Use mid-point to be consistent with boundary conditions flux = bndryval * vpar; } else { // Add flux due to difference in boundary values - flux = s.R * vpar + wave_speed(i, j, k) * (s.R - bndryval); + 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)); + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k)); if (vpar > amax) { // Supersonic flow out of this cell @@ -336,20 +385,20 @@ const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, vpar = 0.5 * (v(i, j, k) + v(i, j - 1, k)); - if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) { + if (is_first_y && (j == mesh->ystart) && !is_periodic_y) { // First point in domain - BoutReal bndryval = 0.5 * (s.c + s.m); + const BoutReal bndryval = 0.5 * (s.c + s.m); if (fixflux) { // Use mid-point to be consistent with boundary conditions flux = bndryval * vpar; } else { // Add flux due to difference in boundary values - flux = s.L * vpar - wave_speed(i, j, k) * (s.L - bndryval); + 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)); + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k)); if (vpar < -amax) { // Supersonic out of this cell @@ -374,20 +423,20 @@ const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, * Div ( n * v ) -- Magnetic drifts * * This uses the expression - * + * * Div( A ) = 1/J * d/di ( J * A^i ) - * + * * Hence the input vector should be contravariant * * Note: Uses to/from FieldAligned * */ template -const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { +Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { ASSERT1(n_in.getLocation() == v.getLocation()); ASSERT1_FIELDS_COMPATIBLE(n_in, v.x); - Mesh* mesh = n_in.getMesh(); + const Mesh* mesh = n_in.getMesh(); CellEdges cellboundary; @@ -406,10 +455,10 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { // Calculate velocities - BoutReal vU = 0.25 * (vz[i.zp()] + vz[i]) * (coord->J[i.zp()] + coord->J[i]); - BoutReal vD = 0.25 * (vz[i.zm()] + vz[i]) * (coord->J[i.zm()] + coord->J[i]); - BoutReal vL = 0.25 * (vx[i.xm()] + vx[i]) * (coord->J[i.xm()] + coord->J[i]); - BoutReal vR = 0.25 * (vx[i.xp()] + vx[i]) * (coord->J[i.xp()] + coord->J[i]); + const BoutReal vU = 0.25 * (vz[i.zp()] + vz[i]) * (coord->J[i.zp()] + coord->J[i]); + const BoutReal vD = 0.25 * (vz[i.zm()] + vz[i]) * (coord->J[i.zm()] + coord->J[i]); + const BoutReal vL = 0.25 * (vx[i.xm()] + vx[i]) * (coord->J[i.xm()] + coord->J[i]); + const BoutReal vR = 0.25 * (vx[i.xp()] + vx[i]) * (coord->J[i.xp()] + coord->J[i]); // X direction Stencil1D s; @@ -424,7 +473,7 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { if ((i.x() == mesh->xend) && (mesh->lastX())) { // At right boundary in X if (bndry_flux) { - BoutReal flux; + BoutReal flux = NAN; if (vR > 0.0) { // Flux to boundary flux = vR * s.R; @@ -439,7 +488,7 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { // Not at a boundary if (vR > 0.0) { // Flux out into next cell - BoutReal flux = vR * s.R; + const BoutReal flux = vR * s.R; result[i] += flux / (coord->dx[i] * coord->J[i]); result[i.xp()] -= flux / (coord->dx[i.xp()] * coord->J[i.xp()]); } @@ -451,7 +500,7 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { // At left boundary in X if (bndry_flux) { - BoutReal flux; + BoutReal flux = NAN; if (vL < 0.0) { // Flux to boundary flux = vL * s.L; @@ -465,7 +514,7 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { } else { // Not at a boundary if (vL < 0.0) { - BoutReal flux = vL * s.L; + const BoutReal flux = vL * s.L; result[i] -= flux / (coord->dx[i] * coord->J[i]); result[i.xm()] += flux / (coord->dx[i.xm()] * coord->J[i.xm()]); } @@ -482,12 +531,12 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { cellboundary(s); if (vU > 0.0) { - BoutReal flux = vU * s.R; + const BoutReal flux = vU * s.R; result[i] += flux / (coord->J[i] * coord->dz[i]); result[i.zp()] -= flux / (coord->J[i.zp()] * coord->dz[i.zp()]); } if (vD < 0.0) { - BoutReal flux = vD * s.L; + const BoutReal flux = vD * s.L; result[i] -= flux / (coord->J[i] * coord->dz[i]); result[i.zm()] += flux / (coord->J[i.zm()] * coord->dz[i.zm()]); } @@ -507,13 +556,13 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { // Y velocities on y boundaries - BoutReal vU = 0.25 * (vy[i] + vy[i.yp()]) * (coord->J[i] + coord->J[i.yp()]); - BoutReal vD = 0.25 * (vy[i] + vy[i.ym()]) * (coord->J[i] + coord->J[i.ym()]); + const BoutReal vU = 0.25 * (vy[i] + vy[i.yp()]) * (coord->J[i] + coord->J[i.yp()]); + const BoutReal vD = 0.25 * (vy[i] + vy[i.ym()]) * (coord->J[i] + coord->J[i.ym()]); // n (advected quantity) on y boundaries // Note: Use unshifted n_in variable - BoutReal nU = 0.5 * (n[i] + n[i.yp()]); - BoutReal nD = 0.5 * (n[i] + n[i.ym()]); + const BoutReal nU = 0.5 * (n[i] + n[i.yp()]); + const BoutReal nD = 0.5 * (n[i] + n[i.ym()]); yresult[i] = (nU * vU - nD * vD) / (coord->J[i] * coord->dy[i]); } @@ -524,5 +573,453 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { * X-Z Finite Volume diffusion operator */ Field3D Div_Perp_Lap(const Field3D& a, const Field3D& f, CELL_LOC outloc = CELL_DEFAULT); + +/// Finite volume parallel divergence +/// +/// NOTE: Modified version, applies limiter to velocity and field +/// Performs better (smaller overshoots) than Div_par +/// +/// Preserves the sum of f*J*dx*dy*dz over the domain +/// +/// @param[in] f_in The field being advected. +/// This will be reconstructed at cell faces +/// using the given CellEdges method +/// @param[in] v_in The advection velocity. +/// This will be interpolated to cell boundaries +/// using linear interpolation +/// @param[in] wave_speed_in Local maximum speed of all waves in the system at each +// point in space +/// @param[in] fixflux Fix the flux at the boundary to be the value at the +/// midpoint (for boundary conditions) +/// +/// @param[out] flow_ylow Flow at the lower Y cell boundary +/// Already includes area factor * flux +template +Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, Field3D& flow_ylow, + bool fixflux = true) { + + Coordinates* coord = f_in.getCoordinates(); + ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); + + if (f_in.isFci()) { + // Use mid-point (cell boundary) averages + + ASSERT1(f_in.hasParallelSlices()); + ASSERT1(v_in.hasParallelSlices()); + + const auto& f_up = f_in.yup(); + const auto& f_down = f_in.ydown(); + + const auto& v_up = v_in.yup(); + const auto& v_down = v_in.ydown(); + + Field3D result{emptyFrom(f_in)}; + BOUT_FOR(i, f_in.getRegion("RGN_NOBNDRY")) { + const auto iyp = i.yp(); + const auto iym = i.ym(); + + result[i] = (0.25 * (f_in[i] + f_up[iyp]) * (v_in[i] + v_up[iyp]) + * (coord->J[i] + coord->J.yup()[iyp]) + / (sqrt(coord->g_22[i]) + sqrt(coord->g_22.yup()[iyp])) + - 0.25 * (f_in[i] + f_down[iym]) * (v_in[i] + v_down[iym]) + * (coord->J[i] + coord->J.ydown()[iym]) + / (sqrt(coord->g_22[i]) + sqrt(coord->g_22.ydown()[iym]))) + / (coord->dy[i] * coord->J[i]); + } + return result; + } + ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in); + + const Mesh* mesh = f_in.getMesh(); + + CellEdges cellboundary; + + ASSERT2(f_in.getDirectionY() == v_in.getDirectionY()); + ASSERT2(f_in.getDirectionY() == wave_speed_in.getDirectionY()); + const bool are_unaligned = + ((f_in.getDirectionY() == YDirectionType::Standard) + and (v_in.getDirectionY() == YDirectionType::Standard) + and (wave_speed_in.getDirectionY() == YDirectionType::Standard)); + + const Field3D f = are_unaligned ? toFieldAligned(f_in, "RGN_NOX") : f_in; + const Field3D v = are_unaligned ? toFieldAligned(v_in, "RGN_NOX") : v_in; + const Field3D wave_speed = + are_unaligned ? toFieldAligned(wave_speed_in, "RGN_NOX") : wave_speed_in; + + Field3D result{zeroFrom(f)}; + flow_ylow = zeroFrom(f); + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + const bool is_periodic_y = mesh->periodicY(i); + const bool is_first_y = mesh->firstY(i); + const bool is_last_y = mesh->lastY(i); + + // Only need one guard cell, so no need to communicate fluxes Instead + // calculate in guard cells to get fluxes consistent between processors, but + // don't include the boundary cell. Note that this implies special handling + // of boundaries later + const int ys = (!is_first_y || is_periodic_y) ? mesh->ystart - 1 : mesh->ystart; + const int ye = (!is_last_y || is_periodic_y) ? mesh->yend + 1 : mesh->yend; + + for (int j = ys; j <= ye; j++) { + // Pre-calculate factors which multiply fluxes +#if not(BOUT_USE_METRIC_3D) + // For right cell boundaries + const BoutReal common_factor_r = + (coord->J(i, j) + coord->J(i, j + 1)) + / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); + + const BoutReal flux_factor_rc = + common_factor_r / (coord->dy(i, j) * coord->J(i, j)); + const BoutReal flux_factor_rp = + common_factor_r / (coord->dy(i, j + 1) * coord->J(i, j + 1)); + + const BoutReal area_rp = + common_factor_r * coord->dx(i, j + 1) * coord->dz(i, j + 1); + + // For left cell boundaries + const BoutReal common_factor_l = + (coord->J(i, j) + coord->J(i, j - 1)) + / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); + + const BoutReal flux_factor_lc = + common_factor_l / (coord->dy(i, j) * coord->J(i, j)); + const BoutReal flux_factor_lm = + common_factor_l / (coord->dy(i, j - 1) * coord->J(i, j - 1)); + + const BoutReal area_lc = common_factor_l * 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 + const BoutReal common_factor_r = + (coord->J(i, j, k) + coord->J(i, j + 1, k)) + / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); + + const BoutReal flux_factor_rc = + common_factor_r / (coord->dy(i, j, k) * coord->J(i, j, k)); + const BoutReal flux_factor_rp = + common_factor_r / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); + + const BoutReal area_rp = + common_factor_r * coord->dx(i, j + 1, k) * coord->dz(i, j + 1, k); + + // For left cell boundaries + const BoutReal common_factor_l = + (coord->J(i, j, k) + coord->J(i, j - 1, k)) + / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); + + const BoutReal flux_factor_lc = + common_factor_l / (coord->dy(i, j, k) * coord->J(i, j, k)); + const BoutReal flux_factor_lm = + common_factor_l / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); + + const BoutReal area_lc = + common_factor_l * 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 + // TODO(peter): We can remove this #ifdef guard after switching to C++20 +#if __cpp_designated_initializers >= 201707L + Stencil1D s{.c = f(i, j, k), .m = f(i, j - 1, k), .p = f(i, j + 1, k)}; +#else + Stencil1D s{f(i, j, k), f(i, j - 1, k), f(i, j + 1, k), BoutNaN, + BoutNaN, BoutNaN, BoutNaN}; +#endif + cellboundary(s); // Calculate s.R and s.L + + //////////////////////////////////////////// + // Reconstruct v at the cell faces + // TODO(peter): We can remove this #ifdef guard after switching to C++20 +#if __cpp_designated_initializers >= 201707L + Stencil1D sv{.c = v(i, j, k), .m = v(i, j - 1, k), .p = v(i, j + 1, k)}; +#else + Stencil1D sv{v(i, j, k), v(i, j - 1, k), v(i, j + 1, k), BoutNaN, + BoutNaN, BoutNaN, BoutNaN}; +#endif + cellboundary(sv); // Calculate sv.R and sv.L + + //////////////////////////////////////////// + // Right boundary + + BoutReal flux = BoutNaN; + + if (is_last_y && (j == mesh->yend) && !is_periodic_y) { + // Last point in domain + + // Calculate velocity at right boundary (y+1/2) + const BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j + 1, k)); + + const BoutReal bndryval = 0.5 * (s.c + s.p); + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = bndryval * vpar; + } else { + // Add flux due to difference in boundary values + flux = (s.R * vpar) + (wave_speed(i, j, k) * (s.R - bndryval)); + } + + } else { + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k), + fabs(v(i, j, k)), fabs(v(i, j + 1, k))); + + flux = s.R * 0.5 * (sv.R + amax); + } + + result(i, j, k) += flux * flux_factor_rc; + result(i, j + 1, k) -= flux * flux_factor_rp; + + flow_ylow(i, j + 1, k) += flux * area_rp; + + //////////////////////////////////////////// + // Calculate at left boundary + + if (is_first_y && (j == mesh->ystart) && !is_periodic_y) { + // First point in domain + const BoutReal bndryval = 0.5 * (s.c + s.m); + const BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j - 1, k)); + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = bndryval * vpar; + } else { + // Add flux due to difference in boundary values + flux = (s.L * vpar) - (wave_speed(i, j, k) * (s.L - bndryval)); + } + } else { + + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k), + fabs(v(i, j, k)), fabs(v(i, j - 1, k))); + + flux = s.L * 0.5 * (sv.L - amax); + } + + result(i, j, k) -= flux * flux_factor_lc; + result(i, j - 1, k) += flux * flux_factor_lm; + + flow_ylow(i, j, k) += flux * area_lc; + } + } + } + if (are_unaligned) { + flow_ylow = fromFieldAligned(flow_ylow, "RGN_NOBNDRY"); + } + return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; +} + +/// This operator calculates Div_par(f v v) +/// It is used primarily (only?) in the parallel momentum equation. +/// +/// This operator is used rather than Div(f fv) so that the values of +/// f and v are consistent with other advection equations: The product +/// fv is not interpolated to cell boundaries. +template +Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true) { + ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); + const Mesh* mesh = f_in.getMesh(); + const Coordinates* coord = f_in.getCoordinates(); + CellEdges cellboundary; + + if (f_in.isFci()) { + // FCI version, using yup/down fields + ASSERT1(f_in.hasParallelSlices()); + ASSERT1(v_in.hasParallelSlices()); + + const auto& B = coord->Bxy; + const auto& B_up = coord->Bxy.yup(); + const auto& B_down = coord->Bxy.ydown(); + + const auto& f_up = f_in.yup(); + const auto& f_down = f_in.ydown(); + + const auto& v_up = v_in.yup(); + const auto& v_down = v_in.ydown(); + + const auto& g_22 = coord->g_22; + const auto& dy = coord->dy; + + Field3D result{emptyFrom(f_in)}; + BOUT_FOR(i, f_in.getRegion("RGN_NOBNDRY")) { + const auto iyp = i.yp(); + const auto iym = i.ym(); + + // Maximum local wave speed + const BoutReal amax = + BOUTMAX(wave_speed_in[i], fabs(v_in[i]), fabs(v_up[iyp]), fabs(v_down[iym])); + + const BoutReal term = (f_up[iyp] * v_up[iyp] * v_up[iyp] / B_up[iyp]) + - (f_down[iym] * v_down[iym] * v_down[iym] / B_down[iym]); + + // Penalty terms. This implementation is very dissipative. + BoutReal penalty = + (amax * (f_in[i] * v_in[i] - f_up[iyp] * v_up[iyp]) / (B[i] + B_up[iyp])) + + (amax * (f_in[i] * v_in[i] - f_down[iym] * v_down[iym]) + / (B[i] + B_down[iym])); + + if (fabs(penalty) > fabs(term) and penalty * v_in[i] > 0) { + if (term * penalty > 0) { + penalty = term; + } else { + penalty = -term; + } + } + + result[i] = B[i] * (term + penalty) / (2 * dy[i] * sqrt(g_22[i])); + +#if CHECK > 0 + if (!std::isfinite(result[i])) { + throw BoutException("Non-finite value in Div_par_fvv at {}\n" + "fup {} vup {} fdown {} vdown {} amax {}\n", + "B {} Bup {} Bdown {} dy {} sqrt(g_22} {}", i, f_up[i], + v_up[i], f_down[i], v_down[i], amax, B[i], B_up[i], B_down[i], + dy[i], sqrt(g_22[i])); + } +#endif + } + return result; + } + + ASSERT1(areFieldsCompatible(f_in, wave_speed_in)); + + /// Ensure that f, v and wave_speed are field aligned + Field3D f = toFieldAligned(f_in, "RGN_NOX"); + Field3D v = toFieldAligned(v_in, "RGN_NOX"); + Field3D wave_speed = toFieldAligned(wave_speed_in, "RGN_NOX"); + + Field3D result{zeroFrom(f)}; + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + const bool is_periodic_y = mesh->periodicY(i); + const bool is_first_y = mesh->firstY(i); + const bool is_last_y = mesh->lastY(i); + + // Only need one guard cell, so no need to communicate fluxes Instead + // calculate in guard cells to get fluxes consistent between processors, but + // don't include the boundary cell. Note that this implies special handling + // of boundaries later + const int ys = (!is_first_y || is_periodic_y) ? mesh->ystart - 1 : mesh->ystart; + const int ye = (!is_last_y || is_periodic_y) ? mesh->yend + 1 : mesh->yend; + + for (int j = ys; j <= ye; j++) { + // Pre-calculate factors which multiply fluxes + + for (int k = 0; k < mesh->LocalNz; k++) { + // For right cell boundaries + const BoutReal common_factor_r = + (coord->J(i, j, k) + coord->J(i, j + 1, k)) + / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); + + const BoutReal flux_factor_rc = + common_factor_r / (coord->dy(i, j, k) * coord->J(i, j, k)); + const BoutReal flux_factor_rp = + common_factor_r / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); + + // For left cell boundaries + const BoutReal common_factor_l = + (coord->J(i, j, k) + coord->J(i, j - 1, k)) + / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); + + const BoutReal flux_factor_lc = + common_factor_l / (coord->dy(i, j, k) * coord->J(i, j, k)); + const BoutReal flux_factor_lm = + common_factor_l / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); + + //////////////////////////////////////////// + // Reconstruct f at the cell faces + // This calculates s.R and s.L for the Right and Left + // face values on this cell + + // Reconstruct f at the cell faces +#if __cpp_designated_initializers >= 201707L + Stencil1D s{.c = f(i, j, k), .m = f(i, j - 1, k), .p = f(i, j + 1, k)}; +#else + Stencil1D s{f(i, j, k), f(i, j - 1, k), f(i, j + 1, k), BoutNaN, + BoutNaN, BoutNaN, BoutNaN}; +#endif + cellboundary(s); // Calculate s.R and s.L + + //////////////////////////////////////////// + // Reconstruct v at the cell faces + // TODO(peter): We can remove this #ifdef guard after switching to C++20 +#if __cpp_designated_initializers >= 201707L + Stencil1D sv{.c = v(i, j, k), .m = v(i, j - 1, k), .p = v(i, j + 1, k)}; +#else + Stencil1D sv{v(i, j, k), v(i, j - 1, k), v(i, j + 1, k), BoutNaN, + BoutNaN, BoutNaN, BoutNaN}; +#endif + cellboundary(sv); + + //////////////////////////////////////////// + // Right boundary + + // Calculate velocity at right boundary (y+1/2) + const BoutReal v_mid_r = 0.5 * (sv.c + sv.p); + // And mid-point density at right boundary + const BoutReal n_mid_r = 0.5 * (s.c + s.p); + BoutReal flux = NAN; + + if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) { + // Last point in domain + + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = n_mid_r * v_mid_r * v_mid_r; + } else { + // Add flux due to difference in boundary values + flux = (s.R * sv.R * sv.R) // Use right cell edge values + + (BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.p)) * n_mid_r + * (sv.R - v_mid_r)); // Damp differences in velocity, not flux + } + } else { + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k), + fabs(sv.c), fabs(sv.p)); + + flux = s.R * 0.5 * (sv.R + amax) * sv.R; + } + + result(i, j, k) += flux * flux_factor_rc; + result(i, j + 1, k) -= flux * flux_factor_rp; + + //////////////////////////////////////////// + // Calculate at left boundary + + const BoutReal v_mid_l = 0.5 * (sv.c + sv.m); + const BoutReal n_mid_l = 0.5 * (s.c + s.m); + + if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) { + // First point in domain + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = n_mid_l * v_mid_l * v_mid_l; + } else { + // Add flux due to difference in boundary values + flux = (s.L * sv.L * sv.L) + - (BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.m)) * n_mid_l + * (sv.L - v_mid_l)); + } + } else { + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k), + fabs(sv.c), fabs(sv.m)); + + flux = s.L * 0.5 * (sv.L - amax) * sv.L; + } + + result(i, j, k) -= flux * flux_factor_lc; + result(i, j - 1, k) += flux * flux_factor_lm; + } + } + } + return fromFieldAligned(result, "RGN_NOBNDRY"); +} } // namespace FV #endif // BOUT_FV_OPS_H diff --git a/include/bout/index_derivs.hxx b/include/bout/index_derivs.hxx index ccce9a7f5e..a8d1bbb1d3 100644 --- a/include/bout/index_derivs.hxx +++ b/include/bout/index_derivs.hxx @@ -148,7 +148,7 @@ struct registerMethod { // removed and we can use nGuard directly in the template statement. const int nGuards = method.meta.nGuards; - auto& derivativeRegister = DerivativeStore::getInstance(); + auto& derivativeRegister = getStore(); switch (method.meta.derivType) { case (DERIV::Standard): diff --git a/include/bout/index_derivs_interface.hxx b/include/bout/index_derivs_interface.hxx index 242492c3f8..cbb1e5ac53 100644 --- a/include/bout/index_derivs_interface.hxx +++ b/include/bout/index_derivs_interface.hxx @@ -29,6 +29,8 @@ #ifndef __INDEX_DERIVS_INTERFACE_HXX__ #define __INDEX_DERIVS_INTERFACE_HXX__ +#include "bout/boutexception.hxx" +#include "bout/field3d.hxx" #include "bout/traits.hxx" #include #include @@ -148,8 +150,8 @@ T standardDerivative(const T& f, CELL_LOC outloc, const std::string& method, } // Lookup the method - auto derivativeMethod = DerivativeStore::getInstance().getStandardDerivative( - method, direction, stagger, derivType); + auto derivativeMethod = + getStore().getStandardDerivative(method, direction, stagger, derivType); // Create the result field T result{emptyFrom(f).setLocation(outloc)}; @@ -176,6 +178,12 @@ T DDX(const T& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "D return standardDerivative(f, outloc, method, region); } +inline Field3D DDX(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT, + const std::string& method = "DEFAULT", + const std::string& region = "RGN_NOBNDRY") { + return DDX(f.asField3D(), outloc, method, region); +} + template T D2DX2(const T& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", const std::string& region = "RGN_NOBNDRY") { @@ -198,17 +206,26 @@ template T DDY(const T& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", const std::string& region = "RGN_NOBNDRY") { - if (f.hasParallelSlices()) { + if (f.isFci() or f.hasParallelSlices()) { ASSERT1(f.getDirectionY() == YDirectionType::Standard); + if (!f.hasParallelSlices()) { + throw BoutException( + "parallel slices needed for parallel derivatives. Make sure to communicate and " + "apply parallel boundary conditions before calling derivative"); + } return standardDerivative(f, outloc, method, region); - } else { - const bool is_unaligned = (f.getDirectionY() == YDirectionType::Standard); - const T f_aligned = is_unaligned ? toFieldAligned(f, "RGN_NOX") : f; - T result = standardDerivative(f_aligned, outloc, - method, region); - return is_unaligned ? fromFieldAligned(result, region) : result; } + const bool is_unaligned = (f.getDirectionY() == YDirectionType::Standard); + const T f_aligned = is_unaligned ? toFieldAligned(f, "RGN_NOX") : f; + T result = standardDerivative(f_aligned, outloc, + method, region); + return is_unaligned ? fromFieldAligned(result, region) : result; +} +inline Field3D DDY(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT, + const std::string& method = "DEFAULT", + const std::string& region = "RGN_NOBNDRY") { + return DDY(f.asField3D(), outloc, method, region); } template @@ -253,6 +270,12 @@ T DDZ(const T& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "D return standardDerivative(f, outloc, method, region); } +inline Field3D DDZ(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT, + const std::string& method = "DEFAULT", + const std::string& region = "RGN_NOBNDRY") { + return DDZ(f.asField3D(), outloc, method, region); +} + template T D2DZ2(const T& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", const std::string& region = "RGN_NOBNDRY") { @@ -356,6 +379,11 @@ T VDDY(const T& vel, const T& f, CELL_LOC outloc = CELL_DEFAULT, return are_unaligned ? fromFieldAligned(result, region) : result; } } +inline Field3D VDDY(const Field3D& v, const Field3DParallel& f, + CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", + const std::string& region = "RGN_NOBNDRY") { + return VDDY(v, f.asField3D(), outloc, method, region); +} template T FDDY(const T& vel, const T& f, CELL_LOC outloc = CELL_DEFAULT, @@ -380,6 +408,11 @@ T FDDY(const T& vel, const T& f, CELL_LOC outloc = CELL_DEFAULT, return are_unaligned ? fromFieldAligned(result, region) : result; } } +inline Field3D FDDY(const Field3D& v, const Field3DParallel& f, + CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", + const std::string& region = "RGN_NOBNDRY") { + return FDDY(v, f.asField3D(), outloc, method, region); +} ////////////// Z DERIVATIVE ///////////////// diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx index 49dea67743..0b425de209 100644 --- a/include/bout/paralleltransform.hxx +++ b/include/bout/paralleltransform.hxx @@ -90,6 +90,10 @@ public: /// require a twist-shift at branch cuts on closed field lines? virtual bool requiresTwistShift(bool twist_shift_enabled, YDirectionType ytype) = 0; + /// Can be implemented to load parallel metrics + /// Needed by FCI + virtual void loadParallelMetrics([[maybe_unused]] Coordinates* coords) {} + protected: /// This method should be called in the constructor to check that if the grid /// has a 'parallel_transform' variable, it has the correct value diff --git a/src/field/field2d.cxx b/src/field/field2d.cxx index c2d35162b3..25e8ca1209 100644 --- a/src/field/field2d.cxx +++ b/src/field/field2d.cxx @@ -44,7 +44,7 @@ #include Field2D::Field2D(Mesh* localmesh, CELL_LOC location_in, DirectionTypes directions_in, - std::optional UNUSED(regionID)) + [[maybe_unused]] std::optional regionID) : Field(localmesh, location_in, directions_in) { if (fieldmesh) { diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 9154f78a4f..38b2682606 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -57,7 +57,7 @@ /// Constructor Field3D::Field3D(Mesh* localmesh, CELL_LOC location_in, DirectionTypes directions_in, std::optional regionID) - : Field(localmesh, location_in, directions_in), regionID{regionID} { + : Field(localmesh, location_in, directions_in), regionID(regionID) { #if BOUT_USE_TRACK name = ""; #endif @@ -96,6 +96,19 @@ Field3D::Field3D(const BoutReal val, Mesh* localmesh) : Field3D(localmesh) { *this = val; } +Field3DParallel::Field3DParallel(const BoutReal val, Mesh* localmesh) + : Field3D(localmesh) { + + *this = val; + if (this->isFci()) { + splitParallelSlices(); + for (size_t i = 0; i < numberParallelSlices(); ++i) { + yup(i) = val; + ydown(i) = val; + } + } +} + Field3D::Field3D(Array data_in, Mesh* localmesh, CELL_LOC datalocation, DirectionTypes directions_in) : Field(localmesh, datalocation, directions_in), data(std::move(data_in)) { @@ -341,6 +354,29 @@ Field3D& Field3D::operator=(const BoutReal val) { allocate(); + BOUT_FOR(i, getRegion("RGN_ALL")) { (*this)[i] = val; } + this->name = "BR"; + + return *this; +} + +Field3DParallel& Field3DParallel::operator=(const BoutReal val) { + TRACE("Field3DParallel = BoutReal"); + track(val, "operator="); + + if (isFci()) { + if (!hasParallelSlices()) { + splitParallelSlices(); + } + for (size_t i = 0; i < numberParallelSlices(); ++i) { + yup(i) = val; + ydown(i) = val; + } + } + resetRegion(); + + allocate(); + BOUT_FOR(i, getRegion("RGN_ALL")) { (*this)[i] = val; } return *this; @@ -847,8 +883,8 @@ void Field3D::setRegion(const std::string& region_name) { regionID = fieldmesh->getRegionID(region_name); } -void Field3D::resetRegionParallel() { - if (isFci()) { +void Field3D::resetRegionParallel(const bool force) { + if (force or isFci()) { for (int i = 0; i < fieldmesh->ystart; ++i) { yup_fields[i].setRegion(fmt::format("RGN_YPAR_{:+d}", i + 1)); ydown_fields[i].setRegion(fmt::format("RGN_YPAR_{:+d}", -i - 1)); @@ -895,6 +931,7 @@ void Field3D::_track(const T& change, std::string operation) { template void Field3D::_track>(const Field3D&, std::string); +template void Field3D::_track(const Field3DParallel&, std::string); template void Field3D::_track(const Field2D&, std::string); template void Field3D::_track<>(const FieldPerp&, std::string); @@ -915,3 +952,35 @@ void Field3D::_track(const BoutReal& change, std::string operation) { {"trace", trace}, }); } + +void Field3DParallel::ensureFieldAligned() { + if (isFci()) { + ASSERT2(hasParallelSlices()); + if (fieldmesh != nullptr) { + for (int i = 0; i < fieldmesh->ystart; ++i) { + ASSERT2(yup_fields[i].getRegionID().has_value()); + ASSERT2(ydown_fields[i].getRegionID().has_value()); + } + } + if (isAllocated()) { + for (int i = 0; i < fieldmesh->ystart; ++i) { + ASSERT2(yup_fields[i].isAllocated()); + ASSERT2(ydown_fields[i].isAllocated()); + } + } + } +} + +Field3DParallel& Field3DParallel::allocate() { + Field3D::allocate(); + if (isFci()) { + ASSERT2(hasParallelSlices()); + if (fieldmesh != nullptr) { + for (int i = 0; i < fieldmesh->ystart; ++i) { + yup_fields[i].allocate(); + ydown_fields[i].allocate(); + } + } + } + return *this; +} diff --git a/src/field/gen_fieldops.jinja b/src/field/gen_fieldops.jinja index 5476242d50..ca171eaf8d 100644 --- a/src/field/gen_fieldops.jinja +++ b/src/field/gen_fieldops.jinja @@ -8,15 +8,31 @@ checkData({{lhs.name}}); checkData({{rhs.name}}); - {% if out == "Field3D" %} - {% if lhs == rhs == "Field3D" %} + {% if out.region_type == "3D" %} + {% if lhs.region_type == rhs.region_type == "3D" %} {{out.name}}.setRegion({{lhs.name}}.getMesh()->getCommonRegion({{lhs.name}}.getRegionID(), {{rhs.name}}.getRegionID())); - {% elif lhs == "Field3D" %} + {% elif lhs.region_type == "3D" %} {{out.name}}.setRegion({{lhs.name}}.getRegionID()); - {% elif rhs == "Field3D" %} + {% elif rhs.region_type == "3D" %} {{out.name}}.setRegion({{rhs.name}}.getRegionID()); {% endif %} + {% if out == "Field3DParallel" %} + if ({{out.name}}.isFci()) { + {{ lhs.assertParallelSlices }} + {{ rhs.assertParallelSlices }} + {{out.name}}.splitParallelSlices(); + {% if lhs.region_type == "3D" %} + for (size_t i{0} ; i < {{lhs.name}}.numberParallelSlices() ; ++i) { + {% else %} + for (size_t i{0} ; i < {{rhs.name}}.numberParallelSlices() ; ++i) { + {% endif %} + {{out.name}}.yup(i) = {{lhs.yup}} {{operator}} {{rhs.yup}}; + {{out.name}}.ydown(i) = {{lhs.ydown}} {{operator}} {{rhs.ydown}}; + } + {{ out.assertParallelSlices }} + } + {% endif %} {% endif %} {% if (out == "Field3D") and ((lhs == "Field2D") or (rhs =="Field2D")) %} @@ -27,17 +43,17 @@ {% else %} {{region_loop}}({{index_var}}, {{rhs.name}}.getRegion({{region_name}})) { {% endif %} - const auto {{mixed_base_ind}} = localmesh->ind2Dto3D({{index_var}}); - {% if (operator == "/") and (rhs == "Field2D") %} + const auto {{mixed_base_ind}} = localmesh->ind2Dto3D({{index_var}}); + {% if (operator == "/") and (rhs == "Field2D") %} const auto tmp = 1.0 / {{rhs.mixed_index}}; - for (int {{jz_var}} = 0; {{jz_var}} < localmesh->LocalNz; ++{{jz_var}}){ - {{out.mixed_index}} = {{lhs.mixed_index}} * tmp; + for (int {{jz_var}} = 0; {{jz_var}} < localmesh->LocalNz; ++{{jz_var}}){ + {{out.mixed_index}} = {{lhs.mixed_index}} * tmp; {% else %} - for (int {{jz_var}} = 0; {{jz_var}} < localmesh->LocalNz; ++{{jz_var}}){ - {{out.mixed_index}} = {{lhs.mixed_index}} {{operator}} {{rhs.mixed_index}}; + for (int {{jz_var}} = 0; {{jz_var}} < localmesh->LocalNz; ++{{jz_var}}){ + {{out.mixed_index}} = {{lhs.mixed_index}} {{operator}} {{rhs.mixed_index}}; {% endif %} - } - } + } + } {% elif out == "FieldPerp" and (lhs == "Field2D" or lhs == "Field3D" or rhs == "Field2D" or rhs == "Field3D")%} Mesh *localmesh = {{lhs.name if lhs.field_type != "BoutReal" else rhs.name}}.getMesh(); @@ -45,11 +61,11 @@ int yind = {{lhs.name if lhs == "FieldPerp" else rhs.name}}.getIndex(); const auto {{mixed_base_ind}} = localmesh->indPerpto3D({{index_var}}, yind); {% if lhs != "FieldPerp" %} - {{out.index}} = {{lhs.base_index}} {{operator}} {{rhs.index}}; + {{out.index}} = {{lhs.base_index}} {{operator}} {{rhs.index}}; {% else %} - {{out.index}} = {{lhs.index}} {{operator}} {{rhs.base_index}}; + {{out.index}} = {{lhs.index}} {{operator}} {{rhs.base_index}}; {% endif %} - } + } {% elif (operator == "/") and (rhs == "BoutReal") %} const auto tmp = 1.0 / {{rhs.index}}; {{region_loop}}({{index_var}}, {{out.name}}.getValidRegionWithDefault({{region_name}})) { @@ -57,51 +73,137 @@ } {% else %} {{region_loop}}({{index_var}}, {{out.name}}.getValidRegionWithDefault({{region_name}})) { - {{out.index}} = {{lhs.index}} {{operator}} {{rhs.index}}; - } + {{out.index}} = {{lhs.index}} {{operator}} {{rhs.index}}; + } {% endif %} checkData({{out.name}}); return {{out.name}}; } +{% if out.field_type == lhs.field_type and lhs == "Field3D" %} +// Provide the C++ operator to update {{lhs}} by {{operator_name}} with {{rhs}} +{{lhs}} &{{lhs}}::update_{{operator_name}}_inplace(const {{rhs.passByReference}}) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + {% if lhs != "BoutReal" and rhs != "BoutReal" %} + ASSERT1_FIELDS_COMPATIBLE(*this, rhs); + {% endif %} + + // Delete existing parallel slices. We don't update parallel slices, so any + // that currently exist will be incorrect. + clearParallelSlices(); + checkData(*this); + checkData({{rhs.name}}); + + {% if lhs.region_type == rhs.region_type == "3D" %} + regionID = fieldmesh->getCommonRegion(regionID, {{rhs.name}}.getRegionID()); + {% endif %} + + + {% if (lhs == "Field3D") and (rhs =="Field2D") %} + {{region_loop}}({{index_var}}, {{rhs.name}}.getRegion({{region_name}})) { + const auto {{mixed_base_ind}} = fieldmesh->ind2Dto3D({{index_var}}); + {% if (operator == "/") and (rhs == "Field2D") %} + const auto tmp = 1.0 / {{rhs.mixed_index}}; + for (int {{jz_var}} = 0; {{jz_var}} < fieldmesh->LocalNz; ++{{jz_var}}){ + (*this)[{{mixed_base_ind}} + {{jz_var}}] *= tmp; + {% else %} + for (int {{jz_var}} = 0; {{jz_var}} < fieldmesh->LocalNz; ++{{jz_var}}){ + (*this)[{{mixed_base_ind}} + {{jz_var}}] {{operator}}= {{rhs.index}}; + {% endif %} + } + } + {% elif rhs == "FieldPerp" and (lhs == "Field3D" or lhs == "Field2D")%} + Mesh *localmesh = this->getMesh(); + + {{region_loop}}({{index_var}}, {{rhs.name}}.getRegion({{region_name}})) { + int yind = {{rhs.name}}.getIndex(); + const auto {{mixed_base_ind}} = localmesh->indPerpto3D({{index_var}}, yind); + (*this)[{{base_ind_var}}] {{operator}}= {{rhs.index}}; + } + {% elif (operator == "/") and (lhs == "Field3D" or lhs == "Field2D") and (rhs =="BoutReal") %} + const auto tmp = 1.0 / {{rhs.index}}; + {{region_loop}}({{index_var}}, this->getRegion({{region_name}})) { + (*this)[{{index_var}}] *= tmp; + } + {% else %} + {{region_loop}}({{index_var}}, this->getRegion({{region_name}})) { + (*this)[{{index_var}}] {{operator}}= {{rhs.index}}; + } + {% endif %} + + {% if lhs.region_type == "3D" %} + track(rhs, "operator{{operator}}="); + {% endif %} +#if BOUT_USE_TRACK + name = fmt::format("{:s} {{operator}}= {:s}", this->name, {{'"BR"' if rhs == "BoutReal" else rhs.name + ".name"}}); +#endif + + checkData(*this); + + return *this; +} +{% endif %} + + {% if out.field_type == lhs.field_type %} // Provide the C++ operator to update {{lhs}} by {{operator_name}} with {{rhs}} {{lhs}} &{{lhs}}::operator{{operator}}=(const {{rhs.passByReference}}) { // only if data is unique we update the field // otherwise just call the non-inplace version +{% if lhs == "Field3DParallel" %} + if (data.unique() or isRef) { +{% else %} if (data.unique()) { +{% endif %} {% if lhs != "BoutReal" and rhs != "BoutReal" %} ASSERT1_FIELDS_COMPATIBLE(*this, rhs); {% endif %} - {% if (lhs == "Field3D") %} - // Delete existing parallel slices. We don't copy parallel slices, so any + {% if lhs == "Field3D" %} + // Delete existing parallel slices. We don't update parallel slices, so any // that currently exist will be incorrect. clearParallelSlices(); - + {% endif %} + {% if lhs == "Field3DParallel" and (rhs.region_type == "3D" or rhs == "BoutReal") %} + if (this->isFci()) { + if (isRef) { + for (size_t i{0} ; i < yup_fields.size() ; ++i) { + yup(i).update_{{operator_name}}_inplace({{rhs.name}}{% if rhs == "Field3D" %}.yup(i){% endif %}); + ydown(i).update_{{operator_name}}_inplace({{rhs.name}}{% if rhs == "Field3D" %}.ydown(i){% endif %}); + } + } else { + for (size_t i{0} ; i < yup_fields.size() ; ++i) { + yup(i) {{operator}}= {{rhs.name}}{% if rhs == "Field3D" %}.yup(i){% endif %}; + ydown(i) {{operator}}= {{rhs.name}}{% if rhs == "Field3D" %}.ydown(i){% endif %}; + } + } + } else { + clearParallelSlices(); + } {% endif %} checkData(*this); checkData({{rhs.name}}); - {% if lhs == rhs == "Field3D" %} - regionID = fieldmesh->getCommonRegion(regionID, {{rhs.name}}.regionID); + {% if lhs.region_type == rhs.region_type == "3D" %} + regionID = fieldmesh->getCommonRegion(regionID, {{rhs.name}}.getRegionID()); {% endif %} {% if (lhs == "Field3D") and (rhs =="Field2D") %} {{region_loop}}({{index_var}}, {{rhs.name}}.getRegion({{region_name}})) { - const auto {{mixed_base_ind}} = fieldmesh->ind2Dto3D({{index_var}}); - {% if (operator == "/") and (rhs == "Field2D") %} + const auto {{mixed_base_ind}} = fieldmesh->ind2Dto3D({{index_var}}); + {% if (operator == "/") and (rhs == "Field2D") %} const auto tmp = 1.0 / {{rhs.mixed_index}}; - for (int {{jz_var}} = 0; {{jz_var}} < fieldmesh->LocalNz; ++{{jz_var}}){ - (*this)[{{mixed_base_ind}} + {{jz_var}}] *= tmp; + for (int {{jz_var}} = 0; {{jz_var}} < fieldmesh->LocalNz; ++{{jz_var}}){ + (*this)[{{mixed_base_ind}} + {{jz_var}}] *= tmp; {% else %} for (int {{jz_var}} = 0; {{jz_var}} < fieldmesh->LocalNz; ++{{jz_var}}){ - (*this)[{{mixed_base_ind}} + {{jz_var}}] {{operator}}= {{rhs.index}}; + (*this)[{{mixed_base_ind}} + {{jz_var}}] {{operator}}= {{rhs.index}}; {% endif %} - } - } + } + } {% elif lhs == "FieldPerp" and (rhs == "Field3D" or rhs == "Field2D")%} Mesh *localmesh = this->getMesh(); @@ -109,7 +211,7 @@ int yind = this->getIndex(); const auto {{mixed_base_ind}} = localmesh->indPerpto3D({{index_var}}, yind); (*this)[{{index_var}}] {{operator}}= {{rhs.base_index}}; - } + } {% elif rhs == "FieldPerp" and (lhs == "Field3D" or lhs == "Field2D")%} Mesh *localmesh = this->getMesh(); @@ -117,7 +219,7 @@ int yind = {{rhs.name}}.getIndex(); const auto {{mixed_base_ind}} = localmesh->indPerpto3D({{index_var}}, yind); (*this)[{{base_ind_var}}] {{operator}}= {{rhs.index}}; - } + } {% elif (operator == "/") and (lhs == "Field3D" or lhs == "Field2D") and (rhs =="BoutReal") %} const auto tmp = 1.0 / {{rhs.index}}; {{region_loop}}({{index_var}}, this->getRegion({{region_name}})) { @@ -129,14 +231,14 @@ } {% endif %} - {% if lhs == "Field3D" %} + {% if lhs.region_type == "3D" %} track(rhs, "operator{{operator}}="); {% endif %} checkData(*this); } else { - {% if lhs == "Field3D" %} + {% if lhs.region_type == "3D" %} track(rhs, "operator{{operator}}="); {% endif %} (*this) = (*this) {{operator}} {{rhs.name}}; diff --git a/src/field/gen_fieldops.py b/src/field/gen_fieldops.py index 3e07d6fec4..91fd4082f0 100755 --- a/src/field/gen_fieldops.py +++ b/src/field/gen_fieldops.py @@ -16,7 +16,6 @@ from builtins import object import argparse -from collections import OrderedDict import contextlib from copy import deepcopy as copy import itertools @@ -53,17 +52,18 @@ def smart_open(filename, mode="r"): # The arthimetic operators -# OrderedDict to (try to) ensure consistency between python 2 & 3 -operators = OrderedDict( - [ - ("*", "multiplication"), - ("/", "division"), - ("+", "addition"), - ("-", "subtraction"), - ] -) +operators = { + "*": "multiplication", + "/": "division", + "+": "addition", + "-": "subtraction", +} + header = """// This file is autogenerated - see gen_fieldops.py +#include "fmt/format.h" +#include "bout/assert.hxx" +#include "bout/build_defines.hxx" #include #include #include @@ -103,7 +103,7 @@ def __init__( self.mixed_base_ind_var = mixed_base_ind_var # Note region_type isn't actually used currently but # may be useful in future. - if self.field_type == "Field3D": + if "Field3D" in self.field_type: self.region_type = "3D" elif self.field_type == "Field2D": self.region_type = "2D" @@ -157,6 +157,35 @@ def base_index(self): else: return "{self.name}[{self.mixed_base_ind_var}]".format(self=self) + @property + def yup(self): + """Returns {{name}}.yup(i) if it is a field with parallel slices. + If it is BoutReal just {{name}}""" + if self.field_type == "BoutReal": + return "{self.name}".format(self=self) + return "{self.name}.yup(i)".format(self=self) + + @property + def ydown(self): + """Returns {{name}}.ydown(i) if it is a field with parallel slices. + If it is BoutReal just {{name}}""" + if self.field_type == "BoutReal": + return "{self.name}".format(self=self) + return "{self.name}.ydown(i)".format(self=self) + + @property + def assertParallelSlices(self): + if self.field_type == "BoutReal": + return "" + + return f""" + ASSERT2({self.name}.hasParallelSlices()); + for (size_t i{{0}} ; i < {self.name}.numberParallelSlices() ; ++i) {{ + ASSERT2({self.ydown}.isAllocated()); + ASSERT2({self.yup}.isAllocated()); + }}\ + """ + def __eq__(self, other): try: return self.field_type == other.field_type @@ -183,6 +212,8 @@ def returnType(f1, f2): return copy(f1) elif f1 == "FieldPerp" or f2 == "FieldPerp": return copy(fieldPerp) + elif f1 == "Field3DParallel" or f2 == "Field3DParallel": + return copy(field3DPar) else: return copy(field3D) @@ -218,7 +249,6 @@ def returnType(f1, f2): region_loop = "BOUT_FOR" # Declare what fields we currently support: - # Field perp is currently missing field3D = Field( "Field3D", ["x", "y", "z"], @@ -226,6 +256,13 @@ def returnType(f1, f2): jz_var=jz_var, mixed_base_ind_var=mixed_base_ind_var, ) + field3DPar = Field( + "Field3DParallel", + ["x", "y", "z"], + index_var=index_var, + jz_var=jz_var, + mixed_base_ind_var=mixed_base_ind_var, + ) field2D = Field( "Field2D", ["x", "y"], @@ -248,7 +285,8 @@ def returnType(f1, f2): mixed_base_ind_var=mixed_base_ind_var, ) - fields = [field3D, field2D, fieldPerp, boutreal] + fields = (field3D, field2D, fieldPerp, boutreal) + fields2 = (field3D, field3DPar, boutreal) with smart_open(args.filename, "w") as f: f.write(header) @@ -258,10 +296,16 @@ def returnType(f1, f2): template = env.get_template("gen_fieldops.jinja") - for lhs, rhs in itertools.product(fields, fields): - # We don't have to define BoutReal BoutReal operations - if lhs == rhs == "BoutReal": + # We don't have to define BoutReal BoutReal operations + done = [(boutreal, boutreal)] + for lhs, rhs in itertools.chain( + itertools.product(fields, fields), + itertools.product(fields2, fields2), + ): + if (lhs, rhs) in done: continue + done.append((lhs, rhs)) + rhs = copy(rhs) lhs = copy(lhs) diff --git a/src/field/generated_fieldops.cxx b/src/field/generated_fieldops.cxx index 24bfde425c..f718d156fe 100644 --- a/src/field/generated_fieldops.cxx +++ b/src/field/generated_fieldops.cxx @@ -1,4 +1,7 @@ // This file is autogenerated - see gen_fieldops.py +#include "fmt/format.h" +#include "bout/assert.hxx" +#include "bout/build_defines.hxx" #include #include #include @@ -24,6 +27,32 @@ Field3D operator*(const Field3D& lhs, const Field3D& rhs) { return result; } +// Provide the C++ operator to update Field3D by multiplication with Field3D +Field3D& Field3D::update_multiplication_inplace(const Field3D& rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + ASSERT1_FIELDS_COMPATIBLE(*this, rhs); + + // Delete existing parallel slices. We don't update parallel slices, so any + // that currently exist will be incorrect. + clearParallelSlices(); + checkData(*this); + checkData(rhs); + + regionID = fieldmesh->getCommonRegion(regionID, rhs.getRegionID()); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] *= rhs[index]; } + + track(rhs, "operator*="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} *= {:s}", this->name, rhs.name); +#endif + + checkData(*this); + + return *this; +} + // Provide the C++ operator to update Field3D by multiplication with Field3D Field3D& Field3D::operator*=(const Field3D& rhs) { // only if data is unique we update the field @@ -31,14 +60,13 @@ Field3D& Field3D::operator*=(const Field3D& rhs) { if (data.unique()) { ASSERT1_FIELDS_COMPATIBLE(*this, rhs); - // Delete existing parallel slices. We don't copy parallel slices, so any + // Delete existing parallel slices. We don't update parallel slices, so any // that currently exist will be incorrect. clearParallelSlices(); - checkData(*this); checkData(rhs); - regionID = fieldmesh->getCommonRegion(regionID, rhs.regionID); + regionID = fieldmesh->getCommonRegion(regionID, rhs.getRegionID()); BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] *= rhs[index]; } @@ -71,6 +99,32 @@ Field3D operator/(const Field3D& lhs, const Field3D& rhs) { return result; } +// Provide the C++ operator to update Field3D by division with Field3D +Field3D& Field3D::update_division_inplace(const Field3D& rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + ASSERT1_FIELDS_COMPATIBLE(*this, rhs); + + // Delete existing parallel slices. We don't update parallel slices, so any + // that currently exist will be incorrect. + clearParallelSlices(); + checkData(*this); + checkData(rhs); + + regionID = fieldmesh->getCommonRegion(regionID, rhs.getRegionID()); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] /= rhs[index]; } + + track(rhs, "operator/="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} /= {:s}", this->name, rhs.name); +#endif + + checkData(*this); + + return *this; +} + // Provide the C++ operator to update Field3D by division with Field3D Field3D& Field3D::operator/=(const Field3D& rhs) { // only if data is unique we update the field @@ -78,14 +132,13 @@ Field3D& Field3D::operator/=(const Field3D& rhs) { if (data.unique()) { ASSERT1_FIELDS_COMPATIBLE(*this, rhs); - // Delete existing parallel slices. We don't copy parallel slices, so any + // Delete existing parallel slices. We don't update parallel slices, so any // that currently exist will be incorrect. clearParallelSlices(); - checkData(*this); checkData(rhs); - regionID = fieldmesh->getCommonRegion(regionID, rhs.regionID); + regionID = fieldmesh->getCommonRegion(regionID, rhs.getRegionID()); BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] /= rhs[index]; } @@ -118,6 +171,32 @@ Field3D operator+(const Field3D& lhs, const Field3D& rhs) { return result; } +// Provide the C++ operator to update Field3D by addition with Field3D +Field3D& Field3D::update_addition_inplace(const Field3D& rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + ASSERT1_FIELDS_COMPATIBLE(*this, rhs); + + // Delete existing parallel slices. We don't update parallel slices, so any + // that currently exist will be incorrect. + clearParallelSlices(); + checkData(*this); + checkData(rhs); + + regionID = fieldmesh->getCommonRegion(regionID, rhs.getRegionID()); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] += rhs[index]; } + + track(rhs, "operator+="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} += {:s}", this->name, rhs.name); +#endif + + checkData(*this); + + return *this; +} + // Provide the C++ operator to update Field3D by addition with Field3D Field3D& Field3D::operator+=(const Field3D& rhs) { // only if data is unique we update the field @@ -125,14 +204,13 @@ Field3D& Field3D::operator+=(const Field3D& rhs) { if (data.unique()) { ASSERT1_FIELDS_COMPATIBLE(*this, rhs); - // Delete existing parallel slices. We don't copy parallel slices, so any + // Delete existing parallel slices. We don't update parallel slices, so any // that currently exist will be incorrect. clearParallelSlices(); - checkData(*this); checkData(rhs); - regionID = fieldmesh->getCommonRegion(regionID, rhs.regionID); + regionID = fieldmesh->getCommonRegion(regionID, rhs.getRegionID()); BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] += rhs[index]; } @@ -165,6 +243,32 @@ Field3D operator-(const Field3D& lhs, const Field3D& rhs) { return result; } +// Provide the C++ operator to update Field3D by subtraction with Field3D +Field3D& Field3D::update_subtraction_inplace(const Field3D& rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + ASSERT1_FIELDS_COMPATIBLE(*this, rhs); + + // Delete existing parallel slices. We don't update parallel slices, so any + // that currently exist will be incorrect. + clearParallelSlices(); + checkData(*this); + checkData(rhs); + + regionID = fieldmesh->getCommonRegion(regionID, rhs.getRegionID()); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] -= rhs[index]; } + + track(rhs, "operator-="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} -= {:s}", this->name, rhs.name); +#endif + + checkData(*this); + + return *this; +} + // Provide the C++ operator to update Field3D by subtraction with Field3D Field3D& Field3D::operator-=(const Field3D& rhs) { // only if data is unique we update the field @@ -172,14 +276,13 @@ Field3D& Field3D::operator-=(const Field3D& rhs) { if (data.unique()) { ASSERT1_FIELDS_COMPATIBLE(*this, rhs); - // Delete existing parallel slices. We don't copy parallel slices, so any + // Delete existing parallel slices. We don't update parallel slices, so any // that currently exist will be incorrect. clearParallelSlices(); - checkData(*this); checkData(rhs); - regionID = fieldmesh->getCommonRegion(regionID, rhs.regionID); + regionID = fieldmesh->getCommonRegion(regionID, rhs.getRegionID()); BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] -= rhs[index]; } @@ -217,6 +320,35 @@ Field3D operator*(const Field3D& lhs, const Field2D& rhs) { return result; } +// Provide the C++ operator to update Field3D by multiplication with Field2D +Field3D& Field3D::update_multiplication_inplace(const Field2D& rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + ASSERT1_FIELDS_COMPATIBLE(*this, rhs); + + // Delete existing parallel slices. We don't update parallel slices, so any + // that currently exist will be incorrect. + clearParallelSlices(); + checkData(*this); + checkData(rhs); + + BOUT_FOR(index, rhs.getRegion("RGN_ALL")) { + const auto base_ind = fieldmesh->ind2Dto3D(index); + for (int jz = 0; jz < fieldmesh->LocalNz; ++jz) { + (*this)[base_ind + jz] *= rhs[index]; + } + } + + track(rhs, "operator*="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} *= {:s}", this->name, rhs.name); +#endif + + checkData(*this); + + return *this; +} + // Provide the C++ operator to update Field3D by multiplication with Field2D Field3D& Field3D::operator*=(const Field2D& rhs) { // only if data is unique we update the field @@ -224,10 +356,9 @@ Field3D& Field3D::operator*=(const Field2D& rhs) { if (data.unique()) { ASSERT1_FIELDS_COMPATIBLE(*this, rhs); - // Delete existing parallel slices. We don't copy parallel slices, so any + // Delete existing parallel slices. We don't update parallel slices, so any // that currently exist will be incorrect. clearParallelSlices(); - checkData(*this); checkData(rhs); @@ -273,6 +404,36 @@ Field3D operator/(const Field3D& lhs, const Field2D& rhs) { return result; } +// Provide the C++ operator to update Field3D by division with Field2D +Field3D& Field3D::update_division_inplace(const Field2D& rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + ASSERT1_FIELDS_COMPATIBLE(*this, rhs); + + // Delete existing parallel slices. We don't update parallel slices, so any + // that currently exist will be incorrect. + clearParallelSlices(); + checkData(*this); + checkData(rhs); + + BOUT_FOR(index, rhs.getRegion("RGN_ALL")) { + const auto base_ind = fieldmesh->ind2Dto3D(index); + const auto tmp = 1.0 / rhs[index]; + for (int jz = 0; jz < fieldmesh->LocalNz; ++jz) { + (*this)[base_ind + jz] *= tmp; + } + } + + track(rhs, "operator/="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} /= {:s}", this->name, rhs.name); +#endif + + checkData(*this); + + return *this; +} + // Provide the C++ operator to update Field3D by division with Field2D Field3D& Field3D::operator/=(const Field2D& rhs) { // only if data is unique we update the field @@ -280,10 +441,9 @@ Field3D& Field3D::operator/=(const Field2D& rhs) { if (data.unique()) { ASSERT1_FIELDS_COMPATIBLE(*this, rhs); - // Delete existing parallel slices. We don't copy parallel slices, so any + // Delete existing parallel slices. We don't update parallel slices, so any // that currently exist will be incorrect. clearParallelSlices(); - checkData(*this); checkData(rhs); @@ -329,6 +489,35 @@ Field3D operator+(const Field3D& lhs, const Field2D& rhs) { return result; } +// Provide the C++ operator to update Field3D by addition with Field2D +Field3D& Field3D::update_addition_inplace(const Field2D& rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + ASSERT1_FIELDS_COMPATIBLE(*this, rhs); + + // Delete existing parallel slices. We don't update parallel slices, so any + // that currently exist will be incorrect. + clearParallelSlices(); + checkData(*this); + checkData(rhs); + + BOUT_FOR(index, rhs.getRegion("RGN_ALL")) { + const auto base_ind = fieldmesh->ind2Dto3D(index); + for (int jz = 0; jz < fieldmesh->LocalNz; ++jz) { + (*this)[base_ind + jz] += rhs[index]; + } + } + + track(rhs, "operator+="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} += {:s}", this->name, rhs.name); +#endif + + checkData(*this); + + return *this; +} + // Provide the C++ operator to update Field3D by addition with Field2D Field3D& Field3D::operator+=(const Field2D& rhs) { // only if data is unique we update the field @@ -336,10 +525,9 @@ Field3D& Field3D::operator+=(const Field2D& rhs) { if (data.unique()) { ASSERT1_FIELDS_COMPATIBLE(*this, rhs); - // Delete existing parallel slices. We don't copy parallel slices, so any + // Delete existing parallel slices. We don't update parallel slices, so any // that currently exist will be incorrect. clearParallelSlices(); - checkData(*this); checkData(rhs); @@ -384,6 +572,35 @@ Field3D operator-(const Field3D& lhs, const Field2D& rhs) { return result; } +// Provide the C++ operator to update Field3D by subtraction with Field2D +Field3D& Field3D::update_subtraction_inplace(const Field2D& rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + ASSERT1_FIELDS_COMPATIBLE(*this, rhs); + + // Delete existing parallel slices. We don't update parallel slices, so any + // that currently exist will be incorrect. + clearParallelSlices(); + checkData(*this); + checkData(rhs); + + BOUT_FOR(index, rhs.getRegion("RGN_ALL")) { + const auto base_ind = fieldmesh->ind2Dto3D(index); + for (int jz = 0; jz < fieldmesh->LocalNz; ++jz) { + (*this)[base_ind + jz] -= rhs[index]; + } + } + + track(rhs, "operator-="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} -= {:s}", this->name, rhs.name); +#endif + + checkData(*this); + + return *this; +} + // Provide the C++ operator to update Field3D by subtraction with Field2D Field3D& Field3D::operator-=(const Field2D& rhs) { // only if data is unique we update the field @@ -391,10 +608,9 @@ Field3D& Field3D::operator-=(const Field2D& rhs) { if (data.unique()) { ASSERT1_FIELDS_COMPATIBLE(*this, rhs); - // Delete existing parallel slices. We don't copy parallel slices, so any + // Delete existing parallel slices. We don't update parallel slices, so any // that currently exist will be incorrect. clearParallelSlices(); - checkData(*this); checkData(rhs); @@ -513,16 +729,38 @@ Field3D operator*(const Field3D& lhs, const BoutReal rhs) { return result; } +// Provide the C++ operator to update Field3D by multiplication with BoutReal +Field3D& Field3D::update_multiplication_inplace(const BoutReal rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + + // Delete existing parallel slices. We don't update parallel slices, so any + // that currently exist will be incorrect. + clearParallelSlices(); + checkData(*this); + checkData(rhs); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] *= rhs; } + + track(rhs, "operator*="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} *= {:s}", this->name, "BR"); +#endif + + checkData(*this); + + return *this; +} + // Provide the C++ operator to update Field3D by multiplication with BoutReal Field3D& Field3D::operator*=(const BoutReal rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { - // Delete existing parallel slices. We don't copy parallel slices, so any + // Delete existing parallel slices. We don't update parallel slices, so any // that currently exist will be incorrect. clearParallelSlices(); - checkData(*this); checkData(rhs); @@ -557,16 +795,39 @@ Field3D operator/(const Field3D& lhs, const BoutReal rhs) { return result; } +// Provide the C++ operator to update Field3D by division with BoutReal +Field3D& Field3D::update_division_inplace(const BoutReal rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + + // Delete existing parallel slices. We don't update parallel slices, so any + // that currently exist will be incorrect. + clearParallelSlices(); + checkData(*this); + checkData(rhs); + + const auto tmp = 1.0 / rhs; + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] *= tmp; } + + track(rhs, "operator/="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} /= {:s}", this->name, "BR"); +#endif + + checkData(*this); + + return *this; +} + // Provide the C++ operator to update Field3D by division with BoutReal Field3D& Field3D::operator/=(const BoutReal rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { - // Delete existing parallel slices. We don't copy parallel slices, so any + // Delete existing parallel slices. We don't update parallel slices, so any // that currently exist will be incorrect. clearParallelSlices(); - checkData(*this); checkData(rhs); @@ -601,16 +862,38 @@ Field3D operator+(const Field3D& lhs, const BoutReal rhs) { return result; } +// Provide the C++ operator to update Field3D by addition with BoutReal +Field3D& Field3D::update_addition_inplace(const BoutReal rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + + // Delete existing parallel slices. We don't update parallel slices, so any + // that currently exist will be incorrect. + clearParallelSlices(); + checkData(*this); + checkData(rhs); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] += rhs; } + + track(rhs, "operator+="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} += {:s}", this->name, "BR"); +#endif + + checkData(*this); + + return *this; +} + // Provide the C++ operator to update Field3D by addition with BoutReal Field3D& Field3D::operator+=(const BoutReal rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { - // Delete existing parallel slices. We don't copy parallel slices, so any + // Delete existing parallel slices. We don't update parallel slices, so any // that currently exist will be incorrect. clearParallelSlices(); - checkData(*this); checkData(rhs); @@ -644,16 +927,38 @@ Field3D operator-(const Field3D& lhs, const BoutReal rhs) { return result; } +// Provide the C++ operator to update Field3D by subtraction with BoutReal +Field3D& Field3D::update_subtraction_inplace(const BoutReal rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + + // Delete existing parallel slices. We don't update parallel slices, so any + // that currently exist will be incorrect. + clearParallelSlices(); + checkData(*this); + checkData(rhs); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] -= rhs; } + + track(rhs, "operator-="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} -= {:s}", this->name, "BR"); +#endif + + checkData(*this); + + return *this; +} + // Provide the C++ operator to update Field3D by subtraction with BoutReal Field3D& Field3D::operator-=(const BoutReal rhs) { // only if data is unique we update the field // otherwise just call the non-inplace version if (data.unique()) { - // Delete existing parallel slices. We don't copy parallel slices, so any + // Delete existing parallel slices. We don't update parallel slices, so any // that currently exist will be incorrect. clearParallelSlices(); - checkData(*this); checkData(rhs); @@ -1960,3 +2265,1280 @@ FieldPerp operator-(const BoutReal lhs, const FieldPerp& rhs) { checkData(result); return result; } + +// Provide the C++ wrapper for multiplication of Field3D and Field3DParallel +Field3DParallel operator*(const Field3D& lhs, const Field3DParallel& rhs) { + ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); + + Field3DParallel result{emptyFrom(rhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); + if (result.isFci()) { + + ASSERT2(lhs.hasParallelSlices()); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + ASSERT2(lhs.ydown(i).isAllocated()); + ASSERT2(lhs.yup(i).isAllocated()); + } + + ASSERT2(rhs.hasParallelSlices()); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + ASSERT2(rhs.ydown(i).isAllocated()); + ASSERT2(rhs.yup(i).isAllocated()); + } + result.splitParallelSlices(); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs.yup(i) * rhs.yup(i); + result.ydown(i) = lhs.ydown(i) * rhs.ydown(i); + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs[index] * rhs[index]; + } + + checkData(result); + return result; +} + +// Provide the C++ wrapper for division of Field3D and Field3DParallel +Field3DParallel operator/(const Field3D& lhs, const Field3DParallel& rhs) { + ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); + + Field3DParallel result{emptyFrom(rhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); + if (result.isFci()) { + + ASSERT2(lhs.hasParallelSlices()); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + ASSERT2(lhs.ydown(i).isAllocated()); + ASSERT2(lhs.yup(i).isAllocated()); + } + + ASSERT2(rhs.hasParallelSlices()); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + ASSERT2(rhs.ydown(i).isAllocated()); + ASSERT2(rhs.yup(i).isAllocated()); + } + result.splitParallelSlices(); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs.yup(i) / rhs.yup(i); + result.ydown(i) = lhs.ydown(i) / rhs.ydown(i); + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs[index] / rhs[index]; + } + + checkData(result); + return result; +} + +// Provide the C++ wrapper for addition of Field3D and Field3DParallel +Field3DParallel operator+(const Field3D& lhs, const Field3DParallel& rhs) { + ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); + + Field3DParallel result{emptyFrom(rhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); + if (result.isFci()) { + + ASSERT2(lhs.hasParallelSlices()); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + ASSERT2(lhs.ydown(i).isAllocated()); + ASSERT2(lhs.yup(i).isAllocated()); + } + + ASSERT2(rhs.hasParallelSlices()); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + ASSERT2(rhs.ydown(i).isAllocated()); + ASSERT2(rhs.yup(i).isAllocated()); + } + result.splitParallelSlices(); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs.yup(i) + rhs.yup(i); + result.ydown(i) = lhs.ydown(i) + rhs.ydown(i); + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs[index] + rhs[index]; + } + + checkData(result); + return result; +} + +// Provide the C++ wrapper for subtraction of Field3D and Field3DParallel +Field3DParallel operator-(const Field3D& lhs, const Field3DParallel& rhs) { + ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); + + Field3DParallel result{emptyFrom(rhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); + if (result.isFci()) { + + ASSERT2(lhs.hasParallelSlices()); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + ASSERT2(lhs.ydown(i).isAllocated()); + ASSERT2(lhs.yup(i).isAllocated()); + } + + ASSERT2(rhs.hasParallelSlices()); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + ASSERT2(rhs.ydown(i).isAllocated()); + ASSERT2(rhs.yup(i).isAllocated()); + } + result.splitParallelSlices(); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs.yup(i) - rhs.yup(i); + result.ydown(i) = lhs.ydown(i) - rhs.ydown(i); + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs[index] - rhs[index]; + } + + checkData(result); + return result; +} + +// Provide the C++ wrapper for multiplication of Field3DParallel and Field3D +Field3DParallel operator*(const Field3DParallel& lhs, const Field3D& rhs) { + ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); + + Field3DParallel result{emptyFrom(lhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); + if (result.isFci()) { + + ASSERT2(lhs.hasParallelSlices()); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + ASSERT2(lhs.ydown(i).isAllocated()); + ASSERT2(lhs.yup(i).isAllocated()); + } + + ASSERT2(rhs.hasParallelSlices()); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + ASSERT2(rhs.ydown(i).isAllocated()); + ASSERT2(rhs.yup(i).isAllocated()); + } + result.splitParallelSlices(); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs.yup(i) * rhs.yup(i); + result.ydown(i) = lhs.ydown(i) * rhs.ydown(i); + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs[index] * rhs[index]; + } + + checkData(result); + return result; +} + +// Provide the C++ operator to update Field3DParallel by multiplication with Field3D +Field3DParallel& Field3DParallel::operator*=(const Field3D& rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + if (data.unique() or isRef) { + ASSERT1_FIELDS_COMPATIBLE(*this, rhs); + + if (this->isFci()) { + if (isRef) { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i).update_multiplication_inplace(rhs.yup(i)); + ydown(i).update_multiplication_inplace(rhs.ydown(i)); + } + } else { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i) *= rhs.yup(i); + ydown(i) *= rhs.ydown(i); + } + } + } else { + clearParallelSlices(); + } + checkData(*this); + checkData(rhs); + + regionID = fieldmesh->getCommonRegion(regionID, rhs.getRegionID()); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] *= rhs[index]; } + + track(rhs, "operator*="); + + checkData(*this); + + } else { + track(rhs, "operator*="); + (*this) = (*this) * rhs; + } + return *this; +} + +// Provide the C++ wrapper for division of Field3DParallel and Field3D +Field3DParallel operator/(const Field3DParallel& lhs, const Field3D& rhs) { + ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); + + Field3DParallel result{emptyFrom(lhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); + if (result.isFci()) { + + ASSERT2(lhs.hasParallelSlices()); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + ASSERT2(lhs.ydown(i).isAllocated()); + ASSERT2(lhs.yup(i).isAllocated()); + } + + ASSERT2(rhs.hasParallelSlices()); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + ASSERT2(rhs.ydown(i).isAllocated()); + ASSERT2(rhs.yup(i).isAllocated()); + } + result.splitParallelSlices(); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs.yup(i) / rhs.yup(i); + result.ydown(i) = lhs.ydown(i) / rhs.ydown(i); + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs[index] / rhs[index]; + } + + checkData(result); + return result; +} + +// Provide the C++ operator to update Field3DParallel by division with Field3D +Field3DParallel& Field3DParallel::operator/=(const Field3D& rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + if (data.unique() or isRef) { + ASSERT1_FIELDS_COMPATIBLE(*this, rhs); + + if (this->isFci()) { + if (isRef) { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i).update_division_inplace(rhs.yup(i)); + ydown(i).update_division_inplace(rhs.ydown(i)); + } + } else { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i) /= rhs.yup(i); + ydown(i) /= rhs.ydown(i); + } + } + } else { + clearParallelSlices(); + } + checkData(*this); + checkData(rhs); + + regionID = fieldmesh->getCommonRegion(regionID, rhs.getRegionID()); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] /= rhs[index]; } + + track(rhs, "operator/="); + + checkData(*this); + + } else { + track(rhs, "operator/="); + (*this) = (*this) / rhs; + } + return *this; +} + +// Provide the C++ wrapper for addition of Field3DParallel and Field3D +Field3DParallel operator+(const Field3DParallel& lhs, const Field3D& rhs) { + ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); + + Field3DParallel result{emptyFrom(lhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); + if (result.isFci()) { + + ASSERT2(lhs.hasParallelSlices()); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + ASSERT2(lhs.ydown(i).isAllocated()); + ASSERT2(lhs.yup(i).isAllocated()); + } + + ASSERT2(rhs.hasParallelSlices()); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + ASSERT2(rhs.ydown(i).isAllocated()); + ASSERT2(rhs.yup(i).isAllocated()); + } + result.splitParallelSlices(); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs.yup(i) + rhs.yup(i); + result.ydown(i) = lhs.ydown(i) + rhs.ydown(i); + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs[index] + rhs[index]; + } + + checkData(result); + return result; +} + +// Provide the C++ operator to update Field3DParallel by addition with Field3D +Field3DParallel& Field3DParallel::operator+=(const Field3D& rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + if (data.unique() or isRef) { + ASSERT1_FIELDS_COMPATIBLE(*this, rhs); + + if (this->isFci()) { + if (isRef) { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i).update_addition_inplace(rhs.yup(i)); + ydown(i).update_addition_inplace(rhs.ydown(i)); + } + } else { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i) += rhs.yup(i); + ydown(i) += rhs.ydown(i); + } + } + } else { + clearParallelSlices(); + } + checkData(*this); + checkData(rhs); + + regionID = fieldmesh->getCommonRegion(regionID, rhs.getRegionID()); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] += rhs[index]; } + + track(rhs, "operator+="); + + checkData(*this); + + } else { + track(rhs, "operator+="); + (*this) = (*this) + rhs; + } + return *this; +} + +// Provide the C++ wrapper for subtraction of Field3DParallel and Field3D +Field3DParallel operator-(const Field3DParallel& lhs, const Field3D& rhs) { + ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); + + Field3DParallel result{emptyFrom(lhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); + if (result.isFci()) { + + ASSERT2(lhs.hasParallelSlices()); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + ASSERT2(lhs.ydown(i).isAllocated()); + ASSERT2(lhs.yup(i).isAllocated()); + } + + ASSERT2(rhs.hasParallelSlices()); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + ASSERT2(rhs.ydown(i).isAllocated()); + ASSERT2(rhs.yup(i).isAllocated()); + } + result.splitParallelSlices(); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs.yup(i) - rhs.yup(i); + result.ydown(i) = lhs.ydown(i) - rhs.ydown(i); + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs[index] - rhs[index]; + } + + checkData(result); + return result; +} + +// Provide the C++ operator to update Field3DParallel by subtraction with Field3D +Field3DParallel& Field3DParallel::operator-=(const Field3D& rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + if (data.unique() or isRef) { + ASSERT1_FIELDS_COMPATIBLE(*this, rhs); + + if (this->isFci()) { + if (isRef) { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i).update_subtraction_inplace(rhs.yup(i)); + ydown(i).update_subtraction_inplace(rhs.ydown(i)); + } + } else { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i) -= rhs.yup(i); + ydown(i) -= rhs.ydown(i); + } + } + } else { + clearParallelSlices(); + } + checkData(*this); + checkData(rhs); + + regionID = fieldmesh->getCommonRegion(regionID, rhs.getRegionID()); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] -= rhs[index]; } + + track(rhs, "operator-="); + + checkData(*this); + + } else { + track(rhs, "operator-="); + (*this) = (*this) - rhs; + } + return *this; +} + +// Provide the C++ wrapper for multiplication of Field3DParallel and Field3DParallel +Field3DParallel operator*(const Field3DParallel& lhs, const Field3DParallel& rhs) { + ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); + + Field3DParallel result{emptyFrom(lhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); + if (result.isFci()) { + + ASSERT2(lhs.hasParallelSlices()); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + ASSERT2(lhs.ydown(i).isAllocated()); + ASSERT2(lhs.yup(i).isAllocated()); + } + + ASSERT2(rhs.hasParallelSlices()); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + ASSERT2(rhs.ydown(i).isAllocated()); + ASSERT2(rhs.yup(i).isAllocated()); + } + result.splitParallelSlices(); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs.yup(i) * rhs.yup(i); + result.ydown(i) = lhs.ydown(i) * rhs.ydown(i); + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs[index] * rhs[index]; + } + + checkData(result); + return result; +} + +// Provide the C++ operator to update Field3DParallel by multiplication with Field3DParallel +Field3DParallel& Field3DParallel::operator*=(const Field3DParallel& rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + if (data.unique() or isRef) { + ASSERT1_FIELDS_COMPATIBLE(*this, rhs); + + if (this->isFci()) { + if (isRef) { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i).update_multiplication_inplace(rhs); + ydown(i).update_multiplication_inplace(rhs); + } + } else { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i) *= rhs; + ydown(i) *= rhs; + } + } + } else { + clearParallelSlices(); + } + checkData(*this); + checkData(rhs); + + regionID = fieldmesh->getCommonRegion(regionID, rhs.getRegionID()); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] *= rhs[index]; } + + track(rhs, "operator*="); + + checkData(*this); + + } else { + track(rhs, "operator*="); + (*this) = (*this) * rhs; + } + return *this; +} + +// Provide the C++ wrapper for division of Field3DParallel and Field3DParallel +Field3DParallel operator/(const Field3DParallel& lhs, const Field3DParallel& rhs) { + ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); + + Field3DParallel result{emptyFrom(lhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); + if (result.isFci()) { + + ASSERT2(lhs.hasParallelSlices()); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + ASSERT2(lhs.ydown(i).isAllocated()); + ASSERT2(lhs.yup(i).isAllocated()); + } + + ASSERT2(rhs.hasParallelSlices()); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + ASSERT2(rhs.ydown(i).isAllocated()); + ASSERT2(rhs.yup(i).isAllocated()); + } + result.splitParallelSlices(); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs.yup(i) / rhs.yup(i); + result.ydown(i) = lhs.ydown(i) / rhs.ydown(i); + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs[index] / rhs[index]; + } + + checkData(result); + return result; +} + +// Provide the C++ operator to update Field3DParallel by division with Field3DParallel +Field3DParallel& Field3DParallel::operator/=(const Field3DParallel& rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + if (data.unique() or isRef) { + ASSERT1_FIELDS_COMPATIBLE(*this, rhs); + + if (this->isFci()) { + if (isRef) { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i).update_division_inplace(rhs); + ydown(i).update_division_inplace(rhs); + } + } else { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i) /= rhs; + ydown(i) /= rhs; + } + } + } else { + clearParallelSlices(); + } + checkData(*this); + checkData(rhs); + + regionID = fieldmesh->getCommonRegion(regionID, rhs.getRegionID()); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] /= rhs[index]; } + + track(rhs, "operator/="); + + checkData(*this); + + } else { + track(rhs, "operator/="); + (*this) = (*this) / rhs; + } + return *this; +} + +// Provide the C++ wrapper for addition of Field3DParallel and Field3DParallel +Field3DParallel operator+(const Field3DParallel& lhs, const Field3DParallel& rhs) { + ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); + + Field3DParallel result{emptyFrom(lhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); + if (result.isFci()) { + + ASSERT2(lhs.hasParallelSlices()); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + ASSERT2(lhs.ydown(i).isAllocated()); + ASSERT2(lhs.yup(i).isAllocated()); + } + + ASSERT2(rhs.hasParallelSlices()); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + ASSERT2(rhs.ydown(i).isAllocated()); + ASSERT2(rhs.yup(i).isAllocated()); + } + result.splitParallelSlices(); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs.yup(i) + rhs.yup(i); + result.ydown(i) = lhs.ydown(i) + rhs.ydown(i); + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs[index] + rhs[index]; + } + + checkData(result); + return result; +} + +// Provide the C++ operator to update Field3DParallel by addition with Field3DParallel +Field3DParallel& Field3DParallel::operator+=(const Field3DParallel& rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + if (data.unique() or isRef) { + ASSERT1_FIELDS_COMPATIBLE(*this, rhs); + + if (this->isFci()) { + if (isRef) { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i).update_addition_inplace(rhs); + ydown(i).update_addition_inplace(rhs); + } + } else { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i) += rhs; + ydown(i) += rhs; + } + } + } else { + clearParallelSlices(); + } + checkData(*this); + checkData(rhs); + + regionID = fieldmesh->getCommonRegion(regionID, rhs.getRegionID()); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] += rhs[index]; } + + track(rhs, "operator+="); + + checkData(*this); + + } else { + track(rhs, "operator+="); + (*this) = (*this) + rhs; + } + return *this; +} + +// Provide the C++ wrapper for subtraction of Field3DParallel and Field3DParallel +Field3DParallel operator-(const Field3DParallel& lhs, const Field3DParallel& rhs) { + ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); + + Field3DParallel result{emptyFrom(lhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); + if (result.isFci()) { + + ASSERT2(lhs.hasParallelSlices()); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + ASSERT2(lhs.ydown(i).isAllocated()); + ASSERT2(lhs.yup(i).isAllocated()); + } + + ASSERT2(rhs.hasParallelSlices()); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + ASSERT2(rhs.ydown(i).isAllocated()); + ASSERT2(rhs.yup(i).isAllocated()); + } + result.splitParallelSlices(); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs.yup(i) - rhs.yup(i); + result.ydown(i) = lhs.ydown(i) - rhs.ydown(i); + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs[index] - rhs[index]; + } + + checkData(result); + return result; +} + +// Provide the C++ operator to update Field3DParallel by subtraction with Field3DParallel +Field3DParallel& Field3DParallel::operator-=(const Field3DParallel& rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + if (data.unique() or isRef) { + ASSERT1_FIELDS_COMPATIBLE(*this, rhs); + + if (this->isFci()) { + if (isRef) { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i).update_subtraction_inplace(rhs); + ydown(i).update_subtraction_inplace(rhs); + } + } else { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i) -= rhs; + ydown(i) -= rhs; + } + } + } else { + clearParallelSlices(); + } + checkData(*this); + checkData(rhs); + + regionID = fieldmesh->getCommonRegion(regionID, rhs.getRegionID()); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] -= rhs[index]; } + + track(rhs, "operator-="); + + checkData(*this); + + } else { + track(rhs, "operator-="); + (*this) = (*this) - rhs; + } + return *this; +} + +// Provide the C++ wrapper for multiplication of Field3DParallel and BoutReal +Field3DParallel operator*(const Field3DParallel& lhs, const BoutReal rhs) { + + Field3DParallel result{emptyFrom(lhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(lhs.getRegionID()); + if (result.isFci()) { + + ASSERT2(lhs.hasParallelSlices()); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + ASSERT2(lhs.ydown(i).isAllocated()); + ASSERT2(lhs.yup(i).isAllocated()); + } + + result.splitParallelSlices(); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs.yup(i) * rhs; + result.ydown(i) = lhs.ydown(i) * rhs; + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs[index] * rhs; + } + + checkData(result); + return result; +} + +// Provide the C++ operator to update Field3DParallel by multiplication with BoutReal +Field3DParallel& Field3DParallel::operator*=(const BoutReal rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + if (data.unique() or isRef) { + + if (this->isFci()) { + if (isRef) { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i).update_multiplication_inplace(rhs); + ydown(i).update_multiplication_inplace(rhs); + } + } else { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i) *= rhs; + ydown(i) *= rhs; + } + } + } else { + clearParallelSlices(); + } + checkData(*this); + checkData(rhs); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] *= rhs; } + + track(rhs, "operator*="); + + checkData(*this); + + } else { + track(rhs, "operator*="); + (*this) = (*this) * rhs; + } + return *this; +} + +// Provide the C++ wrapper for division of Field3DParallel and BoutReal +Field3DParallel operator/(const Field3DParallel& lhs, const BoutReal rhs) { + + Field3DParallel result{emptyFrom(lhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(lhs.getRegionID()); + if (result.isFci()) { + + ASSERT2(lhs.hasParallelSlices()); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + ASSERT2(lhs.ydown(i).isAllocated()); + ASSERT2(lhs.yup(i).isAllocated()); + } + + result.splitParallelSlices(); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs.yup(i) / rhs; + result.ydown(i) = lhs.ydown(i) / rhs; + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + const auto tmp = 1.0 / rhs; + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs[index] * tmp; + } + + checkData(result); + return result; +} + +// Provide the C++ operator to update Field3DParallel by division with BoutReal +Field3DParallel& Field3DParallel::operator/=(const BoutReal rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + if (data.unique() or isRef) { + + if (this->isFci()) { + if (isRef) { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i).update_division_inplace(rhs); + ydown(i).update_division_inplace(rhs); + } + } else { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i) /= rhs; + ydown(i) /= rhs; + } + } + } else { + clearParallelSlices(); + } + checkData(*this); + checkData(rhs); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] /= rhs; } + + track(rhs, "operator/="); + + checkData(*this); + + } else { + track(rhs, "operator/="); + (*this) = (*this) / rhs; + } + return *this; +} + +// Provide the C++ wrapper for addition of Field3DParallel and BoutReal +Field3DParallel operator+(const Field3DParallel& lhs, const BoutReal rhs) { + + Field3DParallel result{emptyFrom(lhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(lhs.getRegionID()); + if (result.isFci()) { + + ASSERT2(lhs.hasParallelSlices()); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + ASSERT2(lhs.ydown(i).isAllocated()); + ASSERT2(lhs.yup(i).isAllocated()); + } + + result.splitParallelSlices(); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs.yup(i) + rhs; + result.ydown(i) = lhs.ydown(i) + rhs; + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs[index] + rhs; + } + + checkData(result); + return result; +} + +// Provide the C++ operator to update Field3DParallel by addition with BoutReal +Field3DParallel& Field3DParallel::operator+=(const BoutReal rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + if (data.unique() or isRef) { + + if (this->isFci()) { + if (isRef) { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i).update_addition_inplace(rhs); + ydown(i).update_addition_inplace(rhs); + } + } else { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i) += rhs; + ydown(i) += rhs; + } + } + } else { + clearParallelSlices(); + } + checkData(*this); + checkData(rhs); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] += rhs; } + + track(rhs, "operator+="); + + checkData(*this); + + } else { + track(rhs, "operator+="); + (*this) = (*this) + rhs; + } + return *this; +} + +// Provide the C++ wrapper for subtraction of Field3DParallel and BoutReal +Field3DParallel operator-(const Field3DParallel& lhs, const BoutReal rhs) { + + Field3DParallel result{emptyFrom(lhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(lhs.getRegionID()); + if (result.isFci()) { + + ASSERT2(lhs.hasParallelSlices()); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + ASSERT2(lhs.ydown(i).isAllocated()); + ASSERT2(lhs.yup(i).isAllocated()); + } + + result.splitParallelSlices(); + for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs.yup(i) - rhs; + result.ydown(i) = lhs.ydown(i) - rhs; + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs[index] - rhs; + } + + checkData(result); + return result; +} + +// Provide the C++ operator to update Field3DParallel by subtraction with BoutReal +Field3DParallel& Field3DParallel::operator-=(const BoutReal rhs) { + // only if data is unique we update the field + // otherwise just call the non-inplace version + if (data.unique() or isRef) { + + if (this->isFci()) { + if (isRef) { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i).update_subtraction_inplace(rhs); + ydown(i).update_subtraction_inplace(rhs); + } + } else { + for (size_t i{0}; i < yup_fields.size(); ++i) { + yup(i) -= rhs; + ydown(i) -= rhs; + } + } + } else { + clearParallelSlices(); + } + checkData(*this); + checkData(rhs); + + BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] -= rhs; } + + track(rhs, "operator-="); + + checkData(*this); + + } else { + track(rhs, "operator-="); + (*this) = (*this) - rhs; + } + return *this; +} + +// Provide the C++ wrapper for multiplication of BoutReal and Field3DParallel +Field3DParallel operator*(const BoutReal lhs, const Field3DParallel& rhs) { + + Field3DParallel result{emptyFrom(rhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(rhs.getRegionID()); + if (result.isFci()) { + + ASSERT2(rhs.hasParallelSlices()); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + ASSERT2(rhs.ydown(i).isAllocated()); + ASSERT2(rhs.yup(i).isAllocated()); + } + result.splitParallelSlices(); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs * rhs.yup(i); + result.ydown(i) = lhs * rhs.ydown(i); + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs * rhs[index]; + } + + checkData(result); + return result; +} + +// Provide the C++ wrapper for division of BoutReal and Field3DParallel +Field3DParallel operator/(const BoutReal lhs, const Field3DParallel& rhs) { + + Field3DParallel result{emptyFrom(rhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(rhs.getRegionID()); + if (result.isFci()) { + + ASSERT2(rhs.hasParallelSlices()); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + ASSERT2(rhs.ydown(i).isAllocated()); + ASSERT2(rhs.yup(i).isAllocated()); + } + result.splitParallelSlices(); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs / rhs.yup(i); + result.ydown(i) = lhs / rhs.ydown(i); + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs / rhs[index]; + } + + checkData(result); + return result; +} + +// Provide the C++ wrapper for addition of BoutReal and Field3DParallel +Field3DParallel operator+(const BoutReal lhs, const Field3DParallel& rhs) { + + Field3DParallel result{emptyFrom(rhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(rhs.getRegionID()); + if (result.isFci()) { + + ASSERT2(rhs.hasParallelSlices()); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + ASSERT2(rhs.ydown(i).isAllocated()); + ASSERT2(rhs.yup(i).isAllocated()); + } + result.splitParallelSlices(); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs + rhs.yup(i); + result.ydown(i) = lhs + rhs.ydown(i); + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs + rhs[index]; + } + + checkData(result); + return result; +} + +// Provide the C++ wrapper for subtraction of BoutReal and Field3DParallel +Field3DParallel operator-(const BoutReal lhs, const Field3DParallel& rhs) { + + Field3DParallel result{emptyFrom(rhs)}; + checkData(lhs); + checkData(rhs); + + result.setRegion(rhs.getRegionID()); + if (result.isFci()) { + + ASSERT2(rhs.hasParallelSlices()); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + ASSERT2(rhs.ydown(i).isAllocated()); + ASSERT2(rhs.yup(i).isAllocated()); + } + result.splitParallelSlices(); + for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { + result.yup(i) = lhs - rhs.yup(i); + result.ydown(i) = lhs - rhs.ydown(i); + } + + ASSERT2(result.hasParallelSlices()); + for (size_t i{0}; i < result.numberParallelSlices(); ++i) { + ASSERT2(result.ydown(i).isAllocated()); + ASSERT2(result.yup(i).isAllocated()); + } + } + + BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { + result[index] = lhs - rhs[index]; + } + + checkData(result); + return result; +} diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index d2634d3f88..1157f3a823 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -591,6 +591,9 @@ Coordinates::Coordinates(Mesh* mesh, Options* options) // IntShiftTorsion will not be used, but set to zero to avoid uninitialized field IntShiftTorsion = 0.; } + + // Allow transform to fix things up + transform->loadParallelMetrics(this); } Coordinates::Coordinates(Mesh* mesh, Options* options, const CELL_LOC loc, @@ -879,6 +882,8 @@ Coordinates::Coordinates(Mesh* mesh, Options* options, const CELL_LOC loc, true, true, false, transform.get()); } } + // Allow transform to fix things up + transform->loadParallelMetrics(this); } void Coordinates::outputVars(Options& output_options) { @@ -957,115 +962,119 @@ int Coordinates::geometry(bool recalculate_staggered, // Note: This calculation is completely general: metric // tensor can be 2D or 3D. For 2D, all DDZ terms are zero - G1_11 = 0.5 * g11 * DDX(g_11) + g12 * (DDX(g_12) - 0.5 * DDY(g_11)) - + g13 * (DDX(g_13) - 0.5 * DDZ(g_11)); - G1_22 = g11 * (DDY(g_12) - 0.5 * DDX(g_22)) + 0.5 * g12 * DDY(g_22) - + g13 * (DDY(g_23) - 0.5 * DDZ(g_22)); - G1_33 = g11 * (DDZ(g_13) - 0.5 * DDX(g_33)) + g12 * (DDZ(g_23) - 0.5 * DDY(g_33)) - + 0.5 * g13 * DDZ(g_33); - G1_12 = 0.5 * g11 * DDY(g_11) + 0.5 * g12 * DDX(g_22) - + 0.5 * g13 * (DDY(g_13) + DDX(g_23) - DDZ(g_12)); - G1_13 = 0.5 * g11 * DDZ(g_11) + 0.5 * g12 * (DDZ(g_12) + DDX(g_23) - DDY(g_13)) - + 0.5 * g13 * DDX(g_33); - G1_23 = 0.5 * g11 * (DDZ(g_12) + DDY(g_13) - DDX(g_23)) - + 0.5 * g12 * (DDZ(g_22) + DDY(g_23) - DDY(g_23)) - // + 0.5 *g13*(DDZ(g_32) + DDY(g_33) - DDZ(g_23)); - // which equals - + 0.5 * g13 * DDY(g_33); - - G2_11 = 0.5 * g12 * DDX(g_11) + g22 * (DDX(g_12) - 0.5 * DDY(g_11)) - + g23 * (DDX(g_13) - 0.5 * DDZ(g_11)); - G2_22 = g12 * (DDY(g_12) - 0.5 * DDX(g_22)) + 0.5 * g22 * DDY(g_22) - + g23 * (DDY(g23) - 0.5 * DDZ(g_22)); - G2_33 = g12 * (DDZ(g_13) - 0.5 * DDX(g_33)) + g22 * (DDZ(g_23) - 0.5 * DDY(g_33)) - + 0.5 * g23 * DDZ(g_33); - G2_12 = 0.5 * g12 * DDY(g_11) + 0.5 * g22 * DDX(g_22) - + 0.5 * g23 * (DDY(g_13) + DDX(g_23) - DDZ(g_12)); - G2_13 = - // 0.5 *g21*(DDZ(g_11) + DDX(g_13) - DDX(g_13)) - // which equals - 0.5 * g12 * (DDZ(g_11) + DDX(g_13) - DDX(g_13)) - // + 0.5 *g22*(DDZ(g_21) + DDX(g_23) - DDY(g_13)) - // which equals - + 0.5 * g22 * (DDZ(g_12) + DDX(g_23) - DDY(g_13)) - // + 0.5 *g23*(DDZ(g_31) + DDX(g_33) - DDZ(g_13)); - // which equals - + 0.5 * g23 * DDX(g_33); - G2_23 = 0.5 * g12 * (DDZ(g_12) + DDY(g_13) - DDX(g_23)) + 0.5 * g22 * DDZ(g_22) - + 0.5 * g23 * DDY(g_33); - - G3_11 = 0.5 * g13 * DDX(g_11) + g23 * (DDX(g_12) - 0.5 * DDY(g_11)) - + g33 * (DDX(g_13) - 0.5 * DDZ(g_11)); - G3_22 = g13 * (DDY(g_12) - 0.5 * DDX(g_22)) + 0.5 * g23 * DDY(g_22) - + g33 * (DDY(g_23) - 0.5 * DDZ(g_22)); - G3_33 = g13 * (DDZ(g_13) - 0.5 * DDX(g_33)) + g23 * (DDZ(g_23) - 0.5 * DDY(g_33)) - + 0.5 * g33 * DDZ(g_33); - G3_12 = - // 0.5 *g31*(DDY(g_11) + DDX(g_12) - DDX(g_12)) - // which equals to - 0.5 * g13 * DDY(g_11) - // + 0.5 *g32*(DDY(g_21) + DDX(g_22) - DDY(g_12)) - // which equals to - + 0.5 * g23 * DDX(g_22) - //+ 0.5 *g33*(DDY(g_31) + DDX(g_32) - DDZ(g_12)); - // which equals to - + 0.5 * g33 * (DDY(g_13) + DDX(g_23) - DDZ(g_12)); - G3_13 = 0.5 * g13 * DDZ(g_11) + 0.5 * g23 * (DDZ(g_12) + DDX(g_23) - DDY(g_13)) - + 0.5 * g33 * DDX(g_33); - G3_23 = 0.5 * g13 * (DDZ(g_12) + DDY(g_13) - DDX(g_23)) + 0.5 * g23 * DDZ(g_22) - + 0.5 * g33 * DDY(g_33); - - auto tmp = J * g12; - localmesh->communicate_no_slices(tmp); - G1 = (DDX(J * g11) + DDY(tmp) + DDZ(J * g13)) / J; - tmp = J * g22; - localmesh->communicate_no_slices(tmp); - G2 = (DDX(J * g12) + DDY(tmp) + DDZ(J * g23)) / J; - tmp = J * g23; - localmesh->communicate_no_slices(tmp); - G3 = (DDX(J * g13) + DDY(tmp) + DDZ(J * g33)) / J; - - // Communicate christoffel symbol terms - output_progress.write("\tCommunicating connection terms\n"); - - localmesh->communicate_no_slices(G1_11, G1_22, G1_33, G1_12, G1_13, G1_23, G2_11, G2_22, - G2_33, G2_12, G2_13, G2_23, G3_11, G3_22, G3_33, G3_12, - G3_13, G3_23, G1, G2, G3); - - // Set boundary guard cells of Christoffel symbol terms - // Ideally, when location is staggered, we would set the upper/outer boundary point - // correctly rather than by extrapolating here: e.g. if location==CELL_YLOW and we are - // at the upper y-boundary the x- and z-derivatives at yend+1 at the boundary can be - // calculated because the guard cells are available, while the y-derivative could be - // calculated from the CELL_CENTRE metric components (which have guard cells available - // past the boundary location). This would avoid the problem that the y-boundary on the - // CELL_YLOW grid is at a 'guard cell' location (yend+1). - // However, the above would require lots of special handling, so just extrapolate for - // now. - G1_11 = interpolateAndExtrapolate(G1_11, location, true, true, true, transform.get()); - G1_22 = interpolateAndExtrapolate(G1_22, location, true, true, true, transform.get()); - G1_33 = interpolateAndExtrapolate(G1_33, location, true, true, true, transform.get()); - G1_12 = interpolateAndExtrapolate(G1_12, location, true, true, true, transform.get()); - G1_13 = interpolateAndExtrapolate(G1_13, location, true, true, true, transform.get()); - G1_23 = interpolateAndExtrapolate(G1_23, location, true, true, true, transform.get()); - - G2_11 = interpolateAndExtrapolate(G2_11, location, true, true, true, transform.get()); - G2_22 = interpolateAndExtrapolate(G2_22, location, true, true, true, transform.get()); - G2_33 = interpolateAndExtrapolate(G2_33, location, true, true, true, transform.get()); - G2_12 = interpolateAndExtrapolate(G2_12, location, true, true, true, transform.get()); - G2_13 = interpolateAndExtrapolate(G2_13, location, true, true, true, transform.get()); - G2_23 = interpolateAndExtrapolate(G2_23, location, true, true, true, transform.get()); - - G3_11 = interpolateAndExtrapolate(G3_11, location, true, true, true, transform.get()); - G3_22 = interpolateAndExtrapolate(G3_22, location, true, true, true, transform.get()); - G3_33 = interpolateAndExtrapolate(G3_33, location, true, true, true, transform.get()); - G3_12 = interpolateAndExtrapolate(G3_12, location, true, true, true, transform.get()); - G3_13 = interpolateAndExtrapolate(G3_13, location, true, true, true, transform.get()); - G3_23 = interpolateAndExtrapolate(G3_23, location, true, true, true, transform.get()); - - G1 = interpolateAndExtrapolate(G1, location, true, true, true, transform.get()); - G2 = interpolateAndExtrapolate(G2, location, true, true, true, transform.get()); - G3 = interpolateAndExtrapolate(G3, location, true, true, true, transform.get()); + if (!g11.isFci()) { + G1_11 = 0.5 * g11 * DDX(g_11) + g12 * (DDX(g_12) - 0.5 * DDY(g_11)) + + g13 * (DDX(g_13) - 0.5 * DDZ(g_11)); + G1_22 = g11 * (DDY(g_12) - 0.5 * DDX(g_22)) + 0.5 * g12 * DDY(g_22) + + g13 * (DDY(g_23) - 0.5 * DDZ(g_22)); + G1_33 = g11 * (DDZ(g_13) - 0.5 * DDX(g_33)) + g12 * (DDZ(g_23) - 0.5 * DDY(g_33)) + + 0.5 * g13 * DDZ(g_33); + G1_12 = 0.5 * g11 * DDY(g_11) + 0.5 * g12 * DDX(g_22) + + 0.5 * g13 * (DDY(g_13) + DDX(g_23) - DDZ(g_12)); + G1_13 = 0.5 * g11 * DDZ(g_11) + 0.5 * g12 * (DDZ(g_12) + DDX(g_23) - DDY(g_13)) + + 0.5 * g13 * DDX(g_33); + G1_23 = 0.5 * g11 * (DDZ(g_12) + DDY(g_13) - DDX(g_23)) + + 0.5 * g12 * (DDZ(g_22) + DDY(g_23) - DDY(g_23)) + // + 0.5 *g13*(DDZ(g_32) + DDY(g_33) - DDZ(g_23)); + // which equals + + 0.5 * g13 * DDY(g_33); + + G2_11 = 0.5 * g12 * DDX(g_11) + g22 * (DDX(g_12) - 0.5 * DDY(g_11)) + + g23 * (DDX(g_13) - 0.5 * DDZ(g_11)); + G2_22 = g12 * (DDY(g_12) - 0.5 * DDX(g_22)) + 0.5 * g22 * DDY(g_22) + + g23 * (DDY(g23) - 0.5 * DDZ(g_22)); + G2_33 = g12 * (DDZ(g_13) - 0.5 * DDX(g_33)) + g22 * (DDZ(g_23) - 0.5 * DDY(g_33)) + + 0.5 * g23 * DDZ(g_33); + G2_12 = 0.5 * g12 * DDY(g_11) + 0.5 * g22 * DDX(g_22) + + 0.5 * g23 * (DDY(g_13) + DDX(g_23) - DDZ(g_12)); + G2_13 = + // 0.5 *g21*(DDZ(g_11) + DDX(g_13) - DDX(g_13)) + // which equals + 0.5 * g12 * (DDZ(g_11) + DDX(g_13) - DDX(g_13)) + // + 0.5 *g22*(DDZ(g_21) + DDX(g_23) - DDY(g_13)) + // which equals + + 0.5 * g22 * (DDZ(g_12) + DDX(g_23) - DDY(g_13)) + // + 0.5 *g23*(DDZ(g_31) + DDX(g_33) - DDZ(g_13)); + // which equals + + 0.5 * g23 * DDX(g_33); + G2_23 = 0.5 * g12 * (DDZ(g_12) + DDY(g_13) - DDX(g_23)) + 0.5 * g22 * DDZ(g_22) + + 0.5 * g23 * DDY(g_33); + + G3_11 = 0.5 * g13 * DDX(g_11) + g23 * (DDX(g_12) - 0.5 * DDY(g_11)) + + g33 * (DDX(g_13) - 0.5 * DDZ(g_11)); + G3_22 = g13 * (DDY(g_12) - 0.5 * DDX(g_22)) + 0.5 * g23 * DDY(g_22) + + g33 * (DDY(g_23) - 0.5 * DDZ(g_22)); + G3_33 = g13 * (DDZ(g_13) - 0.5 * DDX(g_33)) + g23 * (DDZ(g_23) - 0.5 * DDY(g_33)) + + 0.5 * g33 * DDZ(g_33); + G3_12 = + // 0.5 *g31*(DDY(g_11) + DDX(g_12) - DDX(g_12)) + // which equals to + 0.5 * g13 * DDY(g_11) + // + 0.5 *g32*(DDY(g_21) + DDX(g_22) - DDY(g_12)) + // which equals to + + 0.5 * g23 * DDX(g_22) + //+ 0.5 *g33*(DDY(g_31) + DDX(g_32) - DDZ(g_12)); + // which equals to + + 0.5 * g33 * (DDY(g_13) + DDX(g_23) - DDZ(g_12)); + G3_13 = 0.5 * g13 * DDZ(g_11) + 0.5 * g23 * (DDZ(g_12) + DDX(g_23) - DDY(g_13)) + + 0.5 * g33 * DDX(g_33); + G3_23 = 0.5 * g13 * (DDZ(g_12) + DDY(g_13) - DDX(g_23)) + 0.5 * g23 * DDZ(g_22) + + 0.5 * g33 * DDY(g_33); + + G1 = (DDX(J * g11) + DDY(J.asField3DParallel() * g12) + DDZ(J * g13)) / J; + G2 = (DDX(J * g12) + DDY(J.asField3DParallel() * g22) + DDZ(J * g23)) / J; + G3 = (DDX(J * g13) + DDY(J.asField3DParallel() * g23) + DDZ(J * g33)) / J; + + // Communicate christoffel symbol terms + output_progress.write("\tCommunicating connection terms\n"); + + localmesh->communicate_no_slices(G1_11, G1_22, G1_33, G1_12, G1_13, G1_23, G2_11, + G2_22, G2_33, G2_12, G2_13, G2_23, G3_11, G3_22, + G3_33, G3_12, G3_13, G3_23, G1, G2, G3); + + // Set boundary guard cells of Christoffel symbol terms + // Ideally, when location is staggered, we would set the upper/outer boundary point + // correctly rather than by extrapolating here: e.g. if location==CELL_YLOW and we are + // at the upper y-boundary the x- and z-derivatives at yend+1 at the boundary can be + // calculated because the guard cells are available, while the y-derivative could be + // calculated from the CELL_CENTRE metric components (which have guard cells available + // past the boundary location). This would avoid the problem that the y-boundary on the + // CELL_YLOW grid is at a 'guard cell' location (yend+1). + // However, the above would require lots of special handling, so just extrapolate for + // now. + G1_11 = interpolateAndExtrapolate(G1_11, location, true, true, true, transform.get()); + G1_22 = interpolateAndExtrapolate(G1_22, location, true, true, true, transform.get()); + G1_33 = interpolateAndExtrapolate(G1_33, location, true, true, true, transform.get()); + G1_12 = interpolateAndExtrapolate(G1_12, location, true, true, true, transform.get()); + G1_13 = interpolateAndExtrapolate(G1_13, location, true, true, true, transform.get()); + G1_23 = interpolateAndExtrapolate(G1_23, location, true, true, true, transform.get()); + + G2_11 = interpolateAndExtrapolate(G2_11, location, true, true, true, transform.get()); + G2_22 = interpolateAndExtrapolate(G2_22, location, true, true, true, transform.get()); + G2_33 = interpolateAndExtrapolate(G2_33, location, true, true, true, transform.get()); + G2_12 = interpolateAndExtrapolate(G2_12, location, true, true, true, transform.get()); + G2_13 = interpolateAndExtrapolate(G2_13, location, true, true, true, transform.get()); + G2_23 = interpolateAndExtrapolate(G2_23, location, true, true, true, transform.get()); + + G3_11 = interpolateAndExtrapolate(G3_11, location, true, true, true, transform.get()); + G3_22 = interpolateAndExtrapolate(G3_22, location, true, true, true, transform.get()); + G3_33 = interpolateAndExtrapolate(G3_33, location, true, true, true, transform.get()); + G3_12 = interpolateAndExtrapolate(G3_12, location, true, true, true, transform.get()); + G3_13 = interpolateAndExtrapolate(G3_13, location, true, true, true, transform.get()); + G3_23 = interpolateAndExtrapolate(G3_23, location, true, true, true, transform.get()); + + G1 = interpolateAndExtrapolate(G1, location, true, true, true, transform.get()); + G2 = interpolateAndExtrapolate(G2, location, true, true, true, transform.get()); + G3 = interpolateAndExtrapolate(G3, location, true, true, true, transform.get()); + } else { + G1_11 = G1_22 = G1_33 = G1_12 = G1_13 = G1_23 = + + G2_11 = G2_22 = G2_33 = G2_12 = G2_13 = G2_23 = + + G3_11 = G3_22 = G3_33 = G3_12 = G3_13 = G3_23 = + + G1 = G2 = G3 = BoutNaN; + } ////////////////////////////////////////////////////// /// Non-uniform meshes. Need to use DDX, DDY @@ -1102,7 +1111,7 @@ int Coordinates::geometry(bool recalculate_staggered, if (localmesh->get(d2y, "d2y" + suffix, 0.0, false, location)) { output_warn.write( "\tWARNING: differencing quantity 'd2y' not found. Calculating from dy\n"); - d1_dy = DDY(1. / dy); // d/di(1/dy) + d1_dy = DDY(1. / dy.asField3DParallel()); // d/di(1/dy) localmesh->communicate_no_slices(d1_dy); d1_dy = @@ -1485,16 +1494,8 @@ Coordinates::FieldMetric Coordinates::DDY(const Field2D& f, CELL_LOC loc, return bout::derivatives::index::DDY(f, loc, method, region) / dy; } -Field3D Coordinates::DDY(const Field3D& f, CELL_LOC outloc, const std::string& method, - const std::string& region) const { -#if BOUT_USE_METRIC_3D - if (!f.hasParallelSlices() and !transform->canToFromFieldAligned()) { - Field3D f_parallel = f; - transform->calcParallelSlices(f_parallel); - f_parallel.applyParallelBoundary("parallel_neumann_o2"); - return bout::derivatives::index::DDY(f_parallel, outloc, method, region); - } -#endif +Field3D Coordinates::DDY(const Field3DParallel& f, CELL_LOC outloc, + const std::string& method, const std::string& region) const { return bout::derivatives::index::DDY(f, outloc, method, region) / dy; }; @@ -1526,7 +1527,7 @@ Coordinates::FieldMetric Coordinates::Grad_par(const Field2D& var, return DDY(var) * invSg(); } -Field3D Coordinates::Grad_par(const Field3D& var, CELL_LOC outloc, +Field3D Coordinates::Grad_par(const Field3DParallel& var, CELL_LOC outloc, const std::string& method) { ASSERT1(location == outloc || outloc == CELL_DEFAULT); @@ -1546,8 +1547,8 @@ Coordinates::FieldMetric Coordinates::Vpar_Grad_par(const Field2D& v, const Fiel return VDDY(v, f) * invSg(); } -Field3D Coordinates::Vpar_Grad_par(const Field3D& v, const Field3D& f, CELL_LOC outloc, - const std::string& method) { +Field3D Coordinates::Vpar_Grad_par(const Field3D& v, const Field3DParallel& f, + CELL_LOC outloc, const std::string& method) { ASSERT1(location == outloc || outloc == CELL_DEFAULT); return VDDY(v, f, outloc, method) * invSg(); @@ -1568,7 +1569,7 @@ Coordinates::FieldMetric Coordinates::Div_par(const Field2D& f, CELL_LOC outloc, return Bxy * Grad_par(f / Bxy_floc, outloc, method); } -Field3D Coordinates::Div_par(const Field3D& f, CELL_LOC outloc, +Field3D Coordinates::Div_par(const Field3DParallel& f, CELL_LOC outloc, const std::string& method) { ASSERT1(location == outloc || outloc == CELL_DEFAULT); @@ -1577,20 +1578,7 @@ Field3D Coordinates::Div_par(const Field3D& f, CELL_LOC outloc, // Coordinates object const auto& Bxy_floc = f.getCoordinates()->Bxy; - if (!f.hasParallelSlices()) { - // No yup/ydown fields. The Grad_par operator will - // shift to field aligned coordinates - return Bxy * Grad_par(f / Bxy_floc, outloc, method); - } - - // Need to modify yup and ydown fields - Field3D f_B = f / Bxy_floc; - f_B.splitParallelSlices(); - for (int i = 0; i < f.getMesh()->ystart; ++i) { - f_B.yup(i) = f.yup(i) / Bxy_floc.yup(i); - f_B.ydown(i) = f.ydown(i) / Bxy_floc.ydown(i); - } - return Bxy * Grad_par(f_B, outloc, method); + return Bxy * Grad_par(f / Bxy_floc, outloc, method); } ///////////////////////////////////////////////////////// @@ -1608,7 +1596,7 @@ Coordinates::FieldMetric Coordinates::Grad2_par2(const Field2D& f, CELL_LOC outl return result; } -Field3D Coordinates::Grad2_par2(const Field3D& f, CELL_LOC outloc, +Field3D Coordinates::Grad2_par2(const Field3DParallel& f, CELL_LOC outloc, const std::string& method) { if (outloc == CELL_DEFAULT) { @@ -1774,9 +1762,10 @@ Coordinates::FieldMetric Coordinates::Laplace_par(const Field2D& f, CELL_LOC out return D2DY2(f, outloc) / g_22 + DDY(J / g_22, outloc) * DDY(f, outloc) / J; } -Field3D Coordinates::Laplace_par(const Field3D& f, CELL_LOC outloc) { +Field3D Coordinates::Laplace_par(const Field3DParallel& f, CELL_LOC outloc) { ASSERT1(location == outloc || outloc == CELL_DEFAULT); - return D2DY2(f, outloc) / g_22 + DDY(J / g_22, outloc) * ::DDY(f, outloc) / J; + return D2DY2(f, outloc) / g_22 + + DDY(J.asField3DParallel() / g_22, outloc) * ::DDY(f, outloc) / J; } // Full Laplacian operator on scalar field @@ -1796,7 +1785,7 @@ Coordinates::FieldMetric Coordinates::Laplace(const Field2D& f, CELL_LOC outloc, return result; } -Field3D Coordinates::Laplace(const Field3D& f, CELL_LOC outloc, +Field3D Coordinates::Laplace(const Field3DParallel& f, CELL_LOC outloc, const std::string& dfdy_boundary_conditions, const std::string& dfdy_dy_region) { diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 09433b0685..b0b5cc60d5 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -25,19 +25,20 @@ #include "bout/build_defines.hxx" -#include -#include -#include -#include -#include -#include -#include -#include - -#include // Delp2 uses same coefficients as inversion code - -#include -#include +#include "bout/assert.hxx" +#include "bout/derivs.hxx" +#include "bout/difops.hxx" +#include "bout/fft.hxx" +#include "bout/field2d.hxx" +#include "bout/globals.hxx" +#include "bout/interpolation.hxx" +#include "bout/invert_laplace.hxx" // Delp2 uses same coefficients as inversion code +#include "bout/msg_stack.hxx" +#include "bout/region.hxx" +#include "bout/solver.hxx" +#include "bout/unused.hxx" +#include "bout/utils.hxx" +#include "bout/vecops.hxx" #include @@ -366,6 +367,110 @@ Field3D Div_par_K_Grad_par(const Field3D& kY, const Field3D& f, CELL_LOC outloc) + Div_par(kY, outloc) * Grad_par(f, outloc); } +Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, Field3D& flow_ylow, + bool bndry_flux) { + TRACE("FV::Div_par_K_Grad_par_mod"); + + ASSERT2(Kin.getLocation() == fin.getLocation()); + + const Mesh* mesh = Kin.getMesh(); + const Coordinates* coord = fin.getCoordinates(); + + if (Kin.isFci()) { + ASSERT1(Kin.hasParallelSlices()); + ASSERT1(fin.hasParallelSlices()); + // Using parallel slices. + // Note: Y slices may use different coordinate systems + // -> Only B, dy and g_22 can be used in yup/ydown + // Others (e.g J) may not be averaged between y planes. + + const auto& K_up = Kin.yup(); + const auto& K_down = Kin.ydown(); + + const auto& f_up = fin.yup(); + const auto& f_down = fin.ydown(); + + Field3D result{zeroFrom(fin)}; + flow_ylow = zeroFrom(fin); + + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + const auto iyp = i.yp(); + const auto iym = i.ym(); + + // Upper cell edge + const BoutReal c_up = 0.5 * (Kin[i] + K_up[iyp]); // K at the upper boundary + const BoutReal J_up = + 0.5 * (coord->J[i] + coord->J.yup()[iyp]); // Jacobian at boundary + const BoutReal g_22_up = 0.5 * (coord->g_22[i] + coord->g_22.yup()[iyp]); + const BoutReal gradient_up = + 2. * (f_up[iyp] - fin[i]) / (coord->dy[i] + coord->dy.yup()[iyp]); + + const BoutReal flux_up = c_up * J_up * gradient_up / g_22_up; + + // Lower cell edge + const BoutReal c_down = 0.5 * (Kin[i] + K_down[iym]); // K at the lower boundary + const BoutReal J_down = + 0.5 * (coord->J[i] + coord->J.ydown()[iym]); // Jacobian at boundary + const BoutReal g_22_down = 0.5 * (coord->g_22[i] + coord->g_22.ydown()[iym]); + const BoutReal gradient_down = + 2. * (fin[i] - f_down[iym]) / (coord->dy[i] + coord->dy.ydown()[iym]); + + const BoutReal flux_down = c_down * J_down * gradient_down / g_22_down; + + result[i] = (flux_up - flux_down) / (coord->dy[i] * coord->J[i]); + } + + return result; + } + + // Calculate in field-aligned coordinates + const auto& K = toFieldAligned(Kin, "RGN_NOX"); + const auto& f = toFieldAligned(fin, "RGN_NOX"); + + Field3D result{zeroFrom(f)}; + flow_ylow = zeroFrom(f); + + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + // Calculate flux at upper surface + const auto ix = i.x(); + const auto iy = i.y(); + const auto iyp = i.yp(); + const auto iym = i.ym(); + + const bool is_periodic_y = mesh->periodicY(ix); + + if (bndry_flux || is_periodic_y || !mesh->lastY(ix) || (iy != mesh->yend)) { + const BoutReal c = 0.5 * (K[i] + K[iyp]); // K at the upper boundary + const BoutReal J = 0.5 * (coord->J[i] + coord->J[iyp]); // Jacobian at boundary + const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); + const BoutReal gradient = 2. * (f[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); + + const BoutReal flux = c * J * gradient / g_22; + + result[i] += flux / (coord->dy[i] * coord->J[i]); + } + + // Calculate flux at lower surface + if (bndry_flux || is_periodic_y || !mesh->firstY(ix) || (iy != mesh->ystart)) { + const BoutReal c = 0.5 * (K[i] + K[iym]); // K at the lower boundary + const BoutReal J = 0.5 * (coord->J[i] + coord->J[iym]); // Jacobian at boundary + const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); + const BoutReal gradient = 2. * (f[i] - f[iym]) / (coord->dy[i] + coord->dy[iym]); + + const BoutReal flux = c * J * gradient / g_22; + + result[i] -= flux / (coord->dy[i] * coord->J[i]); + flow_ylow[i] = -flux * coord->dx[i] * coord->dz[i]; + } + } + + // Shifted to field aligned coordinates, so need to shift back + result = fromFieldAligned(result, "RGN_NOBNDRY"); + flow_ylow = fromFieldAligned(flow_ylow); + + return result; +} + /******************************************************************************* * Delp2 * perpendicular Laplacian operator diff --git a/src/mesh/fv_ops.cxx b/src/mesh/fv_ops.cxx index fab8beb794..6b8d8a6f21 100644 --- a/src/mesh/fv_ops.cxx +++ b/src/mesh/fv_ops.cxx @@ -1,7 +1,16 @@ -#include -#include -#include -#include +#include "bout/fv_ops.hxx" + +#include "bout/assert.hxx" +#include "bout/bout_types.hxx" +#include "bout/boutexception.hxx" +#include "bout/build_config.hxx" +#include "bout/coordinates.hxx" +#include "bout/field2d.hxx" +#include "bout/field3d.hxx" +#include "bout/globals.hxx" +#include "bout/msg_stack.hxx" +#include "bout/region.hxx" +#include "bout/utils.hxx" namespace { template @@ -33,28 +42,19 @@ Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f) { // Flux in x - int xs = mesh->xstart - 1; - int xe = mesh->xend; - - /* - if(mesh->firstX()) - xs += 1; - */ - /* - if(mesh->lastX()) - xe -= 1; - */ + const int xs = mesh->xstart - 1; + const int xe = mesh->xend; for (int i = xs; i <= xe; i++) { for (int j = mesh->ystart; j <= mesh->yend; j++) { for (int k = mesh->zstart; k <= mesh->zend; k++) { // Calculate flux from i to i+1 - BoutReal fout = 0.5 * (a(i, j, k) + a(i + 1, j, k)) - * (coord->J(i, j, k) * coord->g11(i, j, k) - + coord->J(i + 1, j, k) * coord->g11(i + 1, j, k)) - * (f(i + 1, j, k) - f(i, j, k)) - / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); + const BoutReal fout = 0.5 * (a(i, j, k) + a(i + 1, j, k)) + * (coord->J(i, j, k) * coord->g11(i, j, k) + + coord->J(i + 1, j, k) * coord->g11(i + 1, j, k)) + * (f(i + 1, j, k) - f(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); result(i, j, k) += fout / (coord->dx(i, j, k) * coord->J(i, j, k)); result(i + 1, j, k) -= fout / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); @@ -178,14 +178,13 @@ Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f) { return result; } -const Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, - bool bndry_flux) { +Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, bool bndry_flux) { ASSERT2(Kin.getLocation() == fin.getLocation()); - Mesh* mesh = Kin.getMesh(); + const Mesh* mesh = Kin.getMesh(); - bool use_parallel_slices = (Kin.hasParallelSlices() && fin.hasParallelSlices()); + const bool use_parallel_slices = (Kin.hasParallelSlices() && fin.hasParallelSlices()); const auto& K = use_parallel_slices ? Kin : toFieldAligned(Kin, "RGN_NOX"); const auto& f = use_parallel_slices ? fin : toFieldAligned(fin, "RGN_NOX"); @@ -209,13 +208,13 @@ const Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, if (bndry_flux || mesh->periodicY(i.x()) || !mesh->lastY(i.x()) || (i.y() != mesh->yend)) { - BoutReal c = 0.5 * (K[i] + Kup[iyp]); // K at the upper boundary - BoutReal J = 0.5 * (coord->J[i] + coord->J[iyp]); // Jacobian at boundary - BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); + const BoutReal c = 0.5 * (K[i] + Kup[iyp]); // K at the upper boundary + const BoutReal J = 0.5 * (coord->J[i] + coord->J[iyp]); // Jacobian at boundary + const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); - BoutReal gradient = 2. * (fup[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); + const BoutReal gradient = 2. * (fup[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); - BoutReal flux = c * J * gradient / g_22; + const BoutReal flux = c * J * gradient / g_22; result[i] += flux / (coord->dy[i] * coord->J[i]); } @@ -223,14 +222,15 @@ const Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, // Calculate flux at lower surface if (bndry_flux || mesh->periodicY(i.x()) || !mesh->firstY(i.x()) || (i.y() != mesh->ystart)) { - BoutReal c = 0.5 * (K[i] + Kdown[iym]); // K at the lower boundary - BoutReal J = 0.5 * (coord->J[i] + coord->J[iym]); // Jacobian at boundary + const BoutReal c = 0.5 * (K[i] + Kdown[iym]); // K at the lower boundary + const BoutReal J = 0.5 * (coord->J[i] + coord->J[iym]); // Jacobian at boundary - BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); + const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); - BoutReal gradient = 2. * (f[i] - fdown[iym]) / (coord->dy[i] + coord->dy[iym]); + const BoutReal gradient = + 2. * (f[i] - fdown[iym]) / (coord->dy[i] + coord->dy[iym]); - BoutReal flux = c * J * gradient / g_22; + const BoutReal flux = c * J * gradient / g_22; result[i] -= flux / (coord->dy[i] * coord->J[i]); } @@ -244,10 +244,10 @@ const Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, return result; } -const Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { +Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { ASSERT1_FIELDS_COMPATIBLE(d_in, f_in); - Mesh* mesh = d_in.getMesh(); + const Mesh* mesh = d_in.getMesh(); Coordinates* coord = f_in.getCoordinates(); @@ -263,9 +263,9 @@ const Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { for (int i = mesh->xstart; i <= mesh->xend; i++) { // Check for boundaries - bool yperiodic = mesh->periodicY(i); - bool has_upper_boundary = !yperiodic && mesh->lastY(i); - bool has_lower_boundary = !yperiodic && mesh->firstY(i); + const bool yperiodic = mesh->periodicY(i); + const bool has_upper_boundary = !yperiodic && mesh->lastY(i); + const bool has_lower_boundary = !yperiodic && mesh->firstY(i); // Always calculate fluxes at upper Y cell boundary const int ystart = @@ -281,15 +281,15 @@ const Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { for (int j = ystart; j <= yend; j++) { for (int k = mesh->zstart; k <= mesh->zend; k++) { - BoutReal dy3 = SQ(coord->dy(i, j, k)) * coord->dy(i, j, k); + const BoutReal dy3 = SQ(coord->dy(i, j, k)) * coord->dy(i, j, k); // 3rd derivative at upper boundary - BoutReal d3fdy3 = + const BoutReal d3fdy3 = (f(i, j + 2, k) - 3. * f(i, j + 1, k) + 3. * f(i, j, k) - f(i, j - 1, k)) / dy3; - BoutReal flux = 0.5 * (d(i, j, k) + d(i, j + 1, k)) - * (coord->J(i, j, k) + coord->J(i, j + 1, k)) * d3fdy3; + const BoutReal flux = 0.5 * (d(i, j, k) + d(i, j + 1, k)) + * (coord->J(i, j, k) + coord->J(i, j + 1, k)) * d3fdy3; result(i, j, k) += flux / (coord->J(i, j, k) * coord->dy(i, j, k)); result(i, j + 1, k) -= flux / (coord->J(i, j + 1, k) * coord->dy(i, j + 1, k)); @@ -301,8 +301,8 @@ const Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; } -const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { - Mesh* mesh = f_in.getMesh(); +Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { + const Mesh* mesh = f_in.getMesh(); // Convert to field aligned coordinates const bool is_unaligned = (f_in.getDirectionY() == YDirectionType::Standard); @@ -313,10 +313,10 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { Coordinates* coord = f_in.getCoordinates(); for (int i = mesh->xstart; i <= mesh->xend; i++) { - bool yperiodic = mesh->periodicY(i); + const bool yperiodic = mesh->periodicY(i); - bool has_upper_boundary = !yperiodic && mesh->lastY(i); - bool has_lower_boundary = !yperiodic && mesh->firstY(i); + const bool has_upper_boundary = !yperiodic && mesh->lastY(i); + const bool has_lower_boundary = !yperiodic && mesh->firstY(i); for (int j = mesh->ystart; j <= mesh->yend; j++) { @@ -341,8 +341,8 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { // Not on domain boundary // 3rd derivative at right cell boundary - const BoutReal d3fdx3 = - (f(i, j + 2, k) - 3. * f(i, j + 1, k) + 3. * f(i, j, k) - f(i, j - 1, k)); + const BoutReal d3fdx3 = (f(i, j + 2, k) - (3. * f(i, j + 1, k)) + + (3. * f(i, j, k)) - f(i, j - 1, k)); result(i, j, k) += d3fdx3 * factor_rc; result(i, j + 1, k) -= d3fdx3 * factor_rp; @@ -363,10 +363,10 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { common_factor / (coord->J(i, j + 1, k) * coord->dy(i, j + 1, k)); const BoutReal d3fdx3 = - -((16. / 5) * 0.5 * (f(i, j + 1, k) + f(i, j, k)) // Boundary value f_b - - 6. * f(i, j, k) // f_0 - + 4. * f(i, j - 1, k) // f_1 - - (6. / 5) * f(i, j - 2, k) // f_2 + -(((16. / 5) * 0.5 * (f(i, j + 1, k) + f(i, j, k))) // Boundary value f_b + - (6. * f(i, j, k)) // f_0 + + (4. * f(i, j - 1, k)) // f_1 + - ((6. / 5) * f(i, j - 2, k)) // f_2 ); result(i, j, k) += d3fdx3 * factor_rc; @@ -392,8 +392,8 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { common_factor / (coord->J(i, j - 1, k) * coord->dy(i, j - 1, k)); // Not on a domain boundary - const BoutReal d3fdx3 = - (f(i, j + 1, k) - 3. * f(i, j, k) + 3. * f(i, j - 1, k) - f(i, j - 2, k)); + const BoutReal d3fdx3 = (f(i, j + 1, k) - (3. * f(i, j, k)) + + (3. * f(i, j - 1, k)) - f(i, j - 2, k)); result(i, j, k) -= d3fdx3 * factor_lc; result(i, j - 1, k) += d3fdx3 * factor_lm; @@ -410,10 +410,10 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { const BoutReal factor_lm = common_factor / (coord->J(i, j - 1, k) * coord->dy(i, j - 1, k)); const BoutReal d3fdx3 = - -(-(16. / 5) * 0.5 * (f(i, j - 1, k) + f(i, j, k)) // Boundary value f_b - + 6. * f(i, j, k) // f_0 - - 4. * f(i, j + 1, k) // f_1 - + (6. / 5) * f(i, j + 2, k) // f_2 + -((-(16. / 5) * 0.5 * (f(i, j - 1, k) + f(i, j, k))) // Boundary value f_b + + (6. * f(i, j, k)) // f_0 + - (4. * f(i, j + 1, k)) // f_1 + + ((6. / 5) * f(i, j + 2, k)) // f_2 ); result(i, j, k) -= d3fdx3 * factor_lc; @@ -436,8 +436,9 @@ void communicateFluxes(Field3D& f) { throw BoutException("communicateFluxes: Sorry!"); } - int size = mesh->LocalNy * mesh->LocalNz; - comm_handle xin, xout; + const int size = mesh->LocalNy * mesh->LocalNz; + comm_handle xin = nullptr; + comm_handle xout = nullptr; // Cache results to silence spurious compiler warning about xin, // xout possibly being uninitialised when used const bool not_first = mesh->periodicX || !mesh->firstX(); @@ -496,45 +497,45 @@ Field3D Div_Perp_Lap(const Field3D& a, const Field3D& f, CELL_LOC outloc) { // o --- gD --- o // Coordinates* coords = a.getCoordinates(outloc); - Mesh* mesh = f.getMesh(); + const Mesh* mesh = f.getMesh(); for (int i = mesh->xstart; i <= mesh->xend; i++) { for (int j = mesh->ystart; j <= mesh->yend; j++) { for (int k = 0; k < mesh->LocalNz; k++) { // wrap k-index around as Z is (currently) periodic. - int kp = (k + 1) % (mesh->LocalNz); - int km = (k - 1 + mesh->LocalNz) % (mesh->LocalNz); + const int kp = (k + 1) % (mesh->LocalNz); + const int km = (k - 1 + mesh->LocalNz) % (mesh->LocalNz); // Calculate gradients on cell faces -- assumes constant grid spacing - BoutReal gR = - (coords->g11(i, j, k) + coords->g11(i + 1, j, k)) - * (f(i + 1, j, k) - f(i, j, k)) - / (coords->dx(i + 1, j, k) + coords->dx(i, j, k)) - + 0.5 * (coords->g13(i, j, k) + coords->g13(i + 1, j, k)) - * (f(i + 1, j, kp) - f(i + 1, j, km) + f(i, j, kp) - f(i, j, km)) - / (4. * coords->dz(i, j, k)); - - BoutReal gL = - (coords->g11(i - 1, j, k) + coords->g11(i, j, k)) - * (f(i, j, k) - f(i - 1, j, k)) - / (coords->dx(i - 1, j, k) + coords->dx(i, j, k)) - + 0.5 * (coords->g13(i - 1, j, k) + coords->g13(i, j, k)) - * (f(i - 1, j, kp) - f(i - 1, j, km) + f(i, j, kp) - f(i, j, km)) - / (4 * coords->dz(i, j, k)); - - BoutReal gD = - coords->g13(i, j, k) - * (f(i + 1, j, km) - f(i - 1, j, km) + f(i + 1, j, k) - f(i - 1, j, k)) - / (4. * coords->dx(i, j, k)) - + coords->g33(i, j, k) * (f(i, j, k) - f(i, j, km)) / coords->dz(i, j, k); - - BoutReal gU = - coords->g13(i, j, k) - * (f(i + 1, j, kp) - f(i - 1, j, kp) + f(i + 1, j, k) - f(i - 1, j, k)) - / (4. * coords->dx(i, j, k)) - + coords->g33(i, j, k) * (f(i, j, kp) - f(i, j, k)) / coords->dz(i, j, k); + const BoutReal gR = + ((coords->g11(i, j, k) + coords->g11(i + 1, j, k)) + * (f(i + 1, j, k) - f(i, j, k)) + / (coords->dx(i + 1, j, k) + coords->dx(i, j, k))) + + (0.5 * (coords->g13(i, j, k) + coords->g13(i + 1, j, k)) + * (f(i + 1, j, kp) - f(i + 1, j, km) + f(i, j, kp) - f(i, j, km)) + / (4. * coords->dz(i, j, k))); + + const BoutReal gL = + ((coords->g11(i - 1, j, k) + coords->g11(i, j, k)) + * (f(i, j, k) - f(i - 1, j, k)) + / (coords->dx(i - 1, j, k) + coords->dx(i, j, k))) + + (0.5 * (coords->g13(i - 1, j, k) + coords->g13(i, j, k)) + * (f(i - 1, j, kp) - f(i - 1, j, km) + f(i, j, kp) - f(i, j, km)) + / (4 * coords->dz(i, j, k))); + + const BoutReal gD = + (coords->g13(i, j, k) + * (f(i + 1, j, km) - f(i - 1, j, km) + f(i + 1, j, k) - f(i - 1, j, k)) + / (4. * coords->dx(i, j, k))) + + (coords->g33(i, j, k) * (f(i, j, k) - f(i, j, km)) / coords->dz(i, j, k)); + + const BoutReal gU = + (coords->g13(i, j, k) + * (f(i + 1, j, kp) - f(i - 1, j, kp) + f(i + 1, j, k) - f(i - 1, j, k)) + / (4. * coords->dx(i, j, k))) + + (coords->g33(i, j, k) * (f(i, j, kp) - f(i, j, k)) / coords->dz(i, j, k)); // Flow right BoutReal flux = gR * 0.25 * (coords->J(i + 1, j, k) + coords->J(i, j, k)) diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index 0aa1f2b066..5b9271c93f 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -120,8 +120,12 @@ XZHermiteSplineBase::XZHermiteSplineBase(int y_offset, Mesh if (options == nullptr) { options = &Options::root()["mesh:paralleltransform:xzinterpolation"]; } - abs_fac_monotonic = (*options)["abs_tol"].withDefault(abs_fac_monotonic); - rel_fac_monotonic = (*options)["rel_tol"].withDefault(rel_fac_monotonic); + abs_fac_monotonic = (*options)["atol"] + .doc("Absolute tolerance for clipping overshoot") + .withDefault(abs_fac_monotonic); + rel_fac_monotonic = (*options)["rtol"] + .doc("Relative tolerance for clipping overshoot") + .withDefault(rel_fac_monotonic); } // Index arrays contain guard cells in order to get subscripts right diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index fa965a335a..5b23ce7fc6 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -64,6 +64,56 @@ using namespace std::string_view_literals; +namespace { +// Get a unique name for a field based on the sign/magnitude of the offset +std::string parallel_slice_field_name(std::string field, int offset) { + const std::string direction = (offset > 0) ? "forward" : "backward"; + // We only have a suffix for parallel slices beyond the first + // This is for backwards compatibility + const std::string slice_suffix = + (std::abs(offset) > 1) ? "_" + std::to_string(std::abs(offset)) : ""; + return direction + "_" + field + slice_suffix; +}; + +#if BOUT_USE_METRIC_3D +void load_parallel_metric_component(std::string name, Field3D& component, int offset) { + Mesh* mesh = component.getMesh(); + Field3D tmp{mesh}; + const auto pname = parallel_slice_field_name(name, offset); + if (mesh->get(tmp, pname, 0.0, false) != 0) { + throw BoutException("Could not read {:s} from grid file!\n" + " Fix it up with `zoidberg-update-parallel-metrics `", + pname); + } + if (!component.hasParallelSlices()) { + component.splitParallelSlices(); + component.disallowCalcParallelSlices(); + component.resetRegionParallel(true); + } + auto& pcom = component.ynext(offset); + pcom.allocate(); + BOUT_FOR(i, component.getRegion("RGN_NOBNDRY")) { pcom[i.yp(offset)] = tmp[i]; } +} + +void load_parallel_metric_components(Coordinates* coords, int offset) { +#define LOAD_PAR(var) load_parallel_metric_component(#var, coords->var, offset) + LOAD_PAR(g11); + LOAD_PAR(g22); + LOAD_PAR(g33); + LOAD_PAR(g13); + LOAD_PAR(g_11); + LOAD_PAR(g_22); + LOAD_PAR(g_33); + LOAD_PAR(g_13); + LOAD_PAR(dy); + LOAD_PAR(Bxy); + +#undef LOAD_PAR +} +#endif + +} // namespace + FCIMap::FCIMap(Mesh& mesh, [[maybe_unused]] const Coordinates::FieldMetric& dy, Options& options, int offset, const std::shared_ptr& inner_boundary, @@ -102,38 +152,29 @@ FCIMap::FCIMap(Mesh& mesh, [[maybe_unused]] const Coordinates::FieldMetric& dy, map_mesh->get(R, "R", 0.0, false); map_mesh->get(Z, "Z", 0.0, false); - // Get a unique name for a field based on the sign/magnitude of the offset - const auto parallel_slice_field_name = [&](std::string_view field) -> std::string { - const auto direction = (offset_ > 0) ? "forward"sv : "backward"sv; - // We only have a suffix for parallel slices beyond the first - // This is for backwards compatibility - if (std::abs(offset_) == 1) { - return fmt::format("{}_{}", direction, field); - } - return fmt::format("{}_{}_{}", direction, field, std::abs(offset_)); - }; - // If we can't read in any of these fields, things will silently not // work, so best throw - if (map_mesh->get(xt_prime, parallel_slice_field_name("xt_prime"), 0.0, false) != 0) { + if (map_mesh->get(xt_prime, parallel_slice_field_name("xt_prime", offset), 0.0, false) + != 0) { throw BoutException("Could not read {:s} from grid file!\n" " Either add it to the grid file, or reduce MYG", - parallel_slice_field_name("xt_prime")); + parallel_slice_field_name("xt_prime", offset)); } - if (map_mesh->get(zt_prime, parallel_slice_field_name("zt_prime"), 0.0, false) != 0) { + if (map_mesh->get(zt_prime, parallel_slice_field_name("zt_prime", offset), 0.0, false) + != 0) { throw BoutException("Could not read {:s} from grid file!\n" " Either add it to the grid file, or reduce MYG", - parallel_slice_field_name("zt_prime")); + parallel_slice_field_name("zt_prime", offset)); } - if (map_mesh->get(R_prime, parallel_slice_field_name("R"), 0.0, false) != 0) { + if (map_mesh->get(R_prime, parallel_slice_field_name("R", offset), 0.0, false) != 0) { throw BoutException("Could not read {:s} from grid file!\n" " Either add it to the grid file, or reduce MYG", - parallel_slice_field_name("R")); + parallel_slice_field_name("R", offset)); } - if (map_mesh->get(Z_prime, parallel_slice_field_name("Z"), 0.0, false) != 0) { + if (map_mesh->get(Z_prime, parallel_slice_field_name("Z", offset), 0.0, false) != 0) { throw BoutException("Could not read {:s} from grid file!\n" " Either add it to the grid file, or reduce MYG", - parallel_slice_field_name("Z")); + parallel_slice_field_name("Z", offset)); } // Cell corners @@ -415,3 +456,24 @@ void FCITransform::outputVars(Options& output_options) { output_options["R"].force(R, "FCI"); output_options["Z"].force(Z, "FCI"); } + +void FCITransform::loadParallelMetrics(Coordinates* coords) { +#if BOUT_USE_METRIC_3D + output_info.write("\tLoading parallel metrics\n"); + const auto JB0 = coords->J * coords->Bxy; + coords->J.splitParallelSlices(); + coords->J.disallowCalcParallelSlices(); + coords->J.resetRegionParallel(true); + for (int i = 1; i <= mesh.ystart; ++i) { + load_parallel_metric_components(coords, -i); + load_parallel_metric_components(coords, i); + + coords->J.ynext(i).allocate(); + coords->J.ynext(-i).allocate(); + BOUT_FOR(j, JB0.getRegion("RGN_NOBNDRY")) { + coords->J.ynext(i)[j.yp(i)] = JB0[j] / coords->Bxy.ynext(i)[j.yp(i)]; + coords->J.ynext(-i)[j.yp(-i)] = JB0[j] / coords->Bxy.ynext(-i)[j.yp(-i)]; + } + } +#endif +} diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx index 65529a4c4e..71ac35a192 100644 --- a/src/mesh/parallel/fci.hxx +++ b/src/mesh/parallel/fci.hxx @@ -125,6 +125,8 @@ public: return false; } + void loadParallelMetrics(Coordinates* coords) override; + protected: void checkInputGrid() override; diff --git a/src/sys/derivs.cxx b/src/sys/derivs.cxx index aab75b8f19..2b606b9e5d 100644 --- a/src/sys/derivs.cxx +++ b/src/sys/derivs.cxx @@ -66,7 +66,7 @@ Coordinates::FieldMetric DDX(const Field2D& f, CELL_LOC outloc, const std::strin ////////////// Y DERIVATIVE ///////////////// -Field3D DDY(const Field3D& f, CELL_LOC outloc, const std::string& method, +Field3D DDY(const Field3DParallel& f, CELL_LOC outloc, const std::string& method, const std::string& region) { return bout::derivatives::index::DDY(f, outloc, method, region) / f.getCoordinates(outloc)->dy; @@ -406,7 +406,7 @@ Coordinates::FieldMetric VDDY(const Field2D& v, const Field2D& f, CELL_LOC outlo } // general case -Field3D VDDY(const Field3D& v, const Field3D& f, CELL_LOC outloc, +Field3D VDDY(const Field3D& v, const Field3DParallel& f, CELL_LOC outloc, const std::string& method, const std::string& region) { return bout::derivatives::index::VDDY(v, f, outloc, method, region) / f.getCoordinates(outloc)->dy; @@ -467,7 +467,7 @@ Coordinates::FieldMetric FDDY(const Field2D& v, const Field2D& f, CELL_LOC outlo / f.getCoordinates(outloc)->dy; } -Field3D FDDY(const Field3D& v, const Field3D& f, CELL_LOC outloc, +Field3D FDDY(const Field3D& v, const Field3DParallel& f, CELL_LOC outloc, const std::string& method, const std::string& region) { return bout::derivatives::index::FDDY(v, f, outloc, method, region) / f.getCoordinates(outloc)->dy; diff --git a/tests/MMS/CMakeLists.txt b/tests/MMS/CMakeLists.txt index a6667bcfa5..510385d54c 100644 --- a/tests/MMS/CMakeLists.txt +++ b/tests/MMS/CMakeLists.txt @@ -8,6 +8,7 @@ add_subdirectory(spatial/d2dx2) add_subdirectory(spatial/d2dz2) add_subdirectory(spatial/diffusion) add_subdirectory(spatial/fci) +add_subdirectory(spatial/finite-volume) add_subdirectory(time) add_subdirectory(time-petsc) add_subdirectory(wave-1d) diff --git a/tests/MMS/spatial/fci/CMakeLists.txt b/tests/MMS/spatial/fci/CMakeLists.txt index dcdb183fdd..48721a8199 100644 --- a/tests/MMS/spatial/fci/CMakeLists.txt +++ b/tests/MMS/spatial/fci/CMakeLists.txt @@ -3,5 +3,6 @@ bout_add_mms_test( SOURCES fci_mms.cxx USE_RUNTEST USE_DATA_BOUT_INP REQUIRES zoidberg_FOUND + REQUIRES BOUT_USE_METRIC_3D PROCESSORS 2 ) diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp index 93e2101473..76ac3035c9 100644 --- a/tests/MMS/spatial/fci/data/BOUT.inp +++ b/tests/MMS/spatial/fci/data/BOUT.inp @@ -5,6 +5,9 @@ div_par_solution = (0.01*x + 0.045)*(-12.5663706143592*cos(y - 2*z) - 6.28318530 div_par_K_grad_par_solution = (0.01*x + 0.045)*(6.28318530717959*sin(y - z) - 0.628318530717959*sin(y - z)/(0.01*x + 0.045))*(6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z))/((0.01*x + 0.045)^2 + 1.0) + (6.28318530717959*(0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0))*cos(y - z)/sqrt((0.01*x + 0.045)^2 + 1.0) K = cos(y - z) laplace_par_solution = (0.01*x + 0.045)*(6.28318530717959*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/((0.01*x + 0.045)*sqrt((0.01*x + 0.045)^2 + 1.0)))/sqrt((0.01*x + 0.045)^2 + 1.0) +FV_div_par_mod_solution = (0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*((sin(y - 2*z) + sin(y - z))*sin(y - z)/(0.01*x + 0.045) + (-2*cos(y - 2*z) - cos(y - z))*cos(y - z)/(0.01*x + 0.045)) - 0.628318530717959*(sin(y - 2*z) + sin(y - z))*sin(y - z)/(0.01*x + 0.045) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))*cos(y - z)/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0) +FV_div_par_fvv_solution = (0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(2*(sin(y - 2*z) + sin(y - z))*sin(y - z)*cos(y - z)/(0.01*x + 0.045) + (-2*cos(y - 2*z) - cos(y - z))*cos(y - z)^2/(0.01*x + 0.045)) - 1.25663706143592*(sin(y - 2*z) + sin(y - z))*sin(y - z)*cos(y - z)/(0.01*x + 0.045) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))*cos(y - z)^2/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0) +div_par_K_grad_par_mod_solution = (0.01*x + 0.045)*(6.28318530717959*sin(y - z) - 0.628318530717959*sin(y - z)/(0.01*x + 0.045))*(6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z))/((0.01*x + 0.045)^2 + 1.0) + (6.28318530717959*(0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0))*cos(y - z)/sqrt((0.01*x + 0.045)^2 + 1.0) [mesh] symmetricglobalx = true diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx index 7967452f3d..e75f857a5f 100644 --- a/tests/MMS/spatial/fci/fci_mms.cxx +++ b/tests/MMS/spatial/fci/fci_mms.cxx @@ -4,6 +4,7 @@ #include "bout/field.hxx" #include "bout/field3d.hxx" #include "bout/field_factory.hxx" +#include "bout/fv_ops.hxx" #include "bout/globals.hxx" #include "bout/options.hxx" #include "bout/options_io.hxx" @@ -46,24 +47,27 @@ int main(int argc, char** argv) { Field3D K{FieldFactory::get()->create3D("K", Options::getRoot(), mesh)}; // Communicate to calculate parallel transform. - if constexpr (bout::build::use_metric_3d) { - // Div_par operators require B parallel slices: - // Coordinates::geometry doesn't ensure this (yet) - auto& Bxy = mesh->getCoordinates()->Bxy; - mesh->communicate(Bxy); - } mesh->communicate(input, K); Options dump; // Add mesh geometry variables mesh->outputVars(dump); + // Dummy variable for *_mod overloads + Field3D flow_ylow; + fci_op_test("grad_par", dump, input, Grad_par(input)); fci_op_test("grad2_par2", dump, input, Grad2_par2(input)); fci_op_test("div_par", dump, input, Div_par(input)); fci_op_test("div_par_K_grad_par", dump, input, Div_par_K_Grad_par(K, input)); + fci_op_test("div_par_K_grad_par_mod", dump, input, + Div_par_K_Grad_par_mod(K, input, flow_ylow)); fci_op_test("laplace_par", dump, input, Laplace_par(input)); + // Finite volume methods + fci_op_test("FV_div_par_mod", dump, input, FV::Div_par_mod(input, K, K, flow_ylow)); + fci_op_test("FV_div_par_fvv", dump, input, FV::Div_par_fvv(input, K, K)); + bout::writeDefaultOutputFile(dump); BoutFinalise(); diff --git a/tests/MMS/spatial/fci/mms.py b/tests/MMS/spatial/fci/mms.py index b28e337ac0..801a8d3f26 100755 --- a/tests/MMS/spatial/fci/mms.py +++ b/tests/MMS/spatial/fci/mms.py @@ -16,6 +16,7 @@ f = sin(y - z) + sin(y - 2 * z) K = cos(z - y) + Lx = 0.1 Ly = 10.0 Lz = 1.0 @@ -60,7 +61,10 @@ def FCI_Laplace_par(f: Expr) -> Expr: ("grad2_par2_solution", FCI_grad2_par2(f)), ("div_par_solution", FCI_div_par(f)), ("div_par_K_grad_par_solution", FCI_div_par_K_grad_par(f, K)), + ("div_par_K_grad_par_mod_solution", FCI_div_par_K_grad_par(f, K)), ("laplace_par_solution", FCI_Laplace_par(f)), + ("FV_div_par_mod_solution", FCI_div_par(f * K)), + ("FV_div_par_fvv_solution", FCI_div_par(f * K * K)), ): expr_str = exprToStr(expr) print(f"{name} = {expr_str}") diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 9dd4d1817f..ed22dd70ce 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -22,7 +22,21 @@ from numpy import array, linspace, log, polyfit DIRECTORY = "data" NPROC = 2 MTHREAD = 2 -OPERATORS = ("grad_par", "grad2_par2", "div_par", "div_par_K_grad_par", "laplace_par") +OPERATORS = ( + "grad_par", + "grad2_par2", + "div_par", + "div_par_K_grad_par", + "div_par_K_grad_par_mod", + "laplace_par", + "FV_div_par_mod", + "FV_div_par_fvv", +) + +# div_par_fvv is also tested in ../finite-volume, where 1.5th order is achieved +operator_order = { + "FV_div_par_fvv": 1, +} # Note that we need at least _2_ interior points for hermite spline # interpolation due to an awkwardness with the boundaries NX = 4 @@ -165,7 +179,10 @@ def check_fci_operators(name: str, case: dict) -> bool: for operator in OPERATORS: test_name = f"{operator} {name}" success = assert_convergence( - final_errors[operator]["l_2"], dx, test_name, order + final_errors[operator]["l_2"], + dx, + test_name, + operator_order.get(operator, order), ) if not success: failures.append(test_name) diff --git a/tests/MMS/spatial/finite-volume/CMakeLists.txt b/tests/MMS/spatial/finite-volume/CMakeLists.txt new file mode 100644 index 0000000000..7180eb7d98 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/CMakeLists.txt @@ -0,0 +1,6 @@ +bout_add_mms_test( + MMS-spatial-finite-volume + SOURCES fv_mms.cxx + USE_RUNTEST USE_DATA_BOUT_INP + PROCESSORS 2 +) diff --git a/tests/MMS/spatial/finite-volume/data/BOUT.inp b/tests/MMS/spatial/finite-volume/data/BOUT.inp new file mode 100644 index 0000000000..029011e437 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/data/BOUT.inp @@ -0,0 +1,23 @@ +input_field = 0.1*sin(2.0*y) + 1 +FV_Div_par_mod_solution = -0.188495559215388*sin(3.0*y) +FV_Div_par_fvv_solution = -0.376991118430775*(0.1*cos(3.0*y) + 1)*sin(3.0*y)/(0.1*sin(2.0*y) + 1) - 0.125663706143592*(0.1*cos(3.0*y) + 1)^2*cos(2.0*y)/(0.1*sin(2.0*y) + 1)^2 +FV_Div_par_solution = -0.188495559215388*sin(3.0*y) +FV_Div_par_K_Grad_par_solution = 0.125663706143592*(-0.188495559215388*sin(3.0*y)/(0.1*sin(2.0*y) + 1) - 0.125663706143592*(0.1*cos(3.0*y) + 1)*cos(2.0*y)/(0.1*sin(2.0*y) + 1)^2)*cos(2.0*y) - 0.15791367041743*(0.1*cos(3.0*y) + 1)*sin(2.0*y)/(0.1*sin(2.0*y) + 1) +FV_Div_par_K_Grad_par_mod_solution = 0.125663706143592*(-0.188495559215388*sin(3.0*y)/(0.1*sin(2.0*y) + 1) - 0.125663706143592*(0.1*cos(3.0*y) + 1)*cos(2.0*y)/(0.1*sin(2.0*y) + 1)^2)*cos(2.0*y) - 0.15791367041743*(0.1*cos(3.0*y) + 1)*sin(2.0*y)/(0.1*sin(2.0*y) + 1) +v = (0.1*cos(3.0*y) + 1)/(0.1*sin(2.0*y) + 1) + +[mesh] +MXG = 0 + +nx = 1 +ny = 128 +nz = 1 + +Ly = 10 + +dy = Ly / ny +J = 1 # Identity metric + +[mesh:ddy] +first = C2 +second = C2 diff --git a/tests/MMS/spatial/finite-volume/fv_mms.cxx b/tests/MMS/spatial/finite-volume/fv_mms.cxx new file mode 100644 index 0000000000..19f3f14610 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/fv_mms.cxx @@ -0,0 +1,92 @@ +#include "bout/bout.hxx" +#include "bout/difops.hxx" +#include "bout/field.hxx" +#include "bout/field3d.hxx" +#include "bout/field_factory.hxx" +#include "bout/fv_ops.hxx" +#include "bout/globals.hxx" +#include "bout/options.hxx" +#include "bout/options_io.hxx" +#include "bout/utils.hxx" + +#include + +#include +#include + +namespace { +auto fv_op_test(const std::string& name, Options& dump, const Field3D& input, + const Field3D& result, std::string suffix = "") { + auto* mesh = input.getMesh(); + const Field3D solution{FieldFactory::get()->create3D(fmt::format("{}_solution", name), + Options::getRoot(), mesh)}; + const Field3D error{result - solution}; + + dump[fmt::format("{}{}_l_2", name, suffix)] = + sqrt(mean(SQ(error), true, "RGN_NOBNDRY")); + dump[fmt::format("{}{}_l_inf", name, suffix)] = max(abs(error), true, "RGN_NOBNDRY"); + + dump[fmt::format("{}{}_result", name, suffix)] = result; + dump[fmt::format("{}{}_error", name, suffix)] = error; + dump[fmt::format("{}{}_input", name, suffix)] = input; + dump[fmt::format("{}{}_solution", name, suffix)] = solution; +} +} // namespace + +int main(int argc, char** argv) { + BoutInitialise(argc, argv); + + using bout::globals::mesh; + + Field3D input{FieldFactory::get()->create3D("input_field", Options::getRoot(), mesh)}; + Field3D v{FieldFactory::get()->create3D("v", Options::getRoot(), mesh)}; + + // Communicate to calculate parallel transform. + mesh->communicate(input, v); + + Options dump; + // Add mesh geometry variables + mesh->outputVars(dump); + dump["v"] = v; + + // Dummy variable for *_mod overloads + Field3D flow_ylow; + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), "_MC"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_MC"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), "_MC"); + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), "_Upwind"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_Upwind"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), + "_Upwind"); + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), "_Fromm"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_Fromm"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), + "_Fromm"); + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), "_MinMod"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_MinMod"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), + "_MinMod"); + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), + "_Superbee"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_Superbee"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), + "_Superbee"); + + fv_op_test("FV_Div_par_K_Grad_par", dump, input, FV::Div_par_K_Grad_par(v, input)); + fv_op_test("FV_Div_par_K_Grad_par_mod", dump, input, + Div_par_K_Grad_par_mod(v, input, flow_ylow)); + + bout::writeDefaultOutputFile(dump); + + BoutFinalise(); +} diff --git a/tests/MMS/spatial/finite-volume/makefile b/tests/MMS/spatial/finite-volume/makefile new file mode 100644 index 0000000000..88ba6c77e7 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/makefile @@ -0,0 +1,6 @@ + +BOUT_TOP = ../../../.. + +SOURCEC = fci_mms.cxx + +include $(BOUT_TOP)/make.config diff --git a/tests/MMS/spatial/finite-volume/mms.py b/tests/MMS/spatial/finite-volume/mms.py new file mode 100755 index 0000000000..a95ecc1328 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/mms.py @@ -0,0 +1,64 @@ +#!/usr/bin/env python3 +# +# Generate manufactured solution and sources for FCI test +# + +from math import pi +import warnings + +from boututils.boutwarnings import AlwaysWarning +from boutdata.data import BoutOptionsFile +from boutdata.mms import exprToStr, y, Grad_par, Div_par, Metric +from sympy import sin, cos, Expr + +warnings.simplefilter("ignore", AlwaysWarning) + +# Length of the y domain +Ly = 10.0 + +# Identity +metric = Metric() + +# Define solution in terms of input x,y,z +f = 1 + 0.1 * sin(2 * y) +fv = 1 + 0.1 * cos(3 * y) + + +# Turn solution into real x and z coordinates +replace = [(y, metric.y * 2 * pi / Ly)] + +f = f.subs(replace) +fv = fv.subs(replace) +v = fv / f + +# Substitute back to get input y coordinates +replace = [(metric.y, y * Ly / (2 * pi))] + + +def Grad2_par2(f: Expr) -> Expr: + return Grad_par(Grad_par(f)) + + +def Div_par_K_Grad_par(f: Expr, K: Expr) -> Expr: + return (K * Grad2_par2(f)) + (Div_par(K) * Grad_par(f)) + + +############################################ +# Equations solved + +options = BoutOptionsFile("data/BOUT.inp") + +for name, expr in ( + ("input_field", f), + ("v", v), + ("FV_Div_par_solution", Div_par(f * v)), + ("FV_Div_par_K_Grad_par_solution", Div_par_K_Grad_par(f, v)), + ("FV_Div_par_K_Grad_par_mod_solution", Div_par_K_Grad_par(f, v)), + ("FV_Div_par_mod_solution", Div_par(f * v)), + ("FV_Div_par_fvv_solution", Div_par(f * v * v)), +): + expr_str = exprToStr(expr.subs(replace)) + print(f"{name} = {expr_str}") + options[name] = expr_str + +options.write("data/BOUT.inp", overwrite=True) diff --git a/tests/MMS/spatial/finite-volume/runtest b/tests/MMS/spatial/finite-volume/runtest new file mode 100755 index 0000000000..bcd4672545 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/runtest @@ -0,0 +1,201 @@ +#!/usr/bin/env python3 +# +# Python script to run and analyse MMS test +# + +import argparse +import json +import pathlib +import sys +from time import time + +from boutdata.collect import collect +from boututils.run_wrapper import build_and_log, launch_safe +from numpy import array, log, polyfit + +# Global parameters +DIRECTORY = "data" +NPROC = 2 +MTHREAD = 2 +OPERATORS = { + # Slope-limiters necessarily reduce the accuracy in places + "FV_Div_par_MC": 1.5, + "FV_Div_par_mod_MC": 1.5, + "FV_Div_par_fvv_MC": 1.5, + "FV_Div_par_Upwind": 1, + "FV_Div_par_mod_Upwind": 1, + "FV_Div_par_fvv_Upwind": 1, + "FV_Div_par_Fromm": 1.5, + "FV_Div_par_mod_Fromm": 1.5, + "FV_Div_par_fvv_Fromm": 1.5, + "FV_Div_par_MinMod": 1.5, + "FV_Div_par_mod_MinMod": 1.5, + "FV_Div_par_fvv_MinMod": 1.5, + "FV_Div_par_Superbee": 1.5, + "FV_Div_par_mod_Superbee": 1.5, + "FV_Div_par_fvv_Superbee": 1.5, + "FV_Div_par_K_Grad_par": 2, + "FV_Div_par_K_Grad_par_mod": 2, +} +# Resolution in y and z +NLIST = [8, 16, 32, 64] +dx = 1.0 / array(NLIST) + + +def quiet_collect(name: str) -> float: + # Index to return a plain (numpy) float rather than `BoutArray` + return collect( + name, + tind=[1, 1], + info=False, + path=DIRECTORY, + xguards=False, + yguards=False, + )[()] + + +def assert_convergence(error, dx, name, expected) -> bool: + fit = polyfit(log(dx), log(error), 1) + order = fit[0] + print(f"{name} convergence order = {order:f} (fit)", end="") + + order = log(error[-2] / error[-1]) / log(dx[-2] / dx[-1]) + print(f", {order:f} (small spacing)", end="") + + # Should be close to the expected order + success = order > expected * 0.95 + print(f"\t............ {'PASS' if success else 'FAIL'}") + + return success + + +def run_fv_operators(nz: int) -> dict[str, float]: + # Command to run + cmd = f"./fv_mms MZ={nz} mesh:ny={nz}" + print(f"Running command: {cmd}", end="") + + # Launch using MPI + start = time() + status, out = launch_safe(cmd, nproc=NPROC, mthread=MTHREAD, pipe=True) + print(f" ... done in {time() - start:.3}s") + + # Save output to log file + pathlib.Path(f"run.log.{nz}").write_text(out) + + if status: + print(f"Run failed!\nOutput was:\n{out}") + sys.exit(status) + + return { + operator: { + "l_2": quiet_collect(f"{operator}_l_2"), + "l_inf": quiet_collect(f"{operator}_l_inf"), + } + for operator in OPERATORS + } + + +def transpose( + errors: list[dict[str, dict[str, float]]], +) -> dict[str, dict[str, list[float]]]: + """Turn a list of {operator: error} into a dict of {operator: [errors]}""" + + kinds = ("l_2", "l_inf") + result = {operator: {kind: [] for kind in kinds} for operator in OPERATORS} + for error in errors: + for k, v in error.items(): + for kind in kinds: + result[k][kind].append(v[kind]) + return result + + +def test_fv_operators() -> bool: + failures = [] + + all_errors = [] + + for n in NLIST: + errors = run_fv_operators(n) + all_errors.append(errors) + + for operator in OPERATORS: + l_2 = errors[operator]["l_2"] + l_inf = errors[operator]["l_inf"] + + print(f"{operator} errors: l-2 {l_2:f} l-inf {l_inf:f}") + + final_errors = transpose(all_errors) + for operator, order in OPERATORS.items(): + success = assert_convergence(final_errors[operator]["l_2"], dx, operator, order) + if not success: + failures.append(operator) + + return final_errors, failures + + +def make_plots(cases: dict[str, dict]): + try: + import matplotlib.pyplot as plt + except ImportError: + print("No matplotlib") + return + + num_operators = len(OPERATORS) + fig, axes = plt.subplots(1, num_operators, figsize=(num_operators * 4, 4)) + + for ax, operator in zip(axes, OPERATORS): + ax.loglog(dx, cases[operator]["l_2"], "-", label="$l_2$") + ax.loglog(dx, cases[operator]["l_inf"], "--", label="$l_\\inf$") + ax.legend(loc="upper left") + ax.grid() + ax.set_title(f"Error scaling for {operator}") + ax.set_xlabel(r"Mesh spacing $\delta x$") + ax.set_ylabel("Error norm") + + fig.tight_layout() + fig.savefig("fv_mms.pdf") + print("Plot saved to fv_mms.pdf") + + if args.show_plots: + plt.show() + plt.close() + + +if __name__ == "__main__": + build_and_log("Finite volume MMS test") + + parser = argparse.ArgumentParser("Error scaling test for finite volume operators") + parser.add_argument( + "--make-plots", action="store_true", help="Create plots of error scaling" + ) + parser.add_argument( + "--show-plots", + action="store_true", + help="Stop and show plots, implies --make-plots", + ) + parser.add_argument( + "--dump-errors", + type=str, + help="Output file to dump errors as JSON", + default="fv_operator_errors.json", + ) + + args = parser.parse_args() + + error2, failures = test_fv_operators() + success = len(failures) == 0 + + if args.dump_errors: + pathlib.Path(args.dump_errors).write_text(json.dumps(error2)) + + if args.make_plots or args.show_plots: + make_plots(error2) + + if success: + print("\nAll tests passed") + else: + print("\nSome tests failed:") + for failure in failures: + print(f"\t{failure}") + + sys.exit(0 if success else 1) diff --git a/tests/integrated/test-fci-mpi/CMakeLists.txt b/tests/integrated/test-fci-mpi/CMakeLists.txt index 0abbbce099..ba126d7e76 100644 --- a/tests/integrated/test-fci-mpi/CMakeLists.txt +++ b/tests/integrated/test-fci-mpi/CMakeLists.txt @@ -4,7 +4,7 @@ bout_add_mms_test( USE_RUNTEST USE_DATA_BOUT_INP PROCESSORS 6 DOWNLOAD - https://zenodo.org/records/7614499/files/W7X-conf4-36x8x128.fci.nc?download=1 + https://zenodo.org/records/17571607/files/W7X-conf0-36x8x128.fci.nc?download=1 DOWNLOAD_NAME grid.fci.nc REQUIRES BOUT_HAS_PETSC ) diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx index d80affc3f7..7672ec7dae 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/tests/unit/field/test_field3d.cxx @@ -2450,5 +2450,25 @@ TEST_F(Field3DTest, ZeroFrom) { EXPECT_TRUE(field2.isAllocated()); EXPECT_TRUE(IsFieldEqual(field2, 0.)); } + +TEST_F(Field3DTest, Field3DParallel) { + Field3DParallel field(1.0); + field = 1.0; + + Field3D field2 = field; + + auto& field3 = field.asField3D(); + + field *= 2; + + EXPECT_TRUE(IsFieldEqual(field, 2.0)); + EXPECT_TRUE(IsFieldEqual(field2, 1.0)); + EXPECT_TRUE(IsFieldEqual(field3, 2.0)); + + field3.asField3DParallel() *= 3; + + EXPECT_TRUE(IsFieldEqual(field3, 6.0)); +} + // Restore compiler warnings #pragma GCC diagnostic pop