OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-19 05:53:28 -06:00
Code Issues Packages Projects Releases Wiki Activity
e69d92cca3
opm-simulators/opm/core/linalg/LinearSolverPetsc.cpp

292 lines
9.6 KiB
C++
Raw Normal View History

Add Petsc to OPM.
2014-05-16 01:41:51 -05:00
/*
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"
Removes call_petsc.c and calls the lib from C++ call_petsc.c was really a thin C wrapper around the call to petsc itself and turns out was mostly unnecessary C++ emulation. This removes the file entirely and ports its functionality into LinearSolverPetsc.cpp. All features from the file should now be more readable as well as properly utilising modern C++ features. The patch uses the CHKERRXX macro from petsc to handle errors reported by petsc, and currently does not handle this and give the control back to OPM's error/throw system.
2014-09-22 07:52:51 -05:00
#include <cstring>
Add Petsc to OPM.
2014-05-16 01:41:51 -05:00
#include <opm/core/linalg/LinearSolverPetsc.hpp>
use unordered_map for "string-enum" translations.
2014-05-22 22:09:04 -05:00
#include <unordered_map>
Removes call_petsc.c and calls the lib from C++ call_petsc.c was really a thin C wrapper around the call to petsc itself and turns out was mostly unnecessary C++ emulation. This removes the file entirely and ports its functionality into LinearSolverPetsc.cpp. All features from the file should now be more readable as well as properly utilising modern C++ features. The patch uses the CHKERRXX macro from petsc to handle errors reported by petsc, and currently does not handle this and give the control back to OPM's error/throw system.
2014-09-22 07:52:51 -05:00
#define PETSC_CLANGUAGE_CXX 1 //enable CHKERRXX macro.
Initialize Petsc from constructor. Remove private members for initializing Petsc.
2014-05-16 01:53:49 -05:00
#include <petsc.h>
Add Petsc to OPM.
2014-05-16 01:41:51 -05:00
#include <opm/core/utility/ErrorMacros.hpp>
use unordered_map for "string-enum" translations.
2014-05-22 22:09:04 -05:00
Add Petsc to OPM.
2014-05-16 01:41:51 -05:00
namespace Opm
{
use unordered_map for "string-enum" translations.
2014-05-22 22:09:04 -05:00
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_;
};
Removes call_petsc.c and calls the lib from C++ call_petsc.c was really a thin C wrapper around the call to petsc itself and turns out was mostly unnecessary C++ emulation. This removes the file entirely and ports its functionality into LinearSolverPetsc.cpp. All features from the file should now be more readable as well as properly utilising modern C++ features. The patch uses the CHKERRXX macro from petsc to handle errors reported by petsc, and currently does not handle this and give the control back to OPM's error/throw system.
2014-09-22 07:52:51 -05:00
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 rhs;
Vec solution;
Mat A;
KSP ksp;
PC preconditioner;
OEM_DATA( const int size ) {
Reduces CHKERRXX usage to where necessary. The error checking macro makes it harder to read and harder to write, so instead we now only check for functions that can contain errors. Bounds and range checks are handled by PETSc and not OPM.
2014-10-07 08:50:36 -05:00
VecCreate( PETSC_COMM_WORLD, &solution );
auto err = VecSetSizes( solution, PETSC_DECIDE, size );
CHKERRXX( err );
VecSetFromOptions( solution );
VecDuplicate( solution, &rhs );
MatCreate( PETSC_COMM_WORLD, &A );
err = MatSetSizes( A, PETSC_DECIDE, PETSC_DECIDE, size, size );
CHKERRXX( err );
MatSetFromOptions( A );
MatSetUp( A );
Moved KSPCreate to the Data constructor This is where it was always intended to be called, so this fixes a mistake in earlier development.
2014-11-12 13:40:46 -06:00
KSPCreate( PETSC_COMM_WORLD, &ksp );
Removes call_petsc.c and calls the lib from C++ call_petsc.c was really a thin C wrapper around the call to petsc itself and turns out was mostly unnecessary C++ emulation. This removes the file entirely and ports its functionality into LinearSolverPetsc.cpp. All features from the file should now be more readable as well as properly utilising modern C++ features. The patch uses the CHKERRXX macro from petsc to handle errors reported by petsc, and currently does not handle this and give the control back to OPM's error/throw system.
2014-09-22 07:52:51 -05:00
}
~OEM_DATA() {
Reduces CHKERRXX usage to where necessary. The error checking macro makes it harder to read and harder to write, so instead we now only check for functions that can contain errors. Bounds and range checks are handled by PETSc and not OPM.
2014-10-07 08:50:36 -05:00
VecDestroy( &rhs );
VecDestroy( &solution );
MatDestroy( &A );
KSPDestroy( &ksp );
Removes call_petsc.c and calls the lib from C++ call_petsc.c was really a thin C wrapper around the call to petsc itself and turns out was mostly unnecessary C++ emulation. This removes the file entirely and ports its functionality into LinearSolverPetsc.cpp. All features from the file should now be more readable as well as properly utilising modern C++ features. The patch uses the CHKERRXX macro from petsc to handle errors reported by petsc, and currently does not handle this and give the control back to OPM's error/throw system.
2014-09-22 07:52:51 -05:00
}
};
void to_petsc_vec( const double* x, Vec v ) {
if( !v ) OPM_THROW( std::runtime_error,
"PETSc CopySolution: Invalid PETSc vector." );
PetscScalar* vec;
PetscInt size;
Reduces CHKERRXX usage to where necessary. The error checking macro makes it harder to read and harder to write, so instead we now only check for functions that can contain errors. Bounds and range checks are handled by PETSc and not OPM.
2014-10-07 08:50:36 -05:00
VecGetLocalSize( v, &size );
VecGetArray( v, &vec );
Removes call_petsc.c and calls the lib from C++ call_petsc.c was really a thin C wrapper around the call to petsc itself and turns out was mostly unnecessary C++ emulation. This removes the file entirely and ports its functionality into LinearSolverPetsc.cpp. All features from the file should now be more readable as well as properly utilising modern C++ features. The patch uses the CHKERRXX macro from petsc to handle errors reported by petsc, and currently does not handle this and give the control back to OPM's error/throw system.
2014-09-22 07:52:51 -05:00
std::memcpy( vec, x, size * sizeof( double ) );
Reduces CHKERRXX usage to where necessary. The error checking macro makes it harder to read and harder to write, so instead we now only check for functions that can contain errors. Bounds and range checks are handled by PETSc and not OPM.
2014-10-07 08:50:36 -05:00
VecRestoreArray( v, &vec );
Removes call_petsc.c and calls the lib from C++ call_petsc.c was really a thin C wrapper around the call to petsc itself and turns out was mostly unnecessary C++ emulation. This removes the file entirely and ports its functionality into LinearSolverPetsc.cpp. All features from the file should now be more readable as well as properly utilising modern C++ features. The patch uses the CHKERRXX macro from petsc to handle errors reported by petsc, and currently does not handle this and give the control back to OPM's error/throw system.
2014-09-22 07:52:51 -05:00
}
void from_petsc_vec( double* x, Vec v ) {
if( !v ) OPM_THROW( std::runtime_error,
"PETSc CopySolution: Invalid PETSc vector." );
PetscScalar* vec;
PetscInt size;
Reduces CHKERRXX usage to where necessary. The error checking macro makes it harder to read and harder to write, so instead we now only check for functions that can contain errors. Bounds and range checks are handled by PETSc and not OPM.
2014-10-07 08:50:36 -05:00
VecGetLocalSize( v, &size );
VecGetArray( v, &vec );
Removes call_petsc.c and calls the lib from C++ call_petsc.c was really a thin C wrapper around the call to petsc itself and turns out was mostly unnecessary C++ emulation. This removes the file entirely and ports its functionality into LinearSolverPetsc.cpp. All features from the file should now be more readable as well as properly utilising modern C++ features. The patch uses the CHKERRXX macro from petsc to handle errors reported by petsc, and currently does not handle this and give the control back to OPM's error/throw system.
2014-09-22 07:52:51 -05:00
std::memcpy( x, vec, size * sizeof( double ) );
Reduces CHKERRXX usage to where necessary. The error checking macro makes it harder to read and harder to write, so instead we now only check for functions that can contain errors. Bounds and range checks are handled by PETSc and not OPM.
2014-10-07 08:50:36 -05:00
VecRestoreArray( v, &vec );
Removes call_petsc.c and calls the lib from C++ call_petsc.c was really a thin C wrapper around the call to petsc itself and turns out was mostly unnecessary C++ emulation. This removes the file entirely and ports its functionality into LinearSolverPetsc.cpp. All features from the file should now be more readable as well as properly utilising modern C++ features. The patch uses the CHKERRXX macro from petsc to handle errors reported by petsc, and currently does not handle this and give the control back to OPM's error/throw system.
2014-09-22 07:52:51 -05:00
}
void to_petsc_mat( const int size, const int nonzeros,
const int* ia, const int* ja, const double* sa, Mat A ) {
for( int i = 0; i < size; ++i ) {
int nrows = ia[ i + 1 ] - ia[ i ];
if( nrows == 0 ) continue;
Reduces CHKERRXX usage to where necessary. The error checking macro makes it harder to read and harder to write, so instead we now only check for functions that can contain errors. Bounds and range checks are handled by PETSc and not OPM.
2014-10-07 08:50:36 -05:00
for( int j = ia[ i ]; j < ia[ i + 1 ]; ++j ) {
auto err = MatSetValues( A, 1, &i, 1, &ja[ j ], &sa[ j ], INSERT_VALUES );
CHKERRXX( err );
}
Removes call_petsc.c and calls the lib from C++ call_petsc.c was really a thin C wrapper around the call to petsc itself and turns out was mostly unnecessary C++ emulation. This removes the file entirely and ports its functionality into LinearSolverPetsc.cpp. All features from the file should now be more readable as well as properly utilising modern C++ features. The patch uses the CHKERRXX macro from petsc to handle errors reported by petsc, and currently does not handle this and give the control back to OPM's error/throw system.
2014-09-22 07:52:51 -05:00
Reduces CHKERRXX usage to where necessary. The error checking macro makes it harder to read and harder to write, so instead we now only check for functions that can contain errors. Bounds and range checks are handled by PETSc and not OPM.
2014-10-07 08:50:36 -05:00
auto err = MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY );
CHKERRXX( err );
MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY );
Removes call_petsc.c and calls the lib from C++ call_petsc.c was really a thin C wrapper around the call to petsc itself and turns out was mostly unnecessary C++ emulation. This removes the file entirely and ports its functionality into LinearSolverPetsc.cpp. All features from the file should now be more readable as well as properly utilising modern C++ features. The patch uses the CHKERRXX macro from petsc to handle errors reported by petsc, and currently does not handle this and give the control back to OPM's error/throw system.
2014-09-22 07:52:51 -05:00
}
}
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;
Reduces CHKERRXX usage to where necessary. The error checking macro makes it harder to read and harder to write, so instead we now only check for functions that can contain errors. Bounds and range checks are handled by PETSc and not OPM.
2014-10-07 08:50:36 -05:00
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_TRUE );
KSPSolve( t.ksp, t.rhs, t.solution );
KSPGetConvergedReason( t.ksp, &reason );
KSPGetIterationNumber( t.ksp, &its );
KSPGetResidualNorm( t.ksp, &residual );
Removes call_petsc.c and calls the lib from C++ call_petsc.c was really a thin C wrapper around the call to petsc itself and turns out was mostly unnecessary C++ emulation. This removes the file entirely and ports its functionality into LinearSolverPetsc.cpp. All features from the file should now be more readable as well as properly utilising modern C++ features. The patch uses the CHKERRXX macro from petsc to handle errors reported by petsc, and currently does not handle this and give the control back to OPM's error/throw system.
2014-09-22 07:52:51 -05:00
if( ksp_view )
Reduces CHKERRXX usage to where necessary. The error checking macro makes it harder to read and harder to write, so instead we now only check for functions that can contain errors. Bounds and range checks are handled by PETSc and not OPM.
2014-10-07 08:50:36 -05:00
KSPView( t.ksp, PETSC_VIEWER_STDOUT_WORLD );
Removes call_petsc.c and calls the lib from C++ call_petsc.c was really a thin C wrapper around the call to petsc itself and turns out was mostly unnecessary C++ emulation. This removes the file entirely and ports its functionality into LinearSolverPetsc.cpp. All features from the file should now be more readable as well as properly utilising modern C++ features. The patch uses the CHKERRXX macro from petsc to handle errors reported by petsc, and currently does not handle this and give the control back to OPM's error/throw system.
2014-09-22 07:52:51 -05:00
Reduces CHKERRXX usage to where necessary. The error checking macro makes it harder to read and harder to write, so instead we now only check for functions that can contain errors. Bounds and range checks are handled by PETSc and not OPM.
2014-10-07 08:50:36 -05:00
err = PetscPrintf( PETSC_COMM_WORLD, "KSP Iterations %D, Final Residual %G\n", its, residual );
CHKERRXX( err );
Removes call_petsc.c and calls the lib from C++ call_petsc.c was really a thin C wrapper around the call to petsc itself and turns out was mostly unnecessary C++ emulation. This removes the file entirely and ports its functionality into LinearSolverPetsc.cpp. All features from the file should now be more readable as well as properly utilising modern C++ features. The patch uses the CHKERRXX macro from petsc to handle errors reported by petsc, and currently does not handle this and give the control back to OPM's error/throw system.
2014-09-22 07:52:51 -05:00
}
use unordered_map for "string-enum" translations.
2014-05-22 22:09:04 -05:00
} // anonymous namespace.
Add Petsc to OPM.
2014-05-16 01:41:51 -05:00
LinearSolverPetsc::LinearSolverPetsc(const parameter::ParameterGroup& param)
Petsc constructor now uses intialiser list. The previous implementation set plenty of values in the initialization list and immediately overwrote these values with values looked up from the param group. This patch makes it look up the parameteres from the param group argument, making the constructor simpler.
2014-09-22 08:44:22 -05:00
: 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 ) )
Add Petsc to OPM.
2014-05-16 01:41:51 -05:00
{
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
{
use unordered_map for "string-enum" translations.
2014-05-22 22:09:04 -05:00
KSPTypeMap ksp(ksp_type_);
KSPType ksp_type = ksp.find(ksp_type_);
PCTypeMap pc(pc_type_);
PCType pc_type = pc.find(pc_type_);
Removes call_petsc.c and calls the lib from C++ call_petsc.c was really a thin C wrapper around the call to petsc itself and turns out was mostly unnecessary C++ emulation. This removes the file entirely and ports its functionality into LinearSolverPetsc.cpp. All features from the file should now be more readable as well as properly utilising modern C++ features. The patch uses the CHKERRXX macro from petsc to handle errors reported by petsc, and currently does not handle this and give the control back to OPM's error/throw system.
2014-09-22 07:52:51 -05:00
OEM_DATA t( size );
to_petsc_mat( size, nonzeros, ia, ja, sa, t.A );
to_petsc_vec( rhs, t.rhs );
solve_system( t, ksp_type, pc_type, rtol_, atol_, dtol_, maxits_, ksp_view_ );
from_petsc_vec( solution, t.solution );
use unordered_map for "string-enum" translations.
2014-05-22 22:09:04 -05:00
Add Petsc to OPM.
2014-05-16 01:41:51 -05:00
LinearSolverReport rep = {};
rep.converged = true;
return rep;
}
void LinearSolverPetsc::setTolerance(const double /*tol*/)
{
}
double LinearSolverPetsc::getTolerance() const
{
return -1.;
}
} // namespace Opm
Reference in New Issue Copy Permalink