From d09d8f8ad023532d81a775ac2885e8c47bd5338e Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 15 Mar 2025 15:10:07 -0700 Subject: [PATCH 01/18] difops: Replace Div_par(Field3D, Field3D) When using FCI yup/down fields, each poloidal plane uses a different coordinate system. Quantities like J therefore can't be averaged between planes. Magnetic field strength is a scalar that can be interpolated, and is used here to calculate the divergence of a parallel flow. --- src/mesh/difops.cxx | 69 ++++++++++++++++++++------------------------- 1 file changed, 31 insertions(+), 38 deletions(-) diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 639a5a1612..5d9f579515 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -234,57 +234,50 @@ 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])) + / (dy[i] * sqrt(g_22[i])); + +#if CHECK > 0 + if(!std::isfinite(result[i])) { + output.write("{} {} {} {}\n", f_up[i], v_up[i], f_down[i], v_down[i]); + output.write("{} {} {} {} {}\n", B[i], B_up[i], B_down[i], dy[i], sqrt(g_22[i])); + throw BoutException("Non-finite value in Div_"); } +#endif } return result; From 81ee64eae6298526927bf8ccd02a68c38fb01104 Mon Sep 17 00:00:00 2001 From: bendudson <219233+bendudson@users.noreply.github.com> Date: Sat, 15 Mar 2025 22:16:55 +0000 Subject: [PATCH 02/18] Apply clang-format changes --- src/mesh/difops.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 5d9f579515..f722b85afb 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -267,12 +267,12 @@ Field3D Div_par(const Field3D& f, const Field3D& v) { 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])) - / (dy[i] * sqrt(g_22[i])); + result[i] = B[i] + * ((f_up[i] * v_up[i] / B_up[i]) - (f_down[i] * v_down[i] / B_down[i])) + / (dy[i] * sqrt(g_22[i])); #if CHECK > 0 - if(!std::isfinite(result[i])) { + if (!std::isfinite(result[i])) { output.write("{} {} {} {}\n", f_up[i], v_up[i], f_down[i], v_down[i]); output.write("{} {} {} {} {}\n", B[i], B_up[i], B_down[i], dy[i], sqrt(g_22[i])); throw BoutException("Non-finite value in Div_"); From 3a7743758c7dedca760530601254cc128f640829 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 15 Mar 2025 15:21:57 -0700 Subject: [PATCH 03/18] Coordinates: Fix Div_par(f) for FCI Use Bxy rather than J on neighboring slices: B can be compared between slices, but J cannot (different coordinate system). --- src/mesh/coordinates.cxx | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index d013634644..f39c6c9133 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1600,15 +1600,14 @@ 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); + 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) / 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)); + f_B.yup(i) = f.yup(i) / Bxy_floc.yup(i); + f_B.ydown(i) = f.ydown(i) / Bxy_floc.ydown(i); } - return setName(coords->J / sqrt(coords->g_22) * Grad_par(f_B, outloc, method), + return setName(Bxy * Grad_par(f_B, outloc, method), "Div_par({:s})", f.name); } From 1d785595610ec8c5a2ce2a804c02d2e93b93c00f Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 15 Mar 2025 15:25:42 -0700 Subject: [PATCH 04/18] FCItransform: Only load some parallel metrics Each poloidal slice has a different coordinate system, so metrics can't be directly compared or averaged. Only Bxy and the parallel connection length (dy, g_22) make sense to calculate. --- src/mesh/parallel/fci.cxx | 54 ++++++++++++++++----------------------- 1 file changed, 22 insertions(+), 32 deletions(-) diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index 580897f47a..f7ba0b9aed 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -108,39 +108,29 @@ 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) - LOAD_PAR(g11, false); - LOAD_PAR(g22, false); - LOAD_PAR(g33, false); - LOAD_PAR(g12, false); - LOAD_PAR(g13, false); - LOAD_PAR(g23, false); - - LOAD_PAR(g_11, false); + + // Only some components can be compared between planes + // because each poloidal X-Z plane has a separate coordinate system + LOAD_PAR(dy, false); LOAD_PAR(g_22, false); - LOAD_PAR(g_33, false); - LOAD_PAR(g_12, false); - LOAD_PAR(g_13, false); - LOAD_PAR(g_23, false); - - if (not LOAD_PAR(J, true)) { - auto g = - coords->g11.ynext(offset) * coords->g22.ynext(offset) * coords->g33.ynext(offset) - + 2.0 * coords->g12.ynext(offset) * coords->g13.ynext(offset) - * coords->g23.ynext(offset) - - coords->g11.ynext(offset) * coords->g23.ynext(offset) - * coords->g23.ynext(offset) - - coords->g22.ynext(offset) * coords->g13.ynext(offset) - * coords->g13.ynext(offset) - - coords->g33.ynext(offset) * coords->g12.ynext(offset) - * coords->g12.ynext(offset); - - const auto rgn = fmt::format("RGN_YPAR_{:+d}", offset); - // Check that g is positive - bout::checkPositive(g, "The determinant of g^ij", rgn); - auto J = 1. / sqrt(g); - auto& pcom = coords->J.ynext(offset); - BOUT_FOR(i, J.getRegion(rgn)) { pcom[i] = J[i]; } - } + // LOAD_PAR(Bxy, false); + + // Other components can't have parallel slices + coords->g11.allowParallelSlices(false); + coords->g22.allowParallelSlices(false); + coords->g33.allowParallelSlices(false); + coords->g12.allowParallelSlices(false); + coords->g13.allowParallelSlices(false); + coords->g23.allowParallelSlices(false); + + coords->g_11.allowParallelSlices(false); + coords->g_33.allowParallelSlices(false); + coords->g_12.allowParallelSlices(false); + coords->g_13.allowParallelSlices(false); + coords->g_23.allowParallelSlices(false); + + coords->J.allowParallelSlices(false); + #undef LOAD_PAR #endif } From e5834c18960fa5bcf46201e342793933b3894547 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 15 Mar 2025 15:44:23 -0700 Subject: [PATCH 05/18] DDY: Remove FCI_AUTOMAGIC calcParallelSlices Don't calculate parallel slices in a derivative: The result is unlikely to be correct. --- include/bout/index_derivs_interface.hxx | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) 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); From e058507c293ace2e1ae5c6d1392276369101471d Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 15 Mar 2025 16:14:37 -0700 Subject: [PATCH 06/18] Field: Extend operations to calculate in parallel slices If the argument to `sqrt` has parallel slices, then the slices will be operated on and the result will have parallel slices. To avoid unnecessary work, discard slices before calling. --- include/bout/field.hxx | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/include/bout/field.hxx b/include/bout/field.hxx index 27835ecbd7..ad4fa51b93 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -527,6 +527,10 @@ 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 @@ -540,6 +544,20 @@ T pow(BoutReal lhs, const T& rhs, const std::string& rgn = "RGN_ALL") { /* 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; \ From a46a4f7b0332d70d2382618eff0b6b55fdfd9a05 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 15 Mar 2025 16:36:38 -0700 Subject: [PATCH 07/18] Field: Add withoutParallelSlices() method Returns a shallow copy without parallel slices. Enables user to avoid performing calculations on slices if not needed. --- include/bout/field2d.hxx | 6 ++++-- include/bout/field3d.hxx | 14 +++++++++++++- include/bout/fieldperp.hxx | 9 ++++++--- 3 files changed, 23 insertions(+), 6 deletions(-) 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..96e0722700 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,14 @@ 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 +499,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 { From cedeb2b945aa8fd3eb89b6e76928f9769f407987 Mon Sep 17 00:00:00 2001 From: bendudson <219233+bendudson@users.noreply.github.com> Date: Sat, 15 Mar 2025 23:41:24 +0000 Subject: [PATCH 08/18] Apply clang-format changes --- include/bout/field.hxx | 48 ++++++++++++++++++---------------------- include/bout/field3d.hxx | 3 +-- src/mesh/coordinates.cxx | 3 +-- 3 files changed, 24 insertions(+), 30 deletions(-) diff --git a/include/bout/field.hxx b/include/bout/field.hxx index ad4fa51b93..017fc39e3b 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -535,32 +535,28 @@ T pow(BoutReal lhs, const T& rhs, const std::string& rgn = "RGN_ALL") { #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]); } \ - 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; \ +#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 diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index 96e0722700..c162c12237 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -249,8 +249,7 @@ public: /// Returns a shallow copy without parallel slices Field3D withoutParallelSlices() const { - Field3D result{getMesh(), getLocation(), - getDirections(), getRegionID()}; + Field3D result{getMesh(), getLocation(), getDirections(), getRegionID()}; result.data = data; return result; } diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index f39c6c9133..e2348d5775 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1607,8 +1607,7 @@ Field3D Coordinates::Div_par(const Field3D& f, CELL_LOC outloc, f_B.yup(i) = f.yup(i) / Bxy_floc.yup(i); f_B.ydown(i) = f.ydown(i) / Bxy_floc.ydown(i); } - return setName(Bxy * Grad_par(f_B, outloc, method), - "Div_par({:s})", f.name); + return setName(Bxy * Grad_par(f_B, outloc, method), "Div_par({:s})", f.name); } ///////////////////////////////////////////////////////// From 0b976ec8131ed60321a32eb783061f6843de2dc2 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 15 Mar 2025 22:09:21 -0700 Subject: [PATCH 09/18] field ops: Remove FCI_AUTOMAGIC flag Always perform calculations in the yup/down slices if present in BOTH arguments. i.e. automagic always on for arithmetic operators. To avoid unnecessary work, use `withoutParallelSlices` to pass arguments without parallel slices. --- src/field/gen_fieldops.jinja | 16 +--- src/field/generated_fieldops.cxx | 136 +++++++++++-------------------- 2 files changed, 52 insertions(+), 100 deletions(-) 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..dd2a7add9a 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(); } @@ -329,7 +309,9 @@ Field3D& Field3D::operator*=(const Field2D& rhs) { // Delete existing parallel slices. We don't copy parallel slices, so any // that currently exist will be incorrect. - { clearParallelSlices(); } + { + clearParallelSlices(); + } checkData(*this); checkData(rhs); @@ -391,7 +373,9 @@ Field3D& Field3D::operator/=(const Field2D& rhs) { // Delete existing parallel slices. We don't copy parallel slices, so any // that currently exist will be incorrect. - { clearParallelSlices(); } + { + clearParallelSlices(); + } checkData(*this); checkData(rhs); @@ -453,7 +437,9 @@ Field3D& Field3D::operator+=(const Field2D& rhs) { // Delete existing parallel slices. We don't copy parallel slices, so any // that currently exist will be incorrect. - { clearParallelSlices(); } + { + clearParallelSlices(); + } checkData(*this); checkData(rhs); @@ -514,7 +500,9 @@ Field3D& Field3D::operator-=(const Field2D& rhs) { // Delete existing parallel slices. We don't copy parallel slices, so any // that currently exist will be incorrect. - { clearParallelSlices(); } + { + clearParallelSlices(); + } checkData(*this); checkData(rhs); @@ -640,15 +628,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 +654,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 +691,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 +718,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 +756,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 +782,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 +819,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 +845,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 +2177,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 +2204,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 +2231,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 +2258,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]; From fd93d6901123911f20050d10dd8e43462004d5d5 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 17 Mar 2025 21:12:06 -0700 Subject: [PATCH 10/18] FCI parallel: Revert deletion of coordinate slices Added comment to explain why it's ok to use the yup/down slices of the metric components if they are read from the grid file. --- src/mesh/parallel/fci.cxx | 60 +++++++++++++++++++++++++-------------- 1 file changed, 39 insertions(+), 21 deletions(-) diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index f7ba0b9aed..17a6aaf6fd 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -109,28 +109,45 @@ void load_parallel_metric_components([[maybe_unused]] Coordinates* coords, #define LOAD_PAR(var, doZero) \ load_parallel_metric_component(#var, coords->var, offset, doZero) - // Only some components can be compared between planes - // because each poloidal X-Z plane has a separate coordinate system - LOAD_PAR(dy, false); + // 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(g11, false); + LOAD_PAR(g22, false); + LOAD_PAR(g33, false); + LOAD_PAR(g12, false); + 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(Bxy, false); - - // Other components can't have parallel slices - coords->g11.allowParallelSlices(false); - coords->g22.allowParallelSlices(false); - coords->g33.allowParallelSlices(false); - coords->g12.allowParallelSlices(false); - coords->g13.allowParallelSlices(false); - coords->g23.allowParallelSlices(false); - - coords->g_11.allowParallelSlices(false); - coords->g_33.allowParallelSlices(false); - coords->g_12.allowParallelSlices(false); - coords->g_13.allowParallelSlices(false); - coords->g_23.allowParallelSlices(false); - - coords->J.allowParallelSlices(false); - + LOAD_PAR(g_33, false); + LOAD_PAR(g_12, false); + LOAD_PAR(g_13, false); + LOAD_PAR(g_23, false); + + if (not LOAD_PAR(J, true)) { + auto g = + coords->g11.ynext(offset) * coords->g22.ynext(offset) * coords->g33.ynext(offset) + + 2.0 * coords->g12.ynext(offset) * coords->g13.ynext(offset) + * coords->g23.ynext(offset) + - coords->g11.ynext(offset) * coords->g23.ynext(offset) + * coords->g23.ynext(offset) + - coords->g22.ynext(offset) * coords->g13.ynext(offset) + * coords->g13.ynext(offset) + - coords->g33.ynext(offset) * coords->g12.ynext(offset) + * coords->g12.ynext(offset); + + const auto rgn = fmt::format("RGN_YPAR_{:+d}", offset); + // Check that g is positive + bout::checkPositive(g, "The determinant of g^ij", rgn); + auto J = 1. / sqrt(g); + auto& pcom = coords->J.ynext(offset); + BOUT_FOR(i, J.getRegion(rgn)) { pcom[i] = J[i]; } + } #undef LOAD_PAR #endif } @@ -420,6 +437,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)); } } From 192d10528aa59d8f05253e8e0a60362613c9eaec Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 19 Mar 2025 20:01:58 -0700 Subject: [PATCH 11/18] parallel_boundary_region: Fix ITER macro Iterated for `i < localmesh->ystart - abs_offset()` so when ystart is 1 and the boundary offset is 1, no boundary is applied. --- include/bout/parallel_boundary_region.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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; } From 8d79e504320704600804c58e484968f6bbf4a733 Mon Sep 17 00:00:00 2001 From: bendudson <219233+bendudson@users.noreply.github.com> Date: Thu, 20 Mar 2025 03:16:41 +0000 Subject: [PATCH 12/18] Apply clang-format changes --- src/field/generated_fieldops.cxx | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/src/field/generated_fieldops.cxx b/src/field/generated_fieldops.cxx index dd2a7add9a..1664668e5f 100644 --- a/src/field/generated_fieldops.cxx +++ b/src/field/generated_fieldops.cxx @@ -309,9 +309,7 @@ Field3D& Field3D::operator*=(const Field2D& rhs) { // Delete existing parallel slices. We don't copy parallel slices, so any // that currently exist will be incorrect. - { - clearParallelSlices(); - } + { clearParallelSlices(); } checkData(*this); checkData(rhs); @@ -373,9 +371,7 @@ Field3D& Field3D::operator/=(const Field2D& rhs) { // Delete existing parallel slices. We don't copy parallel slices, so any // that currently exist will be incorrect. - { - clearParallelSlices(); - } + { clearParallelSlices(); } checkData(*this); checkData(rhs); @@ -437,9 +433,7 @@ Field3D& Field3D::operator+=(const Field2D& rhs) { // Delete existing parallel slices. We don't copy parallel slices, so any // that currently exist will be incorrect. - { - clearParallelSlices(); - } + { clearParallelSlices(); } checkData(*this); checkData(rhs); @@ -500,9 +494,7 @@ Field3D& Field3D::operator-=(const Field2D& rhs) { // Delete existing parallel slices. We don't copy parallel slices, so any // that currently exist will be incorrect. - { - clearParallelSlices(); - } + { clearParallelSlices(); } checkData(*this); checkData(rhs); From e23b104c24170fd1f20f1525adb81430deaae601 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 21 Mar 2025 22:02:44 -0700 Subject: [PATCH 13/18] Fixes for FCI: Load dy, improve errors, use J in Div_par Reverting some changes back to David's version. --- src/mesh/coordinates.cxx | 16 +++++++++------- src/mesh/difops.cxx | 22 +++++++++++++++------- src/mesh/parallel/fci.cxx | 2 ++ 3 files changed, 26 insertions(+), 14 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index e2348d5775..bcbd141890 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1601,13 +1601,15 @@ Field3D Coordinates::Div_par(const Field3D& f, CELL_LOC outloc, } // 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 setName(Bxy * Grad_par(f_B, outloc, method), "Div_par({:s})", f.name); + //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 f722b85afb..541fdb8986 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -31,6 +31,7 @@ #include #include #include +#include #include #include #include @@ -267,15 +268,22 @@ Field3D Div_par(const Field3D& f, const Field3D& v) { 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])) - / (dy[i] * sqrt(g_22[i])); + 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 > 0 +#if CHECK > 2 if (!std::isfinite(result[i])) { - output.write("{} {} {} {}\n", f_up[i], v_up[i], f_down[i], v_down[i]); - output.write("{} {} {} {} {}\n", B[i], B_up[i], B_down[i], dy[i], sqrt(g_22[i])); - throw BoutException("Non-finite value in Div_"); + 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 } diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index 17a6aaf6fd..eeccd0e7b1 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -113,6 +113,8 @@ void load_parallel_metric_components([[maybe_unused]] Coordinates* coords, // 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); From f169b8c7de895fd2f2e0c4e78bdc6aa72afd0a5d Mon Sep 17 00:00:00 2001 From: bendudson <219233+bendudson@users.noreply.github.com> Date: Sat, 22 Mar 2025 05:06:27 +0000 Subject: [PATCH 14/18] Apply clang-format changes --- src/mesh/difops.cxx | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 541fdb8986..1537650b0e 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -268,9 +268,9 @@ Field3D Div_par(const Field3D& f, const Field3D& v) { 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])); + 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])) { @@ -279,11 +279,8 @@ Field3D Div_par(const Field3D& f, const Field3D& v) { " 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])); + 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 } From ce70d8684b351240fd72bb4a36d5cb7da577dbc1 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 26 Mar 2025 10:51:00 -0700 Subject: [PATCH 15/18] Remove FCI_AUTOMAGIC compile option Prefer to have consistent behavior so that code works and can be optimized without breaking. The field.withoutParallelSlices() method can be used to strip slices before calling an operation, if slices are not needed. --- CMakeLists.txt | 3 --- include/bout/field.hxx | 10 ++++------ include/bout/fv_ops.hxx | 5 +---- src/field/field3d.cxx | 17 ++++++++--------- 4 files changed, 13 insertions(+), 22 deletions(-) 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 017fc39e3b..81130fcad5 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -684,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 @@ -697,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) { @@ -711,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/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/src/field/field3d.cxx b/src/field/field3d.cxx index 74a6f0853c..72cff777ae 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,14 +360,12 @@ Field3D& Field3D::operator=(const BoutReal val) { TRACE("Field3D = BoutReal"); track(val, "operator="); -#if BOUT_USE_FCI_AUTOMAGIC if (isFci() && 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(); @@ -386,11 +383,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; } From 328880814c81c85edc625b4b3f8ac35b6d018104 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 26 Mar 2025 11:37:08 -0700 Subject: [PATCH 16/18] Remove stray #endif --- src/field/field3d.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 72cff777ae..7878bb46e8 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -369,7 +369,7 @@ Field3D& Field3D::operator=(const BoutReal val) { // Delete existing parallel slices. We don't copy parallel slices, so any // that currently exist will be incorrect. clearParallelSlices(); -#endif + resetRegion(); allocate(); From c840122139c609a7cbdf0e1909fd476fc58fc48a Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 26 Mar 2025 21:20:46 -0700 Subject: [PATCH 17/18] Field3D: Fix operator=(BoutReal) If the field has parallel slices, assign the value also to the slices. --- src/field/field3d.cxx | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 7878bb46e8..70569027c3 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -360,15 +360,12 @@ Field3D& Field3D::operator=(const BoutReal val) { TRACE("Field3D = BoutReal"); track(val, "operator="); - if (isFci() && hasParallelSlices()) { + if (hasParallelSlices()) { for (size_t i = 0; i < numberParallelSlices(); ++i) { yup(i) = val; ydown(i) = val; } } - // Delete existing parallel slices. We don't copy parallel slices, so any - // that currently exist will be incorrect. - clearParallelSlices(); resetRegion(); From 568ad422c094c3a918d69a281b7462c3f4ddeb1d Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 17 Nov 2025 12:54:30 +0100 Subject: [PATCH 18/18] applyParallelBoundary can be a no-op for non-fci --- src/field/field3d.cxx | 35 ++++++++++++++++++++++++++++++----- 1 file changed, 30 insertions(+), 5 deletions(-) diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 70569027c3..a124141634 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -544,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()) { @@ -557,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()) { @@ -570,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(); @@ -589,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(); @@ -611,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();