From fa357c118dd2db162ce8288046c871f835b9b37f Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 19 Mar 2024 15:43:50 +0100 Subject: [PATCH 001/129] Add isFci to check if a field is a FCI field. --- include/bout/field.hxx | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/include/bout/field.hxx b/include/bout/field.hxx index c0693ec0fb..c0ce04dbed 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -178,6 +178,11 @@ inline bool areFieldsCompatible(const Field& field1, const Field& field2) { #define ASSERT1_FIELDS_COMPATIBLE(field1, field2) ; #endif +template +inline bool isFci(const F& f) { + return not f.getCoordinates()->getParallelTransform().canToFromFieldAligned(); +} + /// Return an empty shell field of some type derived from Field, with metadata /// copied and a data array that is allocated but not initialised. template From ca59edb366be164e3583dc6041ea26d28c569ed8 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 3 Nov 2023 14:04:13 +0100 Subject: [PATCH 002/129] Improve isFci Ensure this does not crash if coordinates or transform is not set. In this case no FCI transformation is set, and this returns false. --- include/bout/field.hxx | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/include/bout/field.hxx b/include/bout/field.hxx index c0ce04dbed..4433af6d19 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -180,7 +180,14 @@ inline bool areFieldsCompatible(const Field& field1, const Field& field2) { template inline bool isFci(const F& f) { - return not f.getCoordinates()->getParallelTransform().canToFromFieldAligned(); + const auto coords = f.getCoordinates(); + if (coords == nullptr){ + return false; + } + if (not coords->hasParallelTransform()) { + return false; + } + return not coords->getParallelTransform().canToFromFieldAligned(); } /// Return an empty shell field of some type derived from Field, with metadata From 8b1fbba650fbaa17c572c7036c4a3d949892cece Mon Sep 17 00:00:00 2001 From: dschwoerer Date: Tue, 19 Mar 2024 15:27:28 +0000 Subject: [PATCH 003/129] Apply clang-format changes --- include/bout/field.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/bout/field.hxx b/include/bout/field.hxx index 4433af6d19..61e4af4d4b 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -181,7 +181,7 @@ inline bool areFieldsCompatible(const Field& field1, const Field& field2) { template inline bool isFci(const F& f) { const auto coords = f.getCoordinates(); - if (coords == nullptr){ + if (coords == nullptr) { return false; } if (not coords->hasParallelTransform()) { From 23d309f2532376be77e653535e9e67f96915d20b Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 20 Mar 2024 10:28:57 +0100 Subject: [PATCH 004/129] Make isFci a member function and add missing func --- include/bout/coordinates.hxx | 5 +++-- include/bout/field.hxx | 14 ++------------ src/field/field.cxx | 11 +++++++++++ 3 files changed, 16 insertions(+), 14 deletions(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index 49feffa0a7..5ea2d276fc 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -133,9 +133,10 @@ public: transform = std::move(pt); } + bool hasParallelTransform() const { return transform != nullptr; } /// Return the parallel transform - ParallelTransform& getParallelTransform() { - ASSERT1(transform != nullptr); + ParallelTransform& getParallelTransform() const { + ASSERT1(hasParallelTransform()); return *transform; } diff --git a/include/bout/field.hxx b/include/bout/field.hxx index 61e4af4d4b..d0828e8f6c 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -86,6 +86,8 @@ public: std::string name; + bool isFci() const; + #if CHECK > 0 // Routines to test guard/boundary cells set @@ -178,18 +180,6 @@ inline bool areFieldsCompatible(const Field& field1, const Field& field2) { #define ASSERT1_FIELDS_COMPATIBLE(field1, field2) ; #endif -template -inline bool isFci(const F& f) { - const auto coords = f.getCoordinates(); - if (coords == nullptr) { - return false; - } - if (not coords->hasParallelTransform()) { - return false; - } - return not coords->getParallelTransform().canToFromFieldAligned(); -} - /// Return an empty shell field of some type derived from Field, with metadata /// copied and a data array that is allocated but not initialised. template diff --git a/src/field/field.cxx b/src/field/field.cxx index e48a8f3ef7..c9373454bf 100644 --- a/src/field/field.cxx +++ b/src/field/field.cxx @@ -39,3 +39,14 @@ int Field::getNx() const { return getMesh()->LocalNx; } int Field::getNy() const { return getMesh()->LocalNy; } int Field::getNz() const { return getMesh()->LocalNz; } + +bool Field::isFci() const { + const auto coords = this->getCoordinates(); + if (coords == nullptr) { + return false; + } + if (not coords->hasParallelTransform()) { + return false; + } + return not coords->getParallelTransform().canToFromFieldAligned(); +} From 18fca62f3367981fc173ba58297715e9e3f2c351 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 12 Apr 2024 15:00:21 +0200 Subject: [PATCH 005/129] Make it explicit that this is a pointer --- src/field/field.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/field/field.cxx b/src/field/field.cxx index c9373454bf..797df6c405 100644 --- a/src/field/field.cxx +++ b/src/field/field.cxx @@ -41,7 +41,7 @@ int Field::getNy() const { return getMesh()->LocalNy; } int Field::getNz() const { return getMesh()->LocalNz; } bool Field::isFci() const { - const auto coords = this->getCoordinates(); + const auto* coords = this->getCoordinates(); if (coords == nullptr) { return false; } From 44084cca1fa150f760102aa10d97b4be714ed10e Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 2 Jul 2024 14:38:47 +0200 Subject: [PATCH 006/129] Add isFci also to mesh --- include/bout/mesh.hxx | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index c80716fc12..a3a36ad933 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -828,6 +828,17 @@ public: ASSERT1(RegionID.has_value()); return region3D[RegionID.value()]; } + bool isFci() const { + const auto coords = this->getCoordinatesConst(); + if (coords == nullptr) { + return false; + } + if (not coords->hasParallelTransform()) { + return false; + } + return not coords->getParallelTransform().canToFromFieldAligned(); + } + private: /// Allocates default Coordinates objects From cd288e98d1964d517380ce883fcb83da1106e3f0 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 9 Aug 2024 10:47:16 +0200 Subject: [PATCH 007/129] Add const version for getCoordinates --- include/bout/mesh.hxx | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index a3a36ad933..ccc979987b 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -636,6 +636,19 @@ public: return inserted.first->second; } + std::shared_ptr + getCoordinatesConst(const CELL_LOC location = CELL_CENTRE) const { + ASSERT1(location != CELL_DEFAULT); + ASSERT1(location != CELL_VSHIFT); + + auto found = coords_map.find(location); + if (found != coords_map.end()) { + // True branch most common, returns immediately + return found->second; + } + throw BoutException("Coordinates not yet set. Use non-const version!"); + } + /// Returns the non-CELL_CENTRE location /// allowed as a staggered location CELL_LOC getAllowedStaggerLoc(DIRECTION direction) const { From 061b82ffbd14373f744ad205f6577c3c03120b3a Mon Sep 17 00:00:00 2001 From: dschwoerer Date: Tue, 22 Oct 2024 10:04:02 +0000 Subject: [PATCH 008/129] Apply clang-format changes --- include/bout/mesh.hxx | 1 - 1 file changed, 1 deletion(-) diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index ccc979987b..c1a6a1336d 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -852,7 +852,6 @@ public: return not coords->getParallelTransform().canToFromFieldAligned(); } - private: /// Allocates default Coordinates objects /// By default attempts to read staggered Coordinates from grid data source, From ececfb831dff6ae17aff36924fe01e97c9eeb7d7 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 6 Nov 2024 16:46:26 +0100 Subject: [PATCH 009/129] Apply suggestions from clang-tidy --- include/bout/mesh.hxx | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index c1a6a1336d..ed9f6f8d60 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -651,7 +651,7 @@ public: /// Returns the non-CELL_CENTRE location /// allowed as a staggered location - CELL_LOC getAllowedStaggerLoc(DIRECTION direction) const { + static CELL_LOC getAllowedStaggerLoc(DIRECTION direction) { AUTO_TRACE(); switch (direction) { case (DIRECTION::X): @@ -860,8 +860,7 @@ private: /// (useful if CELL_CENTRE Coordinates have been changed, so reading from file /// would not be correct). std::shared_ptr - createDefaultCoordinates(const CELL_LOC location, - bool force_interpolate_from_centre = false); + createDefaultCoordinates(CELL_LOC location, bool force_interpolate_from_centre = false); //Internal region related information std::map regionMap3D; From aebca7ea5714c1e79297d81022325e59049799fb Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 10 Feb 2025 15:59:34 +0000 Subject: [PATCH 010/129] CI: Bump clang-tidy-review --- .github/workflows/clang-tidy-review.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/clang-tidy-review.yml b/.github/workflows/clang-tidy-review.yml index 8ddcd4ae3f..834caa4d22 100644 --- a/.github/workflows/clang-tidy-review.yml +++ b/.github/workflows/clang-tidy-review.yml @@ -22,7 +22,7 @@ jobs: submodules: true - name: Run clang-tidy - uses: ZedThree/clang-tidy-review@v0.19.0 + uses: ZedThree/clang-tidy-review@v0.20.1 id: review with: annotations: true @@ -48,4 +48,4 @@ jobs: -DBOUT_UPDATE_GIT_SUBMODULE=OFF - name: Upload clang-tidy fixes - uses: ZedThree/clang-tidy-review/upload@v0.19.0 + uses: ZedThree/clang-tidy-review/upload@v0.20.1 From cc8c683a97a098aeb8875d4519ec3b5881e3b495 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 10 Feb 2025 17:28:10 +0000 Subject: [PATCH 011/129] CI: Disable annotations for clang-tidy-review --- .github/workflows/clang-tidy-review.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/clang-tidy-review.yml b/.github/workflows/clang-tidy-review.yml index 834caa4d22..1e2f88d208 100644 --- a/.github/workflows/clang-tidy-review.yml +++ b/.github/workflows/clang-tidy-review.yml @@ -25,7 +25,6 @@ jobs: uses: ZedThree/clang-tidy-review@v0.20.1 id: review with: - annotations: true build_dir: build apt_packages: "libfftw3-dev,libnetcdf-c++4-dev,libopenmpi-dev,petsc-dev,slepc-dev,liblapack-dev,libparpack2-dev,libsundials-dev,uuid-dev" config_file: ".clang-tidy" From 9063bf16186d29fbc73a2cb3fd0f1b8e8af7ce7c Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 11 Feb 2025 09:47:23 +0100 Subject: [PATCH 012/129] CI: Show more output of dnf5 to help debugging --- .ci_fedora.sh | 1 + 1 file changed, 1 insertion(+) diff --git a/.ci_fedora.sh b/.ci_fedora.sh index 77e7f45055..359516ab33 100755 --- a/.ci_fedora.sh +++ b/.ci_fedora.sh @@ -41,6 +41,7 @@ then # Ignore weak depencies echo "install_weak_deps=False" >> /etc/dnf/dnf.conf echo "minrate=10M" >> /etc/dnf/dnf.conf + export FORCE_COLUMNS=200 time dnf -y install dnf5 time dnf5 -y install dnf5-plugins cmake python3-zoidberg python3-natsort # Allow to override packages - see #2073 From d377379e7d5af7373411fb2650d8759f586cd012 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 4 Mar 2025 17:04:21 -0800 Subject: [PATCH 013/129] elm-pb example: Adding relaxing phi boundaries Adds an option phi_boundary_relax. If set to true, then the radial boundary conditions on the potential phi are relaxed over a given timescale towards zero gradient. Adapted from Hermes-3 (https://github.com/boutproject/hermes-3/blob/master/src/vorticity.cxx#L261). --- examples/elm-pb/elm_pb.cxx | 310 ++++++++++++++++++++++++++---- examples/elm-pb/relaxing/BOUT.inp | 299 ++++++++++++++++++++++++++++ 2 files changed, 573 insertions(+), 36 deletions(-) create mode 100644 examples/elm-pb/relaxing/BOUT.inp diff --git a/examples/elm-pb/elm_pb.cxx b/examples/elm-pb/elm_pb.cxx index f830f3d98a..4c34dee0ad 100644 --- a/examples/elm-pb/elm_pb.cxx +++ b/examples/elm-pb/elm_pb.cxx @@ -17,6 +17,9 @@ #include #include +#include +#include + #include #if BOUT_HAS_HYPRE @@ -46,6 +49,7 @@ class ELMpb : public PhysicsModel { Coordinates::FieldMetric U0; // 0th vorticity of equilibrium flow, // radial flux coordinate, normalized radial flux coordinate + bool laplace_perp; // Use Laplace_perp or Delp2? bool constn0; // the total height, average width and center of profile of N0 BoutReal n0_height, n0_ave, n0_width, n0_center, n0_bottom_x, Nbar, Tibar, Tebar; @@ -203,7 +207,11 @@ class ELMpb : public PhysicsModel { bool parallel_lr_diff; // Use left and right shifted stencils for parallel differences - bool phi_constraint; // Solver for phi using a solver constraint + bool phi_constraint; // Solver for phi using a solver constraint + bool phi_boundary_relax; // Relax x boundaries of phi towards Neumann? + bool phi_core_averagey; // Average phi core boundary in Y? + BoutReal phi_boundary_timescale; // Relaxation timescale + BoutReal phi_boundary_last_update; // Time when last updated bool include_rmp; // Include RMP coil perturbation bool simple_rmp; // Just use a simple form for the perturbation @@ -350,7 +358,8 @@ class ELMpb : public PhysicsModel { ////////////////////////////////////////////////////////////// auto& globalOptions = Options::root(); auto& options = globalOptions["highbeta"]; - + laplace_perp = options["laplace_perp"].withDefault(false); + // Use Laplace_perp rather than Delp2 constn0 = options["constn0"].withDefault(true); // use the hyperbolic profile of n0. If both n0_fake_prof and // T0_fake_prof are false, use the profiles from grid file @@ -377,6 +386,9 @@ class ELMpb : public PhysicsModel { phi_constraint = options["phi_constraint"] .doc("Use solver constraint for phi?") .withDefault(false); + phi_boundary_relax = options["phi_boundary_relax"] + .doc("Relax x boundaries of phi towards Neumann?") + .withDefault(false); // Effects to include/exclude include_curvature = options["include_curvature"].withDefault(true); @@ -1148,6 +1160,24 @@ class ELMpb : public PhysicsModel { aparSolver = Laplacian::create(&globalOptions["aparSolver"], loc); + if (phi_boundary_relax) { + // Set the last update time to -1, so it will reset + // the first time RHS function is called + phi_boundary_last_update = -1.; + phi_core_averagey = options["phi_core_averagey"] + .doc("Average phi core boundary in Y?") + .withDefault(false) + and mesh->periodicY(mesh->xstart); + + phi_boundary_timescale = options["phi_boundary_timescale"] + .doc("Timescale for phi boundary relaxation [seconds]") + .withDefault(1e-7) + / Tbar; // Normalise time units to Tbar + + phiSolver->setInnerBoundaryFlags(INVERT_SET); + phiSolver->setOuterBoundaryFlags(INVERT_SET); + } + /////////////// CHECK VACUUM /////////////////////// // In vacuum region, initial vorticity should equal zero @@ -1315,10 +1345,124 @@ class ELMpb : public PhysicsModel { Ctmp.applyBoundary(); Ctmp -= phi; // Now contains error in the boundary - C_phi = Delp2(phi) - U; // Error in the bulk + if (laplace_perp) { + C_phi = Laplace_perp(phi) - U; // Error in the bulk + } else { + C_phi = Delp2(phi) - U; // Error in the bulk + } C_phi.setBoundaryTo(Ctmp); } else { + if (phi_boundary_relax) { + // Update the boundary regions by relaxing towards zero gradient + // on a given timescale. + + if (phi_boundary_last_update < 0.0) { + // First time this has been called. + phi_boundary_last_update = t; + } else if (t > phi_boundary_last_update) { + // Only update if simulation time has advanced + // Uses an exponential decay of the weighting of the value in the boundary + // so that the solution is well behaved for arbitrary steps + BoutReal weight = exp(-(t - phi_boundary_last_update) / phi_boundary_timescale); + phi_boundary_last_update = t; + + if (mesh->firstX()) { + BoutReal phivalue = 0.0; + if (phi_core_averagey) { + // Calculate a single phi boundary value for all Y slices + BoutReal philocal = 0.0; + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + philocal += phi(mesh->xstart, j, k); + } + } + MPI_Comm comm_inner = mesh->getYcomm(0); + int np; + MPI_Comm_size(comm_inner, &np); + MPI_Allreduce(&philocal, &phivalue, 1, MPI_DOUBLE, MPI_SUM, comm_inner); + phivalue /= (np * mesh->LocalNz * mesh->LocalNy); + } + for (int j = mesh->ystart; j <= mesh->yend; j++) { + if (!phi_core_averagey) { + phivalue = 0.0; // Calculate phi boundary for each Y index separately + for (int k = 0; k < mesh->LocalNz; k++) { + phivalue += phi(mesh->xstart, j, k); + } + phivalue /= mesh->LocalNz; // Average in Z of point next to boundary + } + + // Old value of phi at boundary. Note: this is constant in Z + BoutReal oldvalue = + 0.5 * (phi(mesh->xstart - 1, j, 0) + phi(mesh->xstart, j, 0)); + + // New value of phi at boundary, relaxing towards phivalue + BoutReal newvalue = weight * oldvalue + (1. - weight) * phivalue; + + // Set phi at the boundary to this value + for (int k = 0; k < mesh->LocalNz; k++) { + phi(mesh->xstart - 1, j, k) = 2. * newvalue - phi(mesh->xstart, j, k); + phi(mesh->xstart - 2, j, k) = phi(mesh->xstart - 1, j, k); + } + } + } + + if (mesh->lastX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + BoutReal phivalue = 0.0; + BoutReal oldvalue = 0.0; + for (int k = 0; k < mesh->LocalNz; k++) { + phivalue = phi(mesh->xend, j, k); + oldvalue = 0.5 * (phi(mesh->xend + 1, j, k) + phi(mesh->xend, j, k)); + } + + // New value of phi at boundary, relaxing towards phivalue + BoutReal newvalue = weight * oldvalue + (1. - weight) * phivalue; + + // Set phi at the boundary to this value + for (int k = 0; k < mesh->LocalNz; k++) { + phi(mesh->xend + 1, j, k) = 2. * newvalue - phi(mesh->xend, j, k); + phi(mesh->xend + 2, j, k) = phi(mesh->xend + 1, j, k); + } + } + } + } + } + + Field3D phi_shift = phi; + if (constn0 and diamag) { + // Solving for phi + ion pressure term + phi_shift += 0.5 * dnorm * P / B0; + } else { + // Ensure that memory is not shared between phi and phi_shift + phi_shift.allocate(); + } + + // Update boundary conditions. + // The INVERT_SET flag takes the value in the guard (boundary) cell + // and sets the boundary between cells to this value. + // This shift by 1/2 grid cell is important. + + if (mesh->firstX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + // Average phi + Pi at the boundary, and set the boundary cell + // to this value. The phi solver will then put the value back + // onto the cell mid-point + phi_shift(mesh->xstart - 1, j, k) = + 0.5 * (phi_shift(mesh->xstart - 1, j, k) + phi_shift(mesh->xstart, j, k)); + } + } + } + + if (mesh->lastX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + phi_shift(mesh->xend + 1, j, k) = + 0.5 * (phi_shift(mesh->xend + 1, j, k) + phi_shift(mesh->xend, j, k)); + } + } + } if (constn0) { if (split_n0) { @@ -1326,44 +1470,84 @@ class ELMpb : public PhysicsModel { // Boussinesq, split // Split into axisymmetric and non-axisymmetric components Field2D Vort2D = DC(U); // n=0 component + Field2D phi_shift_2d = phi2D; + + if (phi_boundary_relax) { + phi_shift_2d = DC(phi_shift); + } + phi_shift -= phi_shift_2d; // Applies boundary condition for "phi". phi2D.applyBoundary(t); // Solve axisymmetric (n=0) part - phi2D = laplacexy->solve(Vort2D, phi2D); + phi2D = laplacexy->solve(Vort2D, phi_shift_2d); // Solve non-axisymmetric part - phi = phiSolver->solve(U - Vort2D); + phi = phiSolver->solve(U - Vort2D, phi_shift); phi += phi2D; // Add axisymmetric part } else { - phi = phiSolver->solve(U); + if (phi_boundary_relax) { + phi = phiSolver->solve(U, phi_shift); + } else { + phi = phiSolver->solve(U); + } } - if (diamag) { phi -= 0.5 * dnorm * P / B0; } } else { ubyn = U / N0; if (diamag) { - ubyn -= 0.5 * dnorm / (N0 * B0) * Delp2(P); + if (laplace_perp) { + ubyn -= 0.5 * dnorm / (N0 * B0) * Laplace_perp(P); + } else { + ubyn -= 0.5 * dnorm / (N0 * B0) * Delp2(P); + } mesh->communicate(ubyn); } // Invert laplacian for phi phiSolver->setCoefC(N0); phi = phiSolver->solve(ubyn); } - // Apply a boundary condition on phi for target plates - phi.applyBoundary(); mesh->communicate(phi); } + if (mesh->firstX()) { + for (int i = mesh->xstart - 2; i >= 0; --i) { + for (int j = mesh->ystart; j <= mesh->yend; ++j) { + for (int k = 0; k < mesh->LocalNz; ++k) { + phi(i, j, k) = phi(i + 1, j, k); + } + } + } + } + + if (mesh->lastX()) { + for (int i = mesh->xend + 2; i < mesh->LocalNx; ++i) { + for (int j = mesh->ystart; j <= mesh->yend; ++j) { + for (int k = 0; k < mesh->LocalNz; ++k) { + phi(i, j, k) = phi(i - 1, j, k); + } + } + } + } + if (!evolve_jpar) { // Get J from Psi - Jpar = Delp2(Psi); + if (laplace_perp) { + Jpar = Laplace_perp(Psi); + } else { + Jpar = Delp2(Psi); + } + if (include_rmp) { - Jpar += Delp2(rmp_Psi); + if (laplace_perp) { + Jpar += Laplace_perp(rmp_Psi); + } else { + Jpar += Delp2(rmp_Psi); + } } Jpar.applyBoundary(); @@ -1397,8 +1581,11 @@ class ELMpb : public PhysicsModel { } // Get Delp2(J) from J - Jpar2 = Delp2(Jpar); - + if (laplace_perp) { + Jpar2 = Laplace_perp(Jpar); + } else { + Jpar2 = Delp2(Jpar); + } Jpar2.applyBoundary(); mesh->communicate(Jpar2); @@ -1494,7 +1681,11 @@ class ELMpb : public PhysicsModel { // Jpar Field3D B0U = B0 * U; mesh->communicate(B0U); - ddt(Jpar) = -Grad_parP(B0U, loc) / B0 + eta * Delp2(Jpar); + if (laplace_perp) { + ddt(Jpar) = -Grad_parP(B0U, loc) / B0 + eta * Laplace_perp(Jpar); + } else { + ddt(Jpar) = -Grad_parP(B0U, loc) / B0 + eta * Delp2(Jpar); + } if (relax_j_vac) { // Make ddt(Jpar) relax to zero. @@ -1524,11 +1715,19 @@ class ELMpb : public PhysicsModel { } if (hyperresist > 0.0) { // Hyper-resistivity - ddt(Psi) -= eta * hyperresist * Delp2(Jpar); + if (laplace_perp) { + ddt(Psi) -= eta * hyperresist * Laplace_perp(Jpar); + } else { + ddt(Psi) -= eta * hyperresist * Delp2(Jpar); + } } if (ehyperviscos > 0.0) { // electron Hyper-viscosity coefficient - ddt(Psi) -= eta * ehyperviscos * Delp2(Jpar2); + if (laplace_perp) { + ddt(Psi) -= eta * ehyperviscos * Laplace_perp(Jpar2); + } else { + ddt(Psi) -= eta * ehyperviscos * Delp2(Jpar2); + } } // Parallel hyper-viscous diffusion for vector potential @@ -1599,7 +1798,11 @@ class ELMpb : public PhysicsModel { } if (viscos_perp > 0.0) { - ddt(U) += viscos_perp * Delp2(U); // Perpendicular viscosity + if (laplace_perp) { + ddt(U) += viscos_perp * Laplace_perp(U); // Perpendicular viscosity + } else { + ddt(U) += viscos_perp * Delp2(U); // Perpendicular viscosity + } } // Hyper-viscosity @@ -1626,21 +1829,39 @@ class ELMpb : public PhysicsModel { Pi = 0.5 * P; Pi0 = 0.5 * P0; - Dperp2Phi0 = Field3D(Delp2(B0 * phi0)); - Dperp2Phi0.applyBoundary(); - mesh->communicate(Dperp2Phi0); + if (laplace_perp) { + Dperp2Phi0 = Field3D(Laplace_perp(B0 * phi0)); + Dperp2Phi0.applyBoundary(); + mesh->communicate(Dperp2Phi0); + + Dperp2Phi = Laplace_perp(B0 * phi); + Dperp2Phi.applyBoundary(); + mesh->communicate(Dperp2Phi); + + Dperp2Pi0 = Field3D(Laplace_perp(Pi0)); + Dperp2Pi0.applyBoundary(); + mesh->communicate(Dperp2Pi0); + + Dperp2Pi = Laplace_perp(Pi); + Dperp2Pi.applyBoundary(); + mesh->communicate(Dperp2Pi); + } else { + Dperp2Phi0 = Field3D(Delp2(B0 * phi0)); + Dperp2Phi0.applyBoundary(); + mesh->communicate(Dperp2Phi0); - Dperp2Phi = Delp2(B0 * phi); - Dperp2Phi.applyBoundary(); - mesh->communicate(Dperp2Phi); + Dperp2Phi = Delp2(B0 * phi); + Dperp2Phi.applyBoundary(); + mesh->communicate(Dperp2Phi); - Dperp2Pi0 = Field3D(Delp2(Pi0)); - Dperp2Pi0.applyBoundary(); - mesh->communicate(Dperp2Pi0); + Dperp2Pi0 = Field3D(Delp2(Pi0)); + Dperp2Pi0.applyBoundary(); + mesh->communicate(Dperp2Pi0); - Dperp2Pi = Delp2(Pi); - Dperp2Pi.applyBoundary(); - mesh->communicate(Dperp2Pi); + Dperp2Pi = Delp2(Pi); + Dperp2Pi.applyBoundary(); + mesh->communicate(Dperp2Pi); + } bracketPhi0P = bracket(B0 * phi0, Pi, bm_exb); bracketPhi0P.applyBoundary(); @@ -1658,8 +1879,13 @@ class ELMpb : public PhysicsModel { mesh->communicate(B0phi0); ddt(U) += 0.5 * Upara2 * bracket(B0phi, Dperp2Pi0, bm_exb) / B0; ddt(U) += 0.5 * Upara2 * bracket(B0phi0, Dperp2Pi, bm_exb) / B0; - ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhi0P) / B0; - ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP0) / B0; + if (laplace_perp) { + ddt(U) -= 0.5 * Upara2 * Laplace_perp(bracketPhi0P) / B0; + ddt(U) -= 0.5 * Upara2 * Laplace_perp(bracketPhiP0) / B0; + } else { + ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhi0P) / B0; + ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP0) / B0; + } if (nonlinear) { Field3D B0phi = B0 * phi; @@ -1670,7 +1896,11 @@ class ELMpb : public PhysicsModel { ddt(U) -= 0.5 * Upara2 * bracket(Pi, Dperp2Phi, bm_exb) / B0; ddt(U) += 0.5 * Upara2 * bracket(B0phi, Dperp2Pi, bm_exb) / B0; - ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP) / B0; + if (laplace_perp) { + ddt(U) -= 0.5 * Upara2 * Laplace_perp(bracketPhiP) / B0; + } else { + ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP) / B0; + } } } @@ -1808,7 +2038,12 @@ class ELMpb : public PhysicsModel { int precon(BoutReal UNUSED(t), BoutReal gamma, BoutReal UNUSED(delta)) { // First matrix, applying L mesh->communicate(ddt(Psi)); - Field3D Jrhs = Delp2(ddt(Psi)); + Field3D Jrhs; + if (laplace_perp) { + Jrhs = Laplace_perp(ddt(Psi)); + } else { + Jrhs = Delp2(ddt(Psi)); + } Jrhs.applyBoundary("neumann"); if (jpar_bndry_width > 0) { @@ -1880,8 +2115,11 @@ class ELMpb : public PhysicsModel { phi = phiSolver->solve(ddt(U)); - Jpar = Delp2(ddt(Psi)); - + if (laplace_perp) { + Jpar = Laplace_perp(ddt(Psi)); + } else { + Jpar = Delp2(ddt(Psi)); + } mesh->communicate(phi, Jpar); Field3D JP = -b0xGrad_dot_Grad(phi, P0); diff --git a/examples/elm-pb/relaxing/BOUT.inp b/examples/elm-pb/relaxing/BOUT.inp new file mode 100644 index 0000000000..d9489611e6 --- /dev/null +++ b/examples/elm-pb/relaxing/BOUT.inp @@ -0,0 +1,299 @@ +# settings file for BOUT++ +# High-Beta reduced MHD case + +################################################## +# Global settings used by the core code + +nout = 1000 # number of time-steps +timestep = 1 # time between outputs + +zperiod = 5 # Fraction of a torus to simulate +MZ = 64 # Number of points in Z + +grid = "cbm18_dens8.grid_nx68ny64.nc" #"cbm18_dens3_0.5BS_516nx64ny.grid.nc" + +[mesh] + +staggergrids = false # Use staggered grids + +[mesh:paralleltransform] +#type = shifted # Use shifted metric method +type = shiftedinterp + +################################################## +# derivative methods + +[mesh:ddx] +first = C4 # order of first x derivatives +second = C4 # order of second x derivatives +upwind = W3 # order of upwinding method W3 = Weno3 + +[mesh:ddy] +first = C4 +second = C4 +upwind = W3 +flux = C2 + +[mesh:ddz] +first = C4 # Z derivatives can be done using FFT +second = C4 +upwind = W3 + +################################################## +# FFTs + +[fft] + +fft_measurement_flag = measure # If using FFTW, perform tests to determine fastest method + +[output] + +#type = adios +#shiftoutput = true # Put the output into field-aligned coordinates + +[restart_files] + +#type = adios + +[laplace] + +#type = hypre3d +#rtol = 1.e-9 +#atol = 1.e-14 +#pctype = sor # Preconditioner + +#atol = 1e-12 +#rtol = 1e-08 + + +################################################## +# Solver settings + +[solver] + +# mudq, mldq, mukeep, mlkeep preconditioner options +atol = 1.0e-8 # absolute tolerance +rtol = 1.0e-5 # relative tolerance + +use_precon = false # Use preconditioner: User-supplied or BBD + +mxstep = 50000 # Number of internal steps between outputs + +################################################## +# settings for high-beta reduced MHD + +[highbeta] + +density = 1.0e19 # number density of deuterium [m^-3] + # used to produce output normalisations +constn0 = true +n0_fake_prof = false + +sheath_boundaries = false + +evolve_jpar = false # If true, evolve J raher than Psi + +evolve_pressure = true # If false, switch off all pressure evolution + +phi_constraint = false # Solve phi as a constraint (DAE system, needs IDA) + +## Effects to include/exclude + +include_jpar0 = true # determines whether to include jpar0 terms +include_curvature = true # include curvature drive term? + +compress = false # set compressible (evolve Vpar) +nonlinear = true # include non-linear terms? + +diamag = true # Include diamagnetic effects? +diamag_grad_t = false # Include Grad_par(Te) term in Psi equation +diamag_phi0 = true # Balance ExB against Vd for stationary equilibrium + +split_n0 = false + +laplace_perp = true + +################################################## +# BRACKET_METHOD flags: +# 0:BRACKET_STD; derivative methods will be determined +# by the choices C or W in this input file +# 1:BRACKET_SIMPLE; 2:BRACKET_ARAKAWA; 3:BRACKET_CTU. + +bm_exb_flag = 0 +bm_mag_flag = 0 +################################################################## +withflow = false # With flow or not +D_0 = 130000 # differential potential +D_s = 20 # shear parameter +K_H_term = false # Contain K-H term +sign = -1 # flow direction +x0 = 0.855 # peak location +D_min = 3000 # constant +################################################################## + +eHall = false # Include electron pressure effects in Ohm's law? +AA = 2.0 # ion mass in units of proton mass +Zeff = 1.0 +#Zi = 1.0 + +noshear = false # zero all shear + +relax_j_vac = false # Relax to zero-current in the vacuum +relax_j_tconst = 1e-2 # Time constant for vacuum relaxation + +## Toroidal filtering +filter_z = false # remove all except one mode +filter_z_mode = 1 # Specify which harmonic to keep (1 = fundamental) +low_pass_z = 16 # Keep up to and including this harmonic (-1 = keep all) +zonal_flow = true # keep zonal component of vorticity? +zonal_field = false # keep zonal component of Psi? +zonal_bkgd = true # keep zonal component of P? + +## Jpar smoothing +smooth_j_x = true # Filter Jpar in the X direction + +## Magnetic perturbations +include_rmp = false # Read RMP data from grid file + +simple_rmp = false # Enable/disable a simple model of RMP +rmp_n = 3 # Toroidal mode number +rmp_m = 6 # Poloidal mode number +rmp_factor = 1.e-4 # Amplitude of Apar [Tm] +rmp_ramp = 1.e-4 # Timescale [s] of ramp +rmp_polwid = -1.0 # Width of Gaussian factor (< 0 = No Gaussian) +rmp_polpeak = 0.5 # Y location of maximum (fraction) + +## Vacuum region control + +vacuum_pressure = 0.02 # the pressure below which it is considered vacuum + # fraction of peak pressure +vacuum_trans = 0.01 # transition width (fraction of P) + +## Resistivity and Hyper-resistivity + +vac_lund = 1.0e8 # Lundquist number in vacuum (negative -> infinity) +core_lund = 1.0e8 # Lundquist number in core (negative -> infinity) +hyperresist = 1.e-4 # Hyper-resistivity coefficient (like 1 / Lundquist number) + +## Inner boundary damping + +damp_width = -1 # Width of damping region (grid cells) +damp_t_const = 1e-2 # Damping time constant + +## Parallel pressure diffusion + +diffusion_par = 1.0e-2 # Parallel pressure diffusion (< 0 = none) +diffusion_p4 = -1e-05 # parallel hyper-viscous diffusion for pressure (< 0 = none) +diffusion_u4 = -1e-05 # parallel hyper-viscous diffusion for vorticity (< 0 = none) +diffusion_a4 = -1e-05 # parallel hyper-viscous diffusion for vector potential (< 0 = none) + +## heat source in pressure in watts + +heating_P = -1 # heat power in watts (< 0 = none) +hp_width = 0.1 # heat width, in percentage of nx (< 0 = none) +hp_length = 0.3 # heat length in percentage of nx (< 0 = none) + +## sink rate in pressure + +sink_P = -1 # sink rate in pressure (< 0 = none) +sp_width = 0.04 # sink width, in percentage of nx (< 0 = none) +sp_length = 0.15 # sink length in percentage of nx (< 0 = none) + + +## left edge sink rate in vorticity +sink_Ul = 10.0 # left edge sink rate in vorticity (< 0 = none) +su_widthl = 0.06 # left edge sink width, in percentage of nx (< 0 = none) +su_lengthl = 0.1 # left edge sink length in percentage of nx (< 0 = none) + +## right edge sink rate in vorticity +sink_Ur = 10.0 # right edge sink rate in vorticity (< 0 = none) +su_widthr = 0.06 # right edge sink width, in percentage of nx (< 0 = none) +su_lengthr = 0.1 # right edge sink length in percentage of nx (< 0 = none) + +## Viscosity and Hyper-viscosity + +viscos_par = 0.1 # Parallel viscosity (< 0 = none) +viscos_perp = 1.0e-7 # Perpendicular viscosity (< 0 = none) +hyperviscos = -1.0 # Radial hyper viscosity + +## Compressional terms (only when compress = true) +phi_curv = true # Include curvature*Grad(phi) in P equation +# gamma = 1.6666 + +phi_boundary_relax = true # Relax phi at radial boundaries towards zero gradient +phi_core_averagey = false +phi_boundary_timescale = 1.0e-6 # In seconds + +[phiSolver] +# INVERT_DC_GRAD = 1 // Zero-gradient for DC (constant in Z) component. Default is zero value +# INVERT_AC_GRAD = 2 // Zero-gradient for AC (non-constant in Z) component. Default is zero value +# INVERT_AC_LAP = 4 // Use zero-laplacian (decaying solution) to AC component +# INVERT_DC_LAP = 64 // Use zero-laplacian solution for DC component +#type = hypre3d +#rtol = 1.e-9 +#atol = 1.e-14 +#inner_boundary_flags = 0 # 0 for Dirichlet; 2 for Neumann +#outer_boundary_flags = 0 # 0 for Dirichlet; 2 for Neumann +#hypre_print_level = 1 # print information on the matrices, solver settings, and iterations. + +[aparSolver] +#type = hypre3d +#rtol = 1.e-9 +#atol = 1.e-14 +#inner_boundary_flags = 4 # INVERT_AC_LAP +#outer_boundary_flags = 1 + 4 # INVERT_DC_GRAD + INVERT_AC_LAP + +################################################## +# settings for individual variables +# The section "All" defines default settings for all variables +# These can be overridden for individual variables in +# a section of that name. + +[all] +scale = 0.0 # default size of initial perturbations + +# boundary conditions +# ------------------- +# dirichlet - Zero value +# neumann - Zero gradient +# zerolaplace - Laplacian = 0, decaying solution +# constlaplace - Laplacian = const, decaying solution +# +# relax( ) - Make boundary condition relaxing + +bndry_all = dirichlet_o2 # Default to zero-value + +[U] # vorticity +scale = 1e-05 +function = ballooning(gauss(x-0.5,0.1)*gauss(y-pi,0.6*pi)*sin(z),3) + +[P] # pressure +bndry_core = dirichlet +#scale = 1.0e-5 + +[Psi] # Vector potential + +# zero laplacian +bndry_xin = zerolaplace +bndry_xout = zerolaplace + +[Psi_loc] # for staggering + +bndry_xin = zerolaplace +bndry_xout = zerolaplace + +bndry_yup = free_o3 +bndry_ydown = free_o3 + +[J] # parallel current + +# Zero gradient in the core +bndry_core = neumann + +[phi] + +#bndry_core = neumann +bndry_xin = none +bndry_xout = none +#bndry_target = neumann + From 5ce6cae6ea94a76417a4c058fdd848db384d36d1 Mon Sep 17 00:00:00 2001 From: Rob Falgout Date: Tue, 4 Mar 2025 17:04:44 -0800 Subject: [PATCH 014/129] First version of code to eliminate boundary equations to improve HYPRE solves --- include/bout/hypre_interface.hxx | 193 +++++++++++++- .../laplace/impls/hypre3d/hypre3d_laplace.cxx | 6 + .../laplace/impls/hypre3d/hypre_boundary.c | 252 ++++++++++++++++++ 3 files changed, 450 insertions(+), 1 deletion(-) create mode 100644 src/invert/laplace/impls/hypre3d/hypre_boundary.c diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index ae392de4f3..72eec916e2 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -20,6 +20,56 @@ #include "HYPRE_utilities.h" #include "_hypre_utilities.h" +HYPRE_Int +AdjustBCMatrixEquations( + HYPRE_Int nrows, + HYPRE_Int *ncols, + HYPRE_BigInt *rows, + HYPRE_Int **row_indexes_ptr, + HYPRE_BigInt *cols, + HYPRE_Complex *values, + HYPRE_Int nb, // number of boundary equations + HYPRE_Int *bi_array, // row i for each boundary equation + HYPRE_Int **binum_array_ptr, // data index for row i (for each boundary equation) + HYPRE_Int **bjnum_array_ptr, // data index for col j (for each boundary equation) + HYPRE_Complex **bii_array_ptr, // coefficient b_ii (for each boundary equation) + HYPRE_Complex **bij_array_ptr, // coefficient b_ij (for each boundary equation) + HYPRE_Int *na_ptr, // number of interior equations to adjust + HYPRE_Int **aknum_array_ptr, // data index for row k (for each interior equation) + HYPRE_Complex **aki_array_ptr); // coefficient a_ki (for each interior equation) + +HYPRE_Int +AdjustBCRightHandSideEquations( + HYPRE_Complex *rhs, + HYPRE_Int nb, + HYPRE_Int *binum_array, + HYPRE_Complex *bii_array, + HYPRE_Complex *bij_array, + HYPRE_Complex **brhs_array_ptr, + HYPRE_Int na, + HYPRE_Int *aknum_array, + HYPRE_Complex *aki_array); + +HYPRE_Int +AdjustBCSolutionEquations( + HYPRE_Complex *solution, + HYPRE_Int nb, + HYPRE_Int *binum_array, + HYPRE_Int *bjnum_array, + HYPRE_Complex *bii_array, + HYPRE_Complex *bij_array, + HYPRE_Complex *brhs_array); + +HYPRE_Int +AdjustBCEquationsFree( + HYPRE_Int *binum_array, + HYPRE_Int *bjnum_array, + HYPRE_Complex *bii_array, + HYPRE_Complex *bij_array, + HYPRE_Complex *brhs_array, + HYPRE_Int *aknum_array, + HYPRE_Complex *aki_array); + #include // BOUT_ENUM_CLASS does not work inside namespaces @@ -174,6 +224,24 @@ public: HypreMalloc(V, vsize * sizeof(HYPRE_Complex)); } + // Data for eliminating boundary equations (TODO: convert this to a structure) + bool elimBErhs = false; + bool elimBEsol = false; + HYPRE_Complex *rhs; + HYPRE_Int nb; + HYPRE_Int *binum_array; + HYPRE_Int *bjnum_array; + HYPRE_Complex *bii_array; + HYPRE_Complex *bij_array; + HYPRE_Complex *brhs_array; // This is not created from the HypreMatrix class + HYPRE_Int na; + HYPRE_Int *aknum_array; + HYPRE_Complex *aki_array; + + void syncElimBErhs(HypreVector& rhs) { + brhs_array = rhs.brhs_array; + } + void assemble() { CALI_CXX_MARK_FUNCTION; writeCacheToHypre(); @@ -183,11 +251,20 @@ public: } void writeCacheToHypre() { + if (elimBErhs) + { + AdjustBCRightHandSideEquations(V, nb, binum_array, bii_array, bij_array, &brhs_array, + na, aknum_array, aki_array); + } checkHypreError(HYPRE_IJVectorSetValues(hypre_vector, vsize, I, V)); } void readCacheFromHypre() { checkHypreError(HYPRE_IJVectorGetValues(hypre_vector, vsize, I, V)); + if (elimBEsol) + { + AdjustBCSolutionEquations(V, nb, binum_array, bjnum_array, bii_array, bij_array, brhs_array); + } } T toField() { @@ -667,6 +744,43 @@ public: return Element(*this, global_row, global_column, positions, weights); } + // Data for eliminating boundary equations (TODO: convert this to a structure) + bool elimBE = false; + HYPRE_Int nb; + HYPRE_Int *binum_array; + HYPRE_Int *bjnum_array; + HYPRE_Complex *bii_array; + HYPRE_Complex *bij_array; + HYPRE_Int na; + HYPRE_Int *aknum_array; + HYPRE_Complex *aki_array; + + void setElimBE() { + elimBE = true; + } + + void setElimBEVectors(HypreVector& sol, HypreVector& rhs) { + sol.elimBEsol = elimBE; + sol.nb = nb; + sol.binum_array = binum_array; + sol.bjnum_array = bjnum_array; + sol.bii_array = bii_array; + sol.bij_array = bij_array; + sol.na = na; + sol.aknum_array = aknum_array; + sol.aki_array = aki_array; + + rhs.elimBErhs = elimBE; + rhs.nb = nb; + rhs.binum_array = binum_array; + rhs.bjnum_array = bjnum_array; + rhs.bii_array = bii_array; + rhs.bij_array = bij_array; + rhs.na = na; + rhs.aknum_array = aknum_array; + rhs.aki_array = aki_array; + } + void assemble() { CALI_CXX_MARK_FUNCTION; @@ -695,8 +809,85 @@ public: entry++; } } - checkHypreError( + + // Eliminate boundary condition equations in hypre SetValues input arguments + if (elimBE) + { + HYPRE_Int *bi_array; + HYPRE_Int *row_indexes; + // There must be an easier way to get nb + nb = 0; + BOUT_FOR_SERIAL(i, index_converter->getRegionBndry()) { + nb++; + } + HypreMalloc(bi_array, nb * sizeof(HYPRE_Int)); + nb = 0; + BOUT_FOR_SERIAL(i, index_converter->getRegionBndry()) { + bi_array[nb] = index_converter->getGlobal(i); + nb++; + } + AdjustBCMatrixEquations(num_rows, num_cols, rawI, &row_indexes, cols, vals, + nb, bi_array, + &binum_array, &bjnum_array, &bii_array, &bij_array, + &na, &aknum_array, &aki_array); + HypreFree(bi_array); +// { +// FILE *file; +// int i; +// file = fopen("zbout.setvalues.out2", "w"); +// fprintf(file, "nrows %d\n", num_rows); +// for (i = 0; i < num_rows; i++) +// { +// fprintf(file, "ncols %d: %d\n", i, num_cols[i]); +// } +// for (i = 0; i < num_rows; i++) +// { +// fprintf(file, "rows %d: %d\n", i, rawI[i]); +// } +// for (i = 0; i < num_entries; i++) +// { +// fprintf(file, "cols %d: %d\n", i, cols[i]); +// } +// for (i = 0; i < num_entries; i++) +// { +// fprintf(file, "vals %d: %f\n", i, vals[i]); +// } +// fclose(file); +// +// file = fopen("zbout.setvalues.out3", "w"); +// for (i = 0; i < num_rows; i++) +// { +// fprintf(file, "row_indexes %d: %d\n", i, row_indexes[i]); +// } +// fprintf(file, "nb %d\n", nb); +// for (i = 0; i < nb; i++) +// { +// fprintf(file, "bijnum_array %d: %d %d\n", i, binum_array[i], bjnum_array[i]); +// } +// for (i = 0; i < nb; i++) +// { +// fprintf(file, "biiij_array %d: %f %f\n", i, bii_array[i], bij_array[i]); +// } +// fprintf(file, "na %d\n", na); +// for (i = 0; i < na; i++) +// { +// fprintf(file, "aknum_array %d: %d\n", i, aknum_array[i]); +// } +// for (i = 0; i < na; i++) +// { +// fprintf(file, "bki_array %d: %f\n", i, aki_array[i]); +// } +// fclose(file); +// } + checkHypreError( + HYPRE_IJMatrixSetValues2(*hypre_matrix, num_rows, num_cols, rawI, row_indexes, cols, vals)); + HypreFree(row_indexes); + } + else + { + checkHypreError( HYPRE_IJMatrixSetValues(*hypre_matrix, num_rows, num_cols, rawI, cols, vals)); + } checkHypreError(HYPRE_IJMatrixAssemble(*hypre_matrix)); checkHypreError(HYPRE_IJMatrixGetObject(*hypre_matrix, reinterpret_cast(¶llel_matrix))); diff --git a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx index c50be1db85..def54b36b5 100644 --- a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx +++ b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx @@ -28,6 +28,7 @@ #if BOUT_HAS_HYPRE #include "hypre3d_laplace.hxx" +#include "hypre_boundary.c" #include #include @@ -218,11 +219,15 @@ Field3D LaplaceHypre3d::solve(const Field3D& b_in, const Field3D& x0) { CALI_MARK_BEGIN("LaplaceHypre3d_solve:vectorAssemble"); + operator3D.setElimBEVectors(solution, rhs); + rhs.importValuesFromField(b); solution.importValuesFromField(x0); rhs.assemble(); solution.assemble(); + solution.syncElimBErhs(rhs); + CALI_MARK_END("LaplaceHypre3d_solve:vectorAssemble"); CALI_MARK_BEGIN("LaplaceHypre3d_solve:solve"); @@ -411,6 +416,7 @@ void LaplaceHypre3d::updateMatrix3D() { operator3D.ydown(ydown)(l, l.ym().zp()) += -C_d2f_dydz; operator3D.ydown(ydown)(l, l.ym().zm()) += C_d2f_dydz; } + operator3D.setElimBE(); operator3D.assemble(); if (print_matrix) { diff --git a/src/invert/laplace/impls/hypre3d/hypre_boundary.c b/src/invert/laplace/impls/hypre3d/hypre_boundary.c new file mode 100644 index 0000000000..c4e7393c3f --- /dev/null +++ b/src/invert/laplace/impls/hypre3d/hypre_boundary.c @@ -0,0 +1,252 @@ + +/* + * This function modifies the input for the HYPRE_IJMatrixSetValues() routine to + * eliminate the boundary condition equations (see below for details on how the + * equations are adjusted). It modifies the arrays ncols, rows, cols, and + * values. It also returns a row_indexes array. This can then be passed to the + * HYPRE_IJMatrixSetValues2() routine to set up the matrix in hypre. + * + * The arguments nb and bi_array indicate the boundary equations. The routine + * returns info needed to adjust the right-hand-side and solution vector through + * the functions AdjustRightHandSideEquations and AdjustSolutionEquations. + * + * Returned info arrays can be freed with the function AdjustEquationsFree(). + * + * NOTE: It may make sense from an organizational standpoint to collect many of + * these arguments in a structure of some sort. + * + * Notation, assumptions, and other details: + * + * - Boundary equation i is assumed to have two coefficients + * + * b_ii * u_i + b_ij * u_j = rhs_i + * + * - We also assume that each boundary equation has only one interior equation k + * coupled to it (such that k = j) with coupling coefficient a_ki + * + * a_ki * u_i + a_kj * u_j + ... = rhs_k + * + * - Each equation k is adjusted as follows: + * + * a_kj = a_kj - a_ki * b_ij / b_ii + * a_ki = 0 + * + * - Boundary equations are adjusted to be identity equations in the matrix, but + * the boundary coefficients (b_ii, b_ij) are returned for use later + * + * - Right-hand-side equations are adjusted in AdjustRightHandSideEquations() as + * follows: rhs_k = rhs_k - a_ki * rhs_i / b_ii + * + * - Solution unknowns are adjusted at boundaries in AdjustSolutionEquations as + * follows: u_i = (rhs_i - b_ij * u_j) / b_ii + * + * - Naming conventions: Arrays starting with 'b' are boundary equation arrays + * indexed by 'bnum', and arrays starting with 'a' are non-boundary arrays + * (interior matrix equations) indexed by 'anum'. When 'num' is prefixed with + * a row or column number 'i', 'j', or 'k', the array holds the corresponding + * local data index for that row or column (e.g., an index into the local + * solution vector). Matrix coefficients are named as above, e.g., 'bij' is + * the coefficient for b_ij. + */ + +HYPRE_Int +AdjustBCMatrixEquations( + HYPRE_Int nrows, + HYPRE_Int *ncols, + HYPRE_BigInt *rows, + HYPRE_Int **row_indexes_ptr, + HYPRE_BigInt *cols, + HYPRE_Complex *values, + HYPRE_Int nb, // number of boundary equations + HYPRE_Int *bi_array, // row i for each boundary equation + HYPRE_Int **binum_array_ptr, // data index for row i (for each boundary equation) + HYPRE_Int **bjnum_array_ptr, // data index for col j (for each boundary equation) + HYPRE_Complex **bii_array_ptr, // coefficient b_ii (for each boundary equation) + HYPRE_Complex **bij_array_ptr, // coefficient b_ij (for each boundary equation) + HYPRE_Int *na_ptr, // number of interior equations to adjust + HYPRE_Int **aknum_array_ptr, // data index for row k (for each interior equation) + HYPRE_Complex **aki_array_ptr) // coefficient a_ki (for each interior equation) +{ + HYPRE_Int *row_indexes; + HYPRE_Int na, *binum_array, *bjnum_array, *aknum_array; + HYPRE_Complex *bii_array, *bij_array, *aki_array; + HYPRE_Int i, j, k, m, mkj, anum, bnum, acoeffnum, bcoeffnum; + HYPRE_Int binum, aknum; + HYPRE_Complex bii, bij, aki; + + /* Create the row_indexes array */ + row_indexes = (HYPRE_Int *)malloc(sizeof(HYPRE_Int) * nrows); + row_indexes[0] = 0; + for (i = 1; i < nrows; i++) + { + row_indexes[i] = row_indexes[i-1] + ncols[i-1]; + } + + /* Assume just one interior equation coupled to each boundary equation */ + na = nb; + + /* Allocate return arrays */ + HypreMalloc(binum_array, sizeof(HYPRE_Int) * nb); + HypreMalloc(bjnum_array, sizeof(HYPRE_Int) * nb); + HypreMalloc(bii_array, sizeof(HYPRE_Complex) * nb); + HypreMalloc(bij_array, sizeof(HYPRE_Complex) * nb); + HypreMalloc(aknum_array, sizeof(HYPRE_Int) * na); + HypreMalloc(aki_array, sizeof(HYPRE_Complex) * na); + + binum = 0; + aknum = 0; + for (bnum = 0; bnum < nb; bnum++) + { + /* Get boundary equation information and adjust boundary equations */ + /* Find row i in rows array (assume i increases and rows is sorted) */ + i = bi_array[bnum]; + for (; binum < nrows; binum++) + { + if (i == rows[binum]) + { + break; // Found row i in rows array + } + } + bcoeffnum = row_indexes[binum]; + for (m = 0; m < 2; m++) // Assume only two boundary equation coefficients + { + if (cols[bcoeffnum + m] == i) + { + bii = values[bcoeffnum + m]; + values[bcoeffnum + m] = -1.0; // Identity equation (negative definite matrix) + } + else + { + j = cols[bcoeffnum + m]; + bij = values[bcoeffnum + m]; + values[bcoeffnum + m] = 0.0; // Identity equation + } + } + ncols[binum] = 1; // Identity equation + + /* Get interior equation information and adjust interior equations */ + /* Find row k in rows array (assume k increases and rows is sorted) */ + k = j; // Assume equation k = j + for (; aknum < nrows; aknum++) + { + if (k == rows[aknum]) + { + break; // Found row k in rows array + } + } + acoeffnum = row_indexes[aknum]; + for (m = 0; m < ncols[aknum]; m++) + { + if (cols[acoeffnum + m] == j) + { + mkj = m; // Save for update of akj value below + } + if (cols[acoeffnum + m] == i) + { + aki = values[acoeffnum + m]; + values[acoeffnum + m] = 0.0; // Eliminate coupling to boundary equation + } + } + values[acoeffnum + mkj] -= aki * bij / bii; // Update akj value + + /* Update return arrays */ + anum = bnum; // Assume only one interior equation k + binum_array[bnum] = binum; + bjnum_array[bnum] = aknum; // Assume only one interior equation k + bii_array[bnum] = bii; + bij_array[bnum] = bij; + aknum_array[anum] = aknum; + aki_array[anum] = aki; + } + + /* Set return arguments */ + *row_indexes_ptr = row_indexes; + *binum_array_ptr = binum_array; + *bjnum_array_ptr = bjnum_array; + *bii_array_ptr = bii_array; + *bij_array_ptr = bij_array; + *na_ptr = na; + *aknum_array_ptr = aknum_array; + *aki_array_ptr = aki_array; + + return 0; +} + +HYPRE_Int +AdjustBCRightHandSideEquations( + HYPRE_Complex *rhs, + HYPRE_Int nb, + HYPRE_Int *binum_array, + HYPRE_Complex *bii_array, + HYPRE_Complex *bij_array, + HYPRE_Complex **brhs_array_ptr, + HYPRE_Int na, + HYPRE_Int *aknum_array, + HYPRE_Complex *aki_array) +{ + HYPRE_Complex *brhs_array; + HYPRE_Int anum, bnum, binum, aknum; + + HypreMalloc(brhs_array, sizeof(HYPRE_Complex) * nb); + + for (bnum = 0; bnum < nb; bnum++) + { + binum = binum_array[bnum]; + brhs_array[bnum] = rhs[binum]; + } + + for (anum = 0; anum < na; anum++) + { + bnum = anum; // Assume only one interior equation per boundary equation + aknum = aknum_array[anum]; + rhs[aknum] -= aki_array[anum] * bij_array[bnum] / bii_array[bnum]; + } + + *brhs_array_ptr = brhs_array; + + return 0; +} + +HYPRE_Int +AdjustBCSolutionEquations( + HYPRE_Complex *solution, + HYPRE_Int nb, + HYPRE_Int *binum_array, + HYPRE_Int *bjnum_array, + HYPRE_Complex *bii_array, + HYPRE_Complex *bij_array, + HYPRE_Complex *brhs_array) +{ + HYPRE_Int bnum, binum, bjnum; + + for (bnum = 0; bnum < nb; bnum++) + { + binum = binum_array[bnum]; + bjnum = bjnum_array[bnum]; + solution[binum] = (brhs_array[bnum] - bij_array[bnum] * solution[bjnum]) / bii_array[bnum]; + } + + return 0; +} + +HYPRE_Int +AdjustBCEquationsFree( + HYPRE_Int *binum_array, + HYPRE_Int *bjnum_array, + HYPRE_Complex *bii_array, + HYPRE_Complex *bij_array, + HYPRE_Complex *brhs_array, + HYPRE_Int *aknum_array, + HYPRE_Complex *aki_array) +{ + HypreFree(binum_array); + HypreFree(bjnum_array); + HypreFree(bii_array); + HypreFree(bij_array); + HypreFree(brhs_array); + HypreFree(aknum_array); + HypreFree(aki_array); + + return 0; +} + From 14a621b49efc130e9d5851a3138fc1352a719b18 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 4 Mar 2025 17:45:08 -0800 Subject: [PATCH 015/129] elm_pb example: Fix outer boundary phivalue Should average over Z --- examples/elm-pb/elm_pb.cxx | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/examples/elm-pb/elm_pb.cxx b/examples/elm-pb/elm_pb.cxx index 4c34dee0ad..b36bd728d4 100644 --- a/examples/elm-pb/elm_pb.cxx +++ b/examples/elm-pb/elm_pb.cxx @@ -1410,11 +1410,14 @@ class ELMpb : public PhysicsModel { if (mesh->lastX()) { for (int j = mesh->ystart; j <= mesh->yend; j++) { BoutReal phivalue = 0.0; - BoutReal oldvalue = 0.0; for (int k = 0; k < mesh->LocalNz; k++) { - phivalue = phi(mesh->xend, j, k); - oldvalue = 0.5 * (phi(mesh->xend + 1, j, k) + phi(mesh->xend, j, k)); + phivalue += phi(mesh->xend, j, k); } + phivalue /= mesh->LocalNz; // Average in Z of point next to boundary + + // Old value of phi at boundary. Note: this is constant in Z + BoutReal oldvalue = + 0.5 * (phi(mesh->xend + 1, j, 0) + phi(mesh->xend, j, 0)); // New value of phi at boundary, relaxing towards phivalue BoutReal newvalue = weight * oldvalue + (1. - weight) * phivalue; From e3f4a98ac6235b0bef6ffc4d90ec2da2f1ecb668 Mon Sep 17 00:00:00 2001 From: Rob Falgout Date: Wed, 5 Mar 2025 08:35:17 -0800 Subject: [PATCH 016/129] Added HYPRE GMRES SetKDim call. Need to make this an input file parameter. --- include/bout/hypre_interface.hxx | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index 72eec916e2..efb620241e 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -1067,6 +1067,22 @@ public: setMaxIter( options["maxits"].doc("Maximum iterations for Hypre solver").withDefault(10000)); + switch (solver_type) { + case HYPRE_SOLVER_TYPE::gmres: { + HYPRE_ParCSRGMRESSetKDim(solver, 30); // TODO: Make this an input file parameter + break; + } + case HYPRE_SOLVER_TYPE::bicgstab: { + break; + } + case HYPRE_SOLVER_TYPE::pcg: { + break; + } + default: { + throw BoutException("Unsupported hypre_solver_type {}", toString(solver_type)); + } + } + HYPRE_BoomerAMGCreate(&precon); HYPRE_BoomerAMGSetOldDefault(precon); #if BOUT_HAS_CUDA From ef069675e39a6e92355036c45f3f20183a67e628 Mon Sep 17 00:00:00 2001 From: Rob Falgout Date: Wed, 5 Mar 2025 09:53:21 -0800 Subject: [PATCH 017/129] Added code to free up data. The HypreMatrix destroy needs to be updated. --- examples/elm-pb/elm_pb.cxx | 134 ++++++++++++++---- include/bout/hypre_interface.hxx | 41 +++--- .../laplace/impls/hypre3d/hypre_boundary.c | 35 +---- 3 files changed, 129 insertions(+), 81 deletions(-) diff --git a/examples/elm-pb/elm_pb.cxx b/examples/elm-pb/elm_pb.cxx index f830f3d98a..c332edd223 100644 --- a/examples/elm-pb/elm_pb.cxx +++ b/examples/elm-pb/elm_pb.cxx @@ -46,6 +46,7 @@ class ELMpb : public PhysicsModel { Coordinates::FieldMetric U0; // 0th vorticity of equilibrium flow, // radial flux coordinate, normalized radial flux coordinate + bool laplace_perp; // Use Laplace_perp or Delp2? bool constn0; // the total height, average width and center of profile of N0 BoutReal n0_height, n0_ave, n0_width, n0_center, n0_bottom_x, Nbar, Tibar, Tebar; @@ -350,7 +351,8 @@ class ELMpb : public PhysicsModel { ////////////////////////////////////////////////////////////// auto& globalOptions = Options::root(); auto& options = globalOptions["highbeta"]; - + laplace_perp = options["laplace_perp"].withDefault(false); + // Use Laplace_perp rather than Delp2 constn0 = options["constn0"].withDefault(true); // use the hyperbolic profile of n0. If both n0_fake_prof and // T0_fake_prof are false, use the profiles from grid file @@ -1315,7 +1317,11 @@ class ELMpb : public PhysicsModel { Ctmp.applyBoundary(); Ctmp -= phi; // Now contains error in the boundary - C_phi = Delp2(phi) - U; // Error in the bulk + if (laplace_perp) { + C_phi = Laplace_perp(phi) - U; // Error in the bulk + } else { + C_phi = Delp2(phi) - U; // Error in the bulk + } C_phi.setBoundaryTo(Ctmp); } else { @@ -1347,8 +1353,12 @@ class ELMpb : public PhysicsModel { } else { ubyn = U / N0; if (diamag) { - ubyn -= 0.5 * dnorm / (N0 * B0) * Delp2(P); - mesh->communicate(ubyn); + if (laplace_perp) { + ubyn -= 0.5 * dnorm / (N0 * B0) * Laplace_perp(P); + } else { + ubyn -= 0.5 * dnorm / (N0 * B0) * Delp2(P); + } + mesh->communicate(ubyn); } // Invert laplacian for phi phiSolver->setCoefC(N0); @@ -1361,9 +1371,18 @@ class ELMpb : public PhysicsModel { if (!evolve_jpar) { // Get J from Psi - Jpar = Delp2(Psi); + if (laplace_perp) { + Jpar = Laplace_perp(Psi); + } else { + Jpar = Delp2(Psi); + } + if (include_rmp) { - Jpar += Delp2(rmp_Psi); + if (laplace_perp) { + Jpar += Laplace_perp(rmp_Psi); + } else { + Jpar += Delp2(rmp_Psi); + } } Jpar.applyBoundary(); @@ -1397,8 +1416,11 @@ class ELMpb : public PhysicsModel { } // Get Delp2(J) from J - Jpar2 = Delp2(Jpar); - + if (laplace_perp) { + Jpar2 = Laplace_perp(Jpar); + } else { + Jpar2 = Delp2(Jpar); + } Jpar2.applyBoundary(); mesh->communicate(Jpar2); @@ -1494,7 +1516,11 @@ class ELMpb : public PhysicsModel { // Jpar Field3D B0U = B0 * U; mesh->communicate(B0U); - ddt(Jpar) = -Grad_parP(B0U, loc) / B0 + eta * Delp2(Jpar); + if (laplace_perp) { + ddt(Jpar) = -Grad_parP(B0U, loc) / B0 + eta * Laplace_perp(Jpar); + }else { + ddt(Jpar) = -Grad_parP(B0U, loc) / B0 + eta * Delp2(Jpar); + } if (relax_j_vac) { // Make ddt(Jpar) relax to zero. @@ -1524,11 +1550,19 @@ class ELMpb : public PhysicsModel { } if (hyperresist > 0.0) { // Hyper-resistivity - ddt(Psi) -= eta * hyperresist * Delp2(Jpar); + if (laplace_perp) { + ddt(Psi) -= eta * hyperresist * Laplace_perp(Jpar); + } else { + ddt(Psi) -= eta * hyperresist * Delp2(Jpar); + } } if (ehyperviscos > 0.0) { // electron Hyper-viscosity coefficient - ddt(Psi) -= eta * ehyperviscos * Delp2(Jpar2); + if (laplace_perp) { + ddt(Psi) -= eta * ehyperviscos * Laplace_perp(Jpar2); + } else { + ddt(Psi) -= eta * ehyperviscos * Delp2(Jpar2); + } } // Parallel hyper-viscous diffusion for vector potential @@ -1599,7 +1633,11 @@ class ELMpb : public PhysicsModel { } if (viscos_perp > 0.0) { - ddt(U) += viscos_perp * Delp2(U); // Perpendicular viscosity + if (laplace_perp) { + ddt(U) += viscos_perp * Laplace_perp(U); // Perpendicular viscosity + } else { + ddt(U) += viscos_perp * Delp2(U); // Perpendicular viscosity + } } // Hyper-viscosity @@ -1626,21 +1664,40 @@ class ELMpb : public PhysicsModel { Pi = 0.5 * P; Pi0 = 0.5 * P0; - Dperp2Phi0 = Field3D(Delp2(B0 * phi0)); - Dperp2Phi0.applyBoundary(); - mesh->communicate(Dperp2Phi0); + if (laplace_perp) { + Dperp2Phi0 = Field3D(Laplace_perp(B0 * phi0)); + Dperp2Phi0.applyBoundary(); + mesh->communicate(Dperp2Phi0); - Dperp2Phi = Delp2(B0 * phi); - Dperp2Phi.applyBoundary(); - mesh->communicate(Dperp2Phi); + Dperp2Phi = Laplace_perp(B0 * phi); + Dperp2Phi.applyBoundary(); + mesh->communicate(Dperp2Phi); - Dperp2Pi0 = Field3D(Delp2(Pi0)); - Dperp2Pi0.applyBoundary(); - mesh->communicate(Dperp2Pi0); + Dperp2Pi0 = Field3D(Laplace_perp(Pi0)); + Dperp2Pi0.applyBoundary(); + mesh->communicate(Dperp2Pi0); - Dperp2Pi = Delp2(Pi); - Dperp2Pi.applyBoundary(); - mesh->communicate(Dperp2Pi); + Dperp2Pi = Laplace_perp(Pi); + Dperp2Pi.applyBoundary(); + mesh->communicate(Dperp2Pi); + } else { + Dperp2Phi0 = Field3D(Delp2(B0 * phi0)); + Dperp2Phi0.applyBoundary(); + mesh->communicate(Dperp2Phi0); + + Dperp2Phi = Delp2(B0 * phi); + Dperp2Phi.applyBoundary(); + mesh->communicate(Dperp2Phi); + + Dperp2Pi0 = Field3D(Delp2(Pi0)); + Dperp2Pi0.applyBoundary(); + mesh->communicate(Dperp2Pi0); + + Dperp2Pi = Delp2(Pi); + Dperp2Pi.applyBoundary(); + mesh->communicate(Dperp2Pi); + + } bracketPhi0P = bracket(B0 * phi0, Pi, bm_exb); bracketPhi0P.applyBoundary(); @@ -1658,8 +1715,13 @@ class ELMpb : public PhysicsModel { mesh->communicate(B0phi0); ddt(U) += 0.5 * Upara2 * bracket(B0phi, Dperp2Pi0, bm_exb) / B0; ddt(U) += 0.5 * Upara2 * bracket(B0phi0, Dperp2Pi, bm_exb) / B0; - ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhi0P) / B0; - ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP0) / B0; + if (laplace_perp) { + ddt(U) -= 0.5 * Upara2 * Laplace_perp(bracketPhi0P) / B0; + ddt(U) -= 0.5 * Upara2 * Laplace_perp(bracketPhiP0) / B0; + } else { + ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhi0P) / B0; + ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP0) / B0; + } if (nonlinear) { Field3D B0phi = B0 * phi; @@ -1670,7 +1732,11 @@ class ELMpb : public PhysicsModel { ddt(U) -= 0.5 * Upara2 * bracket(Pi, Dperp2Phi, bm_exb) / B0; ddt(U) += 0.5 * Upara2 * bracket(B0phi, Dperp2Pi, bm_exb) / B0; - ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP) / B0; + if (laplace_perp) { + ddt(U) -= 0.5 * Upara2 * Laplace_perp(bracketPhiP) / B0; + } else { + ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP) / B0; + } } } @@ -1808,7 +1874,12 @@ class ELMpb : public PhysicsModel { int precon(BoutReal UNUSED(t), BoutReal gamma, BoutReal UNUSED(delta)) { // First matrix, applying L mesh->communicate(ddt(Psi)); - Field3D Jrhs = Delp2(ddt(Psi)); + Field3D Jrhs; + if (laplace_perp) { + Jrhs = Laplace_perp(ddt(Psi)); + } else { + Jrhs = Delp2(ddt(Psi)); + } Jrhs.applyBoundary("neumann"); if (jpar_bndry_width > 0) { @@ -1880,8 +1951,11 @@ class ELMpb : public PhysicsModel { phi = phiSolver->solve(ddt(U)); - Jpar = Delp2(ddt(Psi)); - + if (laplace_perp) { + Jpar = Laplace_perp(ddt(Psi)); + } else { + Jpar = Delp2(ddt(Psi)); + } mesh->communicate(phi, Jpar); Field3D JP = -b0xGrad_dot_Grad(phi, P0); diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index efb620241e..0025ed062b 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -20,7 +20,7 @@ #include "HYPRE_utilities.h" #include "_hypre_utilities.h" -HYPRE_Int +void AdjustBCMatrixEquations( HYPRE_Int nrows, HYPRE_Int *ncols, @@ -38,7 +38,7 @@ AdjustBCMatrixEquations( HYPRE_Int **aknum_array_ptr, // data index for row k (for each interior equation) HYPRE_Complex **aki_array_ptr); // coefficient a_ki (for each interior equation) -HYPRE_Int +void AdjustBCRightHandSideEquations( HYPRE_Complex *rhs, HYPRE_Int nb, @@ -50,7 +50,7 @@ AdjustBCRightHandSideEquations( HYPRE_Int *aknum_array, HYPRE_Complex *aki_array); -HYPRE_Int +void AdjustBCSolutionEquations( HYPRE_Complex *solution, HYPRE_Int nb, @@ -60,16 +60,6 @@ AdjustBCSolutionEquations( HYPRE_Complex *bij_array, HYPRE_Complex *brhs_array); -HYPRE_Int -AdjustBCEquationsFree( - HYPRE_Int *binum_array, - HYPRE_Int *bjnum_array, - HYPRE_Complex *bii_array, - HYPRE_Complex *bij_array, - HYPRE_Complex *brhs_array, - HYPRE_Int *aknum_array, - HYPRE_Complex *aki_array); - #include // BOUT_ENUM_CLASS does not work inside namespaces @@ -136,6 +126,9 @@ public: checkHypreError(HYPRE_IJVectorDestroy(hypre_vector)); HypreFree(I); HypreFree(V); + if (elimBErhs) { + HypreFree(brhs_array); + } } // Disable copy, at least for now: not clear that HYPRE_IJVector is @@ -251,8 +244,7 @@ public: } void writeCacheToHypre() { - if (elimBErhs) - { + if (elimBErhs) { AdjustBCRightHandSideEquations(V, nb, binum_array, bii_array, bij_array, &brhs_array, na, aknum_array, aki_array); } @@ -261,8 +253,7 @@ public: void readCacheFromHypre() { checkHypreError(HYPRE_IJVectorGetValues(hypre_vector, vsize, I, V)); - if (elimBEsol) - { + if (elimBEsol) { AdjustBCSolutionEquations(V, nb, binum_array, bjnum_array, bii_array, bij_array, brhs_array); } } @@ -419,6 +410,19 @@ public: using ind_type = typename T::ind_type; HypreMatrix() = default; +// The 'if' block below needs to be called when the HypreMatrix is destroyed, +// but I don't know where to put it. The HypreMatrix class is handled +// differently from HypreVector. +// ~HypreVector() { +// if (elimBE) { +// HypreFree(binum_array); +// HypreFree(bjnum_array); +// HypreFree(bii_array); +// HypreFree(bij_array); +// HypreFree(aknum_array); +// HypreFree(aki_array); +// } +// } HypreMatrix(const HypreMatrix&) = delete; HypreMatrix(HypreMatrix&& other) : comm(other.comm), ilower(other.ilower), iupper(other.iupper), @@ -811,8 +815,7 @@ public: } // Eliminate boundary condition equations in hypre SetValues input arguments - if (elimBE) - { + if (elimBE) { HYPRE_Int *bi_array; HYPRE_Int *row_indexes; // There must be an easier way to get nb diff --git a/src/invert/laplace/impls/hypre3d/hypre_boundary.c b/src/invert/laplace/impls/hypre3d/hypre_boundary.c index c4e7393c3f..69e9a2d1b0 100644 --- a/src/invert/laplace/impls/hypre3d/hypre_boundary.c +++ b/src/invert/laplace/impls/hypre3d/hypre_boundary.c @@ -10,8 +10,6 @@ * returns info needed to adjust the right-hand-side and solution vector through * the functions AdjustRightHandSideEquations and AdjustSolutionEquations. * - * Returned info arrays can be freed with the function AdjustEquationsFree(). - * * NOTE: It may make sense from an organizational standpoint to collect many of * these arguments in a structure of some sort. * @@ -49,7 +47,7 @@ * the coefficient for b_ij. */ -HYPRE_Int +void AdjustBCMatrixEquations( HYPRE_Int nrows, HYPRE_Int *ncols, @@ -168,11 +166,9 @@ AdjustBCMatrixEquations( *na_ptr = na; *aknum_array_ptr = aknum_array; *aki_array_ptr = aki_array; - - return 0; } -HYPRE_Int +void AdjustBCRightHandSideEquations( HYPRE_Complex *rhs, HYPRE_Int nb, @@ -203,11 +199,9 @@ AdjustBCRightHandSideEquations( } *brhs_array_ptr = brhs_array; - - return 0; } -HYPRE_Int +void AdjustBCSolutionEquations( HYPRE_Complex *solution, HYPRE_Int nb, @@ -225,28 +219,5 @@ AdjustBCSolutionEquations( bjnum = bjnum_array[bnum]; solution[binum] = (brhs_array[bnum] - bij_array[bnum] * solution[bjnum]) / bii_array[bnum]; } - - return 0; -} - -HYPRE_Int -AdjustBCEquationsFree( - HYPRE_Int *binum_array, - HYPRE_Int *bjnum_array, - HYPRE_Complex *bii_array, - HYPRE_Complex *bij_array, - HYPRE_Complex *brhs_array, - HYPRE_Int *aknum_array, - HYPRE_Complex *aki_array) -{ - HypreFree(binum_array); - HypreFree(bjnum_array); - HypreFree(bii_array); - HypreFree(bij_array); - HypreFree(brhs_array); - HypreFree(aknum_array); - HypreFree(aki_array); - - return 0; } From 67660069a1bd2174bdc36a0750f8a9d79b0a88a9 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 5 Mar 2025 12:11:45 -0800 Subject: [PATCH 018/129] hypre_interface: Wrap BC arrays in structs & shared_ptr Use structs to group together arrays for boundary row elimination. Coefficients from HypreMatrix are stored in one struct, using a shared_ptr to share values with vectors. Eliminated boundary row values need to be passed from RHS to solution vector. This now uses a shared pointer that will manage the HypreFree call. --- include/bout/hypre_interface.hxx | 202 ++++++++++++++----------------- 1 file changed, 93 insertions(+), 109 deletions(-) diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index 0025ed062b..7ad50a008d 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -98,6 +98,82 @@ int checkHypreError(int error) { // TODO: set sizes // TODO: set contiguous blocks at once +/// Wrapper around HYPRE_Complex, that calls HypreFree when destroyed. +struct HypreComplexArray { + HYPRE_Complex *data; + + ~HypreComplexArray() { + HypreFree(data); + } +}; + +/// Shared pointter to a HypreComplexArray. When the last copy is destroyed +/// the HYPRE_Complex array inside will be free'd. +using BCValuesPtr = std::shared_ptr; + +/// Contains information needed to eliminate boundary equations from +/// RHS vectors, and restore boundary values in solution vectors. +struct BCMatrixEquations { + HYPRE_Int nb; + HYPRE_Int *binum_array; + HYPRE_Int *bjnum_array; + HYPRE_Complex *bii_array; + HYPRE_Complex *bij_array; + HYPRE_Int na; + HYPRE_Int *aknum_array; + HYPRE_Complex *aki_array; + + BCMatrixEquations() = delete; + + BCMatrixEquations(HYPRE_Int nrows, + HYPRE_Int *ncols, + HYPRE_BigInt *rows, + HYPRE_Int **row_indexes_ptr, + HYPRE_BigInt *cols, + HYPRE_Complex *values, + HYPRE_Int nb, // number of boundary equations + HYPRE_Int *bi_array) // row i for each boundary equation + : nb(nb) { + AdjustBCMatrixEquations(nrows, ncols, rows, row_indexes_ptr, cols, values, + nb, bi_array, + // Outputs + &binum_array, &bjnum_array, &bii_array, &bij_array, + &na, &aknum_array, &aki_array); + } + + ~BCMatrixEquations() { + // Free arrays + HypreFree(binum_array); + HypreFree(bjnum_array); + HypreFree(bii_array); + HypreFree(bij_array); + HypreFree(aknum_array); + HypreFree(aki_array); + } + + /// Applies in-place modification of the rhs array. + /// + /// Returns an array of boundary values that can be used to apply + /// boundary conditions to a solution vector. + BCValuesPtr adjustBCRightHandSideEquations(HYPRE_Complex *rhs) { + BCValuesPtr brhs = std::make_shared(); + AdjustBCRightHandSideEquations(rhs, nb, binum_array, bii_array, bij_array, &brhs->data, + na, aknum_array, aki_array); + return brhs; + } + + /// Apply boundary conditions to the solution. + /// Uses the BCValuesPtr returned from adjustBCRightHandSideEquations() + void adjustBCSolutionEquations(BCValuesPtr brhs, HYPRE_Complex *solution) { + AdjustBCSolutionEquations(solution, nb, + binum_array, bjnum_array, bii_array, bij_array, brhs->data); + } +}; + +/// A shared pointer to a BCMatrixEquations object +using BCMatrixPtr = std::shared_ptr; + + template class HypreVector { MPI_Comm comm; @@ -126,9 +202,6 @@ public: checkHypreError(HYPRE_IJVectorDestroy(hypre_vector)); HypreFree(I); HypreFree(V); - if (elimBErhs) { - HypreFree(brhs_array); - } } // Disable copy, at least for now: not clear that HYPRE_IJVector is @@ -217,22 +290,14 @@ public: HypreMalloc(V, vsize * sizeof(HYPRE_Complex)); } - // Data for eliminating boundary equations (TODO: convert this to a structure) + // Data for eliminating boundary equation bool elimBErhs = false; bool elimBEsol = false; - HYPRE_Complex *rhs; - HYPRE_Int nb; - HYPRE_Int *binum_array; - HYPRE_Int *bjnum_array; - HYPRE_Complex *bii_array; - HYPRE_Complex *bij_array; - HYPRE_Complex *brhs_array; // This is not created from the HypreMatrix class - HYPRE_Int na; - HYPRE_Int *aknum_array; - HYPRE_Complex *aki_array; + BCMatrixPtr bcmatrix; + BCValuesPtr bcvalues; /// Stores rhs values of BC rows void syncElimBErhs(HypreVector& rhs) { - brhs_array = rhs.brhs_array; + bcvalues = rhs.bcvalues; } void assemble() { @@ -245,8 +310,7 @@ public: void writeCacheToHypre() { if (elimBErhs) { - AdjustBCRightHandSideEquations(V, nb, binum_array, bii_array, bij_array, &brhs_array, - na, aknum_array, aki_array); + bcvalues = bcmatrix->adjustBCRightHandSideEquations(V); } checkHypreError(HYPRE_IJVectorSetValues(hypre_vector, vsize, I, V)); } @@ -254,7 +318,7 @@ public: void readCacheFromHypre() { checkHypreError(HYPRE_IJVectorGetValues(hypre_vector, vsize, I, V)); if (elimBEsol) { - AdjustBCSolutionEquations(V, nb, binum_array, bjnum_array, bii_array, bij_array, brhs_array); + bcmatrix->adjustBCSolutionEquations(bcvalues, V); } } @@ -410,19 +474,6 @@ public: using ind_type = typename T::ind_type; HypreMatrix() = default; -// The 'if' block below needs to be called when the HypreMatrix is destroyed, -// but I don't know where to put it. The HypreMatrix class is handled -// differently from HypreVector. -// ~HypreVector() { -// if (elimBE) { -// HypreFree(binum_array); -// HypreFree(bjnum_array); -// HypreFree(bii_array); -// HypreFree(bij_array); -// HypreFree(aknum_array); -// HypreFree(aki_array); -// } -// } HypreMatrix(const HypreMatrix&) = delete; HypreMatrix(HypreMatrix&& other) : comm(other.comm), ilower(other.ilower), iupper(other.iupper), @@ -748,16 +799,9 @@ public: return Element(*this, global_row, global_column, positions, weights); } - // Data for eliminating boundary equations (TODO: convert this to a structure) + // Data for eliminating boundary equations bool elimBE = false; - HYPRE_Int nb; - HYPRE_Int *binum_array; - HYPRE_Int *bjnum_array; - HYPRE_Complex *bii_array; - HYPRE_Complex *bij_array; - HYPRE_Int na; - HYPRE_Int *aknum_array; - HYPRE_Complex *aki_array; + BCMatrixPtr bcmatrix; // Shared pointer void setElimBE() { elimBE = true; @@ -765,24 +809,10 @@ public: void setElimBEVectors(HypreVector& sol, HypreVector& rhs) { sol.elimBEsol = elimBE; - sol.nb = nb; - sol.binum_array = binum_array; - sol.bjnum_array = bjnum_array; - sol.bii_array = bii_array; - sol.bij_array = bij_array; - sol.na = na; - sol.aknum_array = aknum_array; - sol.aki_array = aki_array; + sol.bcmatrix = bcmatrix; rhs.elimBErhs = elimBE; - rhs.nb = nb; - rhs.binum_array = binum_array; - rhs.bjnum_array = bjnum_array; - rhs.bii_array = bii_array; - rhs.bij_array = bij_array; - rhs.na = na; - rhs.aknum_array = aknum_array; - rhs.aki_array = aki_array; + rhs.bcmatrix = bcmatrix; } void assemble() { @@ -819,7 +849,7 @@ public: HYPRE_Int *bi_array; HYPRE_Int *row_indexes; // There must be an easier way to get nb - nb = 0; + int nb = 0; BOUT_FOR_SERIAL(i, index_converter->getRegionBndry()) { nb++; } @@ -829,59 +859,13 @@ public: bi_array[nb] = index_converter->getGlobal(i); nb++; } - AdjustBCMatrixEquations(num_rows, num_cols, rawI, &row_indexes, cols, vals, - nb, bi_array, - &binum_array, &bjnum_array, &bii_array, &bij_array, - &na, &aknum_array, &aki_array); + + bcmatrix = std::make_shared(num_rows, num_cols, + rawI, &row_indexes, + cols, vals, + nb, bi_array); HypreFree(bi_array); -// { -// FILE *file; -// int i; -// file = fopen("zbout.setvalues.out2", "w"); -// fprintf(file, "nrows %d\n", num_rows); -// for (i = 0; i < num_rows; i++) -// { -// fprintf(file, "ncols %d: %d\n", i, num_cols[i]); -// } -// for (i = 0; i < num_rows; i++) -// { -// fprintf(file, "rows %d: %d\n", i, rawI[i]); -// } -// for (i = 0; i < num_entries; i++) -// { -// fprintf(file, "cols %d: %d\n", i, cols[i]); -// } -// for (i = 0; i < num_entries; i++) -// { -// fprintf(file, "vals %d: %f\n", i, vals[i]); -// } -// fclose(file); -// -// file = fopen("zbout.setvalues.out3", "w"); -// for (i = 0; i < num_rows; i++) -// { -// fprintf(file, "row_indexes %d: %d\n", i, row_indexes[i]); -// } -// fprintf(file, "nb %d\n", nb); -// for (i = 0; i < nb; i++) -// { -// fprintf(file, "bijnum_array %d: %d %d\n", i, binum_array[i], bjnum_array[i]); -// } -// for (i = 0; i < nb; i++) -// { -// fprintf(file, "biiij_array %d: %f %f\n", i, bii_array[i], bij_array[i]); -// } -// fprintf(file, "na %d\n", na); -// for (i = 0; i < na; i++) -// { -// fprintf(file, "aknum_array %d: %d\n", i, aknum_array[i]); -// } -// for (i = 0; i < na; i++) -// { -// fprintf(file, "bki_array %d: %f\n", i, aki_array[i]); -// } -// fclose(file); -// } + checkHypreError( HYPRE_IJMatrixSetValues2(*hypre_matrix, num_rows, num_cols, rawI, row_indexes, cols, vals)); HypreFree(row_indexes); From 29b65c592a229b122da3660fe98a40a23d1f3ed5 Mon Sep 17 00:00:00 2001 From: Rob Falgout Date: Wed, 5 Mar 2025 12:40:04 -0800 Subject: [PATCH 019/129] Fixed a bug in rhs update for boundary elimination code --- include/bout/hypre_interface.hxx | 3 +-- src/invert/laplace/impls/hypre3d/hypre_boundary.c | 3 +-- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index 7ad50a008d..c8377f7f5a 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -44,7 +44,6 @@ AdjustBCRightHandSideEquations( HYPRE_Int nb, HYPRE_Int *binum_array, HYPRE_Complex *bii_array, - HYPRE_Complex *bij_array, HYPRE_Complex **brhs_array_ptr, HYPRE_Int na, HYPRE_Int *aknum_array, @@ -157,7 +156,7 @@ struct BCMatrixEquations { /// boundary conditions to a solution vector. BCValuesPtr adjustBCRightHandSideEquations(HYPRE_Complex *rhs) { BCValuesPtr brhs = std::make_shared(); - AdjustBCRightHandSideEquations(rhs, nb, binum_array, bii_array, bij_array, &brhs->data, + AdjustBCRightHandSideEquations(rhs, nb, binum_array, bii_array, &brhs->data, na, aknum_array, aki_array); return brhs; } diff --git a/src/invert/laplace/impls/hypre3d/hypre_boundary.c b/src/invert/laplace/impls/hypre3d/hypre_boundary.c index 69e9a2d1b0..30f0968765 100644 --- a/src/invert/laplace/impls/hypre3d/hypre_boundary.c +++ b/src/invert/laplace/impls/hypre3d/hypre_boundary.c @@ -174,7 +174,6 @@ AdjustBCRightHandSideEquations( HYPRE_Int nb, HYPRE_Int *binum_array, HYPRE_Complex *bii_array, - HYPRE_Complex *bij_array, HYPRE_Complex **brhs_array_ptr, HYPRE_Int na, HYPRE_Int *aknum_array, @@ -195,7 +194,7 @@ AdjustBCRightHandSideEquations( { bnum = anum; // Assume only one interior equation per boundary equation aknum = aknum_array[anum]; - rhs[aknum] -= aki_array[anum] * bij_array[bnum] / bii_array[bnum]; + rhs[aknum] -= aki_array[anum] * brhs_array[bnum] / bii_array[bnum]; } *brhs_array_ptr = brhs_array; From 62b923b0a40ebdd80b3a25ab0da4ebb76889c72a Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 5 Mar 2025 14:50:31 -0800 Subject: [PATCH 020/129] FCITransform: Save R and Z coordinates to output Reads cell center coordinates from the mesh, and saves to the output. This is to facilitate code coupling and visualisation. --- src/mesh/parallel/fci.cxx | 6 ++++++ src/mesh/parallel/fci.hxx | 13 ++++++++++++- 2 files changed, 18 insertions(+), 1 deletion(-) diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index acd4954f97..659df8600d 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -346,3 +346,9 @@ void FCITransform::integrateParallelSlices(Field3D& f) { f.ynext(map.offset) = map.integrate(f); } } + +void FCITransform::outputVars(Options& output_options) { + // Real-space coordinates of grid points + output_options["R"].force(R, "FCI"); + output_options["Z"].force(Z, "FCI"); +} diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx index 3ec3321a6a..1a02f558e1 100644 --- a/src/mesh/parallel/fci.hxx +++ b/src/mesh/parallel/fci.hxx @@ -73,11 +73,15 @@ public: FCITransform() = delete; FCITransform(Mesh& mesh, const Coordinates::FieldMetric& dy, bool zperiodic = true, Options* opt = nullptr) - : ParallelTransform(mesh, opt) { + : ParallelTransform(mesh, opt), R{&mesh}, Z{&mesh} { // check the coordinate system used for the grid data source FCITransform::checkInputGrid(); + // Real-space coordinates of grid cells + mesh.get(R, "R", 0.0, false); + mesh.get(Z, "Z", 0.0, false); + auto forward_boundary_xin = std::make_shared("FCI_forward", BNDRY_PAR_FWD_XIN, +1, &mesh); auto backward_boundary_xin = std::make_shared( @@ -142,6 +146,10 @@ public: bool canToFromFieldAligned() const override { return false; } + /// Save mesh variables to output + /// If R and Z(x,y,z) coordinates are in the input then these are saved to output. + void outputVars(Options& output_options) override; + bool requiresTwistShift(bool UNUSED(twist_shift_enabled), [[maybe_unused]] YDirectionType ytype) override { // No Field3Ds require twist-shift, because they cannot be field-aligned @@ -156,6 +164,9 @@ protected: private: /// FCI maps for each of the parallel slices std::vector field_line_maps; + + /// Real-space coordinates of grid points + Field3D R, Z; }; #endif // BOUT_FCITRANSFORM_H From 91ddbfd2af7cf85b442c70b1df250d0f282a9c0c Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 5 Mar 2025 14:52:03 -0800 Subject: [PATCH 021/129] Mesh, BoutMesh, Coordinates: update header comment Update copyright statement, contact email. --- include/bout/coordinates.hxx | 4 ++-- include/bout/mesh.hxx | 4 ++-- src/mesh/impls/bout/boutmesh.cxx | 13 ++----------- 3 files changed, 6 insertions(+), 15 deletions(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index fb2fade75d..908c052e73 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -9,9 +9,9 @@ * * ************************************************************************** - * Copyright 2014 B.D.Dudson + * Copyright 2014-2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index a1c88a2634..fc65dd8e0c 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -17,9 +17,9 @@ * * Incorporates code from topology.cpp and Communicator * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010-2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index 31319b87d1..9214c56855 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -2,19 +2,10 @@ * Implementation of the Mesh class, handling input files compatible with * BOUT / BOUT-06. * - * Changelog - * --------- - * - * 2015-01 Ben Dudson - * * - * - * 2010-05 Ben Dudson - * * Initial version, adapted from grid.cpp and topology.cpp - * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010-2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * From 33d78feb1b3a764e08be6156ff7894d6430c9b1b Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 5 Mar 2025 16:05:32 -0800 Subject: [PATCH 022/129] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- examples/elm-pb/elm_pb.cxx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/elm-pb/elm_pb.cxx b/examples/elm-pb/elm_pb.cxx index b36bd728d4..1a2608e62c 100644 --- a/examples/elm-pb/elm_pb.cxx +++ b/examples/elm-pb/elm_pb.cxx @@ -1364,7 +1364,7 @@ class ELMpb : public PhysicsModel { // Only update if simulation time has advanced // Uses an exponential decay of the weighting of the value in the boundary // so that the solution is well behaved for arbitrary steps - BoutReal weight = exp(-(t - phi_boundary_last_update) / phi_boundary_timescale); + BoutReal const weight = exp(-(t - phi_boundary_last_update) / phi_boundary_timescale); phi_boundary_last_update = t; if (mesh->firstX()) { @@ -1378,7 +1378,7 @@ class ELMpb : public PhysicsModel { } } MPI_Comm comm_inner = mesh->getYcomm(0); - int np; + int np = 0; MPI_Comm_size(comm_inner, &np); MPI_Allreduce(&philocal, &phivalue, 1, MPI_DOUBLE, MPI_SUM, comm_inner); phivalue /= (np * mesh->LocalNz * mesh->LocalNy); @@ -1393,11 +1393,11 @@ class ELMpb : public PhysicsModel { } // Old value of phi at boundary. Note: this is constant in Z - BoutReal oldvalue = + BoutReal const oldvalue = 0.5 * (phi(mesh->xstart - 1, j, 0) + phi(mesh->xstart, j, 0)); // New value of phi at boundary, relaxing towards phivalue - BoutReal newvalue = weight * oldvalue + (1. - weight) * phivalue; + BoutReal const newvalue = weight * oldvalue + (1. - weight) * phivalue; // Set phi at the boundary to this value for (int k = 0; k < mesh->LocalNz; k++) { @@ -1420,7 +1420,7 @@ class ELMpb : public PhysicsModel { 0.5 * (phi(mesh->xend + 1, j, 0) + phi(mesh->xend, j, 0)); // New value of phi at boundary, relaxing towards phivalue - BoutReal newvalue = weight * oldvalue + (1. - weight) * phivalue; + BoutReal const newvalue = weight * oldvalue + (1. - weight) * phivalue; // Set phi at the boundary to this value for (int k = 0; k < mesh->LocalNz; k++) { From 08b1261f86e38eb45864265a4d48ec644545256b Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 6 Mar 2025 09:43:05 -0800 Subject: [PATCH 023/129] hypre_interface: Moved BC elimination to src/sys/hypre_interface.cxx Changed C free functions to be members of BCMatrixEquations. Clang formatted, minor style changes. Test results unchanged. --- CMakeLists.txt | 1 + include/bout/hypre_interface.hxx | 187 +++++++-------- .../laplace/impls/hypre3d/hypre3d_laplace.cxx | 1 - .../laplace/impls/hypre3d/hypre_boundary.c | 222 ------------------ src/sys/hypre_interface.cxx | 131 +++++++++++ 5 files changed, 215 insertions(+), 327 deletions(-) delete mode 100644 src/invert/laplace/impls/hypre3d/hypre_boundary.c create mode 100644 src/sys/hypre_interface.cxx diff --git a/CMakeLists.txt b/CMakeLists.txt index c45fca3b72..af628a7db5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -335,6 +335,7 @@ set(BOUT_SOURCES ./src/sys/generator_context.cxx ./include/bout/hyprelib.hxx ./src/sys/hyprelib.cxx + ./src/sys/hypre_interface.cxx ./src/sys/msg_stack.cxx ./src/sys/options.cxx ./src/sys/options/optionparser.hxx diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index c8377f7f5a..33e6218c8f 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -20,45 +20,6 @@ #include "HYPRE_utilities.h" #include "_hypre_utilities.h" -void -AdjustBCMatrixEquations( - HYPRE_Int nrows, - HYPRE_Int *ncols, - HYPRE_BigInt *rows, - HYPRE_Int **row_indexes_ptr, - HYPRE_BigInt *cols, - HYPRE_Complex *values, - HYPRE_Int nb, // number of boundary equations - HYPRE_Int *bi_array, // row i for each boundary equation - HYPRE_Int **binum_array_ptr, // data index for row i (for each boundary equation) - HYPRE_Int **bjnum_array_ptr, // data index for col j (for each boundary equation) - HYPRE_Complex **bii_array_ptr, // coefficient b_ii (for each boundary equation) - HYPRE_Complex **bij_array_ptr, // coefficient b_ij (for each boundary equation) - HYPRE_Int *na_ptr, // number of interior equations to adjust - HYPRE_Int **aknum_array_ptr, // data index for row k (for each interior equation) - HYPRE_Complex **aki_array_ptr); // coefficient a_ki (for each interior equation) - -void -AdjustBCRightHandSideEquations( - HYPRE_Complex *rhs, - HYPRE_Int nb, - HYPRE_Int *binum_array, - HYPRE_Complex *bii_array, - HYPRE_Complex **brhs_array_ptr, - HYPRE_Int na, - HYPRE_Int *aknum_array, - HYPRE_Complex *aki_array); - -void -AdjustBCSolutionEquations( - HYPRE_Complex *solution, - HYPRE_Int nb, - HYPRE_Int *binum_array, - HYPRE_Int *bjnum_array, - HYPRE_Complex *bii_array, - HYPRE_Complex *bij_array, - HYPRE_Complex *brhs_array); - #include // BOUT_ENUM_CLASS does not work inside namespaces @@ -99,46 +60,83 @@ int checkHypreError(int error) { /// Wrapper around HYPRE_Complex, that calls HypreFree when destroyed. struct HypreComplexArray { - HYPRE_Complex *data; + HYPRE_Complex* data; - ~HypreComplexArray() { - HypreFree(data); - } + HypreComplexArray(int n) { HypreMalloc(data, sizeof(HYPRE_Complex) * n); } + + ~HypreComplexArray() { HypreFree(data); } }; /// Shared pointter to a HypreComplexArray. When the last copy is destroyed /// the HYPRE_Complex array inside will be free'd. using BCValuesPtr = std::shared_ptr; -/// Contains information needed to eliminate boundary equations from -/// RHS vectors, and restore boundary values in solution vectors. +/*! + * This function modifies the input for the HYPRE_IJMatrixSetValues() routine to + * eliminate the boundary condition equations (see below for details on how the + * equations are adjusted). It modifies the arrays ncols, rows, cols, and + * values. It also returns a row_indexes array. This can then be passed to the + * HYPRE_IJMatrixSetValues2() routine to set up the matrix in hypre. + * + * The arguments nb and bi_array indicate the boundary equations. The routine + * returns info needed to adjust the right-hand-side and solution vector through + * the functions AdjustRightHandSideEquations and AdjustSolutionEquations. + * + * NOTE: It may make sense from an organizational standpoint to collect many of + * these arguments in a structure of some sort. + * + * Notation, assumptions, and other details: + * + * - Boundary equation i is assumed to have two coefficients + * + * b_ii * u_i + b_ij * u_j = rhs_i + * + * - We also assume that each boundary equation has only one interior equation k + * coupled to it (such that k = j) with coupling coefficient a_ki + * + * a_ki * u_i + a_kj * u_j + ... = rhs_k + * + * - Each equation k is adjusted as follows: + * + * a_kj = a_kj - a_ki * b_ij / b_ii + * a_ki = 0 + * + * - Boundary equations are adjusted to be identity equations in the matrix, but + * the boundary coefficients (b_ii, b_ij) are returned for use later + * + * - Right-hand-side equations are adjusted in AdjustRightHandSideEquations() as + * follows: rhs_k = rhs_k - a_ki * rhs_i / b_ii + * + * - Solution unknowns are adjusted at boundaries in AdjustSolutionEquations as + * follows: u_i = (rhs_i - b_ij * u_j) / b_ii + * + * - Naming conventions: Arrays starting with 'b' are boundary equation arrays + * indexed by 'bnum', and arrays starting with 'a' are non-boundary arrays + * (interior matrix equations) indexed by 'anum'. When 'num' is prefixed with + * a row or column number 'i', 'j', or 'k', the array holds the corresponding + * local data index for that row or column (e.g., an index into the local + * solution vector). Matrix coefficients are named as above, e.g., 'bij' is + * the coefficient for b_ij. + * + * NOTE: Implementation in src/sys/hypre_interface.cxx + */ struct BCMatrixEquations { - HYPRE_Int nb; - HYPRE_Int *binum_array; - HYPRE_Int *bjnum_array; - HYPRE_Complex *bii_array; - HYPRE_Complex *bij_array; - HYPRE_Int na; - HYPRE_Int *aknum_array; - HYPRE_Complex *aki_array; + HYPRE_Int nb; + HYPRE_Int* binum_array; + HYPRE_Int* bjnum_array; + HYPRE_Complex* bii_array; + HYPRE_Complex* bij_array; + HYPRE_Int na; + HYPRE_Int* aknum_array; + HYPRE_Complex* aki_array; BCMatrixEquations() = delete; - BCMatrixEquations(HYPRE_Int nrows, - HYPRE_Int *ncols, - HYPRE_BigInt *rows, - HYPRE_Int **row_indexes_ptr, - HYPRE_BigInt *cols, - HYPRE_Complex *values, - HYPRE_Int nb, // number of boundary equations - HYPRE_Int *bi_array) // row i for each boundary equation - : nb(nb) { - AdjustBCMatrixEquations(nrows, ncols, rows, row_indexes_ptr, cols, values, - nb, bi_array, - // Outputs - &binum_array, &bjnum_array, &bii_array, &bij_array, - &na, &aknum_array, &aki_array); - } + BCMatrixEquations(HYPRE_Int nrows, HYPRE_Int* ncols, HYPRE_BigInt* rows, + HYPRE_Int** row_indexes_ptr, HYPRE_BigInt* cols, + HYPRE_Complex* values, + HYPRE_Int nb, // number of boundary equations + HYPRE_Int* bi_array); // row i for each boundary equation ~BCMatrixEquations() { // Free arrays @@ -154,25 +152,16 @@ struct BCMatrixEquations { /// /// Returns an array of boundary values that can be used to apply /// boundary conditions to a solution vector. - BCValuesPtr adjustBCRightHandSideEquations(HYPRE_Complex *rhs) { - BCValuesPtr brhs = std::make_shared(); - AdjustBCRightHandSideEquations(rhs, nb, binum_array, bii_array, &brhs->data, - na, aknum_array, aki_array); - return brhs; - } + BCValuesPtr adjustBCRightHandSideEquations(HYPRE_Complex* rhs); /// Apply boundary conditions to the solution. /// Uses the BCValuesPtr returned from adjustBCRightHandSideEquations() - void adjustBCSolutionEquations(BCValuesPtr brhs, HYPRE_Complex *solution) { - AdjustBCSolutionEquations(solution, nb, - binum_array, bjnum_array, bii_array, bij_array, brhs->data); - } + void adjustBCSolutionEquations(BCValuesPtr brhs, HYPRE_Complex* solution); }; /// A shared pointer to a BCMatrixEquations object using BCMatrixPtr = std::shared_ptr; - template class HypreVector { MPI_Comm comm; @@ -293,11 +282,9 @@ public: bool elimBErhs = false; bool elimBEsol = false; BCMatrixPtr bcmatrix; - BCValuesPtr bcvalues; /// Stores rhs values of BC rows + BCValuesPtr bcvalues; /// Stores rhs values of BC rows - void syncElimBErhs(HypreVector& rhs) { - bcvalues = rhs.bcvalues; - } + void syncElimBErhs(HypreVector& rhs) { bcvalues = rhs.bcvalues; } void assemble() { CALI_CXX_MARK_FUNCTION; @@ -802,9 +789,7 @@ public: bool elimBE = false; BCMatrixPtr bcmatrix; // Shared pointer - void setElimBE() { - elimBE = true; - } + void setElimBE() { elimBE = true; } void setElimBEVectors(HypreVector& sol, HypreVector& rhs) { sol.elimBEsol = elimBE; @@ -845,34 +830,28 @@ public: // Eliminate boundary condition equations in hypre SetValues input arguments if (elimBE) { - HYPRE_Int *bi_array; - HYPRE_Int *row_indexes; + HYPRE_Int* bi_array; + HYPRE_Int* row_indexes; // There must be an easier way to get nb int nb = 0; - BOUT_FOR_SERIAL(i, index_converter->getRegionBndry()) { - nb++; - } + BOUT_FOR_SERIAL(i, index_converter->getRegionBndry()) { nb++; } HypreMalloc(bi_array, nb * sizeof(HYPRE_Int)); nb = 0; BOUT_FOR_SERIAL(i, index_converter->getRegionBndry()) { - bi_array[nb] = index_converter->getGlobal(i); - nb++; + bi_array[nb] = index_converter->getGlobal(i); + nb++; } - bcmatrix = std::make_shared(num_rows, num_cols, - rawI, &row_indexes, - cols, vals, - nb, bi_array); + bcmatrix = std::make_shared( + num_rows, num_cols, rawI, &row_indexes, cols, vals, nb, bi_array); HypreFree(bi_array); - checkHypreError( - HYPRE_IJMatrixSetValues2(*hypre_matrix, num_rows, num_cols, rawI, row_indexes, cols, vals)); + checkHypreError(HYPRE_IJMatrixSetValues2(*hypre_matrix, num_rows, num_cols, rawI, + row_indexes, cols, vals)); HypreFree(row_indexes); - } - else - { + } else { checkHypreError( - HYPRE_IJMatrixSetValues(*hypre_matrix, num_rows, num_cols, rawI, cols, vals)); + HYPRE_IJMatrixSetValues(*hypre_matrix, num_rows, num_cols, rawI, cols, vals)); } checkHypreError(HYPRE_IJMatrixAssemble(*hypre_matrix)); checkHypreError(HYPRE_IJMatrixGetObject(*hypre_matrix, @@ -1055,7 +1034,7 @@ public: switch (solver_type) { case HYPRE_SOLVER_TYPE::gmres: { - HYPRE_ParCSRGMRESSetKDim(solver, 30); // TODO: Make this an input file parameter + HYPRE_ParCSRGMRESSetKDim(solver, 30); // TODO: Make this an input file parameter break; } case HYPRE_SOLVER_TYPE::bicgstab: { diff --git a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx index def54b36b5..beb83a216d 100644 --- a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx +++ b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx @@ -28,7 +28,6 @@ #if BOUT_HAS_HYPRE #include "hypre3d_laplace.hxx" -#include "hypre_boundary.c" #include #include diff --git a/src/invert/laplace/impls/hypre3d/hypre_boundary.c b/src/invert/laplace/impls/hypre3d/hypre_boundary.c deleted file mode 100644 index 30f0968765..0000000000 --- a/src/invert/laplace/impls/hypre3d/hypre_boundary.c +++ /dev/null @@ -1,222 +0,0 @@ - -/* - * This function modifies the input for the HYPRE_IJMatrixSetValues() routine to - * eliminate the boundary condition equations (see below for details on how the - * equations are adjusted). It modifies the arrays ncols, rows, cols, and - * values. It also returns a row_indexes array. This can then be passed to the - * HYPRE_IJMatrixSetValues2() routine to set up the matrix in hypre. - * - * The arguments nb and bi_array indicate the boundary equations. The routine - * returns info needed to adjust the right-hand-side and solution vector through - * the functions AdjustRightHandSideEquations and AdjustSolutionEquations. - * - * NOTE: It may make sense from an organizational standpoint to collect many of - * these arguments in a structure of some sort. - * - * Notation, assumptions, and other details: - * - * - Boundary equation i is assumed to have two coefficients - * - * b_ii * u_i + b_ij * u_j = rhs_i - * - * - We also assume that each boundary equation has only one interior equation k - * coupled to it (such that k = j) with coupling coefficient a_ki - * - * a_ki * u_i + a_kj * u_j + ... = rhs_k - * - * - Each equation k is adjusted as follows: - * - * a_kj = a_kj - a_ki * b_ij / b_ii - * a_ki = 0 - * - * - Boundary equations are adjusted to be identity equations in the matrix, but - * the boundary coefficients (b_ii, b_ij) are returned for use later - * - * - Right-hand-side equations are adjusted in AdjustRightHandSideEquations() as - * follows: rhs_k = rhs_k - a_ki * rhs_i / b_ii - * - * - Solution unknowns are adjusted at boundaries in AdjustSolutionEquations as - * follows: u_i = (rhs_i - b_ij * u_j) / b_ii - * - * - Naming conventions: Arrays starting with 'b' are boundary equation arrays - * indexed by 'bnum', and arrays starting with 'a' are non-boundary arrays - * (interior matrix equations) indexed by 'anum'. When 'num' is prefixed with - * a row or column number 'i', 'j', or 'k', the array holds the corresponding - * local data index for that row or column (e.g., an index into the local - * solution vector). Matrix coefficients are named as above, e.g., 'bij' is - * the coefficient for b_ij. - */ - -void -AdjustBCMatrixEquations( - HYPRE_Int nrows, - HYPRE_Int *ncols, - HYPRE_BigInt *rows, - HYPRE_Int **row_indexes_ptr, - HYPRE_BigInt *cols, - HYPRE_Complex *values, - HYPRE_Int nb, // number of boundary equations - HYPRE_Int *bi_array, // row i for each boundary equation - HYPRE_Int **binum_array_ptr, // data index for row i (for each boundary equation) - HYPRE_Int **bjnum_array_ptr, // data index for col j (for each boundary equation) - HYPRE_Complex **bii_array_ptr, // coefficient b_ii (for each boundary equation) - HYPRE_Complex **bij_array_ptr, // coefficient b_ij (for each boundary equation) - HYPRE_Int *na_ptr, // number of interior equations to adjust - HYPRE_Int **aknum_array_ptr, // data index for row k (for each interior equation) - HYPRE_Complex **aki_array_ptr) // coefficient a_ki (for each interior equation) -{ - HYPRE_Int *row_indexes; - HYPRE_Int na, *binum_array, *bjnum_array, *aknum_array; - HYPRE_Complex *bii_array, *bij_array, *aki_array; - HYPRE_Int i, j, k, m, mkj, anum, bnum, acoeffnum, bcoeffnum; - HYPRE_Int binum, aknum; - HYPRE_Complex bii, bij, aki; - - /* Create the row_indexes array */ - row_indexes = (HYPRE_Int *)malloc(sizeof(HYPRE_Int) * nrows); - row_indexes[0] = 0; - for (i = 1; i < nrows; i++) - { - row_indexes[i] = row_indexes[i-1] + ncols[i-1]; - } - - /* Assume just one interior equation coupled to each boundary equation */ - na = nb; - - /* Allocate return arrays */ - HypreMalloc(binum_array, sizeof(HYPRE_Int) * nb); - HypreMalloc(bjnum_array, sizeof(HYPRE_Int) * nb); - HypreMalloc(bii_array, sizeof(HYPRE_Complex) * nb); - HypreMalloc(bij_array, sizeof(HYPRE_Complex) * nb); - HypreMalloc(aknum_array, sizeof(HYPRE_Int) * na); - HypreMalloc(aki_array, sizeof(HYPRE_Complex) * na); - - binum = 0; - aknum = 0; - for (bnum = 0; bnum < nb; bnum++) - { - /* Get boundary equation information and adjust boundary equations */ - /* Find row i in rows array (assume i increases and rows is sorted) */ - i = bi_array[bnum]; - for (; binum < nrows; binum++) - { - if (i == rows[binum]) - { - break; // Found row i in rows array - } - } - bcoeffnum = row_indexes[binum]; - for (m = 0; m < 2; m++) // Assume only two boundary equation coefficients - { - if (cols[bcoeffnum + m] == i) - { - bii = values[bcoeffnum + m]; - values[bcoeffnum + m] = -1.0; // Identity equation (negative definite matrix) - } - else - { - j = cols[bcoeffnum + m]; - bij = values[bcoeffnum + m]; - values[bcoeffnum + m] = 0.0; // Identity equation - } - } - ncols[binum] = 1; // Identity equation - - /* Get interior equation information and adjust interior equations */ - /* Find row k in rows array (assume k increases and rows is sorted) */ - k = j; // Assume equation k = j - for (; aknum < nrows; aknum++) - { - if (k == rows[aknum]) - { - break; // Found row k in rows array - } - } - acoeffnum = row_indexes[aknum]; - for (m = 0; m < ncols[aknum]; m++) - { - if (cols[acoeffnum + m] == j) - { - mkj = m; // Save for update of akj value below - } - if (cols[acoeffnum + m] == i) - { - aki = values[acoeffnum + m]; - values[acoeffnum + m] = 0.0; // Eliminate coupling to boundary equation - } - } - values[acoeffnum + mkj] -= aki * bij / bii; // Update akj value - - /* Update return arrays */ - anum = bnum; // Assume only one interior equation k - binum_array[bnum] = binum; - bjnum_array[bnum] = aknum; // Assume only one interior equation k - bii_array[bnum] = bii; - bij_array[bnum] = bij; - aknum_array[anum] = aknum; - aki_array[anum] = aki; - } - - /* Set return arguments */ - *row_indexes_ptr = row_indexes; - *binum_array_ptr = binum_array; - *bjnum_array_ptr = bjnum_array; - *bii_array_ptr = bii_array; - *bij_array_ptr = bij_array; - *na_ptr = na; - *aknum_array_ptr = aknum_array; - *aki_array_ptr = aki_array; -} - -void -AdjustBCRightHandSideEquations( - HYPRE_Complex *rhs, - HYPRE_Int nb, - HYPRE_Int *binum_array, - HYPRE_Complex *bii_array, - HYPRE_Complex **brhs_array_ptr, - HYPRE_Int na, - HYPRE_Int *aknum_array, - HYPRE_Complex *aki_array) -{ - HYPRE_Complex *brhs_array; - HYPRE_Int anum, bnum, binum, aknum; - - HypreMalloc(brhs_array, sizeof(HYPRE_Complex) * nb); - - for (bnum = 0; bnum < nb; bnum++) - { - binum = binum_array[bnum]; - brhs_array[bnum] = rhs[binum]; - } - - for (anum = 0; anum < na; anum++) - { - bnum = anum; // Assume only one interior equation per boundary equation - aknum = aknum_array[anum]; - rhs[aknum] -= aki_array[anum] * brhs_array[bnum] / bii_array[bnum]; - } - - *brhs_array_ptr = brhs_array; -} - -void -AdjustBCSolutionEquations( - HYPRE_Complex *solution, - HYPRE_Int nb, - HYPRE_Int *binum_array, - HYPRE_Int *bjnum_array, - HYPRE_Complex *bii_array, - HYPRE_Complex *bij_array, - HYPRE_Complex *brhs_array) -{ - HYPRE_Int bnum, binum, bjnum; - - for (bnum = 0; bnum < nb; bnum++) - { - binum = binum_array[bnum]; - bjnum = bjnum_array[bnum]; - solution[binum] = (brhs_array[bnum] - bij_array[bnum] * solution[bjnum]) / bii_array[bnum]; - } -} - diff --git a/src/sys/hypre_interface.cxx b/src/sys/hypre_interface.cxx new file mode 100644 index 0000000000..1838bd5158 --- /dev/null +++ b/src/sys/hypre_interface.cxx @@ -0,0 +1,131 @@ + +#include "bout/build_defines.hxx" + +#if BOUT_HAS_HYPRE + +#include "bout/hypre_interface.hxx" + +namespace bout { + +BCMatrixEquations::BCMatrixEquations(HYPRE_Int nrows, HYPRE_Int* ncols, + HYPRE_BigInt* rows, HYPRE_Int** row_indexes_ptr, + HYPRE_BigInt* cols, HYPRE_Complex* values, + HYPRE_Int nb, HYPRE_Int* bi_array) + : nb(nb) { + HYPRE_Int* row_indexes; + + // Create the row_indexes array + row_indexes = (HYPRE_Int*)malloc(sizeof(HYPRE_Int) * nrows); + row_indexes[0] = 0; + for (HYPRE_Int i = 1; i < nrows; i++) { + row_indexes[i] = row_indexes[i - 1] + ncols[i - 1]; + } + + // Assume just one interior equation coupled to each boundary equation + na = nb; + + // Allocate arrays + HypreMalloc(binum_array, sizeof(HYPRE_Int) * nb); + HypreMalloc(bjnum_array, sizeof(HYPRE_Int) * nb); + HypreMalloc(bii_array, sizeof(HYPRE_Complex) * nb); + HypreMalloc(bij_array, sizeof(HYPRE_Complex) * nb); + HypreMalloc(aknum_array, sizeof(HYPRE_Int) * na); + HypreMalloc(aki_array, sizeof(HYPRE_Complex) * na); + + HYPRE_Int binum = 0; + HYPRE_Int aknum = 0; + for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { + // Get boundary equation information and adjust boundary equations + // Find row i in rows array (assume i increases and rows is sorted) + HYPRE_Int i = bi_array[bnum]; + for (; binum < nrows; binum++) { + if (i == rows[binum]) { + break; // Found row i in rows array + } + } + HYPRE_Int bcoeffnum = row_indexes[binum]; + HYPRE_Complex bii{0.0}, bij{0.0}; + HYPRE_Int j = 0; + + for (HYPRE_Int m = 0; m < 2; m++) { // Assume only two boundary equation coefficients + if (cols[bcoeffnum + m] == i) { + bii = values[bcoeffnum + m]; + values[bcoeffnum + m] = -1.0; // Identity equation (negative definite matrix) + } else { + j = cols[bcoeffnum + m]; + bij = values[bcoeffnum + m]; + values[bcoeffnum + m] = 0.0; // Identity equation + } + } + ncols[binum] = 1; // Identity equation + + /* Get interior equation information and adjust interior equations */ + /* Find row k in rows array (assume k increases and rows is sorted) */ + HYPRE_Int k = j; // Assume equation k = j + for (; aknum < nrows; aknum++) { + if (k == rows[aknum]) { + break; // Found row k in rows array + } + } + HYPRE_Int acoeffnum = row_indexes[aknum]; + + HYPRE_Int mkj = 0; + HYPRE_Complex aki{0.0}; + for (HYPRE_Int m = 0; m < ncols[aknum]; m++) { + if (cols[acoeffnum + m] == j) { + mkj = m; // Save for update of akj value below + } + if (cols[acoeffnum + m] == i) { + aki = values[acoeffnum + m]; + values[acoeffnum + m] = 0.0; // Eliminate coupling to boundary equation + } + } + values[acoeffnum + mkj] -= aki * bij / bii; // Update akj value + + // Update arrays + HYPRE_Int anum = bnum; // Assume only one interior equation k + binum_array[bnum] = binum; + bjnum_array[bnum] = aknum; // Assume only one interior equation k + bii_array[bnum] = bii; + bij_array[bnum] = bij; + aknum_array[anum] = aknum; + aki_array[anum] = aki; + } + + // Set return arguments + *row_indexes_ptr = row_indexes; +} + +BCValuesPtr BCMatrixEquations::adjustBCRightHandSideEquations(HYPRE_Complex* rhs) { + + // Allocate array to store boundary row values + BCValuesPtr brhs = std::make_shared(nb); + + for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { + HYPRE_Int binum = binum_array[bnum]; + brhs->data[bnum] = rhs[binum]; + } + + for (HYPRE_Int anum = 0; anum < na; anum++) { + HYPRE_Int bnum = anum; // Assume only one interior equation per boundary equation + HYPRE_Int aknum = aknum_array[anum]; + rhs[aknum] -= aki_array[anum] * brhs->data[bnum] / bii_array[bnum]; + } + + return brhs; +} + +void BCMatrixEquations::adjustBCSolutionEquations(BCValuesPtr brhs, + HYPRE_Complex* solution) { + + for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { + HYPRE_Int binum = binum_array[bnum]; + HYPRE_Int bjnum = bjnum_array[bnum]; + solution[binum] = + (brhs->data[bnum] - bij_array[bnum] * solution[bjnum]) / bii_array[bnum]; + } +} + +} // namespace bout + +#endif // BOUT_HAS_HYPRE From eea57fc2c1a798a80eca0d7687fa8219ffe09f9a Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 6 Mar 2025 09:44:53 -0800 Subject: [PATCH 024/129] test-laplace-hypre3d: Suppress BOUT++ log outputs Difficult to see the test results amongst the log outputs. Now pipe the logs and save to files. --- tests/integrated/test-laplace-hypre3d/runtest | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/tests/integrated/test-laplace-hypre3d/runtest b/tests/integrated/test-laplace-hypre3d/runtest index b50c5993b7..f1f5950547 100755 --- a/tests/integrated/test-laplace-hypre3d/runtest +++ b/tests/integrated/test-laplace-hypre3d/runtest @@ -19,20 +19,24 @@ build_and_log("Laplace 3D with Hypre") success = True for directory, nproc in test_directories: - command = "test-laplace3d -d " + directory + command = "./test-laplace3d -d " + directory print("running on", nproc, "processors:", command) - launch_safe(command, nproc=nproc) + s, out = launch_safe(command, nproc=nproc, pipe=True) + # Save output to log file + with open("run.log." + directory, "w") as f: + f.write(out) error_max = collect("error_max", path=directory, info=False) if error_max > tolerance: - print(directory + " failed with maximum error {}".format(error_max)) + print(" => " + directory + " failed with maximum error {}".format(error_max)) success = False else: - print(directory + " passed with maximum error {}".format(error_max)) + print(" => " + directory + " passed with maximum error {}".format(error_max)) if success: - print("All passed") + print("=> All passed") exit(0) else: + print("=> Some tests failed") exit(1) From ef9fbed3e68ff6ea45294e95cb174457cc5c939b Mon Sep 17 00:00:00 2001 From: bendudson <219233+bendudson@users.noreply.github.com> Date: Thu, 6 Mar 2025 17:46:47 +0000 Subject: [PATCH 025/129] Apply black changes --- tests/integrated/test_suite | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/integrated/test_suite b/tests/integrated/test_suite index 307a8d84b3..77ad7882c4 100755 --- a/tests/integrated/test_suite +++ b/tests/integrated/test_suite @@ -188,7 +188,7 @@ class Test(threading.Thread): self.output += "\n(It is likely that a timeout occured)" else: # ❌ Failed - print("\u274C", end="") # No newline + print("\u274c", end="") # No newline print(" %7.3f s" % (time.time() - self.local.start_time), flush=True) def _cost(self): From af72bd8f83e4ba5770d1800dd25bede21f6f1724 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 6 Mar 2025 19:52:51 +0100 Subject: [PATCH 026/129] Do not fail if clang-format is formatting the code The CI job failed so far, if it had to format the code, as the return code is only zero if nothing needs to be done. --- .github/workflows/clang-format.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/clang-format.yml b/.github/workflows/clang-format.yml index 87f1947802..a6508a2dcd 100644 --- a/.github/workflows/clang-format.yml +++ b/.github/workflows/clang-format.yml @@ -22,7 +22,7 @@ jobs: - name: Run clang-format id: format - run: git clang-format origin/${{ github.base_ref }} + run: git clang-format origin/${{ github.base_ref }} || : - name: Commit to the PR branch uses: stefanzweifel/git-auto-commit-action@v5 From aeb83c6c18cc5fd59e9149d3a507185ab87ad453 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 12 Mar 2025 09:16:05 +0100 Subject: [PATCH 027/129] CI: Avoid issues with special characters --- .github/workflows/clang-format.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/clang-format.yml b/.github/workflows/clang-format.yml index a6508a2dcd..d99b370810 100644 --- a/.github/workflows/clang-format.yml +++ b/.github/workflows/clang-format.yml @@ -22,7 +22,7 @@ jobs: - name: Run clang-format id: format - run: git clang-format origin/${{ github.base_ref }} || : + run: 'git clang-format origin/${{ github.base_ref }} || :' - name: Commit to the PR branch uses: stefanzweifel/git-auto-commit-action@v5 From dfa31b2870962ff86c76679ec5a16f358a4186b1 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 12 Mar 2025 09:50:40 +0100 Subject: [PATCH 028/129] CI: run git-clang-format until there are no more changes That might format more code at once, but should avoid a CI loop. --- .github/workflows/clang-format.yml | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/.github/workflows/clang-format.yml b/.github/workflows/clang-format.yml index d99b370810..5e25154300 100644 --- a/.github/workflows/clang-format.yml +++ b/.github/workflows/clang-format.yml @@ -22,7 +22,11 @@ jobs: - name: Run clang-format id: format - run: 'git clang-format origin/${{ github.base_ref }} || :' + run: + while ! git clang-format origin/${{ github.base_ref }} + do + true + done - name: Commit to the PR branch uses: stefanzweifel/git-auto-commit-action@v5 From e743eb06f5c38e2d8feefc9c652268bbc6b7480f Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 12 Mar 2025 11:06:31 +0100 Subject: [PATCH 029/129] CI: use one line --- .github/workflows/clang-format.yml | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/.github/workflows/clang-format.yml b/.github/workflows/clang-format.yml index 5e25154300..3b8cf6ee50 100644 --- a/.github/workflows/clang-format.yml +++ b/.github/workflows/clang-format.yml @@ -23,10 +23,7 @@ jobs: - name: Run clang-format id: format run: - while ! git clang-format origin/${{ github.base_ref }} - do - true - done + while ! git clang-format origin/${{ github.base_ref }} ; do true ; done - name: Commit to the PR branch uses: stefanzweifel/git-auto-commit-action@v5 From 1cb5277a6d288cc2f2f89c1ca840ede7e11faf52 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 12 Mar 2025 11:16:07 +0100 Subject: [PATCH 030/129] CI: stage before we run git-clang-format again --- .github/workflows/clang-format.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/clang-format.yml b/.github/workflows/clang-format.yml index 3b8cf6ee50..49dfb31e25 100644 --- a/.github/workflows/clang-format.yml +++ b/.github/workflows/clang-format.yml @@ -23,7 +23,7 @@ jobs: - name: Run clang-format id: format run: - while ! git clang-format origin/${{ github.base_ref }} ; do true ; done + while ! git clang-format origin/${{ github.base_ref }} ; do git add . ; done - name: Commit to the PR branch uses: stefanzweifel/git-auto-commit-action@v5 From 9d1b46908fab0c61950905eefa73ad735b4a4c27 Mon Sep 17 00:00:00 2001 From: Rob Falgout Date: Wed, 19 Mar 2025 14:39:19 -0700 Subject: [PATCH 031/129] Turn off test for false convergence in hypre GMRES solver --- include/bout/hypre_interface.hxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index 33e6218c8f..1837d5e275 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -1034,7 +1034,8 @@ public: switch (solver_type) { case HYPRE_SOLVER_TYPE::gmres: { - HYPRE_ParCSRGMRESSetKDim(solver, 30); // TODO: Make this an input file parameter + HYPRE_ParCSRGMRESSetKDim(solver, 30); // TODO: Make this an input file parameter + HYPRE_GMRESSetSkipRealResidualCheck(solver, 1); // TODO: Make this an input file parameter break; } case HYPRE_SOLVER_TYPE::bicgstab: { From f969c011249b60b337817432c026b1bd312291cc Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 19 Mar 2025 14:43:14 -0700 Subject: [PATCH 032/129] hypre_interface: Add kdim and skip_real_residual_check options These options control the GMRES solver. --- include/bout/hypre_interface.hxx | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index 1837d5e275..15ca97cb8a 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -1034,8 +1034,16 @@ public: switch (solver_type) { case HYPRE_SOLVER_TYPE::gmres: { - HYPRE_ParCSRGMRESSetKDim(solver, 30); // TODO: Make this an input file parameter - HYPRE_GMRESSetSkipRealResidualCheck(solver, 1); // TODO: Make this an input file parameter + HYPRE_ParCSRGMRESSetKDim(solver, + options["kdim"] + .doc("Set the maximum size of the Krylov space") + .withDefault(30)); + + if (options["skip_real_residual_check"] + .doc("Skip the evaluation and the check of the actual residual?") + .withDefault(false)) { + HYPRE_GMRESSetSkipRealResidualCheck(solver, 1); + } break; } case HYPRE_SOLVER_TYPE::bicgstab: { From 3aed489c6cbb56137260860c63c4d32b9bd8029b Mon Sep 17 00:00:00 2001 From: Norbert Podhorszki Date: Fri, 26 Apr 2024 09:41:34 -0400 Subject: [PATCH 033/129] Add attributes to ADIOS2 output to "define" dimensions as names. We need this for xarray support that is designed for the NetCDF model, where dimensions are separately defined entities. --- include/bout/adios_object.hxx | 8 +++++++- src/sys/options/options_adios.cxx | 24 ++++++++++++++++++------ 2 files changed, 25 insertions(+), 7 deletions(-) diff --git a/include/bout/adios_object.hxx b/include/bout/adios_object.hxx index b14316f1ba..12da8ce4c7 100755 --- a/include/bout/adios_object.hxx +++ b/include/bout/adios_object.hxx @@ -57,11 +57,17 @@ public: } template - adios2::Variable GetArrayVariable(const std::string& varname, adios2::Dims& shape) { + adios2::Variable GetArrayVariable(const std::string& varname, adios2::Dims& shape, + const std::vector& dimNames, + int rank) { adios2::Variable v = io.InquireVariable(varname); if (!v) { adios2::Dims start(shape.size()); v = io.DefineVariable(varname, shape, start, shape); + if (!rank && dimNames.size()) { + io.DefineAttribute("__xarray_dimensions__", dimNames.data(), + dimNames.size(), varname, "/", true); + } } else { v.SetShape(shape); } diff --git a/src/sys/options/options_adios.cxx b/src/sys/options/options_adios.cxx index b3acbaada6..09797647ed 100644 --- a/src/sys/options/options_adios.cxx +++ b/src/sys/options/options_adios.cxx @@ -307,6 +307,12 @@ void OptionsADIOS::verifyTimesteps() const { return; } +const std::vector DIMS_NONE; +const std::vector DIMS_X = {"x"}; +const std::vector DIMS_XY = {"x", "y"}; +const std::vector DIMS_XZ = {"x", "z"}; +const std::vector DIMS_XYZ = {"x", "y", "z"}; + /// Visit a variant type, and put the data into a NcVar struct ADIOSPutVarVisitor { ADIOSPutVarVisitor(const std::string& name, ADIOSStream& stream) @@ -388,7 +394,8 @@ void ADIOSPutVarVisitor::operator()(const Field2D& value) { adios2::Dims memCount = {static_cast(value.getNx()), static_cast(value.getNy())}; - adios2::Variable var = stream.GetArrayVariable(varname, shape); + adios2::Variable var = + stream.GetArrayVariable(varname, shape, DIMS_XY, BoutComm::rank()); var.SetSelection({start, count}); var.SetMemorySelection({memStart, memCount}); stream.engine.Put(var, &value(0, 0)); @@ -425,7 +432,8 @@ void ADIOSPutVarVisitor::operator()(const Field3D& value) { static_cast(value.getNy()), static_cast(value.getNz())}; - adios2::Variable var = stream.GetArrayVariable(varname, shape); + adios2::Variable var = + stream.GetArrayVariable(varname, shape, DIMS_XYZ, BoutComm::rank()); var.SetSelection({start, count}); var.SetMemorySelection({memStart, memCount}); stream.engine.Put(var, &value(0, 0, 0)); @@ -457,7 +465,8 @@ void ADIOSPutVarVisitor::operator()(const FieldPerp& value) { adios2::Dims memCount = {static_cast(value.getNx()), static_cast(value.getNz())}; - adios2::Variable var = stream.GetArrayVariable(varname, shape); + adios2::Variable var = + stream.GetArrayVariable(varname, shape, DIMS_XZ, BoutComm::rank()); var.SetSelection({start, count}); var.SetMemorySelection({memStart, memCount}); stream.engine.Put(var, &value(0, 0)); @@ -469,7 +478,8 @@ void ADIOSPutVarVisitor::operator()>(const Array& valu adios2::Dims shape = {(size_t)BoutComm::size(), (size_t)value.size()}; adios2::Dims start = {(size_t)BoutComm::rank(), 0}; adios2::Dims count = {1, shape[1]}; - adios2::Variable var = stream.GetArrayVariable(varname, shape); + adios2::Variable var = + stream.GetArrayVariable(varname, shape, DIMS_NONE, BoutComm::rank()); var.SetSelection({start, count}); stream.engine.Put(var, value.begin()); } @@ -482,7 +492,8 @@ void ADIOSPutVarVisitor::operator()>(const Matrix& va (size_t)std::get<1>(s)}; adios2::Dims start = {(size_t)BoutComm::rank(), 0, 0}; adios2::Dims count = {1, shape[1], shape[2]}; - adios2::Variable var = stream.GetArrayVariable(varname, shape); + adios2::Variable var = + stream.GetArrayVariable(varname, shape, DIMS_NONE, BoutComm::rank()); var.SetSelection({start, count}); stream.engine.Put(var, value.begin()); } @@ -495,7 +506,8 @@ void ADIOSPutVarVisitor::operator()>(const Tensor& va (size_t)std::get<1>(s), (size_t)std::get<2>(s)}; adios2::Dims start = {(size_t)BoutComm::rank(), 0, 0, 0}; adios2::Dims count = {1, shape[1], shape[2], shape[3]}; - adios2::Variable var = stream.GetArrayVariable(varname, shape); + adios2::Variable var = + stream.GetArrayVariable(varname, shape, DIMS_NONE, BoutComm::rank()); var.SetSelection({start, count}); stream.engine.Put(var, value.begin()); } From d43cea2ee40a654898a4b4e4801aeb4934bb800d Mon Sep 17 00:00:00 2001 From: Seimon Powell Date: Thu, 24 Apr 2025 14:54:20 +0100 Subject: [PATCH 034/129] Enable alternative stencil shapes in jocobian colouring 4 new options added solver:stencil: = N where