From 358347602cfd43673f4cc64f598440cd5417baeb Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 16 Dec 2021 12:26:52 +0100 Subject: [PATCH 0001/1827] Switch hermite_spline_xz to index Might help with performance --- include/interpolation_xz.hxx | 2 +- src/mesh/interpolation/hermite_spline_xz.cxx | 86 +++++++++---------- .../monotonic_hermite_spline_xz.cxx | 61 ++++++------- 3 files changed, 65 insertions(+), 84 deletions(-) diff --git a/include/interpolation_xz.hxx b/include/interpolation_xz.hxx index 47474ad39c..04269a830d 100644 --- a/include/interpolation_xz.hxx +++ b/include/interpolation_xz.hxx @@ -95,7 +95,7 @@ protected: /// This is protected rather than private so that it can be /// extended and used by HermiteSplineMonotonic - Tensor i_corner; // x-index of bottom-left grid point + Tensor> i_corner; // index of bottom-left grid point Tensor k_corner; // z-index of bottom-left grid point // Basis functions for cubic Hermite spline interpolation diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index f5ce3357cb..9f14d4b836 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -38,7 +38,6 @@ XZHermiteSpline::XZHermiteSpline(int y_offset, Mesh *mesh) // Initialise in order to avoid 'uninitialized value' errors from Valgrind when using // guard-cell values - i_corner = -1; k_corner = -1; // Allocate Field3D members @@ -55,6 +54,8 @@ XZHermiteSpline::XZHermiteSpline(int y_offset, Mesh *mesh) void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z, const std::string& region) { + const int ny = localmesh->LocalNy; + const int nz = localmesh->LocalNz; BOUT_FOR(i, delta_x.getRegion(region)) { const int x = i.x(); const int y = i.y(); @@ -65,29 +66,29 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z // The integer part of xt_prime, zt_prime are the indices of the cell // containing the field line end-point - i_corner(x, y, z) = static_cast(floor(delta_x(x, y, z))); + int i_corn = static_cast(floor(delta_x(x, y, z))); k_corner(x, y, z) = static_cast(floor(delta_z(x, y, z))); // t_x, t_z are the normalised coordinates \in [0,1) within the cell // calculated by taking the remainder of the floating point index - BoutReal t_x = delta_x(x, y, z) - static_cast(i_corner(x, y, z)); + BoutReal t_x = delta_x(x, y, z) - static_cast(i_corn); BoutReal t_z = delta_z(x, y, z) - static_cast(k_corner(x, y, z)); // NOTE: A (small) hack to avoid one-sided differences - if (i_corner(x, y, z) >= localmesh->xend) { - i_corner(x, y, z) = localmesh->xend - 1; + if (i_corn >= localmesh->xend) { + i_corn = localmesh->xend - 1; t_x = 1.0; } - if (i_corner(x, y, z) < localmesh->xstart) { - i_corner(x, y, z) = localmesh->xstart; + if (i_corn < localmesh->xstart) { + i_corn = localmesh->xstart; t_x = 0.0; } // Check that t_x and t_z are in range if ((t_x < 0.0) || (t_x > 1.0)) { throw BoutException( - "t_x={:e} out of range at ({:d},{:d},{:d}) (delta_x={:e}, i_corner={:d})", t_x, - x, y, z, delta_x(x, y, z), i_corner(x, y, z)); + "t_x={:e} out of range at ({:d},{:d},{:d}) (delta_x={:e}, i_corn={:d})", t_x, + x, y, z, delta_x(x, y, z), i_corn); } if ((t_z < 0.0) || (t_z > 1.0)) { @@ -96,17 +97,20 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z x, y, z, delta_z(x, y, z), k_corner(x, y, z)); } - h00_x(x, y, z) = (2. * t_x * t_x * t_x) - (3. * t_x * t_x) + 1.; - h00_z(x, y, z) = (2. * t_z * t_z * t_z) - (3. * t_z * t_z) + 1.; + i_corner[i] = SpecificInd( + (((i_corn * ny) + (y + y_offset)) * nz + k_corner(x, y, z)), ny, nz); - h01_x(x, y, z) = (-2. * t_x * t_x * t_x) + (3. * t_x * t_x); - h01_z(x, y, z) = (-2. * t_z * t_z * t_z) + (3. * t_z * t_z); + h00_x[i] = (2. * t_x * t_x * t_x) - (3. * t_x * t_x) + 1.; + h00_z[i] = (2. * t_z * t_z * t_z) - (3. * t_z * t_z) + 1.; - h10_x(x, y, z) = t_x * (1. - t_x) * (1. - t_x); - h10_z(x, y, z) = t_z * (1. - t_z) * (1. - t_z); + h01_x[i] = (-2. * t_x * t_x * t_x) + (3. * t_x * t_x); + h01_z[i] = (-2. * t_z * t_z * t_z) + (3. * t_z * t_z); - h11_x(x, y, z) = (t_x * t_x * t_x) - (t_x * t_x); - h11_z(x, y, z) = (t_z * t_z * t_z) - (t_z * t_z); + h10_x[i] = t_x * (1. - t_x) * (1. - t_x); + h10_z[i] = t_z * (1. - t_z) * (1. - t_z); + + h11_x[i] = (t_x * t_x * t_x) - (t_x * t_x); + h11_z[i] = (t_z * t_z * t_z) - (t_z * t_z); } } @@ -137,8 +141,8 @@ XZHermiteSpline::getWeightsForYApproximation(int i, int j, int k, int yoffset) { const int ncz = localmesh->LocalNz; const int k_mod = ((k_corner(i, j, k) % ncz) + ncz) % ncz; const int k_mod_m1 = (k_mod > 0) ? (k_mod - 1) : (ncz - 1); - const int k_mod_p1 = (k_mod + 1) % ncz; - const int k_mod_p2 = (k_mod + 2) % ncz; + const int k_mod_p1 = (k_mod == ncz) ? 0 : k_mod + 1; + const int k_mod_p2 = (k_mod_p1 == ncz) ? 0 : k_mod_p1 + 1; return {{i, j + yoffset, k_mod_m1, -0.5 * h10_z(i, j, k)}, {i, j + yoffset, k_mod, h00_z(i, j, k) - 0.5 * h11_z(i, j, k)}, @@ -183,45 +187,35 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region if (skip_mask(x, y, z)) continue; - // Due to lack of guard cells in z-direction, we need to ensure z-index - // wraps around - const int ncz = localmesh->LocalNz; - const int z_mod = ((k_corner(x, y, z) % ncz) + ncz) % ncz; - const int z_mod_p1 = (z_mod + 1) % ncz; + const auto iyp = i.yp(y_offset); - const int y_next = y + y_offset; + const auto ic = i_corner[i]; + const auto iczp = ic.zp(); + const auto icxp = ic.xp(); + const auto icxpzp = iczp.xp(); // Interpolate f in X at Z - const BoutReal f_z = f(i_corner(x, y, z), y_next, z_mod) * h00_x(x, y, z) - + f(i_corner(x, y, z) + 1, y_next, z_mod) * h01_x(x, y, z) - + fx(i_corner(x, y, z), y_next, z_mod) * h10_x(x, y, z) - + fx(i_corner(x, y, z) + 1, y_next, z_mod) * h11_x(x, y, z); + const BoutReal f_z = + f[ic] * h00_x[i] + f[icxp] * h01_x[i] + fx[ic] * h10_x[i] + fx[icxp] * h11_x[i]; // Interpolate f in X at Z+1 - const BoutReal f_zp1 = f(i_corner(x, y, z), y_next, z_mod_p1) * h00_x(x, y, z) - + f(i_corner(x, y, z) + 1, y_next, z_mod_p1) * h01_x(x, y, z) - + fx(i_corner(x, y, z), y_next, z_mod_p1) * h10_x(x, y, z) - + fx(i_corner(x, y, z) + 1, y_next, z_mod_p1) * h11_x(x, y, z); + const BoutReal f_zp1 = f[iczp] * h00_x[i] + f[icxpzp] * h01_x[i] + fx[iczp] * h10_x[i] + + fx[icxpzp] * h11_x[i]; // Interpolate fz in X at Z - const BoutReal fz_z = fz(i_corner(x, y, z), y_next, z_mod) * h00_x(x, y, z) - + fz(i_corner(x, y, z) + 1, y_next, z_mod) * h01_x(x, y, z) - + fxz(i_corner(x, y, z), y_next, z_mod) * h10_x(x, y, z) - + fxz(i_corner(x, y, z) + 1, y_next, z_mod) * h11_x(x, y, z); + const BoutReal fz_z = fz[ic] * h00_x[i] + fz[icxp] * h01_x[i] + fxz[ic] * h10_x[i] + + fxz[icxp] * h11_x[i]; // Interpolate fz in X at Z+1 - const BoutReal fz_zp1 = - fz(i_corner(x, y, z), y_next, z_mod_p1) * h00_x(x, y, z) - + fz(i_corner(x, y, z) + 1, y_next, z_mod_p1) * h01_x(x, y, z) - + fxz(i_corner(x, y, z), y_next, z_mod_p1) * h10_x(x, y, z) - + fxz(i_corner(x, y, z) + 1, y_next, z_mod_p1) * h11_x(x, y, z); + const BoutReal fz_zp1 = fz[iczp] * h00_x[i] + fz[icxpzp] * h01_x[i] + + fxz[iczp] * h10_x[i] + fxz[icxpzp] * h11_x[i]; // Interpolate in Z - f_interp(x, y_next, z) = +f_z * h00_z(x, y, z) + f_zp1 * h01_z(x, y, z) - + fz_z * h10_z(x, y, z) + fz_zp1 * h11_z(x, y, z); + f_interp[iyp] = + +f_z * h00_z[i] + f_zp1 * h01_z[i] + fz_z * h10_z[i] + fz_zp1 * h11_z[i]; - ASSERT2(std::isfinite(f_interp(x, y_next, z)) || x < localmesh->xstart - || x > localmesh->xend); + ASSERT2(std::isfinite(f_interp[iyp]) || i.x() < localmesh->xstart + || i.x() > localmesh->xend); } return f_interp; } diff --git a/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx b/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx index bcf402231b..b2cfdb9515 100644 --- a/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx +++ b/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx @@ -66,44 +66,35 @@ Field3D XZMonotonicHermiteSpline::interpolate(const Field3D& f, if (skip_mask(x, y, z)) continue; - // Due to lack of guard cells in z-direction, we need to ensure z-index - // wraps around - const int ncz = localmesh->LocalNz; - const int z_mod = ((k_corner(x, y, z) % ncz) + ncz) % ncz; - const int z_mod_p1 = (z_mod + 1) % ncz; + const auto iyp = i.yp(y_offset); - const int y_next = y + y_offset; + const auto ic = i_corner[i]; + const auto iczp = ic.zp(); + const auto icxp = ic.xp(); + const auto icxpzp = iczp.xp(); // Interpolate f in X at Z - const BoutReal f_z = f(i_corner(x, y, z), y_next, z_mod) * h00_x(x, y, z) - + f(i_corner(x, y, z) + 1, y_next, z_mod) * h01_x(x, y, z) - + fx(i_corner(x, y, z), y_next, z_mod) * h10_x(x, y, z) - + fx(i_corner(x, y, z) + 1, y_next, z_mod) * h11_x(x, y, z); + const BoutReal f_z = + f[ic] * h00_x[i] + f[icxp] * h01_x[i] + fx[ic] * h10_x[i] + fx[icxp] * h11_x[i]; // Interpolate f in X at Z+1 - const BoutReal f_zp1 = f(i_corner(x, y, z), y_next, z_mod_p1) * h00_x(x, y, z) - + f(i_corner(x, y, z) + 1, y_next, z_mod_p1) * h01_x(x, y, z) - + fx(i_corner(x, y, z), y_next, z_mod_p1) * h10_x(x, y, z) - + fx(i_corner(x, y, z) + 1, y_next, z_mod_p1) * h11_x(x, y, z); + const BoutReal f_zp1 = f[iczp] * h00_x[i] + f[icxpzp] * h01_x[i] + fx[iczp] * h10_x[i] + + fx[icxpzp] * h11_x[i]; // Interpolate fz in X at Z - const BoutReal fz_z = fz(i_corner(x, y, z), y_next, z_mod) * h00_x(x, y, z) - + fz(i_corner(x, y, z) + 1, y_next, z_mod) * h01_x(x, y, z) - + fxz(i_corner(x, y, z), y_next, z_mod) * h10_x(x, y, z) - + fxz(i_corner(x, y, z) + 1, y_next, z_mod) * h11_x(x, y, z); + const BoutReal fz_z = fz[ic] * h00_x[i] + fz[icxp] * h01_x[i] + fxz[ic] * h10_x[i] + + fxz[icxp] * h11_x[i]; // Interpolate fz in X at Z+1 - const BoutReal fz_zp1 = - fz(i_corner(x, y, z), y_next, z_mod_p1) * h00_x(x, y, z) - + fz(i_corner(x, y, z) + 1, y_next, z_mod_p1) * h01_x(x, y, z) - + fxz(i_corner(x, y, z), y_next, z_mod_p1) * h10_x(x, y, z) - + fxz(i_corner(x, y, z) + 1, y_next, z_mod_p1) * h11_x(x, y, z); + const BoutReal fz_zp1 = fz[iczp] * h00_x[i] + fz[icxpzp] * h01_x[i] + + fxz[iczp] * h10_x[i] + fxz[icxpzp] * h11_x[i]; // Interpolate in Z - BoutReal result = +f_z * h00_z(x, y, z) + f_zp1 * h01_z(x, y, z) - + fz_z * h10_z(x, y, z) + fz_zp1 * h11_z(x, y, z); + BoutReal result = + +f_z * h00_z[i] + f_zp1 * h01_z[i] + fz_z * h10_z[i] + fz_zp1 * h11_z[i]; - ASSERT2(std::isfinite(result) || x < localmesh->xstart || x > localmesh->xend); + ASSERT2(std::isfinite(result) || i.x() < localmesh->xstart + || i.x() > localmesh->xend); // Monotonicity // Force the interpolated result to be in the range of the @@ -111,18 +102,14 @@ Field3D XZMonotonicHermiteSpline::interpolate(const Field3D& f, // but also degrades accuracy near maxima and minima. // Perhaps should only impose near boundaries, since that is where // problems most obviously occur. - const BoutReal localmax = BOUTMAX(f(i_corner(x, y, z), y_next, z_mod), - f(i_corner(x, y, z) + 1, y_next, z_mod), - f(i_corner(x, y, z), y_next, z_mod_p1), - f(i_corner(x, y, z) + 1, y_next, z_mod_p1)); + const BoutReal localmax = BOUTMAX(f[ic], f[icxp], f[iczp], f[icxpzp]); - const BoutReal localmin = BOUTMIN(f(i_corner(x, y, z), y_next, z_mod), - f(i_corner(x, y, z) + 1, y_next, z_mod), - f(i_corner(x, y, z), y_next, z_mod_p1), - f(i_corner(x, y, z) + 1, y_next, z_mod_p1)); + const BoutReal localmin = BOUTMIN(f[ic], f[icxp], f[iczp], f[icxpzp]); - ASSERT2(std::isfinite(localmax) || x < localmesh->xstart || x > localmesh->xend); - ASSERT2(std::isfinite(localmin) || x < localmesh->xstart || x > localmesh->xend); + ASSERT2(std::isfinite(localmax) || i.x() < localmesh->xstart + || i.x() > localmesh->xend); + ASSERT2(std::isfinite(localmin) || i.x() < localmesh->xstart + || i.x() > localmesh->xend); if (result > localmax) { result = localmax; @@ -131,7 +118,7 @@ Field3D XZMonotonicHermiteSpline::interpolate(const Field3D& f, result = localmin; } - f_interp(x, y_next, z) = result; + f_interp[iyp] = result; } return f_interp; } From 4ec74efcead5f2aabdac00ea950f166d17c92176 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20Schw=C3=B6rer?= Date: Wed, 22 Sep 2021 11:28:45 +0200 Subject: [PATCH 0002/1827] Switch toward regions --- include/interpolation_xz.hxx | 53 +++++++++++++++---- include/mask.hxx | 10 ++++ src/mesh/interpolation/bilinear_xz.cxx | 12 ++--- src/mesh/interpolation/hermite_spline_xz.cxx | 18 ++----- src/mesh/interpolation/lagrange_4pt_xz.cxx | 15 ++---- .../monotonic_hermite_spline_xz.cxx | 7 +-- 6 files changed, 65 insertions(+), 50 deletions(-) diff --git a/include/interpolation_xz.hxx b/include/interpolation_xz.hxx index 04269a830d..7a58ce3346 100644 --- a/include/interpolation_xz.hxx +++ b/include/interpolation_xz.hxx @@ -37,24 +37,56 @@ const Field3D interpolate(const Field2D &f, const Field3D &delta_x, const Field3D interpolate(const Field2D &f, const Field3D &delta_x); class XZInterpolation { +public: + int y_offset; + protected: Mesh* localmesh{nullptr}; - // 3D vector of points to skip (true -> skip this point) - BoutMask skip_mask; + std::string region_name{""}; + std::shared_ptr> region{nullptr}; public: XZInterpolation(int y_offset = 0, Mesh* localmeshIn = nullptr) - : localmesh(localmeshIn == nullptr ? bout::globals::mesh : localmeshIn), - skip_mask(*localmesh, false), y_offset(y_offset) {} + : y_offset(y_offset), + localmesh(localmeshIn == nullptr ? bout::globals::mesh : localmeshIn), + region_name("RGN_ALL") {} XZInterpolation(const BoutMask &mask, int y_offset = 0, Mesh *mesh = nullptr) : XZInterpolation(y_offset, mesh) { - skip_mask = mask; + region = regionFromMask(mask, localmesh); } + XZInterpolation(const std::string& region_name, int y_offset = 0, Mesh* mesh = nullptr) + : y_offset(y_offset), localmesh(mesh), region_name(region_name) {} + XZInterpolation(std::shared_ptr> region, int y_offset = 0, + Mesh* mesh = nullptr) + : y_offset(y_offset), localmesh(mesh), region(region) {} virtual ~XZInterpolation() = default; - - void setMask(const BoutMask &mask) { skip_mask = mask; } + void setMask(const BoutMask& mask) { + region = regionFromMask(mask, localmesh); + region_name = ""; + } + void setRegion(const std::string& region_name) { + this->region_name = region_name; + this->region = nullptr; + } + void setRegion(const std::shared_ptr>& region) { + this->region_name = ""; + this->region = region; + } + Region getRegion() const { + if (region_name != "") { + return localmesh->getRegion(region_name); + } + ASSERT1(region != nullptr); + return *region; + } + Region getRegion(const std::string& region) const { + if (region != "" and region != "RGN_ALL") { + return getIntersection(localmesh->getRegion(region), getRegion()); + } + return getRegion(); + } virtual void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const std::string& region = "RGN_NOBNDRY") = 0; virtual void calcWeights(const Field3D& delta_x, const Field3D& delta_z, @@ -71,7 +103,6 @@ public: const std::string& region = "RGN_NOBNDRY") = 0; // Interpolate using the field at (x,y+y_offset,z), rather than (x,y,z) - int y_offset; void setYOffset(int offset) { y_offset = offset; } virtual std::vector @@ -119,7 +150,7 @@ public: XZHermiteSpline(int y_offset = 0, Mesh *mesh = nullptr); XZHermiteSpline(const BoutMask &mask, int y_offset = 0, Mesh *mesh = nullptr) : XZHermiteSpline(y_offset, mesh) { - skip_mask = mask; + region = regionFromMask(mask, localmesh); } void calcWeights(const Field3D& delta_x, const Field3D& delta_z, @@ -177,7 +208,7 @@ public: XZLagrange4pt(int y_offset = 0, Mesh *mesh = nullptr); XZLagrange4pt(const BoutMask &mask, int y_offset = 0, Mesh *mesh = nullptr) : XZLagrange4pt(y_offset, mesh) { - skip_mask = mask; + region = regionFromMask(mask, localmesh); } void calcWeights(const Field3D& delta_x, const Field3D& delta_z, @@ -210,7 +241,7 @@ public: XZBilinear(int y_offset = 0, Mesh *mesh = nullptr); XZBilinear(const BoutMask &mask, int y_offset = 0, Mesh *mesh = nullptr) : XZBilinear(y_offset, mesh) { - skip_mask = mask; + region = regionFromMask(mask, localmesh); } void calcWeights(const Field3D& delta_x, const Field3D& delta_z, diff --git a/include/mask.hxx b/include/mask.hxx index 8940edbb16..c26bf31d61 100644 --- a/include/mask.hxx +++ b/include/mask.hxx @@ -72,4 +72,14 @@ public: } }; +inline std::unique_ptr> regionFromMask(const BoutMask& mask, + const Mesh* mesh) { + std::vector indices; + for (auto i : mesh->getRegion("RGN_ALL")) { + if (not mask(i.x(), i.y(), i.z())) { + indices.push_back(i); + } + } + return std::make_unique>(indices); +} #endif //__MASK_H__ diff --git a/src/mesh/interpolation/bilinear_xz.cxx b/src/mesh/interpolation/bilinear_xz.cxx index 7819fafe6f..1869df3218 100644 --- a/src/mesh/interpolation/bilinear_xz.cxx +++ b/src/mesh/interpolation/bilinear_xz.cxx @@ -45,14 +45,11 @@ XZBilinear::XZBilinear(int y_offset, Mesh *mesh) void XZBilinear::calcWeights(const Field3D& delta_x, const Field3D& delta_z, const std::string& region) { - BOUT_FOR(i, delta_x.getRegion(region)) { + BOUT_FOR(i, getRegion(region)) { const int x = i.x(); const int y = i.y(); const int z = i.z(); - if (skip_mask(x, y, z)) - continue; - // The integer part of xt_prime, zt_prime are the indices of the cell // containing the field line end-point i_corner(x, y, z) = static_cast(floor(delta_x(x, y, z))); @@ -87,7 +84,7 @@ void XZBilinear::calcWeights(const Field3D& delta_x, const Field3D& delta_z, void XZBilinear::calcWeights(const Field3D& delta_x, const Field3D& delta_z, const BoutMask& mask, const std::string& region) { - skip_mask = mask; + setMask(mask); calcWeights(delta_x, delta_z, region); } @@ -95,14 +92,11 @@ Field3D XZBilinear::interpolate(const Field3D& f, const std::string& region) con ASSERT1(f.getMesh() == localmesh); Field3D f_interp{emptyFrom(f)}; - BOUT_FOR(i, f.getRegion(region)) { + BOUT_FOR(i, this->getRegion(region)) { const int x = i.x(); const int y = i.y(); const int z = i.z(); - if (skip_mask(x, y, z)) - continue; - const int y_next = y + y_offset; // Due to lack of guard cells in z-direction, we need to ensure z-index // wraps around diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index 9f14d4b836..a682c58839 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -56,14 +56,11 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z const int ny = localmesh->LocalNy; const int nz = localmesh->LocalNz; - BOUT_FOR(i, delta_x.getRegion(region)) { + BOUT_FOR(i, getRegion(region)) { const int x = i.x(); const int y = i.y(); const int z = i.z(); - if (skip_mask(x, y, z)) - continue; - // The integer part of xt_prime, zt_prime are the indices of the cell // containing the field line end-point int i_corn = static_cast(floor(delta_x(x, y, z))); @@ -116,7 +113,7 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z, const BoutMask& mask, const std::string& region) { - skip_mask = mask; + setMask(mask); calcWeights(delta_x, delta_z, region); } @@ -179,16 +176,7 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region localmesh->wait(h); } - BOUT_FOR(i, f.getRegion(region)) { - const int x = i.x(); - const int y = i.y(); - const int z = i.z(); - - if (skip_mask(x, y, z)) - continue; - - const auto iyp = i.yp(y_offset); - + BOUT_FOR(i, getRegion(region)) { const auto ic = i_corner[i]; const auto iczp = ic.zp(); const auto icxp = ic.xp(); diff --git a/src/mesh/interpolation/lagrange_4pt_xz.cxx b/src/mesh/interpolation/lagrange_4pt_xz.cxx index 3a5de28e59..7c79a3f713 100644 --- a/src/mesh/interpolation/lagrange_4pt_xz.cxx +++ b/src/mesh/interpolation/lagrange_4pt_xz.cxx @@ -39,15 +39,12 @@ XZLagrange4pt::XZLagrange4pt(int y_offset, Mesh *mesh) void XZLagrange4pt::calcWeights(const Field3D& delta_x, const Field3D& delta_z, const std::string& region) { - - BOUT_FOR(i, delta_x.getRegion(region)) { + const auto curregion = getRegion(region); + BOUT_FOR(i, curregion) { const int x = i.x(); const int y = i.y(); const int z = i.z(); - if (skip_mask(x, y, z)) - continue; - // The integer part of xt_prime, zt_prime are the indices of the cell // containing the field line end-point i_corner(x, y, z) = static_cast(floor(delta_x(x, y, z))); @@ -80,7 +77,7 @@ void XZLagrange4pt::calcWeights(const Field3D& delta_x, const Field3D& delta_z, void XZLagrange4pt::calcWeights(const Field3D& delta_x, const Field3D& delta_z, const BoutMask& mask, const std::string& region) { - skip_mask = mask; + setMask(mask); calcWeights(delta_x, delta_z, region); } @@ -89,14 +86,12 @@ Field3D XZLagrange4pt::interpolate(const Field3D& f, const std::string& region) ASSERT1(f.getMesh() == localmesh); Field3D f_interp{emptyFrom(f)}; - BOUT_FOR(i, f.getRegion(region)) { + const auto curregion{getRegion(region)}; + BOUT_FOR(i, curregion) { const int x = i.x(); const int y = i.y(); const int z = i.z(); - if (skip_mask(x, y, z)) - continue; - const int jx = i_corner(x, y, z); const int jx2mnew = (jx == 0) ? 0 : (jx - 1); const int jxpnew = jx + 1; diff --git a/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx b/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx index b2cfdb9515..47eeb2df20 100644 --- a/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx +++ b/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx @@ -22,7 +22,7 @@ #include "globals.hxx" #include "interpolation_xz.hxx" -#include "output.hxx" +//#include "output.hxx" #include "bout/index_derivs_interface.hxx" #include "bout/mesh.hxx" @@ -58,14 +58,11 @@ Field3D XZMonotonicHermiteSpline::interpolate(const Field3D& f, localmesh->wait(h); } - BOUT_FOR(i, f.getRegion(region)) { + BOUT_FOR(i, getRegion(region)) { const int x = i.x(); const int y = i.y(); const int z = i.z(); - if (skip_mask(x, y, z)) - continue; - const auto iyp = i.yp(y_offset); const auto ic = i_corner[i]; From d0960204c870dc52cfe5315de22a8e5a4607b4eb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20Schw=C3=B6rer?= Date: Fri, 24 Sep 2021 11:19:50 +0200 Subject: [PATCH 0003/1827] Switch to regions for FCI regions --- include/bout/region.hxx | 23 +++++ include/interpolation_xz.hxx | 14 ++- include/mask.hxx | 1 + include/utils.hxx | 9 ++ src/mesh/interpolation/bilinear_xz.cxx | 6 +- src/mesh/interpolation/lagrange_4pt_xz.cxx | 2 +- .../monotonic_hermite_spline_xz.cxx | 3 +- src/mesh/parallel/fci.cxx | 89 +++++++++---------- src/mesh/parallel/fci.hxx | 7 +- 9 files changed, 98 insertions(+), 56 deletions(-) diff --git a/include/bout/region.hxx b/include/bout/region.hxx index f84c058814..4d0cb51159 100644 --- a/include/bout/region.hxx +++ b/include/bout/region.hxx @@ -51,6 +51,7 @@ #include "bout_types.hxx" #include "bout/assert.hxx" #include "bout/openmpwrap.hxx" +class BoutMask; /// The MAXREGIONBLOCKSIZE value can be tuned to try to optimise /// performance on specific hardware. It determines what the largest @@ -644,6 +645,28 @@ public: return *this; // To allow command chaining }; + /// Return a new region equivalent to *this but with indices contained + /// in mask Region removed + Region mask(const BoutMask& mask) { + // Get the current set of indices that we're going to mask and then + // use to create the result region. + auto currentIndices = getIndices(); + + // Lambda that returns true/false depending if the passed value is in maskIndices + // With C++14 T can be auto instead + auto isInVector = [&](T val) { return mask[val]; }; + + // Erase elements of currentIndices that are in maskIndices + currentIndices.erase( + std::remove_if(std::begin(currentIndices), std::end(currentIndices), isInVector), + std::end(currentIndices)); + + // Update indices + setIndices(currentIndices); + + return *this; // To allow command chaining + }; + /// Returns a new region including only indices contained in both /// this region and the other. Region getIntersection(const Region& otherRegion) { diff --git a/include/interpolation_xz.hxx b/include/interpolation_xz.hxx index 7a58ce3346..915b5a8478 100644 --- a/include/interpolation_xz.hxx +++ b/include/interpolation_xz.hxx @@ -49,8 +49,7 @@ protected: public: XZInterpolation(int y_offset = 0, Mesh* localmeshIn = nullptr) : y_offset(y_offset), - localmesh(localmeshIn == nullptr ? bout::globals::mesh : localmeshIn), - region_name("RGN_ALL") {} + localmesh(localmeshIn == nullptr ? bout::globals::mesh : localmeshIn) {} XZInterpolation(const BoutMask &mask, int y_offset = 0, Mesh *mesh = nullptr) : XZInterpolation(y_offset, mesh) { region = regionFromMask(mask, localmesh); @@ -74,6 +73,10 @@ public: this->region_name = ""; this->region = region; } + void setRegion(const Region& region) { + this->region_name = ""; + this->region = std::make_shared>(region); + } Region getRegion() const { if (region_name != "") { return localmesh->getRegion(region_name); @@ -82,9 +85,14 @@ public: return *region; } Region getRegion(const std::string& region) const { + const bool has_region = region_name != "" or this->region != nullptr; if (region != "" and region != "RGN_ALL") { - return getIntersection(localmesh->getRegion(region), getRegion()); + if (has_region) { + return getIntersection(localmesh->getRegion(region), getRegion()); + } + return localmesh->getRegion(region); } + ASSERT1(has_region); return getRegion(); } virtual void calcWeights(const Field3D& delta_x, const Field3D& delta_z, diff --git a/include/mask.hxx b/include/mask.hxx index c26bf31d61..6113e94d63 100644 --- a/include/mask.hxx +++ b/include/mask.hxx @@ -70,6 +70,7 @@ public: inline const bool& operator()(int jx, int jy, int jz) const { return mask(jx, jy, jz); } + inline const bool& operator[](const Ind3D& i) const { return mask[i]; } }; inline std::unique_ptr> regionFromMask(const BoutMask& mask, diff --git a/include/utils.hxx b/include/utils.hxx index 9de4628358..ea293f7bf8 100644 --- a/include/utils.hxx +++ b/include/utils.hxx @@ -38,6 +38,7 @@ #include "bout/array.hxx" #include "bout/assert.hxx" #include "bout/build_config.hxx" +#include "bout/region.hxx" #include #include @@ -348,6 +349,14 @@ public: return data[(i1*n2+i2)*n3 + i3]; } + const T& operator[](Ind3D i) const { + // ny and nz are private :-( + // ASSERT2(i.nz == n3); + // ASSERT2(i.ny == n2); + ASSERT2(0 <= i.ind && i.ind < n1 * n2 * n3); + return data[i.ind]; + } + Tensor& operator=(const T&val){ for(auto &i: data){ i = val; diff --git a/src/mesh/interpolation/bilinear_xz.cxx b/src/mesh/interpolation/bilinear_xz.cxx index 1869df3218..e36527e765 100644 --- a/src/mesh/interpolation/bilinear_xz.cxx +++ b/src/mesh/interpolation/bilinear_xz.cxx @@ -45,7 +45,8 @@ XZBilinear::XZBilinear(int y_offset, Mesh *mesh) void XZBilinear::calcWeights(const Field3D& delta_x, const Field3D& delta_z, const std::string& region) { - BOUT_FOR(i, getRegion(region)) { + const auto curregion{getRegion(region)}; + BOUT_FOR(i, curregion) { const int x = i.x(); const int y = i.y(); const int z = i.z(); @@ -92,7 +93,8 @@ Field3D XZBilinear::interpolate(const Field3D& f, const std::string& region) con ASSERT1(f.getMesh() == localmesh); Field3D f_interp{emptyFrom(f)}; - BOUT_FOR(i, this->getRegion(region)) { + const auto curregion{getRegion(region)}; + BOUT_FOR(i, curregion) { const int x = i.x(); const int y = i.y(); const int z = i.z(); diff --git a/src/mesh/interpolation/lagrange_4pt_xz.cxx b/src/mesh/interpolation/lagrange_4pt_xz.cxx index 7c79a3f713..caf4ce45eb 100644 --- a/src/mesh/interpolation/lagrange_4pt_xz.cxx +++ b/src/mesh/interpolation/lagrange_4pt_xz.cxx @@ -39,7 +39,7 @@ XZLagrange4pt::XZLagrange4pt(int y_offset, Mesh *mesh) void XZLagrange4pt::calcWeights(const Field3D& delta_x, const Field3D& delta_z, const std::string& region) { - const auto curregion = getRegion(region); + const auto curregion{getRegion(region)}; BOUT_FOR(i, curregion) { const int x = i.x(); const int y = i.y(); diff --git a/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx b/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx index 47eeb2df20..fbffc0b1fd 100644 --- a/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx +++ b/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx @@ -58,7 +58,8 @@ Field3D XZMonotonicHermiteSpline::interpolate(const Field3D& f, localmesh->wait(h); } - BOUT_FOR(i, getRegion(region)) { + const auto curregion{getRegion(region)}; + BOUT_FOR(i, curregion) { const int x = i.x(); const int y = i.y(); const int z = i.z(); diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index c25765852e..77f34fd282 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -50,7 +50,8 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& dy, Options& options, int offset_, BoundaryRegionPar* inner_boundary, BoundaryRegionPar* outer_boundary, bool zperiodic) - : map_mesh(mesh), offset(offset_), boundary_mask(map_mesh), + : map_mesh(mesh), offset(offset_), + region_no_boundary(map_mesh.getRegion("RGN_NOBNDRY")), corner_boundary_mask(map_mesh) { TRACE("Creating FCIMAP for direction {:d}", offset); @@ -156,6 +157,7 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& dy, Options& options, const int ncz = map_mesh.LocalNz; + BoutMask to_remove(map_mesh); // Serial loop because call to BoundaryRegionPar::addPoint // (probably?) can't be done in parallel BOUT_FOR_SERIAL(i, xt_prime.getRegion("RGN_NOBNDRY")) { @@ -185,7 +187,7 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& dy, Options& options, // indices (forward/backward_xt_prime and forward/backward_zt_prime) // are set to -1 - boundary_mask(x, y, z) = true; + to_remove(x, y, z) = true; // Need to specify the index of the boundary intersection, but // this may not be defined in general. @@ -200,13 +202,11 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& dy, Options& options, // and the gradients dR/dx etc. are evaluated at (x,y,z) // Cache the offsets - const auto i_xp = i.xp(); - const auto i_xm = i.xm(); const auto i_zp = i.zp(); const auto i_zm = i.zm(); - const BoutReal dR_dx = 0.5 * (R[i_xp] - R[i_xm]); - const BoutReal dZ_dx = 0.5 * (Z[i_xp] - Z[i_xm]); + const BoutReal dR_dx = 0.5 * (R[i.xp()] - R[i.xm()]); + const BoutReal dZ_dx = 0.5 * (Z[i.xp()] - Z[i.xm()]); BoutReal dR_dz, dZ_dz; // Handle the edge cases in Z @@ -241,8 +241,17 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& dy, Options& options, PI // Right-angle intersection ); } - - interp->setMask(boundary_mask); + region_no_boundary = region_no_boundary.mask(to_remove); + + const auto region = fmt::format("RGN_YPAR_{:+d}", offset); + if (not map_mesh.hasRegion3D(region)) { + // The valid region for this slice + map_mesh.addRegion3D(region, + Region(map_mesh.xstart, map_mesh.xend, + map_mesh.ystart+offset, map_mesh.yend+offset, + 0, map_mesh.LocalNz-1, + map_mesh.LocalNy, map_mesh.LocalNz)); + } } Field3D FCIMap::integrate(Field3D &f) const { @@ -265,45 +274,33 @@ Field3D FCIMap::integrate(Field3D &f) const { int nz = map_mesh.LocalNz; - for(int x = map_mesh.xstart; x <= map_mesh.xend; x++) { - for(int y = map_mesh.ystart; y <= map_mesh.yend; y++) { - - int ynext = y+offset; - - for(int z = 0; z < nz; z++) { - if (boundary_mask(x,y,z)) - continue; - - int zm = z - 1; - if (z == 0) { - zm = nz-1; - } - - BoutReal f_c = centre(x,ynext,z); - - if (corner_boundary_mask(x, y, z) || corner_boundary_mask(x - 1, y, z) || - corner_boundary_mask(x, y, zm) || corner_boundary_mask(x - 1, y, zm) || - (x == map_mesh.xstart)) { - // One of the corners leaves the domain. - // Use the cell centre value, since boundary conditions are not - // currently applied to corners. - result(x, ynext, z) = f_c; - - } else { - BoutReal f_pp = corner(x, ynext, z); // (x+1/2, z+1/2) - BoutReal f_mp = corner(x - 1, ynext, z); // (x-1/2, z+1/2) - BoutReal f_pm = corner(x, ynext, zm); // (x+1/2, z-1/2) - BoutReal f_mm = corner(x - 1, ynext, zm); // (x-1/2, z-1/2) - - // This uses a simple weighted average of centre and corners - // A more sophisticated approach might be to use e.g. Gauss-Lobatto points - // which would include cell edges and corners - result(x, ynext, z) = 0.5 * (f_c + 0.25 * (f_pp + f_mp + f_pm + f_mm)); - - ASSERT2(std::isfinite(result(x,ynext,z))); - } - } + BOUT_FOR(i, region_no_boundary) { + const auto inext = i.yp(offset); + BoutReal f_c = centre[inext]; + const auto izm = i.zm(); + const int x = i.x(); + const int y = i.y(); + const int z = i.z(); + const int zm = izm.z(); + if (corner_boundary_mask(x, y, z) || corner_boundary_mask(x - 1, y, z) + || corner_boundary_mask(x, y, zm) || corner_boundary_mask(x - 1, y, zm) + || (x == map_mesh.xstart)) { + // One of the corners leaves the domain. + // Use the cell centre value, since boundary conditions are not + // currently applied to corners. + result[inext] = f_c; + } else { + BoutReal f_pp = corner[inext]; // (x+1/2, z+1/2) + BoutReal f_mp = corner[inext.xm()]; // (x-1/2, z+1/2) + BoutReal f_pm = corner[inext.zm()]; // (x+1/2, z-1/2) + BoutReal f_mm = corner[inext.xm().zm()]; // (x-1/2, z-1/2) + + // This uses a simple weighted average of centre and corners + // A more sophisticated approach might be to use e.g. Gauss-Lobatto points + // which would include cell edges and corners + result[inext] = 0.5 * (f_c + 0.25 * (f_pp + f_mp + f_pm + f_mm)); } + ASSERT2(finite(result[inext])); } return result; } diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx index 3ecd964bfa..ef7c98693e 100644 --- a/src/mesh/parallel/fci.hxx +++ b/src/mesh/parallel/fci.hxx @@ -54,11 +54,12 @@ public: /// Direction of map const int offset; - /// boundary mask - has the field line left the domain - BoutMask boundary_mask; + /// region containing all points where the field line has not left the + /// domain + Region region_no_boundary; /// If any of the integration area has left the domain BoutMask corner_boundary_mask; - + Field3D interpolate(Field3D& f) const { ASSERT1(&map_mesh == f.getMesh()); return interp->interpolate(f); From d2a81bf142eea4e05fe5bb9c925290586a05fde6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20Schw=C3=B6rer?= Date: Fri, 19 Nov 2021 13:38:21 +0100 Subject: [PATCH 0004/1827] Use region_id in interpolation Otherwise the regions aren't cached. --- include/interpolation_xz.hxx | 61 +++++++++++++++++------------------- include/mask.hxx | 4 +-- 2 files changed, 30 insertions(+), 35 deletions(-) diff --git a/include/interpolation_xz.hxx b/include/interpolation_xz.hxx index 915b5a8478..df58d5b70e 100644 --- a/include/interpolation_xz.hxx +++ b/include/interpolation_xz.hxx @@ -43,8 +43,7 @@ public: protected: Mesh* localmesh{nullptr}; - std::string region_name{""}; - std::shared_ptr> region{nullptr}; + int region_id{-1}; public: XZInterpolation(int y_offset = 0, Mesh* localmeshIn = nullptr) @@ -52,48 +51,44 @@ public: localmesh(localmeshIn == nullptr ? bout::globals::mesh : localmeshIn) {} XZInterpolation(const BoutMask &mask, int y_offset = 0, Mesh *mesh = nullptr) : XZInterpolation(y_offset, mesh) { - region = regionFromMask(mask, localmesh); + setMask(mask); } XZInterpolation(const std::string& region_name, int y_offset = 0, Mesh* mesh = nullptr) - : y_offset(y_offset), localmesh(mesh), region_name(region_name) {} - XZInterpolation(std::shared_ptr> region, int y_offset = 0, + : y_offset(y_offset), localmesh(mesh), region_id(localmesh->getRegionID(region_name)) {} + XZInterpolation(const Region& region, int y_offset = 0, Mesh* mesh = nullptr) - : y_offset(y_offset), localmesh(mesh), region(region) {} + : y_offset(y_offset), localmesh(mesh){ + setRegion(region); + } virtual ~XZInterpolation() = default; void setMask(const BoutMask& mask) { - region = regionFromMask(mask, localmesh); - region_name = ""; + setRegion(regionFromMask(mask, localmesh)); } void setRegion(const std::string& region_name) { - this->region_name = region_name; - this->region = nullptr; - } - void setRegion(const std::shared_ptr>& region) { - this->region_name = ""; - this->region = region; + this->region_id = localmesh->getRegionID(region_name); } void setRegion(const Region& region) { - this->region_name = ""; - this->region = std::make_shared>(region); + std::string name; + int i=0; + do { + name = fmt::format("unsec_reg_xz_interp_{:d}",i++); + } while (localmesh->hasRegion3D(name)); + localmesh->addRegion(name, region); + this->region_id = localmesh->getRegionID(name); } - Region getRegion() const { - if (region_name != "") { - return localmesh->getRegion(region_name); - } - ASSERT1(region != nullptr); - return *region; + const Region& getRegion() const { + ASSERT2(region_id != -1); + return localmesh->getRegion(region_id); } - Region getRegion(const std::string& region) const { - const bool has_region = region_name != "" or this->region != nullptr; - if (region != "" and region != "RGN_ALL") { - if (has_region) { - return getIntersection(localmesh->getRegion(region), getRegion()); - } + const Region& getRegion(const std::string& region) const { + if (region_id == -1) { return localmesh->getRegion(region); } - ASSERT1(has_region); - return getRegion(); + if (region == "" or region == "RGN_ALL"){ + return getRegion(); + } + return localmesh->getRegion(localmesh->getCommonRegion(localmesh->getRegionID(region), region_id)); } virtual void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const std::string& region = "RGN_NOBNDRY") = 0; @@ -158,7 +153,7 @@ public: XZHermiteSpline(int y_offset = 0, Mesh *mesh = nullptr); XZHermiteSpline(const BoutMask &mask, int y_offset = 0, Mesh *mesh = nullptr) : XZHermiteSpline(y_offset, mesh) { - region = regionFromMask(mask, localmesh); + setRegion(regionFromMask(mask, localmesh)); } void calcWeights(const Field3D& delta_x, const Field3D& delta_z, @@ -216,7 +211,7 @@ public: XZLagrange4pt(int y_offset = 0, Mesh *mesh = nullptr); XZLagrange4pt(const BoutMask &mask, int y_offset = 0, Mesh *mesh = nullptr) : XZLagrange4pt(y_offset, mesh) { - region = regionFromMask(mask, localmesh); + setRegion(regionFromMask(mask, localmesh)); } void calcWeights(const Field3D& delta_x, const Field3D& delta_z, @@ -249,7 +244,7 @@ public: XZBilinear(int y_offset = 0, Mesh *mesh = nullptr); XZBilinear(const BoutMask &mask, int y_offset = 0, Mesh *mesh = nullptr) : XZBilinear(y_offset, mesh) { - region = regionFromMask(mask, localmesh); + setRegion(regionFromMask(mask, localmesh)); } void calcWeights(const Field3D& delta_x, const Field3D& delta_z, diff --git a/include/mask.hxx b/include/mask.hxx index 6113e94d63..96d2c99ac3 100644 --- a/include/mask.hxx +++ b/include/mask.hxx @@ -73,7 +73,7 @@ public: inline const bool& operator[](const Ind3D& i) const { return mask[i]; } }; -inline std::unique_ptr> regionFromMask(const BoutMask& mask, +inline Region regionFromMask(const BoutMask& mask, const Mesh* mesh) { std::vector indices; for (auto i : mesh->getRegion("RGN_ALL")) { @@ -81,6 +81,6 @@ inline std::unique_ptr> regionFromMask(const BoutMask& mask, indices.push_back(i); } } - return std::make_unique>(indices); + return Region{indices}; } #endif //__MASK_H__ From ea5640f1b661ab2978deff9084d56f3184b513ad Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 16 Dec 2021 12:26:52 +0100 Subject: [PATCH 0005/1827] Switch hermite_spline_xz to index Might help with performance --- src/mesh/interpolation/hermite_spline_xz.cxx | 2 ++ src/mesh/interpolation/monotonic_hermite_spline_xz.cxx | 4 ---- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index a682c58839..1319f7532d 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -177,6 +177,8 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region } BOUT_FOR(i, getRegion(region)) { + const auto iyp = i.yp(y_offset); + const auto ic = i_corner[i]; const auto iczp = ic.zp(); const auto icxp = ic.xp(); diff --git a/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx b/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx index fbffc0b1fd..e0cdf91ac8 100644 --- a/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx +++ b/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx @@ -60,10 +60,6 @@ Field3D XZMonotonicHermiteSpline::interpolate(const Field3D& f, const auto curregion{getRegion(region)}; BOUT_FOR(i, curregion) { - const int x = i.x(); - const int y = i.y(); - const int z = i.z(); - const auto iyp = i.yp(y_offset); const auto ic = i_corner[i]; From 089ddfa8413948f1fa8b1c759c64b05ce77e2468 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 16 Nov 2022 10:33:57 +0100 Subject: [PATCH 0006/1827] Improve hermitesplinesXZ Precalculate the matrix. This hard-codes the DDX and DDZ derivative to C2. Rather than taking derivatives, this first does the matrix-matrix multiplication and then does only a single matrix-vector operation, rather than 4. Retains a switch to go back to the old version. --- include/interpolation_xz.hxx | 2 + src/mesh/interpolation/hermite_spline_xz.cxx | 161 ++++++++++++++++--- 2 files changed, 142 insertions(+), 21 deletions(-) diff --git a/include/interpolation_xz.hxx b/include/interpolation_xz.hxx index df58d5b70e..e6eb8c3078 100644 --- a/include/interpolation_xz.hxx +++ b/include/interpolation_xz.hxx @@ -147,6 +147,8 @@ protected: Field3D h10_z; Field3D h11_z; + std::vector newWeights; + public: XZHermiteSpline(Mesh *mesh = nullptr) : XZHermiteSpline(0, mesh) {} diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index 1319f7532d..08515c7056 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -49,6 +49,13 @@ XZHermiteSpline::XZHermiteSpline(int y_offset, Mesh *mesh) h01_z.allocate(); h10_z.allocate(); h11_z.allocate(); + + + newWeights.reserve(16); + for (int w=0; w<16;++w){ + newWeights.emplace_back(localmesh); + newWeights[w].allocate(); + } } void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z, @@ -108,6 +115,115 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z h11_x[i] = (t_x * t_x * t_x) - (t_x * t_x); h11_z[i] = (t_z * t_z * t_z) - (t_z * t_z); + +#define USE_NEW_WEIGHTS 1 +#if USE_NEW_WEIGHTS + + for (int w =0; w<16;++w){ + newWeights[w][i]=0; + } + // The distribution of our weights: + // 0 4 8 12 + // 1 5 9 13 + // 2 6 10 14 + // 3 7 11 15 + // e.g. 1 == ic.xm(); 4 == ic.zm(); 5 == ic; 7 == ic.zp(2); + + // f[ic] * h00_x[i] + f[icxp] * h01_x[i] + fx[ic] * h10_x[i] + fx[icxp] * h11_x[i]; + newWeights[5][i] += h00_x[i] * h00_z[i]; + newWeights[9][i] += h01_x[i] * h00_z[i]; + newWeights[9][i] += h10_x[i] * h00_z[i] / 2; + newWeights[1][i] -= h10_x[i] * h00_z[i] / 2; + newWeights[13][i] += h11_x[i] * h00_z[i] / 2; + newWeights[5][i] -= h11_x[i] * h00_z[i] / 2; + + // f[iczp] * h00_x[i] + f[icxpzp] * h01_x[i] + + // fx[iczp] * h10_x[i] + fx[icxpzp] * h11_x[i]; + newWeights[6][i] += h00_x[i] * h01_z[i]; + newWeights[10][i] += h01_x[i] * h01_z[i]; + newWeights[10][i] += h10_x[i] * h01_z[i] / 2; + newWeights[2][i] -= h10_x[i] * h01_z[i] / 2; + newWeights[14][i] += h11_x[i] * h01_z[i] / 2; + newWeights[6][i] -= h11_x[i] * h01_z[i] / 2; + + // fz[ic] * h00_x[i] + fz[icxp] * h01_x[i] + + // fxz[ic] * h10_x[i]+ fxz[icxp] * h11_x[i]; + newWeights[6][i] += h00_x[i] * h10_z[i] / 2; + newWeights[4][i] -= h00_x[i] * h10_z[i] / 2; + newWeights[10][i] += h01_x[i] * h10_z[i] / 2; + newWeights[8][i] -= h01_x[i] * h10_z[i] / 2; + newWeights[10][i] += h10_x[i] * h10_z[i] / 4; + newWeights[8][i] -= h10_x[i] * h10_z[i] / 4; + newWeights[2][i] -= h10_x[i] * h10_z[i] / 4; + newWeights[0][i] += h10_x[i] * h10_z[i] / 4; + newWeights[14][i] += h11_x[i] * h10_z[i] / 4; + newWeights[12][i] -= h11_x[i] * h10_z[i] / 4; + newWeights[6][i] -= h11_x[i] * h10_z[i] / 4; + newWeights[4][i] += h11_x[i] * h10_z[i] / 4; + + // fz[iczp] * h00_x[i] + fz[icxpzp] * h01_x[i] + + // fxz[iczp] * h10_x[i] + fxz[icxpzp] * h11_x[i]; + newWeights[7][i] += h00_x[i] * h11_z[i] / 2; + newWeights[5][i] -= h00_x[i] * h11_z[i] / 2; + newWeights[11][i] += h01_x[i] * h11_z[i] / 2; + newWeights[9][i] -= h01_x[i] * h11_z[i] / 2; + newWeights[11][i] += h10_x[i] * h11_z[i] / 4; + newWeights[9][i] -= h10_x[i] * h11_z[i] / 4; + newWeights[3][i] -= h10_x[i] * h11_z[i] / 4; + newWeights[1][i] += h10_x[i] * h11_z[i] / 4; + newWeights[15][i] += h11_x[i] * h11_z[i] / 4; + newWeights[13][i] -= h11_x[i] * h11_z[i] / 4; + newWeights[7][i] -= h11_x[i] * h11_z[i] / 4; + newWeights[5][i] += h11_x[i] * h11_z[i] / 4; + + + // // f[ic] * h00_x[i] + f[icxp] * h01_x[i] + fx[ic] * h10_x[i] + fx[icxp] * h11_x[i]; + // newWeights[5][i] += h00_x[i] * h00_z[i]; + // newWeights[9][i] += h01_x[i] * h00_z[i]; + // newWeights[9][i] += h10_x[i] * h00_z[i] / 2 / localmesh->dx[ic]; + // newWeights[1][i] -= h10_x[i] * h00_z[i] / 2 / localmesh->dx[ic]; + // newWeights[13][i] += h11_x[i] * h00_z[i] / 2 / localmesh->dx[ic.xp()]; + // newWeights[5][i] -= h11_x[i] * h00_z[i] / 2 / localmesh->dx[ic.xp()]; + + // // f[iczp] * h00_x[i] + f[icxpzp] * h01_x[i] + + // // fx[iczp] * h10_x[i] + fx[icxpzp] * h11_x[i]; + // newWeights[6][i] += h00_x[i] * h01_z[i]; + // newWeights[10][i] += h01_x[i] * h01_z[i]; + // newWeights[10][i] += h10_x[i] * h01_z[i] / 2/ localmesh->dx[ic.zp()]; + // newWeights[2][i] -= h10_x[i] * h01_z[i] / 2/ localmesh->dx[ic.zp()]; + // newWeights[14][i] += h11_x[i] * h01_z[i] / 2/ localmesh->dx[ic.zp().xp()]; + // newWeights[6][i] -= h11_x[i] * h01_z[i] / 2/ localmesh->dx[ic.zp().xp()]; + + // // fz[ic] * h00_x[i] + fz[icxp] * h01_x[i] + + // // fxz[ic] * h10_x[i]+ fxz[icxp] * h11_x[i]; + // newWeights[6][i] += h00_x[i] * h10_z[i] / 2 / localmesh->dz[ic]; + // newWeights[4][i] -= h00_x[i] * h10_z[i] / 2 / localmesh->dz[ic]; + // newWeights[10][i] += h01_x[i] * h10_z[i] / 2 / localmesh->dz[ic.xp()]; + // newWeights[8][i] -= h01_x[i] * h10_z[i] / 2 / localmesh->dz[ic.xp()]; + // newWeights[10][i] += h10_x[i] * h10_z[i] / 4 / localmesh->dz[ic] / localmesh->dx[ic]; + // newWeights[8][i] -= h10_x[i] * h10_z[i] / 4 / localmesh->dz[ic] / localmesh->dx[ic]; + // newWeights[2][i] -= h10_x[i] * h10_z[i] / 4 / localmesh->dz[ic] / localmesh->dx[ic]; + // newWeights[0][i] += h10_x[i] * h10_z[i] / 4 / localmesh->dz[ic] / localmesh->dx[ic]; + // newWeights[14][i] += h11_x[i] * h10_z[i] / 4 / localmesh->dz[ic.xp()] / localmesh->dx[ic.xp()]; + // newWeights[12][i] -= h11_x[i] * h10_z[i] / 4 / localmesh->dz[ic.xp()] / localmesh->dx[ic.xp()]; + // newWeights[6][i] -= h11_x[i] * h10_z[i] / 4 / localmesh->dz[ic.xp()] / localmesh->dx[ic.xp()]; + // newWeights[4][i] += h11_x[i] * h10_z[i] / 4 / localmesh->dz[ic.xp()] / localmesh->dx[ic.xp()]; + + // // fz[iczp] * h00_x[i] + fz[icxpzp] * h01_x[i] + + // // fxz[iczp] * h10_x[i] + fxz[icxpzp] * h11_x[i]; + // newWeights[7][i] += h00_x[i] * h11_z[i] / 2 / localmesh->dz[ic.zp()]; + // newWeights[5][i] -= h00_x[i] * h11_z[i] / 2 / localmesh->dz[ic.zp()]; + // newWeights[11][i] += h01_x[i] * h11_z[i] / 2 / localmesh->dz[ic.zp().xp()]; + // newWeights[9][i] -= h01_x[i] * h11_z[i] / 2 / localmesh->dz[ic.zp().xp()]; + // newWeights[11][i] += h10_x[i] * h11_z[i] / 4 / localmesh->dz[ic.zp()] / localmesh->dx[ic.zp()]; + // newWeights[9][i] -= h10_x[i] * h11_z[i] / 4 / localmesh->dz[ic.zp()] / localmesh->dx[ic.zp()]; + // newWeights[3][i] -= h10_x[i] * h11_z[i] / 4 / localmesh->dz[ic.zp()] / localmesh->dx[ic.zp()]; + // newWeights[1][i] += h10_x[i] * h11_z[i] / 4 / localmesh->dz[ic.zp()] / localmesh->dx[ic.zp()]; + // newWeights[15][i] += h11_x[i] * h11_z[i] / 4 / localmesh->dz[ic.zp().xp()] / localmesh->dx[ic.zp().xp()]; + // newWeights[13][i] -= h11_x[i] * h11_z[i] / 4 / localmesh->dz[ic.zp().xp()] / localmesh->dx[ic.zp().xp()]; + // newWeights[7][i] -= h11_x[i] * h11_z[i] / 4 / localmesh->dz[ic.zp().xp()] / localmesh->dx[ic.zp().xp()]; + // newWeights[5][i] += h11_x[i] * h11_z[i] / 4 / localmesh->dz[ic.zp().xp()] / localmesh->dx[ic.zp().xp()]; +#endif } } @@ -152,29 +268,31 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region ASSERT1(f.getMesh() == localmesh); Field3D f_interp{emptyFrom(f)}; + +#if USE_NEW_WEIGHTS + BOUT_FOR(i, getRegion(region)) { + auto ic = i_corner[i]; + auto iyp = i.yp(y_offset); + + f_interp[iyp]=0; + for (int w = 0; w < 4; ++w){ + f_interp[iyp] += newWeights[w*4+0][i] * f[ic.zm().xp(w-1)]; + f_interp[iyp] += newWeights[w*4+1][i] * f[ic.xp(w-1)]; + f_interp[iyp] += newWeights[w*4+2][i] * f[ic.zp().xp(w-1)]; + f_interp[iyp] += newWeights[w*4+3][i] * f[ic.zp(2).xp(w-1)]; + } + } + return f_interp; +#else // Derivatives are used for tension and need to be on dimensionless // coordinates - Field3D fx = bout::derivatives::index::DDX(f, CELL_DEFAULT, "DEFAULT"); - localmesh->communicateXZ(fx); - // communicate in y, but do not calculate parallel slices - { - auto h = localmesh->sendY(fx); - localmesh->wait(h); - } - Field3D fz = bout::derivatives::index::DDZ(f, CELL_DEFAULT, "DEFAULT", "RGN_ALL"); - localmesh->communicateXZ(fz); - // communicate in y, but do not calculate parallel slices - { - auto h = localmesh->sendY(fz); - localmesh->wait(h); - } - Field3D fxz = bout::derivatives::index::DDX(fz, CELL_DEFAULT, "DEFAULT"); - localmesh->communicateXZ(fxz); - // communicate in y, but do not calculate parallel slices - { - auto h = localmesh->sendY(fxz); - localmesh->wait(h); - } + const auto region2 = fmt::format("RGN_YPAR_{:+d}", y_offset); + // f has been communcated, and thus we can assume that the x-boundaries are + // also valid in the y-boundary. Thus the differentiated field needs no + // extra comms. + Field3D fx = bout::derivatives::index::DDX(f, CELL_DEFAULT, "DEFAULT", region2); + Field3D fz = bout::derivatives::index::DDZ(f, CELL_DEFAULT, "DEFAULT", region2); + Field3D fxz = bout::derivatives::index::DDZ(fx, CELL_DEFAULT, "DEFAULT", region2); BOUT_FOR(i, getRegion(region)) { const auto iyp = i.yp(y_offset); @@ -208,6 +326,7 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region || i.x() > localmesh->xend); } return f_interp; +# endif } Field3D XZHermiteSpline::interpolate(const Field3D& f, const Field3D& delta_x, From de7d1deab6a419bf974fb8df5074fba4514175e3 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 17 Nov 2022 10:36:29 +0100 Subject: [PATCH 0007/1827] Cleanup --- src/mesh/interpolation/hermite_spline_xz.cxx | 50 +------------------- 1 file changed, 1 insertion(+), 49 deletions(-) diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index 08515c7056..22da552125 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -158,7 +158,7 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z newWeights[0][i] += h10_x[i] * h10_z[i] / 4; newWeights[14][i] += h11_x[i] * h10_z[i] / 4; newWeights[12][i] -= h11_x[i] * h10_z[i] / 4; - newWeights[6][i] -= h11_x[i] * h10_z[i] / 4; + newWeights[6][i] -= h11_x[i] * h10_z[i] / 4; newWeights[4][i] += h11_x[i] * h10_z[i] / 4; // fz[iczp] * h00_x[i] + fz[icxpzp] * h01_x[i] + @@ -175,54 +175,6 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z newWeights[13][i] -= h11_x[i] * h11_z[i] / 4; newWeights[7][i] -= h11_x[i] * h11_z[i] / 4; newWeights[5][i] += h11_x[i] * h11_z[i] / 4; - - - // // f[ic] * h00_x[i] + f[icxp] * h01_x[i] + fx[ic] * h10_x[i] + fx[icxp] * h11_x[i]; - // newWeights[5][i] += h00_x[i] * h00_z[i]; - // newWeights[9][i] += h01_x[i] * h00_z[i]; - // newWeights[9][i] += h10_x[i] * h00_z[i] / 2 / localmesh->dx[ic]; - // newWeights[1][i] -= h10_x[i] * h00_z[i] / 2 / localmesh->dx[ic]; - // newWeights[13][i] += h11_x[i] * h00_z[i] / 2 / localmesh->dx[ic.xp()]; - // newWeights[5][i] -= h11_x[i] * h00_z[i] / 2 / localmesh->dx[ic.xp()]; - - // // f[iczp] * h00_x[i] + f[icxpzp] * h01_x[i] + - // // fx[iczp] * h10_x[i] + fx[icxpzp] * h11_x[i]; - // newWeights[6][i] += h00_x[i] * h01_z[i]; - // newWeights[10][i] += h01_x[i] * h01_z[i]; - // newWeights[10][i] += h10_x[i] * h01_z[i] / 2/ localmesh->dx[ic.zp()]; - // newWeights[2][i] -= h10_x[i] * h01_z[i] / 2/ localmesh->dx[ic.zp()]; - // newWeights[14][i] += h11_x[i] * h01_z[i] / 2/ localmesh->dx[ic.zp().xp()]; - // newWeights[6][i] -= h11_x[i] * h01_z[i] / 2/ localmesh->dx[ic.zp().xp()]; - - // // fz[ic] * h00_x[i] + fz[icxp] * h01_x[i] + - // // fxz[ic] * h10_x[i]+ fxz[icxp] * h11_x[i]; - // newWeights[6][i] += h00_x[i] * h10_z[i] / 2 / localmesh->dz[ic]; - // newWeights[4][i] -= h00_x[i] * h10_z[i] / 2 / localmesh->dz[ic]; - // newWeights[10][i] += h01_x[i] * h10_z[i] / 2 / localmesh->dz[ic.xp()]; - // newWeights[8][i] -= h01_x[i] * h10_z[i] / 2 / localmesh->dz[ic.xp()]; - // newWeights[10][i] += h10_x[i] * h10_z[i] / 4 / localmesh->dz[ic] / localmesh->dx[ic]; - // newWeights[8][i] -= h10_x[i] * h10_z[i] / 4 / localmesh->dz[ic] / localmesh->dx[ic]; - // newWeights[2][i] -= h10_x[i] * h10_z[i] / 4 / localmesh->dz[ic] / localmesh->dx[ic]; - // newWeights[0][i] += h10_x[i] * h10_z[i] / 4 / localmesh->dz[ic] / localmesh->dx[ic]; - // newWeights[14][i] += h11_x[i] * h10_z[i] / 4 / localmesh->dz[ic.xp()] / localmesh->dx[ic.xp()]; - // newWeights[12][i] -= h11_x[i] * h10_z[i] / 4 / localmesh->dz[ic.xp()] / localmesh->dx[ic.xp()]; - // newWeights[6][i] -= h11_x[i] * h10_z[i] / 4 / localmesh->dz[ic.xp()] / localmesh->dx[ic.xp()]; - // newWeights[4][i] += h11_x[i] * h10_z[i] / 4 / localmesh->dz[ic.xp()] / localmesh->dx[ic.xp()]; - - // // fz[iczp] * h00_x[i] + fz[icxpzp] * h01_x[i] + - // // fxz[iczp] * h10_x[i] + fxz[icxpzp] * h11_x[i]; - // newWeights[7][i] += h00_x[i] * h11_z[i] / 2 / localmesh->dz[ic.zp()]; - // newWeights[5][i] -= h00_x[i] * h11_z[i] / 2 / localmesh->dz[ic.zp()]; - // newWeights[11][i] += h01_x[i] * h11_z[i] / 2 / localmesh->dz[ic.zp().xp()]; - // newWeights[9][i] -= h01_x[i] * h11_z[i] / 2 / localmesh->dz[ic.zp().xp()]; - // newWeights[11][i] += h10_x[i] * h11_z[i] / 4 / localmesh->dz[ic.zp()] / localmesh->dx[ic.zp()]; - // newWeights[9][i] -= h10_x[i] * h11_z[i] / 4 / localmesh->dz[ic.zp()] / localmesh->dx[ic.zp()]; - // newWeights[3][i] -= h10_x[i] * h11_z[i] / 4 / localmesh->dz[ic.zp()] / localmesh->dx[ic.zp()]; - // newWeights[1][i] += h10_x[i] * h11_z[i] / 4 / localmesh->dz[ic.zp()] / localmesh->dx[ic.zp()]; - // newWeights[15][i] += h11_x[i] * h11_z[i] / 4 / localmesh->dz[ic.zp().xp()] / localmesh->dx[ic.zp().xp()]; - // newWeights[13][i] -= h11_x[i] * h11_z[i] / 4 / localmesh->dz[ic.zp().xp()] / localmesh->dx[ic.zp().xp()]; - // newWeights[7][i] -= h11_x[i] * h11_z[i] / 4 / localmesh->dz[ic.zp().xp()] / localmesh->dx[ic.zp().xp()]; - // newWeights[5][i] += h11_x[i] * h11_z[i] / 4 / localmesh->dz[ic.zp().xp()] / localmesh->dx[ic.zp().xp()]; #endif } } From 8a1c2306b9496a09d768fe4fe990352f5fb34936 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 2 Feb 2023 10:55:01 +0100 Subject: [PATCH 0008/1827] Enable splitting in X using PETSc --- include/interpolation_xz.hxx | 27 ++++ src/mesh/interpolation/hermite_spline_xz.cxx | 140 ++++++++++++++++++- 2 files changed, 161 insertions(+), 6 deletions(-) diff --git a/include/interpolation_xz.hxx b/include/interpolation_xz.hxx index e6eb8c3078..620d146edf 100644 --- a/include/interpolation_xz.hxx +++ b/include/interpolation_xz.hxx @@ -26,6 +26,15 @@ #include "mask.hxx" +#define USE_NEW_WEIGHTS 1 +#if BOUT_HAS_PETSC +#define HS_USE_PETSC 1 +#endif + +#ifdef HS_USE_PETSC +#include "bout/petsclib.hxx" +#endif + class Options; /// Interpolate a field onto a perturbed set of points @@ -149,6 +158,13 @@ protected: std::vector newWeights; +#if HS_USE_PETSC + PetscLib* petsclib; + bool isInit{false}; + Mat petscWeights; + Vec rhs, result; +#endif + public: XZHermiteSpline(Mesh *mesh = nullptr) : XZHermiteSpline(0, mesh) {} @@ -157,6 +173,17 @@ public: : XZHermiteSpline(y_offset, mesh) { setRegion(regionFromMask(mask, localmesh)); } + ~XZHermiteSpline() { +#if HS_USE_PETSC + if (isInit) { + MatDestroy(&petscWeights); + VecDestroy(&rhs); + VecDestroy(&result); + isInit = false; + delete petsclib; + } +#endif + } void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const std::string& region = "RGN_NOBNDRY") override; diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index 22da552125..1b63c84231 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -23,10 +23,84 @@ #include "globals.hxx" #include "interpolation_xz.hxx" #include "bout/index_derivs_interface.hxx" -#include "bout/mesh.hxx" +#include "../impls/bout/boutmesh.hxx" #include +class IndConverter { +public: + IndConverter(Mesh* mesh) + : mesh(dynamic_cast(mesh)), nxpe(mesh->getNXPE()), nype(mesh->getNYPE()), + xstart(mesh->xstart), ystart(mesh->ystart), zstart(0), + lnx(mesh->LocalNx - 2 * xstart), lny(mesh->LocalNy - 2 * ystart), + lnz(mesh->LocalNz - 2 * zstart) {} + // ix and iy are global indices + // iy is local + int fromMeshToGlobal(int ix, int iy, int iz) { + const int xstart = mesh->xstart; + const int lnx = mesh->LocalNx - xstart * 2; + // x-proc-id + int pex = divToNeg(ix - xstart, lnx); + if (pex < 0) { + pex = 0; + } + if (pex >= nxpe) { + pex = nxpe - 1; + } + const int zstart = 0; + const int lnz = mesh->LocalNz - zstart * 2; + // z-proc-id + // pez only for wrapping around ; later needs similar treatment than pey + const int pez = divToNeg(iz - zstart, lnz); + // y proc-id - y is already local + const int ystart = mesh->ystart; + const int lny = mesh->LocalNy - ystart * 2; + const int pey_offset = divToNeg(iy - ystart, lny); + int pey = pey_offset + mesh->getYProcIndex(); + while (pey < 0) { + pey += nype; + } + while (pey >= nype) { + pey -= nype; + } + ASSERT2(pex >= 0); + ASSERT2(pex < nxpe); + ASSERT2(pey >= 0); + ASSERT2(pey < nype); + return fromLocalToGlobal(ix - pex * lnx, iy - pey_offset * lny, iz - pez * lnz, pex, + pey, 0); + } + int fromLocalToGlobal(const int ilocalx, const int ilocaly, const int ilocalz) { + return fromLocalToGlobal(ilocalx, ilocaly, ilocalz, mesh->getXProcIndex(), + mesh->getYProcIndex(), 0); + } + int fromLocalToGlobal(const int ilocalx, const int ilocaly, const int ilocalz, + const int pex, const int pey, const int pez) { + ASSERT3(ilocalx >= 0); + ASSERT3(ilocaly >= 0); + ASSERT3(ilocalz >= 0); + const int ilocal = ((ilocalx * mesh->LocalNy) + ilocaly) * mesh->LocalNz + ilocalz; + const int ret = ilocal + + mesh->LocalNx * mesh->LocalNy * mesh->LocalNz + * ((pey * nxpe + pex) * nzpe + pez); + ASSERT3(ret >= 0); + ASSERT3(ret < nxpe * nype * mesh->LocalNx * mesh->LocalNy * mesh->LocalNz); + return ret; + } + +private: + // number of procs + BoutMesh* mesh; + const int nxpe; + const int nype; + const int nzpe{1}; + const int xstart, ystart, zstart; + const int lnx, lny, lnz; + static int divToNeg(const int n, const int d) { + return (n < 0) ? ((n - d + 1) / d) : (n / d); + } +}; + XZHermiteSpline::XZHermiteSpline(int y_offset, Mesh *mesh) : XZInterpolation(y_offset, mesh), h00_x(localmesh), h01_x(localmesh), h10_x(localmesh), h11_x(localmesh), @@ -50,12 +124,27 @@ XZHermiteSpline::XZHermiteSpline(int y_offset, Mesh *mesh) h10_z.allocate(); h11_z.allocate(); - +#if USE_NEW_WEIGHTS newWeights.reserve(16); for (int w=0; w<16;++w){ newWeights.emplace_back(localmesh); newWeights[w].allocate(); } +#ifdef HS_USE_PETSC + petsclib = new PetscLib( + &Options::root()["mesh:paralleltransform:xzinterpolation:hermitespline"]); + // MatCreate(MPI_Comm comm,Mat *A) + // MatCreate(MPI_COMM_WORLD, &petscWeights); + // MatSetSizes(petscWeights, m, m, M, M); + // PetscErrorCode MatCreateAIJ(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, + // PetscInt N, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt + //o_nnz[], Mat *A) + // MatSetSizes(Mat A,PetscInt m,PetscInt n,PetscInt M,PetscInt N) + const int m = mesh->LocalNx * mesh->LocalNy * mesh->LocalNz; + const int M = m * mesh->getNXPE() * mesh->getNYPE(); + MatCreateAIJ(MPI_COMM_WORLD, m, m, M, M, 16, nullptr, 16, nullptr, &petscWeights); +#endif +#endif } void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z, @@ -63,6 +152,11 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z const int ny = localmesh->LocalNy; const int nz = localmesh->LocalNz; + const int xend = (localmesh->xend - localmesh->xstart + 1) * localmesh->getNXPE() + + localmesh->xstart - 1; +#ifdef HS_USE_PETSC + IndConverter conv{localmesh}; +#endif BOUT_FOR(i, getRegion(region)) { const int x = i.x(); const int y = i.y(); @@ -79,8 +173,8 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z BoutReal t_z = delta_z(x, y, z) - static_cast(k_corner(x, y, z)); // NOTE: A (small) hack to avoid one-sided differences - if (i_corn >= localmesh->xend) { - i_corn = localmesh->xend - 1; + if (i_corn >= xend) { + i_corn = xend - 1; t_x = 1.0; } if (i_corn < localmesh->xstart) { @@ -116,7 +210,6 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z h11_x[i] = (t_x * t_x * t_x) - (t_x * t_x); h11_z[i] = (t_z * t_z * t_z) - (t_z * t_z); -#define USE_NEW_WEIGHTS 1 #if USE_NEW_WEIGHTS for (int w =0; w<16;++w){ @@ -175,8 +268,28 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z newWeights[13][i] -= h11_x[i] * h11_z[i] / 4; newWeights[7][i] -= h11_x[i] * h11_z[i] / 4; newWeights[5][i] += h11_x[i] * h11_z[i] / 4; +#ifdef HS_USE_PETSC + PetscInt idxn[1] = {/* ; idxn[0] = */ conv.fromLocalToGlobal(x, y + y_offset, z)}; + // ixstep = mesh->LocalNx * mesh->LocalNz; + for (int j = 0; j < 4; ++j) { + PetscInt idxm[4]; + PetscScalar vals[4]; + for (int k = 0; k < 4; ++k) { + idxm[k] = conv.fromMeshToGlobal(i_corn - 1 + j, y + y_offset, + k_corner(x, y, z) - 1 + k); + vals[k] = newWeights[j * 4 + k][i]; + } + MatSetValues(petscWeights, 1, idxn, 4, idxm, vals, INSERT_VALUES); + } +#endif #endif } +#ifdef HS_USE_PETSC + isInit = true; + MatAssemblyBegin(petscWeights, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(petscWeights, MAT_FINAL_ASSEMBLY); + MatCreateVecs(petscWeights, &rhs, &result); +#endif } void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z, @@ -220,8 +333,22 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region ASSERT1(f.getMesh() == localmesh); Field3D f_interp{emptyFrom(f)}; - #if USE_NEW_WEIGHTS +#ifdef HS_USE_PETSC + BoutReal* ptr; + const BoutReal* cptr; + VecGetArray(rhs, &ptr); + BOUT_FOR(i, f.getRegion("RGN_NOY")) { ptr[int(i)] = f[i]; } + VecRestoreArray(rhs, &ptr); + MatMult(petscWeights, rhs, result); + VecGetArrayRead(result, &cptr); + const auto region2 = fmt::format("RGN_YPAR_{:+d}", y_offset); + BOUT_FOR(i, f.getRegion(region2)) { + f_interp[i] = cptr[int(i)]; + ASSERT2(std::isfinite(cptr[int(i)])); + } + VecRestoreArrayRead(result, &cptr); +#else BOUT_FOR(i, getRegion(region)) { auto ic = i_corner[i]; auto iyp = i.yp(y_offset); @@ -234,6 +361,7 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region f_interp[iyp] += newWeights[w*4+3][i] * f[ic.zp(2).xp(w-1)]; } } +#endif return f_interp; #else // Derivatives are used for tension and need to be on dimensionless From 8234ccf92a8725f05453b22e874ec989c4e52121 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 2 Feb 2023 10:56:08 +0100 Subject: [PATCH 0009/1827] Test parallised interpolation if PETSc is found --- tests/MMS/spatial/fci/data/BOUT.inp | 2 +- tests/MMS/spatial/fci/runtest | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp index 5f377bbb56..b845e22012 100644 --- a/tests/MMS/spatial/fci/data/BOUT.inp +++ b/tests/MMS/spatial/fci/data/BOUT.inp @@ -20,4 +20,4 @@ y_periodic = true z_periodic = true [mesh:paralleltransform:xzinterpolation] -type = lagrange4pt +type = hermitespline diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 1613155ed2..afff928087 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -19,13 +19,13 @@ from sys import stdout import zoidberg as zb -nx = 3 # Not changed for these tests +nx = 4 # Not changed for these tests # Resolution in y and z nlist = [8, 16, 32, 64, 128] # Number of parallel slices (in each direction) -nslices = [1, 2] +nslices = [1] # , 2] directory = "data" @@ -59,7 +59,7 @@ for nslice in nslices: # Note that the Bz and Bzprime parameters here must be the same as in mms.py field = zb.field.Slab(Bz=0.05, Bzprime=0.1) # Create rectangular poloidal grids - poloidal_grid = zb.poloidal_grid.RectangularPoloidalGrid(nx, n, 0.1, 1.0) + poloidal_grid = zb.poloidal_grid.RectangularPoloidalGrid(nx, n, 0.1, 1.0, MXG=1) # Set the ylength and y locations ylength = 10.0 @@ -72,12 +72,12 @@ for nslice in nslices: # Create the grid grid = zb.grid.Grid(poloidal_grid, ycoords, ylength, yperiodic=yperiodic) # Make and write maps - maps = zb.make_maps(grid, field, nslice=nslice, quiet=True) + maps = zb.make_maps(grid, field, nslice=nslice, quiet=True, MXG=1) zb.write_maps( grid, field, maps, new_names=False, metric2d=conf.isMetric2D(), quiet=True ) - args = " MZ={} MYG={} mesh:paralleltransform:y_periodic={} mesh:ddy:first={}".format( + args = " MZ={} MYG={} mesh:paralleltransform:y_periodic={} mesh:ddy:first={} NXPE=2".format( n, nslice, yperiodic, method_orders[nslice]["name"] ) From 9525e2cb7ccdbfe4d888e9f0e0f019e29936d1de Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 2 Feb 2023 12:43:13 +0100 Subject: [PATCH 0010/1827] Fall back to region if not shifted --- src/mesh/interpolation/hermite_spline_xz.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index 1b63c84231..37dc48662f 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -342,7 +342,7 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region VecRestoreArray(rhs, &ptr); MatMult(petscWeights, rhs, result); VecGetArrayRead(result, &cptr); - const auto region2 = fmt::format("RGN_YPAR_{:+d}", y_offset); + const auto region2 = y_offset == 0 ? region : fmt::format("RGN_YPAR_{:+d}", y_offset); BOUT_FOR(i, f.getRegion(region2)) { f_interp[i] = cptr[int(i)]; ASSERT2(std::isfinite(cptr[int(i)])); From b62082c8cf680deb391fbeaab2184a56d178177b Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 2 Feb 2023 14:05:56 +0100 Subject: [PATCH 0011/1827] Split in X only if we have PETSc --- tests/MMS/spatial/fci/runtest | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index afff928087..7a54b87ccf 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -77,8 +77,12 @@ for nslice in nslices: grid, field, maps, new_names=False, metric2d=conf.isMetric2D(), quiet=True ) - args = " MZ={} MYG={} mesh:paralleltransform:y_periodic={} mesh:ddy:first={} NXPE=2".format( - n, nslice, yperiodic, method_orders[nslice]["name"] + args = " MZ={} MYG={} mesh:paralleltransform:y_periodic={} mesh:ddy:first={} NXPE={}".format( + n, + nslice, + yperiodic, + method_orders[nslice]["name"], + 2 if conf.has["petsc"] else 1, ) # Command to run From e147bc5fa3151726a1b9281993c40ba19b6cfbdc Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 2 Feb 2023 18:32:55 +0100 Subject: [PATCH 0012/1827] Add test-interpolate for splitting in X --- .../integrated/test-interpolate/data/BOUT.inp | 1 - tests/integrated/test-interpolate/runtest | 22 ++++++++++++++----- .../test-interpolate/test_interpolate.cxx | 3 ++- 3 files changed, 18 insertions(+), 8 deletions(-) diff --git a/tests/integrated/test-interpolate/data/BOUT.inp b/tests/integrated/test-interpolate/data/BOUT.inp index 101c63f3c7..804e780bbe 100644 --- a/tests/integrated/test-interpolate/data/BOUT.inp +++ b/tests/integrated/test-interpolate/data/BOUT.inp @@ -4,7 +4,6 @@ # MZ = 4 # Z size -NXPE = 1 ZMAX = 1 MXG = 2 diff --git a/tests/integrated/test-interpolate/runtest b/tests/integrated/test-interpolate/runtest index 08975cfd33..3e8f2c5bc6 100755 --- a/tests/integrated/test-interpolate/runtest +++ b/tests/integrated/test-interpolate/runtest @@ -6,6 +6,7 @@ from boututils.run_wrapper import build_and_log, shell, launch_safe from boutdata import collect +import boutconfig from numpy import sqrt, max, abs, mean, array, log, polyfit from sys import stdout, exit @@ -16,7 +17,7 @@ show_plot = False nxlist = [16, 32, 64, 128] # Only testing 2D (x, z) slices, so only need one processor -nproc = 1 +nproc = 2 # Variables to compare varlist = ["a", "b", "c"] @@ -48,11 +49,9 @@ for method in methods: for nx in nxlist: dx = 1.0 / (nx) - args = ( - " mesh:nx={nx4} mesh:dx={dx} MZ={nx} xzinterpolation:type={method}".format( - nx4=nx + 4, dx=dx, nx=nx, method=method - ) - ) + args = f" mesh:nx={nx + 4} mesh:dx={dx} MZ={nx} xzinterpolation:type={method}" + NXPE = 2 if method == "hermitespline" and boutconfig.has["petsc"] else 1 + args += f" NXPE={NXPE}" cmd = "./test_interpolate" + args @@ -71,6 +70,17 @@ for method in methods: E = interp - solution + if False: + import matplotlib.pyplot as plt + + def myplot(f, lbl=None): + plt.plot(f[:, 0, 6], label=lbl) + + myplot(interp, "interp") + myplot(solution, "sol") + plt.legend() + plt.show() + l2 = float(sqrt(mean(E**2))) linf = float(max(abs(E))) diff --git a/tests/integrated/test-interpolate/test_interpolate.cxx b/tests/integrated/test-interpolate/test_interpolate.cxx index 958409bbc1..8f19cd2a28 100644 --- a/tests/integrated/test-interpolate/test_interpolate.cxx +++ b/tests/integrated/test-interpolate/test_interpolate.cxx @@ -72,7 +72,7 @@ int main(int argc, char **argv) { BoutReal dz = index.z() + dice(); // For the last point, put the displacement inwards // Otherwise we try to interpolate in the guard cells, which doesn't work so well - if (index.x() >= mesh->xend) { + if (index.x() >= mesh->xend && mesh->getNXPE() - 1 == mesh->getXProcIndex()) { dx = index.x() - dice(); } deltax[index] = dx; @@ -87,6 +87,7 @@ int main(int argc, char **argv) { c_solution[index] = c_gen->generate(pos); } + deltax += (mesh->LocalNx - mesh->xstart * 2) * mesh->getXProcIndex(); // Create the interpolation object from the input options auto interp = XZInterpolationFactory::getInstance().create(); From 62d5bbec211127ffa6aeb6c06a8abaeee237b85d Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 2 Feb 2023 18:41:27 +0100 Subject: [PATCH 0013/1827] Cleanup --- src/mesh/interpolation/hermite_spline_xz.cxx | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index 37dc48662f..e41dfa4d03 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -269,7 +269,10 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z newWeights[7][i] -= h11_x[i] * h11_z[i] / 4; newWeights[5][i] += h11_x[i] * h11_z[i] / 4; #ifdef HS_USE_PETSC - PetscInt idxn[1] = {/* ; idxn[0] = */ conv.fromLocalToGlobal(x, y + y_offset, z)}; + PetscInt idxn[1] = {conv.fromLocalToGlobal(x, y + y_offset, z)}; + // output.write("debug: {:d} -> {:d}: {:d}:{:d} -> {:d}:{:d}\n", conv.fromLocalToGlobal(x, y + y_offset, z), + // conv.fromMeshToGlobal(i_corn, y + y_offset, k_corner(x, y, z)), + // x, z, i_corn, k_corner(x, y, z)); // ixstep = mesh->LocalNx * mesh->LocalNz; for (int j = 0; j < 4; ++j) { PetscInt idxm[4]; From 5c324157ca1d2f3eb241eac68e5f55feaef889ae Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 2 Feb 2023 19:44:17 +0100 Subject: [PATCH 0014/1827] Only run in parallel if we split in X --- tests/integrated/test-interpolate/runtest | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/tests/integrated/test-interpolate/runtest b/tests/integrated/test-interpolate/runtest index 3e8f2c5bc6..f5460aff2a 100755 --- a/tests/integrated/test-interpolate/runtest +++ b/tests/integrated/test-interpolate/runtest @@ -16,9 +16,6 @@ show_plot = False # List of NX values to use nxlist = [16, 32, 64, 128] -# Only testing 2D (x, z) slices, so only need one processor -nproc = 2 - # Variables to compare varlist = ["a", "b", "c"] markers = ["bo", "r^", "kx"] @@ -50,8 +47,8 @@ for method in methods: dx = 1.0 / (nx) args = f" mesh:nx={nx + 4} mesh:dx={dx} MZ={nx} xzinterpolation:type={method}" - NXPE = 2 if method == "hermitespline" and boutconfig.has["petsc"] else 1 - args += f" NXPE={NXPE}" + nproc = 2 if method == "hermitespline" and boutconfig.has["petsc"] else 1 + args += f" NXPE={nproc}" cmd = "./test_interpolate" + args From a4a28c6e3c4025500f7a8bbbc0b0e474720703e0 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 3 Feb 2023 13:10:21 +0100 Subject: [PATCH 0015/1827] Delete object release leaks the pointer, reset free's the object. --- tests/integrated/test-interpolate/test_interpolate.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/integrated/test-interpolate/test_interpolate.cxx b/tests/integrated/test-interpolate/test_interpolate.cxx index 8f19cd2a28..ed8e80f43f 100644 --- a/tests/integrated/test-interpolate/test_interpolate.cxx +++ b/tests/integrated/test-interpolate/test_interpolate.cxx @@ -113,6 +113,7 @@ int main(int argc, char **argv) { bout::writeDefaultOutputFile(dump); bout::checkForUnusedOptions(); + interp.reset(); BoutFinalise(); return 0; From 0e28dc1c5dbaa8c5e73bb09020ac2d01924aa850 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 3 Feb 2023 14:05:38 +0100 Subject: [PATCH 0016/1827] Create PetscVecs only once --- src/mesh/interpolation/hermite_spline_xz.cxx | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index e41dfa4d03..4e6fc8320c 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -288,10 +288,12 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z #endif } #ifdef HS_USE_PETSC - isInit = true; MatAssemblyBegin(petscWeights, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(petscWeights, MAT_FINAL_ASSEMBLY); - MatCreateVecs(petscWeights, &rhs, &result); + if (!isInit) { + MatCreateVecs(petscWeights, &rhs, &result); + } + isInit = true; #endif } From 66bde042f6553e372465eb47be0ede532a815bde Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 3 Feb 2023 14:06:09 +0100 Subject: [PATCH 0017/1827] Be more general about cleaning up before BoutFinialise --- tests/integrated/test-interpolate/test_interpolate.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/integrated/test-interpolate/test_interpolate.cxx b/tests/integrated/test-interpolate/test_interpolate.cxx index ed8e80f43f..517d9c2445 100644 --- a/tests/integrated/test-interpolate/test_interpolate.cxx +++ b/tests/integrated/test-interpolate/test_interpolate.cxx @@ -30,7 +30,7 @@ std::shared_ptr getGeneratorFromOptions(const std::string& varna int main(int argc, char **argv) { BoutInitialise(argc, argv); - + { // Random number generator std::default_random_engine generator; // Uniform distribution of BoutReals from 0 to 1 @@ -113,7 +113,7 @@ int main(int argc, char **argv) { bout::writeDefaultOutputFile(dump); bout::checkForUnusedOptions(); - interp.reset(); + } BoutFinalise(); return 0; From adba8774b274073e07eb13cd8165a7594dd99bd0 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 3 Feb 2023 14:56:13 +0100 Subject: [PATCH 0018/1827] Run different interpolations in the fci test --- tests/MMS/spatial/fci/runtest | 197 +++++++++++++++++++--------------- 1 file changed, 108 insertions(+), 89 deletions(-) diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 7a54b87ccf..d68b6c6ca1 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -47,96 +47,115 @@ failures = [] build_and_log("FCI MMS test") for nslice in nslices: - error_2[nslice] = [] - error_inf[nslice] = [] - - # Which central difference scheme to use and its expected order - order = nslice * 2 - method_orders[nslice] = {"name": "C{}".format(order), "order": order} - - for n in nlist: - # Define the magnetic field using new poloidal gridding method - # Note that the Bz and Bzprime parameters here must be the same as in mms.py - field = zb.field.Slab(Bz=0.05, Bzprime=0.1) - # Create rectangular poloidal grids - poloidal_grid = zb.poloidal_grid.RectangularPoloidalGrid(nx, n, 0.1, 1.0, MXG=1) - # Set the ylength and y locations - ylength = 10.0 - - if yperiodic: - ycoords = linspace(0.0, ylength, n, endpoint=False) + for method in [ + "hermitespline", + "lagrange4pt", + "bilinear", + # "monotonichermitespline", + ]: + error_2[nslice] = [] + error_inf[nslice] = [] + + # Which central difference scheme to use and its expected order + order = nslice * 2 + method_orders[nslice] = {"name": "C{}".format(order), "order": order} + + for n in nlist: + # Define the magnetic field using new poloidal gridding method + # Note that the Bz and Bzprime parameters here must be the same as in mms.py + field = zb.field.Slab(Bz=0.05, Bzprime=0.1) + # Create rectangular poloidal grids + poloidal_grid = zb.poloidal_grid.RectangularPoloidalGrid( + nx, n, 0.1, 1.0, MXG=1 + ) + # Set the ylength and y locations + ylength = 10.0 + + if yperiodic: + ycoords = linspace(0.0, ylength, n, endpoint=False) + else: + # Doesn't include the end points + ycoords = (arange(n) + 0.5) * ylength / float(n) + + # Create the grid + grid = zb.grid.Grid(poloidal_grid, ycoords, ylength, yperiodic=yperiodic) + # Make and write maps + maps = zb.make_maps(grid, field, nslice=nslice, quiet=True, MXG=1) + zb.write_maps( + grid, + field, + maps, + new_names=False, + metric2d=conf.isMetric2D(), + quiet=True, + ) + + args = " MZ={} MYG={} mesh:paralleltransform:y_periodic={} mesh:ddy:first={} NXPE={}".format( + n, + nslice, + yperiodic, + method_orders[nslice]["name"], + 2 if conf.has["petsc"] and method == "hermitespline" else 1, + ) + args += f" mesh:paralleltransform:xzinterpolation:type={method}" + + # Command to run + cmd = "./fci_mms " + args + + print("Running command: " + cmd) + + # Launch using MPI + s, out = launch_safe(cmd, nproc=nproc, mthread=mthread, pipe=True) + + # Save output to log file + with open("run.log." + str(n), "w") as f: + f.write(out) + + if s: + print("Run failed!\nOutput was:\n") + print(out) + exit(s) + + # Collect data + l_2 = collect( + "l_2", + tind=[1, 1], + info=False, + path=directory, + xguards=False, + yguards=False, + ) + l_inf = collect( + "l_inf", + tind=[1, 1], + info=False, + path=directory, + xguards=False, + yguards=False, + ) + + error_2[nslice].append(l_2) + error_inf[nslice].append(l_inf) + + print("Errors : l-2 {:f} l-inf {:f}".format(l_2, l_inf)) + + dx = 1.0 / array(nlist) + + # Calculate convergence order + fit = polyfit(log(dx), log(error_2[nslice]), 1) + order = fit[0] + stdout.write("Convergence order = {:f} (fit)".format(order)) + + order = log(error_2[nslice][-2] / error_2[nslice][-1]) / log(dx[-2] / dx[-1]) + stdout.write(", {:f} (small spacing)".format(order)) + + # Should be close to the expected order + if order > method_orders[nslice]["order"] * 0.95: + print("............ PASS\n") else: - # Doesn't include the end points - ycoords = (arange(n) + 0.5) * ylength / float(n) - - # Create the grid - grid = zb.grid.Grid(poloidal_grid, ycoords, ylength, yperiodic=yperiodic) - # Make and write maps - maps = zb.make_maps(grid, field, nslice=nslice, quiet=True, MXG=1) - zb.write_maps( - grid, field, maps, new_names=False, metric2d=conf.isMetric2D(), quiet=True - ) - - args = " MZ={} MYG={} mesh:paralleltransform:y_periodic={} mesh:ddy:first={} NXPE={}".format( - n, - nslice, - yperiodic, - method_orders[nslice]["name"], - 2 if conf.has["petsc"] else 1, - ) - - # Command to run - cmd = "./fci_mms " + args - - print("Running command: " + cmd) - - # Launch using MPI - s, out = launch_safe(cmd, nproc=nproc, mthread=mthread, pipe=True) - - # Save output to log file - with open("run.log." + str(n), "w") as f: - f.write(out) - - if s: - print("Run failed!\nOutput was:\n") - print(out) - exit(s) - - # Collect data - l_2 = collect( - "l_2", tind=[1, 1], info=False, path=directory, xguards=False, yguards=False - ) - l_inf = collect( - "l_inf", - tind=[1, 1], - info=False, - path=directory, - xguards=False, - yguards=False, - ) - - error_2[nslice].append(l_2) - error_inf[nslice].append(l_inf) - - print("Errors : l-2 {:f} l-inf {:f}".format(l_2, l_inf)) - - dx = 1.0 / array(nlist) - - # Calculate convergence order - fit = polyfit(log(dx), log(error_2[nslice]), 1) - order = fit[0] - stdout.write("Convergence order = {:f} (fit)".format(order)) - - order = log(error_2[nslice][-2] / error_2[nslice][-1]) / log(dx[-2] / dx[-1]) - stdout.write(", {:f} (small spacing)".format(order)) - - # Should be close to the expected order - if order > method_orders[nslice]["order"] * 0.95: - print("............ PASS\n") - else: - print("............ FAIL\n") - success = False - failures.append(method_orders[nslice]["name"]) + print("............ FAIL\n") + success = False + failures.append(method_orders[nslice]["name"]) with open("fci_mms.pkl", "wb") as output: From cbaf894b2b77b09111f0b38d17653476b3263c35 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 7 Feb 2023 14:06:13 +0100 Subject: [PATCH 0019/1827] Fix parallel boundary region with x splitting --- src/mesh/parallel/fci.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index 77f34fd282..baf0f3fc6b 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -158,6 +158,7 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& dy, Options& options, const int ncz = map_mesh.LocalNz; BoutMask to_remove(map_mesh); + const int xend = map_mesh.xstart + (map_mesh.xend - map_mesh.xstart + 1) * map_mesh.getNXPE() - 1; // Serial loop because call to BoundaryRegionPar::addPoint // (probably?) can't be done in parallel BOUT_FOR_SERIAL(i, xt_prime.getRegion("RGN_NOBNDRY")) { @@ -171,7 +172,7 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& dy, Options& options, } } - if ((xt_prime[i] >= map_mesh.xstart) and (xt_prime[i] <= map_mesh.xend)) { + if ((xt_prime[i] >= map_mesh.xstart) and (xt_prime[i] <= xend)) { // Not a boundary continue; } From 5673f0c2ea8ecef929b99fb473cc25bc393a76e7 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 7 Feb 2023 14:12:52 +0100 Subject: [PATCH 0020/1827] Add integrated test for FCI X splitting * helped finding the bug for the boundary * rather slow (around 1 minute) * needs internet connectivity --- cmake/BOUT++functions.cmake | 16 +++++- tests/integrated/CMakeLists.txt | 2 + tests/integrated/test-fci-mpi/CMakeLists.txt | 8 +++ tests/integrated/test-fci-mpi/data/BOUT.inp | 28 ++++++++++ tests/integrated/test-fci-mpi/fci_mpi.cxx | 37 +++++++++++++ tests/integrated/test-fci-mpi/runtest | 58 ++++++++++++++++++++ 6 files changed, 148 insertions(+), 1 deletion(-) create mode 100644 tests/integrated/test-fci-mpi/CMakeLists.txt create mode 100644 tests/integrated/test-fci-mpi/data/BOUT.inp create mode 100644 tests/integrated/test-fci-mpi/fci_mpi.cxx create mode 100755 tests/integrated/test-fci-mpi/runtest diff --git a/cmake/BOUT++functions.cmake b/cmake/BOUT++functions.cmake index 40e45f99be..77279dfd4b 100644 --- a/cmake/BOUT++functions.cmake +++ b/cmake/BOUT++functions.cmake @@ -162,7 +162,7 @@ endfunction() # function(bout_add_integrated_or_mms_test BUILD_CHECK_TARGET TESTNAME) set(options USE_RUNTEST USE_DATA_BOUT_INP) - set(oneValueArgs EXECUTABLE_NAME PROCESSORS) + set(oneValueArgs EXECUTABLE_NAME PROCESSORS DOWNLOAD DOWNLOAD_NAME) set(multiValueArgs SOURCES EXTRA_FILES REQUIRES CONFLICTS TESTARGS EXTRA_DEPENDS) cmake_parse_arguments(BOUT_TEST_OPTIONS "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN}) @@ -202,6 +202,20 @@ function(bout_add_integrated_or_mms_test BUILD_CHECK_TARGET TESTNAME) add_custom_target(${TESTNAME}) endif() + if (BOUT_TEST_OPTIONS_DOWNLOAD) + if (NOT BOUT_TEST_OPTIONS_DOWNLOAD_NAME) + message(FATAL_ERROR "We need DOWNLOAD_NAME if we should DOWNLOAD!") + endif() + set(output ) + add_custom_command(OUTPUT ${BOUT_TEST_OPTIONS_DOWNLOAD_NAME} + COMMAND wget ${BOUT_TEST_OPTIONS_DOWNLOAD} -O ${BOUT_TEST_OPTIONS_DOWNLOAD_NAME} + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + COMMENT "Downloading ${BOUT_TEST_OPTIONS_DOWNLOAD_NAME}" + ) + add_custom_target(download_test_data DEPENDS ${BOUT_TEST_OPTIONS_DOWNLOAD_NAME}) + add_dependencies(${TESTNAME} download_test_data) + endif() + if (BOUT_TEST_OPTIONS_EXTRA_DEPENDS) add_dependencies(${TESTNAME} ${BOUT_TEST_OPTIONS_EXTRA_DEPENDS}) endif() diff --git a/tests/integrated/CMakeLists.txt b/tests/integrated/CMakeLists.txt index 2fe72dfe2d..89cbf6ffe6 100644 --- a/tests/integrated/CMakeLists.txt +++ b/tests/integrated/CMakeLists.txt @@ -9,6 +9,8 @@ add_subdirectory(test-cyclic) add_subdirectory(test-delp2) add_subdirectory(test-drift-instability) add_subdirectory(test-drift-instability-staggered) +add_subdirectory(test-fci-boundary) +add_subdirectory(test-fci-mpi) add_subdirectory(test-fieldgroupComm) add_subdirectory(test-griddata) add_subdirectory(test-griddata-yboundary-guards) diff --git a/tests/integrated/test-fci-mpi/CMakeLists.txt b/tests/integrated/test-fci-mpi/CMakeLists.txt new file mode 100644 index 0000000000..6a1ec33ac6 --- /dev/null +++ b/tests/integrated/test-fci-mpi/CMakeLists.txt @@ -0,0 +1,8 @@ +bout_add_mms_test(test-fci-mpi + SOURCES fci_mpi.cxx + USE_RUNTEST + USE_DATA_BOUT_INP + PROCESSORS 6 + DOWNLOAD https://zenodo.org/record/7614499/files/W7X-conf4-36x8x128.fci.nc?download=1 + DOWNLOAD_NAME grid.fci.nc +) diff --git a/tests/integrated/test-fci-mpi/data/BOUT.inp b/tests/integrated/test-fci-mpi/data/BOUT.inp new file mode 100644 index 0000000000..47272dab61 --- /dev/null +++ b/tests/integrated/test-fci-mpi/data/BOUT.inp @@ -0,0 +1,28 @@ +grid = grid.fci.nc + +[mesh] +symmetricglobalx = true + +[mesh:ddy] +first = C2 +second = C2 + +[mesh:paralleltransform] +type = fci +y_periodic = true +z_periodic = true + +[mesh:paralleltransform:xzinterpolation] +type = hermitespline + +[input_0] +function = sin(z) + +[input_1] +function = cos(y) + +[input_2] +function = sin(x) + +[input_3] +function = sin(x) * sin(z) * cos(y) diff --git a/tests/integrated/test-fci-mpi/fci_mpi.cxx b/tests/integrated/test-fci-mpi/fci_mpi.cxx new file mode 100644 index 0000000000..b353493dda --- /dev/null +++ b/tests/integrated/test-fci-mpi/fci_mpi.cxx @@ -0,0 +1,37 @@ +#include "bout.hxx" +#include "derivs.hxx" +#include "field_factory.hxx" + +int main(int argc, char** argv) { + BoutInitialise(argc, argv); + { + using bout::globals::mesh; + Options *options = Options::getRoot(); + int i=0; + std::string default_str {"not_set"}; + Options dump; + while (true) { + std::string temp_str; + options->get(fmt::format("input_{:d}:function", i), temp_str, default_str); + if (temp_str == default_str) { + break; + } + Field3D input{FieldFactory::get()->create3D(fmt::format("input_{:d}:function", i), Options::getRoot(), mesh)}; + //options->get(fmt::format("input_{:d}:boundary_perp", i), temp_str, s"free_o3"); + mesh->communicate(input); + input.applyParallelBoundary("parallel_neumann_o2"); + for (int slice = -mesh->ystart; slice <= mesh->ystart; ++slice) { + if (slice) { + Field3D tmp{0.}; + BOUT_FOR(i, tmp.getRegion("RGN_NOBNDRY")) { + tmp[i] = input.ynext(slice)[i.yp(slice)]; + } + dump[fmt::format("output_{:d}_{:+d}", i, slice)] = tmp; + } + } + ++i; + } + bout::writeDefaultOutputFile(dump); + } + BoutFinalise(); +} diff --git a/tests/integrated/test-fci-mpi/runtest b/tests/integrated/test-fci-mpi/runtest new file mode 100755 index 0000000000..4ac0e43460 --- /dev/null +++ b/tests/integrated/test-fci-mpi/runtest @@ -0,0 +1,58 @@ +#!/usr/bin/env python3 +# +# Python script to run and analyse MMS test +# + +# Cores: 8 +# requires: metric_3d + +from boututils.run_wrapper import build_and_log, launch_safe, shell_safe +from boutdata.collect import collect +import boutconfig as conf +import itertools + +import numpy as np + +# Resolution in x and y +nlist = [1, 2, 4] + +maxcores = 8 + +nslices = [1] + +success = True + +build_and_log("FCI MMS test") + +for nslice in nslices: + for NXPE, NYPE in itertools.product(nlist, nlist): + + if NXPE * NYPE > maxcores: + continue + + args = f"NXPE={NXPE} NYPE={NYPE}" + # Command to run + cmd = f"./fci_mpi {args}" + + print(f"Running command: {cmd}") + + mthread = maxcores // (NXPE * NYPE) + # Launch using MPI + _, out = launch_safe(cmd, nproc=NXPE * NYPE, mthread=mthread, pipe=True) + + # Save output to log file + with open("run.log.{NXPE}.{NYPE}.{nslice}.log", "w") as f: + f.write(out) + + collect_kw = dict(info=False, xguards=False, yguards=False, path="data") + if NXPE == NYPE == 1: + # reference data! + ref = {} + for i in range(4): + for yp in range(1, nslice + 1): + for y in [-yp, yp]: + name = f"output_{i}_{y:+d}" + ref[name] = collect(name, **collect_kw) + else: + for name, val in ref.items(): + assert np.allclose(val, collect(name, **collect_kw)) From 9149bf406fe01a98966095455757f71c0305272c Mon Sep 17 00:00:00 2001 From: dschwoerer Date: Wed, 8 Feb 2023 09:20:18 +0000 Subject: [PATCH 0021/1827] Apply black changes --- bin/bout-v5-xzinterpolation-upgrader.py | 3 --- bin/bout_3to4.py | 1 - src/field/gen_fieldops.py | 2 -- tests/MMS/GBS/circle.py | 1 - tests/MMS/GBS/mms-slab3d.py | 1 - tests/MMS/GBS/runtest-slab3d | 1 - tests/integrated/test-drift-instability/runtest | 1 - tests/integrated/test-fci-mpi/runtest | 1 - tests/integrated/test-multigrid_laplace/runtest | 1 - .../test-multigrid_laplace/runtest_multiple_grids | 1 - tests/integrated/test-multigrid_laplace/runtest_unsheared | 1 - tests/integrated/test-naulin-laplace/runtest | 1 - .../integrated/test-naulin-laplace/runtest_multiple_grids | 1 - tests/integrated/test-naulin-laplace/runtest_unsheared | 1 - tests/integrated/test-petsc_laplace_MAST-grid/runtest | 1 - tests/integrated/test-twistshift-staggered/runtest | 1 + tests/integrated/test-twistshift/runtest | 2 ++ tests/integrated/test-yupdown-weights/runtest | 1 - tests/integrated/test-yupdown/runtest | 1 - tests/integrated/test_suite | 1 + tools/pylib/post_bout/__init__.py | 1 - tools/pylib/post_bout/basic_info.py | 3 --- tools/pylib/post_bout/grate2.py | 2 -- tools/pylib/post_bout/pb_corral.py | 4 ---- tools/pylib/post_bout/pb_draw.py | 7 ------- tools/pylib/post_bout/pb_nonlinear.py | 1 - tools/pylib/post_bout/pb_present.py | 1 - tools/pylib/post_bout/read_cxx.py | 1 - tools/pylib/post_bout/read_inp.py | 3 --- tools/pylib/post_bout/rms.py | 1 - tools/tokamak_grids/elite/elite2nc | 1 + tools/tokamak_grids/gato/gato2nc | 1 + 32 files changed, 6 insertions(+), 44 deletions(-) diff --git a/bin/bout-v5-xzinterpolation-upgrader.py b/bin/bout-v5-xzinterpolation-upgrader.py index 1e19c6a034..37c79e0de8 100755 --- a/bin/bout-v5-xzinterpolation-upgrader.py +++ b/bin/bout-v5-xzinterpolation-upgrader.py @@ -64,7 +64,6 @@ def fix_header_includes(old_header, new_header, source): def fix_interpolations(old_interpolation, new_interpolation, source): - return re.sub( r""" \b{}\b @@ -118,7 +117,6 @@ def clang_fix_interpolation(old_interpolation, new_interpolation, node, source): def fix_factories(old_factory, new_factory, source): - return re.sub( r""" \b{}\b @@ -186,7 +184,6 @@ def apply_fixes(headers, interpolations, factories, source): def clang_apply_fixes(headers, interpolations, factories, filename, source): - # translation unit tu = clang_parse(filename, source) diff --git a/bin/bout_3to4.py b/bin/bout_3to4.py index 02481fcf6e..d41db36fbb 100755 --- a/bin/bout_3to4.py +++ b/bin/bout_3to4.py @@ -195,7 +195,6 @@ def throw_warnings(line_text, filename, line_num): if __name__ == "__main__": - epilog = """ Currently bout_3to4 can detect the following transformations are needed: - Triple square brackets instead of round brackets for subscripts diff --git a/src/field/gen_fieldops.py b/src/field/gen_fieldops.py index 646580559f..68a3a1b059 100755 --- a/src/field/gen_fieldops.py +++ b/src/field/gen_fieldops.py @@ -189,7 +189,6 @@ def returnType(f1, f2): if __name__ == "__main__": - parser = argparse.ArgumentParser( description="Generate code for the Field arithmetic operators" ) @@ -274,7 +273,6 @@ def returnType(f1, f2): rhs.name = "rhs" for operator, operator_name in operators.items(): - template_args = { "operator": operator, "operator_name": operator_name, diff --git a/tests/MMS/GBS/circle.py b/tests/MMS/GBS/circle.py index b2553f47c3..b33cc77d8c 100644 --- a/tests/MMS/GBS/circle.py +++ b/tests/MMS/GBS/circle.py @@ -30,7 +30,6 @@ def generate( mxg=2, file="circle.nc", ): - # q = rBt / RBp Bp = r * Bt / (R * q) diff --git a/tests/MMS/GBS/mms-slab3d.py b/tests/MMS/GBS/mms-slab3d.py index f88a897c98..2f958cfa69 100644 --- a/tests/MMS/GBS/mms-slab3d.py +++ b/tests/MMS/GBS/mms-slab3d.py @@ -277,7 +277,6 @@ def C(f): # print "\n\nDelp2 phi = ", Delp2(phi, metric).subs(replace) if not estatic: - Spsi = Delp2(psi, metric) - 0.5 * Ne * mi_me * beta_e * psi - Ne * (Vi - VePsi) print("\n[psi]") print("\nsolution = " + exprToStr(psi.subs(replace))) diff --git a/tests/MMS/GBS/runtest-slab3d b/tests/MMS/GBS/runtest-slab3d index cdcde3e1b3..f193583e47 100755 --- a/tests/MMS/GBS/runtest-slab3d +++ b/tests/MMS/GBS/runtest-slab3d @@ -95,7 +95,6 @@ with open("mms-slab3d.pkl", "wb") as output: # plot errors for var, mark in zip(varlist, markers): - order = log(error_2[var][-1] / error_2[var][-2]) / log(dx[-1] / dx[-2]) print("%s Convergence order = %f" % (var, order)) if 1.9 < order < 2.1: diff --git a/tests/integrated/test-drift-instability/runtest b/tests/integrated/test-drift-instability/runtest index 16b35d69d9..e8ddddc11d 100755 --- a/tests/integrated/test-drift-instability/runtest +++ b/tests/integrated/test-drift-instability/runtest @@ -197,7 +197,6 @@ def run_zeff_case(zeff): if __name__ == "__main__": - parser = argparse.ArgumentParser("Run drift-instability test") parser.add_argument( "-Z", diff --git a/tests/integrated/test-fci-mpi/runtest b/tests/integrated/test-fci-mpi/runtest index 4ac0e43460..e12b330326 100755 --- a/tests/integrated/test-fci-mpi/runtest +++ b/tests/integrated/test-fci-mpi/runtest @@ -26,7 +26,6 @@ build_and_log("FCI MMS test") for nslice in nslices: for NXPE, NYPE in itertools.product(nlist, nlist): - if NXPE * NYPE > maxcores: continue diff --git a/tests/integrated/test-multigrid_laplace/runtest b/tests/integrated/test-multigrid_laplace/runtest index 76914d19e4..4a7455f80b 100755 --- a/tests/integrated/test-multigrid_laplace/runtest +++ b/tests/integrated/test-multigrid_laplace/runtest @@ -29,7 +29,6 @@ print("Running multigrid Laplacian inversion test") success = True for nproc in [1, 3]: - # Make sure we don't use too many cores: # Reduce number of OpenMP threads when using multiple MPI processes mthread = 2 diff --git a/tests/integrated/test-multigrid_laplace/runtest_multiple_grids b/tests/integrated/test-multigrid_laplace/runtest_multiple_grids index 26cdd96ff6..6817120b13 100755 --- a/tests/integrated/test-multigrid_laplace/runtest_multiple_grids +++ b/tests/integrated/test-multigrid_laplace/runtest_multiple_grids @@ -26,7 +26,6 @@ success = True for nproc in [1, 2, 4]: for inputfile in ["BOUT_jy4.inp", "BOUT_jy63.inp", "BOUT_jy127.inp"]: - # set nxpe on the command line as we only use solution from one point in y, so splitting in y-direction is redundant (and also doesn't help test the multigrid solver) cmd = "./test_multigrid_laplace -f " + inputfile + " NXPE=" + str(nproc) diff --git a/tests/integrated/test-multigrid_laplace/runtest_unsheared b/tests/integrated/test-multigrid_laplace/runtest_unsheared index fd0e11fbbf..cda68f2167 100755 --- a/tests/integrated/test-multigrid_laplace/runtest_unsheared +++ b/tests/integrated/test-multigrid_laplace/runtest_unsheared @@ -25,7 +25,6 @@ print("Running multigrid Laplacian inversion test") success = True for nproc in [1, 3]: - # Make sure we don't use too many cores: # Reduce number of OpenMP threads when using multiple MPI processes mthread = 2 diff --git a/tests/integrated/test-naulin-laplace/runtest b/tests/integrated/test-naulin-laplace/runtest index 82dd22776e..f972eab6cc 100755 --- a/tests/integrated/test-naulin-laplace/runtest +++ b/tests/integrated/test-naulin-laplace/runtest @@ -29,7 +29,6 @@ print("Running LaplaceNaulin inversion test") success = True for nproc in [1, 3]: - # Make sure we don't use too many cores: # Reduce number of OpenMP threads when using multiple MPI processes mthread = 2 diff --git a/tests/integrated/test-naulin-laplace/runtest_multiple_grids b/tests/integrated/test-naulin-laplace/runtest_multiple_grids index 2d583fd1b8..c0281c3a4e 100755 --- a/tests/integrated/test-naulin-laplace/runtest_multiple_grids +++ b/tests/integrated/test-naulin-laplace/runtest_multiple_grids @@ -26,7 +26,6 @@ success = True for nproc in [1, 2, 4]: for inputfile in ["BOUT_jy4.inp", "BOUT_jy63.inp", "BOUT_jy127.inp"]: - # set nxpe on the command line as we only use solution from one point in y, so splitting in y-direction is redundant (and also doesn't help test the solver) cmd = "./test_naulin_laplace -f " + inputfile + " NXPE=" + str(nproc) diff --git a/tests/integrated/test-naulin-laplace/runtest_unsheared b/tests/integrated/test-naulin-laplace/runtest_unsheared index ec956686ef..8f47f33026 100755 --- a/tests/integrated/test-naulin-laplace/runtest_unsheared +++ b/tests/integrated/test-naulin-laplace/runtest_unsheared @@ -25,7 +25,6 @@ print("Running LaplaceNaulin inversion test") success = True for nproc in [1, 3]: - # Make sure we don't use too many cores: # Reduce number of OpenMP threads when using multiple MPI processes mthread = 2 diff --git a/tests/integrated/test-petsc_laplace_MAST-grid/runtest b/tests/integrated/test-petsc_laplace_MAST-grid/runtest index c1c973e3b2..bf6656fbf2 100755 --- a/tests/integrated/test-petsc_laplace_MAST-grid/runtest +++ b/tests/integrated/test-petsc_laplace_MAST-grid/runtest @@ -43,7 +43,6 @@ for nproc in [1, 2, 4]: # if nproc > 2: # nxpe = 2 for jy in [2, 34, 65, 81, 113]: - cmd = ( "./test_petsc_laplace_MAST_grid grid=grids/grid_MAST_SOL_jyis{}.nc".format( jy diff --git a/tests/integrated/test-twistshift-staggered/runtest b/tests/integrated/test-twistshift-staggered/runtest index 2976a80ecd..5bb9963db7 100755 --- a/tests/integrated/test-twistshift-staggered/runtest +++ b/tests/integrated/test-twistshift-staggered/runtest @@ -24,6 +24,7 @@ check = collect("check", path=datapath, yguards=True, info=False) success = True + # Check test_aligned is *not* periodic in y def test1(ylower, yupper): global success diff --git a/tests/integrated/test-twistshift/runtest b/tests/integrated/test-twistshift/runtest index 26b5c2b135..c4858970c4 100755 --- a/tests/integrated/test-twistshift/runtest +++ b/tests/integrated/test-twistshift/runtest @@ -24,6 +24,7 @@ result = collect("result", path=datapath, yguards=True, info=False) success = True + # Check test_aligned is *not* periodic in y def test1(ylower, yupper): global success @@ -49,6 +50,7 @@ if numpy.any(numpy.abs(result - test) > tol): print("Fail - result has not been communicated correctly - is different from input") success = False + # Check result is periodic in y def test2(ylower, yupper): global success diff --git a/tests/integrated/test-yupdown-weights/runtest b/tests/integrated/test-yupdown-weights/runtest index e40b81368f..4b59d08cb0 100755 --- a/tests/integrated/test-yupdown-weights/runtest +++ b/tests/integrated/test-yupdown-weights/runtest @@ -10,7 +10,6 @@ build_and_log("parallel slices and weights test") failed = False for shifttype in ["shiftedinterp"]: - s, out = launch_safe( "./test_yupdown_weights mesh:paralleltransform:type=" + shifttype, nproc=1, diff --git a/tests/integrated/test-yupdown/runtest b/tests/integrated/test-yupdown/runtest index 1b24f86327..34fcd36496 100755 --- a/tests/integrated/test-yupdown/runtest +++ b/tests/integrated/test-yupdown/runtest @@ -12,7 +12,6 @@ build_and_log("parallel slices test") failed = False for shifttype in ["shifted", "shiftedinterp"]: - s, out = launch_safe( "./test_yupdown mesh:paralleltransform:type=" + shifttype, nproc=1, diff --git a/tests/integrated/test_suite b/tests/integrated/test_suite index e6012e3869..307a8d84b3 100755 --- a/tests/integrated/test_suite +++ b/tests/integrated/test_suite @@ -241,6 +241,7 @@ if args.get_list: print(test) sys.exit(0) + # A function to get more threads from the job server def get_threads(): global js_read diff --git a/tools/pylib/post_bout/__init__.py b/tools/pylib/post_bout/__init__.py index a06792ed41..069ae3e85b 100644 --- a/tools/pylib/post_bout/__init__.py +++ b/tools/pylib/post_bout/__init__.py @@ -15,7 +15,6 @@ import os try: - boutpath = os.environ["BOUT_TOP"] pylibpath = boutpath + "/tools/pylib" boutdatapath = pylibpath + "/boutdata" diff --git a/tools/pylib/post_bout/basic_info.py b/tools/pylib/post_bout/basic_info.py index 28660c39bc..563bfe6e98 100644 --- a/tools/pylib/post_bout/basic_info.py +++ b/tools/pylib/post_bout/basic_info.py @@ -10,7 +10,6 @@ def basic_info(data, meta, rescale=True, rotate=False, user_peak=0, nonlinear=None): - print("in basic_info") # from . import read_grid,parse_inp,read_inp,show @@ -227,7 +226,6 @@ def fft_info( jumps = np.where(abs(phase_r) > old_div(np.pi, 32)) # print jumps if len(jumps[0]) != 0: - all_pts = np.array(list(range(0, nt))) good_pts = (np.where(abs(phase_r) < old_div(np.pi, 3)))[0] # print good_pts,good_pts @@ -363,7 +361,6 @@ def fft_info( # return a 2d array fof boolean values, a very simple boolian filter def local_maxima(array2d, user_peak, index=False, count=4, floor=0, bug=False): - from operator import itemgetter, attrgetter if user_peak == 0: diff --git a/tools/pylib/post_bout/grate2.py b/tools/pylib/post_bout/grate2.py index 58f5a47e2b..5157863cc2 100644 --- a/tools/pylib/post_bout/grate2.py +++ b/tools/pylib/post_bout/grate2.py @@ -13,7 +13,6 @@ def avgrate(p, y=None, tind=None): - if tind is None: tind = 0 @@ -25,7 +24,6 @@ def avgrate(p, y=None, tind=None): growth = np.zeros((ni, nj)) with np.errstate(divide="ignore"): - for i in range(ni): for j in range(nj): growth[i, j] = np.gradient(np.log(rmsp_f[tind::, i, j]))[-1] diff --git a/tools/pylib/post_bout/pb_corral.py b/tools/pylib/post_bout/pb_corral.py index 412a677921..df9a9b00b6 100644 --- a/tools/pylib/post_bout/pb_corral.py +++ b/tools/pylib/post_bout/pb_corral.py @@ -40,7 +40,6 @@ def corral( cached=True, refresh=False, debug=False, IConly=1, logname="status.log", skew=False ): - print("in corral") log = read_log(logname=logname) # done = log['done'] @@ -53,7 +52,6 @@ def corral( print("current:", current) if refresh == True: - for i, path in enumerate(runs): print(i, path) a = post_bout.save(path=path, IConly=IConly) # re post-process a run @@ -61,7 +59,6 @@ def corral( elif ( cached == False ): # if all the ind. simulation pkl files are in place skip this part - a = post_bout.save(path=current) # save to current dir # here is really where you shoudl write to status.log # write_log('status.log', @@ -111,7 +108,6 @@ def islist(input): class LinRes(object): def __init__(self, all_modes): - self.mode_db = all_modes self.db = all_modes # self.ave_db = all_ave diff --git a/tools/pylib/post_bout/pb_draw.py b/tools/pylib/post_bout/pb_draw.py index 75bca03a15..272aab9c35 100644 --- a/tools/pylib/post_bout/pb_draw.py +++ b/tools/pylib/post_bout/pb_draw.py @@ -140,7 +140,6 @@ def plottheory( fig1.savefig(pp, format="pdf") plt.close(fig1) else: # if not plot its probably plotted iwth sim data, print chi somewhere - for i, m in enumerate(s.models): textstr = r"$\chi^2$" + "$=%.2f$" % (m.chi[comp].sum()) print(textstr) @@ -173,7 +172,6 @@ def plotomega( trans=False, infobox=True, ): - colors = [ "b.", "r.", @@ -1064,7 +1062,6 @@ def plotmodes( linestyle="-", summary=True, ): - Nplots = self.nrun colors = ["b", "g", "r", "c", "m", "y", "k", "b", "g", "r", "c", "m", "y", "k"] @@ -1265,7 +1262,6 @@ def plotmodes( def plotradeigen( self, pp, field="Ni", comp="amp", yscale="linear", xscale="linear" ): - Nplots = self.nrun colors = ["b", "g", "r", "c", "m", "y", "k", "b", "g", "r", "c", "m", "y", "k"] fig1 = plt.figure() @@ -1368,7 +1364,6 @@ def plotmodes2( xrange=1, debug=False, ): - Nplots = self.nrun Modes = subset(self.db, "field", [field]) # pick field colors = ["b", "g", "r", "c", "m", "y", "k", "b", "g", "r", "c", "m", "y", "k"] @@ -1464,7 +1459,6 @@ def plotMacroDep( def savemovie( self, field="Ni", yscale="log", xscale="log", moviename="spectrum.avi" ): - print("Making movie animation.mpg - this make take a while") files = [] @@ -1504,7 +1498,6 @@ def savemovie( os.system("rm *png") def printmeta(self, pp, filename="output2.pdf", debug=False): - import os from pyPdf import PdfFileWriter, PdfFileReader diff --git a/tools/pylib/post_bout/pb_nonlinear.py b/tools/pylib/post_bout/pb_nonlinear.py index c14072cff5..3fb726d4f8 100644 --- a/tools/pylib/post_bout/pb_nonlinear.py +++ b/tools/pylib/post_bout/pb_nonlinear.py @@ -46,7 +46,6 @@ def plotnlrhs( xscale="linear", xrange=1, ): - colors = ["b", "g", "r", "c", "m", "y", "k", "b", "g", "r", "c", "m", "y", "k"] Modes = subset(self.db, "field", [field]) # pick field diff --git a/tools/pylib/post_bout/pb_present.py b/tools/pylib/post_bout/pb_present.py index e1fc92ea6d..64fcaf1ec6 100644 --- a/tools/pylib/post_bout/pb_present.py +++ b/tools/pylib/post_bout/pb_present.py @@ -96,7 +96,6 @@ def show( for j in list( set(s.dz).union() ): # looping over runs, over unique 'dz' key values - ss = subset(s.db, "dz", [j]) # subset where dz = j plt.scatter(ss.MN[:, 1], ss.MN[:, 0], c=colors[i]) plt.annotate(str(j), (ss.MN[0, 1], ss.MN[0, 0])) diff --git a/tools/pylib/post_bout/read_cxx.py b/tools/pylib/post_bout/read_cxx.py index 907edc48a7..eda88aac5b 100644 --- a/tools/pylib/post_bout/read_cxx.py +++ b/tools/pylib/post_bout/read_cxx.py @@ -98,7 +98,6 @@ def get_evolved_cxx(cxxfile=None): def read_cxx(path=".", boutcxx="physics_code.cxx.ref", evolved=""): - # print path, boutcxx boutcxx = path + "/" + boutcxx # boutcxx = open(boutcxx,'r').readlines() diff --git a/tools/pylib/post_bout/read_inp.py b/tools/pylib/post_bout/read_inp.py index 40a60774af..87de7ddf3a 100644 --- a/tools/pylib/post_bout/read_inp.py +++ b/tools/pylib/post_bout/read_inp.py @@ -12,7 +12,6 @@ def read_inp(path="", boutinp="BOUT.inp"): - boutfile = path + "/" + boutinp boutinp = open(boutfile, "r").readlines() @@ -29,7 +28,6 @@ def read_inp(path="", boutinp="BOUT.inp"): def parse_inp(boutlist): - import re from ordereddict import OrderedDict @@ -67,7 +65,6 @@ def parse_inp(boutlist): def read_log(path=".", logname="status.log"): - print("in read_log") import re from ordereddict import OrderedDict diff --git a/tools/pylib/post_bout/rms.py b/tools/pylib/post_bout/rms.py index 9ec23d9f90..6a9bdb1929 100644 --- a/tools/pylib/post_bout/rms.py +++ b/tools/pylib/post_bout/rms.py @@ -14,7 +14,6 @@ def rms(f): - nt = f.shape[0] ns = f.shape[1] diff --git a/tools/tokamak_grids/elite/elite2nc b/tools/tokamak_grids/elite/elite2nc index eb17c13bd9..669c36aef9 100755 --- a/tools/tokamak_grids/elite/elite2nc +++ b/tools/tokamak_grids/elite/elite2nc @@ -57,6 +57,7 @@ if not desc: print("Description: " + desc) + # Define a generator to get the next token from the file def file_tokens(fp): toklist = [] diff --git a/tools/tokamak_grids/gato/gato2nc b/tools/tokamak_grids/gato/gato2nc index ba4cdc6e69..4ed2b2d632 100755 --- a/tools/tokamak_grids/gato/gato2nc +++ b/tools/tokamak_grids/gato/gato2nc @@ -68,6 +68,7 @@ print("Date: " + date) desc = f.readline() print("Description: " + desc) + # Define a generator to get the next token from the file def file_tokens(fp): """Generator to get numbers from a text file""" From a088700cd86a4a1244a9159b0f15dd4f62fa40da Mon Sep 17 00:00:00 2001 From: dschwoerer Date: Wed, 8 Feb 2023 10:19:40 +0000 Subject: [PATCH 0022/1827] Apply clang-format changes --- include/interpolation_xz.hxx | 21 ++- include/mask.hxx | 3 +- src/mesh/interpolation/hermite_spline_xz.cxx | 84 +++++----- src/mesh/parallel/fci.cxx | 12 +- tests/integrated/test-fci-mpi/fci_mpi.cxx | 27 ++-- .../test-interpolate/test_interpolate.cxx | 146 +++++++++--------- 6 files changed, 147 insertions(+), 146 deletions(-) diff --git a/include/interpolation_xz.hxx b/include/interpolation_xz.hxx index 620d146edf..df4c7fc61a 100644 --- a/include/interpolation_xz.hxx +++ b/include/interpolation_xz.hxx @@ -63,25 +63,23 @@ public: setMask(mask); } XZInterpolation(const std::string& region_name, int y_offset = 0, Mesh* mesh = nullptr) - : y_offset(y_offset), localmesh(mesh), region_id(localmesh->getRegionID(region_name)) {} - XZInterpolation(const Region& region, int y_offset = 0, - Mesh* mesh = nullptr) - : y_offset(y_offset), localmesh(mesh){ + : y_offset(y_offset), localmesh(mesh), + region_id(localmesh->getRegionID(region_name)) {} + XZInterpolation(const Region& region, int y_offset = 0, Mesh* mesh = nullptr) + : y_offset(y_offset), localmesh(mesh) { setRegion(region); } virtual ~XZInterpolation() = default; - void setMask(const BoutMask& mask) { - setRegion(regionFromMask(mask, localmesh)); - } + void setMask(const BoutMask& mask) { setRegion(regionFromMask(mask, localmesh)); } void setRegion(const std::string& region_name) { this->region_id = localmesh->getRegionID(region_name); } void setRegion(const Region& region) { std::string name; - int i=0; + int i = 0; do { - name = fmt::format("unsec_reg_xz_interp_{:d}",i++); + name = fmt::format("unsec_reg_xz_interp_{:d}", i++); } while (localmesh->hasRegion3D(name)); localmesh->addRegion(name, region); this->region_id = localmesh->getRegionID(name); @@ -94,10 +92,11 @@ public: if (region_id == -1) { return localmesh->getRegion(region); } - if (region == "" or region == "RGN_ALL"){ + if (region == "" or region == "RGN_ALL") { return getRegion(); } - return localmesh->getRegion(localmesh->getCommonRegion(localmesh->getRegionID(region), region_id)); + return localmesh->getRegion( + localmesh->getCommonRegion(localmesh->getRegionID(region), region_id)); } virtual void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const std::string& region = "RGN_NOBNDRY") = 0; diff --git a/include/mask.hxx b/include/mask.hxx index 96d2c99ac3..20211b5d02 100644 --- a/include/mask.hxx +++ b/include/mask.hxx @@ -73,8 +73,7 @@ public: inline const bool& operator[](const Ind3D& i) const { return mask[i]; } }; -inline Region regionFromMask(const BoutMask& mask, - const Mesh* mesh) { +inline Region regionFromMask(const BoutMask& mask, const Mesh* mesh) { std::vector indices; for (auto i : mesh->getRegion("RGN_ALL")) { if (not mask(i.x(), i.y(), i.z())) { diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index 4e6fc8320c..a5b9c8bd05 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -20,10 +20,10 @@ * **************************************************************************/ +#include "../impls/bout/boutmesh.hxx" #include "globals.hxx" #include "interpolation_xz.hxx" #include "bout/index_derivs_interface.hxx" -#include "../impls/bout/boutmesh.hxx" #include @@ -126,7 +126,7 @@ XZHermiteSpline::XZHermiteSpline(int y_offset, Mesh *mesh) #if USE_NEW_WEIGHTS newWeights.reserve(16); - for (int w=0; w<16;++w){ + for (int w = 0; w < 16; ++w) { newWeights.emplace_back(localmesh); newWeights[w].allocate(); } @@ -137,8 +137,9 @@ XZHermiteSpline::XZHermiteSpline(int y_offset, Mesh *mesh) // MatCreate(MPI_COMM_WORLD, &petscWeights); // MatSetSizes(petscWeights, m, m, M, M); // PetscErrorCode MatCreateAIJ(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, - // PetscInt N, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt - //o_nnz[], Mat *A) + // PetscInt N, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, + // const PetscInt + // o_nnz[], Mat *A) // MatSetSizes(Mat A,PetscInt m,PetscInt n,PetscInt M,PetscInt N) const int m = mesh->LocalNx * mesh->LocalNy * mesh->LocalNz; const int M = m * mesh->getNXPE() * mesh->getNYPE(); @@ -185,8 +186,8 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z // Check that t_x and t_z are in range if ((t_x < 0.0) || (t_x > 1.0)) { throw BoutException( - "t_x={:e} out of range at ({:d},{:d},{:d}) (delta_x={:e}, i_corn={:d})", t_x, - x, y, z, delta_x(x, y, z), i_corn); + "t_x={:e} out of range at ({:d},{:d},{:d}) (delta_x={:e}, i_corn={:d})", t_x, x, + y, z, delta_x(x, y, z), i_corn); } if ((t_z < 0.0) || (t_z > 1.0)) { @@ -212,8 +213,8 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z #if USE_NEW_WEIGHTS - for (int w =0; w<16;++w){ - newWeights[w][i]=0; + for (int w = 0; w < 16; ++w) { + newWeights[w][i] = 0; } // The distribution of our weights: // 0 4 8 12 @@ -223,54 +224,55 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z // e.g. 1 == ic.xm(); 4 == ic.zm(); 5 == ic; 7 == ic.zp(2); // f[ic] * h00_x[i] + f[icxp] * h01_x[i] + fx[ic] * h10_x[i] + fx[icxp] * h11_x[i]; - newWeights[5][i] += h00_x[i] * h00_z[i]; - newWeights[9][i] += h01_x[i] * h00_z[i]; - newWeights[9][i] += h10_x[i] * h00_z[i] / 2; - newWeights[1][i] -= h10_x[i] * h00_z[i] / 2; + newWeights[5][i] += h00_x[i] * h00_z[i]; + newWeights[9][i] += h01_x[i] * h00_z[i]; + newWeights[9][i] += h10_x[i] * h00_z[i] / 2; + newWeights[1][i] -= h10_x[i] * h00_z[i] / 2; newWeights[13][i] += h11_x[i] * h00_z[i] / 2; - newWeights[5][i] -= h11_x[i] * h00_z[i] / 2; + newWeights[5][i] -= h11_x[i] * h00_z[i] / 2; // f[iczp] * h00_x[i] + f[icxpzp] * h01_x[i] + // fx[iczp] * h10_x[i] + fx[icxpzp] * h11_x[i]; - newWeights[6][i] += h00_x[i] * h01_z[i]; + newWeights[6][i] += h00_x[i] * h01_z[i]; newWeights[10][i] += h01_x[i] * h01_z[i]; newWeights[10][i] += h10_x[i] * h01_z[i] / 2; - newWeights[2][i] -= h10_x[i] * h01_z[i] / 2; + newWeights[2][i] -= h10_x[i] * h01_z[i] / 2; newWeights[14][i] += h11_x[i] * h01_z[i] / 2; - newWeights[6][i] -= h11_x[i] * h01_z[i] / 2; + newWeights[6][i] -= h11_x[i] * h01_z[i] / 2; // fz[ic] * h00_x[i] + fz[icxp] * h01_x[i] + // fxz[ic] * h10_x[i]+ fxz[icxp] * h11_x[i]; - newWeights[6][i] += h00_x[i] * h10_z[i] / 2; - newWeights[4][i] -= h00_x[i] * h10_z[i] / 2; + newWeights[6][i] += h00_x[i] * h10_z[i] / 2; + newWeights[4][i] -= h00_x[i] * h10_z[i] / 2; newWeights[10][i] += h01_x[i] * h10_z[i] / 2; - newWeights[8][i] -= h01_x[i] * h10_z[i] / 2; + newWeights[8][i] -= h01_x[i] * h10_z[i] / 2; newWeights[10][i] += h10_x[i] * h10_z[i] / 4; - newWeights[8][i] -= h10_x[i] * h10_z[i] / 4; - newWeights[2][i] -= h10_x[i] * h10_z[i] / 4; - newWeights[0][i] += h10_x[i] * h10_z[i] / 4; + newWeights[8][i] -= h10_x[i] * h10_z[i] / 4; + newWeights[2][i] -= h10_x[i] * h10_z[i] / 4; + newWeights[0][i] += h10_x[i] * h10_z[i] / 4; newWeights[14][i] += h11_x[i] * h10_z[i] / 4; newWeights[12][i] -= h11_x[i] * h10_z[i] / 4; - newWeights[6][i] -= h11_x[i] * h10_z[i] / 4; - newWeights[4][i] += h11_x[i] * h10_z[i] / 4; + newWeights[6][i] -= h11_x[i] * h10_z[i] / 4; + newWeights[4][i] += h11_x[i] * h10_z[i] / 4; // fz[iczp] * h00_x[i] + fz[icxpzp] * h01_x[i] + // fxz[iczp] * h10_x[i] + fxz[icxpzp] * h11_x[i]; - newWeights[7][i] += h00_x[i] * h11_z[i] / 2; - newWeights[5][i] -= h00_x[i] * h11_z[i] / 2; + newWeights[7][i] += h00_x[i] * h11_z[i] / 2; + newWeights[5][i] -= h00_x[i] * h11_z[i] / 2; newWeights[11][i] += h01_x[i] * h11_z[i] / 2; - newWeights[9][i] -= h01_x[i] * h11_z[i] / 2; + newWeights[9][i] -= h01_x[i] * h11_z[i] / 2; newWeights[11][i] += h10_x[i] * h11_z[i] / 4; - newWeights[9][i] -= h10_x[i] * h11_z[i] / 4; - newWeights[3][i] -= h10_x[i] * h11_z[i] / 4; - newWeights[1][i] += h10_x[i] * h11_z[i] / 4; + newWeights[9][i] -= h10_x[i] * h11_z[i] / 4; + newWeights[3][i] -= h10_x[i] * h11_z[i] / 4; + newWeights[1][i] += h10_x[i] * h11_z[i] / 4; newWeights[15][i] += h11_x[i] * h11_z[i] / 4; newWeights[13][i] -= h11_x[i] * h11_z[i] / 4; - newWeights[7][i] -= h11_x[i] * h11_z[i] / 4; - newWeights[5][i] += h11_x[i] * h11_z[i] / 4; + newWeights[7][i] -= h11_x[i] * h11_z[i] / 4; + newWeights[5][i] += h11_x[i] * h11_z[i] / 4; #ifdef HS_USE_PETSC PetscInt idxn[1] = {conv.fromLocalToGlobal(x, y + y_offset, z)}; - // output.write("debug: {:d} -> {:d}: {:d}:{:d} -> {:d}:{:d}\n", conv.fromLocalToGlobal(x, y + y_offset, z), + // output.write("debug: {:d} -> {:d}: {:d}:{:d} -> {:d}:{:d}\n", + // conv.fromLocalToGlobal(x, y + y_offset, z), // conv.fromMeshToGlobal(i_corn, y + y_offset, k_corner(x, y, z)), // x, z, i_corn, k_corner(x, y, z)); // ixstep = mesh->LocalNx * mesh->LocalNz; @@ -355,15 +357,15 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region VecRestoreArrayRead(result, &cptr); #else BOUT_FOR(i, getRegion(region)) { - auto ic = i_corner[i]; + auto ic = i_corner[i]; auto iyp = i.yp(y_offset); - f_interp[iyp]=0; - for (int w = 0; w < 4; ++w){ - f_interp[iyp] += newWeights[w*4+0][i] * f[ic.zm().xp(w-1)]; - f_interp[iyp] += newWeights[w*4+1][i] * f[ic.xp(w-1)]; - f_interp[iyp] += newWeights[w*4+2][i] * f[ic.zp().xp(w-1)]; - f_interp[iyp] += newWeights[w*4+3][i] * f[ic.zp(2).xp(w-1)]; + f_interp[iyp] = 0; + for (int w = 0; w < 4; ++w) { + f_interp[iyp] += newWeights[w * 4 + 0][i] * f[ic.zm().xp(w - 1)]; + f_interp[iyp] += newWeights[w * 4 + 1][i] * f[ic.xp(w - 1)]; + f_interp[iyp] += newWeights[w * 4 + 2][i] * f[ic.zp().xp(w - 1)]; + f_interp[iyp] += newWeights[w * 4 + 3][i] * f[ic.zp(2).xp(w - 1)]; } } #endif @@ -411,7 +413,7 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region || i.x() > localmesh->xend); } return f_interp; -# endif +#endif } Field3D XZHermiteSpline::interpolate(const Field3D& f, const Field3D& delta_x, diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index baf0f3fc6b..4b964ae0fa 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -158,7 +158,8 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& dy, Options& options, const int ncz = map_mesh.LocalNz; BoutMask to_remove(map_mesh); - const int xend = map_mesh.xstart + (map_mesh.xend - map_mesh.xstart + 1) * map_mesh.getNXPE() - 1; + const int xend = + map_mesh.xstart + (map_mesh.xend - map_mesh.xstart + 1) * map_mesh.getNXPE() - 1; // Serial loop because call to BoundaryRegionPar::addPoint // (probably?) can't be done in parallel BOUT_FOR_SERIAL(i, xt_prime.getRegion("RGN_NOBNDRY")) { @@ -247,11 +248,10 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& dy, Options& options, const auto region = fmt::format("RGN_YPAR_{:+d}", offset); if (not map_mesh.hasRegion3D(region)) { // The valid region for this slice - map_mesh.addRegion3D(region, - Region(map_mesh.xstart, map_mesh.xend, - map_mesh.ystart+offset, map_mesh.yend+offset, - 0, map_mesh.LocalNz-1, - map_mesh.LocalNy, map_mesh.LocalNz)); + map_mesh.addRegion3D( + region, Region(map_mesh.xstart, map_mesh.xend, map_mesh.ystart + offset, + map_mesh.yend + offset, 0, map_mesh.LocalNz - 1, + map_mesh.LocalNy, map_mesh.LocalNz)); } } diff --git a/tests/integrated/test-fci-mpi/fci_mpi.cxx b/tests/integrated/test-fci-mpi/fci_mpi.cxx index b353493dda..6ae711351e 100644 --- a/tests/integrated/test-fci-mpi/fci_mpi.cxx +++ b/tests/integrated/test-fci-mpi/fci_mpi.cxx @@ -6,28 +6,29 @@ int main(int argc, char** argv) { BoutInitialise(argc, argv); { using bout::globals::mesh; - Options *options = Options::getRoot(); - int i=0; - std::string default_str {"not_set"}; + Options* options = Options::getRoot(); + int i = 0; + std::string default_str{"not_set"}; Options dump; while (true) { std::string temp_str; options->get(fmt::format("input_{:d}:function", i), temp_str, default_str); if (temp_str == default_str) { - break; + break; } - Field3D input{FieldFactory::get()->create3D(fmt::format("input_{:d}:function", i), Options::getRoot(), mesh)}; - //options->get(fmt::format("input_{:d}:boundary_perp", i), temp_str, s"free_o3"); + Field3D input{FieldFactory::get()->create3D(fmt::format("input_{:d}:function", i), + Options::getRoot(), mesh)}; + // options->get(fmt::format("input_{:d}:boundary_perp", i), temp_str, s"free_o3"); mesh->communicate(input); input.applyParallelBoundary("parallel_neumann_o2"); for (int slice = -mesh->ystart; slice <= mesh->ystart; ++slice) { - if (slice) { - Field3D tmp{0.}; - BOUT_FOR(i, tmp.getRegion("RGN_NOBNDRY")) { - tmp[i] = input.ynext(slice)[i.yp(slice)]; - } - dump[fmt::format("output_{:d}_{:+d}", i, slice)] = tmp; - } + if (slice) { + Field3D tmp{0.}; + BOUT_FOR(i, tmp.getRegion("RGN_NOBNDRY")) { + tmp[i] = input.ynext(slice)[i.yp(slice)]; + } + dump[fmt::format("output_{:d}_{:+d}", i, slice)] = tmp; + } } ++i; } diff --git a/tests/integrated/test-interpolate/test_interpolate.cxx b/tests/integrated/test-interpolate/test_interpolate.cxx index 517d9c2445..a14208e7b3 100644 --- a/tests/integrated/test-interpolate/test_interpolate.cxx +++ b/tests/integrated/test-interpolate/test_interpolate.cxx @@ -31,88 +31,88 @@ std::shared_ptr getGeneratorFromOptions(const std::string& varna int main(int argc, char **argv) { BoutInitialise(argc, argv); { - // Random number generator - std::default_random_engine generator; - // Uniform distribution of BoutReals from 0 to 1 - std::uniform_real_distribution distribution{0.0, 1.0}; - - using bout::globals::mesh; - - FieldFactory f(mesh); - - // Set up generators and solutions for three different analtyic functions - std::string a_func; - auto a_gen = getGeneratorFromOptions("a", a_func); - Field3D a = f.create3D(a_func); - Field3D a_solution = 0.0; - Field3D a_interp = 0.0; - - std::string b_func; - auto b_gen = getGeneratorFromOptions("b", b_func); - Field3D b = f.create3D(b_func); - Field3D b_solution = 0.0; - Field3D b_interp = 0.0; - - std::string c_func; - auto c_gen = getGeneratorFromOptions("c", c_func); - Field3D c = f.create3D(c_func); - Field3D c_solution = 0.0; - Field3D c_interp = 0.0; - - // x and z displacements - Field3D deltax = 0.0; - Field3D deltaz = 0.0; - - // Bind the random number generator and distribution into a single function - auto dice = std::bind(distribution, generator); - - for (const auto &index : deltax) { - // Get some random displacements - BoutReal dx = index.x() + dice(); - BoutReal dz = index.z() + dice(); - // For the last point, put the displacement inwards - // Otherwise we try to interpolate in the guard cells, which doesn't work so well - if (index.x() >= mesh->xend && mesh->getNXPE() - 1 == mesh->getXProcIndex()) { - dx = index.x() - dice(); + // Random number generator + std::default_random_engine generator; + // Uniform distribution of BoutReals from 0 to 1 + std::uniform_real_distribution distribution{0.0, 1.0}; + + using bout::globals::mesh; + + FieldFactory f(mesh); + + // Set up generators and solutions for three different analtyic functions + std::string a_func; + auto a_gen = getGeneratorFromOptions("a", a_func); + Field3D a = f.create3D(a_func); + Field3D a_solution = 0.0; + Field3D a_interp = 0.0; + + std::string b_func; + auto b_gen = getGeneratorFromOptions("b", b_func); + Field3D b = f.create3D(b_func); + Field3D b_solution = 0.0; + Field3D b_interp = 0.0; + + std::string c_func; + auto c_gen = getGeneratorFromOptions("c", c_func); + Field3D c = f.create3D(c_func); + Field3D c_solution = 0.0; + Field3D c_interp = 0.0; + + // x and z displacements + Field3D deltax = 0.0; + Field3D deltaz = 0.0; + + // Bind the random number generator and distribution into a single function + auto dice = std::bind(distribution, generator); + + for (const auto& index : deltax) { + // Get some random displacements + BoutReal dx = index.x() + dice(); + BoutReal dz = index.z() + dice(); + // For the last point, put the displacement inwards + // Otherwise we try to interpolate in the guard cells, which doesn't work so well + if (index.x() >= mesh->xend && mesh->getNXPE() - 1 == mesh->getXProcIndex()) { + dx = index.x() - dice(); + } + deltax[index] = dx; + deltaz[index] = dz; + // Get the global indices + bout::generator::Context pos{index, CELL_CENTRE, deltax.getMesh(), 0.0}; + pos.set("x", mesh->GlobalX(dx), "z", + TWOPI * static_cast(dz) / static_cast(mesh->LocalNz)); + // Generate the analytic solution at the displacements + a_solution[index] = a_gen->generate(pos); + b_solution[index] = b_gen->generate(pos); + c_solution[index] = c_gen->generate(pos); } - deltax[index] = dx; - deltaz[index] = dz; - // Get the global indices - bout::generator::Context pos{index, CELL_CENTRE, deltax.getMesh(), 0.0}; - pos.set("x", mesh->GlobalX(dx), - "z", TWOPI * static_cast(dz) / static_cast(mesh->LocalNz)); - // Generate the analytic solution at the displacements - a_solution[index] = a_gen->generate(pos); - b_solution[index] = b_gen->generate(pos); - c_solution[index] = c_gen->generate(pos); - } - deltax += (mesh->LocalNx - mesh->xstart * 2) * mesh->getXProcIndex(); - // Create the interpolation object from the input options - auto interp = XZInterpolationFactory::getInstance().create(); + deltax += (mesh->LocalNx - mesh->xstart * 2) * mesh->getXProcIndex(); + // Create the interpolation object from the input options + auto interp = XZInterpolationFactory::getInstance().create(); - // Interpolate the analytic functions at the displacements - a_interp = interp->interpolate(a, deltax, deltaz); - b_interp = interp->interpolate(b, deltax, deltaz); - c_interp = interp->interpolate(c, deltax, deltaz); + // Interpolate the analytic functions at the displacements + a_interp = interp->interpolate(a, deltax, deltaz); + b_interp = interp->interpolate(b, deltax, deltaz); + c_interp = interp->interpolate(c, deltax, deltaz); - Options dump; + Options dump; - dump["a"] = a; - dump["a_interp"] = a_interp; - dump["a_solution"] = a_solution; + dump["a"] = a; + dump["a_interp"] = a_interp; + dump["a_solution"] = a_solution; - dump["b"] = b; - dump["b_interp"] = b_interp; - dump["b_solution"] = b_solution; + dump["b"] = b; + dump["b_interp"] = b_interp; + dump["b_solution"] = b_solution; - dump["c"] = c; - dump["c_interp"] = c_interp; - dump["c_solution"] = c_solution; + dump["c"] = c; + dump["c_interp"] = c_interp; + dump["c_solution"] = c_solution; - bout::writeDefaultOutputFile(dump); + bout::writeDefaultOutputFile(dump); - bout::checkForUnusedOptions(); + bout::checkForUnusedOptions(); } BoutFinalise(); From 32ea2fdcf7a10efb86e9c8541d6dcb47aaa1f538 Mon Sep 17 00:00:00 2001 From: dschwoerer Date: Wed, 8 Feb 2023 13:10:31 +0000 Subject: [PATCH 0023/1827] Apply clang-format changes --- src/mesh/interpolation/hermite_spline_xz.cxx | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index a5b9c8bd05..2d1649a7c1 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -137,9 +137,8 @@ XZHermiteSpline::XZHermiteSpline(int y_offset, Mesh *mesh) // MatCreate(MPI_COMM_WORLD, &petscWeights); // MatSetSizes(petscWeights, m, m, M, M); // PetscErrorCode MatCreateAIJ(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, - // PetscInt N, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, - // const PetscInt - // o_nnz[], Mat *A) + // PetscInt N, PetscInt d_nz, const PetscInt d_nnz[], PetscInt + // o_nz, const PetscInt o_nnz[], Mat *A) // MatSetSizes(Mat A,PetscInt m,PetscInt n,PetscInt M,PetscInt N) const int m = mesh->LocalNx * mesh->LocalNy * mesh->LocalNz; const int M = m * mesh->getNXPE() * mesh->getNYPE(); From faac69c321c957d542b7631e900d34e690238f02 Mon Sep 17 00:00:00 2001 From: dschwoerer Date: Wed, 8 Feb 2023 13:10:52 +0000 Subject: [PATCH 0024/1827] Apply clang-format changes --- src/mesh/interpolation/hermite_spline_xz.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index 2d1649a7c1..93f1c326e8 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -137,8 +137,8 @@ XZHermiteSpline::XZHermiteSpline(int y_offset, Mesh *mesh) // MatCreate(MPI_COMM_WORLD, &petscWeights); // MatSetSizes(petscWeights, m, m, M, M); // PetscErrorCode MatCreateAIJ(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, - // PetscInt N, PetscInt d_nz, const PetscInt d_nnz[], PetscInt - // o_nz, const PetscInt o_nnz[], Mat *A) + // PetscInt N, PetscInt d_nz, const PetscInt d_nnz[], + // PetscInt o_nz, const PetscInt o_nnz[], Mat *A) // MatSetSizes(Mat A,PetscInt m,PetscInt n,PetscInt M,PetscInt N) const int m = mesh->LocalNx * mesh->LocalNy * mesh->LocalNz; const int M = m * mesh->getNXPE() * mesh->getNYPE(); From cc1672d1b9e28af98089fcdff59d6c2f8aa8594d Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 14 Jun 2023 12:06:27 -0700 Subject: [PATCH 0025/1827] Add Array to Options Enable 1D arrays of integers to be read from NetCDF files, and made available via `Mesh::get()`. Mainly useful for mesh construction. --- include/bout/options.hxx | 6 +++++- include/bout/utils.hxx | 4 ++-- src/mesh/data/gridfromfile.cxx | 20 +++++++++++++------ src/sys/options.cxx | 31 ++++++++++++++++++++++++++++++ src/sys/options/options_netcdf.cxx | 4 ++++ 5 files changed, 56 insertions(+), 9 deletions(-) diff --git a/include/bout/options.hxx b/include/bout/options.hxx index bf1704100f..fae05b8ecc 100644 --- a/include/bout/options.hxx +++ b/include/bout/options.hxx @@ -196,7 +196,7 @@ public: /// The type used to store values using ValueType = bout::utils::variant, Matrix, Tensor>; + Array, Array, Matrix, Tensor>; /// The type used to store attributes /// Extends the variant class so that cast operator can be implemented @@ -874,6 +874,8 @@ void Options::assign<>(FieldPerp val, std::string source); template <> void Options::assign<>(Array val, std::string source); template <> +void Options::assign<>(Array val, std::string source); +template <> void Options::assign<>(Matrix val, std::string source); template <> void Options::assign<>(Tensor val, std::string source); @@ -902,6 +904,8 @@ FieldPerp Options::as(const FieldPerp& similar_to) const; template <> Array Options::as>(const Array& similar_to) const; template <> +Array Options::as>(const Array& similar_to) const; +template <> Matrix Options::as>(const Matrix& similar_to) const; template <> Tensor Options::as>(const Tensor& similar_to) const; diff --git a/include/bout/utils.hxx b/include/bout/utils.hxx index 450d3eaf87..e2ec5b24db 100644 --- a/include/bout/utils.hxx +++ b/include/bout/utils.hxx @@ -543,8 +543,8 @@ std::string toString(const T& val) { /// where the type may be std::string. inline std::string toString(const std::string& val) { return val; } -template <> -inline std::string toString<>(const Array& UNUSED(val)) { +template +inline std::string toString(const Array& UNUSED(val)) { return ""; } diff --git a/src/mesh/data/gridfromfile.cxx b/src/mesh/data/gridfromfile.cxx index 7e875bd109..7fab486e54 100644 --- a/src/mesh/data/gridfromfile.cxx +++ b/src/mesh/data/gridfromfile.cxx @@ -138,7 +138,8 @@ struct GetDimensions { std::vector operator()(MAYBE_UNUSED(int value)) { return {1}; } std::vector operator()(MAYBE_UNUSED(BoutReal value)) { return {1}; } std::vector operator()(MAYBE_UNUSED(const std::string& value)) { return {1}; } - std::vector operator()(const Array& array) { return {array.size()}; } + template + std::vector operator()(const Array& array) { return {array.size()}; } std::vector operator()(const Matrix& array) { const auto shape = array.shape(); return {std::get<0>(shape), std::get<1>(shape)}; @@ -471,13 +472,20 @@ void GridFile::readField(Mesh* m, const std::string& name, int UNUSED(ys), int U } } -bool GridFile::get(MAYBE_UNUSED(Mesh* m), MAYBE_UNUSED(std::vector& var), - MAYBE_UNUSED(const std::string& name), MAYBE_UNUSED(int len), - MAYBE_UNUSED(int offset), - MAYBE_UNUSED(GridDataSource::Direction dir)) { +bool GridFile::get(Mesh* UNUSED(m), std::vector& var, const std::string& name, + int len, int offset, GridDataSource::Direction UNUSED(dir)) { TRACE("GridFile::get(vector)"); - return false; + if (not data.isSet(name)) { + return false; + } + + const auto full_var = data[name].as>(); + const auto* it = std::begin(full_var); + std::advance(it, offset); + std::copy_n(it, len, std::begin(var)); + + return true; } bool GridFile::get(Mesh* UNUSED(m), std::vector& var, const std::string& name, diff --git a/src/sys/options.cxx b/src/sys/options.cxx index 8b49b1f3f1..fec4073682 100644 --- a/src/sys/options.cxx +++ b/src/sys/options.cxx @@ -288,6 +288,10 @@ void Options::assign<>(Array val, std::string source) { _set_no_check(std::move(val), std::move(source)); } template <> +void Options::assign<>(Array val, std::string source) { + _set_no_check(std::move(val), std::move(source)); +} +template <> void Options::assign<>(Matrix val, std::string source) { _set_no_check(std::move(val), std::move(source)); } @@ -723,6 +727,33 @@ Array Options::as>(const Array& similar_to) return result; } +template <> +Array Options::as>(const Array& similar_to) const { + if (is_section) { + throw BoutException(_("Option {:s} has no value"), full_name); + } + + Array result = bout::utils::visit( + ConvertContainer>{ + fmt::format( + _("Value for option {:s} cannot be converted to an Array"), + full_name), + similar_to}, + value); + + // Mark this option as used + value_used = true; + + output_info << _("\tOption ") << full_name << " = Array"; + if (hasAttribute("source")) { + // Specify the source of the setting + output_info << " (" << bout::utils::variantToString(attributes.at("source")) << ")"; + } + output_info << endl; + + return result; +} + template <> Matrix Options::as>(const Matrix& similar_to) const { if (is_section) { diff --git a/src/sys/options/options_netcdf.cxx b/src/sys/options/options_netcdf.cxx index d7ceeaea60..7c21b309bc 100644 --- a/src/sys/options/options_netcdf.cxx +++ b/src/sys/options/options_netcdf.cxx @@ -88,6 +88,10 @@ void readGroup(const std::string& filename, const NcGroup& group, Options& resul Array value(static_cast(dims[0].getSize())); var.getVar(value.begin()); result[var_name] = value; + } else if (var_type == ncInt or var_type == ncShort) { + Array value(static_cast(dims[0].getSize())); + var.getVar(value.begin()); + result[var_name] = value; } else if ((var_type == ncString) or (var_type == ncChar)) { std::string value; value.resize(dims[0].getSize()); From e654af65c6ea32a0ee4f6f07659e246ac70ac88b Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 14 Jun 2023 18:06:49 -0700 Subject: [PATCH 0026/1827] GridFile: Check sizes, allocate vector ints When reading 1D vector of ints, check that the input is long enough, and resize the output vector. --- src/mesh/data/gridfromfile.cxx | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/mesh/data/gridfromfile.cxx b/src/mesh/data/gridfromfile.cxx index 7fab486e54..7c02cf4ab7 100644 --- a/src/mesh/data/gridfromfile.cxx +++ b/src/mesh/data/gridfromfile.cxx @@ -481,6 +481,16 @@ bool GridFile::get(Mesh* UNUSED(m), std::vector& var, const std::string& na } const auto full_var = data[name].as>(); + + // Check size + if (full_var.size() < len + offset) { + throw BoutException("{} has length {}. Expected {} elements + {} offset", + name, full_var.size(), len, offset); + } + + // Ensure that output variable has the correct size + var.resize(len); + const auto* it = std::begin(full_var); std::advance(it, offset); std::copy_n(it, len, std::begin(var)); From 359766313a8b878f056008358192cb4653130384 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 14 Jun 2023 18:08:32 -0700 Subject: [PATCH 0027/1827] Add limiters in input mesh A short-term fix, allowing extra limiters to be added with inputs: limiter_count : int limiter_yinds : [int], length limiter_count limiter_xstarts : [int], length limiter_count limiter_xends : [int], length limiter_count The limiter(s) are added between yinds[i] and yinds[i] + 1 --- src/mesh/impls/bout/boutmesh.cxx | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index a802d3f5b3..448832ea57 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -2073,6 +2073,32 @@ void BoutMesh::topology() { add_target(ny_inner - 1, 0, nx); } + // Additional limiters + // Each limiter needs 3 indices: A Y index, start and end X indices + int limiter_count = 0; + Mesh::get(limiter_count, "limiter_count", 0); + if (limiter_count > 0) { + std::vector limiter_yinds, limiter_xstarts, limiter_xends; + if (!source->get(this, limiter_yinds, "limiter_yinds", limiter_count)) { + throw BoutException("Couldn't read limiter_yinds vector of length {} from mesh", limiter_count); + } + if (!source->get(this, limiter_xstarts, "limiter_xstarts", limiter_count)) { + throw BoutException("Couldn't read limiter_xstarts vector of length {} from mesh", limiter_count); + } + if (!source->get(this, limiter_xends, "limiter_xends", limiter_count)) { + throw BoutException("Couldn't read limiter_xend vector of length {} from mesh", limiter_count); + } + + for (int i = 0; i < limiter_count; ++i) { + int yind = limiter_yinds[i]; + int xstart = limiter_xstarts[i]; + int xend = limiter_xends[i]; + output_info.write("Adding a limiter between y={} and {}. X indices {} to {}\n", + yind, yind+1, xstart, xend); + add_target(yind, xstart, xend); + } + } + if ((ixseps_inner > 0) && (((PE_YIND * MYSUB > jyseps1_1) && (PE_YIND * MYSUB <= jyseps2_1)) || ((PE_YIND * MYSUB > jyseps1_2) && (PE_YIND * MYSUB <= jyseps2_2)))) { From 1c1d834c822857239efa5f8189b187845cebe073 Mon Sep 17 00:00:00 2001 From: bendudson Date: Thu, 15 Jun 2023 01:17:21 +0000 Subject: [PATCH 0028/1827] Apply clang-format changes --- include/bout/options.hxx | 6 +++--- src/mesh/data/gridfromfile.cxx | 10 ++++++---- src/mesh/impls/bout/boutmesh.cxx | 11 +++++++---- src/sys/options.cxx | 5 ++--- 4 files changed, 18 insertions(+), 14 deletions(-) diff --git a/include/bout/options.hxx b/include/bout/options.hxx index fae05b8ecc..e790b405b5 100644 --- a/include/bout/options.hxx +++ b/include/bout/options.hxx @@ -194,9 +194,9 @@ public: static void cleanup(); /// The type used to store values - using ValueType = - bout::utils::variant, Array, Matrix, Tensor>; + using ValueType = bout::utils::variant, Array, + Matrix, Tensor>; /// The type used to store attributes /// Extends the variant class so that cast operator can be implemented diff --git a/src/mesh/data/gridfromfile.cxx b/src/mesh/data/gridfromfile.cxx index 7c02cf4ab7..eed70776b3 100644 --- a/src/mesh/data/gridfromfile.cxx +++ b/src/mesh/data/gridfromfile.cxx @@ -138,8 +138,10 @@ struct GetDimensions { std::vector operator()(MAYBE_UNUSED(int value)) { return {1}; } std::vector operator()(MAYBE_UNUSED(BoutReal value)) { return {1}; } std::vector operator()(MAYBE_UNUSED(const std::string& value)) { return {1}; } - template - std::vector operator()(const Array& array) { return {array.size()}; } + template + std::vector operator()(const Array& array) { + return {array.size()}; + } std::vector operator()(const Matrix& array) { const auto shape = array.shape(); return {std::get<0>(shape), std::get<1>(shape)}; @@ -484,8 +486,8 @@ bool GridFile::get(Mesh* UNUSED(m), std::vector& var, const std::string& na // Check size if (full_var.size() < len + offset) { - throw BoutException("{} has length {}. Expected {} elements + {} offset", - name, full_var.size(), len, offset); + throw BoutException("{} has length {}. Expected {} elements + {} offset", name, + full_var.size(), len, offset); } // Ensure that output variable has the correct size diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index 448832ea57..c2c9b51ba3 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -2080,13 +2080,16 @@ void BoutMesh::topology() { if (limiter_count > 0) { std::vector limiter_yinds, limiter_xstarts, limiter_xends; if (!source->get(this, limiter_yinds, "limiter_yinds", limiter_count)) { - throw BoutException("Couldn't read limiter_yinds vector of length {} from mesh", limiter_count); + throw BoutException("Couldn't read limiter_yinds vector of length {} from mesh", + limiter_count); } if (!source->get(this, limiter_xstarts, "limiter_xstarts", limiter_count)) { - throw BoutException("Couldn't read limiter_xstarts vector of length {} from mesh", limiter_count); + throw BoutException("Couldn't read limiter_xstarts vector of length {} from mesh", + limiter_count); } if (!source->get(this, limiter_xends, "limiter_xends", limiter_count)) { - throw BoutException("Couldn't read limiter_xend vector of length {} from mesh", limiter_count); + throw BoutException("Couldn't read limiter_xend vector of length {} from mesh", + limiter_count); } for (int i = 0; i < limiter_count; ++i) { @@ -2094,7 +2097,7 @@ void BoutMesh::topology() { int xstart = limiter_xstarts[i]; int xend = limiter_xends[i]; output_info.write("Adding a limiter between y={} and {}. X indices {} to {}\n", - yind, yind+1, xstart, xend); + yind, yind + 1, xstart, xend); add_target(yind, xstart, xend); } } diff --git a/src/sys/options.cxx b/src/sys/options.cxx index fec4073682..1c50db5896 100644 --- a/src/sys/options.cxx +++ b/src/sys/options.cxx @@ -735,9 +735,8 @@ Array Options::as>(const Array& similar_to) const { Array result = bout::utils::visit( ConvertContainer>{ - fmt::format( - _("Value for option {:s} cannot be converted to an Array"), - full_name), + fmt::format(_("Value for option {:s} cannot be converted to an Array"), + full_name), similar_to}, value); From 3f311e0896799974a154638ef0ff62b85401becb Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 22 Jan 2024 14:04:00 +0000 Subject: [PATCH 0029/1827] Docs: Fix some double-backtick issues --- manual/sphinx/user_docs/bout_options.rst | 30 ++++++++++++------------ 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/manual/sphinx/user_docs/bout_options.rst b/manual/sphinx/user_docs/bout_options.rst index 85a8a17d59..f63e65b6df 100644 --- a/manual/sphinx/user_docs/bout_options.rst +++ b/manual/sphinx/user_docs/bout_options.rst @@ -139,10 +139,10 @@ Boolean expressions Boolean values must be "true", "false", "True", "False", "1" or "0". All lowercase ("true"/"false") is preferred, but the uppercase versions are allowed to support Python string conversions. Booleans -can be combined into expressions using binary operators `&` (logical -AND), `|` (logical OR), and unary operator `!` (logical NOT). For -example "true & false" evaluates to `false`; "!false" evaluates to -`true`. Like real values and integers, boolean expressions can refer +can be combined into expressions using binary operators ``&`` (logical +AND), ``|`` (logical OR), and unary operator ``!`` (logical NOT). For +example "true & false" evaluates to ``false``; "!false" evaluates to +``true``. Like real values and integers, boolean expressions can refer to other variables: .. code-block:: cfg @@ -150,8 +150,8 @@ to other variables: switch = true other_switch = !switch -Boolean expressions can be formed by comparing real values using -`>` and `<` comparison operators: +Boolean expressions can be formed by comparing real values using ``>`` +and ``<`` comparison operators: .. code-block:: cfg @@ -160,15 +160,15 @@ Boolean expressions can be formed by comparing real values using is_false = value < 2 .. note:: - Previous BOUT++ versions (v5.1.0 and earlier) were case - insensitive when reading boolean values, so would read "True" or - "yEs" as `true`, and "False" or "No" as `false`. These earlier - versions did not allow boolean expressions. - -Internally, booleans are evaluated as real values, with `true` being 1 -and `false` being 0. Logical operators (`&`, `|`, `!`) check that -their left and right arguments are either close to 0 or close to 1 -(like integers, "close to" is within 1e-3). + Previous BOUT++ versions (v5.1.0 and earlier) were case insensitive + when reading boolean values, so would read "True" or "yEs" as + ``true``, and "False" or "No" as ``false``. These earlier versions + did not allow boolean expressions. + +Internally, booleans are evaluated as real values, with ``true`` being +1 and ``false`` being 0. Logical operators (``&``, ``|``, ``!``) check +that their left and right arguments are either close to 0 or close to +1 (like integers, "close to" is within 1e-3). Special symbols in Option names ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ From 0bc5206ed28aa41ab43e9f32707aa524318386eb Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 22 Jan 2024 15:07:54 +0000 Subject: [PATCH 0030/1827] Docs: Update code layout Prefer linking to sections rather than specific files, as these have changed a fair bit --- manual/sphinx/developer_docs/code_layout.rst | 294 ++++--------------- manual/sphinx/developer_docs/mesh.rst | 4 + manual/sphinx/user_docs/variable_init.rst | 2 + 3 files changed, 58 insertions(+), 242 deletions(-) diff --git a/manual/sphinx/developer_docs/code_layout.rst b/manual/sphinx/developer_docs/code_layout.rst index 7605355d7f..e79245cde3 100644 --- a/manual/sphinx/developer_docs/code_layout.rst +++ b/manual/sphinx/developer_docs/code_layout.rst @@ -49,275 +49,85 @@ Directories ----------- The source code for the core of BOUT++ is divided into include files -(which can be used in physics models) in ``bout++/include``, and source -code and low-level includes in ``bout++/src``. Many parts of the code -are defined by their interface, and can have multiple different -implementations. An example is the time-integration solvers: many -different implementations are available, some of which use external -libraries, but all have the same interface and can be used +(which can be used in physics models) in ``bout++/include/bout``, and +source code and low-level includes in ``bout++/src``. Many parts of +the code are defined by their interface, and can have multiple +different implementations. An example is the time-integration solvers: +many different implementations are available, some of which use +external libraries, but all have the same interface and can be used interchangeably. This is reflected in the directory structure inside -``bout++/src``. A common pattern is to store individual implementations -of an interface in a subdirectory called ``impls``. +``bout++/src``. A common pattern is to store individual +implementations of an interface in a subdirectory called ``impls``. :: - include/foo.hxx + include/bout/foo.hxx src/.../foo.cxx - src/.../foo_factory.hxx - src/.../foo_factory.cxx src/.../impls/one/one.hxx src/.../impls/one/one.cxx where ``foo.hxx`` defines the interface, ``foo.cxx`` implements common -functions used in several implementations. ``foo_factory`` creates new -implementations, and is the only file which includes all the -implementations. Individual implementations are stored in their own -subdirectories of ``impls``. Components which follow this pattern -include ``fileio`` formats, ``invert/laplace`` and ``invert/parderiv`` -inversion codes, ``mesh``, and ``solver``. +functions used in several implementations. Individual implementations +are stored in their own subdirectories of ``impls``. Components which +follow this pattern include ``invert/laplace`` and ``invert/parderiv`` +inversions, ``mesh``, and ``solver``. -The current source code files are: +The layout of the ``src/`` directory is as follows: -- :doc:`bout++.cxx<../_breathe_autogen/file/bout_09_09_8cxx>`: Main - file which initialises, runs and finalises BOUT++. Currently - contains a `main()` function, though this is being removed shortly. +- :doc:`src/bout++.cxx<../_breathe_autogen/file/bout_09_09_8cxx>`: Main + file which initialises, runs and finalises BOUT++. -- field +- ``src/field`` - - :doc:`field2d.cxx<../_breathe_autogen/file/field2d_8cxx>` - implements the `Field2D` class. This is a scalar field which - varies only in :math:`x` and :math:`y` and is used for things - like metric tensor components and initial profiles. It supplies - lots of overloaded operators and functions on these objects. + - Implementations of "fields" (the scalars `Field2D`, `Field3D`, + `FieldPerp`, and vectors `Vector2D`, `Vector3D`), as well as basic + operations (arithmetic, :ref:`sec-algebraic-ops`, and vector + calculus), and :ref:`initialisation `. - - :doc:`field3d.cxx<../_breathe_autogen/file/field3d_8cxx>` - implements the `Field3D` class, which varies in :math:`x`, - :math:`y` and :math:`z`. Since these handle a lot more memory - than Field2D objects, the memory management is more complicated - and includes reference counting. See section - :ref:`sec-memorymanage` for more details. +- ``src/invert`` - - :doc:`field_data.cxx<../_breathe_autogen/file/field__data_8cxx>` - Implements some functions in the `FieldData` class. This is a - mainly pure virtual interface class which is inherited by - `Field2D` and `Field3D`. + - Implementations of different inverse operators, including Fourier + transforms (via an interface to `FFTW `_, see + `bout::fft::rfft` and `bout::fft::irfft`). - - :doc:`fieldperp.cxx<../_breathe_autogen/file/fieldperp_8cxx>` - implements a `FieldPerp` class to store slices perpendicular to - the magnetic field i.e. they are a function of :math:`x` and - :math:`z` only. This is mainly used for Laplacian inversion - routines, and needs to be integrated with the other fields - better. +- ``src/invert/laplace`` - - :doc:`initialprofiles.cxx<../_breathe_autogen/file/initialprofiles_8cxx>` - routines to set the initial values of fields when a simulation - first starts. Reads settings from the option file based on the name - of the variable. + - Implementations of some inverse generalised Laplacian operator + variations, where some directions may be constant. See + :ref:`sec-laplacian` for more details. - - :doc:`vecops.cxx<../_breathe_autogen/file/vecops_8cxx>` a - collection of function to operate on vectors. Contains things - like ``Grad``, ``Div`` and ``Curl``, and uses a combination of - field differential operators (in - :doc:`difops.cxx<../_breathe_autogen/file/difops_8cxx>`) and - metric tensor components (in `Mesh`). +- ``src/invert/parderiv`` - - :doc:`vector2d.cxx<../_breathe_autogen/file/vector2d_8cxx>` - implements the `Vector2D` class, which uses a `Field2D` object - for each of its 3 components. Overloads operators to supply - things like dot and cross products. + - Inversion of parallel derivatives, intended for use in + preconditioners. See `InvertPar` - - :doc:`vector3d.cxx<../_breathe_autogen/file/vector3d_8cxx>` - implements `Vector3D` by using a `Field3D` - object for each component. +- ``src/invert/pardiv`` - - :doc:`where.cxx<../_breathe_autogen/file/where_8cxx>` supplies - functions for choosing between values based on selection - criteria. + - Inversion of parallel divergence, intended for use in + :ref:`sec-nonlocal-heatflux`. -- invert +- ``src/mesh`` - - :doc:`fft_fftw.cxx<../_breathe_autogen/file/fft__fftw_8cxx>` - implements the :doc:`fft.hxx<../_breathe_autogen/file/fft_8hxx>` - interface by calling the Fastest Fourier Transform in the West - (FFTW) library. + - Implementations of low-level numerical routines. This includes + things like the :ref:`inter-process communications + `, :ref:`boundary conditions `, + :ref:`derivative operators `, interpolation, the + coordinate system (including :ref:`sec-parallel-transforms`), and + the :ref:`sec-mesh` itself. -- invert / laplace +- ``src/physics`` - - :doc:`invert_laplace.cxx<../_breathe_autogen/file/invert__laplace_8cxx>` uses Fourier - decomposition in :math:`z` combined with tri- and band-diagonal - solvers in :math:`x` to solve Laplacian problems. + - This contains some specialised physics operators and routines, + such as gyro-averaging and :ref:`sec-nonlocal-heatflux`. - - impls +- ``src/solver`` - - serial\_tri + - Implementations of :ref:`time integration solvers + ` - - :doc:`serial_tri.hxx<../_breathe_autogen/file/serial__tri_8hxx>` +- ``src/sys`` - - :doc:`serial_tri.cxx<../_breathe_autogen/file/serial__tri_8cxx>` - - - serial\_band - - - :doc:`serial_band.hxx<../_breathe_autogen/file/serial__band_8hxx>` - - - :doc:`serial_band.cxx<../_breathe_autogen/file/serial__band_8cxx>` - - - spt - - - :doc:`spt.hxx<../_breathe_autogen/file/spt_8hxx>` - - - :doc:`spt.cxx<../_breathe_autogen/file/spt_8cxx>` - -- invert / parderiv - - - - :doc:`invert_parderiv.cxx<../_breathe_autogen/file/invert__parderiv_8cxx>` - inverts a problem involving only parallel :math:`y` - derivatives. Intended for use in some preconditioners. - - - impls - - - cyclic - - - :doc:`cyclic.cxx<../_breathe_autogen/file/cyclic_8cxx>` - - - :doc:`cyclic.hxx<../_breathe_autogen/file/cyclic_8hxx>` - -- :doc:`lapack_routines.cxx<../_breathe_autogen/file/lapack__routines_8cxx>` supplies an - interface to the LAPACK linear solvers, which are used by the - ``invert_laplace`` routines. - -- mesh - - - :doc:`boundary_factory.cxx<../_breathe_autogen/file/boundary__factory_8cxx>` creates boundary - condition operators which can then be applied to - fields. Described in section :ref:`sec-BoundaryFactory`. - - - :doc:`boundary_region.cxx<../_breathe_autogen/file/boundary__region_8cxx>` implements a way - to describe and iterate over boundary regions. Created by the - mesh, and then used by boundary conditions. See - section :ref:`sec-BoundaryRegion` for more details. - - - :doc:`boundary_standard.cxx<../_breathe_autogen/file/boundary__standard_8cxx>` implements some - standard boundary operations and modifiers such as ``Neumann`` - and ``Dirichlet``. - - - :doc:`difops.cxx<../_breathe_autogen/file/difops_8cxx>` is a - collection of differential operators on scalar fields. It uses - the differential methods in :doc:`derivs.cxx<../_breathe_autogen/file/derivs_8cxx>` and the metric tensor - components in `Mesh` to compute operators. - - - :doc:`interpolation.cxx<../_breathe_autogen/file/interpolation_8cxx>` contains functions - for interpolating fields - - - :doc:`mesh.cxx<../_breathe_autogen/file/mesh_8cxx>` is the base - class for the `Mesh` object. Contains routines useful - for all `Mesh` implementations. - - - impls - - - bout - - - :doc:`boutmesh.cxx<../_breathe_autogen/file/boutmesh_8cxx>` - implements a mesh interface which is compatible with BOUT - grid files. - - - :doc:`boutmesh.hxx<../_breathe_autogen/file/boutmesh_8hxx>` - -- physics - - - :doc:`gyro_average.cxx<../_breathe_autogen/file/gyro__average_8cxx>` - gyro-averaging operators - - - :doc:`smoothing.cxx<../_breathe_autogen/file/smoothing_8cxx>` - provides smoothing routines on scalar fields - - - :doc:`sourcex.cxx<../_breathe_autogen/file/sourcex_8cxx>` contains - some useful routines for creating sources and sinks in physics - equations. - -- solver - - - :doc:`solver.cxx<../_breathe_autogen/file/solver_8cxx>` is the - interface for all solvers - - - impls - - - cvode - - - :doc:`cvode.cxx<../_breathe_autogen/file/cvode_8cxx>` is the - implementation of `Solver` which interfaces with - the SUNDIALS CVODE library. - - - :doc:`cvode.hxx<../_breathe_autogen/file/cvode_8hxx>` - - - ida - - - :doc:`ida.cxx<../_breathe_autogen/file/ida_8cxx>` is the - implementation which interfaces with the SUNDIALS IDA - library - - - :doc:`ida.hxx<../_breathe_autogen/file/ida_8hxx>` - - - petsc - - - :doc:`petsc.cxx<../_breathe_autogen/file/petsc_8cxx>` is the - interface to the PETSc time integration routines - - - :doc:`petsc.hxx<../_breathe_autogen/file/petsc_8hxx>` - - - pvode - - - :doc:`pvode.cxx<../_breathe_autogen/file/pvode_8cxx>` - interfaces with the 1998 (pre-SUNDIALS) version of PVODE - (which became CVODE). - - - :doc:`pvode.hxx<../_breathe_autogen/file/pvode_8hxx>` - -- sys - - - :doc:`boutcomm.cxx<../_breathe_autogen/file/boutcomm_8cxx>` - - - :doc:`boutexception.cxx<../_breathe_autogen/file/boutexception_8cxx>` - is an exception class which are used for error handling - - - :doc:`derivs.cxx<../_breathe_autogen/file/derivs_8cxx>` contains - basic derivative methods such as upwinding, central difference - and WENO methods. These are then used by - :doc:`difops.cxx<../_breathe_autogen/file/difops_8cxx>`. Details are - given in section :ref:`sec-diffops`. - - - :doc:`msg_stack.cxx<../_breathe_autogen/file/msg__stack_8cxx>` is - part of the error handling system. It maintains a stack of - messages which can be pushed onto the stack at the start of a - function, then removed (popped) at the end. If an error occurs or - a segmentation fault is caught then this stack is printed out and - can help to find errors. - - - :doc:`options.cxx<../_breathe_autogen/file/options_8cxx>` provides - an interface to the BOUT.inp option file and the command-line - options. - - - :doc:`optionsreader.cxx<../_breathe_autogen/file/optionsreader_8cxx>` - - - :doc:`output.cxx<../_breathe_autogen/file/output_8cxx>` - - - :doc:`range.cxx<../_breathe_autogen/file/range_8cxx>` Provides the - RangeIterator class, used to iterate over a set of - ranges. Described in section :ref:`sec-rangeiterator` - - - :doc:`timer.cxx<../_breathe_autogen/file/timer_8cxx>` a class for - timing parts of the code like communications and file - I/O. Described in section :ref:`sec-timerclass` - - - :doc:`utils.cxx<../_breathe_autogen/file/utils_8cxx>` contains - miscellaneous small useful routines such as allocating and - freeing arrays. - - - options - - - :doc:`optionparser.hxx<../_breathe_autogen/file/optionparser_8hxx>` - - - :doc:`options_ini.cxx<../_breathe_autogen/file/options__ini_8cxx>` - - - :doc:`options_ini.hxx<../_breathe_autogen/file/options__ini_8hxx>` + - General purpose utilities used throughout the library, such as + `BoutException`, wrappers for C libraries like ``PETSc`` and + ``HYPRE``, screen and file input and output. diff --git a/manual/sphinx/developer_docs/mesh.rst b/manual/sphinx/developer_docs/mesh.rst index ff3e40d4b0..c592773324 100644 --- a/manual/sphinx/developer_docs/mesh.rst +++ b/manual/sphinx/developer_docs/mesh.rst @@ -1,3 +1,5 @@ +.. _sec-mesh: + Mesh ==== @@ -99,6 +101,8 @@ it:: To read 2D and 3D fields, the branch-cuts need to be taken into account. +.. _sec-communications: + Communications -------------- diff --git a/manual/sphinx/user_docs/variable_init.rst b/manual/sphinx/user_docs/variable_init.rst index 4ac13e0ead..ebe989c233 100644 --- a/manual/sphinx/user_docs/variable_init.rst +++ b/manual/sphinx/user_docs/variable_init.rst @@ -1,3 +1,5 @@ +.. _sec-variable-init: + Variable initialisation ======================= From 297b28430988b871815947071e4bdd7ae461eca2 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 22 Jan 2024 15:08:31 +0000 Subject: [PATCH 0031/1827] Docs: Fix underline too short --- manual/sphinx/user_docs/installing.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/manual/sphinx/user_docs/installing.rst b/manual/sphinx/user_docs/installing.rst index eb155909bf..22fc3e724f 100644 --- a/manual/sphinx/user_docs/installing.rst +++ b/manual/sphinx/user_docs/installing.rst @@ -397,7 +397,7 @@ installation for that. Working with an active ``conda`` environment -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ When ``conda`` is used, it installs separate versions of several libraries. These can cause warnings or even failures when linking BOUT++ executables. There are From e6bc5c5d0e85b5e895b124b066bbf62aed6c45b1 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 22 Jan 2024 15:12:09 +0000 Subject: [PATCH 0032/1827] Docs: Don't include boutdata/boututils in docs --- manual/sphinx/conf.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/manual/sphinx/conf.py b/manual/sphinx/conf.py index 29c0985841..d8ab76cf23 100755 --- a/manual/sphinx/conf.py +++ b/manual/sphinx/conf.py @@ -31,7 +31,8 @@ import subprocess import sys -sys.path.append("../../tools/pylib") +sys.path.append("../../tools/pylib/boutpp") +sys.path.append("../../tools/pylib/boutconfig") # Are we running on readthedocs? on_readthedocs = os.environ.get("READTHEDOCS") == "True" From 41da8a816f3217c608ca3739330f9f3fc9186902 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 22 Jan 2024 15:35:11 +0000 Subject: [PATCH 0033/1827] Docs: Fix reference link --- manual/sphinx/user_docs/physics_models.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/manual/sphinx/user_docs/physics_models.rst b/manual/sphinx/user_docs/physics_models.rst index 56e3b719b4..fa62541f3d 100644 --- a/manual/sphinx/user_docs/physics_models.rst +++ b/manual/sphinx/user_docs/physics_models.rst @@ -28,7 +28,7 @@ either physics models (e.g. ``test-delp2`` and Building Physics Models ----------------------- -After building the library (see :ref:`sec-cmake`), you can build a +After building the library (see :ref:`sec-config-bout`), you can build a physics model in several different ways. For the bundled examples, perhaps the easiest is to build it directly From 96723d088369ad73e8e74329a860a1775395da6c Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 22 Jan 2024 15:47:00 +0000 Subject: [PATCH 0034/1827] Docs: Fix bullet points (require blank line before start) --- manual/sphinx/user_docs/installing.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/manual/sphinx/user_docs/installing.rst b/manual/sphinx/user_docs/installing.rst index 22fc3e724f..4f47255d2e 100644 --- a/manual/sphinx/user_docs/installing.rst +++ b/manual/sphinx/user_docs/installing.rst @@ -402,14 +402,17 @@ Working with an active ``conda`` environment When ``conda`` is used, it installs separate versions of several libraries. These can cause warnings or even failures when linking BOUT++ executables. There are several alternatives to deal with this problem: + * The simplest but least convenient option is to use ``conda deactivate`` before configuring, compiling, or running any BOUT++ program. + * You might sometimes want to link to the conda-installed libraries. This is probably not ideal for production runs on an HPC system (as conda downloads binary packages that will not be optimized for specific hardware), but can be a simple way to get packages for testing or on a personal computer. In this case just keep your ``conda`` environment active, and with luck the libraries should be picked up by the standard search mechanisms. + * In case you do want a fully optimized and as-stable-as-possible build for production runs, it is probably best not to depend on any conda packages for compiling or running BOUT++ executables (restrict ``conda`` to providing Python From 82c0708b6be5a2db608b5a1b884c76178404c89a Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 22 Jan 2024 15:47:46 +0000 Subject: [PATCH 0035/1827] Docs: Fix wrong start of literal block --- manual/sphinx/user_docs/running_bout.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/manual/sphinx/user_docs/running_bout.rst b/manual/sphinx/user_docs/running_bout.rst index 59c5b53828..5ee7a9bfb3 100644 --- a/manual/sphinx/user_docs/running_bout.rst +++ b/manual/sphinx/user_docs/running_bout.rst @@ -457,7 +457,7 @@ values loaded, the solver can be started:: 3d fields = 2, 2d fields = 0 neq=100, local_N=100 This last line gives the number of equations being evolved (in this case -100), and the number of these on this processor (here 100).:: +100), and the number of these on this processor (here 100). The absolute and relative tolerances come next:: From e1591bcda65b93788e42392c470461a690307726 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 22 Jan 2024 15:55:24 +0000 Subject: [PATCH 0036/1827] Docs: Fix snippet for Scorep --- .../developer_docs/performance_profiling.rst | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/manual/sphinx/developer_docs/performance_profiling.rst b/manual/sphinx/developer_docs/performance_profiling.rst index d93c1e904b..de14d695da 100644 --- a/manual/sphinx/developer_docs/performance_profiling.rst +++ b/manual/sphinx/developer_docs/performance_profiling.rst @@ -85,9 +85,7 @@ Configure and build Configure with ``-BOUT_USE_SCOREP`` to enable Scorep instrumentation, then build as normal. This option can be combined with other options, but it is usually desirable to profile the optimized code, configuring -with the flags ````. Build the code with ``make`` as normal. - -With CMake: +with the flags: .. code-block:: bash @@ -95,14 +93,13 @@ With CMake: -DCMAKE_C_COMPILER=scorep-mpicc \ -DCMAKE_CXX_COMPILER=scorep-mpicxx \ -DCMAKE_CXX_FLAGS=-O3 -DCHECK=0 \ + -DCMAKE_BUILD_TYPE=RelWithDebInfo \ -This will turn off the instrumentation during the configure -step. Please be aware that if you change ``CMakeLists.txt``, CMake -will try to automatically reconfigure the build, which the Score-P -wrappers interfere with. In this case you will need to restart the -configure step from scratch (i.e. remove the build directory and start -again). +``SCOREP_WRAPPER=off`` is used turn off the instrumentation during the +configure step. You probably want to use a fresh build directory to +prevent any issues from cached results when setting the compilers like +this. Run and analysis ~~~~~~~~~~~~~~~~ From dcb33f44a215a267f98a8dcb69a3056e9a4dac5a Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 22 Jan 2024 15:58:41 +0000 Subject: [PATCH 0037/1827] Docs: Fix unknown role --- manual/sphinx/developer_docs/file_io.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/manual/sphinx/developer_docs/file_io.rst b/manual/sphinx/developer_docs/file_io.rst index aece8cfb1e..0d13e70cfc 100644 --- a/manual/sphinx/developer_docs/file_io.rst +++ b/manual/sphinx/developer_docs/file_io.rst @@ -114,5 +114,5 @@ in the first output file, which `collect` uses to get its type and dimensions. .. [2] Actually, the C++ I/O code should work fine even if a `FieldPerp` object is defined with different y-indices on different processors. This may be useful for diagnostic or debugging purposes. However, Python routines like `collect` and - :py:`boutdata.restart.redistribute` will fail because they find inconsistent + `boutdata.restart.redistribute` will fail because they find inconsistent `yindex_global` values. From aa4a6553394f9420564d12ff8df681a8ba6d0618 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 22 Jan 2024 16:07:44 +0000 Subject: [PATCH 0038/1827] Docs: Fix syntax for footnote --- manual/sphinx/user_docs/boundary_options.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/manual/sphinx/user_docs/boundary_options.rst b/manual/sphinx/user_docs/boundary_options.rst index 57c6658891..d5e1b9199c 100644 --- a/manual/sphinx/user_docs/boundary_options.rst +++ b/manual/sphinx/user_docs/boundary_options.rst @@ -108,7 +108,7 @@ Boundary conditions for non-orthogonal grids If non-orthogonal grids are used (meaning that the x- and y-directions are not orthogonal, so ``g12 != 0.``), then corner cells may be required. The boundary conditions are applied -in corner cells[#disablecorners]_ by applying the y-boundary condition using x-boundary +in corner cells [#disablecorners]_ by applying the y-boundary condition using x-boundary values. This requires that x-boundary conditions are applied before y-boundary conditions. The ordering is taken care of by the methods described in this section, but also needs to be respected by any custom boundary conditions in user code (e.g. sheath boundary From 6640e940e350c3636d83026f2ee2a60322fb9440 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 22 Jan 2024 17:00:55 +0000 Subject: [PATCH 0039/1827] Fix a couple of "most vexing parse" that doxygen doesn't like --- include/bout/dcomplex.hxx | 4 +++- include/bout/sys/uuid.h | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/include/bout/dcomplex.hxx b/include/bout/dcomplex.hxx index 569b5f2c13..3bb360b055 100644 --- a/include/bout/dcomplex.hxx +++ b/include/bout/dcomplex.hxx @@ -37,7 +37,9 @@ using dcomplex = std::complex; -const dcomplex Im(0, 1); // 1i +/// Imaginary unit, ``1i`` +//NOLINTNEXTLINE(readability-identifier-length) +constexpr dcomplex Im{0, 1}; /// Complex type for passing data to/from FORTRAN struct fcmplx { diff --git a/include/bout/sys/uuid.h b/include/bout/sys/uuid.h index d043ea780e..edc252bb44 100644 --- a/include/bout/sys/uuid.h +++ b/include/bout/sys/uuid.h @@ -261,7 +261,7 @@ class sha1 { size_t m_byteCount; }; -static std::mt19937 clock_gen(std::random_device{}()); +static std::mt19937 clock_gen{std::random_device{}()}; static std::uniform_int_distribution clock_dis{-32768, 32767}; static std::atomic_short clock_sequence{clock_dis(clock_gen)}; } // namespace detail From 534fe4db38604369f111e1f2020f93e3f005edd6 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 22 Jan 2024 17:07:50 +0000 Subject: [PATCH 0040/1827] Docs: Don't keep separate doxygen config for readthedocs --- manual/doxygen/Doxyfile | 1390 +-------------------------- manual/doxygen/Doxyfile_readthedocs | 1119 --------------------- manual/sphinx/conf.py | 5 +- 3 files changed, 19 insertions(+), 2495 deletions(-) delete mode 100644 manual/doxygen/Doxyfile_readthedocs diff --git a/manual/doxygen/Doxyfile b/manual/doxygen/Doxyfile index 216ef60915..e27df586ec 100644 --- a/manual/doxygen/Doxyfile +++ b/manual/doxygen/Doxyfile @@ -246,34 +246,6 @@ ALIASES = TCL_SUBST = -# Set the OPTIMIZE_OUTPUT_FOR_C tag to YES if your project consists of C sources -# only. Doxygen will then generate output that is more tailored for C. For -# instance, some of the names that are used will be different. The list of all -# members will be omitted, etc. -# The default value is: NO. - -OPTIMIZE_OUTPUT_FOR_C = NO - -# Set the OPTIMIZE_OUTPUT_JAVA tag to YES if your project consists of Java or -# Python sources only. Doxygen will then generate output that is more tailored -# for that language. For instance, namespaces will be presented as packages, -# qualified scopes will look different, etc. -# The default value is: NO. - -OPTIMIZE_OUTPUT_JAVA = NO - -# Set the OPTIMIZE_FOR_FORTRAN tag to YES if your project consists of Fortran -# sources. Doxygen will then generate output that is tailored for Fortran. -# The default value is: NO. - -OPTIMIZE_FOR_FORTRAN = NO - -# Set the OPTIMIZE_OUTPUT_VHDL tag to YES if your project consists of VHDL -# sources. Doxygen will then generate output that is tailored for VHDL. -# The default value is: NO. - -OPTIMIZE_OUTPUT_VHDL = NO - # Doxygen selects the parser to use depending on the extension of the files it # parses. With this tag you can assign which parser to use for a given # extension. Doxygen has a built-in mapping, but you can override or extend it @@ -328,21 +300,7 @@ AUTOLINK_SUPPORT = YES # diagrams that involve STL classes more complete and accurate. # The default value is: NO. -BUILTIN_STL_SUPPORT = NO - -# If you use Microsoft's C++/CLI language, you should set this option to YES to -# enable parsing support. -# The default value is: NO. - -CPP_CLI_SUPPORT = NO - -# Set the SIP_SUPPORT tag to YES if your project consists of sip (see: -# http://www.riverbankcomputing.co.uk/software/sip/intro) sources only. Doxygen -# will parse them like normal C++ but will assume all classes use public instead -# of private inheritance when no explicit protection keyword is present. -# The default value is: NO. - -SIP_SUPPORT = NO +BUILTIN_STL_SUPPORT = YES # For Microsoft's IDL there are propget and propput attributes to indicate # getter and setter methods for a property. Setting this option to YES will make @@ -352,7 +310,7 @@ SIP_SUPPORT = NO # should set this option to NO. # The default value is: YES. -IDL_PROPERTY_SUPPORT = YES +IDL_PROPERTY_SUPPORT = NO # If member grouping is used in the documentation and the DISTRIBUTE_GROUP_DOC # tag is set to YES then doxygen will reuse the documentation of the first @@ -435,13 +393,13 @@ LOOKUP_CACHE_SIZE = 0 # normally produced when WARNINGS is set to YES. # The default value is: NO. -EXTRACT_ALL = YES +EXTRACT_ALL = NO # If the EXTRACT_PRIVATE tag is set to YES, all private members of a class will # be included in the documentation. # The default value is: NO. -EXTRACT_PRIVATE = YES +EXTRACT_PRIVATE = NO # If the EXTRACT_PACKAGE tag is set to YES, all members with package or internal # scope will be included in the documentation. @@ -463,14 +421,6 @@ EXTRACT_STATIC = YES EXTRACT_LOCAL_CLASSES = YES -# This flag is only useful for Objective-C code. If set to YES, local methods, -# which are defined in the implementation section but not in the interface are -# included in the documentation. If set to NO, only methods in the interface are -# included. -# The default value is: NO. - -EXTRACT_LOCAL_METHODS = YES - # If this flag is set to YES, the members of anonymous namespaces will be # extracted and appear in the documentation as a namespace called # 'anonymous_namespace{file}', where file will be replaced with the base name of @@ -508,7 +458,7 @@ HIDE_FRIEND_COMPOUNDS = NO # blocks will be appended to the function's detailed documentation block. # The default value is: NO. -HIDE_IN_BODY_DOCS = NO +HIDE_IN_BODY_DOCS = YES # The INTERNAL_DOCS tag determines if documentation that is typed after a # \internal command is included. If the tag is set to NO then the documentation @@ -590,7 +540,7 @@ SORT_BRIEF_DOCS = NO # detailed member documentation. # The default value is: NO. -SORT_MEMBERS_CTORS_1ST = NO +SORT_MEMBERS_CTORS_1ST = YES # If the SORT_GROUP_NAMES tag is set to YES then doxygen will sort the hierarchy # of group names into alphabetical order. If set to NO the group names will @@ -623,19 +573,19 @@ STRICT_PROTO_MATCHING = NO # list. This list is created by putting \todo commands in the documentation. # The default value is: YES. -GENERATE_TODOLIST = YES +GENERATE_TODOLIST = NO # The GENERATE_TESTLIST tag can be used to enable (YES) or disable (NO) the test # list. This list is created by putting \test commands in the documentation. # The default value is: YES. -GENERATE_TESTLIST = YES +GENERATE_TESTLIST = NO # The GENERATE_BUGLIST tag can be used to enable (YES) or disable (NO) the bug # list. This list is created by putting \bug commands in the documentation. # The default value is: YES. -GENERATE_BUGLIST = YES +GENERATE_BUGLIST = NO # The GENERATE_DEPRECATEDLIST tag can be used to enable (YES) or disable (NO) # the deprecated list. This list is created by putting \deprecated commands in @@ -673,7 +623,7 @@ SHOW_USED_FILES = YES # (if specified). # The default value is: YES. -SHOW_FILES = YES +SHOW_FILES = NO # Set the SHOW_NAMESPACES tag to NO to disable the generation of the Namespaces # page. This will remove the Namespaces entry from the Quick Index and from the @@ -791,7 +741,7 @@ WARN_LOGFILE = # Note: If this tag is empty the current directory is searched. INPUT = ../../include/bout \ - ../../src/ + ../../src # This tag can be used to specify the character encoding of the source files # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses @@ -858,99 +808,6 @@ EXCLUDE_SYMLINKS = NO EXCLUDE_PATTERNS = -# The EXCLUDE_SYMBOLS tag can be used to specify one or more symbol names -# (namespaces, classes, functions, etc.) that should be excluded from the -# output. The symbol name can be a fully qualified name, a word, or if the -# wildcard * is used, a substring. Examples: ANamespace, AClass, -# AClass::ANamespace, ANamespace::*Test -# -# Note that the wildcards are matched against the file with absolute path, so to -# exclude all test directories use the pattern */test/* - -EXCLUDE_SYMBOLS = - -# The EXAMPLE_PATH tag can be used to specify one or more files or directories -# that contain example code fragments that are included (see the \include -# command). - -EXAMPLE_PATH = - -# If the value of the EXAMPLE_PATH tag contains directories, you can use the -# EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp and -# *.h) to filter out the source-files in the directories. If left blank all -# files are included. - -EXAMPLE_PATTERNS = * - -# If the EXAMPLE_RECURSIVE tag is set to YES then subdirectories will be -# searched for input files to be used with the \include or \dontinclude commands -# irrespective of the value of the RECURSIVE tag. -# The default value is: NO. - -EXAMPLE_RECURSIVE = NO - -# The IMAGE_PATH tag can be used to specify one or more files or directories -# that contain images that are to be included in the documentation (see the -# \image command). - -IMAGE_PATH = - -# The INPUT_FILTER tag can be used to specify a program that doxygen should -# invoke to filter for each input file. Doxygen will invoke the filter program -# by executing (via popen()) the command: -# -# -# -# where is the value of the INPUT_FILTER tag, and is the -# name of an input file. Doxygen will then use the output that the filter -# program writes to standard output. If FILTER_PATTERNS is specified, this tag -# will be ignored. -# -# Note that the filter must not add or remove lines; it is applied before the -# code is scanned, but not when the output code is generated. If lines are added -# or removed, the anchors will not be placed correctly. -# -# Note that for custom extensions or not directly supported extensions you also -# need to set EXTENSION_MAPPING for the extension otherwise the files are not -# properly processed by doxygen. - -INPUT_FILTER = - -# The FILTER_PATTERNS tag can be used to specify filters on a per file pattern -# basis. Doxygen will compare the file name with each pattern and apply the -# filter if there is a match. The filters are a list of the form: pattern=filter -# (like *.cpp=my_cpp_filter). See INPUT_FILTER for further information on how -# filters are used. If the FILTER_PATTERNS tag is empty or if none of the -# patterns match the file name, INPUT_FILTER is applied. -# -# Note that for custom extensions or not directly supported extensions you also -# need to set EXTENSION_MAPPING for the extension otherwise the files are not -# properly processed by doxygen. - -FILTER_PATTERNS = - -# If the FILTER_SOURCE_FILES tag is set to YES, the input filter (if set using -# INPUT_FILTER) will also be used to filter the input files that are used for -# producing the source files to browse (i.e. when SOURCE_BROWSER is set to YES). -# The default value is: NO. - -FILTER_SOURCE_FILES = NO - -# The FILTER_SOURCE_PATTERNS tag can be used to specify source filters per file -# pattern. A pattern will override the setting for FILTER_PATTERN (if any) and -# it is also possible to disable source filtering for a specific pattern using -# *.ext= (so without naming a filter). -# This tag requires that the tag FILTER_SOURCE_FILES is set to YES. - -FILTER_SOURCE_PATTERNS = - -# If the USE_MDFILE_AS_MAINPAGE tag refers to the name of a markdown file that -# is part of the input, its contents will be placed on the main page -# (index.html). This can be useful if you have a project on for instance GitHub -# and want to reuse the introduction page also for the doxygen output. - -USE_MDFILE_AS_MAINPAGE = - #--------------------------------------------------------------------------- # Configuration options related to source browsing #--------------------------------------------------------------------------- @@ -1005,7 +862,7 @@ REFERENCES_LINK_SOURCE = YES # The default value is: YES. # This tag requires that the tag SOURCE_BROWSER is set to YES. -SOURCE_TOOLTIPS = YES +SOURCE_TOOLTIPS = NO # If the USE_HTAGS tag is set to YES then the references to source code will # point to the HTML generated by the htags(1) tool instead of doxygen built-in @@ -1048,21 +905,6 @@ VERBATIM_HEADERS = NO ALPHABETICAL_INDEX = NO -# The COLS_IN_ALPHA_INDEX tag can be used to specify the number of columns in -# which the alphabetical index list will be split. -# Minimum value: 1, maximum value: 20, default value: 5. -# This tag requires that the tag ALPHABETICAL_INDEX is set to YES. - -COLS_IN_ALPHA_INDEX = 5 - -# In case all classes in a project start with a common prefix, all classes will -# be put under the same header in the alphabetical index. The IGNORE_PREFIX tag -# can be used to specify a prefix (or a list of prefixes) that should be ignored -# while generating the index headers. -# This tag requires that the tag ALPHABETICAL_INDEX is set to YES. - -IGNORE_PREFIX = - #--------------------------------------------------------------------------- # Configuration options related to the HTML output #--------------------------------------------------------------------------- @@ -1070,542 +912,7 @@ IGNORE_PREFIX = # If the GENERATE_HTML tag is set to YES, doxygen will generate HTML output # The default value is: YES. -GENERATE_HTML = YES - -# The HTML_OUTPUT tag is used to specify where the HTML docs will be put. If a -# relative path is entered the value of OUTPUT_DIRECTORY will be put in front of -# it. -# The default directory is: html. -# This tag requires that the tag GENERATE_HTML is set to YES. - -HTML_OUTPUT = html - -# The HTML_FILE_EXTENSION tag can be used to specify the file extension for each -# generated HTML page (for example: .htm, .php, .asp). -# The default value is: .html. -# This tag requires that the tag GENERATE_HTML is set to YES. - -HTML_FILE_EXTENSION = .html - -# The HTML_HEADER tag can be used to specify a user-defined HTML header file for -# each generated HTML page. If the tag is left blank doxygen will generate a -# standard header. -# -# To get valid HTML the header file that includes any scripts and style sheets -# that doxygen needs, which is dependent on the configuration options used (e.g. -# the setting GENERATE_TREEVIEW). It is highly recommended to start with a -# default header using -# doxygen -w html new_header.html new_footer.html new_stylesheet.css -# YourConfigFile -# and then modify the file new_header.html. See also section "Doxygen usage" -# for information on how to generate the default header that doxygen normally -# uses. -# Note: The header is subject to change so you typically have to regenerate the -# default header when upgrading to a newer version of doxygen. For a description -# of the possible markers and block names see the documentation. -# This tag requires that the tag GENERATE_HTML is set to YES. - -HTML_HEADER = - -# The HTML_FOOTER tag can be used to specify a user-defined HTML footer for each -# generated HTML page. If the tag is left blank doxygen will generate a standard -# footer. See HTML_HEADER for more information on how to generate a default -# footer and what special commands can be used inside the footer. See also -# section "Doxygen usage" for information on how to generate the default footer -# that doxygen normally uses. -# This tag requires that the tag GENERATE_HTML is set to YES. - -HTML_FOOTER = - -# The HTML_STYLESHEET tag can be used to specify a user-defined cascading style -# sheet that is used by each HTML page. It can be used to fine-tune the look of -# the HTML output. If left blank doxygen will generate a default style sheet. -# See also section "Doxygen usage" for information on how to generate the style -# sheet that doxygen normally uses. -# Note: It is recommended to use HTML_EXTRA_STYLESHEET instead of this tag, as -# it is more robust and this tag (HTML_STYLESHEET) will in the future become -# obsolete. -# This tag requires that the tag GENERATE_HTML is set to YES. - -HTML_STYLESHEET = - -# The HTML_EXTRA_STYLESHEET tag can be used to specify additional user-defined -# cascading style sheets that are included after the standard style sheets -# created by doxygen. Using this option one can overrule certain style aspects. -# This is preferred over using HTML_STYLESHEET since it does not replace the -# standard style sheet and is therefore more robust against future updates. -# Doxygen will copy the style sheet files to the output directory. -# Note: The order of the extra style sheet files is of importance (e.g. the last -# style sheet in the list overrules the setting of the previous ones in the -# list). For an example see the documentation. -# This tag requires that the tag GENERATE_HTML is set to YES. - -HTML_EXTRA_STYLESHEET = - -# The HTML_EXTRA_FILES tag can be used to specify one or more extra images or -# other source files which should be copied to the HTML output directory. Note -# that these files will be copied to the base HTML output directory. Use the -# $relpath^ marker in the HTML_HEADER and/or HTML_FOOTER files to load these -# files. In the HTML_STYLESHEET file, use the file name only. Also note that the -# files will be copied as-is; there are no commands or markers available. -# This tag requires that the tag GENERATE_HTML is set to YES. - -HTML_EXTRA_FILES = - -# The HTML_COLORSTYLE_HUE tag controls the color of the HTML output. Doxygen -# will adjust the colors in the style sheet and background images according to -# this color. Hue is specified as an angle on a colorwheel, see -# http://en.wikipedia.org/wiki/Hue for more information. For instance the value -# 0 represents red, 60 is yellow, 120 is green, 180 is cyan, 240 is blue, 300 -# purple, and 360 is red again. -# Minimum value: 0, maximum value: 359, default value: 220. -# This tag requires that the tag GENERATE_HTML is set to YES. - -HTML_COLORSTYLE_HUE = 220 - -# The HTML_COLORSTYLE_SAT tag controls the purity (or saturation) of the colors -# in the HTML output. For a value of 0 the output will use grayscales only. A -# value of 255 will produce the most vivid colors. -# Minimum value: 0, maximum value: 255, default value: 100. -# This tag requires that the tag GENERATE_HTML is set to YES. - -HTML_COLORSTYLE_SAT = 100 - -# The HTML_COLORSTYLE_GAMMA tag controls the gamma correction applied to the -# luminance component of the colors in the HTML output. Values below 100 -# gradually make the output lighter, whereas values above 100 make the output -# darker. The value divided by 100 is the actual gamma applied, so 80 represents -# a gamma of 0.8, The value 220 represents a gamma of 2.2, and 100 does not -# change the gamma. -# Minimum value: 40, maximum value: 240, default value: 80. -# This tag requires that the tag GENERATE_HTML is set to YES. - -HTML_COLORSTYLE_GAMMA = 80 - -# If the HTML_TIMESTAMP tag is set to YES then the footer of each generated HTML -# page will contain the date and time when the page was generated. Setting this -# to YES can help to show when doxygen was last run and thus if the -# documentation is up to date. -# The default value is: NO. -# This tag requires that the tag GENERATE_HTML is set to YES. - -HTML_TIMESTAMP = NO - -# If the HTML_DYNAMIC_SECTIONS tag is set to YES then the generated HTML -# documentation will contain sections that can be hidden and shown after the -# page has loaded. -# The default value is: NO. -# This tag requires that the tag GENERATE_HTML is set to YES. - -HTML_DYNAMIC_SECTIONS = NO - -# With HTML_INDEX_NUM_ENTRIES one can control the preferred number of entries -# shown in the various tree structured indices initially; the user can expand -# and collapse entries dynamically later on. Doxygen will expand the tree to -# such a level that at most the specified number of entries are visible (unless -# a fully collapsed tree already exceeds this amount). So setting the number of -# entries 1 will produce a full collapsed tree by default. 0 is a special value -# representing an infinite number of entries and will result in a full expanded -# tree by default. -# Minimum value: 0, maximum value: 9999, default value: 100. -# This tag requires that the tag GENERATE_HTML is set to YES. - -HTML_INDEX_NUM_ENTRIES = 100 - -# If the GENERATE_DOCSET tag is set to YES, additional index files will be -# generated that can be used as input for Apple's Xcode 3 integrated development -# environment (see: http://developer.apple.com/tools/xcode/), introduced with -# OSX 10.5 (Leopard). To create a documentation set, doxygen will generate a -# Makefile in the HTML output directory. Running make will produce the docset in -# that directory and running make install will install the docset in -# ~/Library/Developer/Shared/Documentation/DocSets so that Xcode will find it at -# startup. See http://developer.apple.com/tools/creatingdocsetswithdoxygen.html -# for more information. -# The default value is: NO. -# This tag requires that the tag GENERATE_HTML is set to YES. - -GENERATE_DOCSET = NO - -# This tag determines the name of the docset feed. A documentation feed provides -# an umbrella under which multiple documentation sets from a single provider -# (such as a company or product suite) can be grouped. -# The default value is: Doxygen generated docs. -# This tag requires that the tag GENERATE_DOCSET is set to YES. - -DOCSET_FEEDNAME = "Doxygen generated docs" - -# This tag specifies a string that should uniquely identify the documentation -# set bundle. This should be a reverse domain-name style string, e.g. -# com.mycompany.MyDocSet. Doxygen will append .docset to the name. -# The default value is: org.doxygen.Project. -# This tag requires that the tag GENERATE_DOCSET is set to YES. - -DOCSET_BUNDLE_ID = org.doxygen.Project - -# The DOCSET_PUBLISHER_ID tag specifies a string that should uniquely identify -# the documentation publisher. This should be a reverse domain-name style -# string, e.g. com.mycompany.MyDocSet.documentation. -# The default value is: org.doxygen.Publisher. -# This tag requires that the tag GENERATE_DOCSET is set to YES. - -DOCSET_PUBLISHER_ID = org.doxygen.Publisher - -# The DOCSET_PUBLISHER_NAME tag identifies the documentation publisher. -# The default value is: Publisher. -# This tag requires that the tag GENERATE_DOCSET is set to YES. - -DOCSET_PUBLISHER_NAME = Publisher - -# If the GENERATE_HTMLHELP tag is set to YES then doxygen generates three -# additional HTML index files: index.hhp, index.hhc, and index.hhk. The -# index.hhp is a project file that can be read by Microsoft's HTML Help Workshop -# (see: http://www.microsoft.com/en-us/download/details.aspx?id=21138) on -# Windows. -# -# The HTML Help Workshop contains a compiler that can convert all HTML output -# generated by doxygen into a single compiled HTML file (.chm). Compiled HTML -# files are now used as the Windows 98 help format, and will replace the old -# Windows help format (.hlp) on all Windows platforms in the future. Compressed -# HTML files also contain an index, a table of contents, and you can search for -# words in the documentation. The HTML workshop also contains a viewer for -# compressed HTML files. -# The default value is: NO. -# This tag requires that the tag GENERATE_HTML is set to YES. - -GENERATE_HTMLHELP = NO - -# The CHM_FILE tag can be used to specify the file name of the resulting .chm -# file. You can add a path in front of the file if the result should not be -# written to the html output directory. -# This tag requires that the tag GENERATE_HTMLHELP is set to YES. - -CHM_FILE = - -# The HHC_LOCATION tag can be used to specify the location (absolute path -# including file name) of the HTML help compiler (hhc.exe). If non-empty, -# doxygen will try to run the HTML help compiler on the generated index.hhp. -# The file has to be specified with full path. -# This tag requires that the tag GENERATE_HTMLHELP is set to YES. - -HHC_LOCATION = - -# The GENERATE_CHI flag controls if a separate .chi index file is generated -# (YES) or that it should be included in the master .chm file (NO). -# The default value is: NO. -# This tag requires that the tag GENERATE_HTMLHELP is set to YES. - -GENERATE_CHI = NO - -# The CHM_INDEX_ENCODING is used to encode HtmlHelp index (hhk), content (hhc) -# and project file content. -# This tag requires that the tag GENERATE_HTMLHELP is set to YES. - -CHM_INDEX_ENCODING = - -# The BINARY_TOC flag controls whether a binary table of contents is generated -# (YES) or a normal table of contents (NO) in the .chm file. Furthermore it -# enables the Previous and Next buttons. -# The default value is: NO. -# This tag requires that the tag GENERATE_HTMLHELP is set to YES. - -BINARY_TOC = NO - -# The TOC_EXPAND flag can be set to YES to add extra items for group members to -# the table of contents of the HTML help documentation and to the tree view. -# The default value is: NO. -# This tag requires that the tag GENERATE_HTMLHELP is set to YES. - -TOC_EXPAND = NO - -# If the GENERATE_QHP tag is set to YES and both QHP_NAMESPACE and -# QHP_VIRTUAL_FOLDER are set, an additional index file will be generated that -# can be used as input for Qt's qhelpgenerator to generate a Qt Compressed Help -# (.qch) of the generated HTML documentation. -# The default value is: NO. -# This tag requires that the tag GENERATE_HTML is set to YES. - -GENERATE_QHP = NO - -# If the QHG_LOCATION tag is specified, the QCH_FILE tag can be used to specify -# the file name of the resulting .qch file. The path specified is relative to -# the HTML output folder. -# This tag requires that the tag GENERATE_QHP is set to YES. - -QCH_FILE = - -# The QHP_NAMESPACE tag specifies the namespace to use when generating Qt Help -# Project output. For more information please see Qt Help Project / Namespace -# (see: http://qt-project.org/doc/qt-4.8/qthelpproject.html#namespace). -# The default value is: org.doxygen.Project. -# This tag requires that the tag GENERATE_QHP is set to YES. - -QHP_NAMESPACE = org.doxygen.Project - -# The QHP_VIRTUAL_FOLDER tag specifies the namespace to use when generating Qt -# Help Project output. For more information please see Qt Help Project / Virtual -# Folders (see: http://qt-project.org/doc/qt-4.8/qthelpproject.html#virtual- -# folders). -# The default value is: doc. -# This tag requires that the tag GENERATE_QHP is set to YES. - -QHP_VIRTUAL_FOLDER = doc - -# If the QHP_CUST_FILTER_NAME tag is set, it specifies the name of a custom -# filter to add. For more information please see Qt Help Project / Custom -# Filters (see: http://qt-project.org/doc/qt-4.8/qthelpproject.html#custom- -# filters). -# This tag requires that the tag GENERATE_QHP is set to YES. - -QHP_CUST_FILTER_NAME = - -# The QHP_CUST_FILTER_ATTRS tag specifies the list of the attributes of the -# custom filter to add. For more information please see Qt Help Project / Custom -# Filters (see: http://qt-project.org/doc/qt-4.8/qthelpproject.html#custom- -# filters). -# This tag requires that the tag GENERATE_QHP is set to YES. - -QHP_CUST_FILTER_ATTRS = - -# The QHP_SECT_FILTER_ATTRS tag specifies the list of the attributes this -# project's filter section matches. Qt Help Project / Filter Attributes (see: -# http://qt-project.org/doc/qt-4.8/qthelpproject.html#filter-attributes). -# This tag requires that the tag GENERATE_QHP is set to YES. - -QHP_SECT_FILTER_ATTRS = - -# The QHG_LOCATION tag can be used to specify the location of Qt's -# qhelpgenerator. If non-empty doxygen will try to run qhelpgenerator on the -# generated .qhp file. -# This tag requires that the tag GENERATE_QHP is set to YES. - -QHG_LOCATION = - -# If the GENERATE_ECLIPSEHELP tag is set to YES, additional index files will be -# generated, together with the HTML files, they form an Eclipse help plugin. To -# install this plugin and make it available under the help contents menu in -# Eclipse, the contents of the directory containing the HTML and XML files needs -# to be copied into the plugins directory of eclipse. The name of the directory -# within the plugins directory should be the same as the ECLIPSE_DOC_ID value. -# After copying Eclipse needs to be restarted before the help appears. -# The default value is: NO. -# This tag requires that the tag GENERATE_HTML is set to YES. - -GENERATE_ECLIPSEHELP = NO - -# A unique identifier for the Eclipse help plugin. When installing the plugin -# the directory name containing the HTML and XML files should also have this -# name. Each documentation set should have its own identifier. -# The default value is: org.doxygen.Project. -# This tag requires that the tag GENERATE_ECLIPSEHELP is set to YES. - -ECLIPSE_DOC_ID = org.doxygen.Project - -# If you want full control over the layout of the generated HTML pages it might -# be necessary to disable the index and replace it with your own. The -# DISABLE_INDEX tag can be used to turn on/off the condensed index (tabs) at top -# of each HTML page. A value of NO enables the index and the value YES disables -# it. Since the tabs in the index contain the same information as the navigation -# tree, you can set this option to YES if you also set GENERATE_TREEVIEW to YES. -# The default value is: NO. -# This tag requires that the tag GENERATE_HTML is set to YES. - -DISABLE_INDEX = NO - -# The GENERATE_TREEVIEW tag is used to specify whether a tree-like index -# structure should be generated to display hierarchical information. If the tag -# value is set to YES, a side panel will be generated containing a tree-like -# index structure (just like the one that is generated for HTML Help). For this -# to work a browser that supports JavaScript, DHTML, CSS and frames is required -# (i.e. any modern browser). Windows users are probably better off using the -# HTML help feature. Via custom style sheets (see HTML_EXTRA_STYLESHEET) one can -# further fine-tune the look of the index. As an example, the default style -# sheet generated by doxygen has an example that shows how to put an image at -# the root of the tree instead of the PROJECT_NAME. Since the tree basically has -# the same information as the tab index, you could consider setting -# DISABLE_INDEX to YES when enabling this option. -# The default value is: NO. -# This tag requires that the tag GENERATE_HTML is set to YES. - -GENERATE_TREEVIEW = NO - -# The ENUM_VALUES_PER_LINE tag can be used to set the number of enum values that -# doxygen will group on one line in the generated HTML documentation. -# -# Note that a value of 0 will completely suppress the enum values from appearing -# in the overview section. -# Minimum value: 0, maximum value: 20, default value: 4. -# This tag requires that the tag GENERATE_HTML is set to YES. - -ENUM_VALUES_PER_LINE = 4 - -# If the treeview is enabled (see GENERATE_TREEVIEW) then this tag can be used -# to set the initial width (in pixels) of the frame in which the tree is shown. -# Minimum value: 0, maximum value: 1500, default value: 250. -# This tag requires that the tag GENERATE_HTML is set to YES. - -TREEVIEW_WIDTH = 250 - -# If the EXT_LINKS_IN_WINDOW option is set to YES, doxygen will open links to -# external symbols imported via tag files in a separate window. -# The default value is: NO. -# This tag requires that the tag GENERATE_HTML is set to YES. - -EXT_LINKS_IN_WINDOW = NO - -# Use this tag to change the font size of LaTeX formulas included as images in -# the HTML documentation. When you change the font size after a successful -# doxygen run you need to manually remove any form_*.png images from the HTML -# output directory to force them to be regenerated. -# Minimum value: 8, maximum value: 50, default value: 10. -# This tag requires that the tag GENERATE_HTML is set to YES. - -FORMULA_FONTSIZE = 10 - -# Use the FORMULA_TRANPARENT tag to determine whether or not the images -# generated for formulas are transparent PNGs. Transparent PNGs are not -# supported properly for IE 6.0, but are supported on all modern browsers. -# -# Note that when changing this option you need to delete any form_*.png files in -# the HTML output directory before the changes have effect. -# The default value is: YES. -# This tag requires that the tag GENERATE_HTML is set to YES. - -FORMULA_TRANSPARENT = YES - -# Enable the USE_MATHJAX option to render LaTeX formulas using MathJax (see -# http://www.mathjax.org) which uses client side Javascript for the rendering -# instead of using pre-rendered bitmaps. Use this if you do not have LaTeX -# installed or if you want to formulas look prettier in the HTML output. When -# enabled you may also need to install MathJax separately and configure the path -# to it using the MATHJAX_RELPATH option. -# The default value is: NO. -# This tag requires that the tag GENERATE_HTML is set to YES. - -USE_MATHJAX = NO - -# When MathJax is enabled you can set the default output format to be used for -# the MathJax output. See the MathJax site (see: -# http://docs.mathjax.org/en/latest/output.html) for more details. -# Possible values are: HTML-CSS (which is slower, but has the best -# compatibility), NativeMML (i.e. MathML) and SVG. -# The default value is: HTML-CSS. -# This tag requires that the tag USE_MATHJAX is set to YES. - -MATHJAX_FORMAT = HTML-CSS - -# When MathJax is enabled you need to specify the location relative to the HTML -# output directory using the MATHJAX_RELPATH option. The destination directory -# should contain the MathJax.js script. For instance, if the mathjax directory -# is located at the same level as the HTML output directory, then -# MATHJAX_RELPATH should be ../mathjax. The default value points to the MathJax -# Content Delivery Network so you can quickly see the result without installing -# MathJax. However, it is strongly recommended to install a local copy of -# MathJax from http://www.mathjax.org before deployment. -# The default value is: http://cdn.mathjax.org/mathjax/latest. -# This tag requires that the tag USE_MATHJAX is set to YES. - -MATHJAX_RELPATH = http://cdn.mathjax.org/mathjax/latest - -# The MATHJAX_EXTENSIONS tag can be used to specify one or more MathJax -# extension names that should be enabled during MathJax rendering. For example -# MATHJAX_EXTENSIONS = TeX/AMSmath TeX/AMSsymbols -# This tag requires that the tag USE_MATHJAX is set to YES. - -MATHJAX_EXTENSIONS = - -# The MATHJAX_CODEFILE tag can be used to specify a file with javascript pieces -# of code that will be used on startup of the MathJax code. See the MathJax site -# (see: http://docs.mathjax.org/en/latest/output.html) for more details. For an -# example see the documentation. -# This tag requires that the tag USE_MATHJAX is set to YES. - -MATHJAX_CODEFILE = - -# When the SEARCHENGINE tag is enabled doxygen will generate a search box for -# the HTML output. The underlying search engine uses javascript and DHTML and -# should work on any modern browser. Note that when using HTML help -# (GENERATE_HTMLHELP), Qt help (GENERATE_QHP), or docsets (GENERATE_DOCSET) -# there is already a search function so this one should typically be disabled. -# For large projects the javascript based search engine can be slow, then -# enabling SERVER_BASED_SEARCH may provide a better solution. It is possible to -# search using the keyboard; to jump to the search box use + S -# (what the is depends on the OS and browser, but it is typically -# , /