OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #702 from jorgekva/linsolver_petsc

Browse Source 
PETSc support in OPM
This commit is contained in:
Atgeirr Flø Rasmussen
2014-12-19 11:25:58 +01:00
parent c2116149c9 263e1a6928
commit dea6e8c044
5 changed files with 400 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files

13
opm/core/linalg/LinearSolverFactory.cpp
View File

@@ -31,6 +31,10 @@
#include <opm/core/linalg/LinearSolverIstl.hpp>
#endif
#if HAVE_PETSC
#include <opm/core/linalg/LinearSolverPetsc.hpp>
#endif
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <string>
@@ -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.");

9
opm/core/linalg/LinearSolverFactory.hpp
View File

@@ -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.

277
opm/core/linalg/LinearSolverPetsc.cpp Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <cstring>
#include <opm/core/linalg/LinearSolverPetsc.hpp>
#include <unordered_map>
#define PETSC_CLANGUAGE_CXX 1 //enable CHKERRXX macro.
#include <petsc.h>
#include <opm/core/utility/ErrorMacros.hpp>
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<std::string, KSPType> 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<std::string, PCType> 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

92
opm/core/linalg/LinearSolverPetsc.hpp Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_LINEARSOLVERPETSC_HEADER_INCLUDED
#define OPM_LINEARSOLVERPETSC_HEADER_INCLUDED
#include <opm/core/linalg/LinearSolverInterface.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <string>
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

13
tests/test_linearsolver.cpp
View File

@@ -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