mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
NewtonBlackølInterationInterleaved: NewtonBlackølInterationInterleavedImpl implements
the Iteration for a fixed number of cell variabled which is then used by the NewtonBlackølInterationInterleaved class via a switch-case over the actually existing numbers.
This commit is contained in:
Robert Kloefkorn
parent
62511215d4
commit
c485e3fdc7
@ -3,6 +3,7 @@
|
|||||||
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
|
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
|
||||||
Copyright 2015 NTNU
|
Copyright 2015 NTNU
|
||||||
Copyright 2015 Statoil AS
|
Copyright 2015 Statoil AS
|
||||||
|
Copyright 2015 IRIS AS
|
||||||
|
|
||||||
This file is part of the Open Porous Media project (OPM).
|
This file is part of the Open Porous Media project (OPM).
|
||||||
|
|
||||||
@ -30,6 +31,13 @@
|
|||||||
#include <opm/core/utility/Exceptions.hpp>
|
#include <opm/core/utility/Exceptions.hpp>
|
||||||
#include <opm/core/linalg/ParallelIstlInformation.hpp>
|
#include <opm/core/linalg/ParallelIstlInformation.hpp>
|
||||||
|
|
||||||
|
#include <dune/istl/scalarproducts.hh>
|
||||||
|
#include <dune/istl/operators.hh>
|
||||||
|
#include <dune/istl/preconditioners.hh>
|
||||||
|
#include <dune/istl/solvers.hh>
|
||||||
|
#include <dune/istl/owneroverlapcopy.hh>
|
||||||
|
#include <dune/istl/paamg/amg.hh>
|
||||||
|
|
||||||
#include <opm/common/utility/platform_dependent/disable_warnings.h>
|
#include <opm/common/utility/platform_dependent/disable_warnings.h>
|
||||||
#if HAVE_UMFPACK
|
#if HAVE_UMFPACK
|
||||||
#include <Eigen/UmfPackSupport>
|
#include <Eigen/UmfPackSupport>
|
||||||
@ -48,156 +56,6 @@ namespace Opm
|
|||||||
typedef ADB::M M;
|
typedef ADB::M M;
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/// Construct a system solver.
|
|
||||||
NewtonIterationBlackoilInterleaved::NewtonIterationBlackoilInterleaved(const parameter::ParameterGroup& param,
|
|
||||||
const boost::any& parallelInformation_arg)
|
|
||||||
: iterations_( 0 ),
|
|
||||||
parallelInformation_(parallelInformation_arg),
|
|
||||||
newton_use_gmres_( param.getDefault("newton_use_gmres", false ) ),
|
|
||||||
linear_solver_reduction_( param.getDefault("linear_solver_reduction", 1e-2 ) ),
|
|
||||||
linear_solver_maxiter_( param.getDefault("linear_solver_maxiter", 50 ) ),
|
|
||||||
linear_solver_restart_( param.getDefault("linear_solver_restart", 40 ) ),
|
|
||||||
linear_solver_verbosity_( param.getDefault("linear_solver_verbosity", 0 ))
|
|
||||||
{
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/// Solve the linear system Ax = b, with A being the
|
|
||||||
/// combined derivative matrix of the residual and b
|
|
||||||
/// being the residual itself.
|
|
||||||
/// \param[in] residual residual object containing A and b.
|
|
||||||
/// \return the solution x
|
|
||||||
NewtonIterationBlackoilInterleaved::SolutionVector
|
|
||||||
NewtonIterationBlackoilInterleaved::computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const
|
|
||||||
{
|
|
||||||
// Build the vector of equations.
|
|
||||||
const int np = residual.material_balance_eq.size();
|
|
||||||
std::vector<ADB> eqs;
|
|
||||||
eqs.reserve(np + 2);
|
|
||||||
for (int phase = 0; phase < np; ++phase) {
|
|
||||||
eqs.push_back(residual.material_balance_eq[phase]);
|
|
||||||
}
|
|
||||||
|
|
||||||
// check if wells are present
|
|
||||||
const bool hasWells = residual.well_flux_eq.size() > 0 ;
|
|
||||||
std::vector<ADB> elim_eqs;
|
|
||||||
if( hasWells )
|
|
||||||
{
|
|
||||||
eqs.push_back(residual.well_flux_eq);
|
|
||||||
eqs.push_back(residual.well_eq);
|
|
||||||
|
|
||||||
// Eliminate the well-related unknowns, and corresponding equations.
|
|
||||||
elim_eqs.reserve(2);
|
|
||||||
elim_eqs.push_back(eqs[np]);
|
|
||||||
eqs = eliminateVariable(eqs, np); // Eliminate well flux unknowns.
|
|
||||||
elim_eqs.push_back(eqs[np]);
|
|
||||||
eqs = eliminateVariable(eqs, np); // Eliminate well bhp unknowns.
|
|
||||||
assert(int(eqs.size()) == np);
|
|
||||||
}
|
|
||||||
|
|
||||||
// Scale material balance equations.
|
|
||||||
for (int phase = 0; phase < np; ++phase) {
|
|
||||||
eqs[phase] = eqs[phase] * residual.matbalscale[phase];
|
|
||||||
}
|
|
||||||
|
|
||||||
// calculating the size for b
|
|
||||||
int size_b = 0;
|
|
||||||
for (int elem = 0; elem < np; ++elem) {
|
|
||||||
const int loc_size = eqs[elem].size();
|
|
||||||
size_b += loc_size;
|
|
||||||
}
|
|
||||||
|
|
||||||
V b(size_b);
|
|
||||||
|
|
||||||
int pos = 0;
|
|
||||||
for (int elem = 0; elem < np; ++elem) {
|
|
||||||
const int loc_size = eqs[elem].size();
|
|
||||||
b.segment(pos, loc_size) = eqs[elem].value();
|
|
||||||
pos += loc_size;
|
|
||||||
}
|
|
||||||
assert(pos == size_b);
|
|
||||||
|
|
||||||
// Create ISTL matrix with interleaved rows and columns (block structured).
|
|
||||||
Mat istlA;
|
|
||||||
formInterleavedSystem(eqs, istlA);
|
|
||||||
|
|
||||||
// Solve reduced system.
|
|
||||||
SolutionVector dx(SolutionVector::Zero(b.size()));
|
|
||||||
|
|
||||||
// Right hand side.
|
|
||||||
const int size = istlA.N();
|
|
||||||
Vector istlb(size);
|
|
||||||
for (int i = 0; i < size; ++i) {
|
|
||||||
istlb[i][0] = b(i);
|
|
||||||
istlb[i][1] = b(size + i);
|
|
||||||
istlb[i][2] = b(2*size + i);
|
|
||||||
}
|
|
||||||
|
|
||||||
// System solution
|
|
||||||
Vector x(istlA.M());
|
|
||||||
x = 0.0;
|
|
||||||
|
|
||||||
Dune::InverseOperatorResult result;
|
|
||||||
// Parallel version is deactivated until we figure out how to do it properly.
|
|
||||||
#if HAVE_MPI
|
|
||||||
if (parallelInformation_.type() == typeid(ParallelISTLInformation))
|
|
||||||
{
|
|
||||||
typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
|
|
||||||
const ParallelISTLInformation& info =
|
|
||||||
boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);
|
|
||||||
Comm istlComm(info.communicator());
|
|
||||||
// As we use a dune-istl with block size np the number of components
|
|
||||||
// per parallel is only one.
|
|
||||||
info.copyValuesTo(istlComm.indexSet(), istlComm.remoteIndices(),
|
|
||||||
size, 1);
|
|
||||||
// Construct operator, scalar product and vectors needed.
|
|
||||||
typedef Dune::OverlappingSchwarzOperator<Mat,Vector,Vector,Comm> Operator;
|
|
||||||
Operator opA(istlA, istlComm);
|
|
||||||
constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, x, istlb, istlComm, result);
|
|
||||||
}
|
|
||||||
else
|
|
||||||
#endif
|
|
||||||
{
|
|
||||||
// Construct operator, scalar product and vectors needed.
|
|
||||||
typedef Dune::MatrixAdapter<Mat,Vector,Vector> Operator;
|
|
||||||
Operator opA(istlA);
|
|
||||||
Dune::Amg::SequentialInformation info;
|
|
||||||
constructPreconditionerAndSolve(opA, x, istlb, info, result);
|
|
||||||
}
|
|
||||||
|
|
||||||
// store number of iterations
|
|
||||||
iterations_ = result.iterations;
|
|
||||||
|
|
||||||
// Check for failure of linear solver.
|
|
||||||
if (!result.converged) {
|
|
||||||
OPM_THROW(LinearSolverProblem, "Convergence failure for linear solver.");
|
|
||||||
}
|
|
||||||
|
|
||||||
// Copy solver output to dx.
|
|
||||||
for (int i = 0; i < size; ++i) {
|
|
||||||
dx(i) = x[i][0];
|
|
||||||
dx(size + i) = x[i][1];
|
|
||||||
dx(2*size + i) = x[i][2];
|
|
||||||
}
|
|
||||||
|
|
||||||
if ( hasWells ) {
|
|
||||||
// Compute full solution using the eliminated equations.
|
|
||||||
// Recovery in inverse order of elimination.
|
|
||||||
dx = recoverVariable(elim_eqs[1], dx, np);
|
|
||||||
dx = recoverVariable(elim_eqs[0], dx, np);
|
|
||||||
}
|
|
||||||
return dx;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
namespace detail {
|
namespace detail {
|
||||||
/**
|
/**
|
||||||
* Simple binary operator that always returns 0.1
|
* Simple binary operator that always returns 0.1
|
||||||
@ -212,79 +70,429 @@ namespace Opm
|
|||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
void NewtonIterationBlackoilInterleaved::formInterleavedSystem(const std::vector<ADB>& eqs,
|
/// This class solves the fully implicit black-oil system by
|
||||||
Mat& istlA) const
|
/// solving the reduced system (after eliminating well variables)
|
||||||
|
/// as a block-structured matrix (one block for all cell variables) for a fixed
|
||||||
|
/// number of cell variables np .
|
||||||
|
template <int np>
|
||||||
|
class NewtonIterationBlackoilInterleavedImpl : public NewtonIterationBlackoilInterface
|
||||||
{
|
{
|
||||||
const int np = eqs.size();
|
typedef Dune::FieldVector<double, np > VectorBlockType;
|
||||||
// Find sparsity structure as union of basic block sparsity structures,
|
typedef Dune::FieldMatrix<double, np, np> MatrixBlockType;
|
||||||
// corresponding to the jacobians with respect to pressure.
|
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
|
||||||
// Use our custom PointOneOp to get to the union structure.
|
typedef Dune::BlockVector<VectorBlockType> Vector;
|
||||||
// Note that we only iterate over the pressure derivatives on purpose.
|
|
||||||
Eigen::SparseMatrix<double, Eigen::ColMajor> col_major = eqs[0].derivative()[0].getSparse();
|
public:
|
||||||
detail::PointOneOp<double> point_one;
|
typedef NewtonIterationBlackoilInterface :: SolutionVector SolutionVector;
|
||||||
for (int phase = 1; phase < np; ++phase) {
|
/// Construct a system solver.
|
||||||
const AutoDiffMatrix::SparseRep& mat = eqs[phase].derivative()[0].getSparse();
|
/// \param[in] param parameters controlling the behaviour of
|
||||||
col_major = col_major.binaryExpr(mat, point_one);
|
/// the preconditioning and choice of
|
||||||
|
/// linear solvers.
|
||||||
|
/// Parameters:
|
||||||
|
/// cpr_relax (default 1.0) relaxation for the preconditioner
|
||||||
|
/// cpr_ilu_n (default 0) use ILU(n) for preconditioning of the linear system
|
||||||
|
/// cpr_use_amg (default false) if true, use AMG preconditioner for elliptic part
|
||||||
|
/// cpr_use_bicgstab (default true) if true, use BiCGStab (else use CG) for elliptic part
|
||||||
|
/// \param[in] parallelInformation In the case of a parallel run
|
||||||
|
/// with dune-istl the information about the parallelization.
|
||||||
|
NewtonIterationBlackoilInterleavedImpl(const parameter::ParameterGroup& param,
|
||||||
|
const boost::any& parallelInformation_arg=boost::any())
|
||||||
|
: iterations_( 0 ),
|
||||||
|
parallelInformation_(parallelInformation_arg),
|
||||||
|
newton_use_gmres_( param.getDefault("newton_use_gmres", false ) ),
|
||||||
|
linear_solver_reduction_( param.getDefault("linear_solver_reduction", 1e-2 ) ),
|
||||||
|
linear_solver_maxiter_( param.getDefault("linear_solver_maxiter", 50 ) ),
|
||||||
|
linear_solver_restart_( param.getDefault("linear_solver_restart", 40 ) ),
|
||||||
|
linear_solver_verbosity_( param.getDefault("linear_solver_verbosity", 0 ))
|
||||||
|
{
|
||||||
}
|
}
|
||||||
|
|
||||||
// Automatically convert the column major structure to a row-major structure
|
/// Solve the system of linear equations Ax = b, with A being the
|
||||||
Eigen::SparseMatrix<double, Eigen::RowMajor> row_major = col_major;
|
/// combined derivative matrix of the residual and b
|
||||||
|
/// being the residual itself.
|
||||||
|
/// \param[in] residual residual object containing A and b.
|
||||||
|
/// \return the solution x
|
||||||
|
|
||||||
const int size = row_major.rows();
|
/// \copydoc NewtonIterationBlackoilInterface::iterations
|
||||||
assert(size == row_major.cols());
|
int iterations () const { return iterations_; }
|
||||||
|
|
||||||
// Create ISTL matrix with interleaved rows and columns (block structured).
|
/// \copydoc NewtonIterationBlackoilInterface::parallelInformation
|
||||||
assert(np == 3);
|
const boost::any& parallelInformation() const { return parallelInformation_; }
|
||||||
istlA.setSize(row_major.rows(), row_major.cols(), row_major.nonZeros());
|
|
||||||
istlA.setBuildMode(Mat::row_wise);
|
public:
|
||||||
const int* ia = row_major.outerIndexPtr();
|
/// \brief construct the CPR preconditioner and the solver.
|
||||||
const int* ja = row_major.innerIndexPtr();
|
/// \tparam P The type of the parallel information.
|
||||||
for (Mat::CreateIterator row = istlA.createbegin(); row != istlA.createend(); ++row) {
|
/// \param parallelInformation the information about the parallelization.
|
||||||
const int ri = row.index();
|
template<int category=Dune::SolverCategory::sequential, class O, class POrComm>
|
||||||
for (int i = ia[ri]; i < ia[ri + 1]; ++i) {
|
void constructPreconditionerAndSolve(O& opA,
|
||||||
row.insert(ja[i]);
|
Vector& x, Vector& istlb,
|
||||||
|
const POrComm& parallelInformation_arg,
|
||||||
|
Dune::InverseOperatorResult& result) const
|
||||||
|
{
|
||||||
|
// Construct scalar product.
|
||||||
|
typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
|
||||||
|
typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
|
||||||
|
SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg));
|
||||||
|
|
||||||
|
// Construct preconditioner.
|
||||||
|
auto precond = constructPrecond(opA, parallelInformation_arg);
|
||||||
|
|
||||||
|
// Communicate if parallel.
|
||||||
|
parallelInformation_arg.copyOwnerToAll(istlb, istlb);
|
||||||
|
|
||||||
|
// Solve.
|
||||||
|
solve(opA, x, istlb, *sp, *precond, result);
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
typedef Dune::SeqILU0<Mat, Vector, Vector> SeqPreconditioner;
|
||||||
|
|
||||||
|
template <class Operator>
|
||||||
|
std::unique_ptr<SeqPreconditioner> constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const
|
||||||
|
{
|
||||||
|
const double relax = 1.0;
|
||||||
|
std::unique_ptr<SeqPreconditioner> precond(new SeqPreconditioner(opA.getmat(), relax));
|
||||||
|
return precond;
|
||||||
|
}
|
||||||
|
|
||||||
|
#if HAVE_MPI
|
||||||
|
typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
|
||||||
|
typedef Dune::BlockPreconditioner<Vector, Vector, Comm, SeqPreconditioner> ParPreconditioner;
|
||||||
|
|
||||||
|
template <class Operator>
|
||||||
|
std::unique_ptr<ParPreconditioner,
|
||||||
|
AdditionalObjectDeleter<SeqPreconditioner> >
|
||||||
|
constructPrecond(Operator& opA, const Comm& comm) const
|
||||||
|
{
|
||||||
|
typedef AdditionalObjectDeleter<SeqPreconditioner> Deleter;
|
||||||
|
typedef std::unique_ptr<ParPreconditioner, Deleter> Pointer;
|
||||||
|
int ilu_setup_successful = 1;
|
||||||
|
std::string message;
|
||||||
|
const double relax = 1.0;
|
||||||
|
SeqPreconditioner* seq_precond = nullptr;
|
||||||
|
try {
|
||||||
|
seq_precond = new SeqPreconditioner(opA.getmat(), relax);
|
||||||
|
}
|
||||||
|
catch ( Dune::MatrixBlockError error )
|
||||||
|
{
|
||||||
|
message = error.what();
|
||||||
|
std::cerr<<"Exception occured on process " <<
|
||||||
|
comm.communicator().rank() << " during " <<
|
||||||
|
"setup of ILU0 preconditioner with message: " <<
|
||||||
|
message<<std::endl;
|
||||||
|
ilu_setup_successful = 0;
|
||||||
|
}
|
||||||
|
// Check whether there was a problem on some process
|
||||||
|
if ( comm.communicator().min(ilu_setup_successful) == 0 )
|
||||||
|
{
|
||||||
|
if ( seq_precond ) // not null if constructor succeeded
|
||||||
|
{
|
||||||
|
// prevent memory leak
|
||||||
|
delete seq_precond;
|
||||||
|
throw Dune::MatrixBlockError();
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return Pointer(new ParPreconditioner(*seq_precond, comm),
|
||||||
|
Deleter(*seq_precond));
|
||||||
|
}
|
||||||
|
#endif
|
||||||
|
|
||||||
|
/// \brief Solve the system using the given preconditioner and scalar product.
|
||||||
|
template <class Operator, class ScalarProd, class Precond>
|
||||||
|
void solve(Operator& opA, Vector& x, Vector& istlb, ScalarProd& sp, Precond& precond, Dune::InverseOperatorResult& result) const
|
||||||
|
{
|
||||||
|
// TODO: Revise when linear solvers interface opm-core is done
|
||||||
|
// Construct linear solver.
|
||||||
|
// GMRes solver
|
||||||
|
if ( newton_use_gmres_ ) {
|
||||||
|
Dune::RestartedGMResSolver<Vector> linsolve(opA, sp, precond,
|
||||||
|
linear_solver_reduction_, linear_solver_restart_, linear_solver_maxiter_, linear_solver_verbosity_);
|
||||||
|
// Solve system.
|
||||||
|
linsolve.apply(x, istlb, result);
|
||||||
|
}
|
||||||
|
else { // BiCGstab solver
|
||||||
|
Dune::BiCGSTABSolver<Vector> linsolve(opA, sp, precond,
|
||||||
|
linear_solver_reduction_, linear_solver_maxiter_, linear_solver_verbosity_);
|
||||||
|
// Solve system.
|
||||||
|
linsolve.apply(x, istlb, result);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
// Set all blocks to zero.
|
void formInterleavedSystem(const std::vector<LinearisedBlackoilResidual::ADB>& eqs,
|
||||||
for (int row = 0; row < size; ++row) {
|
Mat& istlA) const
|
||||||
for (int col_ix = ia[row]; col_ix < ia[row + 1]; ++col_ix) {
|
{
|
||||||
const int col = ja[col_ix];
|
assert( np == int(eqs.size()) );
|
||||||
istlA[row][col] = 0.0;
|
// Find sparsity structure as union of basic block sparsity structures,
|
||||||
|
// corresponding to the jacobians with respect to pressure.
|
||||||
|
// Use our custom PointOneOp to get to the union structure.
|
||||||
|
// Note that we only iterate over the pressure derivatives on purpose.
|
||||||
|
Eigen::SparseMatrix<double, Eigen::ColMajor> col_major = eqs[0].derivative()[0].getSparse();
|
||||||
|
detail::PointOneOp<double> point_one;
|
||||||
|
for (int phase = 1; phase < np; ++phase) {
|
||||||
|
const AutoDiffMatrix::SparseRep& mat = eqs[phase].derivative()[0].getSparse();
|
||||||
|
col_major = col_major.binaryExpr(mat, point_one);
|
||||||
}
|
}
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
// Automatically convert the column major structure to a row-major structure
|
||||||
* Go through all jacobians, and insert in correct spot
|
Eigen::SparseMatrix<double, Eigen::RowMajor> row_major = col_major;
|
||||||
*
|
|
||||||
* The straight forward way to do this would be to run through each
|
const int size = row_major.rows();
|
||||||
* element in the output matrix, and set all block entries by gathering
|
assert(size == row_major.cols());
|
||||||
* from all "input matrices" (derivatives).
|
|
||||||
*
|
// Create ISTL matrix with interleaved rows and columns (block structured).
|
||||||
* A faster alternative is to instead run through each "input matrix" and
|
istlA.setSize(row_major.rows(), row_major.cols(), row_major.nonZeros());
|
||||||
* insert its elements in the correct spot in the output matrix.
|
istlA.setBuildMode(Mat::row_wise);
|
||||||
*
|
const int* ia = row_major.outerIndexPtr();
|
||||||
*/
|
const int* ja = row_major.innerIndexPtr();
|
||||||
for (int col = 0; col < size; ++col) {
|
const typename Mat::CreateIterator endrow = istlA.createend();
|
||||||
for (int p1 = 0; p1 < np; ++p1) {
|
for (typename Mat::CreateIterator row = istlA.createbegin(); row != endrow; ++row) {
|
||||||
for (int p2 = 0; p2 < np; ++p2) {
|
const int ri = row.index();
|
||||||
// Note that that since these are CSC and not CSR matrices,
|
for (int i = ia[ri]; i < ia[ri + 1]; ++i) {
|
||||||
// ja contains row numbers instead of column numbers.
|
row.insert(ja[i]);
|
||||||
const AutoDiffMatrix::SparseRep& s = eqs[p1].derivative()[p2].getSparse();
|
}
|
||||||
const int* ia = s.outerIndexPtr();
|
}
|
||||||
const int* ja = s.innerIndexPtr();
|
|
||||||
const double* sa = s.valuePtr();
|
// Set all blocks to zero.
|
||||||
for (int elem_ix = ia[col]; elem_ix < ia[col + 1]; ++elem_ix) {
|
for (int row = 0; row < size; ++row) {
|
||||||
const int row = ja[elem_ix];
|
for (int col_ix = ia[row]; col_ix < ia[row + 1]; ++col_ix) {
|
||||||
istlA[row][col][p1][p2] = sa[elem_ix];
|
const int col = ja[col_ix];
|
||||||
|
istlA[row][col] = 0.0;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Go through all jacobians, and insert in correct spot
|
||||||
|
*
|
||||||
|
* The straight forward way to do this would be to run through each
|
||||||
|
* element in the output matrix, and set all block entries by gathering
|
||||||
|
* from all "input matrices" (derivatives).
|
||||||
|
*
|
||||||
|
* A faster alternative is to instead run through each "input matrix" and
|
||||||
|
* insert its elements in the correct spot in the output matrix.
|
||||||
|
*
|
||||||
|
*/
|
||||||
|
for (int col = 0; col < size; ++col) {
|
||||||
|
for (int p1 = 0; p1 < np; ++p1) {
|
||||||
|
for (int p2 = 0; p2 < np; ++p2) {
|
||||||
|
// Note that that since these are CSC and not CSR matrices,
|
||||||
|
// ja contains row numbers instead of column numbers.
|
||||||
|
const AutoDiffMatrix::SparseRep& s = eqs[p1].derivative()[p2].getSparse();
|
||||||
|
const int* ia = s.outerIndexPtr();
|
||||||
|
const int* ja = s.innerIndexPtr();
|
||||||
|
const double* sa = s.valuePtr();
|
||||||
|
for (int elem_ix = ia[col]; elem_ix < ia[col + 1]; ++elem_ix) {
|
||||||
|
const int row = ja[elem_ix];
|
||||||
|
istlA[row][col][p1][p2] = sa[elem_ix];
|
||||||
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/// Solve the linear system Ax = b, with A being the
|
||||||
|
/// combined derivative matrix of the residual and b
|
||||||
|
/// being the residual itself.
|
||||||
|
/// \param[in] residual residual object containing A and b.
|
||||||
|
/// \return the solution x
|
||||||
|
SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const
|
||||||
|
{
|
||||||
|
// Build the vector of equations.
|
||||||
|
//const int np = residual.material_balance_eq.size();
|
||||||
|
assert( np == int(residual.material_balance_eq.size()) );
|
||||||
|
std::vector<ADB> eqs;
|
||||||
|
eqs.reserve(np + 2);
|
||||||
|
for (int phase = 0; phase < np; ++phase) {
|
||||||
|
eqs.push_back(residual.material_balance_eq[phase]);
|
||||||
|
}
|
||||||
|
|
||||||
|
// check if wells are present
|
||||||
|
const bool hasWells = residual.well_flux_eq.size() > 0 ;
|
||||||
|
std::vector<ADB> elim_eqs;
|
||||||
|
if( hasWells )
|
||||||
|
{
|
||||||
|
eqs.push_back(residual.well_flux_eq);
|
||||||
|
eqs.push_back(residual.well_eq);
|
||||||
|
|
||||||
|
// Eliminate the well-related unknowns, and corresponding equations.
|
||||||
|
elim_eqs.reserve(2);
|
||||||
|
elim_eqs.push_back(eqs[np]);
|
||||||
|
eqs = eliminateVariable(eqs, np); // Eliminate well flux unknowns.
|
||||||
|
elim_eqs.push_back(eqs[np]);
|
||||||
|
eqs = eliminateVariable(eqs, np); // Eliminate well bhp unknowns.
|
||||||
|
assert(int(eqs.size()) == np);
|
||||||
|
}
|
||||||
|
|
||||||
|
// Scale material balance equations.
|
||||||
|
for (int phase = 0; phase < np; ++phase) {
|
||||||
|
eqs[phase] = eqs[phase] * residual.matbalscale[phase];
|
||||||
|
}
|
||||||
|
|
||||||
|
// calculating the size for b
|
||||||
|
int size_b = 0;
|
||||||
|
for (int elem = 0; elem < np; ++elem) {
|
||||||
|
const int loc_size = eqs[elem].size();
|
||||||
|
size_b += loc_size;
|
||||||
|
}
|
||||||
|
|
||||||
|
V b(size_b);
|
||||||
|
|
||||||
|
int pos = 0;
|
||||||
|
for (int elem = 0; elem < np; ++elem) {
|
||||||
|
const int loc_size = eqs[elem].size();
|
||||||
|
b.segment(pos, loc_size) = eqs[elem].value();
|
||||||
|
pos += loc_size;
|
||||||
|
}
|
||||||
|
assert(pos == size_b);
|
||||||
|
|
||||||
|
// Create ISTL matrix with interleaved rows and columns (block structured).
|
||||||
|
Mat istlA;
|
||||||
|
formInterleavedSystem(eqs, istlA);
|
||||||
|
|
||||||
|
// Solve reduced system.
|
||||||
|
SolutionVector dx(SolutionVector::Zero(b.size()));
|
||||||
|
|
||||||
|
// Right hand side.
|
||||||
|
const int size = istlA.N();
|
||||||
|
Vector istlb(size);
|
||||||
|
for (int i = 0; i < size; ++i) {
|
||||||
|
for( int p = 0, idx = i; p<np; ++p, idx += size ) {
|
||||||
|
istlb[i][p] = b(idx);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// System solution
|
||||||
|
Vector x(istlA.M());
|
||||||
|
x = 0.0;
|
||||||
|
|
||||||
|
Dune::InverseOperatorResult result;
|
||||||
|
// Parallel version is deactivated until we figure out how to do it properly.
|
||||||
|
#if HAVE_MPI
|
||||||
|
if (parallelInformation_.type() == typeid(ParallelISTLInformation))
|
||||||
|
{
|
||||||
|
typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
|
||||||
|
const ParallelISTLInformation& info =
|
||||||
|
boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);
|
||||||
|
Comm istlComm(info.communicator());
|
||||||
|
// As we use a dune-istl with block size np the number of components
|
||||||
|
// per parallel is only one.
|
||||||
|
info.copyValuesTo(istlComm.indexSet(), istlComm.remoteIndices(),
|
||||||
|
size, 1);
|
||||||
|
// Construct operator, scalar product and vectors needed.
|
||||||
|
typedef Dune::OverlappingSchwarzOperator<Mat,Vector,Vector,Comm> Operator;
|
||||||
|
Operator opA(istlA, istlComm);
|
||||||
|
constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, x, istlb, istlComm, result);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
#endif
|
||||||
|
{
|
||||||
|
// Construct operator, scalar product and vectors needed.
|
||||||
|
typedef Dune::MatrixAdapter<Mat,Vector,Vector> Operator;
|
||||||
|
Operator opA(istlA);
|
||||||
|
Dune::Amg::SequentialInformation info;
|
||||||
|
constructPreconditionerAndSolve(opA, x, istlb, info, result);
|
||||||
|
}
|
||||||
|
|
||||||
|
// store number of iterations
|
||||||
|
iterations_ = result.iterations;
|
||||||
|
|
||||||
|
// Check for failure of linear solver.
|
||||||
|
if (!result.converged) {
|
||||||
|
OPM_THROW(LinearSolverProblem, "Convergence failure for linear solver.");
|
||||||
|
}
|
||||||
|
|
||||||
|
// Copy solver output to dx.
|
||||||
|
for (int i = 0; i < size; ++i) {
|
||||||
|
for( int p=0, idx = i; p<np; ++p, idx += size ) {
|
||||||
|
dx(idx) = x[i][p];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
if ( hasWells ) {
|
||||||
|
// Compute full solution using the eliminated equations.
|
||||||
|
// Recovery in inverse order of elimination.
|
||||||
|
dx = recoverVariable(elim_eqs[1], dx, np);
|
||||||
|
dx = recoverVariable(elim_eqs[0], dx, np);
|
||||||
|
}
|
||||||
|
return dx;
|
||||||
|
}
|
||||||
|
|
||||||
|
protected:
|
||||||
|
mutable int iterations_;
|
||||||
|
boost::any parallelInformation_;
|
||||||
|
|
||||||
|
const bool newton_use_gmres_;
|
||||||
|
const double linear_solver_reduction_;
|
||||||
|
const int linear_solver_maxiter_;
|
||||||
|
const int linear_solver_restart_;
|
||||||
|
const int linear_solver_verbosity_;
|
||||||
|
|
||||||
|
}; // end NewtonIterationBlackoilInterleavedImpl
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
/// Construct a system solver.
|
||||||
|
NewtonIterationBlackoilInterleaved::NewtonIterationBlackoilInterleaved(const parameter::ParameterGroup& param,
|
||||||
|
const boost::any& parallelInformation_arg)
|
||||||
|
: newtonIncrement_( 6 ),
|
||||||
|
param_( param ),
|
||||||
|
parallelInformation_(parallelInformation_arg),
|
||||||
|
iterations_( 0 )
|
||||||
|
{
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
/// Solve the linear system Ax = b, with A being the
|
||||||
|
/// combined derivative matrix of the residual and b
|
||||||
|
/// being the residual itself.
|
||||||
|
/// \param[in] residual residual object containing A and b.
|
||||||
|
/// \return the solution x
|
||||||
|
NewtonIterationBlackoilInterleaved::SolutionVector
|
||||||
|
NewtonIterationBlackoilInterleaved::computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const
|
||||||
|
{
|
||||||
|
// get np and call appropriate template method
|
||||||
|
const int np = residual.material_balance_eq.size();
|
||||||
|
switch( np )
|
||||||
|
{
|
||||||
|
case 2:
|
||||||
|
return computeNewtonIncrementImpl< 2 >( residual );
|
||||||
|
case 3:
|
||||||
|
return computeNewtonIncrementImpl< 3 >( residual );
|
||||||
|
case 4:
|
||||||
|
return computeNewtonIncrementImpl< 4 >( residual );
|
||||||
|
case 5:
|
||||||
|
return computeNewtonIncrementImpl< 5 >( residual );
|
||||||
|
case 6:
|
||||||
|
return computeNewtonIncrementImpl< 6 >( residual );
|
||||||
|
default:
|
||||||
|
OPM_THROW(std::runtime_error,"computeNewtonIncrement: number of variables not supported yet. Consult the code to add a case for np = " << np);
|
||||||
|
return SolutionVector();
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/// Solve the linear system Ax = b, with A being the
|
||||||
|
/// combined derivative matrix of the residual and b
|
||||||
|
/// being the residual itself.
|
||||||
|
/// \param[in] residual residual object containing A and b.
|
||||||
|
/// \return the solution x
|
||||||
|
template <int np>
|
||||||
|
NewtonIterationBlackoilInterleaved::SolutionVector
|
||||||
|
NewtonIterationBlackoilInterleaved::computeNewtonIncrementImpl(const LinearisedBlackoilResidual& residual) const
|
||||||
|
{
|
||||||
|
assert( np < int(newtonIncrement_.size()) );
|
||||||
|
// create NewtonIncrement with fixed np
|
||||||
|
if( ! newtonIncrement_[ np ] )
|
||||||
|
newtonIncrement_[ np ].reset( new NewtonIterationBlackoilInterleavedImpl< np >( param_, parallelInformation_ ) );
|
||||||
|
|
||||||
|
// compute newton increment
|
||||||
|
SolutionVector dx = newtonIncrement_[ np ]->computeNewtonIncrement( residual );
|
||||||
|
// get number of linear iterations
|
||||||
|
iterations_ = newtonIncrement_[ np ]->iterations();
|
||||||
|
return dx;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
const boost::any& NewtonIterationBlackoilInterleaved::parallelInformation() const
|
const boost::any& NewtonIterationBlackoilInterleaved::parallelInformation() const
|
||||||
@ -294,7 +502,5 @@ namespace Opm
|
|||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
} // namespace Opm
|
} // namespace Opm
|
||||||
|
|
||||||
|
|||||||
@ -4,6 +4,7 @@
|
|||||||
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
|
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
|
||||||
Copyright 2015 NTNU
|
Copyright 2015 NTNU
|
||||||
Copyright 2015 Statoil AS
|
Copyright 2015 Statoil AS
|
||||||
|
Copyright 2015 IRIS AS
|
||||||
|
|
||||||
This file is part of the Open Porous Media project (OPM).
|
This file is part of the Open Porous Media project (OPM).
|
||||||
|
|
||||||
@ -32,13 +33,6 @@
|
|||||||
|
|
||||||
#include <opm/common/utility/platform_dependent/disable_warnings.h>
|
#include <opm/common/utility/platform_dependent/disable_warnings.h>
|
||||||
|
|
||||||
#include <dune/istl/scalarproducts.hh>
|
|
||||||
#include <dune/istl/operators.hh>
|
|
||||||
#include <dune/istl/preconditioners.hh>
|
|
||||||
#include <dune/istl/solvers.hh>
|
|
||||||
#include <dune/istl/owneroverlapcopy.hh>
|
|
||||||
#include <dune/istl/paamg/amg.hh>
|
|
||||||
|
|
||||||
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
|
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
|
||||||
|
|
||||||
#include <memory>
|
#include <memory>
|
||||||
@ -51,11 +45,6 @@ namespace Opm
|
|||||||
/// as a block-structured matrix (one block for all cell variables).
|
/// as a block-structured matrix (one block for all cell variables).
|
||||||
class NewtonIterationBlackoilInterleaved : public NewtonIterationBlackoilInterface
|
class NewtonIterationBlackoilInterleaved : public NewtonIterationBlackoilInterface
|
||||||
{
|
{
|
||||||
typedef Dune::FieldVector<double, 3 > VectorBlockType;
|
|
||||||
typedef Dune::FieldMatrix<double, 3, 3> MatrixBlockType;
|
|
||||||
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
|
|
||||||
typedef Dune::BlockVector<VectorBlockType> Vector;
|
|
||||||
|
|
||||||
public:
|
public:
|
||||||
|
|
||||||
/// Construct a system solver.
|
/// Construct a system solver.
|
||||||
@ -86,118 +75,13 @@ namespace Opm
|
|||||||
virtual const boost::any& parallelInformation() const;
|
virtual const boost::any& parallelInformation() const;
|
||||||
|
|
||||||
private:
|
private:
|
||||||
|
template <int np>
|
||||||
|
SolutionVector computeNewtonIncrementImpl( const LinearisedBlackoilResidual& residual ) const;
|
||||||
|
|
||||||
/// \brief construct the CPR preconditioner and the solver.
|
mutable std::vector< std::unique_ptr< NewtonIterationBlackoilInterface > > newtonIncrement_;
|
||||||
/// \tparam P The type of the parallel information.
|
const parameter::ParameterGroup& param_;
|
||||||
/// \param parallelInformation the information about the parallelization.
|
|
||||||
template<int category=Dune::SolverCategory::sequential, class O, class POrComm>
|
|
||||||
void constructPreconditionerAndSolve(O& opA,
|
|
||||||
Vector& x, Vector& istlb,
|
|
||||||
const POrComm& parallelInformation_arg,
|
|
||||||
Dune::InverseOperatorResult& result) const
|
|
||||||
{
|
|
||||||
// Construct scalar product.
|
|
||||||
typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
|
|
||||||
typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
|
|
||||||
SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg));
|
|
||||||
|
|
||||||
// Construct preconditioner.
|
|
||||||
auto precond = constructPrecond(opA, parallelInformation_arg);
|
|
||||||
|
|
||||||
// Communicate if parallel.
|
|
||||||
parallelInformation_arg.copyOwnerToAll(istlb, istlb);
|
|
||||||
|
|
||||||
// Solve.
|
|
||||||
solve(opA, x, istlb, *sp, *precond, result);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
typedef Dune::SeqILU0<Mat, Vector, Vector> SeqPreconditioner;
|
|
||||||
|
|
||||||
template <class Operator>
|
|
||||||
std::unique_ptr<SeqPreconditioner> constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const
|
|
||||||
{
|
|
||||||
const double relax = 1.0;
|
|
||||||
std::unique_ptr<SeqPreconditioner>
|
|
||||||
precond(new SeqPreconditioner(opA.getmat(), relax));
|
|
||||||
return precond;
|
|
||||||
}
|
|
||||||
|
|
||||||
#if HAVE_MPI
|
|
||||||
typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
|
|
||||||
typedef Dune::BlockPreconditioner<Vector, Vector, Comm, SeqPreconditioner> ParPreconditioner;
|
|
||||||
|
|
||||||
template <class Operator>
|
|
||||||
std::unique_ptr<ParPreconditioner,
|
|
||||||
AdditionalObjectDeleter<SeqPreconditioner> >
|
|
||||||
constructPrecond(Operator& opA, const Comm& comm) const
|
|
||||||
{
|
|
||||||
typedef AdditionalObjectDeleter<SeqPreconditioner> Deleter;
|
|
||||||
typedef std::unique_ptr<ParPreconditioner, Deleter> Pointer;
|
|
||||||
int ilu_setup_successful = 1;
|
|
||||||
std::string message;
|
|
||||||
const double relax = 1.0;
|
|
||||||
SeqPreconditioner* seq_precond = nullptr;
|
|
||||||
try {
|
|
||||||
seq_precond = new SeqPreconditioner(opA.getmat(),
|
|
||||||
relax);
|
|
||||||
}
|
|
||||||
catch ( Dune::MatrixBlockError error )
|
|
||||||
{
|
|
||||||
message = error.what();
|
|
||||||
std::cerr<<"Exception occured on process " <<
|
|
||||||
comm.communicator().rank() << " during " <<
|
|
||||||
"setup of ILU0 preconditioner with message: " <<
|
|
||||||
message<<std::endl;
|
|
||||||
ilu_setup_successful = 0;
|
|
||||||
}
|
|
||||||
// Check whether there was a problem on some process
|
|
||||||
if ( comm.communicator().min(ilu_setup_successful) == 0 )
|
|
||||||
{
|
|
||||||
if ( seq_precond ) // not null if constructor succeeded
|
|
||||||
{
|
|
||||||
// prevent memory leak
|
|
||||||
delete seq_precond;
|
|
||||||
}
|
|
||||||
throw Dune::MatrixBlockError();
|
|
||||||
}
|
|
||||||
return Pointer(new ParPreconditioner(*seq_precond, comm),
|
|
||||||
Deleter(*seq_precond));
|
|
||||||
}
|
|
||||||
#endif
|
|
||||||
|
|
||||||
/// \brief Solve the system using the given preconditioner and scalar product.
|
|
||||||
template <class Operator, class ScalarProd, class Precond>
|
|
||||||
void solve(Operator& opA, Vector& x, Vector& istlb, ScalarProd& sp, Precond& precond, Dune::InverseOperatorResult& result) const
|
|
||||||
{
|
|
||||||
// TODO: Revise when linear solvers interface opm-core is done
|
|
||||||
// Construct linear solver.
|
|
||||||
// GMRes solver
|
|
||||||
if ( newton_use_gmres_ ) {
|
|
||||||
Dune::RestartedGMResSolver<Vector> linsolve(opA, sp, precond,
|
|
||||||
linear_solver_reduction_, linear_solver_restart_, linear_solver_maxiter_, linear_solver_verbosity_);
|
|
||||||
// Solve system.
|
|
||||||
linsolve.apply(x, istlb, result);
|
|
||||||
}
|
|
||||||
else { // BiCGstab solver
|
|
||||||
Dune::BiCGSTABSolver<Vector> linsolve(opA, sp, precond,
|
|
||||||
linear_solver_reduction_, linear_solver_maxiter_, linear_solver_verbosity_);
|
|
||||||
// Solve system.
|
|
||||||
linsolve.apply(x, istlb, result);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
void formInterleavedSystem(const std::vector<LinearisedBlackoilResidual::ADB>& eqs,
|
|
||||||
Mat& istlA) const;
|
|
||||||
|
|
||||||
mutable int iterations_;
|
|
||||||
boost::any parallelInformation_;
|
boost::any parallelInformation_;
|
||||||
|
mutable int iterations_;
|
||||||
const bool newton_use_gmres_;
|
|
||||||
const double linear_solver_reduction_;
|
|
||||||
const int linear_solver_maxiter_;
|
|
||||||
const int linear_solver_restart_;
|
|
||||||
const int linear_solver_verbosity_;
|
|
||||||
};
|
};
|
||||||
|
|
||||||
} // namespace Opm
|
} // namespace Opm
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user