Add Petsc to OPM.

This commit is contained in:
Liu Ming 2014-05-16 14:41:51 +08:00
parent d50b7ab4aa
commit 34dfb85919
5 changed files with 514 additions and 0 deletions

View File

@ -31,6 +31,7 @@
#include <opm/core/linalg/LinearSolverIstl.hpp>
#endif
#include <opm/core/linalg/LinearSolverPetsc.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <string>
@ -46,6 +47,8 @@ namespace Opm
#elif HAVE_DUNE_ISTL
solver_.reset(new LinearSolverIstl);
#else
solver_.reset(new LinearSolverPetsc);
OPM_THROW(std::runtime_error, "No linear solver available, you must have UMFPACK or dune-istl installed to use LinearSolverFactory.");
#endif
}
@ -69,6 +72,9 @@ namespace Opm
solver_.reset(new LinearSolverIstl(param));
#endif
}
else if (ls == "petsc"){
solver_.reset(new LinearSolverPetsc(param));
}
else {
OPM_THROW(std::runtime_error, "Linear solver " << ls << " is unknown.");

View File

@ -0,0 +1,110 @@
/*
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 <opm/core/linalg/LinearSolverPetsc.hpp>
#include <opm/core/linalg/call_petsc.h>
#include <opm/core/utility/ErrorMacros.hpp>
namespace Opm
{
LinearSolverPetsc::LinearSolverPetsc()
{
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")
, view_ksp_(false)
, rtol_(1e-5)
, atol_(1e-50)
, dtol_(1e5)
, maxits_(1e5)
{
int argc = 0;
char** argv = NULL;
PetscInitialize(&argc, &argv, (char*)0, "Petsc interface for OPM!\n");
ksp_type_ = (param.getDefault("ksp_type", std::string(ksp_type_)));
pc_type_ = (param.getDefault("pc_type", std::string(pc_type_)));
view_ksp_ = (param.getDefault("ksp_view", int(view_ksp_)));
rtol_ = param.getDefault("ksp_rtol", rtol_);
atol_ = param.getDefault("ksp_atol", atol_);
dtol_ = param.getDefault("ksp_dtol", dtol_);
maxits_ = param.getDefault("ksp_max_it", maxits_);
}
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
{
int ksp_type=4;
int pc_type=0;
const std::string ksp[]={"richardson", "chebyshev", "cg", "bicg",
"gmres", "fgmres", "dgmres", "gcr", "bcgs",
"cgs", "tfqmr", "tcqmr", "cr", "lsqr", "preonly"};
const std::string pc[] = {"jacobi", "bjacobi", "sor", "eisenstat",
"icc", "ilu", "pcasm", "gamg", "ksp",
"composit", "lu", "cholesky", "none"};
for (int i = 0; i < 15; ++i){
if (ksp[i] == ksp_type_){
ksp_type=i;
break;
}
}
for (int i = 0; i < 13; ++i){
if (pc[i] == pc_type_){
pc_type=i;
break;
}
}
call_Petsc(size, nonzeros, ia, ja, sa, rhs, solution, argc_, argv_, ksp_type, pc_type, rtol_, atol_, dtol_, maxits_, view_ksp_);
LinearSolverReport rep = {};
rep.converged = true;
return rep;
}
void LinearSolverPetsc::setTolerance(const double /*tol*/)
{
}
double LinearSolverPetsc::getTolerance() const
{
return -1.;
}
} // namespace Opm

View File

@ -0,0 +1,96 @@
/*
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.
/// All parameters controlling the solver are defaulted:
/// ksp_type gmres
/// pc_type jacobi
/// ksp_view 0
/// ksp_rtol 1e-5
/// ksp_atol 1e-50
/// ksp_dtol 1e5
/// ksp_max_it 1e5
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 view_ksp_;
double rtol_;
double atol_;
double dtol_;
int maxits_;
};
} // namespace Opm
#endif // OPM_LINEARSOLVERPETSC_HEADER_INCLUDED

View File

@ -0,0 +1,263 @@
/*===========================================================================
//
// 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 <petsc.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; /*weather view ksp detail information*/
double rtol;
double atol;
double dtol;
int maxits;
} KSP_OPT;
static const KSPType ksp_list[] = {
KSPRICHARDSON, /*richardson*/
KSPCHEBYSHEV, /*chebyshev*/
KSPCG, /*cg*/
KSPBICG, /*bicgs*/
KSPGMRES, /*gmres*/
KSPFGMRES, /*fgmres*/
KSPDGMRES, /*dgmres*/
KSPGCR, /*gcr*/
KSPBCGS, /*bcgs*/
KSPCGS, /*cgs*/
KSPTFQMR, /*tfqmr*/
KSPTCQMR, /*tcqmr*/
KSPCR, /*cr*/
KSPLSQR, /*lsqr*/
KSPPREONLY /*preonly*/
};
static const PCType pc_list[] = {
PCJACOBI, /*jacobi*/
PCBJACOBI, /*bjacobi*/
PCSOR, /*sor*/
PCEISENSTAT, /*eisenstat*/
PCICC, /*icc*/
PCILU, /*ilu*/
PCASM, /*asm*/
PCGAMG, /*gamg*/
PCKSP, /*ksp*/
PCCOMPOSITE, /*composite*/
PCLU, /*lu*/
PCCHOLESKY, /*cholesky*/
PCNONE /*none*/
};
static int
init(OEM_DATA* t, KSP_OPT* opts)
{
if (t == NULL)
t = calloc(1, sizeof(OEM_DATA));
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 = PCSOR;
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(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)
printf("PETSc CopySolution: invalid PETSc vector.\n");
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)
printf("PETSc CopySolution: invalid PETSc vector.\n");
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(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)
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);
if (opts == NULL)
return 0;
free(opts);
return 0;
}
static int
set_ksp_opts(const int ksp_type, const int pc_type, const double rtol, const double atol, const double dtol, const int maxits, const int view_ksp, KSP_OPT* opts)
{
if(opts == NULL) {
opts = calloc(1, sizeof(KSP_OPT));
}
opts->method = ksp_list[ksp_type];
opts->pcname = pc_list[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, int argc, char** argv, const int ksp_type, const int 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 = calloc(1, sizeof(OEM_DATA));
opts = calloc(1, sizeof(KSP_OPT));
init(t, opts);
create(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(t, opts);
from_petsc_vec(x, t->u);
destory(t, opts);
return 0;
}

View File

@ -0,0 +1,39 @@
/*===========================================================================
//
// 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
int
call_Petsc(const int size, const int nonzeros, const int* ia, const int* ja, const double* sa, const double* b, double* x, int argc, char** argv, const int ksp_type, const int 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 */