From 4491eb7a6ae7d9564898a87c8661b366c9957cdb Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Fri, 16 May 2014 14:41:51 +0800 Subject: [PATCH 01/19] Add Petsc to OPM. --- opm/core/linalg/LinearSolverFactory.cpp | 6 ++ opm/core/linalg/LinearSolverPetsc.cpp | 110 ++++++++++++++++++++++++ opm/core/linalg/LinearSolverPetsc.hpp | 96 +++++++++++++++++++++ 3 files changed, 212 insertions(+) create mode 100644 opm/core/linalg/LinearSolverPetsc.cpp create mode 100644 opm/core/linalg/LinearSolverPetsc.hpp diff --git a/opm/core/linalg/LinearSolverFactory.cpp b/opm/core/linalg/LinearSolverFactory.cpp index 207c4ea54..923feec69 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 000000000..4e5108541 --- /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 000000000..61ad4e27e --- /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 From ccb83a8503bc8690f43f4aad1e82ba787e852d29 Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Fri, 16 May 2014 14:53:49 +0800 Subject: [PATCH 02/19] Initialize Petsc from constructor. Remove private members for initializing Petsc. --- opm/core/linalg/LinearSolverPetsc.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/opm/core/linalg/LinearSolverPetsc.cpp b/opm/core/linalg/LinearSolverPetsc.cpp index 4e5108541..fcb9b1270 100644 --- a/opm/core/linalg/LinearSolverPetsc.cpp +++ b/opm/core/linalg/LinearSolverPetsc.cpp @@ -21,6 +21,7 @@ #include "config.h" #include #include +#include #include namespace Opm { @@ -32,7 +33,7 @@ namespace Opm LinearSolverPetsc::LinearSolverPetsc(const parameter::ParameterGroup& param) - , ksp_type_("gmres") + : ksp_type_("gmres") , pc_type_("sor") , view_ksp_(false) , rtol_(1e-5) @@ -90,7 +91,7 @@ namespace Opm break; } } - call_Petsc(size, nonzeros, ia, ja, sa, rhs, solution, argc_, argv_, ksp_type, pc_type, rtol_, atol_, dtol_, maxits_, view_ksp_); + call_Petsc(size, nonzeros, ia, ja, sa, rhs, solution, ksp_type, pc_type, rtol_, atol_, dtol_, maxits_, view_ksp_); LinearSolverReport rep = {}; rep.converged = true; return rep; From 559003ef1e7849866f5d69e0c723d597384a2d83 Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Fri, 16 May 2014 14:55:23 +0800 Subject: [PATCH 03/19] Test Petsc. --- tests/test_linearsolver.cpp | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/tests/test_linearsolver.cpp b/tests/test_linearsolver.cpp index 184b7765d..2143fe009 100644 --- a/tests/test_linearsolver.cpp +++ b/tests/test_linearsolver.cpp @@ -184,3 +184,14 @@ BOOST_AUTO_TEST_CASE(KAMGTest) } #endif #endif + +BOOST_AUTO_TEST_CASE(PETScTest) +{ + Opm::parameter::ParameterGroup param; + param.insertParameter(std::string("linsolver"), std::string("petsc")); + param.insertParameter(std::string("ksp_type"), std::string("cg")); + param.insertParameter(std::string("pc_type"), std::string("jacobi")); + param.insertParameter(std::string("ksp_rtol"), std::string("1e-10")); + param.insertParameter(std::string("ksp_view"), std::string("0")); + run_test(param); +} From 6f2d28d08abbfb7e496996509264bc7c4950e4e4 Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Fri, 16 May 2014 15:02:49 +0800 Subject: [PATCH 04/19] Throw information for Petsc. --- opm/core/linalg/LinearSolverFactory.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/linalg/LinearSolverFactory.cpp b/opm/core/linalg/LinearSolverFactory.cpp index 923feec69..034123a22 100644 --- a/opm/core/linalg/LinearSolverFactory.cpp +++ b/opm/core/linalg/LinearSolverFactory.cpp @@ -49,7 +49,7 @@ namespace Opm #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."); + OPM_THROW(std::runtime_error, "No linear solver available, you must have UMFPACK , dune-istl or Petsc installed to use LinearSolverFactory."); #endif } From 9d985f7421adf3567f1a99b95340c4486bc49177 Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Fri, 23 May 2014 11:09:04 +0800 Subject: [PATCH 05/19] use unordered_map for "string-enum" translations. --- opm/core/linalg/LinearSolverPetsc.cpp | 125 +++++++++++++++++++++----- 1 file changed, 101 insertions(+), 24 deletions(-) diff --git a/opm/core/linalg/LinearSolverPetsc.cpp b/opm/core/linalg/LinearSolverPetsc.cpp index fcb9b1270..c0eaa8896 100644 --- a/opm/core/linalg/LinearSolverPetsc.cpp +++ b/opm/core/linalg/LinearSolverPetsc.cpp @@ -21,10 +21,102 @@ #include "config.h" #include #include +#include #include #include + namespace Opm { + +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 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 Map; + + std::string default_type_; + Map type_map_; + }; + +} // anonymous namespace. LinearSolverPetsc::LinearSolverPetsc() { @@ -35,7 +127,7 @@ namespace Opm LinearSolverPetsc::LinearSolverPetsc(const parameter::ParameterGroup& param) : ksp_type_("gmres") , pc_type_("sor") - , view_ksp_(false) + , ksp_view_(false) , rtol_(1e-5) , atol_(1e-50) , dtol_(1e5) @@ -46,7 +138,7 @@ namespace Opm 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_))); + ksp_view_ = (param.getDefault("ksp_view", int(ksp_view_))); rtol_ = param.getDefault("ksp_rtol", rtol_); atol_ = param.getDefault("ksp_atol", atol_); dtol_ = param.getDefault("ksp_dtol", dtol_); @@ -60,7 +152,6 @@ namespace Opm } - LinearSolverInterface::LinearSolverReport LinearSolverPetsc::solve(const int size, const int nonzeros, @@ -71,27 +162,13 @@ namespace Opm 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, ksp_type, pc_type, rtol_, atol_, dtol_, maxits_, view_ksp_); + KSPTypeMap ksp(ksp_type_); + 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_); + LinearSolverReport rep = {}; rep.converged = true; return rep; From 8a118f45fe666c152264377af712576b20835047 Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Fri, 23 May 2014 11:10:22 +0800 Subject: [PATCH 06/19] rename the private member. --- opm/core/linalg/LinearSolverPetsc.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/linalg/LinearSolverPetsc.hpp b/opm/core/linalg/LinearSolverPetsc.hpp index 61ad4e27e..94dab44d6 100644 --- a/opm/core/linalg/LinearSolverPetsc.hpp +++ b/opm/core/linalg/LinearSolverPetsc.hpp @@ -81,7 +81,7 @@ namespace Opm private: std::string ksp_type_; std::string pc_type_; - int view_ksp_; + int ksp_view_; double rtol_; double atol_; double dtol_; From 28b354c6252d28a213a20669a86f84e0dd23c95a Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Tue, 8 Jul 2014 10:58:39 +0800 Subject: [PATCH 07/19] Add HAVE_PETSC macro. --- opm/core/linalg/LinearSolverFactory.cpp | 9 +++++++-- opm/core/linalg/LinearSolverFactory.hpp | 9 ++++++--- tests/test_linearsolver.cpp | 2 ++ 3 files changed, 15 insertions(+), 5 deletions(-) diff --git a/opm/core/linalg/LinearSolverFactory.cpp b/opm/core/linalg/LinearSolverFactory.cpp index 034123a22..0e22910b2 100644 --- a/opm/core/linalg/LinearSolverFactory.cpp +++ b/opm/core/linalg/LinearSolverFactory.cpp @@ -31,7 +31,10 @@ #include #endif +#if HAVE_PETSC #include +#endif + #include #include #include @@ -46,9 +49,9 @@ namespace Opm solver_.reset(new LinearSolverUmfpack); #elif HAVE_DUNE_ISTL solver_.reset(new LinearSolverIstl); -#else +#elif HAVE_PETSC solver_.reset(new LinearSolverPetsc); - +#else OPM_THROW(std::runtime_error, "No linear solver available, you must have UMFPACK , dune-istl or Petsc installed to use LinearSolverFactory."); #endif } @@ -73,7 +76,9 @@ namespace Opm #endif } else if (ls == "petsc"){ +#if HAVE_PETSC solver_.reset(new LinearSolverPetsc(param)); +#endif } else { diff --git a/opm/core/linalg/LinearSolverFactory.hpp b/opm/core/linalg/LinearSolverFactory.hpp index f27140bc0..92f70c337 100644 --- a/opm/core/linalg/LinearSolverFactory.hpp +++ b/opm/core/linalg/LinearSolverFactory.hpp @@ -42,16 +42,19 @@ namespace Opm /// Construct from parameters. /// The accepted parameters are (default) (allowed values): - /// linsolver ("umfpack") ("umfpack", "istl") + /// linsolver ("umfpack") ("umfpack", "istl", "petsc") /// For the umfpack solver to be available, this class must be /// compiled with UMFPACK support, as indicated by the /// variable HAVE_SUITESPARSE_UMFPACK_H in config.h. /// For the istl solver to be available, this class must be /// compiled with dune-istl support, as indicated by the /// variable HAVE_DUNE_ISTL in config.h. + /// For the petsc solver to be available, this class must be + /// compiled with petsc support, as indicated by the + /// variable HAVE_PETSC in config.h. /// Any further parameters are passed on to the constructors - /// of the actual solver used, see LinearSolverUmfpack - /// and LinearSolverIstl for details. + /// of the actual solver used, see LinearSolverUmfpack, + /// LinearSolverIstl and LinearSolverPetsc for details. LinearSolverFactory(const parameter::ParameterGroup& param); /// Destructor. diff --git a/tests/test_linearsolver.cpp b/tests/test_linearsolver.cpp index 2143fe009..6021218e6 100644 --- a/tests/test_linearsolver.cpp +++ b/tests/test_linearsolver.cpp @@ -185,6 +185,7 @@ BOOST_AUTO_TEST_CASE(KAMGTest) #endif #endif +#if HAVE_PETSC BOOST_AUTO_TEST_CASE(PETScTest) { Opm::parameter::ParameterGroup param; @@ -195,3 +196,4 @@ BOOST_AUTO_TEST_CASE(PETScTest) param.insertParameter(std::string("ksp_view"), std::string("0")); run_test(param); } +#endif From ade7aa658effae51ad2c26ec291868ccfac07bef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Kvalsvik?= Date: Mon, 22 Sep 2014 14:52:51 +0200 Subject: [PATCH 08/19] 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. --- opm/core/linalg/LinearSolverPetsc.cpp | 114 ++++++++++++++++++++++++-- 1 file changed, 109 insertions(+), 5 deletions(-) diff --git a/opm/core/linalg/LinearSolverPetsc.cpp b/opm/core/linalg/LinearSolverPetsc.cpp index c0eaa8896..9c8528994 100644 --- a/opm/core/linalg/LinearSolverPetsc.cpp +++ b/opm/core/linalg/LinearSolverPetsc.cpp @@ -19,9 +19,10 @@ */ #include "config.h" +#include #include -#include #include +#define PETSC_CLANGUAGE_CXX 1 //enable CHKERRXX macro. #include #include @@ -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 From ce3981a55edb8d4d82302bf2cd6c1f9f2f35083c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Kvalsvik?= Date: Mon, 22 Sep 2014 15:35:57 +0200 Subject: [PATCH 09/19] Makes using wrong constructor a compile-time error Petsc only supports initialisation through the ParameterGroup constructor. Calling the default, non-arg constructor is a static error, and not implementing it makes using it break compiles. --- opm/core/linalg/LinearSolverPetsc.cpp | 5 ----- opm/core/linalg/LinearSolverPetsc.hpp | 12 ++++-------- 2 files changed, 4 insertions(+), 13 deletions(-) diff --git a/opm/core/linalg/LinearSolverPetsc.cpp b/opm/core/linalg/LinearSolverPetsc.cpp index 9c8528994..68397b99c 100644 --- a/opm/core/linalg/LinearSolverPetsc.cpp +++ b/opm/core/linalg/LinearSolverPetsc.cpp @@ -219,11 +219,6 @@ namespace{ } // anonymous namespace. - 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") diff --git a/opm/core/linalg/LinearSolverPetsc.hpp b/opm/core/linalg/LinearSolverPetsc.hpp index 94dab44d6..9750ca3d7 100644 --- a/opm/core/linalg/LinearSolverPetsc.hpp +++ b/opm/core/linalg/LinearSolverPetsc.hpp @@ -33,14 +33,10 @@ namespace Opm { 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 + /// Declared, but not implemented. Petsc can only be created through + /// the ParameterGroup constructor, everything else is an error. This way + /// the error is caught compile time and not rune time, which is nice as + /// it is a static error. LinearSolverPetsc(); /// Construct from parameters From 93a84303394af39caf00818864e65933a37a45b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Kvalsvik?= Date: Mon, 22 Sep 2014 15:44:22 +0200 Subject: [PATCH 10/19] 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. --- opm/core/linalg/LinearSolverPetsc.cpp | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/opm/core/linalg/LinearSolverPetsc.cpp b/opm/core/linalg/LinearSolverPetsc.cpp index 68397b99c..8c7b3dd3c 100644 --- a/opm/core/linalg/LinearSolverPetsc.cpp +++ b/opm/core/linalg/LinearSolverPetsc.cpp @@ -220,24 +220,17 @@ namespace{ } // anonymous namespace. LinearSolverPetsc::LinearSolverPetsc(const parameter::ParameterGroup& param) - : ksp_type_("gmres") - , pc_type_("sor") - , ksp_view_(false) - , rtol_(1e-5) - , atol_(1e-50) - , dtol_(1e5) - , maxits_(1e5) + : 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"); - ksp_type_ = (param.getDefault("ksp_type", std::string(ksp_type_))); - pc_type_ = (param.getDefault("pc_type", std::string(pc_type_))); - ksp_view_ = (param.getDefault("ksp_view", int(ksp_view_))); - 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_); } From 7c78afa13d629d43e689012cdd35f370988263d4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Kvalsvik?= Date: Tue, 7 Oct 2014 15:50:36 +0200 Subject: [PATCH 11/19] 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. --- opm/core/linalg/LinearSolverPetsc.cpp | 82 +++++++++++++++------------ 1 file changed, 46 insertions(+), 36 deletions(-) diff --git a/opm/core/linalg/LinearSolverPetsc.cpp b/opm/core/linalg/LinearSolverPetsc.cpp index 8c7b3dd3c..e5fd3a614 100644 --- a/opm/core/linalg/LinearSolverPetsc.cpp +++ b/opm/core/linalg/LinearSolverPetsc.cpp @@ -128,22 +128,24 @@ namespace{ 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 ) ); + VecCreate( PETSC_COMM_WORLD, &solution ); + auto err = VecSetSizes( solution, PETSC_DECIDE, size ); + CHKERRXX( err ); + VecSetFromOptions( solution ); + VecDuplicate( solution, &rhs ); - CHKERRXX( MatCreate( PETSC_COMM_WORLD, &A ) ); - CHKERRXX( MatSetSizes( A, PETSC_DECIDE, PETSC_DECIDE, size, size ) ); - CHKERRXX( MatSetFromOptions( A ) ); - CHKERRXX( MatSetUp( A ) ); + MatCreate( PETSC_COMM_WORLD, &A ); + err = MatSetSizes( A, PETSC_DECIDE, PETSC_DECIDE, size, size ); + CHKERRXX( err ); + MatSetFromOptions( A ); + MatSetUp( A ); } ~OEM_DATA() { - CHKERRXX( VecDestroy( &rhs ) ); - CHKERRXX( VecDestroy( &solution ) ); - CHKERRXX( MatDestroy( &A ) ); - CHKERRXX( KSPDestroy( &ksp ) ); + VecDestroy( &rhs ); + VecDestroy( &solution ); + MatDestroy( &A ); + KSPDestroy( &ksp ); } }; @@ -154,11 +156,11 @@ namespace{ PetscScalar* vec; PetscInt size; - CHKERRXX( VecGetLocalSize( v, &size ) ); - CHKERRXX( VecGetArray( v, &vec ) ); + VecGetLocalSize( v, &size ); + VecGetArray( v, &vec ); std::memcpy( vec, x, size * sizeof( double ) ); - CHKERRXX( VecRestoreArray( v, &vec ) ); + VecRestoreArray( v, &vec ); } void from_petsc_vec( double* x, Vec v ) { @@ -168,11 +170,11 @@ namespace{ PetscScalar* vec; PetscInt size; - CHKERRXX( VecGetLocalSize( v, &size ) ); - CHKERRXX( VecGetArray( v, &vec ) ); + VecGetLocalSize( v, &size ); + VecGetArray( v, &vec ); std::memcpy( x, vec, size * sizeof( double ) ); - CHKERRXX( VecRestoreArray( v, &vec ) ); + VecRestoreArray( v, &vec ); } void to_petsc_mat( const int size, const int nonzeros, @@ -183,11 +185,14 @@ namespace{ 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 ) ); + 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 ); + } - CHKERRXX( MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY ) ); - CHKERRXX( MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY ) ); + auto err = MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY ); + CHKERRXX( err ); + MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY ); } } @@ -198,23 +203,28 @@ namespace{ 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 ) ); + KSPCreate( PETSC_COMM_WORLD, &t.ksp ); + 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 ); if( ksp_view ) - CHKERRXX( KSPView( t.ksp, PETSC_VIEWER_STDOUT_WORLD ) ); + KSPView( t.ksp, PETSC_VIEWER_STDOUT_WORLD ); - CHKERRXX( PetscPrintf( PETSC_COMM_WORLD, "KSP Iterations %D, Final Residual %G\n", its, residual ) ); + err = PetscPrintf( PETSC_COMM_WORLD, "KSP Iterations %D, Final Residual %G\n", its, residual ); + CHKERRXX( err ); } } // anonymous namespace. From e69d92cca35639446287de932dc3e155fadf5639 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Kvalsvik?= Date: Wed, 12 Nov 2014 20:40:46 +0100 Subject: [PATCH 12/19] Moved KSPCreate to the Data constructor This is where it was always intended to be called, so this fixes a mistake in earlier development. --- opm/core/linalg/LinearSolverPetsc.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opm/core/linalg/LinearSolverPetsc.cpp b/opm/core/linalg/LinearSolverPetsc.cpp index e5fd3a614..1adf55986 100644 --- a/opm/core/linalg/LinearSolverPetsc.cpp +++ b/opm/core/linalg/LinearSolverPetsc.cpp @@ -139,6 +139,8 @@ namespace{ CHKERRXX( err ); MatSetFromOptions( A ); MatSetUp( A ); + + KSPCreate( PETSC_COMM_WORLD, &ksp ); } ~OEM_DATA() { @@ -203,7 +205,6 @@ namespace{ PetscReal residual; KSPConvergedReason reason; - KSPCreate( PETSC_COMM_WORLD, &t.ksp ); KSPSetOperators( t.ksp, t.A, t.A, DIFFERENT_NONZERO_PATTERN ); KSPGetPC( t.ksp, &t.preconditioner ); auto err = KSPSetType( t.ksp, method ); From 43d82a232cda4b511616bf5f98705a402958266e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Kvalsvik?= Date: Wed, 10 Dec 2014 18:05:04 +0100 Subject: [PATCH 13/19] Fixes performance bug in matrix construction. At the same time changes the construction to be value-oriented, not constructing it in OEM_DATA, but rather in to_petsc_mat --- opm/core/linalg/LinearSolverPetsc.cpp | 30 +++++++-------------------- 1 file changed, 8 insertions(+), 22 deletions(-) diff --git a/opm/core/linalg/LinearSolverPetsc.cpp b/opm/core/linalg/LinearSolverPetsc.cpp index 1adf55986..a0fb5771f 100644 --- a/opm/core/linalg/LinearSolverPetsc.cpp +++ b/opm/core/linalg/LinearSolverPetsc.cpp @@ -134,11 +134,7 @@ namespace{ 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 ); + VecSetFromOptions( b ); KSPCreate( PETSC_COMM_WORLD, &ksp ); } @@ -179,23 +175,13 @@ namespace{ VecRestoreArray( v, &vec ); } - void to_petsc_mat( const int size, const int nonzeros, - const int* ia, const int* ja, const double* sa, Mat A ) { + Mat to_petsc_mat( const int size, const int nonzeros, + const int* ia, const int* ja, const double* sa ) { - 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 ) { - auto err = MatSetValues( A, 1, &i, 1, &ja[ j ], &sa[ j ], INSERT_VALUES ); - CHKERRXX( err ); - } - - auto err = MatAssemblyBegin( A, MAT_FINAL_ASSEMBLY ); - CHKERRXX( err ); - MatAssemblyEnd( A, MAT_FINAL_ASSEMBLY ); - } + Mat A; + auto err = MatCreateSeqAIJWithArrays( PETSC_COMM_WORLD, size, size, (int*)ia, (int*)ja, (double*)sa, &A ); + CHKERRXX( err ); + return A; } @@ -267,7 +253,7 @@ namespace{ PCType pc_type = pc.find(pc_type_); OEM_DATA t( size ); - to_petsc_mat( size, nonzeros, ia, ja, sa, t.A ); + t.A = to_petsc_mat( size, nonzeros, ia, ja, sa ); to_petsc_vec( rhs, t.rhs ); solve_system( t, ksp_type, pc_type, rtol_, atol_, dtol_, maxits_, ksp_view_ ); From 6c70d2938dc33acac3227fc9719122338dfd8f40 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Kvalsvik?= Date: Wed, 10 Dec 2014 18:06:28 +0100 Subject: [PATCH 14/19] Removes unecessary whitespace. --- opm/core/linalg/LinearSolverPetsc.cpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/opm/core/linalg/LinearSolverPetsc.cpp b/opm/core/linalg/LinearSolverPetsc.cpp index a0fb5771f..4d6528831 100644 --- a/opm/core/linalg/LinearSolverPetsc.cpp +++ b/opm/core/linalg/LinearSolverPetsc.cpp @@ -28,7 +28,7 @@ namespace Opm { - + namespace{ class KSPTypeMap { @@ -52,8 +52,8 @@ namespace{ type_map_.insert(std::make_pair("cr", KSPCR)); type_map_.insert(std::make_pair("preonly", KSPPREONLY)); } - - KSPType + + KSPType find(const std::string& type) const { Map::const_iterator it = type_map_.find(type); @@ -61,7 +61,7 @@ namespace{ if (it == type_map_.end()) { it = type_map_.find(default_type_); } - + if (it == type_map_.end()) { OPM_THROW(std::runtime_error, "Unknown KSPType: '" << type << "'"); } @@ -95,8 +95,8 @@ namespace{ type_map_.insert(std::make_pair("cholesky", PCCHOLESKY)); type_map_.insert(std::make_pair("none", PCNONE)); } - - PCType + + PCType find(const std::string& type) const { Map::const_iterator it = type_map_.find(type); @@ -104,10 +104,11 @@ namespace{ 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: From 37a47711e52e678868e3815d32da016a64949178 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Kvalsvik?= Date: Wed, 10 Dec 2014 18:07:45 +0100 Subject: [PATCH 15/19] Vector construction is now value oriented. No longer initialised and constructed in two places. --- opm/core/linalg/LinearSolverPetsc.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/opm/core/linalg/LinearSolverPetsc.cpp b/opm/core/linalg/LinearSolverPetsc.cpp index 4d6528831..58e659bd5 100644 --- a/opm/core/linalg/LinearSolverPetsc.cpp +++ b/opm/core/linalg/LinearSolverPetsc.cpp @@ -148,18 +148,18 @@ namespace{ } }; - void to_petsc_vec( const double* x, Vec v ) { - if( !v ) OPM_THROW( std::runtime_error, - "PETSc CopySolution: Invalid PETSc vector." ); - + Vec to_petsc_vec( const double* x, int size ) { PetscScalar* vec; - PetscInt size; + Vec v; + + VecSetSizes( v, PETSC_DECIDE, size ); + VecSetFromOptions( v ); - VecGetLocalSize( v, &size ); VecGetArray( v, &vec ); - std::memcpy( vec, x, size * sizeof( double ) ); + VecRestoreArray( v, &vec ); + return v; } void from_petsc_vec( double* x, Vec v ) { @@ -255,7 +255,7 @@ namespace{ OEM_DATA t( size ); t.A = to_petsc_mat( size, nonzeros, ia, ja, sa ); - to_petsc_vec( rhs, t.rhs ); + t.rhs = to_petsc_vec( rhs, size ); solve_system( t, ksp_type, pc_type, rtol_, atol_, dtol_, maxits_, ksp_view_ ); from_petsc_vec( solution, t.solution ); From 8fbbdad722a24e16ca9da8c94f343f296e3295ff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Kvalsvik?= Date: Wed, 10 Dec 2014 18:09:01 +0100 Subject: [PATCH 16/19] Renamed variables to match Ax = b --- opm/core/linalg/LinearSolverPetsc.cpp | 21 +++++++++------------ 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/opm/core/linalg/LinearSolverPetsc.cpp b/opm/core/linalg/LinearSolverPetsc.cpp index 58e659bd5..202db19b1 100644 --- a/opm/core/linalg/LinearSolverPetsc.cpp +++ b/opm/core/linalg/LinearSolverPetsc.cpp @@ -122,27 +122,24 @@ namespace{ /* 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; + Vec x; + Vec b; Mat A; KSP ksp; PC preconditioner; OEM_DATA( const int size ) { - VecCreate( PETSC_COMM_WORLD, &solution ); - auto err = VecSetSizes( solution, PETSC_DECIDE, size ); + VecCreate( PETSC_COMM_WORLD, &b ); + auto err = VecSetSizes( b, PETSC_DECIDE, size ); CHKERRXX( err ); - VecSetFromOptions( solution ); - VecDuplicate( solution, &rhs ); - VecSetFromOptions( b ); KSPCreate( PETSC_COMM_WORLD, &ksp ); } ~OEM_DATA() { - VecDestroy( &rhs ); - VecDestroy( &solution ); + VecDestroy( &x ); + VecDestroy( &b ); MatDestroy( &A ); KSPDestroy( &ksp ); } @@ -203,7 +200,7 @@ namespace{ err = KSPSetFromOptions( t.ksp ); CHKERRXX( err ); KSPSetInitialGuessNonzero( t.ksp, PETSC_TRUE ); - KSPSolve( t.ksp, t.rhs, t.solution ); + KSPSolve( t.ksp, t.x, t.b ); KSPGetConvergedReason( t.ksp, &reason ); KSPGetIterationNumber( t.ksp, &its ); KSPGetResidualNorm( t.ksp, &residual ); @@ -255,10 +252,10 @@ namespace{ OEM_DATA t( size ); t.A = to_petsc_mat( size, nonzeros, ia, ja, sa ); - t.rhs = to_petsc_vec( rhs, size ); + 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.solution ); + from_petsc_vec( solution, t.b ); LinearSolverReport rep = {}; rep.converged = true; From b3ebe04e141e13f60ed23f508b5337b2faf248a7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Kvalsvik?= Date: Wed, 10 Dec 2014 18:09:13 +0100 Subject: [PATCH 17/19] Adds missing preconditoner destruction. --- opm/core/linalg/LinearSolverPetsc.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/opm/core/linalg/LinearSolverPetsc.cpp b/opm/core/linalg/LinearSolverPetsc.cpp index 202db19b1..daaa71447 100644 --- a/opm/core/linalg/LinearSolverPetsc.cpp +++ b/opm/core/linalg/LinearSolverPetsc.cpp @@ -142,6 +142,7 @@ namespace{ VecDestroy( &b ); MatDestroy( &A ); KSPDestroy( &ksp ); + PCDestroy( &preconditioner ); } }; From 294ec58bed3f3e0439a5d5807462ab8504a6a60f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Kvalsvik?= Date: Wed, 10 Dec 2014 18:15:13 +0100 Subject: [PATCH 18/19] Adds missing VecCreate. --- opm/core/linalg/LinearSolverPetsc.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/opm/core/linalg/LinearSolverPetsc.cpp b/opm/core/linalg/LinearSolverPetsc.cpp index daaa71447..a6e7f5ee5 100644 --- a/opm/core/linalg/LinearSolverPetsc.cpp +++ b/opm/core/linalg/LinearSolverPetsc.cpp @@ -150,6 +150,7 @@ namespace{ PetscScalar* vec; Vec v; + VecCreate( PETSC_COMM_WORLD, &v ); VecSetSizes( v, PETSC_DECIDE, size ); VecSetFromOptions( v ); From 263e1a692835a5fa97b6e939bf0ca0b842e6f6b4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Kvalsvik?= Date: Wed, 10 Dec 2014 18:15:33 +0100 Subject: [PATCH 19/19] Sets initial guess to zero. Setting it to nonzero means starting to approximate from memory garbage values. --- opm/core/linalg/LinearSolverPetsc.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/linalg/LinearSolverPetsc.cpp b/opm/core/linalg/LinearSolverPetsc.cpp index a6e7f5ee5..1b4fe4fb1 100644 --- a/opm/core/linalg/LinearSolverPetsc.cpp +++ b/opm/core/linalg/LinearSolverPetsc.cpp @@ -201,7 +201,7 @@ namespace{ CHKERRXX( err ); err = KSPSetFromOptions( t.ksp ); CHKERRXX( err ); - KSPSetInitialGuessNonzero( t.ksp, PETSC_TRUE ); + KSPSetInitialGuessNonzero( t.ksp, PETSC_FALSE ); KSPSolve( t.ksp, t.x, t.b ); KSPGetConvergedReason( t.ksp, &reason ); KSPGetIterationNumber( t.ksp, &its );