From 34dfb85919d91c4d417266198fcdb898b074dc56 Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Fri, 16 May 2014 14:41:51 +0800 Subject: [PATCH] Add Petsc to OPM. --- opm/core/linalg/LinearSolverFactory.cpp | 6 + opm/core/linalg/LinearSolverPetsc.cpp | 110 ++++++++++ opm/core/linalg/LinearSolverPetsc.hpp | 96 +++++++++ opm/core/linalg/call_petsc.c | 263 ++++++++++++++++++++++++ opm/core/linalg/call_petsc.h | 39 ++++ 5 files changed, 514 insertions(+) create mode 100644 opm/core/linalg/LinearSolverPetsc.cpp create mode 100644 opm/core/linalg/LinearSolverPetsc.hpp create mode 100644 opm/core/linalg/call_petsc.c create mode 100644 opm/core/linalg/call_petsc.h diff --git a/opm/core/linalg/LinearSolverFactory.cpp b/opm/core/linalg/LinearSolverFactory.cpp index 207c4ea5..923feec6 100644 --- a/opm/core/linalg/LinearSolverFactory.cpp +++ b/opm/core/linalg/LinearSolverFactory.cpp @@ -31,6 +31,7 @@ #include #endif +#include #include #include #include @@ -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."); diff --git a/opm/core/linalg/LinearSolverPetsc.cpp b/opm/core/linalg/LinearSolverPetsc.cpp new file mode 100644 index 00000000..4e510854 --- /dev/null +++ b/opm/core/linalg/LinearSolverPetsc.cpp @@ -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 . +*/ + +#include "config.h" +#include +#include +#include +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 + diff --git a/opm/core/linalg/LinearSolverPetsc.hpp b/opm/core/linalg/LinearSolverPetsc.hpp new file mode 100644 index 00000000..61ad4e27 --- /dev/null +++ b/opm/core/linalg/LinearSolverPetsc.hpp @@ -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 . +*/ + +#ifndef OPM_LINEARSOLVERPETSC_HEADER_INCLUDED +#define OPM_LINEARSOLVERPETSC_HEADER_INCLUDED +#include +#include +#include + +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 diff --git a/opm/core/linalg/call_petsc.c b/opm/core/linalg/call_petsc.c new file mode 100644 index 00000000..9b48757f --- /dev/null +++ b/opm/core/linalg/call_petsc.c @@ -0,0 +1,263 @@ +/*=========================================================================== +// +// File: call_petsc.c +// +// Created: 2014-05-07 11:20:08 CST +// +// Authors: Ming Liu +//==========================================================================*/ +/* + 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 . +*/ + +#include "config.h" +#include +#include +#include +#include +#include + +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; +} + diff --git a/opm/core/linalg/call_petsc.h b/opm/core/linalg/call_petsc.h new file mode 100644 index 00000000..b69a49b7 --- /dev/null +++ b/opm/core/linalg/call_petsc.h @@ -0,0 +1,39 @@ +/*=========================================================================== +// +// File: call_petsc.h +// +// Created: 2014-05-07 10:21:21 CST +// +// Authors: Ming Liu +//==========================================================================*/ +/* + 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 . +*/ +#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 */