diff --git a/opm/core/linalg/LinearSolverFactory.cpp b/opm/core/linalg/LinearSolverFactory.cpp index 207c4ea54..0e22910b2 100644 --- a/opm/core/linalg/LinearSolverFactory.cpp +++ b/opm/core/linalg/LinearSolverFactory.cpp @@ -31,6 +31,10 @@ #include #endif +#if HAVE_PETSC +#include +#endif + #include #include #include @@ -45,8 +49,10 @@ namespace Opm solver_.reset(new LinearSolverUmfpack); #elif HAVE_DUNE_ISTL solver_.reset(new LinearSolverIstl); +#elif HAVE_PETSC + solver_.reset(new LinearSolverPetsc); #else - OPM_THROW(std::runtime_error, "No linear solver available, you must have UMFPACK or dune-istl installed to use LinearSolverFactory."); + OPM_THROW(std::runtime_error, "No linear solver available, you must have UMFPACK , dune-istl or Petsc installed to use LinearSolverFactory."); #endif } @@ -69,6 +75,11 @@ namespace Opm solver_.reset(new LinearSolverIstl(param)); #endif } + else if (ls == "petsc"){ +#if HAVE_PETSC + solver_.reset(new LinearSolverPetsc(param)); +#endif + } else { OPM_THROW(std::runtime_error, "Linear solver " << ls << " is unknown."); diff --git a/opm/core/linalg/LinearSolverFactory.hpp b/opm/core/linalg/LinearSolverFactory.hpp index f27140bc0..92f70c337 100644 --- a/opm/core/linalg/LinearSolverFactory.hpp +++ b/opm/core/linalg/LinearSolverFactory.hpp @@ -42,16 +42,19 @@ namespace Opm /// Construct from parameters. /// The accepted parameters are (default) (allowed values): - /// linsolver ("umfpack") ("umfpack", "istl") + /// linsolver ("umfpack") ("umfpack", "istl", "petsc") /// For the umfpack solver to be available, this class must be /// compiled with UMFPACK support, as indicated by the /// variable HAVE_SUITESPARSE_UMFPACK_H in config.h. /// For the istl solver to be available, this class must be /// compiled with dune-istl support, as indicated by the /// variable HAVE_DUNE_ISTL in config.h. + /// For the petsc solver to be available, this class must be + /// compiled with petsc support, as indicated by the + /// variable HAVE_PETSC in config.h. /// Any further parameters are passed on to the constructors - /// of the actual solver used, see LinearSolverUmfpack - /// and LinearSolverIstl for details. + /// of the actual solver used, see LinearSolverUmfpack, + /// LinearSolverIstl and LinearSolverPetsc for details. LinearSolverFactory(const parameter::ParameterGroup& param); /// Destructor. diff --git a/opm/core/linalg/LinearSolverPetsc.cpp b/opm/core/linalg/LinearSolverPetsc.cpp new file mode 100644 index 000000000..1b4fe4fb1 --- /dev/null +++ b/opm/core/linalg/LinearSolverPetsc.cpp @@ -0,0 +1,277 @@ +/* + Copyright 2014 SINTEF ICT, Applied Mathematics. + Copyright 2014 STATOIL ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include "config.h" +#include +#include +#include +#define PETSC_CLANGUAGE_CXX 1 //enable CHKERRXX macro. +#include +#include + +namespace Opm +{ + +namespace{ + + class KSPTypeMap { + public: + explicit + KSPTypeMap(const std::string& default_type = "gmres") + : default_type_(default_type) + { + type_map_.insert(std::make_pair("richardson", KSPRICHARDSON)); + type_map_.insert(std::make_pair("chebyshev", KSPCHEBYSHEV)); + type_map_.insert(std::make_pair("cg", KSPCG)); + type_map_.insert(std::make_pair("bicgs", KSPBICG)); + type_map_.insert(std::make_pair("gmres", KSPGMRES)); + type_map_.insert(std::make_pair("fgmres", KSPFGMRES)); + type_map_.insert(std::make_pair("dgmres", KSPDGMRES)); + type_map_.insert(std::make_pair("gcr", KSPGCR)); + type_map_.insert(std::make_pair("bcgs", KSPBCGS)); + type_map_.insert(std::make_pair("cgs", KSPCGS)); + type_map_.insert(std::make_pair("tfqmr", KSPTFQMR)); + type_map_.insert(std::make_pair("tcqmr", KSPTCQMR)); + type_map_.insert(std::make_pair("cr", KSPCR)); + type_map_.insert(std::make_pair("preonly", KSPPREONLY)); + } + + KSPType + find(const std::string& type) const + { + Map::const_iterator it = type_map_.find(type); + + if (it == type_map_.end()) { + it = type_map_.find(default_type_); + } + + if (it == type_map_.end()) { + OPM_THROW(std::runtime_error, "Unknown KSPType: '" << type << "'"); + } + return it->second; + } + private: + typedef std::unordered_map Map; + + std::string default_type_; + Map type_map_; + }; + + + class PCTypeMap { + public: + explicit + PCTypeMap(const std::string& default_type = "jacobi") + : default_type_(default_type) + { + type_map_.insert(std::make_pair("jacobi", PCJACOBI)); + type_map_.insert(std::make_pair("bjacobi", PCBJACOBI)); + type_map_.insert(std::make_pair("sor", PCSOR)); + type_map_.insert(std::make_pair("eisenstat", PCEISENSTAT)); + type_map_.insert(std::make_pair("icc", PCICC)); + type_map_.insert(std::make_pair("ilu", PCILU)); + type_map_.insert(std::make_pair("asm", PCASM)); + type_map_.insert(std::make_pair("gamg", PCGAMG)); + type_map_.insert(std::make_pair("ksp", PCKSP)); + type_map_.insert(std::make_pair("composite", PCCOMPOSITE)); + type_map_.insert(std::make_pair("lu", PCLU)); + type_map_.insert(std::make_pair("cholesky", PCCHOLESKY)); + type_map_.insert(std::make_pair("none", PCNONE)); + } + + PCType + find(const std::string& type) const + { + Map::const_iterator it = type_map_.find(type); + + if (it == type_map_.end()) { + it = type_map_.find(default_type_); + } + + if (it == type_map_.end()) { + OPM_THROW(std::runtime_error, "Unknown PCType: '" << type << "'"); + } + + return it->second; + } + private: + typedef std::unordered_map Map; + + std::string default_type_; + Map type_map_; + }; + + struct OEM_DATA { + /* Convenience struct to handle automatic (de)allocation of some useful + * variables, as well as group them up for easier parameter passing + */ + Vec x; + Vec b; + Mat A; + KSP ksp; + PC preconditioner; + + OEM_DATA( const int size ) { + VecCreate( PETSC_COMM_WORLD, &b ); + auto err = VecSetSizes( b, PETSC_DECIDE, size ); + CHKERRXX( err ); + VecSetFromOptions( b ); + + KSPCreate( PETSC_COMM_WORLD, &ksp ); + } + + ~OEM_DATA() { + VecDestroy( &x ); + VecDestroy( &b ); + MatDestroy( &A ); + KSPDestroy( &ksp ); + PCDestroy( &preconditioner ); + } + }; + + Vec to_petsc_vec( const double* x, int size ) { + PetscScalar* vec; + Vec v; + + VecCreate( PETSC_COMM_WORLD, &v ); + VecSetSizes( v, PETSC_DECIDE, size ); + VecSetFromOptions( v ); + + VecGetArray( v, &vec ); + std::memcpy( vec, x, size * sizeof( double ) ); + + VecRestoreArray( v, &vec ); + return v; + } + + void from_petsc_vec( double* x, Vec v ) { + if( !v ) OPM_THROW( std::runtime_error, + "PETSc CopySolution: Invalid PETSc vector." ); + + PetscScalar* vec; + PetscInt size; + + VecGetLocalSize( v, &size ); + VecGetArray( v, &vec ); + + std::memcpy( x, vec, size * sizeof( double ) ); + VecRestoreArray( v, &vec ); + } + + Mat to_petsc_mat( const int size, const int nonzeros, + const int* ia, const int* ja, const double* sa ) { + + Mat A; + auto err = MatCreateSeqAIJWithArrays( PETSC_COMM_WORLD, size, size, (int*)ia, (int*)ja, (double*)sa, &A ); + CHKERRXX( err ); + return A; + } + + + void solve_system( OEM_DATA& t, KSPType method, PCType pcname, + double rtol, double atol, double dtol, int maxits, int ksp_view ) { + PetscInt its; + PetscReal residual; + KSPConvergedReason reason; + + KSPSetOperators( t.ksp, t.A, t.A, DIFFERENT_NONZERO_PATTERN ); + KSPGetPC( t.ksp, &t.preconditioner ); + auto err = KSPSetType( t.ksp, method ); + CHKERRXX( err ); + err = PCSetType( t.preconditioner, pcname ); + CHKERRXX( err ); + err = KSPSetTolerances( t.ksp, rtol, atol, dtol, maxits ); + CHKERRXX( err ); + err = KSPSetFromOptions( t.ksp ); + CHKERRXX( err ); + KSPSetInitialGuessNonzero( t.ksp, PETSC_FALSE ); + KSPSolve( t.ksp, t.x, t.b ); + KSPGetConvergedReason( t.ksp, &reason ); + KSPGetIterationNumber( t.ksp, &its ); + KSPGetResidualNorm( t.ksp, &residual ); + + if( ksp_view ) + KSPView( t.ksp, PETSC_VIEWER_STDOUT_WORLD ); + + err = PetscPrintf( PETSC_COMM_WORLD, "KSP Iterations %D, Final Residual %G\n", its, residual ); + CHKERRXX( err ); + } + +} // anonymous namespace. + + LinearSolverPetsc::LinearSolverPetsc(const parameter::ParameterGroup& param) + : ksp_type_( param.getDefault( std::string( "ksp_type" ), std::string( "gmres" ) ) ) + , pc_type_( param.getDefault( std::string( "pc_type" ), std::string( "sor" ) ) ) + , ksp_view_( param.getDefault( std::string( "ksp_view" ), int( false ) ) ) + , rtol_( param.getDefault( std::string( "ksp_rtol" ), 1e-5 ) ) + , atol_( param.getDefault( std::string( "ksp_atol" ), 1e-50 ) ) + , dtol_( param.getDefault( std::string( "ksp_dtol" ), 1e5 ) ) + , maxits_( param.getDefault( std::string( "ksp_max_it" ), 1e5 ) ) + { + int argc = 0; + char** argv = NULL; + PetscInitialize(&argc, &argv, (char*)0, "Petsc interface for OPM!\n"); + } + + + LinearSolverPetsc::~LinearSolverPetsc() + { + PetscFinalize(); + } + + + LinearSolverInterface::LinearSolverReport + LinearSolverPetsc::solve(const int size, + const int nonzeros, + const int* ia, + const int* ja, + const double* sa, + const double* rhs, + double* solution, + const boost::any&) const + { + KSPTypeMap ksp(ksp_type_); + KSPType ksp_type = ksp.find(ksp_type_); + PCTypeMap pc(pc_type_); + PCType pc_type = pc.find(pc_type_); + + OEM_DATA t( size ); + t.A = to_petsc_mat( size, nonzeros, ia, ja, sa ); + t.x = to_petsc_vec( rhs, size ); + + solve_system( t, ksp_type, pc_type, rtol_, atol_, dtol_, maxits_, ksp_view_ ); + from_petsc_vec( solution, t.b ); + + LinearSolverReport rep = {}; + rep.converged = true; + return rep; + } + + void LinearSolverPetsc::setTolerance(const double /*tol*/) + { + } + + double LinearSolverPetsc::getTolerance() const + { + return -1.; + } + +} // namespace Opm + diff --git a/opm/core/linalg/LinearSolverPetsc.hpp b/opm/core/linalg/LinearSolverPetsc.hpp new file mode 100644 index 000000000..9750ca3d7 --- /dev/null +++ b/opm/core/linalg/LinearSolverPetsc.hpp @@ -0,0 +1,92 @@ +/* + Copyright 2014 SINTEF ICT, Applied Mathematics. + Copyright 2014 STATOIL ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_LINEARSOLVERPETSC_HEADER_INCLUDED +#define OPM_LINEARSOLVERPETSC_HEADER_INCLUDED +#include +#include +#include + +namespace Opm +{ + + + /// Concrete class encapsulating some Petsc linear solvers. + class LinearSolverPetsc : public LinearSolverInterface + { + public: + /// Default constructor. + /// Declared, but not implemented. Petsc can only be created through + /// the ParameterGroup constructor, everything else is an error. This way + /// the error is caught compile time and not rune time, which is nice as + /// it is a static error. + LinearSolverPetsc(); + + /// Construct from parameters + /// Accepted parameters are, with defaults, listed in the + /// default constructor. + LinearSolverPetsc(const parameter::ParameterGroup& param); + + /// Destructor. + virtual ~LinearSolverPetsc(); + + using LinearSolverInterface::solve; + + /// Solve a linear system, with a matrix given in compressed sparse row format. + /// \param[in] size # of rows in matrix + /// \param[in] nonzeros # of nonzeros elements in matrix + /// \param[in] ia array of length (size + 1) containing start and end indices for each row + /// \param[in] ja array of length nonzeros containing column numbers for the nonzero elements + /// \param[in] sa array of length nonzeros containing the values of the nonzero elements + /// \param[in] rhs array of length size containing the right hand side + /// \param[inout] solution array of length size to which the solution will be written, may also be used + /// as initial guess by iterative solvers. + virtual LinearSolverReport solve(const int size, + const int nonzeros, + const int* ia, + const int* ja, + const double* sa, + const double* rhs, + double* solution, + const boost::any&) const; + + /// Set tolerance for the residual in dune istl linear solver. + /// \param[in] tol tolerance value + virtual void setTolerance(const double tol); + + /// Get tolerance ofthe linear solver. + /// \param[out] tolerance value + virtual double getTolerance() const; + private: + std::string ksp_type_; + std::string pc_type_; + int ksp_view_; + double rtol_; + double atol_; + double dtol_; + int maxits_; + }; + + +} // namespace Opm + + + +#endif // OPM_LINEARSOLVERPETSC_HEADER_INCLUDED diff --git a/tests/test_linearsolver.cpp b/tests/test_linearsolver.cpp index 184b7765d..6021218e6 100644 --- a/tests/test_linearsolver.cpp +++ b/tests/test_linearsolver.cpp @@ -184,3 +184,16 @@ BOOST_AUTO_TEST_CASE(KAMGTest) } #endif #endif + +#if HAVE_PETSC +BOOST_AUTO_TEST_CASE(PETScTest) +{ + Opm::parameter::ParameterGroup param; + param.insertParameter(std::string("linsolver"), std::string("petsc")); + param.insertParameter(std::string("ksp_type"), std::string("cg")); + param.insertParameter(std::string("pc_type"), std::string("jacobi")); + param.insertParameter(std::string("ksp_rtol"), std::string("1e-10")); + param.insertParameter(std::string("ksp_view"), std::string("0")); + run_test(param); +} +#endif