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.
This commit is contained in:
Jørgen Kvalsvik
@@ -51,7 +51,6 @@ list (APPEND MAIN_SOURCE_FILES
|
||||
opm/core/linalg/LinearSolverUmfpack.cpp
|
||||
opm/core/linalg/LinearSolverPetsc.cpp
|
||||
opm/core/linalg/call_umfpack.c
|
||||
opm/core/linalg/call_petsc.c
|
||||
opm/core/linalg/sparse_sys.c
|
||||
opm/core/pressure/CompressibleTpfa.cpp
|
||||
opm/core/pressure/FlowBCManager.cpp
|
||||
|
||||
@@ -19,9 +19,10 @@
|
||||
*/
|
||||
|
||||
#include "config.h"
|
||||
#include <cstring>
|
||||
#include <opm/core/linalg/LinearSolverPetsc.hpp>
|
||||
#include <opm/core/linalg/call_petsc.h>
|
||||
#include <unordered_map>
|
||||
#define PETSC_CLANGUAGE_CXX 1 //enable CHKERRXX macro.
|
||||
#include <petsc.h>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
|
||||
@@ -116,6 +117,106 @@ namespace{
|
||||
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 rhs;
|
||||
Vec solution;
|
||||
Mat A;
|
||||
KSP ksp;
|
||||
PC preconditioner;
|
||||
|
||||
OEM_DATA( const int size ) {
|
||||
CHKERRXX( VecCreate( PETSC_COMM_WORLD, &solution ) );
|
||||
CHKERRXX( VecSetSizes( solution, PETSC_DECIDE, size ) );
|
||||
CHKERRXX( VecSetFromOptions( solution ) );
|
||||
CHKERRXX( VecDuplicate( solution, &rhs ) );
|
||||
|
||||
CHKERRXX( MatCreate( PETSC_COMM_WORLD, &A ) );
|
||||
CHKERRXX( MatSetSizes( A, PETSC_DECIDE, PETSC_DECIDE, size, size ) );
|
||||
CHKERRXX( MatSetFromOptions( A ) );
|
||||
CHKERRXX( MatSetUp( A ) );
|
||||
}
|
||||
|
||||
~OEM_DATA() {
|
||||
CHKERRXX( VecDestroy( &rhs ) );
|
||||
CHKERRXX( VecDestroy( &solution ) );
|
||||
CHKERRXX( MatDestroy( &A ) );
|
||||
CHKERRXX( KSPDestroy( &ksp ) );
|
||||
}
|
||||
};
|
||||
|
||||
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;
|
||||
|
||||
CHKERRXX( VecGetLocalSize( v, &size ) );
|
||||
CHKERRXX( VecGetArray( v, &vec ) );
|
||||
|
||||
std::memcpy( vec, x, size * sizeof( double ) );
|
||||
CHKERRXX( VecRestoreArray( v, &vec ) );
|
||||
}
|
||||
|
||||
void from_petsc_vec( double* x, Vec v ) {
|
||||
if( !v ) OPM_THROW( std::runtime_error,
|
||||
"PETSc CopySolution: Invalid PETSc vector." );
|
||||
|
||||
PetscScalar* vec;
|
||||
PetscInt size;
|
||||
|
||||
CHKERRXX( VecGetLocalSize( v, &size ) );
|
||||
CHKERRXX( VecGetArray( v, &vec ) );
|
||||
|
||||
std::memcpy( x, vec, size * sizeof( double ) );
|
||||
CHKERRXX( VecRestoreArray( v, &vec ) );
|
||||
}
|
||||
|
||||
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;
|
||||
|
||||
for( int j = ia[ i ]; j < ia[ i + 1 ]; ++j )
|
||||
CHKERRXX( MatSetValues( A, 1, &i, 1, &ja[ j ], &sa[ j ], INSERT_VALUES ) );
|
||||
|
||||
CHKERRXX( MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY ) );
|
||||
CHKERRXX( MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY ) );
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
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;
|
||||
|
||||
CHKERRXX( KSPCreate( PETSC_COMM_WORLD, &t.ksp ) );
|
||||
CHKERRXX( KSPSetOperators( t.ksp, t.A, t.A, DIFFERENT_NONZERO_PATTERN ) );
|
||||
CHKERRXX( KSPGetPC( t.ksp, &t.preconditioner ) );
|
||||
CHKERRXX( KSPSetType( t.ksp, method ) );
|
||||
CHKERRXX( PCSetType( t.preconditioner, pcname ) );
|
||||
CHKERRXX( KSPSetTolerances( t.ksp, rtol, atol, dtol, maxits ) );
|
||||
CHKERRXX( KSPSetFromOptions( t.ksp ) );
|
||||
CHKERRXX( KSPSetInitialGuessNonzero( t.ksp, PETSC_TRUE ) );
|
||||
CHKERRXX( KSPSolve( t.ksp, t.rhs, t.solution ) );
|
||||
CHKERRXX( KSPGetConvergedReason( t.ksp, &reason ) );
|
||||
CHKERRXX( KSPGetIterationNumber( t.ksp, &its ) );
|
||||
CHKERRXX( KSPGetResidualNorm( t.ksp, &residual ) );
|
||||
|
||||
if( ksp_view )
|
||||
CHKERRXX( KSPView( t.ksp, PETSC_VIEWER_STDOUT_WORLD ) );
|
||||
|
||||
CHKERRXX( PetscPrintf( PETSC_COMM_WORLD, "KSP Iterations %D, Final Residual %G\n", its, residual ) );
|
||||
}
|
||||
|
||||
} // anonymous namespace.
|
||||
|
||||
LinearSolverPetsc::LinearSolverPetsc()
|
||||
@@ -123,7 +224,6 @@ namespace{
|
||||
OPM_THROW(std::runtime_error, "Pestc just can be called through paramers.\n");
|
||||
}
|
||||
|
||||
|
||||
LinearSolverPetsc::LinearSolverPetsc(const parameter::ParameterGroup& param)
|
||||
: ksp_type_("gmres")
|
||||
, pc_type_("sor")
|
||||
@@ -166,8 +266,13 @@ namespace{
|
||||
KSPType ksp_type = ksp.find(ksp_type_);
|
||||
PCTypeMap pc(pc_type_);
|
||||
PCType pc_type = pc.find(pc_type_);
|
||||
|
||||
call_Petsc(size, nonzeros, ia, ja, sa, rhs, solution, ksp_type, pc_type, rtol_, atol_, dtol_, maxits_, ksp_view_);
|
||||
|
||||
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 );
|
||||
|
||||
LinearSolverReport rep = {};
|
||||
rep.converged = true;
|
||||
@@ -183,6 +288,5 @@ namespace{
|
||||
return -1.;
|
||||
}
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
|
||||
@@ -1,244 +0,0 @@
|
||||
/*===========================================================================
|
||||
//
|
||||
// File: call_petsc.c
|
||||
//
|
||||
// Created: 2014-05-07 11:20:08 CST
|
||||
//
|
||||
// Authors: Ming Liu <miliu@statoil.com>
|
||||
//==========================================================================*/
|
||||
/*
|
||||
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 <assert.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <opm/core/linalg/call_petsc.h>
|
||||
|
||||
static PetscErrorCode code;
|
||||
|
||||
typedef struct {
|
||||
Vec b; /*b = rhs*/
|
||||
Vec u; /*u = solution*/
|
||||
Mat A; /*A = matrix*/
|
||||
KSP ksp; /*ksp = solver*/
|
||||
PC pc; /*pc = preconditioner*/
|
||||
} OEM_DATA;
|
||||
|
||||
typedef struct {
|
||||
KSPType method; /*ksp method*/
|
||||
PCType pcname; /*pc method*/
|
||||
int view_ksp; /*whether view ksp detail information*/
|
||||
double rtol;
|
||||
double atol;
|
||||
double dtol;
|
||||
int maxits;
|
||||
} KSP_OPT;
|
||||
|
||||
static int
|
||||
init(OEM_DATA* t, KSP_OPT* opts)
|
||||
{
|
||||
t->b = PETSC_NULL;
|
||||
t->u = PETSC_NULL;
|
||||
t->A = PETSC_NULL;
|
||||
t->ksp = PETSC_NULL;
|
||||
t->pc = PETSC_NULL;
|
||||
|
||||
/*set default options for ksp solvers*/
|
||||
opts->method = KSPGMRES;
|
||||
opts->pcname = PCJACOBI;
|
||||
opts->view_ksp = 0;
|
||||
opts->rtol = PETSC_DEFAULT;
|
||||
opts->atol = PETSC_DEFAULT;
|
||||
opts->dtol = PETSC_DEFAULT;
|
||||
opts->maxits = PETSC_DEFAULT;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
static int
|
||||
create_mat_vec(const int size, OEM_DATA* t)
|
||||
{
|
||||
code = VecCreate(PETSC_COMM_WORLD,&t->u);CHKERRQ(code);
|
||||
code = VecSetSizes(t->u,PETSC_DECIDE, size);CHKERRQ(code);
|
||||
code = VecSetFromOptions(t->u);CHKERRQ(code);
|
||||
code = VecDuplicate(t->u, &t->b);CHKERRQ(code);
|
||||
|
||||
code = MatCreate(PETSC_COMM_WORLD, &t->A);CHKERRQ(code);
|
||||
code = MatSetSizes(t->A,PETSC_DECIDE,PETSC_DECIDE, size, size);CHKERRQ(code);
|
||||
code = MatSetFromOptions(t->A);CHKERRQ(code);
|
||||
code = MatSetUp(t->A);CHKERRQ(code);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
static int
|
||||
to_petsc_vec(const double *x, Vec v)
|
||||
{
|
||||
PetscScalar *vec;
|
||||
PetscInt size;
|
||||
int i;
|
||||
|
||||
if (v == PETSC_NULL) {
|
||||
PetscPrintf(PETSC_COMM_WORLD,"PETSc CopySolution: invalid PETSc vector.\n");
|
||||
CHKERRQ(code);
|
||||
}
|
||||
code = VecGetLocalSize(v, &size); CHKERRQ(code);
|
||||
code = VecGetArray(v, &vec); CHKERRQ(code);
|
||||
for (i = 0; i < size; ++i) {
|
||||
vec[i] = x[i];
|
||||
}
|
||||
code = VecRestoreArray(v, &vec); CHKERRQ(code);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
static int
|
||||
from_petsc_vec(double *x, Vec v)
|
||||
{
|
||||
PetscScalar *vec;
|
||||
PetscInt size;
|
||||
int i;
|
||||
|
||||
if (v == PETSC_NULL) {
|
||||
PetscPrintf(PETSC_COMM_WORLD,"PETSc CopySolution: invalid PETSc vector.\n");
|
||||
CHKERRQ(code);
|
||||
}
|
||||
code = VecGetLocalSize(v, &size); CHKERRQ(code);
|
||||
code = VecGetArray(v, &vec); CHKERRQ(code);
|
||||
for (i = 0; i < size; ++i) {
|
||||
x[i] = vec[i];
|
||||
}
|
||||
code = VecRestoreArray(v, &vec); CHKERRQ(code);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
static int
|
||||
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;
|
||||
}
|
||||
for (int j = ia[i]; j < ia[i+1]; ++j) {
|
||||
code = MatSetValues(A,1,&i,1,&ja[j],&sa[j],INSERT_VALUES);
|
||||
CHKERRQ(code);
|
||||
}
|
||||
}
|
||||
code = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(code);
|
||||
code = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(code);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
static int
|
||||
solve_system(OEM_DATA* t, KSP_OPT* opts)
|
||||
{
|
||||
PetscInt its;
|
||||
PetscReal residual;
|
||||
KSPConvergedReason reason;
|
||||
code = KSPCreate(PETSC_COMM_WORLD,&t->ksp);CHKERRQ(code);
|
||||
code = KSPSetOperators(t->ksp,t->A,t->A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(code);
|
||||
code = KSPGetPC(t->ksp,&t->pc);CHKERRQ(code);
|
||||
code = KSPSetType(t->ksp, opts->method);
|
||||
code = PCSetType(t->pc, opts->pcname);CHKERRQ(code);
|
||||
code = KSPSetTolerances(t->ksp,opts->rtol,opts->atol,opts->dtol,opts->maxits);CHKERRQ(code);
|
||||
code = KSPSetFromOptions(t->ksp);CHKERRQ(code);
|
||||
code = KSPSetInitialGuessNonzero(t->ksp,PETSC_TRUE);CHKERRQ(code);
|
||||
code = KSPSolve(t->ksp,t->b,t->u);CHKERRQ(code);
|
||||
code = KSPGetConvergedReason(t->ksp, &reason); CHKERRQ(code);
|
||||
code = KSPGetIterationNumber(t->ksp, &its); CHKERRQ(code);
|
||||
code = KSPGetResidualNorm(t->ksp, &residual); CHKERRQ(code);
|
||||
if (opts->view_ksp) {
|
||||
code = KSPView(t->ksp,PETSC_VIEWER_STDOUT_WORLD);
|
||||
CHKERRQ(code);
|
||||
}
|
||||
|
||||
code = PetscPrintf(PETSC_COMM_WORLD,"KSP Iterations %D, Final Residual %G\n",its,residual);CHKERRQ(code);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
static int
|
||||
destory(OEM_DATA* t, KSP_OPT* opts)
|
||||
{
|
||||
if(t == NULL && opts == NULL) {
|
||||
return 0;
|
||||
}
|
||||
if(t->u != PETSC_NULL) {
|
||||
code = VecDestroy(&t->u);
|
||||
CHKERRQ(code);
|
||||
}
|
||||
if(t->b != PETSC_NULL) {
|
||||
code = VecDestroy(&t->b);
|
||||
CHKERRQ(code);
|
||||
}
|
||||
if(t->A != PETSC_NULL) {
|
||||
code = MatDestroy(&t->A);
|
||||
CHKERRQ(code);
|
||||
}
|
||||
if(t->ksp != PETSC_NULL) {
|
||||
code = KSPDestroy(&t->ksp);
|
||||
CHKERRQ(code);
|
||||
}
|
||||
|
||||
free(t);
|
||||
free(opts);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
static int
|
||||
set_ksp_opts(const KSPType ksp_type, const PCType pc_type, const double rtol, const double atol, const double dtol, const int maxits, const int view_ksp, KSP_OPT* opts)
|
||||
{
|
||||
opts->method = ksp_type;
|
||||
opts->pcname = pc_type;
|
||||
opts->view_ksp = view_ksp;
|
||||
opts->rtol = rtol;
|
||||
opts->atol = atol;
|
||||
opts->dtol = dtol;
|
||||
opts->maxits = maxits;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
int
|
||||
call_Petsc(const int size, const int nonzeros, const int* ia, const int* ja, const double* sa, const double* b, double* x, const KSPType ksp_type, const PCType pc_type, const double rtol, const double atol, const double dtol, const int maxits, const int view_ksp)
|
||||
{
|
||||
OEM_DATA* t;
|
||||
KSP_OPT* opts;
|
||||
t = malloc(sizeof(OEM_DATA));
|
||||
opts = malloc(sizeof(KSP_OPT));
|
||||
assert(t);
|
||||
assert(opts);
|
||||
init(t, opts);
|
||||
create_mat_vec(size, t);
|
||||
to_petsc_mat(size, nonzeros, ia, ja, sa, t->A);
|
||||
to_petsc_vec(b, t->b);
|
||||
set_ksp_opts(ksp_type, pc_type, rtol, atol, dtol, maxits, view_ksp, opts);
|
||||
solve_system(t, opts);
|
||||
from_petsc_vec(x, t->u);
|
||||
destory(t, opts);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
@@ -1,40 +0,0 @@
|
||||
/*===========================================================================
|
||||
//
|
||||
// File: call_petsc.h
|
||||
//
|
||||
// Created: 2014-05-07 10:21:21 CST
|
||||
//
|
||||
// Authors: Ming Liu <miliu@statoil.com>
|
||||
//==========================================================================*/
|
||||
/*
|
||||
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_CALL_PETSC_H_HEADER
|
||||
#define OPM_CALL_PETSC_H_HEADER
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
#include <petsc.h>
|
||||
int
|
||||
call_Petsc(const int size, const int nonzeros, const int* ia, const int* ja, const double* sa, const double* b, double* x, const KSPType ksp_type, const PCType pc_type, const double rtol, const double atol, const double dtol, const int maxits, const int view_ksp);
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
#endif /* OPM_CALL_PETSC_H_HEADER */
|
||||
Reference in New Issue
Block a user