OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #702 from jorgekva/linsolver_petsc

Browse Source 
PETSc support in OPM
This commit is contained in:
Atgeirr Flø Rasmussen 2014-12-19 11:25:58 +01:00
commit 060a2fea0f
9 changed files with 563 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
CMakeLists.txt
View File

@ -71,6 +71,12 @@ macro (sources_hook)
)
endif (NOT SuiteSparse_FOUND)
if (NOT PETSC_FOUND)
list (REMOVE_ITEM opm-core_SOURCES
${PROJECT_SOURCE_DIR}/${opm-core_DIR}/core/linalg/call_petsc.c
${PROJECT_SOURCE_DIR}/${opm-core_DIR}/core/linalg/LinearSolverPetsc.cpp
)
endif (NOT PETSC_FOUND)
if ((NOT MPI_FOUND) OR (NOT DUNE_ISTL_FOUND))
list (REMOVE_ITEM tests_SOURCES
${PROJECT_SOURCE_DIR}/tests/test_parallel_linearsolver.cpp

2
CMakeLists_files.cmake
View File

@ -49,6 +49,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/linalg/LinearSolverInterface.cpp
opm/core/linalg/LinearSolverIstl.cpp
opm/core/linalg/LinearSolverUmfpack.cpp
opm/core/linalg/LinearSolverPetsc.cpp
opm/core/linalg/call_umfpack.c
opm/core/linalg/sparse_sys.c
opm/core/pressure/CompressibleTpfa.cpp
@ -286,6 +287,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/linalg/LinearSolverInterface.hpp
opm/core/linalg/LinearSolverIstl.hpp
opm/core/linalg/LinearSolverUmfpack.hpp
opm/core/linalg/LinearSolverPetsc.hpp
opm/core/linalg/blas_lapack.h
opm/core/linalg/call_umfpack.h
opm/core/linalg/sparse_sys.h

153
cmake/Modules/FindPetsc.cmake Normal file
View File

@ -0,0 +1,153 @@
# - Try to find Petsc lib
#
# This module supports requiring a minimum version, e.g. you can do
# find_package(Petsc)
#
# Once done this will define
#
# PETSC_FOUND - system has Petsc lib with correct version
# PETSC_INCLUDE_DIRS - the Petsc include directory
# PETSC_LIBRARIES - the Petsc library.
# Copyright (c) 2006, 2007 Montel Laurent, <montel@kde.org>
# Copyright (c) 2008, 2009 Gael Guennebaud, <g.gael@free.fr>
# Copyright (c) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
# Redistribution and use is allowed according to the terms of the 2-clause BSD license.
# find out the size of a pointer. this is required to only search for
# libraries in the directories relevant for the architecture
if (CMAKE_SIZEOF_VOID_P)
math (EXPR _BITS "8 * ${CMAKE_SIZEOF_VOID_P}")
endif (CMAKE_SIZEOF_VOID_P)
# if PETSC_ROOT is set, then this is the only place searched for petsc headers
# and includes
set(_no_default_path "")
if(PETSC_ROOT)
set (_no_default_path "NO_DEFAULT_PATH")
endif()
# look for a system-wide BLAS library
set(PETSC_BLAS_LIBRARY "")
find_package(BLAS QUIET)
list(APPEND PETSC_BLAS_LIBRARY "${BLAS_LIBRARIES}")
# if BLAS wasn't found, look for it in PETSC_ROOT. Won't search if
# PETSC_BLAS_LIBRARY is set.
find_library(PETSC_BLAS_LIBRARY
NAME "blas"
PATH ${PETSC_ROOT}
PATH_SUFFIXES "lib" "lib${_BITS}" "lib/${CMAKE_LIBRARY_ARCHITECTURE}"
${_no_default_path}
)
# print message if there was still no blas found!
if(NOT BLAS_FOUND AND NOT PETSC_BLAS_LIBRARY)
message(STATUS "BLAS not found but required for PETSC")
return()
endif()
list(APPEND CMAKE_REQUIRED_LIBRARIES "${PETSC_BLAS_LIBRARY}")
set(PETSC_LAPACK_LIBRARY "")
find_package(LAPACK QUIET)
list(APPEND PETSC_LAPACK_LIBRARY "${LAPACK_LIBRARIES}")
# if LAPACK wasn't found, look for it in PETSC_ROOT
find_library(PETSC_LAPACK_LIBRARY
NAME "lapack"
PATH ${PETSC_ROOT}
PATH_SUFFIXES "lib" "lib${_BITS}" "lib/${CMAKE_LIBRARY_ARCHITECTURE}"
${_no_default_path}
)
# print message if there was still no blas found!
if(NOT LAPACK_FOUND AND NOT PETSC_LAPACK_LIBRARY)
message(STATUS "LAPACK not found but required for PETSC")
return()
endif()
list(APPEND CMAKE_REQUIRED_LIBRARIES "${PETSC_LAPACK_LIBRARY}")
find_package(X11 QUIET)
if (X11_FOUND)
list(APPEND PETSC_X11_LIBRARY "${X11_LIBRARIES}")
endif()
list(APPEND CMAKE_REQUIRED_LIBRARIES "${PETSC_X11_LIBRARY}")
# only probe if we haven't a path in our cache
if (Petsc_ROOT)
set (PETSC_ROOT "${Petsc_ROOT}")
endif (Petsc_ROOT)
find_path (PETSC_NORMAL_INCLUDE_DIR
NAMES "petsc.h"
PATHS ${PETSC_ROOT}
PATH_SUFFIXES "petsc-3.4.4" "include" "petsc"
${_no_default_path}
)
# if parallel computing is explicitly enabled - reuse the paths and links from
# OpmMainLib + OpmFind
# this needs to be called explicitly as FindPetsc runs before OpmMainLib starts
# looking for MPI (for some reason). Ideally this isn't necessary
if(USE_MPI)
find_package(MPI)
endif()
set(PETSC_MPI_FOUND ${MPI_FOUND})
set(PETSC_MPI_INCLUDE_DIRS ${MPI_INCLUDE_PATH})
set(PETSC_MPI_LIBRARIES ${MPI_LIBRARIES})
# fallback - use the petsc provided implementation of serial MPI.
if (NOT PETSC_MPI_INCLUDE_DIRS)
message(STATUS "Building without MPI support - looking for PETSc provided serial implementation")
find_path (PETSC_MPI_INCLUDE_DIRS
NAMES "mpi.h"
PATHS ${PETSC_ROOT}/include
PATH_SUFFIXES "mpiuni"
${_no_default_path}
)
if(PETSC_MPI_INCLUDE_DIRS)
set(PETSC_MPI_FOUND 1)
endif(PETSC_MPI_INCLUDE_DIRS)
endif(NOT PETSC_MPI_INCLUDE_DIRS)
# couldn't find any usable mpi implementation - abort
if(NOT PETSC_MPI_FOUND)
message(STATUS "Could not find any suitable MPI implementation.")
return()
endif()
# look for actual Petsc library
find_library(PETSC_LIBRARY
NAMES "petsc-3.4.3" "petsc-3.4.4" "petsc"
PATHS ${PETSC_ROOT}
PATH_SUFFIXES "lib" "lib${_BITS}" "lib/${CMAKE_LIBRARY_ARCHITECTURE}"
${_no_default_path}
)
if(NOT PETSC_LIBRARY)
message(STATUS "Could not find the PETSc library")
return()
endif()
list(APPEND CMAKE_REQUIRED_LIBRARIES "${PETSC_LIBRARY}" "${PETSC_MPI_LIBRARIES}")
if (PETSC_MPI_INCLUDE_DIRS AND PETSC_NORMAL_INCLUDE_DIR)
list (APPEND PETSC_INCLUDE_DIR ${PETSC_MPI_INCLUDE_DIRS} ${PETSC_NORMAL_INCLUDE_DIR})
endif()
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(Petsc DEFAULT_MSG PETSC_INCLUDE_DIR PETSC_LIBRARY)
mark_as_advanced(PETSC_INCLUDE_DIR PETSC_LIBRARY)
# if both headers and library are found, store results
if(PETSC_FOUND)
set(PETSC_INCLUDE_DIRS ${PETSC_INCLUDE_DIR})
set(PETSC_LIBRARIES ${PETSC_LIBRARY})
list(APPEND PETSC_LIBRARIES ${PETSC_BLAS_LIBRARY})
list(APPEND PETSC_LIBRARIES ${PETSC_LAPACK_LIBRARY})
list(APPEND PETSC_LIBRARIES ${PETSC_X11_LIBRARY})
endif()

2
cmake/Modules/opm-core-prereqs.cmake
View File

@ -7,6 +7,7 @@ set (opm-core_CONFIG_VAR
HAVE_SUITESPARSE_UMFPACK_H
HAVE_DUNE_ISTL
HAVE_MPI
HAVE_PETSC
)
# dependencies
@ -25,6 +26,7 @@ set (opm-core_DEPS
"SuiteSparse COMPONENTS umfpack"
# solver
"SuperLU"
"Petsc"
# xml processing (for config parsing)
"TinyXML"
# Ensembles-based Reservoir Tools (ERT)

13
opm/core/linalg/LinearSolverFactory.cpp
View File

@ -31,6 +31,10 @@
#include <opm/core/linalg/LinearSolverIstl.hpp>
#endif
#if HAVE_PETSC
#include <opm/core/linalg/LinearSolverPetsc.hpp>
#endif
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <string>
@ -45,8 +49,10 @@ namespace Opm
solver_.reset(new LinearSolverUmfpack);
#elif HAVE_DUNE_ISTL
solver_.reset(new LinearSolverIstl);
#elif HAVE_PETSC
solver_.reset(new LinearSolverPetsc);
#else
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
}
@ -69,6 +75,11 @@ namespace Opm
solver_.reset(new LinearSolverIstl(param));
#endif
}
else if (ls == "petsc"){
#if HAVE_PETSC
solver_.reset(new LinearSolverPetsc(param));
#endif
}
else {
OPM_THROW(std::runtime_error, "Linear solver " << ls << " is unknown.");

9
opm/core/linalg/LinearSolverFactory.hpp
View File

@ -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.

277
opm/core/linalg/LinearSolverPetsc.cpp Normal file
View File

@ -0,0 +1,277 @@
/*
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 <cstring>
#include <opm/core/linalg/LinearSolverPetsc.hpp>
#include <unordered_map>
#define PETSC_CLANGUAGE_CXX 1 //enable CHKERRXX macro.
#include <petsc.h>
#include <opm/core/utility/ErrorMacros.hpp>
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<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_;
};
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 x;
Vec b;
Mat A;
KSP ksp;
PC preconditioner;
OEM_DATA( const int size ) {
VecCreate( PETSC_COMM_WORLD, &b );
auto err = VecSetSizes( b, PETSC_DECIDE, size );
CHKERRXX( err );
VecSetFromOptions( b );
KSPCreate( PETSC_COMM_WORLD, &ksp );
}
~OEM_DATA() {
VecDestroy( &x );
VecDestroy( &b );
MatDestroy( &A );
KSPDestroy( &ksp );
PCDestroy( &preconditioner );
}
};
Vec to_petsc_vec( const double* x, int size ) {
PetscScalar* vec;
Vec v;
VecCreate( PETSC_COMM_WORLD, &v );
VecSetSizes( v, PETSC_DECIDE, size );
VecSetFromOptions( v );
VecGetArray( v, &vec );
std::memcpy( vec, x, size * sizeof( double ) );
VecRestoreArray( v, &vec );
return v;
}
void from_petsc_vec( double* x, Vec v ) {
if( !v ) OPM_THROW( std::runtime_error,
"PETSc CopySolution: Invalid PETSc vector." );
PetscScalar* vec;
PetscInt size;
VecGetLocalSize( v, &size );
VecGetArray( v, &vec );
std::memcpy( x, vec, size * sizeof( double ) );
VecRestoreArray( v, &vec );
}
Mat to_petsc_mat( const int size, const int nonzeros,
const int* ia, const int* ja, const double* sa ) {
Mat A;
auto err = MatCreateSeqAIJWithArrays( PETSC_COMM_WORLD, size, size, (int*)ia, (int*)ja, (double*)sa, &A );
CHKERRXX( err );
return A;
}
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;
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_FALSE );
KSPSolve( t.ksp, t.x, t.b );
KSPGetConvergedReason( t.ksp, &reason );
KSPGetIterationNumber( t.ksp, &its );
KSPGetResidualNorm( t.ksp, &residual );
if( ksp_view )
KSPView( t.ksp, PETSC_VIEWER_STDOUT_WORLD );
err = PetscPrintf( PETSC_COMM_WORLD, "KSP Iterations %D, Final Residual %G\n", its, residual );
CHKERRXX( err );
}
} // anonymous namespace.
LinearSolverPetsc::LinearSolverPetsc(const parameter::ParameterGroup& param)
: 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");
}
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
{
KSPTypeMap ksp(ksp_type_);
KSPType ksp_type = ksp.find(ksp_type_);
PCTypeMap pc(pc_type_);
PCType pc_type = pc.find(pc_type_);
OEM_DATA t( size );
t.A = to_petsc_mat( size, nonzeros, ia, ja, sa );
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.b );
LinearSolverReport rep = {};
rep.converged = true;
return rep;
}
void LinearSolverPetsc::setTolerance(const double /*tol*/)
{
}
double LinearSolverPetsc::getTolerance() const
{
return -1.;
}
} // namespace Opm

92
opm/core/linalg/LinearSolverPetsc.hpp Normal file
View File

@ -0,0 +1,92 @@
/*
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.
/// 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
/// 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 ksp_view_;
double rtol_;
double atol_;
double dtol_;
int maxits_;
};
} // namespace Opm
#endif // OPM_LINEARSOLVERPETSC_HEADER_INCLUDED

13
tests/test_linearsolver.cpp
View File

@ -184,3 +184,16 @@ BOOST_AUTO_TEST_CASE(KAMGTest)
}
#endif
#endif
#if HAVE_PETSC
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);
}
#endif