1515#include " bout/field2d.hxx"
1616#include " bout/globalindexer.hxx"
1717#include " bout/hyprelib.hxx"
18+ #include " bout/msg_stack.hxx"
1819#include " bout/options.hxx"
1920#include " bout/output.hxx"
2021#include " bout/paralleltransform.hxx"
2122#include " bout/region.hxx"
23+ #include " bout/unused.hxx"
2224#include " bout/utils.hxx"
2325
2426#include " HYPRE.h"
2527#include " HYPRE_IJ_mv.h"
2628#include " HYPRE_parcsr_ls.h"
2729#include " HYPRE_parcsr_mv.h"
2830#include " HYPRE_utilities.h"
29- #include " _hypre_utilities.h"
3031
3132#include < algorithm>
3233#include < cmath>
3334#include < iterator>
3435#include < memory>
36+ #include < utility>
3537#include < vector>
3638
3739// BOUT_ENUM_CLASS does not work inside namespaces
@@ -48,7 +50,7 @@ namespace bout {
4850
4951namespace {
5052int checkHypreError (int error) {
51- if (error) {
53+ if (error != 0 ) {
5254 // Aaron Fisher warns that it is possible that the error code is non-zero even
5355 // though a solve was successful. If this occurs, we might want to make this a
5456 // warning, or otherwise investigate further.
@@ -217,7 +219,7 @@ class HypreVector {
217219 HYPRE_Complex* V{nullptr };
218220 HYPRE_Complex* workV{nullptr };
219221 bool cache_current{false };
220- HypreLib hyprelib{} ;
222+ HypreLib hyprelib;
221223
222224public:
223225 static_assert (bout::utils::is_Field_v<T>, " HypreVector only works with Fields" );
@@ -241,7 +243,7 @@ public:
241243 HypreVector (const HypreVector<T>&) = delete ;
242244 auto operator =(const HypreVector<T>&) = delete ;
243245
244- HypreVector (HypreVector<T>&& other)
246+ HypreVector (HypreVector<T>&& other) noexcept
245247 : comm(other.comm), jlower(other.jlower), jupper(other.jupper), vsize(other.vsize),
246248 indexConverter (other.indexConverter), location(other.location),
247249 initialised(other.initialised), have_indices(other.have_indices), I(other.I),
@@ -312,7 +314,8 @@ public:
312314 }
313315
314316 // / Construct a vector with given index set, but don't set any values
315- explicit HypreVector (IndexerPtr<T> indConverter) : indexConverter(indConverter) {
317+ explicit HypreVector (IndexerPtr<T> indConverter)
318+ : indexConverter(indConverter), location(CELL_LOC ::centre) {
316319 Mesh& mesh = *indConverter->getMesh ();
317320 const MPI_Comm comm =
318321 std::is_same_v<T, FieldPerp> ? mesh.getXZcomm () : BoutComm::get ();
@@ -325,7 +328,6 @@ public:
325328 checkHypreError (HYPRE_IJVectorSetObjectType (hypre_vector, HYPRE_PARCSR ));
326329 checkHypreError (HYPRE_IJVectorInitialize (hypre_vector));
327330 initialised = true ;
328- location = CELL_LOC ::centre;
329331 HypreMalloc (I, vsize * sizeof (HYPRE_BigInt));
330332 HypreMalloc (V, vsize * sizeof (HYPRE_Complex));
331333 HypreMalloc (workV, vsize * sizeof (HYPRE_Complex));
@@ -546,13 +548,13 @@ class HypreMatrix {
546548 ParallelTransform* parallel_transform{nullptr };
547549 bool assembled{false };
548550 HYPRE_BigInt num_rows;
549- std::vector<HYPRE_BigInt>* I;
550- std::vector<std::vector<HYPRE_BigInt>>* J;
551- std::vector<std::vector<HYPRE_Complex>>* V;
551+ std::shared_ptr<std:: vector<HYPRE_BigInt>> I;
552+ std::shared_ptr<std:: vector<std::vector<HYPRE_BigInt>>> J;
553+ std::shared_ptr<std:: vector<std::vector<HYPRE_Complex>>> V;
552554 // / Enable reduction of boundary equations before assembling the Hypre matrix.
553555 bool use_boundary_elimination{false };
554556 BoundaryEliminationPtr boundary_elimination;
555- HypreLib hyprelib{} ;
557+ HypreLib hyprelib;
556558
557559 // todo also take care of I,J,V
558560 struct MatrixDeleter {
@@ -575,8 +577,8 @@ public:
575577 index_converter (other.index_converter), location(other.location),
576578 initialised(other.initialised), yoffset(other.yoffset),
577579 parallel_transform(other.parallel_transform), assembled(other.assembled),
578- num_rows(other.num_rows), I(other.I), J(other.J), V(other.V ),
579- use_boundary_elimination(other.use_boundary_elimination),
580+ num_rows(other.num_rows), I(std::move( other.I)) , J(std::move( other.J)),
581+ V(std::move(other.V)), use_boundary_elimination(other.use_boundary_elimination),
580582 boundary_elimination(other.boundary_elimination) {
581583 std::swap (hypre_matrix, other.hypre_matrix );
582584 std::swap (parallel_matrix, other.parallel_matrix );
@@ -595,9 +597,9 @@ public:
595597 parallel_transform = other.parallel_transform ;
596598 assembled = other.assembled ;
597599 num_rows = other.num_rows ;
598- I = other.I ;
599- J = other.J ;
600- V = other.V ;
600+ I = std::move ( other.I ) ;
601+ J = std::move ( other.J ) ;
602+ V = std::move ( other.V ) ;
601603 use_boundary_elimination = other.use_boundary_elimination ;
602604 boundary_elimination = other.boundary_elimination ;
603605 return *this ;
@@ -608,27 +610,21 @@ public:
608610 // /
609611 // / note: preallocate not currently used, but here to match PetscMatrix interface
610612 explicit HypreMatrix (IndexerPtr<T> indConverter, bool UNUSED (preallocate) = true)
611- : hypre_matrix(new HYPRE_IJMatrix, MatrixDeleter{}), index_converter(indConverter) {
613+ : ilower(indConverter->getGlobalStart ()),
614+ iupper(ilower + indConverter->size () - 1), // inclusive end
615+ hypre_matrix(new HYPRE_IJMatrix, MatrixDeleter{}), index_converter(indConverter),
616+ num_rows(iupper - ilower + 1 ),
617+ I(std::make_shared<std::vector<HYPRE_BigInt>>(num_rows)),
618+ J(std::make_shared<std::vector<std::vector<HYPRE_BigInt>>>(num_rows)),
619+ V(std::make_shared<std::vector<std::vector<HYPRE_Complex>>>(num_rows)) {
612620 Mesh* mesh = indConverter->getMesh ();
613621 const MPI_Comm comm =
614622 std::is_same_v<T, FieldPerp> ? mesh->getXZcomm () : BoutComm::get ();
615623 parallel_transform = &mesh->getCoordinates ()->getParallelTransform ();
616624
617- ilower = indConverter->getGlobalStart ();
618- iupper = ilower + indConverter->size () - 1 ; // inclusive end
619- num_rows = iupper - ilower + 1 ;
620-
621- I = new std::vector<HYPRE_BigInt>;
622- J = new std::vector<std::vector<HYPRE_BigInt>>;
623- V = new std::vector<std::vector<HYPRE_Complex>>;
624- (*I).resize (num_rows);
625- (*J).resize (num_rows);
626- (*V).resize (num_rows);
627625 for (HYPRE_BigInt i = 0 ; i < num_rows; ++i) {
628626 (*I)[i] = ilower + i;
629- (*J)[i].resize (0 );
630627 (*J)[i].reserve (10 );
631- (*V)[i].resize (0 );
632628 (*V)[i].reserve (10 );
633629 }
634630
@@ -659,8 +655,8 @@ public:
659655 Element (HypreMatrix<T>& matrix_, HYPRE_BigInt row_, HYPRE_BigInt column_,
660656 std::vector<HYPRE_BigInt> positions_ = {},
661657 std::vector<BoutReal> weights_ = {})
662- : hypre_matrix(matrix_.get()), row(row_), column(column_), positions(positions_),
663- weights (weights_) {
658+ : hypre_matrix(matrix_.get()), row(row_), column(column_),
659+ positions (std::move(positions_)), weights(std::move( weights_) ) {
664660#if CHECK > 2
665661 for (const auto val : weights) {
666662 ASSERT3 (std::isfinite (val));
@@ -1032,7 +1028,7 @@ public:
10321028 const HYPRE_Complex operator_value = elimination_state.interior_values ->data [i];
10331029 elimination_state.interior_values ->data [i] = alpha * operator_value;
10341030 elimination_state.boundary_values ->data [i] =
1035- alpha * operator_value
1031+ ( alpha * operator_value)
10361032 + (y_boundary_values != nullptr ? beta * y_boundary_values->data [i] : 0.0 );
10371033 }
10381034 }
@@ -1073,7 +1069,7 @@ private:
10731069 HYPRE_Solver precon;
10741070 bool solver_setup{false };
10751071 HYPRE_SOLVER_TYPE solver_type = HYPRE_SOLVER_TYPE ::gmres;
1076- HypreLib hyprelib{} ;
1072+ HypreLib hyprelib;
10771073
10781074 // Basically a manual vtable for the various implementations of the
10791075 // HYPRE solver methods. By storing function pointers and setting
@@ -1092,14 +1088,13 @@ private:
10921088 decltype (HYPRE_ParCSRGMRESSolve)* solverSolve{nullptr };
10931089
10941090public:
1095- HypreSystem (Mesh& mesh, Options& options) {
1091+ HypreSystem (Mesh& mesh, Options& options)
1092+ : comm(std::is_same_v<T, FieldPerp> ? mesh.getXZcomm() : BoutComm::get()) {
10961093 solver_type = options[" hypre_solver_type" ]
10971094 .doc (" Type of solver to use when solving Hypre system. Possible "
10981095 " values are: gmres, bicgstab, pcg" )
10991096 .withDefault (HYPRE_SOLVER_TYPE ::bicgstab);
11001097
1101- comm = std::is_same_v<T, FieldPerp> ? mesh.getXZcomm () : BoutComm::get ();
1102-
11031098 auto print_level =
11041099 options[" hypre_print_level" ]
11051100 .doc (" Verbosity for Hypre solver. Integer from 0 (silent) to 4 "
0 commit comments