diff --git a/opm/models/discretization/common/fvbasediscretization.hh b/opm/models/discretization/common/fvbasediscretization.hh index e2b00b6e9b..b17326a4b5 100644 --- a/opm/models/discretization/common/fvbasediscretization.hh +++ b/opm/models/discretization/common/fvbasediscretization.hh @@ -1787,6 +1787,8 @@ public: } return *adaptationManager_; } + + const DiscreteFunctionSpace& space() const { return space_; } #endif const Opm::Timer& prePostProcessTimer() const diff --git a/opm/models/discretization/common/fvbaselinearizer.hh b/opm/models/discretization/common/fvbaselinearizer.hh index 07690d169f..9295dfd60c 100644 --- a/opm/models/discretization/common/fvbaselinearizer.hh +++ b/opm/models/discretization/common/fvbaselinearizer.hh @@ -95,8 +95,6 @@ class FvBaseLinearizer typedef GlobalEqVector Vector; - typedef typename SparseMatrixAdapter::IstlMatrix IstlMatrix; - enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) }; enum { historySize = GET_PROP_VALUE(TypeTag, TimeDiscHistorySize) }; diff --git a/opm/simulators/linalg/amgxsolverbackend.hh b/opm/simulators/linalg/amgxsolverbackend.hh new file mode 100644 index 0000000000..58acca222d --- /dev/null +++ b/opm/simulators/linalg/amgxsolverbackend.hh @@ -0,0 +1,333 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Ewoms::Linear::AmgXSolverBackend + */ +#ifndef EWOMS_AMGX_SOLVER_BACKEND_HH +#define EWOMS_AMGX_SOLVER_BACKEND_HH + +#include + +#if USE_AMGX_SOLVERS + +#if !HAVE_PETSC || !HAVE_AMGXSOLVER +#error "PETSc and AmgXSolver is needed for the AMGX solver backend" +#endif + +#define DISABLE_AMG_DIRECTSOLVER 1 +#include +#include +#include +#include + +#include + +#include + + +#include +#include + +#include + +#include +#include +#include + +namespace Ewoms { +namespace Linear { +template +class AmgXSolverBackend; +}} // namespace Linear, Ewoms + + +BEGIN_PROPERTIES + +NEW_TYPE_TAG(AmgXSolverBackend); + +SET_TYPE_PROP(AmgXSolverBackend, + LinearSolverBackend, + Ewoms::Linear::AmgXSolverBackend); + +NEW_PROP_TAG(LinearSolverTolerance); +NEW_PROP_TAG(LinearSolverAbsTolerance); +NEW_PROP_TAG(LinearSolverMaxIterations); +NEW_PROP_TAG(LinearSolverVerbosity); +NEW_PROP_TAG(LinearSolverMaxError); +NEW_PROP_TAG(LinearSolverOverlapSize); +//! The order of the sequential preconditioner +NEW_PROP_TAG(PreconditionerOrder); + +//! The relaxation factor of the preconditioner +NEW_PROP_TAG(PreconditionerRelaxation); + +//! Solver mode for AMGX solver +NEW_PROP_TAG(AmgxSolverMode); + +//! Filename for AMGX solver configuration +NEW_PROP_TAG(AmgxSolverConfigFileName); + +//! Filename for DUNE-FEM solver parameters +NEW_PROP_TAG(FemSolverParameterFileName); + +//! make the linear solver shut up by default +SET_INT_PROP(AmgXSolverBackend, LinearSolverVerbosity, 0); + +//! set the default number of maximum iterations for the linear solver +SET_INT_PROP(AmgXSolverBackend, LinearSolverMaxIterations, 1000); + +SET_SCALAR_PROP(AmgXSolverBackend, LinearSolverMaxError, 1e7); + +//! set the default overlap size to 2 +SET_INT_PROP(AmgXSolverBackend, LinearSolverOverlapSize, 2); + +//! set the preconditioner order to 0 by default +SET_INT_PROP(AmgXSolverBackend, PreconditionerOrder, 0); + +//! set the preconditioner relaxation parameter to 1.0 by default +SET_SCALAR_PROP(AmgXSolverBackend, PreconditionerRelaxation, 1.0); + +//! make the linear solver shut up by default +SET_SCALAR_PROP(AmgXSolverBackend, LinearSolverTolerance, 0.01); + +//! make the linear solver shut up by default +SET_SCALAR_PROP(AmgXSolverBackend, LinearSolverAbsTolerance, 0.01); + +//! set default solver mode for AMGX solver configuratrion +SET_STRING_PROP(AmgXSolverBackend, AmgxSolverMode, "dDDI"); + +//! set default filename for AMGX solver configuratrion +SET_STRING_PROP(AmgXSolverBackend, AmgxSolverConfigFileName, "amgxconfig.json"); + +//! set the preconditioner relaxation parameter to 1.0 by default +SET_STRING_PROP(AmgXSolverBackend, FemSolverParameterFileName, ""); + +END_PROPERTIES + +namespace Ewoms { +namespace Linear { +/*! + * \ingroup Linear + * + * \brief Provides a linear solver backend that utilizes the AMG-X package. + */ +template +class AmgXSolverBackend +{ +protected: + typedef typename GET_PROP_TYPE(TypeTag, LinearSolverBackend) Implementation; + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) LinearOperator; + typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + + typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunctionSpace) DiscreteFunctionSpace; + typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunction) DiscreteFunction; + + // discrete function to wrap what is used as Vector in eWoms + typedef Dune::Fem::ISTLBlockVectorDiscreteFunction + VectorWrapperDiscreteFunction; + typedef Dune::Fem::PetscDiscreteFunction + PetscDiscreteFunctionType; + + enum { dimWorld = GridView::dimensionworld }; + +public: + AmgXSolverBackend(Simulator& simulator) + : simulator_(simulator) + , space_(simulator.vanguard().gridPart()) + , amgxSolver_() + , rhs_(nullptr) + , iterations_(0) + { + std::string paramFileName = EWOMS_GET_PARAM(TypeTag, std::string, FemSolverParameterFileName); + if (paramFileName != "") + Dune::Fem::Parameter::append(paramFileName); + } + + ~AmgXSolverBackend() + { + cleanup_(); + if (A_) { + ::Dune::Petsc::MatDestroy(A_.operator ->()); + A_.reset(); + } + } + + /*! + * \brief Register all run-time parameters for the linear solver. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LinearSolverTolerance, + "The maximum allowed error between of the linear solver"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LinearSolverAbsTolerance, + "The maximum allowed error between of the linear solver"); + EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverMaxIterations, + "The maximum number of iterations of the linear solver"); + EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverVerbosity, + "The verbosity level of the linear solver"); + EWOMS_REGISTER_PARAM(TypeTag, unsigned, LinearSolverOverlapSize, + "The size of the algebraic overlap for the linear solver"); + + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LinearSolverMaxError, + "The maximum residual error which the linear solver tolerates" + " without giving up"); + + EWOMS_REGISTER_PARAM(TypeTag, int, PreconditionerOrder, + "The order of the preconditioner"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, PreconditionerRelaxation, + "The relaxation factor of the preconditioner"); + + EWOMS_REGISTER_PARAM(TypeTag, std::string, AmgxSolverMode, + "The name of the solver mode for AMGX"); + EWOMS_REGISTER_PARAM(TypeTag, std::string, AmgxSolverConfigFileName, + "The name of the file which contains the AMGX solver configuration"); + + EWOMS_REGISTER_PARAM(TypeTag, std::string, FemSolverParameterFileName, + "The name of the file which contains the parameters for the DUNE-FEM solvers"); + + } + + /*! + * \brief Causes the solve() method to discared the structure of the linear system of + * equations the next time it is called. + */ + void eraseMatrix() + { cleanup_(); } + + void prepare(const LinearOperator& op, Vector& b) + { + //Scalar linearSolverTolerance = EWOMS_GET_PARAM(TypeTag, Scalar, LinearSolverTolerance); + //Scalar linearSolverAbsTolerance = EWOMS_GET_PARAM(TypeTag, Scalar, LinearSolverAbsTolerance); + + if (!amgxSolver_) { + // reset linear solver + std::string mode = EWOMS_GET_PARAM(TypeTag, std::string, AmgxSolverMode); + std::string solverconfig = EWOMS_GET_PARAM(TypeTag, std::string, AmgxSolverConfigFileName); + + amgxSolver_.reset(new AmgXSolver()); + amgxSolver_->initialize(PETSC_COMM_WORLD, mode, solverconfig); + + if (!A_) { + A_.reset(new Mat()); + } + + // convert MATBAIJ to MATAIJ which is needed by AmgXSolver + ::Dune::Petsc::ErrorCheck(::MatConvert(op.petscMatrix(), MATAIJ, MAT_INITIAL_MATRIX, A_.operator ->())); + // set up the matrix used by AmgX + amgxSolver_->setA(*A_); + } + + // store pointer to right hand side + rhs_ = &b; + } + + /*! + * \brief Actually solve the linear system of equations. + * + * \return true if the residual reduction could be achieved, else false. + */ + bool solve(Vector& x) + { + // wrap x into discrete function X (no copy) + VectorWrapperDiscreteFunction X("FSB::x", space_, x); + VectorWrapperDiscreteFunction B("FSB::rhs", space_, *rhs_); + + if (!petscRhs_) + petscRhs_.reset(new PetscDiscreteFunctionType("AMGX::rhs", space_)); + + if (!petscX_) + petscX_.reset(new PetscDiscreteFunctionType("AMGX::X", space_)); + + petscRhs_->assign(B); + petscX_->assign(X); + + // solve with right hand side rhs and store in x + Vec& p = *(petscX_->petscVec()); + Vec& rhs = *(petscRhs_->petscVec()); + + try { + amgxSolver_->solve(p, rhs); + } + catch (...) { + OPM_THROW(Opm::NumericalIssue, "AmgXSolver: no convergence of linear solver!"); + } + + // copy result to ewoms solution + X.assign(*petscX_); + + amgxSolver_->getIters(iterations_); + + // return the result of the solver + return true; + } + + /*! + * \brief Return number of iterations used during last solve. + */ + size_t iterations () const { + return iterations_; + } + +protected: + Implementation& asImp_() + { return *static_cast(this); } + + const Implementation& asImp_() const + { return *static_cast(this); } + + void cleanup_() + { + if (amgxSolver_) { + amgxSolver_->finalize(); + amgxSolver_.reset(); + } + + rhs_ = nullptr; + + petscRhs_.reset(); + petscX_.reset(); + } + + const Simulator& simulator_; + + DiscreteFunctionSpace space_; + std::unique_ptr A_; + + std::unique_ptr petscRhs_; + std::unique_ptr petscX_; + + std::unique_ptr amgxSolver_; + + Vector* rhs_; + int iterations_; +}; +}} // namespace Linear, Ewoms + +#endif // HAVE_DUNE_FEM + +#endif diff --git a/opm/simulators/linalg/femsolverbackend.hh b/opm/simulators/linalg/femsolverbackend.hh new file mode 100644 index 0000000000..d2a001453d --- /dev/null +++ b/opm/simulators/linalg/femsolverbackend.hh @@ -0,0 +1,361 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::Linear::FemSolverBackend + */ +#ifndef OPM_FEM_SOLVER_BACKEND_HH +#define OPM_FEM_SOLVER_BACKEND_HH + +#include + +#define DISABLE_AMG_DIRECTSOLVER 1 +#include +#include +#if HAVE_PETSC +#include +#endif + +#if HAVE_VIENNACL +#include +#endif + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +namespace Opm { +namespace Linear { +template +class FemSolverBackend; +}} // namespace Linear, Opm + + +BEGIN_PROPERTIES + +NEW_TYPE_TAG(FemSolverBackend); + +SET_TYPE_PROP(FemSolverBackend, + LinearSolverBackend, + Opm::Linear::FemSolverBackend); + +NEW_PROP_TAG(DiscreteFunction); +NEW_PROP_TAG(SparseMatrixAdapter); + +//NEW_PROP_TAG(LinearSolverTolerance); +NEW_PROP_TAG(LinearSolverMaxIterations); +NEW_PROP_TAG(LinearSolverVerbosity); +NEW_PROP_TAG(LinearSolverMaxError); +NEW_PROP_TAG(LinearSolverOverlapSize); +//! The order of the sequential preconditioner +NEW_PROP_TAG(PreconditionerOrder); + +//! The relaxation factor of the preconditioner +NEW_PROP_TAG(PreconditionerRelaxation); + +//! Filename for DUNE-FEM solver parameters +NEW_PROP_TAG(FemSolverParameterFileName); + +//! make the linear solver shut up by default +SET_INT_PROP(FemSolverBackend, LinearSolverVerbosity, 0); + +//! set the default number of maximum iterations for the linear solver +SET_INT_PROP(FemSolverBackend, LinearSolverMaxIterations, 1000); + +SET_SCALAR_PROP(FemSolverBackend, LinearSolverMaxError, 1e7); + +//! set the default overlap size to 2 +SET_INT_PROP(FemSolverBackend, LinearSolverOverlapSize, 2); + +//! set the preconditioner order to 0 by default +SET_INT_PROP(FemSolverBackend, PreconditionerOrder, 0); + +//! set the preconditioner relaxation parameter to 1.0 by default +SET_SCALAR_PROP(FemSolverBackend, PreconditionerRelaxation, 1.0); + +//! set the preconditioner relaxation parameter to 1.0 by default +SET_STRING_PROP(FemSolverBackend, FemSolverParameterFileName, ""); + +//! make the linear solver shut up by default +//SET_SCALAR_PROP(FemSolverBackend, LinearSolverTolerance, 0.01); + +SET_PROP(FemSolverBackend, DiscreteFunction) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunctionSpace) DiscreteFunctionSpace; + + //typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; +public: + // discrete function storing solution data + typedef Dune::Fem::ISTLBlockVectorDiscreteFunction type; +}; + +//! Set the type of a global jacobian matrix for linear solvers that are based on +//! dune-istl. +SET_PROP(FemSolverBackend, SparseMatrixAdapter) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunctionSpace) DiscreteFunctionSpace; + //typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + // discrete function storing solution data + typedef Dune::Fem::ISTLBlockVectorDiscreteFunction DiscreteFunction; + +public: + typedef Opm::Linear::FemSparseRowMatrixAdapter< DiscreteFunction > type; +}; + + +END_PROPERTIES + +namespace Opm { +namespace Linear { +/*! + * \ingroup Linear + * + * \brief Uses the dune-fem infrastructure to solve the linear system of equations. + */ +template +class FemSolverBackend +{ +protected: + typedef typename GET_PROP_TYPE(TypeTag, LinearSolverBackend) Implementation; + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter; + typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + + typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunctionSpace) DiscreteFunctionSpace; + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + // discrete function storing solution data + typedef Dune::Fem::ISTLBlockVectorDiscreteFunction DiscreteFunction; + + // discrete function to wrap what is used as Vector in eWoms + typedef Dune::Fem::ISTLBlockVectorDiscreteFunction + VectorWrapperDiscreteFunction; + + template + struct SolverSelector; + + template + struct SolverSelector< d, FemSparseRowMatrixAdapter< DF > > + { + typedef Dune::Fem::KrylovInverseOperator type; + }; + +#if HAVE_PETSC + template + struct SolverSelector > + { + typedef Dune::Fem::PetscInverseOperator type; + }; +#endif + + template + struct SolverSelector > + { + typedef Dune::Fem::ISTLBICGSTABOp type; + }; + + // select solver type depending on linear operator type + typedef typename SolverSelector<0, SparseMatrixAdapter>::type InverseLinearOperator; + + enum { dimWorld = GridView::dimensionworld }; + +public: + FemSolverBackend(Simulator& simulator) + : simulator_(simulator) + , invOp_() + , rhs_(nullptr) + , space_(simulator.vanguard().gridPart()) + { + std::string paramFileName;// = EWOMS_GET_PARAM(TypeTag, std::string, FemSolverParameterFileName); + if (paramFileName != "") + Dune::Fem::Parameter::append(paramFileName); + else { + // default parameters + Dune::Fem::Parameter::append("fem.solver.errormeasure", "residualreduction"); + Dune::Fem::Parameter::append("fem.solver.verbose", "true"); + Dune::Fem::Parameter::append("fem.verboserank", "0"); + + // Krylov solver + Dune::Fem::Parameter::append("fem.solver.method", "bicgstab"); + + // Preconditioner + Dune::Fem::Parameter::append("fem.solver.preconditioning.method", "ilu"); + Dune::Fem::Parameter::append("fem.solver.preconditioning.level", "0"); + Dune::Fem::Parameter::append("fem.solver.preconditioning.relaxation", "0.9"); + } + } + + ~FemSolverBackend() + { cleanup_(); } + + /*! + * \brief Register all run-time parameters for the linear solver. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LinearSolverTolerance, + "The maximum allowed error between of the linear solver"); + /* + EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverMaxIterations, + "The maximum number of iterations of the linear solver"); + EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverVerbosity, + "The verbosity level of the linear solver"); + EWOMS_REGISTER_PARAM(TypeTag, unsigned, LinearSolverOverlapSize, + "The size of the algebraic overlap for the linear solver"); + + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LinearSolverMaxError, + "The maximum residual error which the linear solver tolerates" + " without giving up"); + + EWOMS_REGISTER_PARAM(TypeTag, int, PreconditionerOrder, + "The order of the preconditioner"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, PreconditionerRelaxation, + "The relaxation factor of the preconditioner"); + + EWOMS_REGISTER_PARAM(TypeTag, std::string, FemSolverParameterFileName, + "The name of the file which contains the parameters for the DUNE-FEM solvers"); + + */ + } + + /*! + * \brief Causes the solve() method to discared the structure of the linear system of + * equations the next time it is called. + */ + void eraseMatrix() + { cleanup_(); } + + void prepare(const SparseMatrixAdapter& op, Vector& b) + { + Scalar linearSolverTolerance = EWOMS_GET_PARAM(TypeTag, Scalar, LinearSolverTolerance); + Scalar linearSolverAbsTolerance = this->simulator_.model().newtonMethod().tolerance() / 100000.0; + //Scalar linearSolverTolerance = 1e-3;//EWOMS_GET_PARAM(TypeTag, Scalar, LinearSolverTolerance); + //Scalar linearSolverAbsTolerance = this->simulator_.model().newtonMethod().tolerance() / 100000.0; + + // reset linear solver + if (!invOp_) + invOp_.reset(new InverseLinearOperator(linearSolverTolerance, linearSolverAbsTolerance)); + + setMatrix(op); + setResidual(b); + + // not needed + asImp_().rescale_(); + } + + /*! + * \brief Assign values to the internal data structure for the residual vector. + * + * This method also cares about synchronizing that vector with the peer processes. + */ + void setResidual(const Vector& b) + { rhs_ = &b ; } + + /*! + * \brief Retrieve the synchronized internal residual vector. + * + * This only deals with entries which are local to the current process. + */ + void getResidual(Vector& b) const + { + assert(rhs_); + b = *rhs_; + } + + /*! + * \brief Sets the values of the residual's Jacobian matrix. + * + * This method also synchronizes the data structure across the processes which are + * involved in the simulation run. + */ + void setMatrix(const SparseMatrixAdapter& op) + { + invOp_->bind(op); + } + + + /*! + * \brief Actually solve the linear system of equations. + * + * \return true if the residual reduction could be achieved, else false. + */ + bool solve(Vector& x) + { + // wrap x into discrete function X (no copy) + VectorWrapperDiscreteFunction X("FSB::x", space_, x); + assert(rhs_); + VectorWrapperDiscreteFunction B("FSB::rhs", space_, *rhs_); + + // solve with right hand side rhs and store in x + (*invOp_)(B, X); + + // return the result of the solver + return invOp_->iterations() < 0 ? false : true; + } + + /*! + * \brief Return number of iterations used during last solve. + */ + size_t iterations () const + { + assert(invOp_); + return invOp_->iterations(); + } + +protected: + Implementation& asImp_() + { return *static_cast(this); } + + const Implementation& asImp_() const + { return *static_cast(this); } + + void rescale_() + { } + + void cleanup_() + { + invOp_.reset(); + rhs_ = nullptr; + } + + const Simulator& simulator_; + std::unique_ptr invOp_; + const Vector* rhs_; + DiscreteFunctionSpace space_; +}; +}} // namespace Linear, Opm + +#endif diff --git a/opm/simulators/linalg/femsparsematrixadapter.hh b/opm/simulators/linalg/femsparsematrixadapter.hh new file mode 100644 index 0000000000..e76af5d380 --- /dev/null +++ b/opm/simulators/linalg/femsparsematrixadapter.hh @@ -0,0 +1,101 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::Linear::FemSparseMatrixAdapter + */ +#ifndef OPM_FEM_SPARSE_MATRIX_ADAPTER_HH +#define OPM_FEM_SPARSE_MATRIX_ADAPTER_HH + +// this code only works with dune-fem available +#if HAVE_DUNE_FEM + +// the following implementation of FemSparseMatrixAdapter only works for +// dune-fem version 2.7 or higher +#if DUNE_VERSION_NEWER(DUNE_FEM, 2, 7) +#include + +#if HAVE_PETSC +#include +#endif + +#include +#include + + +namespace Opm { +namespace Linear { + +/*! + * \ingroup Linear + * \brief A sparse matrix interface backend for linear operators from dune-fem. + * + * \note LinearOperators from dune-fem implement most methods needed for SparseMatrixAdapter + * and here we simply add a few forwarding methods. + */ +template +struct FemSparseMatrixAdapter : public LinearOperator +{ + typedef LinearOperator ParentType; + typedef typename LinearOperator :: MatrixType Matrix; + typedef typename ParentType :: MatrixBlockType MatrixBlock; + + typedef typename LinearOperator :: RangeFunctionType :: RangeFieldType Scalar; + + template + FemSparseMatrixAdapter( const Simulator& simulator ) + : LinearOperator("Opm::Jacobian", simulator.model().space(), simulator.model().space() ) + {} + + void commit() + { + this->flushAssembly(); + } + + template< class LocalBlock > + void addToBlock ( const size_t row, const size_t col, const LocalBlock& block ) + { + this->addBlock( row, col, block ); + } + + void clearRow( const size_t row, const Scalar diag = 1.0 ) { this->unitRow( row ); } +}; + +template +using FemSparseRowMatrixAdapter = FemSparseMatrixAdapter< Dune::Fem::SparseRowLinearOperator< DiscreteFunction, DiscreteFunction > >; + +#if HAVE_PETSC +template +using FemPetscMatrixAdapter = FemSparseMatrixAdapter< Dune::Fem::PetscLinearOperator< DiscreteFunction, DiscreteFunction > >; +#endif + +#if HAVE_DUNE_ISTL +template +using FemISTLMatrixAdapter = FemSparseMatrixAdapter< Dune::Fem::ISTLLinearOperator< DiscreteFunction, DiscreteFunction > >; +#endif + +}} // namespace Linear, Opm + +#endif // DUNE_VERSION_NEWER(DUNE_FEM, 2, 7) +#endif // HAVE_DUNE_FEM +#endif