adapt to the eWoms interface for abstracting sparse matrices

This commit is contained in:
Andreas Lauser 2018-11-12 12:40:11 +01:00
parent 46641d5ace
commit f9104ca3d7
8 changed files with 29 additions and 18 deletions

View File

@ -19,6 +19,7 @@
#ifndef OPM_AMG_HEADER_INCLUDED #ifndef OPM_AMG_HEADER_INCLUDED
#define OPM_AMG_HEADER_INCLUDED #define OPM_AMG_HEADER_INCLUDED
#include <ewoms/linear/matrixblock.hh>
#include <opm/autodiff/ParallelOverlappingILU0.hpp> #include <opm/autodiff/ParallelOverlappingILU0.hpp>
#include <opm/autodiff/CPRPreconditioner.hpp> #include <opm/autodiff/CPRPreconditioner.hpp>
#include <dune/istl/paamg/twolevelmethod.hh> #include <dune/istl/paamg/twolevelmethod.hh>
@ -168,6 +169,12 @@ struct ScalarType<Dune::MatrixBlock<FieldType, ROWS, COLS> >
typedef Dune::MatrixBlock<FieldType, 1, 1> value; typedef Dune::MatrixBlock<FieldType, 1, 1> value;
}; };
template<typename FieldType, int ROWS, int COLS>
struct ScalarType<Ewoms::MatrixBlock<FieldType, ROWS, COLS> >
{
typedef Ewoms::MatrixBlock<FieldType, 1, 1> value;
};
template<typename BlockType, typename Allocator> template<typename BlockType, typename Allocator>
struct ScalarType<Dune::BCRSMatrix<BlockType, Allocator> > struct ScalarType<Dune::BCRSMatrix<BlockType, Allocator> >
{ {

View File

@ -109,6 +109,7 @@ namespace Opm {
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector ; typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector ;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables ; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables ;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
@ -126,8 +127,8 @@ namespace Opm {
static const int temperatureIdx = Indices::temperatureIdx; static const int temperatureIdx = Indices::temperatureIdx;
typedef Dune::FieldVector<Scalar, numEq > VectorBlockType; typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType; typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat; typedef typename SparseMatrixAdapter::IstlMatrix Mat;
typedef Dune::BlockVector<VectorBlockType> BVector; typedef Dune::BlockVector<VectorBlockType> BVector;
typedef ISTLSolverEbos<TypeTag> ISTLSolverType; typedef ISTLSolverEbos<TypeTag> ISTLSolverType;
@ -358,13 +359,13 @@ namespace Opm {
ebosSimulator_.model().linearizer().linearize(); ebosSimulator_.model().linearizer().linearize();
ebosSimulator_.problem().endIteration(); ebosSimulator_.problem().endIteration();
auto& ebosJac = ebosSimulator_.model().linearizer().matrix(); auto& ebosJac = ebosSimulator_.model().linearizer().jacobian();
if (param_.matrix_add_well_contributions_) { if (param_.matrix_add_well_contributions_) {
wellModel().addWellContributions(ebosJac); wellModel().addWellContributions(ebosJac.istlMatrix());
} }
if ( param_.preconditioner_add_well_contributions_ && if ( param_.preconditioner_add_well_contributions_ &&
! param_.matrix_add_well_contributions_ ) { ! param_.matrix_add_well_contributions_ ) {
matrix_for_preconditioner_ .reset(new Mat(ebosJac)); matrix_for_preconditioner_ .reset(new Mat(ebosJac.istlMatrix()));
wellModel().addWellContributions(*matrix_for_preconditioner_); wellModel().addWellContributions(*matrix_for_preconditioner_);
} }
@ -494,7 +495,7 @@ namespace Opm {
// r = [r, r_well], where r is the residual and r_well the well residual. // r = [r, r_well], where r is the residual and r_well the well residual.
// r -= B^T * D^-1 r_well // r -= B^T * D^-1 r_well
auto& ebosJac = ebosSimulator_.model().linearizer().matrix(); auto& ebosJac = ebosSimulator_.model().linearizer().jacobian();
auto& ebosResid = ebosSimulator_.model().linearizer().residual(); auto& ebosResid = ebosSimulator_.model().linearizer().residual();
wellModel().apply(ebosResid); wellModel().apply(ebosResid);
@ -502,13 +503,13 @@ namespace Opm {
// set initial guess // set initial guess
x = 0.0; x = 0.0;
const Mat& actual_mat_for_prec = matrix_for_preconditioner_ ? *matrix_for_preconditioner_.get() : ebosJac; const Mat& actual_mat_for_prec = matrix_for_preconditioner_ ? *matrix_for_preconditioner_.get() : ebosJac.istlMatrix();
// Solve system. // Solve system.
if( isParallel() ) if( isParallel() )
{ {
typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, true > Operator; typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, true > Operator;
auto ebosJacIgnoreOverlap = Mat(ebosJac); auto ebosJacIgnoreOverlap = Mat(ebosJac.istlMatrix());
//remove ghost rows in local matrix //remove ghost rows in local matrix
makeOverlapRowsInvalid(ebosJacIgnoreOverlap); makeOverlapRowsInvalid(ebosJacIgnoreOverlap);
@ -522,7 +523,7 @@ namespace Opm {
else else
{ {
typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, false > Operator; typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, false > Operator;
Operator opA(ebosJac, actual_mat_for_prec, wellModel()); Operator opA(ebosJac.istlMatrix(), actual_mat_for_prec, wellModel());
istlSolver().solve( opA, x, ebosResid ); istlSolver().solve( opA, x, ebosResid );
} }
} }

View File

@ -86,7 +86,7 @@ namespace Opm {
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector; typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix; typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter;
typedef typename Ewoms::BaseAuxiliaryModule<TypeTag>::NeighborSet NeighborSet; typedef typename Ewoms::BaseAuxiliaryModule<TypeTag>::NeighborSet NeighborSet;
@ -105,7 +105,7 @@ namespace Opm {
#else #else
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType; typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
#endif #endif
typedef Dune::BCRSMatrix <MatrixBlockType> Mat; typedef typename SparseMatrixAdapter::IstlMatrix Mat;
typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule; typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
@ -129,7 +129,7 @@ namespace Opm {
void applyInitial() void applyInitial()
{} {}
void linearize(JacobianMatrix& mat , GlobalEqVector& res); void linearize(SparseMatrixAdapter& mat , GlobalEqVector& res);
void postSolve(GlobalEqVector& deltaX) void postSolve(GlobalEqVector& deltaX)
{ {

View File

@ -94,7 +94,7 @@ namespace Opm {
template<typename TypeTag> template<typename TypeTag>
void void
BlackoilWellModel<TypeTag>:: BlackoilWellModel<TypeTag>::
linearize(JacobianMatrix& mat , GlobalEqVector& res) linearize(SparseMatrixAdapter& mat , GlobalEqVector& res)
{ {
if (!localWellsActive()) if (!localWellsActive())
return; return;

View File

@ -71,10 +71,12 @@ namespace Opm
class ISTLSolverEbos : public NewtonIterationBlackoilInterface class ISTLSolverEbos : public NewtonIterationBlackoilInterface
{ {
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) Matrix; typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter;
typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector; typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename SparseMatrixAdapter::IstlMatrix Matrix;
enum { pressureIndex = Indices::pressureSwitchIdx }; enum { pressureIndex = Indices::pressureSwitchIdx };
public: public:

View File

@ -26,13 +26,13 @@
// Define making clear that the simulator supports AMG // Define making clear that the simulator supports AMG
#define FLOW_SUPPORT_AMG !defined(HAVE_UMFPACK) #define FLOW_SUPPORT_AMG !defined(HAVE_UMFPACK)
#include <opm/autodiff/DuneMatrix.hpp>
#include <opm/autodiff/CPRPreconditioner.hpp> #include <opm/autodiff/CPRPreconditioner.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp> #include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
#include <opm/autodiff/NewtonIterationUtilities.hpp> #include <opm/autodiff/NewtonIterationUtilities.hpp>
#include <opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp> #include <opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp>
#include <opm/autodiff/ParallelOverlappingILU0.hpp> #include <opm/autodiff/ParallelOverlappingILU0.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp> #include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/DuneMatrix.hpp>
#include <opm/common/Exceptions.hpp> #include <opm/common/Exceptions.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp> #include <opm/core/linalg/ParallelIstlInformation.hpp>

View File

@ -36,7 +36,7 @@ class WellConnectionAuxiliaryModule
: public Ewoms::BaseAuxiliaryModule<TypeTag> : public Ewoms::BaseAuxiliaryModule<TypeTag>
{ {
typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector; typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix; typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter;
public: public:
@ -116,7 +116,7 @@ public:
void applyInitial() void applyInitial()
{} {}
void linearize(JacobianMatrix& , GlobalEqVector&) void linearize(SparseMatrixAdapter& , GlobalEqVector&)
{ {
// Linearization is done in StandardDenseWells // Linearization is done in StandardDenseWells
} }

View File

@ -80,6 +80,7 @@ namespace Opm
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter;
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
static const int numEq = Indices::numEq; static const int numEq = Indices::numEq;
@ -87,7 +88,7 @@ namespace Opm
typedef Dune::FieldVector<Scalar, numEq > VectorBlockType; typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType; typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat; typedef typename SparseMatrixAdapter::IstlMatrix Mat;
typedef Dune::BlockVector<VectorBlockType> BVector; typedef Dune::BlockVector<VectorBlockType> BVector;
typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval; typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;