diff --git a/cajita/src/CMakeLists.txt b/cajita/src/CMakeLists.txt index 547697312..1d2b7dd93 100644 --- a/cajita/src/CMakeLists.txt +++ b/cajita/src/CMakeLists.txt @@ -44,7 +44,9 @@ set(HEADERS_PUBLIC if(Cabana_ENABLE_HYPRE) list(APPEND HEADERS_PUBLIC + Cajita_Hypre.hpp Cajita_HypreStructuredSolver.hpp + Cajita_HypreSemiStructuredSolver.hpp ) endif() diff --git a/cajita/src/Cajita.hpp b/cajita/src/Cajita.hpp index 870fd57e1..8262d343d 100644 --- a/cajita/src/Cajita.hpp +++ b/cajita/src/Cajita.hpp @@ -46,6 +46,8 @@ #include #ifdef Cabana_ENABLE_HYPRE +#include +#include #include #endif diff --git a/cajita/src/Cajita_Hypre.hpp b/cajita/src/Cajita_Hypre.hpp new file mode 100644 index 000000000..d54b467a3 --- /dev/null +++ b/cajita/src/Cajita_Hypre.hpp @@ -0,0 +1,85 @@ +/**************************************************************************** + * Copyright (c) 2018-2022 by the Cabana authors * + * All rights reserved. * + * * + * This file is part of the Cabana library. Cabana is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +/*! + \file Cajita_Hypre.hpp + \brief HYPRE memory space handling +*/ +#ifndef CAJITA_HYPRE_HPP +#define CAJITA_HYPRE_HPP + +#include +#include +#include + +#include + +#include + +namespace Cajita +{ +//---------------------------------------------------------------------------// +// Hypre memory space selection. Don't compile if HYPRE wasn't configured to +// use the device. +// ---------------------------------------------------------------------------// + +//! Hypre device compatibility check. +template +struct HypreIsCompatibleWithMemorySpace : std::false_type +{ +}; + +// FIXME: This is currently written in this structure because HYPRE only has +// compile-time switches for backends and hence only one can be used at a +// time. Once they have a run-time switch we can use that instead. +#ifdef HYPRE_USING_CUDA +#ifdef KOKKOS_ENABLE_CUDA +#ifdef HYPRE_USING_DEVICE_MEMORY +//! Hypre device compatibility check - CUDA memory. +template <> +struct HypreIsCompatibleWithMemorySpace : std::true_type +{ +}; +#endif // end HYPRE_USING_DEVICE_MEMORY + +//! Hypre device compatibility check - CUDA UVM memory. +#ifdef HYPRE_USING_UNIFIED_MEMORY +template <> +struct HypreIsCompatibleWithMemorySpace : std::true_type +{ +}; +#endif // end HYPRE_USING_UNIFIED_MEMORY +#endif // end KOKKOS_ENABLE_CUDA +#endif // end HYPRE_USING_CUDA + +#ifdef HYPRE_USING_HIP +#ifdef KOKKOS_ENABLE_HIP +//! Hypre device compatibility check - HIP memory. FIXME - make this true when +//! the HYPRE CMake includes HIP +template <> +struct HypreIsCompatibleWithMemorySpace + : std::false_type +{ +}; +#endif // end KOKKOS_ENABLE_HIP +#endif // end HYPRE_USING_HIP + +#ifndef HYPRE_USING_GPU +//! Hypre device compatibility check - host memory. +template <> +struct HypreIsCompatibleWithMemorySpace : std::true_type +{ +}; +#endif // end HYPRE_USING_GPU + +} // namespace Cajita + +#endif // end HYPRE_HPP diff --git a/cajita/src/Cajita_HypreSemiStructuredSolver.hpp b/cajita/src/Cajita_HypreSemiStructuredSolver.hpp new file mode 100644 index 000000000..26f8fdb55 --- /dev/null +++ b/cajita/src/Cajita_HypreSemiStructuredSolver.hpp @@ -0,0 +1,1230 @@ +/**************************************************************************** + * Copyright (c) 2018-2022 by the Cabana authors * + * All rights reserved. * + * * + * This file is part of the Cabana library. Cabana is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +/*! + \file Cajita_HypreSemiStructuredSolver.hpp + \brief HYPRE semi-structured solver interface +*/ +#ifndef CAJITA_HYPRESEMISTRUCTUREDSOLVER_HPP +#define CAJITA_HYPRESEMISTRUCTUREDSOLVER_HPP + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace Cajita +{ + +//---------------------------------------------------------------------------// +//! Hypre semi-structured solver interface for scalar fields. +template +class HypreSemiStructuredSolver +{ + public: + //! Entity type. + using entity_type = EntityType; + //! Kokkos memory space.. + using memory_space = MemorySpace; + //! Scalar value type. + using value_type = Scalar; + //! Object Type for SStruct + const int object_type = HYPRE_SSTRUCT; + //! Hypre memory space compatibility check. + static_assert( HypreIsCompatibleWithMemorySpace::value, + "HYPRE not compatible with solver memory space" ); + + /*! + \brief Constructor. + \param layout The array layout defining the vector space of the solver. + \param n_vars Number of variable in the system domain + \param is_preconditioner Flag indicating if this solver will be used as + a preconditioner. + */ + template + HypreSemiStructuredSolver( const ArrayLayout_t& layout, int n_vars, + const bool is_preconditioner = false ) + : _comm( layout.localGrid()->globalGrid().comm() ) + , _is_preconditioner( is_preconditioner ) + { + HYPRE_Init(); + + static_assert( is_array_layout::value, + "Must use an array layout" ); + static_assert( + std::is_same::value, + "Array layout entity type mush match solver entity type" ); + + // Spatial dimension. + const std::size_t num_space_dim = ArrayLayout_t::num_space_dim; + + // Only a single part grid is supported initially + int n_parts = 1; + int part = 0; + + // Only create data structures if this is not a preconditioner. + if ( !_is_preconditioner ) + { + // Create the grid. + auto error = HYPRE_SStructGridCreate( _comm, num_space_dim, n_parts, + &_grid ); + checkHypreError( error ); + + // Get the global index space spanned by the local grid on this + // rank. Note that the upper bound is not a bound but rather the + // last index as this is what Hypre wants. Note that we reordered + // this to KJI from IJK to be consistent with HYPRE ordering. By + // setting up the grid like this, HYPRE will then want layout-right + // data indexed as (i,j,k) or (i,j,k,l) which will allow us to + // directly use Kokkos::deep_copy to move data between Cajita arrays + // and HYPRE data structures. + auto global_space = layout.indexSpace( Own(), Global() ); + _lower.resize( num_space_dim ); + _upper.resize( num_space_dim ); + for ( std::size_t d = 0; d < num_space_dim; ++d ) + { + _lower[d] = static_cast( + global_space.min( num_space_dim - d - 1 ) ); + _upper[d] = static_cast( + global_space.max( num_space_dim - d - 1 ) - 1 ); + } + error = HYPRE_SStructGridSetExtents( _grid, part, _lower.data(), + _upper.data() ); + checkHypreError( error ); + + // Get periodicity. Note we invert the order of this to KJI as well. + const auto& global_grid = layout.localGrid()->globalGrid(); + HYPRE_Int periodic[num_space_dim]; + for ( std::size_t d = 0; d < num_space_dim; ++d ) + periodic[num_space_dim - 1 - d] = + global_grid.isPeriodic( d ) + ? layout.localGrid()->globalGrid().globalNumEntity( + EntityType(), d ) + : 0; + error = HYPRE_SStructGridSetPeriodic( _grid, part, periodic ); + checkHypreError( error ); + + // Set the variables on the HYPRE grid + std::vector vartypes; + vartypes.resize( n_vars ); + for ( int i = 0; i < n_vars; ++i ) + { + vartypes[i] = HYPRE_SSTRUCT_VARIABLE_CELL; + } + error = HYPRE_SStructGridSetVariables( _grid, part, n_vars, + vartypes.data() ); + + // Assemble the grid. + error = HYPRE_SStructGridAssemble( _grid ); + checkHypreError( error ); + + // Allocate LHS and RHS vectors and initialize to zero. Note that we + // are fixing the views under these vectors to layout-right. + + std::array reorder_size; + for ( std::size_t d = 0; d < num_space_dim; ++d ) + { + reorder_size[d] = global_space.extent( d ); + } + reorder_size.back() = n_vars; + IndexSpace reorder_space( reorder_size ); + auto vector_values = + createView( + "vector_values0", reorder_space ); + Kokkos::deep_copy( vector_values, 0.0 ); + + _stencils.resize( n_vars ); + _stencil_size.resize( n_vars ); + _stencil_index.resize( n_vars, + std::vector( n_vars + 1 ) ); + + error = HYPRE_SStructVectorCreate( _comm, _grid, &_b ); + checkHypreError( error ); + error = HYPRE_SStructVectorSetObjectType( _b, object_type ); + checkHypreError( error ); + error = HYPRE_SStructVectorInitialize( _b ); + checkHypreError( error ); + for ( int i = 0; i < n_vars; ++i ) + { + error = HYPRE_SStructVectorSetBoxValues( + _b, part, _lower.data(), _upper.data(), i, + vector_values.data() ); + checkHypreError( error ); + } + error = HYPRE_SStructVectorAssemble( _b ); + checkHypreError( error ); + + error = HYPRE_SStructVectorCreate( _comm, _grid, &_x ); + checkHypreError( error ); + error = HYPRE_SStructVectorSetObjectType( _x, object_type ); + checkHypreError( error ); + error = HYPRE_SStructVectorInitialize( _x ); + checkHypreError( error ); + for ( int i = 0; i < n_vars; ++i ) + { + error = HYPRE_SStructVectorSetBoxValues( + _x, part, _lower.data(), _upper.data(), i, + vector_values.data() ); + checkHypreError( error ); + } + checkHypreError( error ); + error = HYPRE_SStructVectorAssemble( _x ); + checkHypreError( error ); + } + } + + // Destructor. + virtual ~HypreSemiStructuredSolver() + { + // We only make data if this is not a preconditioner. + if ( !_is_preconditioner ) + { + HYPRE_SStructVectorDestroy( _x ); + HYPRE_SStructVectorDestroy( _b ); + HYPRE_SStructMatrixDestroy( _A ); + for ( std::size_t i = 0; i < _stencils.size(); ++i ) + { + HYPRE_SStructStencilDestroy( _stencils[i] ); + } + HYPRE_SStructGridDestroy( _grid ); + HYPRE_SStructGraphDestroy( _graph ); + + HYPRE_Finalize(); + } + } + + //! Return if this solver is a preconditioner. + bool isPreconditioner() const { return _is_preconditioner; } + + /*! + \brief Create the operator stencil to be filled by setMatrixStencil + \param NumSpaceDim The number of spatial dimensions in the linear system + being solved. + \param var The variable number that the stencil corresponds to, in essence + which equation number in the linear system + \param n_vars number of variables in the linear system + \param stencil_length A vector containing the length of the stencil for + variable `var` for each variable in the system to be created for HYPRE + */ + void createMatrixStencil( int NumSpaceDim, int var = 0, int n_vars = 3, + std::vector stencil_length = { 7, 7, + 7 } ) + { + // This function is only valid for non-preconditioners. + if ( _is_preconditioner ) + throw std::logic_error( + "Cannot call createMatrixStencil() on preconditioners" ); + + // Generate the stencil indexing + unsigned index = 0; + for ( int i = 0; i < n_vars; ++i ) + { + _stencil_index[var][i] = index; + index += stencil_length[i]; + } + _stencil_index[var][n_vars] = index; + + // Create the stencil. + _stencil_size[var] = index; + auto error = HYPRE_SStructStencilCreate( + NumSpaceDim, _stencil_size[var], &_stencils[var] ); + checkHypreError( error ); + } + + /*! + \brief Set the operator stencil. + \param stencil The (i,j,k) offsets describing the structured matrix + entries at each grid point. Offsets are defined relative to an index. + \param var The variable number that the stencil corresponds to, in essence + which equation number in the linear system + \param dep The integer for the independent variable in the linear system + that is currently being set + */ + template + void + setMatrixStencil( const std::vector>& stencil, + int var = 0, int dep = 0 ) + { + // This function is only valid for non-preconditioners. + if ( _is_preconditioner ) + throw std::logic_error( + "Cannot call setMatrixStencil() on preconditioners" ); + + std::array offset; + + auto index = _stencil_index[var][dep]; + for ( unsigned n = index; n < index + stencil.size(); ++n ) + { + for ( std::size_t d = 0; d < NumSpaceDim; ++d ) + offset[d] = stencil[n - index][d]; + auto error = HYPRE_SStructStencilSetEntry( _stencils[var], n, + offset.data(), dep ); + checkHypreError( error ); + } + } + + /*! + \brief Set the solver graph + \param n_vars The number of variables (and equations) in the + specified problem domain + */ + void setSolverGraph( int n_vars ) + { + // This function is only valid for non-preconditioners. + if ( _is_preconditioner ) + throw std::logic_error( + "Cannot call setSolverGraph() on preconditioners" ); + + int part = 0; + + // Setup the Graph for the non-zero structure of the matrix + // Create the graph with hypre + auto error = HYPRE_SStructGraphCreate( _comm, _grid, &_graph ); + checkHypreError( error ); + + // Set up the object type + error = HYPRE_SStructGraphSetObjectType( _graph, object_type ); + checkHypreError( error ); + + // Set the stencil to the graph + for ( int i = 0; i < n_vars; ++i ) + { + error = + HYPRE_SStructGraphSetStencil( _graph, part, i, _stencils[i] ); + checkHypreError( error ); + } + + // Assemble the graph + error = HYPRE_SStructGraphAssemble( _graph ); + checkHypreError( error ); + + // Create the matrix. + error = HYPRE_SStructMatrixCreate( _comm, _graph, &_A ); + checkHypreError( error ); + + // Set the SStruct matrix object type + error = HYPRE_SStructMatrixSetObjectType( _A, object_type ); + checkHypreError( error ); + } + + /*! + \brief Prepare the hypre matrix to have it's values set + */ + void initializeHypreMatrix() + { + // Initialize the matrix. + auto error = HYPRE_SStructMatrixInitialize( _A ); + checkHypreError( error ); + } + + /*! + \brief Set the matrix values. + \param values The matrix entry values. For each entity over which the + vector space is defined an entry for each stencil element is + required. The order of the stencil elements is that same as that in the + stencil definition. Note that values corresponding to stencil entries + outside of the domain should be set to zero. + \param v_x The variable index for the independent variable (column) + that is being set by the current call to setMatrixValues + \param v_h The variable index for the equation (row) that is being + set by the current call to setMatrixValues + */ + template + void setMatrixValues( const Array_t& values, int v_x, int v_h ) + { + static_assert( is_array::value, "Must use an array" ); + static_assert( + std::is_same::value, + "Array entity type mush match solver entity type" ); + static_assert( + std::is_same::value, + "Array device type and solver device type are different." ); + + static_assert( + std::is_same::value, + "Array value type and solver value type are different." ); + + // This function is only valid for non-preconditioners. + if ( _is_preconditioner ) + throw std::logic_error( + "Cannot call setMatrixValues() on preconditioners" ); + + int index_size = + _stencil_index[v_h][v_x + 1] - _stencil_index[v_h][v_x]; + + // Ensure the values array matches up in dimension with the stencil size + if ( values.layout()->dofsPerEntity() != + static_cast( index_size ) ) + throw std::runtime_error( + "Number of matrix values does not match stencil size" ); + + // Spatial dimension. + const std::size_t num_space_dim = Array_t::num_space_dim; + + int part = 0; + + // Copy the matrix entries into HYPRE. The HYPRE layout is fixed as + // layout-right. + auto owned_space = values.layout()->indexSpace( Own(), Local() ); + std::array reorder_size; + for ( std::size_t d = 0; d < num_space_dim; ++d ) + { + reorder_size[d] = owned_space.extent( d ); + } + + reorder_size.back() = index_size; + IndexSpace reorder_space( reorder_size ); + auto a_values = + createView( + "a_values", reorder_space ); + + auto values_subv = createSubview( values.view(), owned_space ); + Kokkos::deep_copy( a_values, values_subv ); + + // Insert values into the HYPRE matrix. + std::vector indices( index_size ); + int start = _stencil_index[v_h][v_x]; + std::iota( indices.begin(), indices.end(), start ); + auto error = HYPRE_SStructMatrixSetBoxValues( + _A, part, _lower.data(), _upper.data(), v_h, indices.size(), + indices.data(), a_values.data() ); + checkHypreError( error ); + } + + /*! + \brief Print the hypre matrix to ouput file + \param prefix File prefix for where hypre output is written + */ + void printMatrix( const char* prefix ) + { + HYPRE_SStructMatrixPrint( prefix, _A, 0 ); + } + + /*! + \brief Print the hypre LHS to ouput file + \param prefix File prefix for where hypre output is written + */ + void printLHS( const char* prefix ) + { + HYPRE_SStructVectorPrint( prefix, _x, 0 ); + } + + /*! + \brief Print the hypre RHS to ouput file + \param prefix File prefix for where hypre output is written + */ + void printRHS( const char* prefix ) + { + HYPRE_SStructVectorPrint( prefix, _b, 0 ); + } + + //! Set convergence tolerance implementation. + void setTolerance( const double tol ) { this->setToleranceImpl( tol ); } + + //! Set maximum iteration implementation. + void setMaxIter( const int max_iter ) { this->setMaxIterImpl( max_iter ); } + + //! Set the output level. + void setPrintLevel( const int print_level ) + { + this->setPrintLevelImpl( print_level ); + } + + //! Set a preconditioner. + void + setPreconditioner( const std::shared_ptr>& preconditioner ) + { + // This function is only valid for non-preconditioners. + if ( _is_preconditioner ) + throw std::logic_error( + "Cannot call setPreconditioner() on a preconditioner" ); + + // Only a preconditioner can be used as a preconditioner. + if ( !preconditioner->isPreconditioner() ) + throw std::logic_error( "Not a preconditioner" ); + + _preconditioner = preconditioner; + this->setPreconditionerImpl( *_preconditioner ); + } + + //! Setup the problem. + void setup() + { + // This function is only valid for non-preconditioners. + if ( _is_preconditioner ) + throw std::logic_error( "Cannot call setup() on preconditioners" ); + + auto error = HYPRE_SStructMatrixAssemble( _A ); + checkHypreError( error ); + + this->setupImpl( _A, _b, _x ); + } + + /*! + \brief Solve the problem Ax = b for x. + \param b The forcing term. + \param x The solution. + \param n_vars Number of variables in the solution domain + */ + template + void solve( const Array_t& b, Array_t& x, int n_vars = 3 ) + { + Kokkos::Profiling::pushRegion( + "Cajita::HypreSemiStructuredSolver::solve" ); + + static_assert( is_array::value, "Must use an array" ); + static_assert( + std::is_same::value, + "Array entity type mush match solver entity type" ); + static_assert( + std::is_same::value, + "Array device type and solver device type are different." ); + + static_assert( + std::is_same::value, + "Array value type and solver value type are different." ); + + // This function is only valid for non-preconditioners. + if ( _is_preconditioner ) + throw std::logic_error( "Cannot call solve() on preconditioners" ); + + // Spatial dimension. + const std::size_t num_space_dim = Array_t::num_space_dim; + + int part = 0; + + // Initialize the RHS. + auto error = HYPRE_SStructVectorInitialize( _b ); + checkHypreError( error ); + + // Copy the RHS into HYPRE. The HYPRE layout is fixed as layout-right. + auto owned_space = b.layout()->indexSpace( Own(), Local() ); + std::array reorder_min; + std::array reorder_max; + for ( std::size_t d = 0; d < num_space_dim; ++d ) + { + reorder_min[d] = owned_space.min( d ); + reorder_max[d] = owned_space.max( d ); + } + + // Insert b values into the HYPRE vector. + // The process of creating the view and then deep copying each + // variable is functional, but we should avoid this process + // for performance if possible + for ( int var = 0; var < n_vars; ++var ) + { + reorder_min.back() = var; + reorder_max.back() = var + 1; + + IndexSpace reorder_space( reorder_min, + reorder_max ); + auto b_values = + createView( + "vector_values", reorder_space ); + // Extract one variable at at time. + auto b_subv = createSubview( b.view(), reorder_space ); + + Kokkos::deep_copy( b_values, b_subv ); + error = HYPRE_SStructVectorSetBoxValues( + _b, part, _lower.data(), _upper.data(), var, b_values.data() ); + checkHypreError( error ); + } + + error = HYPRE_SStructVectorAssemble( _b ); + checkHypreError( error ); + + // Solve the problem + this->solveImpl( _A, _b, _x ); + + // Extract the solution from the LHS + for ( int var = 0; var < n_vars; ++var ) + { + reorder_min.back() = var; + reorder_max.back() = var + 1; + + IndexSpace reorder_space( reorder_min, + reorder_max ); + auto x_values = + createView( + "vector_values", reorder_space ); + + // Extract one variable at at time. + // Use a pair here to retain the view rank. + + error = HYPRE_SStructVectorGetBoxValues( + _x, part, _lower.data(), _upper.data(), var, x_values.data() ); + checkHypreError( error ); + + // Copy the HYPRE solution to the LHS. + auto x_subv = createSubview( x.view(), reorder_space ); + Kokkos::deep_copy( x_subv, x_values ); + } + + Kokkos::Profiling::popRegion(); + } + + //! Get the number of iterations taken on the last solve. + int getNumIter() { return this->getNumIterImpl(); } + + //! Get the relative residual norm achieved on the last solve. + double getFinalRelativeResidualNorm() + { + return this->getFinalRelativeResidualNormImpl(); + } + + //! Get the preconditioner. + virtual HYPRE_SStructSolver getHypreSolver() const = 0; + //! Get the preconditioner setup function. + virtual HYPRE_PtrToSStructSolverFcn getHypreSetupFunction() const = 0; + //! Get the preconditioner solve function. + virtual HYPRE_PtrToSStructSolverFcn getHypreSolveFunction() const = 0; + + protected: + //! Set convergence tolerance implementation. + virtual void setToleranceImpl( const double tol ) = 0; + + //! Set maximum iteration implementation. + virtual void setMaxIterImpl( const int max_iter ) = 0; + + //! Set the output level. + virtual void setPrintLevelImpl( const int print_level ) = 0; + + //! Setup implementation. + virtual void setupImpl( HYPRE_SStructMatrix A, HYPRE_SStructVector b, + HYPRE_SStructVector x ) = 0; + + //! Solver implementation. + virtual void solveImpl( HYPRE_SStructMatrix A, HYPRE_SStructVector b, + HYPRE_SStructVector x ) = 0; + + //! Get the number of iterations taken on the last solve. + virtual int getNumIterImpl() = 0; + + //! Get the relative residual norm achieved on the last solve. + virtual double getFinalRelativeResidualNormImpl() = 0; + + //! Set a preconditioner. + virtual void setPreconditionerImpl( + const HypreSemiStructuredSolver& + preconditioner ) = 0; + + //! Check a hypre error. + void checkHypreError( const int error ) const + { + if ( error > 0 ) + { + char error_msg[256]; + HYPRE_DescribeError( error, error_msg ); + std::stringstream out; + out << "HYPRE semi-structured solver error: "; + out << error << " " << error_msg; + HYPRE_ClearError( error ); + throw std::runtime_error( out.str() ); + } + } + + private: + MPI_Comm _comm; + bool _is_preconditioner; + HYPRE_SStructGrid _grid; + std::vector _lower; + std::vector _upper; + std::vector _stencils; + HYPRE_SStructGraph _graph; + std::vector _stencil_size; + std::vector> _stencil_index; + HYPRE_SStructMatrix _A; + HYPRE_SStructVector _b; + HYPRE_SStructVector _x; + std::shared_ptr> + _preconditioner; +}; + +//---------------------------------------------------------------------------// +//! PCG solver. +template +class HypreSemiStructPCG + : public HypreSemiStructuredSolver +{ + public: + //! Base HYPRE semi-structured solver type. + using Base = HypreSemiStructuredSolver; + //! Constructor + template + HypreSemiStructPCG( const ArrayLayout_t& layout, int n_vars, + const bool is_preconditioner = false ) + : Base( layout, n_vars, is_preconditioner ) + { + if ( is_preconditioner ) + throw std::logic_error( + "HYPRE PCG cannot be used as a preconditioner" ); + + auto error = HYPRE_SStructPCGCreate( + layout.localGrid()->globalGrid().comm(), &_solver ); + this->checkHypreError( error ); + + HYPRE_SStructPCGSetTwoNorm( _solver, 1 ); + } + + ~HypreSemiStructPCG() { HYPRE_SStructPCGDestroy( _solver ); } + + // PCG SETTINGS + + //! Set the absolute tolerance + void setAbsoluteTol( const double tol ) + { + auto error = HYPRE_SStructPCGSetAbsoluteTol( _solver, tol ); + this->checkHypreError( error ); + } + + //! Additionally require that the relative difference in successive + //! iterates be small. + void setRelChange( const int rel_change ) + { + auto error = HYPRE_SStructPCGSetRelChange( _solver, rel_change ); + this->checkHypreError( error ); + } + + //! Set the amount of logging to do. + void setLogging( const int logging ) + { + auto error = HYPRE_SStructPCGSetLogging( _solver, logging ); + this->checkHypreError( error ); + } + + HYPRE_SStructSolver getHypreSolver() const override { return _solver; } + HYPRE_PtrToSStructSolverFcn getHypreSetupFunction() const override + { + return HYPRE_SStructPCGSetup; + } + HYPRE_PtrToSStructSolverFcn getHypreSolveFunction() const override + { + return HYPRE_SStructPCGSolve; + } + + protected: + void setToleranceImpl( const double tol ) override + { + auto error = HYPRE_SStructPCGSetTol( _solver, tol ); + this->checkHypreError( error ); + } + + void setMaxIterImpl( const int max_iter ) override + { + auto error = HYPRE_SStructPCGSetMaxIter( _solver, max_iter ); + this->checkHypreError( error ); + } + + void setPrintLevelImpl( const int print_level ) override + { + auto error = HYPRE_SStructPCGSetPrintLevel( _solver, print_level ); + this->checkHypreError( error ); + } + + void setupImpl( HYPRE_SStructMatrix A, HYPRE_SStructVector b, + HYPRE_SStructVector x ) override + { + auto error = HYPRE_SStructPCGSetup( _solver, A, b, x ); + this->checkHypreError( error ); + } + + void solveImpl( HYPRE_SStructMatrix A, HYPRE_SStructVector b, + HYPRE_SStructVector x ) override + { + auto error = HYPRE_SStructPCGSolve( _solver, A, b, x ); + this->checkHypreError( error ); + } + + int getNumIterImpl() override + { + HYPRE_Int num_iter; + auto error = HYPRE_SStructPCGGetNumIterations( _solver, &num_iter ); + this->checkHypreError( error ); + return num_iter; + } + + double getFinalRelativeResidualNormImpl() override + { + HYPRE_Real norm; + auto error = + HYPRE_SStructPCGGetFinalRelativeResidualNorm( _solver, &norm ); + this->checkHypreError( error ); + return norm; + } + + void setPreconditionerImpl( + const HypreSemiStructuredSolver& + preconditioner ) override + { + auto error = HYPRE_SStructPCGSetPrecond( + _solver, preconditioner.getHypreSolveFunction(), + preconditioner.getHypreSetupFunction(), + preconditioner.getHypreSolver() ); + this->checkHypreError( error ); + } + + private: + HYPRE_SStructSolver _solver; +}; + +//---------------------------------------------------------------------------// +//! GMRES solver. +template +class HypreSemiStructGMRES + : public HypreSemiStructuredSolver +{ + public: + //! Base HYPRE semi-structured solver type. + using Base = HypreSemiStructuredSolver; + //! Constructor + template + HypreSemiStructGMRES( const ArrayLayout_t& layout, int n_vars, + const bool is_preconditioner = false ) + : Base( layout, n_vars, is_preconditioner ) + { + if ( is_preconditioner ) + throw std::logic_error( + "HYPRE GMRES cannot be used as a preconditioner" ); + + auto error = HYPRE_SStructGMRESCreate( + layout.localGrid()->globalGrid().comm(), &_solver ); + this->checkHypreError( error ); + } + + ~HypreSemiStructGMRES() { HYPRE_SStructGMRESDestroy( _solver ); } + + // GMRES SETTINGS + + //! Set the absolute tolerance + void setAbsoluteTol( const double tol ) + { + auto error = HYPRE_SStructGMRESSetAbsoluteTol( _solver, tol ); + this->checkHypreError( error ); + } + + //! Set the max size of the Krylov space. + void setKDim( const int k_dim ) + { + auto error = HYPRE_SStructGMRESSetKDim( _solver, k_dim ); + this->checkHypreError( error ); + } + + //! Set the amount of logging to do. + void setLogging( const int logging ) + { + auto error = HYPRE_SStructGMRESSetLogging( _solver, logging ); + this->checkHypreError( error ); + } + + HYPRE_SStructSolver getHypreSolver() const override { return _solver; } + HYPRE_PtrToSStructSolverFcn getHypreSetupFunction() const override + { + return HYPRE_SStructGMRESSetup; + } + HYPRE_PtrToSStructSolverFcn getHypreSolveFunction() const override + { + return HYPRE_SStructGMRESSolve; + } + + protected: + void setToleranceImpl( const double tol ) override + { + auto error = HYPRE_SStructGMRESSetTol( _solver, tol ); + this->checkHypreError( error ); + } + + void setMaxIterImpl( const int max_iter ) override + { + auto error = HYPRE_SStructGMRESSetMaxIter( _solver, max_iter ); + this->checkHypreError( error ); + } + + void setPrintLevelImpl( const int print_level ) override + { + auto error = HYPRE_SStructGMRESSetPrintLevel( _solver, print_level ); + this->checkHypreError( error ); + } + + void setupImpl( HYPRE_SStructMatrix A, HYPRE_SStructVector b, + HYPRE_SStructVector x ) override + { + auto error = HYPRE_SStructGMRESSetup( _solver, A, b, x ); + this->checkHypreError( error ); + } + + void solveImpl( HYPRE_SStructMatrix A, HYPRE_SStructVector b, + HYPRE_SStructVector x ) override + { + auto error = HYPRE_SStructGMRESSolve( _solver, A, b, x ); + this->checkHypreError( error ); + } + + int getNumIterImpl() override + { + HYPRE_Int num_iter; + auto error = HYPRE_SStructGMRESGetNumIterations( _solver, &num_iter ); + this->checkHypreError( error ); + return num_iter; + } + + double getFinalRelativeResidualNormImpl() override + { + HYPRE_Real norm; + auto error = + HYPRE_SStructGMRESGetFinalRelativeResidualNorm( _solver, &norm ); + this->checkHypreError( error ); + return norm; + } + + void setPreconditionerImpl( + const HypreSemiStructuredSolver& + preconditioner ) override + { + auto error = HYPRE_SStructGMRESSetPrecond( + _solver, preconditioner.getHypreSolveFunction(), + preconditioner.getHypreSetupFunction(), + preconditioner.getHypreSolver() ); + this->checkHypreError( error ); + } + + private: + HYPRE_SStructSolver _solver; +}; + +//---------------------------------------------------------------------------// +//! BiCGSTAB solver. +template +class HypreSemiStructBiCGSTAB + : public HypreSemiStructuredSolver +{ + public: + //! Base HYPRE semi-structured solver type. + using Base = HypreSemiStructuredSolver; + //! Constructor + template + HypreSemiStructBiCGSTAB( const ArrayLayout_t& layout, + const bool is_preconditioner = false, + int n_vars = 3 ) + : Base( layout, n_vars, is_preconditioner ) + { + if ( is_preconditioner ) + throw std::logic_error( + "HYPRE BiCGSTAB cannot be used as a preconditioner" ); + + auto error = HYPRE_SStructBiCGSTABCreate( + layout.localGrid()->globalGrid().comm(), &_solver ); + this->checkHypreError( error ); + } + + ~HypreSemiStructBiCGSTAB() { HYPRE_SStructBiCGSTABDestroy( _solver ); } + + // BiCGSTAB SETTINGS + + //! Set the absolute tolerance + void setAbsoluteTol( const double tol ) + { + auto error = HYPRE_SStructBiCGSTABSetAbsoluteTol( _solver, tol ); + this->checkHypreError( error ); + } + + //! Set the amount of logging to do. + void setLogging( const int logging ) + { + auto error = HYPRE_SStructBiCGSTABSetLogging( _solver, logging ); + this->checkHypreError( error ); + } + + HYPRE_SStructSolver getHypreSolver() const override { return _solver; } + HYPRE_PtrToSStructSolverFcn getHypreSetupFunction() const override + { + return HYPRE_SStructBiCGSTABSetup; + } + HYPRE_PtrToSStructSolverFcn getHypreSolveFunction() const override + { + return HYPRE_SStructBiCGSTABSolve; + } + + protected: + void setToleranceImpl( const double tol ) override + { + auto error = HYPRE_SStructBiCGSTABSetTol( _solver, tol ); + this->checkHypreError( error ); + } + + void setMaxIterImpl( const int max_iter ) override + { + auto error = HYPRE_SStructBiCGSTABSetMaxIter( _solver, max_iter ); + this->checkHypreError( error ); + } + + void setPrintLevelImpl( const int print_level ) override + { + auto error = HYPRE_SStructBiCGSTABSetPrintLevel( _solver, print_level ); + this->checkHypreError( error ); + } + + void setupImpl( HYPRE_SStructMatrix A, HYPRE_SStructVector b, + HYPRE_SStructVector x ) override + { + auto error = HYPRE_SStructBiCGSTABSetup( _solver, A, b, x ); + this->checkHypreError( error ); + } + + void solveImpl( HYPRE_SStructMatrix A, HYPRE_SStructVector b, + HYPRE_SStructVector x ) override + { + auto error = HYPRE_SStructBiCGSTABSolve( _solver, A, b, x ); + this->checkHypreError( error ); + } + + int getNumIterImpl() override + { + HYPRE_Int num_iter; + auto error = + HYPRE_SStructBiCGSTABGetNumIterations( _solver, &num_iter ); + this->checkHypreError( error ); + return num_iter; + } + + double getFinalRelativeResidualNormImpl() override + { + HYPRE_Real norm; + auto error = + HYPRE_SStructBiCGSTABGetFinalRelativeResidualNorm( _solver, &norm ); + this->checkHypreError( error ); + return norm; + } + + void setPreconditionerImpl( + const HypreSemiStructuredSolver& + preconditioner ) override + { + auto error = HYPRE_SStructBiCGSTABSetPrecond( + _solver, preconditioner.getHypreSolveFunction(), + preconditioner.getHypreSetupFunction(), + preconditioner.getHypreSolver() ); + this->checkHypreError( error ); + } + + private: + HYPRE_SStructSolver _solver; +}; + +//---------------------------------------------------------------------------// +//! Diagonal preconditioner. +template +class HypreSemiStructDiagonal + : public HypreSemiStructuredSolver +{ + public: + //! Base HYPRE semi-structured solver type. + using Base = HypreSemiStructuredSolver; + //! Constructor + template + HypreSemiStructDiagonal( const ArrayLayout_t& layout, + const bool is_preconditioner = false, + int n_vars = 3 ) + : Base( layout, n_vars, is_preconditioner ) + { + if ( !is_preconditioner ) + throw std::logic_error( + "Diagonal preconditioner cannot be used as a solver" ); + } + + HYPRE_SStructSolver getHypreSolver() const override { return nullptr; } + HYPRE_PtrToSStructSolverFcn getHypreSetupFunction() const override + { + return HYPRE_SStructDiagScaleSetup; + } + HYPRE_PtrToSStructSolverFcn getHypreSolveFunction() const override + { + return HYPRE_SStructDiagScale; + } + + protected: + void setToleranceImpl( const double ) override + { + throw std::logic_error( + "Diagonal preconditioner cannot be used as a solver" ); + } + + void setMaxIterImpl( const int ) override + { + throw std::logic_error( + "Diagonal preconditioner cannot be used as a solver" ); + } + + void setPrintLevelImpl( const int ) override + { + throw std::logic_error( + "Diagonal preconditioner cannot be used as a solver" ); + } + + void setupImpl( HYPRE_SStructMatrix, HYPRE_SStructVector, + HYPRE_SStructVector ) override + { + throw std::logic_error( + "Diagonal preconditioner cannot be used as a solver" ); + } + + void solveImpl( HYPRE_SStructMatrix, HYPRE_SStructVector, + HYPRE_SStructVector ) override + { + throw std::logic_error( + "Diagonal preconditioner cannot be used as a solver" ); + } + + int getNumIterImpl() override + { + throw std::logic_error( + "Diagonal preconditioner cannot be used as a solver" ); + } + + double getFinalRelativeResidualNormImpl() override + { + throw std::logic_error( + "Diagonal preconditioner cannot be used as a solver" ); + } + + void setPreconditionerImpl( + const HypreSemiStructuredSolver& ) + override + { + throw std::logic_error( + "Diagonal preconditioner does not support preconditioning." ); + } +}; + +//---------------------------------------------------------------------------// +// Builders +//---------------------------------------------------------------------------// +//! Create a HYPRE PCG semi-structured solver. +template +std::shared_ptr> +createHypreSemiStructPCG( const ArrayLayout_t& layout, + const bool is_preconditioner = false, int n_vars = 3 ) +{ + static_assert( is_array_layout::value, + "Must use an array layout" ); + return std::make_shared>( + layout, n_vars, is_preconditioner ); +} + +//! Create a HYPRE GMRES semi-structured solver. +template +std::shared_ptr> +createHypreSemiStructGMRES( const ArrayLayout_t& layout, + const bool is_preconditioner = false, + int n_vars = 3 ) +{ + static_assert( is_array_layout::value, + "Must use an array layout" ); + return std::make_shared>( + layout, n_vars, is_preconditioner ); +} + +//! Create a HYPRE BiCGSTAB semi-structured solver. +template +std::shared_ptr> +createHypreSemiStructBiCGSTAB( const ArrayLayout_t& layout, + const bool is_preconditioner = false, + int n_vars = 3 ) +{ + static_assert( is_array_layout::value, + "Must use an array layout" ); + return std::make_shared>( + layout, is_preconditioner, n_vars ); +} + +//! Create a HYPRE Diagonal semi-structured solver. +template +std::shared_ptr> +createHypreSemiStructDiagonal( const ArrayLayout_t& layout, + const bool is_preconditioner = false, + int n_vars = 3 ) +{ + static_assert( is_array_layout::value, + "Must use an array layout" ); + return std::make_shared>( + layout, is_preconditioner, n_vars ); +} + +//---------------------------------------------------------------------------// +// Factory +//---------------------------------------------------------------------------// +/*! + \brief Create a HYPRE semi-structured solver. + + \param solver_type Solver name. + \param layout The ArrayLayout defining the vector space of the solver. + \param is_preconditioner Use as a preconditioner. + \param n_vars Number of variables in the solver +*/ +template +std::shared_ptr> +createHypreSemiStructuredSolver( const std::string& solver_type, + const ArrayLayout_t& layout, + const bool is_preconditioner = false, + int n_vars = 3 ) +{ + static_assert( is_array_layout::value, + "Must use an array layout" ); + + if ( "PCG" == solver_type ) + return createHypreSemiStructPCG( + layout, is_preconditioner, n_vars ); + else if ( "GMRES" == solver_type ) + return createHypreSemiStructGMRES( + layout, is_preconditioner, n_vars ); + else if ( "BiCGSTAB" == solver_type ) + return createHypreSemiStructBiCGSTAB( + layout, is_preconditioner, n_vars ); + else if ( "Diagonal" == solver_type ) + return createHypreSemiStructDiagonal( + layout, is_preconditioner, n_vars ); + else + throw std::runtime_error( "Invalid solver type" ); +} + +//---------------------------------------------------------------------------// + +} // end namespace Cajita + +#endif // end CAJITA_HypreSemiStRUCTUREDSOLVER_HPP diff --git a/cajita/src/Cajita_HypreStructuredSolver.hpp b/cajita/src/Cajita_HypreStructuredSolver.hpp index 269306546..c7e7a4b60 100644 --- a/cajita/src/Cajita_HypreStructuredSolver.hpp +++ b/cajita/src/Cajita_HypreStructuredSolver.hpp @@ -18,6 +18,7 @@ #include #include +#include #include #include #include @@ -38,60 +39,6 @@ namespace Cajita { -//---------------------------------------------------------------------------// -// Hypre memory space selection. Don't compile if HYPRE wasn't configured to -// use the device. -// ---------------------------------------------------------------------------// - -//! Hypre device compatibility check. -template -struct HypreIsCompatibleWithMemorySpace : std::false_type -{ -}; - -// FIXME: This is currently written in this structure because HYPRE only has -// compile-time switches for backends and hence only one can be used at a -// time. Once they have a run-time switch we can use that instead. -#ifdef HYPRE_USING_CUDA -#ifdef KOKKOS_ENABLE_CUDA -#ifdef HYPRE_USING_DEVICE_MEMORY -//! Hypre device compatibility check - CUDA memory. -template <> -struct HypreIsCompatibleWithMemorySpace : std::true_type -{ -}; -#endif // end HYPRE_USING_DEVICE_MEMORY - -//! Hypre device compatibility check - CUDA UVM memory. -#ifdef HYPRE_USING_UNIFIED_MEMORY -template <> -struct HypreIsCompatibleWithMemorySpace : std::true_type -{ -}; -#endif // end HYPRE_USING_UNIFIED_MEMORY -#endif // end KOKKOS_ENABLE_CUDA -#endif // end HYPRE_USING_CUDA - -#ifdef HYPRE_USING_HIP -#ifdef KOKKOS_ENABLE_HIP -//! Hypre device compatibility check - HIP memory. FIXME - make this true when -//! the HYPRE CMake includes HIP -template <> -struct HypreIsCompatibleWithMemorySpace - : std::false_type -{ -}; -#endif // end KOKKOS_ENABLE_HIP -#endif // end HYPRE_USING_HIP - -#ifndef HYPRE_USING_GPU -//! Hypre device compatibility check - host memory. -template <> -struct HypreIsCompatibleWithMemorySpace : std::true_type -{ -}; -#endif // end HYPRE_USING_GPU - //---------------------------------------------------------------------------// //! Hypre structured solver interface for scalar fields. template @@ -333,6 +280,33 @@ class HypreStructuredSolver checkHypreError( error ); } + /*! + \brief Print the hypre matrix to ouput file + \param prefix File prefix for where hypre output is written + */ + void printMatrix( const char* prefix ) + { + HYPRE_StructMatrixPrint( prefix, _A, 0 ); + } + + /*! + \brief Print the hypre LHS to ouput file + \param prefix File prefix for where hypre output is written + */ + void printLHS( const char* prefix ) + { + HYPRE_StructVectorPrint( prefix, _x, 0 ); + } + + /*! + \brief Print the hypre RHS to ouput file + \param prefix File prefix for where hypre output is written + */ + void printRHS( const char* prefix ) + { + HYPRE_StructVectorPrint( prefix, _b, 0 ); + } + //! Set convergence tolerance implementation. void setTolerance( const double tol ) { this->setToleranceImpl( tol ); } diff --git a/cajita/unit_test/CMakeLists.txt b/cajita/unit_test/CMakeLists.txt index f0e1f4d27..750187e5b 100644 --- a/cajita/unit_test/CMakeLists.txt +++ b/cajita/unit_test/CMakeLists.txt @@ -56,6 +56,8 @@ if(Cabana_ENABLE_HYPRE) list(APPEND MPI_TESTS HypreStructuredSolver3d HypreStructuredSolver2d + HypreSemiStructuredSolver + HypreSemiStructuredSolverMulti ) endif() diff --git a/cajita/unit_test/tstHypreSemiStructuredSolver.hpp b/cajita/unit_test/tstHypreSemiStructuredSolver.hpp new file mode 100644 index 000000000..085cdecac --- /dev/null +++ b/cajita/unit_test/tstHypreSemiStructuredSolver.hpp @@ -0,0 +1,260 @@ +/**************************************************************************** + * Copyright (c) 2018-2022 by the Cabana authors * + * All rights reserved. * + * * + * This file is part of the Cabana library. Cabana is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include +#include + +using namespace Cajita; + +namespace Test +{ + +//---------------------------------------------------------------------------// +// FIXME: Only run test if HYPRE is compatible with the memory space. This +// is currently written in this structure because HYPRE only has +// compile-time switches for backends and hence only one can be used at a +// time. Once they have a run-time switch we can use that instead. +template +std::enable_if_t::value, void> +poissonTest( const std::string&, const std::string&, MemorySpace ) +{ +} + +template +std::enable_if_t::value, void> +poissonTest( const std::string& solver_type, const std::string& precond_type, + MemorySpace ) +{ + // Create the global grid. + double cell_size = 0.25; + std::array is_dim_periodic = { false, false, false }; + std::array global_low_corner = { -1.0, -2.0, -1.0 }; + std::array global_high_corner = { 1.0, 1.0, 0.5 }; + auto global_mesh = createUniformGlobalMesh( global_low_corner, + global_high_corner, cell_size ); + + // Create the global grid. + DimBlockPartitioner<3> partitioner; + auto global_grid = createGlobalGrid( MPI_COMM_WORLD, global_mesh, + is_dim_periodic, partitioner ); + + // Create a local grid. + auto local_mesh = createLocalGrid( global_grid, 1 ); + auto owned_space = local_mesh->indexSpace( Own(), Cell(), Local() ); + + // Create the RHS. + auto vector_layout = createArrayLayout( local_mesh, 1, Cell() ); + auto rhs = createArray( "rhs", vector_layout ); + ArrayOp::assign( *rhs, 1.0, Own() ); + + // Create the LHS. + auto lhs = createArray( "lhs", vector_layout ); + ArrayOp::assign( *lhs, 0.0, Own() ); + + HYPRE_Init(); + + // Create a solver. + auto solver = createHypreSemiStructuredSolver( + solver_type, *vector_layout, false, 1 ); + + // Create a 7-point 3d laplacian stencil. + std::vector> stencil = { + { 0, 0, 0 }, { -1, 0, 0 }, { 1, 0, 0 }, { 0, -1, 0 }, + { 0, 1, 0 }, { 0, 0, -1 }, { 0, 0, 1 } }; + solver->createMatrixStencil( 3, 0, 1, { 7 } ); + solver->setMatrixStencil( stencil, 0, 0 ); + + solver->setSolverGraph( 1 ); + + // Create the matrix entries. The stencil is defined over cells. + auto matrix_entry_layout = createArrayLayout( local_mesh, 7, Cell() ); + auto matrix_entries = createArray( + "matrix_entries", matrix_entry_layout ); + auto entry_view = matrix_entries->view(); + Kokkos::parallel_for( + "fill_matrix_entries", + createExecutionPolicy( owned_space, TEST_EXECSPACE() ), + KOKKOS_LAMBDA( const int i, const int j, const int k ) { + entry_view( i, j, k, 0 ) = 6.0; + entry_view( i, j, k, 1 ) = -1.0; + entry_view( i, j, k, 2 ) = -1.0; + entry_view( i, j, k, 3 ) = -1.0; + entry_view( i, j, k, 4 ) = -1.0; + entry_view( i, j, k, 5 ) = -1.0; + entry_view( i, j, k, 6 ) = -1.0; + } ); + + solver->initializeHypreMatrix(); + + solver->setMatrixValues( *matrix_entries, 0, 0 ); + + // Set the tolerance. + solver->setTolerance( 1.0e-9 ); + + // Set the maximum iterations. + solver->setMaxIter( 2000 ); + + // Set the print level. + solver->setPrintLevel( 2 ); + + // Create a preconditioner. + if ( "none" != precond_type ) + { + auto preconditioner = + createHypreSemiStructuredSolver( + precond_type, *vector_layout, true ); + solver->setPreconditioner( preconditioner ); + } + + // Setup the problem. + solver->setup(); + + // Solve the problem. + solver->solve( *rhs, *lhs, 1 ); + + // Create a solver reference for comparison. + auto lhs_ref = createArray( "lhs_ref", vector_layout ); + ArrayOp::assign( *lhs_ref, 0.0, Own() ); + + auto ref_solver = + createReferenceConjugateGradient( *vector_layout ); + ref_solver->setMatrixStencil( stencil ); + const auto& ref_entries = ref_solver->getMatrixValues(); + auto matrix_view = ref_entries.view(); + auto global_space = local_mesh->indexSpace( Own(), Cell(), Global() ); + int ncell_i = global_grid->globalNumEntity( Cell(), Dim::I ); + int ncell_j = global_grid->globalNumEntity( Cell(), Dim::J ); + int ncell_k = global_grid->globalNumEntity( Cell(), Dim::K ); + Kokkos::parallel_for( + "fill_ref_entries", + createExecutionPolicy( owned_space, TEST_EXECSPACE() ), + KOKKOS_LAMBDA( const int i, const int j, const int k ) { + int gi = i + global_space.min( Dim::I ) - owned_space.min( Dim::I ); + int gj = j + global_space.min( Dim::J ) - owned_space.min( Dim::J ); + int gk = k + global_space.min( Dim::K ) - owned_space.min( Dim::K ); + matrix_view( i, j, k, 0 ) = 6.0; + matrix_view( i, j, k, 1 ) = ( gi - 1 >= 0 ) ? -1.0 : 0.0; + matrix_view( i, j, k, 2 ) = ( gi + 1 < ncell_i ) ? -1.0 : 0.0; + matrix_view( i, j, k, 3 ) = ( gj - 1 >= 0 ) ? -1.0 : 0.0; + matrix_view( i, j, k, 4 ) = ( gj + 1 < ncell_j ) ? -1.0 : 0.0; + matrix_view( i, j, k, 5 ) = ( gk - 1 >= 0 ) ? -1.0 : 0.0; + matrix_view( i, j, k, 6 ) = ( gk + 1 < ncell_k ) ? -1.0 : 0.0; + } ); + + std::vector> diag_stencil = { { 0, 0, 0 } }; + ref_solver->setPreconditionerStencil( diag_stencil ); + const auto& preconditioner_entries = ref_solver->getPreconditionerValues(); + auto preconditioner_view = preconditioner_entries.view(); + Kokkos::parallel_for( + "fill_preconditioner_entries", + createExecutionPolicy( owned_space, TEST_EXECSPACE() ), + KOKKOS_LAMBDA( const int i, const int j, const int k ) { + preconditioner_view( i, j, k, 0 ) = 1.0 / 6.0; + } ); + + ref_solver->setTolerance( 1.0e-11 ); + ref_solver->setPrintLevel( 1 ); + ref_solver->setup(); + ref_solver->solve( *rhs, *lhs_ref ); + + // Check the results. + auto lhs_host = + Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace(), lhs->view() ); + auto lhs_ref_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), lhs_ref->view() ); + for ( int i = owned_space.min( Dim::I ); i < owned_space.max( Dim::I ); + ++i ) + for ( int j = owned_space.min( Dim::J ); j < owned_space.max( Dim::J ); + ++j ) + for ( int k = owned_space.min( Dim::K ); + k < owned_space.max( Dim::K ); ++k ) + EXPECT_FLOAT_EQ( lhs_host( i, j, k, 0 ), + lhs_ref_host( i, j, k, 0 ) ); + + // Setup the problem again. We would need to do this if we changed the + // matrix entries. + solver->setup(); + + // Solve the problem again + ArrayOp::assign( *rhs, 2.0, Own() ); + ArrayOp::assign( *lhs, 0.0, Own() ); + solver->solve( *rhs, *lhs, 1 ); + + // Compute another reference solution. + ArrayOp::assign( *lhs_ref, 0.0, Own() ); + ref_solver->solve( *rhs, *lhs_ref ); + + // Check the results again + lhs_host = + Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace(), lhs->view() ); + lhs_ref_host = Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace(), + lhs_ref->view() ); + for ( int i = owned_space.min( Dim::I ); i < owned_space.max( Dim::I ); + ++i ) + for ( int j = owned_space.min( Dim::J ); j < owned_space.max( Dim::J ); + ++j ) + for ( int k = owned_space.min( Dim::K ); + k < owned_space.max( Dim::K ); ++k ) + EXPECT_FLOAT_EQ( lhs_host( i, j, k, 0 ), + lhs_ref_host( i, j, k, 0 ) ); + + HYPRE_Finalize(); +} + +//---------------------------------------------------------------------------// +// RUN TESTS +//---------------------------------------------------------------------------// +TEST( semi_structured_solver, pcg_none_test ) +{ + poissonTest( "PCG", "none", TEST_MEMSPACE{} ); +} + +TEST( semi_structured_solver, gmres_none_test ) +{ + poissonTest( "GMRES", "none", TEST_MEMSPACE{} ); +} + +TEST( semi_structured_solver, bicgstab_none_test ) +{ + poissonTest( "BiCGSTAB", "none", TEST_MEMSPACE{} ); +} + +TEST( semi_structured_solver, pcg_diag_test ) +{ + poissonTest( "PCG", "Diagonal", TEST_MEMSPACE{} ); +} + +TEST( semi_structured_solver, gmres_diag_test ) +{ + poissonTest( "GMRES", "Diagonal", TEST_MEMSPACE{} ); +} + +TEST( semi_structured_solver, bicgstab_diag_test ) +{ + poissonTest( "BiCGSTAB", "Diagonal", TEST_MEMSPACE{} ); +} + +//---------------------------------------------------------------------------// + +} // end namespace Test diff --git a/cajita/unit_test/tstHypreSemiStructuredSolverMulti.hpp b/cajita/unit_test/tstHypreSemiStructuredSolverMulti.hpp new file mode 100644 index 000000000..455fde2dd --- /dev/null +++ b/cajita/unit_test/tstHypreSemiStructuredSolverMulti.hpp @@ -0,0 +1,282 @@ + + +/**************************************************************************** + * Copyright (c) 2018-2022 by the Cabana authors * + * All rights reserved. * + * * + * This file is part of the Cabana library. Cabana is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include +#include + +using namespace Cajita; + +namespace Test +{ + +//---------------------------------------------------------------------------// +// FIXME: Only run test if HYPRE is compatible with the memory space. This +// is currently written in this structure because HYPRE only has +// compile-time switches for backends and hence only one can be used at a +// time. Once they have a run-time switch we can use that instead. +template +std::enable_if_t::value, void> +poissonTest( const std::string&, const std::string&, MemorySpace ) +{ +} + +template +std::enable_if_t::value, void> +poissonTest( const std::string& solver_type, const std::string& precond_type, + MemorySpace ) +{ + // Create the global grid. + double cell_size = 0.25; + std::array is_dim_periodic = { false, false, false }; + std::array global_low_corner = { -1.0, -2.0, -1.0 }; + std::array global_high_corner = { 1.0, 1.0, 0.5 }; + auto global_mesh = createUniformGlobalMesh( global_low_corner, + global_high_corner, cell_size ); + + // Create the global grid. + DimBlockPartitioner<3> partitioner; + auto global_grid = createGlobalGrid( MPI_COMM_WORLD, global_mesh, + is_dim_periodic, partitioner ); + + // Create a local grid. + auto local_mesh = createLocalGrid( global_grid, 1 ); + auto owned_space = local_mesh->indexSpace( Own(), Cell(), Local() ); + + // Create the RHS. + auto vector_layout = createArrayLayout( local_mesh, 3, Cell() ); + auto rhs = createArray( "rhs", vector_layout ); + // auto rhs = createArray( + // "rhs", vector_layout ); + ArrayOp::assign( *rhs, 1.0, Own() ); + + // Create the LHS. + auto lhs = createArray( "lhs", vector_layout ); + // auto lhs = createArray( + // "lhs", vector_layout ); + ArrayOp::assign( *lhs, 0.0, Own() ); + + HYPRE_Init(); + + // Create a solver. + auto solver = createHypreSemiStructuredSolver( + solver_type, *vector_layout, false, 3 ); + + // Create a 7-point 3d laplacian stencil. + std::vector> stencil = { + { 0, 0, 0 }, { -1, 0, 0 }, { 1, 0, 0 }, { 0, -1, 0 }, + { 0, 1, 0 }, { 0, 0, -1 }, { 0, 0, 1 } }; + + solver->createMatrixStencil( 3, 0, 3, { 7, 0, 0 } ); + solver->createMatrixStencil( 3, 1, 3, { 0, 7, 0 } ); + solver->createMatrixStencil( 3, 2, 3, { 0, 0, 7 } ); + for ( int v_h = 0; v_h < 3; ++v_h ) + { + solver->setMatrixStencil( stencil, v_h, v_h ); + } + + solver->setSolverGraph( 3 ); + + // Create the matrix entries. The stencil is defined over cells. + auto matrix_entry_layout = createArrayLayout( local_mesh, 7, Cell() ); + auto matrix_entries = createArray( + "matrix_entries", matrix_entry_layout ); + auto entry_view = matrix_entries->view(); + Kokkos::parallel_for( + "fill_matrix_entries", + createExecutionPolicy( owned_space, TEST_EXECSPACE() ), + KOKKOS_LAMBDA( const int i, const int j, const int k ) { + entry_view( i, j, k, 0 ) = 6.0; + entry_view( i, j, k, 1 ) = -1.0; + entry_view( i, j, k, 2 ) = -1.0; + entry_view( i, j, k, 3 ) = -1.0; + entry_view( i, j, k, 4 ) = -1.0; + entry_view( i, j, k, 5 ) = -1.0; + entry_view( i, j, k, 6 ) = -1.0; + } ); + + solver->initializeHypreMatrix(); + + for ( int v_h = 0; v_h < 3; ++v_h ) + { + solver->setMatrixValues( *matrix_entries, v_h, v_h ); + } + + // Set the tolerance. + solver->setTolerance( 1.0e-9 ); + + // Set the maximum iterations. + solver->setMaxIter( 2000 ); + + // Set the print level. + solver->setPrintLevel( 2 ); + + // Create a preconditioner. + if ( "none" != precond_type ) + { + auto preconditioner = + createHypreSemiStructuredSolver( + precond_type, *vector_layout, true ); + solver->setPreconditioner( preconditioner ); + } + + // Setup the problem. + solver->setup(); + + // Solve the problem. + solver->solve( *rhs, *lhs, 3 ); + + // Create a solver reference for comparison. + vector_layout = createArrayLayout( local_mesh, 1, Cell() ); + + auto lhs_ref = createArray( "lhs_ref", vector_layout ); + ArrayOp::assign( *lhs_ref, 0.0, Own() ); + + auto rhs_ref = createArray( "rhs_ref", vector_layout ); + ArrayOp::assign( *rhs_ref, 1.0, Own() ); + + auto ref_solver = + createReferenceConjugateGradient( *vector_layout ); + ref_solver->setMatrixStencil( stencil ); + const auto& ref_entries = ref_solver->getMatrixValues(); + auto matrix_view = ref_entries.view(); + auto global_space = local_mesh->indexSpace( Own(), Cell(), Global() ); + int ncell_i = global_grid->globalNumEntity( Cell(), Dim::I ); + int ncell_j = global_grid->globalNumEntity( Cell(), Dim::J ); + int ncell_k = global_grid->globalNumEntity( Cell(), Dim::K ); + Kokkos::parallel_for( + "fill_ref_entries", + createExecutionPolicy( owned_space, TEST_EXECSPACE() ), + KOKKOS_LAMBDA( const int i, const int j, const int k ) { + int gi = i + global_space.min( Dim::I ) - owned_space.min( Dim::I ); + int gj = j + global_space.min( Dim::J ) - owned_space.min( Dim::J ); + int gk = k + global_space.min( Dim::K ) - owned_space.min( Dim::K ); + matrix_view( i, j, k, 0 ) = 6.0; + matrix_view( i, j, k, 1 ) = ( gi - 1 >= 0 ) ? -1.0 : 0.0; + matrix_view( i, j, k, 2 ) = ( gi + 1 < ncell_i ) ? -1.0 : 0.0; + matrix_view( i, j, k, 3 ) = ( gj - 1 >= 0 ) ? -1.0 : 0.0; + matrix_view( i, j, k, 4 ) = ( gj + 1 < ncell_j ) ? -1.0 : 0.0; + matrix_view( i, j, k, 5 ) = ( gk - 1 >= 0 ) ? -1.0 : 0.0; + matrix_view( i, j, k, 6 ) = ( gk + 1 < ncell_k ) ? -1.0 : 0.0; + } ); + + std::vector> diag_stencil = { { 0, 0, 0 } }; + ref_solver->setPreconditionerStencil( diag_stencil ); + const auto& preconditioner_entries = ref_solver->getPreconditionerValues(); + auto preconditioner_view = preconditioner_entries.view(); + Kokkos::parallel_for( + "fill_preconditioner_entries", + createExecutionPolicy( owned_space, TEST_EXECSPACE() ), + KOKKOS_LAMBDA( const int i, const int j, const int k ) { + preconditioner_view( i, j, k, 0 ) = 1.0 / 6.0; + } ); + + ref_solver->setTolerance( 1.0e-11 ); + ref_solver->setPrintLevel( 1 ); + ref_solver->setup(); + ref_solver->solve( *rhs_ref, *lhs_ref ); + + // Check the results. + auto lhs_host = + Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace(), lhs->view() ); + auto lhs_ref_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), lhs_ref->view() ); + for ( int i = owned_space.min( Dim::I ); i < owned_space.max( Dim::I ); + ++i ) + for ( int j = owned_space.min( Dim::J ); j < owned_space.max( Dim::J ); + ++j ) + for ( int k = owned_space.min( Dim::K ); + k < owned_space.max( Dim::K ); ++k ) + for ( int var = 0; var < 3; ++var ) + EXPECT_FLOAT_EQ( lhs_host( i, j, k, var ), + lhs_ref_host( i, j, k, 0 ) ); + + // Setup the problem again. We would need to do this if we changed the + // matrix entries. + solver->setup(); + + // Solve the problem again + ArrayOp::assign( *rhs, 2.0, Own() ); + ArrayOp::assign( *lhs, 0.0, Own() ); + solver->solve( *rhs, *lhs, 3 ); + + // Compute another reference solution. + ArrayOp::assign( *lhs_ref, 0.0, Own() ); + ArrayOp::assign( *rhs_ref, 2.0, Own() ); + ref_solver->solve( *rhs_ref, *lhs_ref ); + + // Check the results again + lhs_host = + Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace(), lhs->view() ); + lhs_ref_host = Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace(), + lhs_ref->view() ); + for ( int i = owned_space.min( Dim::I ); i < owned_space.max( Dim::I ); + ++i ) + for ( int j = owned_space.min( Dim::J ); j < owned_space.max( Dim::J ); + ++j ) + for ( int k = owned_space.min( Dim::K ); + k < owned_space.max( Dim::K ); ++k ) + EXPECT_FLOAT_EQ( lhs_host( i, j, k, 1 ), + lhs_ref_host( i, j, k, 0 ) ); + + HYPRE_Finalize(); +} + +//---------------------------------------------------------------------------// +// RUN TESTS +//---------------------------------------------------------------------------// +TEST( semi_structured_solver, pcg_none_test ) +{ + poissonTest( "PCG", "none", TEST_MEMSPACE{} ); +} + +TEST( semi_structured_solver, gmres_none_test ) +{ + poissonTest( "GMRES", "none", TEST_MEMSPACE{} ); +} + +TEST( semi_structured_solver, bicgstab_none_test ) +{ + poissonTest( "BiCGSTAB", "none", TEST_MEMSPACE{} ); +} + +TEST( semi_structured_solver, pcg_diag_test ) +{ + poissonTest( "PCG", "Diagonal", TEST_MEMSPACE{} ); +} + +TEST( semi_structured_solver, gmres_diag_test ) +{ + poissonTest( "GMRES", "Diagonal", TEST_MEMSPACE{} ); +} + +TEST( semi_structured_solver, bicgstab_diag_test ) +{ + poissonTest( "BiCGSTAB", "Diagonal", TEST_MEMSPACE{} ); +} + +//---------------------------------------------------------------------------// + +} // end namespace Test diff --git a/example/cajita_tutorial/11_semi_structured_solver_multi_variate/CMakeLists.txt b/example/cajita_tutorial/11_semi_structured_solver_multi_variate/CMakeLists.txt new file mode 100644 index 000000000..93d592793 --- /dev/null +++ b/example/cajita_tutorial/11_semi_structured_solver_multi_variate/CMakeLists.txt @@ -0,0 +1,15 @@ +############################################################################ +# Copyright (c) 2018-2022 by the Cabana authors # +# All rights reserved. # +# # +# This file is part of the Cabana library. Cabana is distributed under a # +# BSD 3-clause license. For the licensing terms see the LICENSE file in # +# the top-level directory. # +# # +# SPDX-License-Identifier: BSD-3-Clause # +############################################################################ + +add_executable(HypreSemiStructuredSolverMulti hypre_semi_structured_solver_multi_example.cpp) +target_link_libraries(HypreSemiStructuredSolverMulti Cajita) +add_test(NAME Cajita_tutorial_11_hypre_SS COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_MAX_NUMPROCS} ${MPIEXEC_PREFLAGS} $ ${MPIEXEC_POSTFLAGS}) +set_tests_properties(Cajita_tutorial_11_hypre_SS PROPERTIES PROCESSORS ${MPIEXEC_MAX_NUMPROCS}) diff --git a/example/cajita_tutorial/11_semi_structured_solver_multi_variate/hypre_semi_structured_solver_multi_example.cpp b/example/cajita_tutorial/11_semi_structured_solver_multi_variate/hypre_semi_structured_solver_multi_example.cpp new file mode 100644 index 000000000..5a2c3132a --- /dev/null +++ b/example/cajita_tutorial/11_semi_structured_solver_multi_variate/hypre_semi_structured_solver_multi_example.cpp @@ -0,0 +1,200 @@ +/**************************************************************************** + * Copyright (c) 2018-2022 by the Cabana authors * + * All rights reserved. * + * * + * This file is part of the Cabana library. Cabana is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +#include + +#include + +#include + +#include +#include + +//---------------------------------------------------------------------------// +// HYPRE Semi-Structured Solver Example +//---------------------------------------------------------------------------// +void hypreSemiStructuredSolverExample() +{ + /* + In this example we will demonstrate building a HYPRE semi-structured + solver that solves 3, independent, Poisson equations with designated + solution tolerance, + + Laplacian( lhs ) = rhs, + + This is discretized at {i,j,k} + + Laplacian( lhs )_{i,j,k} = rhs_{i,j,k}, + + which includes 7 stencils at current {i,j,k} + + { 0, 0, 0 }, { -1, 0, 0 }, { 1, 0, 0 }, { 0, -1, 0 }, + { 0, 1, 0 }, { 0, 0, -1 }, { 0, 0, 1 } + + You can try one of the following solver type and preconditioner type + + solver type : PCG, GMRES, BiCGSTAB, + preconditioner type : none, Diagonal + */ + + std::cout << "Cajita HYPRE Semi-Structured Solver Example\n" << std::endl; + + /* + As with all Cajita examples, we start by defining everything from the + global mesh to the local grid. + */ + using MemorySpace = Kokkos::HostSpace; + using ExecutionSpace = Kokkos::DefaultHostExecutionSpace; + + // Create the global grid. + double cell_size = 0.25; + std::array is_dim_periodic = { false, false, false }; + std::array global_low_corner = { -1.0, -2.0, -1.0 }; + std::array global_high_corner = { 1.0, 1.0, 0.5 }; + auto global_mesh = Cajita::createUniformGlobalMesh( + global_low_corner, global_high_corner, cell_size ); + + // Create the global grid. + Cajita::DimBlockPartitioner<3> partitioner; + auto global_grid = Cajita::createGlobalGrid( MPI_COMM_WORLD, global_mesh, + is_dim_periodic, partitioner ); + + // Create a local grid. + auto local_mesh = createLocalGrid( global_grid, 1 ); + auto owned_space = local_mesh->indexSpace( Cajita::Own(), Cajita::Cell(), + Cajita::Local() ); + + /************************************************************************/ + + // Create the RHS. + auto vector_layout = createArrayLayout( local_mesh, 3, Cajita::Cell() ); + auto rhs = Cajita::createArray( "rhs", vector_layout ); + Cajita::ArrayOp::assign( *rhs, 1.0, Cajita::Own() ); + + // Create the LHS. + auto lhs = Cajita::createArray( "lhs", vector_layout ); + Cajita::ArrayOp::assign( *lhs, 0.0, Cajita::Own() ); + + // Create a solver. + auto solver = Cajita::createHypreSemiStructuredSolver( + "PCG", *vector_layout, false, 3 ); + + // Create a 7-point 3d laplacian stencil. + std::vector> stencil = { + { 0, 0, 0 }, { -1, 0, 0 }, { 1, 0, 0 }, { 0, -1, 0 }, + { 0, 1, 0 }, { 0, 0, -1 }, { 0, 0, 1 } }; + + solver->createMatrixStencil( 3, 0, 3, { 7, 0, 0 } ); + solver->createMatrixStencil( 3, 1, 3, { 0, 7, 0 } ); + solver->createMatrixStencil( 3, 2, 3, { 0, 0, 7 } ); + for ( int v = 0; v < 3; ++v ) + { + solver->setMatrixStencil( stencil, v, v ); + } + + solver->setSolverGraph( 3 ); + + // Create the matrix entries. The stencil is defined over cells. + auto matrix_entry_layout = + createArrayLayout( local_mesh, 7, Cajita::Cell() ); + auto matrix_entries = Cajita::createArray( + "matrix_entries", matrix_entry_layout ); + auto entry_view = matrix_entries->view(); + Kokkos::parallel_for( + "fill_matrix_entries", + createExecutionPolicy( owned_space, ExecutionSpace() ), + KOKKOS_LAMBDA( const int i, const int j, const int k ) { + entry_view( i, j, k, 0 ) = 6.0; + entry_view( i, j, k, 1 ) = -1.0; + entry_view( i, j, k, 2 ) = -1.0; + entry_view( i, j, k, 3 ) = -1.0; + entry_view( i, j, k, 4 ) = -1.0; + entry_view( i, j, k, 5 ) = -1.0; + entry_view( i, j, k, 6 ) = -1.0; + } ); + + solver->initializeHypreMatrix(); + + for ( int v_h = 0; v_h < 3; ++v_h ) + { + solver->setMatrixValues( *matrix_entries, v_h, v_h ); + } + + // The desired tolerance must be set for each solve. + solver->setTolerance( 1.0e-9 ); + + // Set the maximum iterations. + solver->setMaxIter( 2000 ); + + /* + The print level defines the information output from HYPRE during the solve + */ + solver->setPrintLevel( 2 ); + + /* + Create a preconditioner - in this case we use Diagonal + */ + std::string precond_type = "Diagonal"; + auto preconditioner = + Cajita::createHypreSemiStructuredSolver( + precond_type, *vector_layout, true, 3 ); + solver->setPreconditioner( preconditioner ); + + // Setup the problem - this is necessary before solving. + solver->setup(); + + // Now solve the problem. + solver->solve( *rhs, *lhs, 3 ); + + /* + Setup the problem again. We would need to do this if we changed the matrix + entries, but in this case we just leave it unchanged. + */ + solver->setup(); + // Reset to the same initial condition and solve the problem again. + Cajita::ArrayOp::assign( *rhs, 2.0, Cajita::Own() ); + Cajita::ArrayOp::assign( *lhs, 0.0, Cajita::Own() ); + solver->solve( *rhs, *lhs, 3 ); +} + +//---------------------------------------------------------------------------// +// Main. +//---------------------------------------------------------------------------// +int main( int argc, char* argv[] ) +{ + // MPI only needed to create the grid/mesh. Not intended to be run with + // multiple ranks. + MPI_Init( &argc, &argv ); + { + Kokkos::ScopeGuard scope_guard( argc, argv ); + + /* + The hypre solver capabilities used by Cabana must be initialized and + finalized. HYPRE_Init() initializes hypre. A call to HYPRE_Init() must + be included before any hypre calls occur + */ + HYPRE_Init(); + + hypreSemiStructuredSolverExample(); + + /* + The hypre solver capabilities used by Cabana must be initialized and + finalized. HYPRE_Finalize() finalizes hypre. A call to + HYPRE_Finalize() should not occur before all calls to hypre + capabilites are finished. + */ + HYPRE_Finalize(); + } + MPI_Finalize(); + + return 0; +} +//---------------------------------------------------------------------------// diff --git a/example/cajita_tutorial/11_structured_solver_hypre/hypre_structured_solver_example.cpp b/example/cajita_tutorial/11_structured_solver_hypre/hypre_structured_solver_example.cpp index aeaa499eb..92a84fbe8 100644 --- a/example/cajita_tutorial/11_structured_solver_hypre/hypre_structured_solver_example.cpp +++ b/example/cajita_tutorial/11_structured_solver_hypre/hypre_structured_solver_example.cpp @@ -113,6 +113,8 @@ void hypreStructuredSolverExample() solver->setMatrixValues( *matrix_entries ); + solver->printMatrix( "Struct.mat" ); + // The desired tolerance must be set for each solve. solver->setTolerance( 1.0e-9 ); diff --git a/example/cajita_tutorial/CMakeLists.txt b/example/cajita_tutorial/CMakeLists.txt index 75c13da54..b96346390 100644 --- a/example/cajita_tutorial/CMakeLists.txt +++ b/example/cajita_tutorial/CMakeLists.txt @@ -24,6 +24,7 @@ endif() add_subdirectory(11_structured_solver) if(Cabana_ENABLE_HYPRE AND (NOT Kokkos_ENABLE_CUDA AND NOT Kokkos_ENABLE_HIP AND NOT Kokkos_ENABLE_SYCL)) add_subdirectory(11_structured_solver_hypre) + add_subdirectory(11_semi_structured_solver_multi_variate) endif() add_subdirectory(12_halo) if(Cabana_ENABLE_ALL)