diff --git a/CMakeLists.txt b/CMakeLists.txt index 9e57495885..73a40ccd66 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -593,9 +593,6 @@ else() endif() set(BOUT_USE_METRIC_3D ${BOUT_ENABLE_METRIC_3D}) -option(BOUT_ENABLE_FCI_AUTOMAGIC "Enable (slow?) automatic features for FCI" ON) -set(BOUT_USE_FCI_AUTOMAGIC ${BOUT_ENABLE_FCI_AUTOMAGIC}) - include(CheckCXXSourceCompiles) check_cxx_source_compiles("int main() { const char* name = __PRETTY_FUNCTION__; }" HAS_PRETTY_FUNCTION) diff --git a/include/bout/field.hxx b/include/bout/field.hxx index 27835ecbd7..81130fcad5 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -527,22 +527,36 @@ T pow(BoutReal lhs, const T& rhs, const std::string& rgn = "RGN_ALL") { * and uses checkData() to, if CHECK >= 3, check * result for non-finite numbers * + * If the input field has parallel slices, then those will also be + * operated on. To avoid this, use `withoutParallelSlices()` to + * discard slices before calling. + * */ #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") { \ - AUTO_TRACE(); \ - /* 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]); } \ - result.name = std::string(#_name "(") + f.name + std::string(")"); \ - checkData(result); \ - return result; \ +#define FIELD_FUNC(_name, func) \ + template > \ + inline T _name(const T& f, const std::string& rgn = "RGN_ALL") { \ + AUTO_TRACE(); \ + /* 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 (f.hasParallelSlices()) { \ + /* Operate on parallel slices */ \ + result.splitParallelSlicesAndAllocate(); \ + for (size_t i{0}; i != f.numberParallelSlices(); ++i) { \ + BOUT_FOR(d, result.getRegion(rgn)) { result.yup(i)[d] = func(f.yup(i)[d]); } \ + checkData(result.yup(i)); \ + BOUT_FOR(d, result.getRegion(rgn)) { result.ydown(i)[d] = func(f.ydown(i)[d]); } \ + checkData(result.ydown(i)); \ + } \ + } \ + result.name = std::string(#_name "(") + f.name + std::string(")"); \ + checkData(result); \ + return result; \ } #endif @@ -670,6 +684,8 @@ T copy(const T& f) { /// Apply a floor value \p f to a field \p var. Any value lower than /// the floor is set to the floor. /// +/// If the field has parallel slices then they will also be floored. +/// /// @param[in] var Variable to apply floor to /// @param[in] f The floor value /// @param[in] rgn The region to calculate the result over @@ -683,8 +699,8 @@ inline T floor(const T& var, BoutReal f, const std::string& rgn = "RGN_ALL") { result[d] = f; } } -#if BOUT_USE_FCI_AUTOMAGIC - if (var.isFci()) { + + if (var.hasParallelSlices()) { for (size_t i = 0; i < result.numberParallelSlices(); ++i) { BOUT_FOR(d, result.yup(i).getRegion(rgn)) { if (result.yup(i)[d] < f) { @@ -697,10 +713,6 @@ inline T floor(const T& var, BoutReal f, const std::string& rgn = "RGN_ALL") { } } } - } else -#endif - { - result.clearParallelSlices(); } return result; } diff --git a/include/bout/field2d.hxx b/include/bout/field2d.hxx index 5eab330e8e..eaf5e21d72 100644 --- a/include/bout/field2d.hxx +++ b/include/bout/field2d.hxx @@ -134,12 +134,14 @@ public: return *this; } - /// Dummy functions to increase portability + /// Dummy functions to replicate Field3D interface bool hasParallelSlices() const { return true; } void calcParallelSlices() const {} void splitParallelSlices() const {} + void splitParallelSlicesAndAllocate() const {} void clearParallelSlices() const {} - int numberParallelSlices() const { return 0; } + Field2D withoutParallelSlices() const { return *this; } + size_t numberParallelSlices() const { return 0; } Field2D& yup(std::vector::size_type UNUSED(index) = 0) { return *this; } const Field2D& yup(std::vector::size_type UNUSED(index) = 0) const { diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index 0643efbe0a..c162c12237 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -237,6 +237,9 @@ public: */ void splitParallelSlices(); + /// Separate fields for yup and ydown, with memory allocated. Note: + /// After this the parallel slices will be allocated and unique, but + /// may contain uninitialised values. void splitParallelSlicesAndAllocate(); /*! @@ -244,6 +247,13 @@ public: */ void clearParallelSlices(); + /// Returns a shallow copy without parallel slices + Field3D withoutParallelSlices() const { + Field3D result{getMesh(), getLocation(), getDirections(), getRegionID()}; + result.data = data; + return result; + } + /// Check if this field has yup and ydown fields bool hasParallelSlices() const { #if CHECK > 2 @@ -488,10 +498,11 @@ public: friend class Vector2D; Field3D& calcParallelSlices(); - void allowParallelSlices([[maybe_unused]] bool allow) { + Field3D& allowParallelSlices([[maybe_unused]] bool allow) { #if CHECK > 0 allowCalcParallelSlices = allow; #endif + return *this; } void applyBoundary(bool init = false) override; diff --git a/include/bout/fieldperp.hxx b/include/bout/fieldperp.hxx index b50eef1991..376d333d1c 100644 --- a/include/bout/fieldperp.hxx +++ b/include/bout/fieldperp.hxx @@ -158,11 +158,14 @@ public: return *this; } - /// Dummy functions to increase portability + /// Dummy functions to match Field3D interface bool hasParallelSlices() const { return true; } + void splitParallelSlices() const {} + void splitParallelSlicesAndAllocate() const {} void calcParallelSlices() const {} - void clearParallelSlices() {} - int numberParallelSlices() { return 0; } + void clearParallelSlices() const {} + FieldPerp withoutParallelSlices() const { return *this; } + size_t numberParallelSlices() const { return 0; } FieldPerp& yup(std::vector::size_type UNUSED(index) = 0) { return *this; } const FieldPerp& yup(std::vector::size_type UNUSED(index) = 0) const { diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 8a9baaf3e7..454818d8b5 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -188,16 +188,13 @@ void communicateFluxes(Field3D& f); /// @param[in] fixflux Fix the flux at the boundary to be the value at the /// midpoint (for boundary conditions) /// -/// NB: Uses to/from FieldAligned coordinates +/// Handles FCI by calling ::Div_par(f_in, v_in) template const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, bool fixflux = true) { - -#if BOUT_USE_FCI_AUTOMAGIC if (f_in.isFci()) { return ::Div_par(f_in, v_in); } -#endif ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in); diff --git a/include/bout/index_derivs_interface.hxx b/include/bout/index_derivs_interface.hxx index bc9a687b34..fb1e7ef020 100644 --- a/include/bout/index_derivs_interface.hxx +++ b/include/bout/index_derivs_interface.hxx @@ -4,10 +4,9 @@ * Definition of main derivative kernels * ************************************************************************** - * Copyright 2018 - * D.Dickinson + * Copyright 2018 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -204,13 +203,9 @@ T DDY(const T& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "D ASSERT1(f.getDirectionY() == YDirectionType::Standard); T f_tmp = f; if (!f.hasParallelSlices()) { -#if BOUT_USE_FCI_AUTOMAGIC - f_tmp.calcParallelSlices(); -#else throw BoutException( "parallel slices needed for parallel derivatives. Make sure to communicate and " "apply parallel boundary conditions before calling derivative"); -#endif } return standardDerivative(f_tmp, outloc, method, region); diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 849bb0ffe3..453f7099e4 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -150,7 +150,7 @@ public: return bndry_position != rhs.bndry_position; } -#define ITER() for (int i = 0; i < localmesh->ystart - abs_offset(); ++i) +#define ITER() for (int i = 0; i <= localmesh->ystart - abs_offset(); ++i) // dirichlet boundary code void dirichlet_o1(Field3D& f, BoutReal value) const { ITER() { getAt(f, i) = value; } diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 74a6f0853c..a124141634 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -93,7 +93,7 @@ Field3D::Field3D(const BoutReal val, Mesh* localmesh) : Field3D(localmesh) { TRACE("Field3D: Copy constructor from value"); *this = val; -#if BOUT_USE_FCI_AUTOMAGIC + if (this->isFci()) { splitParallelSlices(); for (size_t i = 0; i < numberParallelSlices(); ++i) { @@ -101,7 +101,6 @@ Field3D::Field3D(const BoutReal val, Mesh* localmesh) : Field3D(localmesh) { ydown(i) = val; } } -#endif } Field3D::Field3D(Array data_in, Mesh* localmesh, CELL_LOC datalocation, @@ -361,18 +360,13 @@ Field3D& Field3D::operator=(const BoutReal val) { TRACE("Field3D = BoutReal"); track(val, "operator="); -#if BOUT_USE_FCI_AUTOMAGIC - if (isFci() && hasParallelSlices()) { + if (hasParallelSlices()) { for (size_t i = 0; i < numberParallelSlices(); ++i) { yup(i) = val; ydown(i) = val; } } -#else - // Delete existing parallel slices. We don't copy parallel slices, so any - // that currently exist will be incorrect. - clearParallelSlices(); -#endif + resetRegion(); allocate(); @@ -386,11 +380,13 @@ Field3D& Field3D::operator=(const BoutReal val) { Field3D& Field3D::calcParallelSlices() { ASSERT2(allowCalcParallelSlices); getCoordinates()->getParallelTransform().calcParallelSlices(*this); -#if BOUT_USE_FCI_AUTOMAGIC - if (this->isFci()) { - this->applyParallelBoundary("parallel_neumann_o2"); - } -#endif + + // This seems like it could work but give wrong results. + // Better to insert NaNs into boundary cells + // if (this->isFci()) { + // this->applyParallelBoundary("parallel_neumann_o2"); + // } + return *this; } @@ -548,7 +544,12 @@ void Field3D::applyParallelBoundary() { TRACE("Field3D::applyParallelBoundary()"); checkData(*this); - ASSERT1(hasParallelSlices()); + if (isFci()) { + ASSERT1(hasParallelSlices()); + } + if (!hasParallelSlices()) { + return; + } // Apply boundary to this field for (const auto& bndry : getBoundaryOpPars()) { @@ -561,7 +562,12 @@ void Field3D::applyParallelBoundary(BoutReal t) { TRACE("Field3D::applyParallelBoundary(t)"); checkData(*this); - ASSERT1(hasParallelSlices()); + if (isFci()) { + ASSERT1(hasParallelSlices()); + } + if (!hasParallelSlices()) { + return; + } // Apply boundary to this field for (const auto& bndry : getBoundaryOpPars()) { @@ -574,7 +580,12 @@ void Field3D::applyParallelBoundary(const std::string& condition) { TRACE("Field3D::applyParallelBoundary(condition)"); checkData(*this); - ASSERT1(hasParallelSlices()); + if (isFci()) { + ASSERT1(hasParallelSlices()); + } + if (!hasParallelSlices()) { + return; + } /// Get the boundary factory (singleton) BoundaryFactory* bfact = BoundaryFactory::getInstance(); @@ -593,7 +604,12 @@ void Field3D::applyParallelBoundary(const std::string& region, TRACE("Field3D::applyParallelBoundary(region, condition)"); checkData(*this); - ASSERT1(hasParallelSlices()); + if (isFci()) { + ASSERT1(hasParallelSlices()); + } + if (!hasParallelSlices()) { + return; + } /// Get the boundary factory (singleton) BoundaryFactory* bfact = BoundaryFactory::getInstance(); @@ -615,7 +631,12 @@ void Field3D::applyParallelBoundary(const std::string& region, TRACE("Field3D::applyParallelBoundary(region, condition, f)"); checkData(*this); - ASSERT1(hasParallelSlices()); + if (isFci()) { + ASSERT1(hasParallelSlices()); + } + if (!hasParallelSlices()) { + return; + } /// Get the boundary factory (singleton) BoundaryFactory* bfact = BoundaryFactory::getInstance(); diff --git a/src/field/gen_fieldops.jinja b/src/field/gen_fieldops.jinja index 88e877c197..9d1c3ee68d 100644 --- a/src/field/gen_fieldops.jinja +++ b/src/field/gen_fieldops.jinja @@ -12,40 +12,34 @@ {% if lhs == rhs == "Field3D" %} {{out.name}}.setRegion({{lhs.name}}.getMesh()->getCommonRegion({{lhs.name}}.getRegionID(), {{rhs.name}}.getRegionID())); -#if BOUT_USE_FCI_AUTOMAGIC - if ({{lhs.name}}.isFci() and {{lhs.name}}.hasParallelSlices() and {{rhs.name}}.hasParallelSlices()) { + if ({{lhs.name}}.hasParallelSlices() and {{rhs.name}}.hasParallelSlices()) { {{out.name}}.splitParallelSlices(); for (size_t i{0} ; i < {{lhs.name}}.numberParallelSlices() ; ++i) { {{out.name}}.yup(i) = {{lhs.name}}.yup(i) {{operator}} {{rhs.name}}.yup(i); {{out.name}}.ydown(i) = {{lhs.name}}.ydown(i) {{operator}} {{rhs.name}}.ydown(i); } } -#endif {% elif lhs == "Field3D" %} {{out.name}}.setRegion({{lhs.name}}.getRegionID()); {% if rhs == "BoutReal" %} -#if BOUT_USE_FCI_AUTOMAGIC - if ({{lhs.name}}.isFci() and {{lhs.name}}.hasParallelSlices()) { + if ({{lhs.name}}.hasParallelSlices()) { {{out.name}}.splitParallelSlices(); for (size_t i{0} ; i < {{lhs.name}}.numberParallelSlices() ; ++i) { {{out.name}}.yup(i) = {{lhs.name}}.yup(i) {{operator}} {{rhs.name}}; {{out.name}}.ydown(i) = {{lhs.name}}.ydown(i) {{operator}} {{rhs.name}}; } } -#endif {% endif %} {% elif rhs == "Field3D" %} {{out.name}}.setRegion({{rhs.name}}.getRegionID()); {% if lhs == "BoutReal" %} -#if BOUT_USE_FCI_AUTOMAGIC - if ({{rhs.name}}.isFci() and {{rhs.name}}.hasParallelSlices()) { + if ({{rhs.name}}.hasParallelSlices()) { {{out.name}}.splitParallelSlices(); for (size_t i{0} ; i < {{rhs.name}}.numberParallelSlices() ; ++i) { {{out.name}}.yup(i) = {{lhs.name}} {{operator}} {{rhs.name}}.yup(i); {{out.name}}.ydown(i) = {{lhs.name}} {{operator}} {{rhs.name}}.ydown(i); } } -#endif {% endif %} {% endif %} {% endif %} @@ -114,14 +108,12 @@ // Delete existing parallel slices. We don't copy parallel slices, so any // that currently exist will be incorrect. {% if (rhs == "Field3D" or rhs == "BoutReal") %} -#if BOUT_USE_FCI_AUTOMAGIC - if (this->isFci() and this->hasParallelSlices() {% if rhs == "Field3D" %} and {{rhs.name}}.hasParallelSlices() {% endif %}) { + if (this->hasParallelSlices() {% if rhs == "Field3D" %} and {{rhs.name}}.hasParallelSlices() {% endif %}) { 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 -#endif {% endif %} { clearParallelSlices(); diff --git a/src/field/generated_fieldops.cxx b/src/field/generated_fieldops.cxx index 75d2ede82d..1664668e5f 100644 --- a/src/field/generated_fieldops.cxx +++ b/src/field/generated_fieldops.cxx @@ -15,15 +15,13 @@ Field3D operator*(const Field3D& lhs, const Field3D& rhs) { checkData(rhs); result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); -#if BOUT_USE_FCI_AUTOMAGIC - if (lhs.isFci() and lhs.hasParallelSlices() and rhs.hasParallelSlices()) { + if (lhs.hasParallelSlices() and rhs.hasParallelSlices()) { 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); } } -#endif BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { result[index] = lhs[index] * rhs[index]; @@ -44,16 +42,13 @@ Field3D& Field3D::operator*=(const Field3D& rhs) { ASSERT1_FIELDS_COMPATIBLE(*this, rhs); // Delete existing parallel slices. We don't copy parallel slices, so any -// that currently exist will be incorrect. -#if BOUT_USE_FCI_AUTOMAGIC - if (this->isFci() and this->hasParallelSlices() and rhs.hasParallelSlices()) { + // that currently exist will be incorrect. + if (this->hasParallelSlices() and rhs.hasParallelSlices()) { for (size_t i{0}; i < yup_fields.size(); ++i) { yup(i) *= rhs.yup(i); ydown(i) *= rhs.ydown(i); } - } else -#endif - { + } else { clearParallelSlices(); } @@ -87,15 +82,13 @@ Field3D operator/(const Field3D& lhs, const Field3D& rhs) { checkData(rhs); result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); -#if BOUT_USE_FCI_AUTOMAGIC - if (lhs.isFci() and lhs.hasParallelSlices() and rhs.hasParallelSlices()) { + if (lhs.hasParallelSlices() and rhs.hasParallelSlices()) { 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); } } -#endif BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { result[index] = lhs[index] / rhs[index]; @@ -116,16 +109,13 @@ Field3D& Field3D::operator/=(const Field3D& rhs) { ASSERT1_FIELDS_COMPATIBLE(*this, rhs); // Delete existing parallel slices. We don't copy parallel slices, so any -// that currently exist will be incorrect. -#if BOUT_USE_FCI_AUTOMAGIC - if (this->isFci() and this->hasParallelSlices() and rhs.hasParallelSlices()) { + // that currently exist will be incorrect. + if (this->hasParallelSlices() and rhs.hasParallelSlices()) { for (size_t i{0}; i < yup_fields.size(); ++i) { yup(i) /= rhs.yup(i); ydown(i) /= rhs.ydown(i); } - } else -#endif - { + } else { clearParallelSlices(); } @@ -159,15 +149,13 @@ Field3D operator+(const Field3D& lhs, const Field3D& rhs) { checkData(rhs); result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); -#if BOUT_USE_FCI_AUTOMAGIC - if (lhs.isFci() and lhs.hasParallelSlices() and rhs.hasParallelSlices()) { + if (lhs.hasParallelSlices() and rhs.hasParallelSlices()) { 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); } } -#endif BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { result[index] = lhs[index] + rhs[index]; @@ -188,16 +176,13 @@ Field3D& Field3D::operator+=(const Field3D& rhs) { ASSERT1_FIELDS_COMPATIBLE(*this, rhs); // Delete existing parallel slices. We don't copy parallel slices, so any -// that currently exist will be incorrect. -#if BOUT_USE_FCI_AUTOMAGIC - if (this->isFci() and this->hasParallelSlices() and rhs.hasParallelSlices()) { + // that currently exist will be incorrect. + if (this->hasParallelSlices() and rhs.hasParallelSlices()) { for (size_t i{0}; i < yup_fields.size(); ++i) { yup(i) += rhs.yup(i); ydown(i) += rhs.ydown(i); } - } else -#endif - { + } else { clearParallelSlices(); } @@ -231,15 +216,13 @@ Field3D operator-(const Field3D& lhs, const Field3D& rhs) { checkData(rhs); result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); -#if BOUT_USE_FCI_AUTOMAGIC - if (lhs.isFci() and lhs.hasParallelSlices() and rhs.hasParallelSlices()) { + if (lhs.hasParallelSlices() and rhs.hasParallelSlices()) { 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); } } -#endif BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { result[index] = lhs[index] - rhs[index]; @@ -260,16 +243,13 @@ Field3D& Field3D::operator-=(const Field3D& rhs) { ASSERT1_FIELDS_COMPATIBLE(*this, rhs); // Delete existing parallel slices. We don't copy parallel slices, so any -// that currently exist will be incorrect. -#if BOUT_USE_FCI_AUTOMAGIC - if (this->isFci() and this->hasParallelSlices() and rhs.hasParallelSlices()) { + // that currently exist will be incorrect. + if (this->hasParallelSlices() and rhs.hasParallelSlices()) { for (size_t i{0}; i < yup_fields.size(); ++i) { yup(i) -= rhs.yup(i); ydown(i) -= rhs.ydown(i); } - } else -#endif - { + } else { clearParallelSlices(); } @@ -640,15 +620,13 @@ Field3D operator*(const Field3D& lhs, const BoutReal rhs) { checkData(rhs); result.setRegion(lhs.getRegionID()); -#if BOUT_USE_FCI_AUTOMAGIC - if (lhs.isFci() and lhs.hasParallelSlices()) { + if (lhs.hasParallelSlices()) { 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; } } -#endif BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { result[index] = lhs[index] * rhs; @@ -668,16 +646,13 @@ Field3D& Field3D::operator*=(const BoutReal rhs) { if (data.unique()) { // Delete existing parallel slices. We don't copy parallel slices, so any -// that currently exist will be incorrect. -#if BOUT_USE_FCI_AUTOMAGIC - if (this->isFci() and this->hasParallelSlices()) { + // that currently exist will be incorrect. + if (this->hasParallelSlices()) { for (size_t i{0}; i < yup_fields.size(); ++i) { yup(i) *= rhs; ydown(i) *= rhs; } - } else -#endif - { + } else { clearParallelSlices(); } @@ -708,15 +683,13 @@ Field3D operator/(const Field3D& lhs, const BoutReal rhs) { checkData(rhs); result.setRegion(lhs.getRegionID()); -#if BOUT_USE_FCI_AUTOMAGIC - if (lhs.isFci() and lhs.hasParallelSlices()) { + if (lhs.hasParallelSlices()) { 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; } } -#endif const auto tmp = 1.0 / rhs; BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { @@ -737,16 +710,13 @@ Field3D& Field3D::operator/=(const BoutReal rhs) { if (data.unique()) { // Delete existing parallel slices. We don't copy parallel slices, so any -// that currently exist will be incorrect. -#if BOUT_USE_FCI_AUTOMAGIC - if (this->isFci() and this->hasParallelSlices()) { + // that currently exist will be incorrect. + if (this->hasParallelSlices()) { for (size_t i{0}; i < yup_fields.size(); ++i) { yup(i) /= rhs; ydown(i) /= rhs; } - } else -#endif - { + } else { clearParallelSlices(); } @@ -778,15 +748,13 @@ Field3D operator+(const Field3D& lhs, const BoutReal rhs) { checkData(rhs); result.setRegion(lhs.getRegionID()); -#if BOUT_USE_FCI_AUTOMAGIC - if (lhs.isFci() and lhs.hasParallelSlices()) { + if (lhs.hasParallelSlices()) { 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; } } -#endif BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { result[index] = lhs[index] + rhs; @@ -806,16 +774,13 @@ Field3D& Field3D::operator+=(const BoutReal rhs) { if (data.unique()) { // Delete existing parallel slices. We don't copy parallel slices, so any -// that currently exist will be incorrect. -#if BOUT_USE_FCI_AUTOMAGIC - if (this->isFci() and this->hasParallelSlices()) { + // that currently exist will be incorrect. + if (this->hasParallelSlices()) { for (size_t i{0}; i < yup_fields.size(); ++i) { yup(i) += rhs; ydown(i) += rhs; } - } else -#endif - { + } else { clearParallelSlices(); } @@ -846,15 +811,13 @@ Field3D operator-(const Field3D& lhs, const BoutReal rhs) { checkData(rhs); result.setRegion(lhs.getRegionID()); -#if BOUT_USE_FCI_AUTOMAGIC - if (lhs.isFci() and lhs.hasParallelSlices()) { + if (lhs.hasParallelSlices()) { 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; } } -#endif BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { result[index] = lhs[index] - rhs; @@ -874,16 +837,13 @@ Field3D& Field3D::operator-=(const BoutReal rhs) { if (data.unique()) { // Delete existing parallel slices. We don't copy parallel slices, so any -// that currently exist will be incorrect. -#if BOUT_USE_FCI_AUTOMAGIC - if (this->isFci() and this->hasParallelSlices()) { + // that currently exist will be incorrect. + if (this->hasParallelSlices()) { for (size_t i{0}; i < yup_fields.size(); ++i) { yup(i) -= rhs; ydown(i) -= rhs; } - } else -#endif - { + } else { clearParallelSlices(); } @@ -2209,15 +2169,13 @@ Field3D operator*(const BoutReal lhs, const Field3D& rhs) { checkData(rhs); result.setRegion(rhs.getRegionID()); -#if BOUT_USE_FCI_AUTOMAGIC - if (rhs.isFci() and rhs.hasParallelSlices()) { + if (rhs.hasParallelSlices()) { 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); } } -#endif BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { result[index] = lhs * rhs[index]; @@ -2238,15 +2196,13 @@ Field3D operator/(const BoutReal lhs, const Field3D& rhs) { checkData(rhs); result.setRegion(rhs.getRegionID()); -#if BOUT_USE_FCI_AUTOMAGIC - if (rhs.isFci() and rhs.hasParallelSlices()) { + if (rhs.hasParallelSlices()) { 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); } } -#endif BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { result[index] = lhs / rhs[index]; @@ -2267,15 +2223,13 @@ Field3D operator+(const BoutReal lhs, const Field3D& rhs) { checkData(rhs); result.setRegion(rhs.getRegionID()); -#if BOUT_USE_FCI_AUTOMAGIC - if (rhs.isFci() and rhs.hasParallelSlices()) { + if (rhs.hasParallelSlices()) { 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); } } -#endif BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { result[index] = lhs + rhs[index]; @@ -2296,15 +2250,13 @@ Field3D operator-(const BoutReal lhs, const Field3D& rhs) { checkData(rhs); result.setRegion(rhs.getRegionID()); -#if BOUT_USE_FCI_AUTOMAGIC - if (rhs.isFci() and rhs.hasParallelSlices()) { + if (rhs.hasParallelSlices()) { 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); } } -#endif BOUT_FOR(index, result.getValidRegionWithDefault("RGN_ALL")) { result[index] = lhs - rhs[index]; diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index d013634644..bcbd141890 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1600,15 +1600,15 @@ Field3D Coordinates::Div_par(const Field3D& f, CELL_LOC outloc, return Bxy * Grad_par(f / Bxy_floc, outloc, method); } - auto coords = f.getCoordinates(); // Need to modify yup and ydown fields - Field3D f_B = f / coords->J * sqrt(coords->g_22); - f_B.splitParallelSlices(); - for (int i = 0; i < f.getMesh()->ystart; ++i) { - f_B.yup(i) = f.yup(i) / coords->J.yup(i) * sqrt(coords->g_22.yup(i)); - f_B.ydown(i) = f.ydown(i) / coords->J.ydown(i) * sqrt(coords->g_22.ydown(i)); - } - return setName(coords->J / sqrt(coords->g_22) * Grad_par(f_B, outloc, method), + //Field3D f_B = f / Bxy_floc; + //return setName(Bxy * Grad_par(f_B, outloc, method), + // "Div_par({:s})", f.name); + + auto coords = f.getCoordinates(); + // Note: Arithmetic operators iterate over yup/ydown slices + Field3D f_B = f * coords->J / sqrt(coords->g_22); + return setName(sqrt(coords->g_22) * Grad_par(f_B, outloc, method) / coords->J, "Div_par({:s})", f.name); } diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 639a5a1612..1537650b0e 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -31,6 +31,7 @@ #include #include #include +#include #include #include #include @@ -234,57 +235,54 @@ Field3D Div_par(const Field3D& f, const std::string& method, CELL_LOC outloc) { return f.getCoordinates(outloc)->Div_par(f, outloc, method); } -Field3D Div_par(const Field3D& f_in, const Field3D& v_in) { -#if BOUT_USE_FCI_AUTOMAGIC - auto f{f_in}; - auto v{v_in}; +Field3D Div_par(const Field3D& f, const Field3D& v) { + AUTO_TRACE(); + ASSERT1_FIELDS_COMPATIBLE(f, v); + + // Either both have parallel slices or neither if (!f.hasParallelSlices()) { - f.calcParallelSlices(); - } - if (!v.hasParallelSlices()) { - v.calcParallelSlices(); + // No parallel slices + ASSERT1(!v.hasParallelSlices()); + + return Div_par(f * v); } -#else - const auto& f{f_in}; - const auto& v{v_in}; -#endif - ASSERT1_FIELDS_COMPATIBLE(f, v); + // Using parallel slices ASSERT1(f.hasParallelSlices()); ASSERT1(v.hasParallelSlices()); - // Parallel divergence, using velocities at cell boundaries - // Note: Not guaranteed to be flux conservative - Mesh* mesh = f.getMesh(); - - Field3D result{emptyFrom(f)}; - Coordinates* coord = f.getCoordinates(); - for (int i = mesh->xstart; i <= mesh->xend; i++) { - for (int j = mesh->ystart; j <= mesh->yend; j++) { - for (int k = mesh->zstart; k <= mesh->zend; k++) { - // Value of f and v at left cell face - BoutReal fL = 0.5 * (f(i, j, k) + f.ydown()(i, j - 1, k)); - BoutReal vL = 0.5 * (v(i, j, k) + v.ydown()(i, j - 1, k)); + auto B = coord->Bxy; + auto B_up = coord->Bxy.yup(); + auto B_down = coord->Bxy.ydown(); - BoutReal fR = 0.5 * (f(i, j, k) + f.yup()(i, j + 1, k)); - BoutReal vR = 0.5 * (v(i, j, k) + v.yup()(i, j + 1, k)); + auto f_up = f.yup(); + auto f_down = f.ydown(); - // Calculate flux at right boundary (y+1/2) - BoutReal fluxRight = - fR * vR * (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))); + auto v_up = v.yup(); + auto v_down = v.ydown(); - // Calculate at left boundary (y-1/2) - BoutReal fluxLeft = - fL * vL * (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))); + auto g_22 = coord->g_22; + auto dy = coord->dy; - result(i, j, k) = - (fluxRight - fluxLeft) / (coord->dy(i, j, k) * coord->J(i, j, k)); - } + Field3D result{emptyFrom(f)}; + BOUT_FOR(i, f.getRegion("RGN_NOBNDRY")) { + result[i] = B[i] + * ((f_up[i] * v_up[i] / B_up[i]) - (f_down[i] * v_down[i] / B_down[i])) + / (2 * dy[i] * sqrt(g_22[i])); + +#if CHECK > 2 + if (!std::isfinite(result[i])) { + throw BoutException("Non-finite value in Div_par(f, v) at {}\n" + " f down {} up {}\n" + " v down {} up {}\n" + " B down {} central {} up {}\n" + " dy {} sqrt(g_22) {}\n", + i, f_down[i], f_up[i], v_down[i], v_up[i], B_down[i], B[i], + B_up[i], dy[i], sqrt(g_22[i])); } +#endif } return result; diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index 580897f47a..eeccd0e7b1 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -108,6 +108,13 @@ void load_parallel_metric_components([[maybe_unused]] Coordinates* coords, #if BOUT_USE_METRIC_3D #define LOAD_PAR(var, doZero) \ load_parallel_metric_component(#var, coords->var, offset, doZero) + + // Parallel slices of metric components must NOT be calculated by + // interpolation. The X and Z coordinates are defined independently on + // on each poloidal plane. Instead, these yup/ydown fields are calculated + // by mapping coordinate points, and calculating metrics on the mapped points. + LOAD_PAR(dy, false); + LOAD_PAR(g11, false); LOAD_PAR(g22, false); LOAD_PAR(g33, false); @@ -115,6 +122,8 @@ void load_parallel_metric_components([[maybe_unused]] Coordinates* coords, LOAD_PAR(g13, false); LOAD_PAR(g23, false); + // LOAD_PAR(Bxy, false); // Not yet written to mesh file + LOAD_PAR(g_11, false); LOAD_PAR(g_22, false); LOAD_PAR(g_33, false); @@ -430,6 +439,7 @@ void FCITransform::calcParallelSlices(Field3D& f) { // Interpolate f onto yup and ydown fields for (const auto& map : field_line_maps) { f.ynext(map.offset) = map.interpolate(f); + //f.ynext(map.offset) = map.integrate(f); f.ynext(map.offset).setRegion(fmt::format("RGN_YPAR_{:+d}", map.offset)); } }