mirror of
https://github.com/OPM/opm-simulators.git
synced 2024-12-25 16:51:00 -06:00
290 lines
9.0 KiB
C++
290 lines
9.0 KiB
C++
/*
|
|
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 <opm/common/utility/platform_dependent/disable_warnings.h>
|
|
#include <petsc.h>
|
|
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
|
|
#include <opm/common/ErrorMacros.hpp>
|
|
|
|
namespace Opm
|
|
{
|
|
|
|
namespace{
|
|
|
|
class KSPTypeMap {
|
|
public:
|
|
explicit
|
|
KSPTypeMap(const std::string& default_type = "gmres")
|
|
: default_type_(default_type)
|
|
{
|
|
// g++-4.4 has problems converting const char* to char*
|
|
// The problem is caused by the mapped type being PCType
|
|
// which (at least in PETSc 3.2) is char* because of C
|
|
// (in the header there is "#define PCType character*(80)").
|
|
// and the KSP... defines being const char* (because of C++).
|
|
type_map_["richardson"] = KSPRICHARDSON;
|
|
// Not available in PETSC 3.2 on Debian
|
|
//type_map_["chebyshev"] = KSPCHEBYSHEV;
|
|
type_map_["cg"] = KSPCG;
|
|
type_map_["bicgs"] = KSPBICG;
|
|
type_map_["gmres"] = KSPGMRES;
|
|
type_map_["fgmres"] = KSPFGMRES;
|
|
type_map_["dgmres"] = KSPDGMRES;
|
|
type_map_["gcr"] = KSPGCR;
|
|
type_map_["bcgs"] = KSPBCGS;
|
|
type_map_["cgs"] = KSPCGS;
|
|
type_map_["tfqmr"] = KSPTFQMR;
|
|
type_map_["tcqmr"] = KSPTCQMR;
|
|
type_map_["cr"] = KSPCR;
|
|
type_map_["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_["jacobi"] = PCJACOBI;
|
|
type_map_["bjacobi"] = PCBJACOBI;
|
|
type_map_["sor"] = PCSOR;
|
|
type_map_["eisenstat"] = PCEISENSTAT;
|
|
type_map_["icc"] = PCICC;
|
|
type_map_["ilu"] = PCILU;
|
|
type_map_["asm"] = PCASM;
|
|
type_map_["gamg"] = PCGAMG;
|
|
type_map_["ksp"] = PCKSP;
|
|
type_map_["composite"] = PCCOMPOSITE;
|
|
type_map_["lu"] = PCLU;
|
|
type_map_["cholesky"] = PCCHOLESKY;
|
|
type_map_["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 );
|
|
}
|
|
};
|
|
|
|
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, const_cast<int*>(ia), const_cast<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;
|
|
|
|
#if PETSC_VERSION_MAJOR <= 3 && PETSC_VERSION_MINOR < 5
|
|
KSPSetOperators( t.ksp, t.A, t.A, DIFFERENT_NONZERO_PATTERN );
|
|
#else
|
|
KSPSetOperators( t.ksp, t.A, t.A );
|
|
KSPSetReusePreconditioner(t.ksp, PETSC_FALSE);
|
|
#endif
|
|
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, (double)residual );
|
|
CHKERRXX( err );
|
|
}
|
|
|
|
} // anonymous namespace.
|
|
|
|
LinearSolverPetsc::LinearSolverPetsc(const 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
|