mirror of
https://github.com/OPM/opm-simulators.git
synced 2026-09-05 04:40:19 -05:00
Merge pull request #2848 from atgeirr/unify-linsolve-only
Unify linear solvers
This commit is contained in:
Markus Blatt
@@ -32,6 +32,7 @@ list (APPEND MAIN_SOURCE_FILES
|
||||
opm/simulators/linalg/FlexibleSolver2.cpp
|
||||
opm/simulators/linalg/FlexibleSolver3.cpp
|
||||
opm/simulators/linalg/FlexibleSolver4.cpp
|
||||
opm/simulators/linalg/setupPropertyTree.cpp
|
||||
opm/simulators/utils/readDeck.cpp
|
||||
opm/simulators/timestepping/TimeStepControl.cpp
|
||||
opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp
|
||||
@@ -170,10 +171,8 @@ list (APPEND PUBLIC_HEADER_FILES
|
||||
opm/simulators/linalg/bda/openclSolverBackend.hpp
|
||||
opm/simulators/linalg/bda/MultisegmentWellContribution.hpp
|
||||
opm/simulators/linalg/bda/WellContributions.hpp
|
||||
opm/simulators/linalg/BlackoilAmg.hpp
|
||||
opm/simulators/linalg/amgcpr.hh
|
||||
opm/simulators/linalg/twolevelmethodcpr.hh
|
||||
opm/simulators/linalg/CPRPreconditioner.hpp
|
||||
opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp
|
||||
opm/simulators/linalg/FlexibleSolver.hpp
|
||||
opm/simulators/linalg/FlexibleSolver_impl.hpp
|
||||
|
||||
@@ -50,14 +50,6 @@ namespace Opm {
|
||||
static constexpr int value = 100;
|
||||
};
|
||||
template<class TypeTag>
|
||||
struct UseAmg<TypeTag, TTag::EclFlowProblemSimple> { // probably not used
|
||||
static constexpr bool value = true;
|
||||
};
|
||||
template<class TypeTag>
|
||||
struct UseCpr<TypeTag, TTag::EclFlowProblemSimple> {
|
||||
static constexpr bool value = true;
|
||||
};
|
||||
template<class TypeTag>
|
||||
struct CprMaxEllIter<TypeTag,TTag::EclFlowProblemSimple> {
|
||||
static constexpr int value = 1;
|
||||
};
|
||||
@@ -70,17 +62,9 @@ namespace Opm {
|
||||
static constexpr int value = 3;
|
||||
};
|
||||
template<class TypeTag>
|
||||
struct CprSolverVerbose<TypeTag, TTag::EclFlowProblemSimple> {
|
||||
static constexpr int value = 0;
|
||||
};
|
||||
template<class TypeTag>
|
||||
struct LinearSolverConfiguration<TypeTag, TTag::EclFlowProblemSimple> {
|
||||
struct Linsolver<TypeTag, TTag::EclFlowProblemSimple> {
|
||||
static constexpr auto value = "ilu0";
|
||||
};
|
||||
template<class TypeTag>
|
||||
struct SystemStrategy<TypeTag, TTag::EclFlowProblemSimple> {
|
||||
static constexpr auto value = "quasiimpes";
|
||||
};
|
||||
|
||||
template<class TypeTag>
|
||||
struct FluidSystem<TypeTag, TTag::EclFlowProblemSimple>
|
||||
|
||||
@@ -16,9 +16,6 @@
|
||||
*/
|
||||
#include "config.h"
|
||||
|
||||
// Define making clear that the simulator supports AMG
|
||||
#define FLOW_SUPPORT_AMG 1
|
||||
|
||||
#include <flow/flow_ebos_blackoil.hpp>
|
||||
|
||||
#include <opm/material/common/ResetLocale.hpp>
|
||||
|
||||
@@ -16,9 +16,6 @@
|
||||
*/
|
||||
#include "config.h"
|
||||
|
||||
// Define making clear that the simulator supports AMG
|
||||
#define FLOW_SUPPORT_AMG 1
|
||||
|
||||
#include <flow/flow_ebos_gasoil.hpp>
|
||||
|
||||
#include <opm/material/common/ResetLocale.hpp>
|
||||
|
||||
@@ -16,9 +16,6 @@
|
||||
*/
|
||||
#include "config.h"
|
||||
|
||||
// Define making clear that the simulator supports AMG
|
||||
#define FLOW_SUPPORT_AMG 1
|
||||
|
||||
#include <flow/flow_ebos_oilwater.hpp>
|
||||
|
||||
#include <opm/material/common/ResetLocale.hpp>
|
||||
|
||||
@@ -16,9 +16,6 @@
|
||||
*/
|
||||
#include "config.h"
|
||||
|
||||
// Define making clear that the simulator supports AMG
|
||||
#define FLOW_SUPPORT_AMG 1
|
||||
|
||||
#include <flow/flow_ebos_oilwater_brine.hpp>
|
||||
|
||||
#include <opm/material/common/ResetLocale.hpp>
|
||||
|
||||
@@ -16,9 +16,6 @@
|
||||
*/
|
||||
#include "config.h"
|
||||
|
||||
// Define making clear that the simulator supports AMG
|
||||
#define FLOW_SUPPORT_AMG 1
|
||||
|
||||
#include <flow/flow_ebos_oilwater_polymer.hpp>
|
||||
|
||||
#include <opm/material/common/ResetLocale.hpp>
|
||||
|
||||
@@ -16,9 +16,6 @@
|
||||
*/
|
||||
#include "config.h"
|
||||
|
||||
// Define making clear that the simulator supports AMG
|
||||
#define FLOW_SUPPORT_AMG 1
|
||||
|
||||
#include <flow/flow_ebos_oilwater_polymer_injectivity.hpp>
|
||||
|
||||
#include <opm/material/common/ResetLocale.hpp>
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@@ -1,528 +0,0 @@
|
||||
/*
|
||||
Copyright 2014 SINTEF ICT, Applied Mathematics.
|
||||
Copyright 2014 IRIS AS.
|
||||
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
|
||||
Copyright 2015 NTNU
|
||||
Copyright 2015 Statoil AS
|
||||
|
||||
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_CPRPRECONDITIONER_HEADER_INCLUDED
|
||||
#define OPM_CPRPRECONDITIONER_HEADER_INCLUDED
|
||||
|
||||
#include <memory>
|
||||
#include <type_traits>
|
||||
|
||||
#include <opm/common/utility/platform_dependent/disable_warnings.h>
|
||||
#include <opm/common/utility/parameters/ParameterGroup.hpp>
|
||||
#include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
|
||||
#include <opm/simulators/linalg/FlowLinearSolverParameters.hpp>
|
||||
|
||||
#include <dune/istl/bvector.hh>
|
||||
#include <dune/istl/bcrsmatrix.hh>
|
||||
#include <dune/istl/operators.hh>
|
||||
#include <dune/istl/io.hh>
|
||||
#include <dune/istl/owneroverlapcopy.hh>
|
||||
#include <dune/istl/preconditioners.hh>
|
||||
#include <dune/istl/schwarz.hh>
|
||||
#include <dune/istl/solvers.hh>
|
||||
#include <dune/istl/paamg/amg.hh>
|
||||
#include <dune/istl/paamg/kamg.hh>
|
||||
#include <dune/istl/paamg/pinfo.hh>
|
||||
|
||||
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
|
||||
|
||||
#include <opm/common/ErrorMacros.hpp>
|
||||
#include <opm/common/Exceptions.hpp>
|
||||
|
||||
#include <opm/simulators/linalg/ParallelRestrictedAdditiveSchwarz.hpp>
|
||||
#include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
template<typename O, typename S, typename C,
|
||||
typename P, std::size_t COMPONENT_INDEX, std::size_t VARIABLE_INDEX>
|
||||
class BlackoilAmg;
|
||||
|
||||
namespace Amg
|
||||
{
|
||||
template<int Row, int Column>
|
||||
class Element;
|
||||
}
|
||||
|
||||
namespace ISTLUtility
|
||||
{
|
||||
|
||||
template<class T>
|
||||
void setILUParameters(Opm::ParallelOverlappingILU0Args<T>& args,
|
||||
const CPRParameter& params)
|
||||
{
|
||||
args.setN(params.cpr_ilu_n_);
|
||||
args.setMilu(params.cpr_ilu_milu_);
|
||||
}
|
||||
|
||||
template<class T>
|
||||
void setILUParameters(Opm::ParallelOverlappingILU0Args<T>& args,
|
||||
MILU_VARIANT milu, int n=0)
|
||||
{
|
||||
args.setN(n);
|
||||
args.setMilu(milu);
|
||||
}
|
||||
|
||||
template<class S, class P>
|
||||
void setILUParameters(S&, const P&)
|
||||
{}
|
||||
|
||||
template<class S, class P>
|
||||
void setILUParameters(S&, bool, int)
|
||||
{}
|
||||
|
||||
///
|
||||
/// \brief A traits class for selecting the types of the preconditioner.
|
||||
///
|
||||
/// \tparam M The type of the matrix.
|
||||
/// \tparam X The type of the domain of the linear problem.
|
||||
/// \tparam Y The type of the range of the linear problem.
|
||||
/// \tparam P The type of the parallel information.
|
||||
////
|
||||
template<class M, class X, class Y, class P>
|
||||
struct CPRSelector
|
||||
{
|
||||
/// \brief The information about the parallelization and communication
|
||||
using ParallelInformation = P;
|
||||
/// \brief The operator type;
|
||||
typedef Dune::OverlappingSchwarzOperator<M,X,X,ParallelInformation> Operator;
|
||||
typedef ParallelOverlappingILU0<M,X, X, ParallelInformation>
|
||||
EllipticPreconditioner;
|
||||
/// \brief The type of the unique pointer to the preconditioner of the elliptic part.
|
||||
typedef std::unique_ptr<EllipticPreconditioner>
|
||||
EllipticPreconditionerPointer;
|
||||
|
||||
typedef EllipticPreconditioner Smoother;
|
||||
typedef Dune::Amg::AMG<Operator, X, Smoother, ParallelInformation> AMG;
|
||||
|
||||
/// \brief creates an Operator from the matrix
|
||||
/// \param M The matrix to use.
|
||||
/// \param p The parallel information to use.
|
||||
static Operator* makeOperator(const M& m, const ParallelInformation& p)
|
||||
{
|
||||
return new Operator(m, p);
|
||||
}
|
||||
};
|
||||
|
||||
template<class M, class X, class Y>
|
||||
struct CPRSelector<M,X,Y,Dune::Amg::SequentialInformation>
|
||||
{
|
||||
/// \brief The information about the parallelization and communication
|
||||
typedef Dune::Amg::SequentialInformation ParallelInformation;
|
||||
/// \brief The operator type;
|
||||
typedef Dune::MatrixAdapter<M, X, Y> Operator;
|
||||
/// \brief The type of the preconditioner used for the elliptic part.
|
||||
typedef ParallelOverlappingILU0<M,X, X, ParallelInformation>
|
||||
EllipticPreconditioner;
|
||||
/// \brief The type of the unique pointer to the preconditioner of the elliptic part.
|
||||
typedef std::unique_ptr<EllipticPreconditioner> EllipticPreconditionerPointer;
|
||||
|
||||
/// \brief type of AMG used to precondition the elliptic system.
|
||||
typedef EllipticPreconditioner Smoother;
|
||||
typedef Dune::Amg::AMG<Operator, X, Smoother, ParallelInformation> AMG;
|
||||
|
||||
/// \brief creates an Operator from the matrix
|
||||
/// \param M The matrix to use.
|
||||
/// \param p The parallel information to use.
|
||||
static Operator* makeOperator(const M& m, const ParallelInformation&)
|
||||
{
|
||||
return new Operator(m);
|
||||
}
|
||||
|
||||
};
|
||||
//! \brief Creates and initializes a unique pointer to an sequential ILU0 preconditioner.
|
||||
//! \param A The matrix of the linear system to solve.
|
||||
//! \param relax The relaxation factor to use.
|
||||
template<class M, class X, class C>
|
||||
std::shared_ptr<ParallelOverlappingILU0<M,X,X,C> >
|
||||
createILU0Ptr(const M& A, const C& comm, double relax, MILU_VARIANT milu)
|
||||
{
|
||||
return std::make_shared<ParallelOverlappingILU0<M,X,X,C> >(A, comm, relax, milu);
|
||||
}
|
||||
//! \brief Creates and initializes a shared pointer to an ILUn preconditioner.
|
||||
//! \param A The matrix of the linear system to solve.
|
||||
//! \param ilu_n The n parameter for the extension of the nonzero pattern.
|
||||
//! \param relax The relaxation factor to use.
|
||||
template<class M, class X, class C>
|
||||
std::shared_ptr<ParallelOverlappingILU0<M,X,X,C> >
|
||||
createILUnPtr(const M& A, const C& comm, int ilu_n, double relax, MILU_VARIANT milu)
|
||||
{
|
||||
return std::make_shared<ParallelOverlappingILU0<M,X,X,C> >( A, comm, ilu_n, relax, milu );
|
||||
}
|
||||
/// \brief Creates the elliptic preconditioner (ILU0)
|
||||
/// \param Ae The matrix of the elliptic system.
|
||||
/// \param relax The relaxation parameter for ILU0.
|
||||
/// \param milu If true, the modified ilu approach is used. Dropped elements
|
||||
/// will get added to the diagonal of U to preserve the row sum
|
||||
/// for constant vectors (Ae = LUe).
|
||||
/// \param comm The object describing the parallelization information and communication.
|
||||
template<class M, class X=typename M::range_type, class P>
|
||||
typename CPRSelector<M,X,X,P>::EllipticPreconditionerPointer
|
||||
createEllipticPreconditionerPointer(const M& Ae, double relax,
|
||||
MILU_VARIANT milu, const P& comm)
|
||||
{
|
||||
typedef typename CPRSelector<M,X,X,P >
|
||||
::EllipticPreconditioner ParallelPreconditioner;
|
||||
|
||||
typedef typename CPRSelector<M,X,X,P>
|
||||
::EllipticPreconditionerPointer EllipticPreconditionerPointer;
|
||||
return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax, milu));
|
||||
}
|
||||
|
||||
template < class C, class Op, class P, class S, std::size_t PressureEqnIndex, std::size_t PressureVarIndex, class Vector>
|
||||
inline void
|
||||
createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm,
|
||||
std::unique_ptr< BlackoilAmg<Op,S,C,P,PressureEqnIndex,PressureVarIndex> >& amgPtr,
|
||||
const CPRParameter& params,
|
||||
const Vector& weights)
|
||||
{
|
||||
using AMG = BlackoilAmg<Op,S,C,P,PressureEqnIndex,PressureVarIndex>;
|
||||
int verbosity = 0;
|
||||
if (comm.communicator().rank() == 0) {
|
||||
verbosity = params.cpr_solver_verbose_;
|
||||
}
|
||||
// TODO: revise choice of parameters
|
||||
int coarsenTarget=1200;
|
||||
using Criterion = C;
|
||||
Criterion criterion(15, coarsenTarget);
|
||||
criterion.setDebugLevel( verbosity ); // no debug information, 1 for printing hierarchy information
|
||||
criterion.setDefaultValuesIsotropic(2);
|
||||
criterion.setNoPostSmoothSteps( 1 );
|
||||
criterion.setNoPreSmoothSteps( 1 );
|
||||
|
||||
// Since DUNE 2.2 we also need to pass the smoother args instead of steps directly
|
||||
typedef typename AMG::Smoother Smoother;
|
||||
typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
|
||||
SmootherArgs smootherArgs;
|
||||
smootherArgs.iterations = 1;
|
||||
smootherArgs.relaxationFactor = relax;
|
||||
setILUParameters(smootherArgs, params);
|
||||
|
||||
amgPtr.reset( new AMG( params, weights, opA, criterion, smootherArgs, comm ) );
|
||||
}
|
||||
|
||||
template < class C, class Op, class P, class AMG >
|
||||
inline void
|
||||
createAMGPreconditionerPointer(Op& opA, const double relax, const MILU_VARIANT milu, const P& comm, std::unique_ptr< AMG >& amgPtr)
|
||||
{
|
||||
// TODO: revise choice of parameters
|
||||
int coarsenTarget=1200;
|
||||
using Criterion = C;
|
||||
Criterion criterion(15, coarsenTarget);
|
||||
criterion.setDebugLevel( 0 ); // no debug information, 1 for printing hierarchy information
|
||||
criterion.setDefaultValuesIsotropic(2);
|
||||
criterion.setNoPostSmoothSteps( 1 );
|
||||
criterion.setNoPreSmoothSteps( 1 );
|
||||
|
||||
// for DUNE 2.2 we also need to pass the smoother args
|
||||
typedef typename AMG::Smoother Smoother;
|
||||
typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
|
||||
SmootherArgs smootherArgs;
|
||||
smootherArgs.iterations = 1;
|
||||
smootherArgs.relaxationFactor = relax;
|
||||
setILUParameters(smootherArgs, milu);
|
||||
|
||||
amgPtr.reset( new AMG(opA, criterion, smootherArgs, comm ) );
|
||||
}
|
||||
|
||||
/// \brief Creates the elliptic preconditioner (ILU0)
|
||||
/// \param opA The operator representing the matrix of the system.
|
||||
/// \param relax The relaxation parameter for ILU0.
|
||||
/// \param comm The object describing the parallelization information and communication.
|
||||
// \param amgPtr The unique_ptr to be filled (return)
|
||||
template < int PressureEqnIndex, int PressureVarIndex, class Op, class P, class AMG >
|
||||
inline void
|
||||
createAMGPreconditionerPointer( Op& opA, const double relax, const MILU_VARIANT milu, const P& comm, std::unique_ptr< AMG >& amgPtr )
|
||||
{
|
||||
// type of matrix
|
||||
typedef typename Op::matrix_type M;
|
||||
|
||||
// The coupling metric used in the AMG
|
||||
typedef Opm::Amg::Element<PressureEqnIndex, PressureVarIndex> CouplingMetric;
|
||||
|
||||
// The coupling criterion used in the AMG
|
||||
typedef Dune::Amg::SymmetricCriterion<M, CouplingMetric> CritBase;
|
||||
|
||||
// The coarsening criterion used in the AMG
|
||||
typedef Dune::Amg::CoarsenCriterion<CritBase> Criterion;
|
||||
|
||||
createAMGPreconditionerPointer<Criterion>(opA, relax, milu, comm, amgPtr);
|
||||
}
|
||||
|
||||
} // end namespace ISTLUtility
|
||||
|
||||
/*!
|
||||
\brief CPR preconditioner.
|
||||
|
||||
This is a two-stage preconditioner, combining an elliptic-type
|
||||
partial solution with ILU0 for the whole system.
|
||||
|
||||
\tparam M The matrix type to operate on
|
||||
\tparam X Type of the update
|
||||
\tparam Y Type of the defect
|
||||
\tparam P Type of the parallel information. If not provided
|
||||
this will be Dune::Amg::SequentialInformation.
|
||||
The preconditioner is parallel if this is
|
||||
Dune::OwnerOverlapCopyCommunication<int,int>
|
||||
*/
|
||||
template<class M, class X, class Y,
|
||||
class P=Dune::Amg::SequentialInformation>
|
||||
class CPRPreconditioner : public Dune::Preconditioner<X,Y>
|
||||
{
|
||||
// prohibit copying for now
|
||||
CPRPreconditioner( const CPRPreconditioner& );
|
||||
|
||||
public:
|
||||
//! \brief The type describing the parallel information
|
||||
typedef P ParallelInformation;
|
||||
//! \brief The matrix type the preconditioner is for.
|
||||
typedef typename std::remove_const<M>::type matrix_type;
|
||||
//! \brief The domain type of the preconditioner.
|
||||
typedef X domain_type;
|
||||
//! \brief The range type of the preconditioner.
|
||||
typedef Y range_type;
|
||||
//! \brief The field type of the preconditioner.
|
||||
typedef typename X::field_type field_type;
|
||||
|
||||
// define the category
|
||||
Dune::SolverCategory::Category category() const override
|
||||
{
|
||||
return std::is_same<P,Dune::Amg::SequentialInformation>::value ?
|
||||
Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping;
|
||||
}
|
||||
|
||||
typedef ISTLUtility::CPRSelector<M,X,X,P> CPRSelectorType ;
|
||||
|
||||
//! \brief Elliptic Operator
|
||||
typedef typename CPRSelectorType::Operator Operator;
|
||||
|
||||
//! \brief preconditioner for the whole system (here either ILU(0) or ILU(n)
|
||||
typedef Dune::Preconditioner<X,X> WholeSystemPreconditioner;
|
||||
|
||||
//! \brief type of the unique pointer to the ilu-0 preconditioner
|
||||
//! used the for the elliptic system
|
||||
typedef typename CPRSelectorType::EllipticPreconditionerPointer
|
||||
EllipticPreconditionerPointer;
|
||||
|
||||
//! \brief amg preconditioner for the elliptic system
|
||||
typedef typename CPRSelectorType::AMG AMG;
|
||||
|
||||
/*! \brief Constructor.
|
||||
|
||||
Constructor gets all parameters to operate the prec.
|
||||
\param A The matrix to operate on.
|
||||
\param Ae The top-left elliptic part of A.
|
||||
\param relax The ILU0 relaxation factor.
|
||||
\param useAMG if true, AMG is used as a preconditioner for the elliptic sub-system, otherwise ilu-0 (default)
|
||||
\param useBiCG if true, BiCG solver is used (default), otherwise CG solver
|
||||
\param paralleInformation The information about the parallelization, if this is a
|
||||
parallel run
|
||||
*/
|
||||
CPRPreconditioner (const CPRParameter& param, const M& A, const M& Ae,
|
||||
const ParallelInformation& comm=ParallelInformation(),
|
||||
const ParallelInformation& commAe=ParallelInformation())
|
||||
: param_( param ),
|
||||
A_(A),
|
||||
Ae_(Ae),
|
||||
de_( Ae_.N() ),
|
||||
ve_( Ae_.M() ),
|
||||
dmodified_( A_.N() ),
|
||||
opAe_(CPRSelectorType::makeOperator(Ae_, commAe)),
|
||||
precond_(), // ilu0 preconditioner for elliptic system
|
||||
amg_(), // amg preconditioner for elliptic system
|
||||
pre_(), // copy A will be made be the preconditioner
|
||||
vilu_( A_.N() ),
|
||||
comm_(comm),
|
||||
commAe_(commAe)
|
||||
{
|
||||
// create appropriate preconditioner for elliptic system
|
||||
createEllipticPreconditioner( param_.cpr_use_amg_, commAe_ );
|
||||
|
||||
// create the preconditioner for the whole system.
|
||||
if( param_.cpr_ilu_n_ == 0 ) {
|
||||
pre_ = ISTLUtility::createILU0Ptr<M,X>( A_, comm, param_.cpr_relax_, param_.cpr_ilu_milu_ );
|
||||
}
|
||||
else {
|
||||
pre_ = ISTLUtility::createILUnPtr<M,X>( A_, comm, param_.cpr_ilu_n_, param_.cpr_relax_, param_.cpr_ilu_milu_);
|
||||
}
|
||||
}
|
||||
|
||||
/*!
|
||||
\brief Prepare the preconditioner.
|
||||
|
||||
\copydoc Preconditioner::pre(X&,Y&)
|
||||
*/
|
||||
virtual void pre (X& /*x*/, Y& /*b*/)
|
||||
{
|
||||
}
|
||||
|
||||
/*!
|
||||
\brief Apply the preconditoner.
|
||||
|
||||
\copydoc Preconditioner::apply(X&,const Y&)
|
||||
*/
|
||||
virtual void apply (X& v, const Y& d)
|
||||
{
|
||||
// Extract part of d corresponding to elliptic part.
|
||||
// Note: Assumes that the elliptic part comes first.
|
||||
std::copy_n(d.begin(), de_.size(), de_.begin());
|
||||
|
||||
// Solve elliptic part, extend solution to full.
|
||||
// reset result
|
||||
ve_ = 0;
|
||||
solveElliptic( ve_, de_ );
|
||||
|
||||
//reset return value
|
||||
v = 0.0;
|
||||
// Again assuming that the elliptic part comes first.
|
||||
std::copy(ve_.begin(), ve_.end(), v.begin());
|
||||
|
||||
// Subtract elliptic residual from initial residual.
|
||||
// dmodified = d - A * vfull
|
||||
dmodified_ = d;
|
||||
A_.mmv(v, dmodified_);
|
||||
// A is not parallel, do communication manually.
|
||||
comm_.copyOwnerToAll(dmodified_, dmodified_);
|
||||
|
||||
// Apply Preconditioner for whole system (relax will be applied already)
|
||||
pre_->apply( vilu_, dmodified_);
|
||||
|
||||
// don't apply relaxation if relax_ == 1
|
||||
if( std::abs( param_.cpr_relax_ - 1.0 ) < 1e-12 ) {
|
||||
v += vilu_;
|
||||
}
|
||||
else {
|
||||
v *= param_.cpr_relax_;
|
||||
v += vilu_;
|
||||
}
|
||||
}
|
||||
|
||||
/*!
|
||||
\brief Clean up.
|
||||
|
||||
\copydoc Preconditioner::post(X&)
|
||||
*/
|
||||
virtual void post (X& /*x*/)
|
||||
{
|
||||
}
|
||||
|
||||
protected:
|
||||
void solveElliptic(Y& x, Y& de)
|
||||
{
|
||||
// Linear solver parameters
|
||||
const double tolerance = param_.cpr_solver_tol_;
|
||||
const int maxit = param_.cpr_max_ell_iter_;
|
||||
const int verbosity = ( param_.cpr_solver_verbose_ &&
|
||||
comm_.communicator().rank()==0 ) ? 1 : 0;
|
||||
|
||||
// operator result containing iterations etc.
|
||||
Dune::InverseOperatorResult result;
|
||||
|
||||
// the scalar product chooser
|
||||
auto sp = Dune::createScalarProduct<X,ParallelInformation>(commAe_, category());
|
||||
|
||||
if( amg_ )
|
||||
{
|
||||
// Solve system with AMG
|
||||
if( param_.cpr_use_bicgstab_ ) {
|
||||
Dune::BiCGSTABSolver<X> linsolve(*opAe_, *sp, (*amg_), tolerance, maxit, verbosity);
|
||||
linsolve.apply(x, de, result);
|
||||
}
|
||||
else {
|
||||
Dune::CGSolver<X> linsolve(*opAe_, *sp, (*amg_), tolerance, maxit, verbosity);
|
||||
linsolve.apply(x, de, result);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
assert( precond_ );
|
||||
// Solve system with ILU-0
|
||||
if( param_.cpr_use_bicgstab_ ) {
|
||||
Dune::BiCGSTABSolver<X> linsolve(*opAe_, *sp, (*precond_), tolerance, maxit, verbosity);
|
||||
linsolve.apply(x, de, result);
|
||||
}
|
||||
else {
|
||||
Dune::CGSolver<X> linsolve(*opAe_, *sp, (*precond_), tolerance, maxit, verbosity);
|
||||
linsolve.apply(x, de, result);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
if (!result.converged) {
|
||||
OPM_THROW(LinearSolverProblem, "CPRPreconditioner failed to solve elliptic subsystem.");
|
||||
}
|
||||
}
|
||||
|
||||
//! \brief Parameter collection for CPR
|
||||
const CPRParameter& param_;
|
||||
|
||||
//! \brief The matrix for the full linear problem.
|
||||
const matrix_type& A_;
|
||||
//! \brief The elliptic part of the matrix.
|
||||
const matrix_type& Ae_;
|
||||
|
||||
//! \brief temporary variables for elliptic solve
|
||||
Y de_, ve_, dmodified_;
|
||||
|
||||
//! \brief elliptic operator
|
||||
std::unique_ptr<Operator> opAe_;
|
||||
|
||||
//! \brief ILU0 preconditioner for the elliptic system
|
||||
EllipticPreconditionerPointer precond_;
|
||||
//! \brief AMG preconditioner with ILU0 smoother
|
||||
std::unique_ptr< AMG > amg_;
|
||||
|
||||
//! \brief The preconditioner for the whole system
|
||||
//!
|
||||
//! We have to use a shared_ptr instead of a unique_ptr
|
||||
//! as we need to use a custom allocator based on dynamic
|
||||
//! information. But for unique_ptr the type of this deleter
|
||||
//! has to be available at coompile time.
|
||||
std::shared_ptr< WholeSystemPreconditioner > pre_;
|
||||
|
||||
//! \brief temporary variables for ILU solve
|
||||
Y vilu_;
|
||||
|
||||
//! \brief The information about the parallelization of the whole system.
|
||||
const P& comm_;
|
||||
//! \brief The information about the parallelization of the elliptic part
|
||||
//! of the system
|
||||
const P& commAe_;
|
||||
protected:
|
||||
void createEllipticPreconditioner( const bool amg, const P& comm )
|
||||
{
|
||||
if( amg )
|
||||
{
|
||||
ISTLUtility::createAMGPreconditionerPointer( *opAe_ , param_.cpr_relax_, param_.cpr_ilu_milu_, comm, amg_ );
|
||||
}
|
||||
else
|
||||
{
|
||||
precond_ = ISTLUtility::createEllipticPreconditionerPointer<M,X>( Ae_, param_.cpr_relax_, param_.cpr_ilu_milu_, comm);
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif // OPM_CPRPRECONDITIONER_HEADER_INCLUDED
|
||||
@@ -1,5 +1,5 @@
|
||||
/*
|
||||
Copyright 2015 SINTEF ICT, Applied Mathematics.
|
||||
Copyright 2015, 2020 SINTEF Digital, Mathematics and Cybernetics.
|
||||
Copyright 2015 IRIS AS
|
||||
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
|
||||
Copyright 2015 NTNU
|
||||
@@ -94,34 +94,14 @@ struct LinearSolverIgnoreConvergenceFailure{
|
||||
using type = UndefinedProperty;
|
||||
};
|
||||
template<class TypeTag, class MyTypeTag>
|
||||
struct UseAmg {
|
||||
using type = UndefinedProperty;
|
||||
};
|
||||
template<class TypeTag, class MyTypeTag>
|
||||
struct UseCpr {
|
||||
using type = UndefinedProperty;
|
||||
};
|
||||
template<class TypeTag, class MyTypeTag>
|
||||
struct PreconditionerAddWellContributions {
|
||||
using type = UndefinedProperty;
|
||||
};
|
||||
template<class TypeTag, class MyTypeTag>
|
||||
struct SystemStrategy {
|
||||
using type = UndefinedProperty;
|
||||
};
|
||||
template<class TypeTag, class MyTypeTag>
|
||||
struct ScaleLinearSystem {
|
||||
using type = UndefinedProperty;
|
||||
};
|
||||
template<class TypeTag, class MyTypeTag>
|
||||
struct CprSolverVerbose {
|
||||
using type = UndefinedProperty;
|
||||
};
|
||||
template<class TypeTag, class MyTypeTag>
|
||||
struct CprUseDrs {
|
||||
using type = UndefinedProperty;
|
||||
};
|
||||
template<class TypeTag, class MyTypeTag>
|
||||
struct CprMaxEllIter {
|
||||
using type = UndefinedProperty;
|
||||
};
|
||||
@@ -134,11 +114,7 @@ struct CprReuseSetup {
|
||||
using type = UndefinedProperty;
|
||||
};
|
||||
template<class TypeTag, class MyTypeTag>
|
||||
struct LinearSolverConfiguration {
|
||||
using type = UndefinedProperty;
|
||||
};
|
||||
template<class TypeTag, class MyTypeTag>
|
||||
struct LinearSolverConfigurationJsonFile {
|
||||
struct Linsolver {
|
||||
using type = UndefinedProperty;
|
||||
};
|
||||
template<class TypeTag, class MyTypeTag>
|
||||
@@ -205,14 +181,6 @@ struct LinearSolverIgnoreConvergenceFailure<TypeTag, TTag::FlowIstlSolverParams>
|
||||
static constexpr bool value = false;
|
||||
};
|
||||
template<class TypeTag>
|
||||
struct UseAmg<TypeTag, TTag::FlowIstlSolverParams> {
|
||||
static constexpr bool value = false;
|
||||
};
|
||||
template<class TypeTag>
|
||||
struct UseCpr<TypeTag, TTag::FlowIstlSolverParams> {
|
||||
static constexpr bool value = false;
|
||||
};
|
||||
template<class TypeTag>
|
||||
struct LinearSolverBackend<TypeTag, TTag::FlowIstlSolverParams> {
|
||||
using type = Opm::ISTLSolverEbos<TypeTag>;
|
||||
};
|
||||
@@ -221,22 +189,10 @@ struct PreconditionerAddWellContributions<TypeTag, TTag::FlowIstlSolverParams> {
|
||||
static constexpr bool value = false;
|
||||
};
|
||||
template<class TypeTag>
|
||||
struct SystemStrategy<TypeTag, TTag::FlowIstlSolverParams> {
|
||||
static constexpr auto value = "none";
|
||||
};
|
||||
template<class TypeTag>
|
||||
struct ScaleLinearSystem<TypeTag, TTag::FlowIstlSolverParams> {
|
||||
static constexpr bool value = false;
|
||||
};
|
||||
template<class TypeTag>
|
||||
struct CprSolverVerbose<TypeTag, TTag::FlowIstlSolverParams> {
|
||||
static constexpr int value = 0;
|
||||
};
|
||||
template<class TypeTag>
|
||||
struct CprUseDrs<TypeTag, TTag::FlowIstlSolverParams> {
|
||||
static constexpr bool value = false;
|
||||
};
|
||||
template<class TypeTag>
|
||||
struct CprMaxEllIter<TypeTag, TTag::FlowIstlSolverParams> {
|
||||
static constexpr int value = 20;
|
||||
};
|
||||
@@ -249,14 +205,10 @@ struct CprReuseSetup<TypeTag, TTag::FlowIstlSolverParams> {
|
||||
static constexpr int value = 3;
|
||||
};
|
||||
template<class TypeTag>
|
||||
struct LinearSolverConfiguration<TypeTag, TTag::FlowIstlSolverParams> {
|
||||
struct Linsolver<TypeTag, TTag::FlowIstlSolverParams> {
|
||||
static constexpr auto value = "ilu0";
|
||||
};
|
||||
template<class TypeTag>
|
||||
struct LinearSolverConfigurationJsonFile<TypeTag, TTag::FlowIstlSolverParams> {
|
||||
static constexpr auto value = "none";
|
||||
};
|
||||
template<class TypeTag>
|
||||
struct GpuMode<TypeTag, TTag::FlowIstlSolverParams> {
|
||||
static constexpr auto value = "none";
|
||||
};
|
||||
@@ -274,51 +226,9 @@ struct OpenclPlatformId<TypeTag, TTag::FlowIstlSolverParams> {
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
/**
|
||||
* \brief Parameters used to configure the CPRPreconditioner.
|
||||
*/
|
||||
struct CPRParameter
|
||||
{
|
||||
double cpr_relax_;
|
||||
double cpr_solver_tol_;
|
||||
int cpr_ilu_n_;
|
||||
MILU_VARIANT cpr_ilu_milu_;
|
||||
bool cpr_ilu_redblack_;
|
||||
bool cpr_ilu_reorder_sphere_;
|
||||
bool cpr_use_drs_;
|
||||
int cpr_max_ell_iter_;
|
||||
int cpr_ell_solvetype_;
|
||||
bool cpr_use_amg_;
|
||||
bool cpr_use_bicgstab_;
|
||||
int cpr_solver_verbose_;
|
||||
bool cpr_pressure_aggregation_;
|
||||
int cpr_reuse_setup_;
|
||||
CPRParameter() { reset(); }
|
||||
|
||||
void reset()
|
||||
{
|
||||
cpr_solver_tol_ = 1e-2;
|
||||
cpr_ilu_n_ = 0;
|
||||
cpr_ilu_milu_ = MILU_VARIANT::ILU;
|
||||
cpr_ilu_redblack_ = false;
|
||||
cpr_ilu_reorder_sphere_ = true;
|
||||
cpr_max_ell_iter_ = 25;
|
||||
cpr_ell_solvetype_ = 0;
|
||||
cpr_use_drs_ = false;
|
||||
cpr_use_amg_ = true;
|
||||
cpr_use_bicgstab_ = true;
|
||||
cpr_solver_verbose_ = 0;
|
||||
cpr_pressure_aggregation_ = false;
|
||||
cpr_reuse_setup_ = 0;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
/// This class carries all parameters for the NewtonIterationBlackoilInterleaved class.
|
||||
struct FlowLinearSolverParameters
|
||||
: public CPRParameter
|
||||
{
|
||||
double linear_solver_reduction_;
|
||||
double ilu_relaxation_;
|
||||
@@ -332,15 +242,14 @@ namespace Opm
|
||||
bool newton_use_gmres_;
|
||||
bool require_full_sparsity_pattern_;
|
||||
bool ignoreConvergenceFailure_;
|
||||
bool linear_solver_use_amg_;
|
||||
bool use_cpr_;
|
||||
std::string system_strategy_;
|
||||
bool scale_linear_system_;
|
||||
std::string linear_solver_configuration_;
|
||||
std::string linear_solver_configuration_json_file_;
|
||||
std::string linsolver_;
|
||||
std::string gpu_mode_;
|
||||
int bda_device_id_;
|
||||
int opencl_platform_id_;
|
||||
int cpr_max_ell_iter_ = 20;
|
||||
int cpr_reuse_setup_ = 0;
|
||||
bool use_gpu_;
|
||||
|
||||
template <class TypeTag>
|
||||
void init()
|
||||
@@ -358,17 +267,10 @@ namespace Opm
|
||||
newton_use_gmres_ = EWOMS_GET_PARAM(TypeTag, bool, UseGmres);
|
||||
require_full_sparsity_pattern_ = EWOMS_GET_PARAM(TypeTag, bool, LinearSolverRequireFullSparsityPattern);
|
||||
ignoreConvergenceFailure_ = EWOMS_GET_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure);
|
||||
linear_solver_use_amg_ = EWOMS_GET_PARAM(TypeTag, bool, UseAmg);
|
||||
use_cpr_ = EWOMS_GET_PARAM(TypeTag, bool, UseCpr);
|
||||
system_strategy_ = EWOMS_GET_PARAM(TypeTag, std::string, SystemStrategy);
|
||||
scale_linear_system_ = EWOMS_GET_PARAM(TypeTag, bool, ScaleLinearSystem);
|
||||
cpr_solver_verbose_ = EWOMS_GET_PARAM(TypeTag, int, CprSolverVerbose);
|
||||
cpr_use_drs_ = EWOMS_GET_PARAM(TypeTag, bool, CprUseDrs);
|
||||
cpr_max_ell_iter_ = EWOMS_GET_PARAM(TypeTag, int, CprMaxEllIter);
|
||||
cpr_ell_solvetype_ = EWOMS_GET_PARAM(TypeTag, int, CprEllSolvetype);
|
||||
cpr_reuse_setup_ = EWOMS_GET_PARAM(TypeTag, int, CprReuseSetup);
|
||||
linear_solver_configuration_ = EWOMS_GET_PARAM(TypeTag, std::string, LinearSolverConfiguration);
|
||||
linear_solver_configuration_json_file_ = EWOMS_GET_PARAM(TypeTag, std::string, LinearSolverConfigurationJsonFile);
|
||||
linsolver_ = EWOMS_GET_PARAM(TypeTag, std::string, Linsolver);
|
||||
gpu_mode_ = EWOMS_GET_PARAM(TypeTag, std::string, GpuMode);
|
||||
bda_device_id_ = EWOMS_GET_PARAM(TypeTag, int, BdaDeviceId);
|
||||
opencl_platform_id_ = EWOMS_GET_PARAM(TypeTag, int, OpenclPlatformId);
|
||||
@@ -389,17 +291,10 @@ namespace Opm
|
||||
EWOMS_REGISTER_PARAM(TypeTag, bool, UseGmres, "Use GMRES as the linear solver");
|
||||
EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverRequireFullSparsityPattern, "Produce the full sparsity pattern for the linear solver");
|
||||
EWOMS_REGISTER_PARAM(TypeTag, bool, LinearSolverIgnoreConvergenceFailure, "Continue with the simulation like nothing happened after the linear solver did not converge");
|
||||
EWOMS_REGISTER_PARAM(TypeTag, bool, UseAmg, "Use AMG as the linear solver's preconditioner");
|
||||
EWOMS_REGISTER_PARAM(TypeTag, bool, UseCpr, "Use CPR as the linear solver's preconditioner");
|
||||
EWOMS_REGISTER_PARAM(TypeTag, std::string, SystemStrategy, "Strategy for reformulating and scaling linear system (none: no scaling -- should not be used with CPR, original: use weights that are equivalent to no scaling -- should not be used with CPR, simple: form pressure equation as simple sum of conservation equations, quasiimpes: form pressure equation based on diagonal block, trueimpes: form pressure equation based on linearization of accumulation term)");
|
||||
EWOMS_REGISTER_PARAM(TypeTag, bool, ScaleLinearSystem, "Scale linear system according to equation scale and primary variable types");
|
||||
EWOMS_REGISTER_PARAM(TypeTag, int, CprSolverVerbose, "Verbosity of cpr solver (0: silent, 1: print summary of inner linear solver, 2: print extensive information about inner linear solve, including setup information)");
|
||||
EWOMS_REGISTER_PARAM(TypeTag, bool, CprUseDrs, "Use dynamic row sum using weights");
|
||||
EWOMS_REGISTER_PARAM(TypeTag, int, CprMaxEllIter, "MaxIterations of the elliptic pressure part of the cpr solver");
|
||||
EWOMS_REGISTER_PARAM(TypeTag, int, CprEllSolvetype, "Solver type of elliptic pressure solve (0: bicgstab, 1: cg, 2: only amg preconditioner)");
|
||||
EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseSetup, "Reuse Amg Setup");
|
||||
EWOMS_REGISTER_PARAM(TypeTag, std::string, LinearSolverConfiguration, "Configuration of solver valid is: ilu0 (default), cpr_quasiimpes, cpr_trueimpes or file (specified in LinearSolverConfigurationJsonFile) ");
|
||||
EWOMS_REGISTER_PARAM(TypeTag, std::string, LinearSolverConfigurationJsonFile, "Filename of JSON configuration for flexible linear solver system.");
|
||||
EWOMS_REGISTER_PARAM(TypeTag, int, CprReuseSetup, "Reuse preconditioner setup. Valid options are 0: recreate the preconditioner for every linear solve, 1: recreate once every timestep, 2: recreate if last linear solve took more than 10 iterations, 3: never recreate");
|
||||
EWOMS_REGISTER_PARAM(TypeTag, std::string, Linsolver, "Configuration of solver. Valid options are: ilu0 (default), cpr (an alias for cpr_trueimpes), cpr_quasiimpes, cpr_trueimpes or amg. Alternatively, you can request a configuration to be read from a JSON file by giving the filename here, ending with '.json.'");
|
||||
EWOMS_REGISTER_PARAM(TypeTag, std::string, GpuMode, "Use GPU cusparseSolver or openclSolver as the linear solver, usage: '--gpu-mode=[none|cusparse|opencl]'");
|
||||
EWOMS_REGISTER_PARAM(TypeTag, int, BdaDeviceId, "Choose device ID for cusparseSolver or openclSolver, use 'nvidia-smi' or 'clinfo' to determine valid IDs");
|
||||
EWOMS_REGISTER_PARAM(TypeTag, int, OpenclPlatformId, "Choose platform ID for openclSolver, use 'clinfo' to determine valid platform IDs");
|
||||
@@ -410,7 +305,6 @@ namespace Opm
|
||||
// set default values
|
||||
void reset()
|
||||
{
|
||||
use_cpr_ = false;
|
||||
newton_use_gmres_ = false;
|
||||
linear_solver_reduction_ = 1e-2;
|
||||
linear_solver_maxiter_ = 150;
|
||||
@@ -418,7 +312,6 @@ namespace Opm
|
||||
linear_solver_verbosity_ = 0;
|
||||
require_full_sparsity_pattern_ = false;
|
||||
ignoreConvergenceFailure_ = false;
|
||||
linear_solver_use_amg_ = false;
|
||||
ilu_fillin_level_ = 0;
|
||||
ilu_relaxation_ = 0.9;
|
||||
ilu_milu_ = MILU_VARIANT::ILU;
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,108 @@
|
||||
/*
|
||||
Copyright 2019, 2020 SINTEF Digital, Mathematics and Cybernetics.
|
||||
|
||||
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 <opm/simulators/linalg/setupPropertyTree.hpp>
|
||||
|
||||
#include <boost/property_tree/json_parser.hpp>
|
||||
#include <boost/version.hpp>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
boost::property_tree::ptree
|
||||
setupCPR(const std::string& conf, const FlowLinearSolverParameters& p)
|
||||
{
|
||||
boost::property_tree::ptree prm;
|
||||
prm.put("maxiter", p.linear_solver_maxiter_);
|
||||
prm.put("tol", p.linear_solver_reduction_);
|
||||
prm.put("verbosity", p.linear_solver_verbosity_);
|
||||
prm.put("solver", "bicgstab");
|
||||
prm.put("preconditioner.type", "cpr");
|
||||
if (conf == "cpr_quasiimpes") {
|
||||
prm.put("preconditioner.weight_type", "quasiimpes");
|
||||
} else {
|
||||
prm.put("preconditioner.weight_type", "trueimpes");
|
||||
}
|
||||
prm.put("preconditioner.finesmoother.type", "ParOverILU0");
|
||||
prm.put("preconditioner.finesmoother.relaxation", 1.0);
|
||||
prm.put("preconditioner.pressure_var_index", 1);
|
||||
prm.put("preconditioner.verbosity", 0);
|
||||
prm.put("preconditioner.coarsesolver.maxiter", 1);
|
||||
prm.put("preconditioner.coarsesolver.tol", 1e-1);
|
||||
prm.put("preconditioner.coarsesolver.solver", "loopsolver");
|
||||
prm.put("preconditioner.coarsesolver.verbosity", 0);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.type", "amg");
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.alpha", 0.333333333333);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.relaxation", 1.0);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.iterations", p.cpr_max_ell_iter_);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.coarsenTarget", 1200);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.pre_smooth", 1);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.post_smooth", 1);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.beta", 1e-5);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.smoother", "ILU0");
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.verbosity", 0);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.maxlevel", 15);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.skip_isolated", 0);
|
||||
return prm;
|
||||
}
|
||||
|
||||
|
||||
boost::property_tree::ptree
|
||||
setupAMG([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p)
|
||||
{
|
||||
boost::property_tree::ptree prm;
|
||||
prm.put("tol", p.linear_solver_reduction_);
|
||||
prm.put("maxiter", p.linear_solver_maxiter_);
|
||||
prm.put("verbosity", p.linear_solver_verbosity_);
|
||||
prm.put("solver", "bicgstab");
|
||||
prm.put("preconditioner.type", "amg");
|
||||
prm.put("preconditioner.alpha", 0.333333333333);
|
||||
prm.put("preconditioner.relaxation", 1.0);
|
||||
prm.put("preconditioner.iterations", 20);
|
||||
prm.put("preconditioner.coarsenTarget", 1200);
|
||||
prm.put("preconditioner.pre_smooth", 1);
|
||||
prm.put("preconditioner.post_smooth", 1);
|
||||
prm.put("preconditioner.beta", 1e-5);
|
||||
prm.put("preconditioner.smoother", "ILU0");
|
||||
prm.put("preconditioner.verbosity", 0);
|
||||
prm.put("preconditioner.maxlevel", 15);
|
||||
prm.put("preconditioner.skip_isolated", 0);
|
||||
return prm;
|
||||
}
|
||||
|
||||
|
||||
boost::property_tree::ptree
|
||||
setupILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p)
|
||||
{
|
||||
boost::property_tree::ptree prm;
|
||||
prm.put("tol", p.linear_solver_reduction_);
|
||||
prm.put("maxiter", p.linear_solver_maxiter_);
|
||||
prm.put("verbosity", p.linear_solver_verbosity_);
|
||||
prm.put("solver", "bicgstab");
|
||||
prm.put("preconditioner.type", "ParOverILU0");
|
||||
prm.put("preconditioner.relaxation", p.ilu_relaxation_);
|
||||
prm.put("preconditioner.ilulevel", p.ilu_fillin_level_);
|
||||
return prm;
|
||||
}
|
||||
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
@@ -28,7 +28,11 @@ namespace Opm
|
||||
{
|
||||
|
||||
template<class TypeTag>
|
||||
boost::property_tree::ptree setupPropertyTree(const FlowLinearSolverParameters& p);
|
||||
boost::property_tree::ptree setupPropertyTree(FlowLinearSolverParameters p);
|
||||
|
||||
boost::property_tree::ptree setupCPR(const std::string& conf, const FlowLinearSolverParameters& p);
|
||||
boost::property_tree::ptree setupAMG(const std::string& conf, const FlowLinearSolverParameters& p);
|
||||
boost::property_tree::ptree setupILU(const std::string& conf, const FlowLinearSolverParameters& p);
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
|
||||
@@ -34,89 +34,60 @@ namespace Opm
|
||||
/// from file the data in the JSON file will override any other options.
|
||||
template<class TypeTag>
|
||||
boost::property_tree::ptree
|
||||
setupPropertyTree(const FlowLinearSolverParameters& p)
|
||||
setupPropertyTree(FlowLinearSolverParameters p) // Note: copying the parameters to potentially override.
|
||||
{
|
||||
boost::property_tree::ptree prm;
|
||||
if (p.linear_solver_configuration_ == "file") {
|
||||
std::string conf = p.linsolver_;
|
||||
|
||||
// Get configuration from file.
|
||||
if (conf.size() > 5 && conf.substr(conf.size() - 5, 5) == ".json") { // the ends_with() method is not available until C++20
|
||||
#if BOOST_VERSION / 100 % 1000 > 48
|
||||
if (p.linear_solver_configuration_json_file_ == "none"){
|
||||
OPM_THROW(std::invalid_argument,
|
||||
"--linear-solver-configuration=file requires that a filename "
|
||||
<< "is passed with "
|
||||
<< "--linear-solver-configuration-json-file=filename.");
|
||||
}else{
|
||||
boost::property_tree::read_json(p.linear_solver_configuration_json_file_, prm);
|
||||
try {
|
||||
boost::property_tree::ptree prm;
|
||||
boost::property_tree::read_json(conf, prm);
|
||||
return prm;
|
||||
}
|
||||
catch (...) {
|
||||
OPM_THROW(std::invalid_argument, "Failed reading linear solver configuration from JSON file " << conf);
|
||||
}
|
||||
#else
|
||||
OPM_THROW(std::invalid_argument,
|
||||
"--linear-solver-configuration=file not supported with "
|
||||
<< "boost version. Needs versoin > 1.48.");
|
||||
OPM_THROW(std::invalid_argument,
|
||||
"--linear-solver-configuration=file.json not supported with "
|
||||
<< "boost version. Needs version > 1.48.");
|
||||
#endif
|
||||
}
|
||||
else
|
||||
{
|
||||
std::string conf = p.linear_solver_configuration_;
|
||||
// Support old UseCpr if not configuration was set
|
||||
if (!EWOMS_PARAM_IS_SET(TypeTag, std::string, LinearSolverConfiguration) && p.use_cpr_)
|
||||
{
|
||||
|
||||
// Use CPR configuration.
|
||||
if ((conf == "cpr") || (conf == "cpr_trueimpes") || (conf == "cpr_quasiimpes")) {
|
||||
if (conf == "cpr") {
|
||||
// Treat "cpr" as short cut for the true IMPES variant.
|
||||
conf = "cpr_trueimpes";
|
||||
}
|
||||
|
||||
if((conf == "cpr_trueimpes") || (conf == "cpr_quasiimpes")){
|
||||
prm.put("tol", p.linear_solver_reduction_);
|
||||
if (EWOMS_PARAM_IS_SET(TypeTag, int, LinearSolverMaxIter))
|
||||
prm.put("maxiter", p.linear_solver_maxiter_);// Trust that the user knows what he does
|
||||
else
|
||||
prm.put("maxiter", 20); // Use our own default.
|
||||
prm.put("verbosity", p.linear_solver_verbosity_);
|
||||
prm.put("solver", "bicgstab");
|
||||
prm.put("preconditioner.type", "cpr");
|
||||
if (conf == "cpr_quasiimpes" )
|
||||
{
|
||||
prm.put("preconditioner.weight_type","quasiimpes");
|
||||
}
|
||||
else
|
||||
{
|
||||
prm.put("preconditioner.weight_type","trueimpes");
|
||||
}
|
||||
prm.put("preconditioner.finesmoother.type", "ParOverILU0");
|
||||
prm.put("preconditioner.finesmoother.relaxation", 1.0);
|
||||
prm.put("preconditioner.pressure_var_index",1);
|
||||
prm.put("preconditioner.verbosity",0);
|
||||
prm.put("preconditioner.coarsesolver.maxiter",1);
|
||||
prm.put("preconditioner.coarsesolver.tol",1e-1);
|
||||
prm.put("preconditioner.coarsesolver.solver","loopsolver");
|
||||
prm.put("preconditioner.coarsesolver.verbosity",0);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.type","amg");
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.alpha",0.333333333333);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.relaxation",1.0);
|
||||
if (EWOMS_PARAM_IS_SET(TypeTag, int, CprMaxEllIter))
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.iterations", p.cpr_max_ell_iter_);
|
||||
else
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.iterations",1);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.coarsenTarget",1200);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.pre_smooth",1);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.post_smooth",1);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.beta",1e-5);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.smoother","ILU0");
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.verbosity",0);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.maxlevel",15);
|
||||
prm.put("preconditioner.coarsesolver.preconditioner.skip_isolated",0);
|
||||
} else {
|
||||
if(conf != "ilu0"){
|
||||
OPM_THROW(std::invalid_argument, conf << "is not a valid setting for --linear-solver-configuration."
|
||||
<< " Please use ilu0, cpr_trueimpes, or cpr_quasiimpes");
|
||||
}
|
||||
prm.put("tol", p.linear_solver_reduction_);
|
||||
prm.put("maxiter", p.linear_solver_maxiter_);
|
||||
prm.put("verbosity", p.linear_solver_verbosity_);
|
||||
prm.put("solver", "bicgstab");
|
||||
prm.put("preconditioner.type", "ParOverILU0");
|
||||
prm.put("preconditioner.relaxation", p.ilu_relaxation_);
|
||||
prm.put("preconditioner.ilulevel", p.ilu_fillin_level_);
|
||||
if (!EWOMS_PARAM_IS_SET(TypeTag, int, LinearSolverMaxIter)) {
|
||||
// Use our own default unless it was explicitly overridden by user.
|
||||
p.linear_solver_maxiter_ = 20;
|
||||
}
|
||||
if (!EWOMS_PARAM_IS_SET(TypeTag, int, CprMaxEllIter)) {
|
||||
// Use our own default unless it was explicitly overridden by user.
|
||||
p.cpr_max_ell_iter_ = 1;
|
||||
}
|
||||
return setupCPR(conf, p);
|
||||
}
|
||||
return prm;
|
||||
|
||||
if (conf == "amg") {
|
||||
return setupAMG(conf, p);
|
||||
}
|
||||
|
||||
// Use ILU0 configuration.
|
||||
if (conf == "ilu0") {
|
||||
return setupILU(conf, p);
|
||||
}
|
||||
|
||||
// No valid configuration option found.
|
||||
OPM_THROW(std::invalid_argument,
|
||||
conf << " is not a valid setting for --linear-solver-configuration."
|
||||
<< " Please use ilu0, cpr, cpr_trueimpes, or cpr_quasiimpes");
|
||||
}
|
||||
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
+13
-21
@@ -21,13 +21,15 @@
|
||||
#define BOOST_TEST_MODULE BlackoilAmgTest
|
||||
#define BOOST_TEST_NO_MAIN
|
||||
#include <boost/test/unit_test.hpp>
|
||||
#include <opm/simulators/linalg/BlackoilAmg.hpp>
|
||||
#include <opm/simulators/linalg/PreconditionerFactory.hpp>
|
||||
#include <opm/simulators/linalg/FlexibleSolver.hpp>
|
||||
|
||||
#include <dune/common/parallel/mpihelper.hh>
|
||||
#include <dune/common/fmatrix.hh>
|
||||
#include <dune/common/unused.hh>
|
||||
#include <dune/common/parallel/indexset.hh>
|
||||
#include <dune/common/parallel/plocalindex.hh>
|
||||
#include <dune/common/version.hh>
|
||||
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 7)
|
||||
#include <dune/common/parallel/communication.hh>
|
||||
#else
|
||||
@@ -35,6 +37,7 @@
|
||||
#endif
|
||||
#include <dune/istl/bcrsmatrix.hh>
|
||||
#include <dune/istl/owneroverlapcopy.hh>
|
||||
#include <dune/istl/schwarz.hh>
|
||||
|
||||
class MPIError {
|
||||
public:
|
||||
@@ -260,6 +263,7 @@ M setupAnisotropic2d(int N, Dune::ParallelIndexSet<G,L,s>& indices, const Dune::
|
||||
return mat;
|
||||
}
|
||||
|
||||
|
||||
//BOOST_AUTO_TEST_CASE(runBlackoilAmgLaplace)
|
||||
void runBlackoilAmgLaplace()
|
||||
{
|
||||
@@ -289,30 +293,18 @@ void runBlackoilAmgLaplace()
|
||||
setBoundary(x, b, N, comm.indexSet());
|
||||
|
||||
Operator fop(mat, comm);
|
||||
typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<BCRSMat,Dune::Amg::FirstDiagonal> >
|
||||
Criterion;
|
||||
typedef Dune::SeqSSOR<BCRSMat,Vector,Vector> Smoother;
|
||||
//typedef Dune::SeqJac<BCRSMat,Vector,Vector> Smoother;
|
||||
//typedef Dune::SeqILU0<BCRSMat,Vector,Vector> Smoother;
|
||||
//typedef Dune::SeqILUn<BCRSMat,Vector,Vector> Smoother;
|
||||
typedef Dune::BlockPreconditioner<Vector,Vector,Communication,Smoother> ParSmoother;
|
||||
typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
|
||||
|
||||
Dune::OverlappingSchwarzScalarProduct<Vector,Communication> sp(comm);
|
||||
|
||||
Dune::InverseOperatorResult r;
|
||||
SmootherArgs smootherArgs;
|
||||
Criterion criterion;
|
||||
|
||||
smootherArgs.iterations = 1;
|
||||
Opm::CPRParameter param;
|
||||
boost::property_tree::ptree prm;
|
||||
prm.put("type", "amg");
|
||||
std::function<Vector()> weights = [&mat]() {
|
||||
return Opm::Amg::getQuasiImpesWeights<BCRSMat, Vector>(mat, 0, false);
|
||||
};
|
||||
auto amg = Opm::PreconditionerFactory<Operator, Communication>::create(fop, prm, weights, comm);
|
||||
|
||||
Dune::CGSolver<Vector> amgCG(fop, sp, *amg, 10e-8, 300, (ccomm.rank()==0) ? 2 : 0);
|
||||
|
||||
Opm::BlackoilAmg<Operator,ParSmoother,Criterion,Communication,0,0> amg(param,
|
||||
{},
|
||||
fop, criterion,
|
||||
smootherArgs,
|
||||
comm);
|
||||
Dune::CGSolver<Vector> amgCG(fop, sp, amg, 10e-8, 300, (ccomm.rank()==0) ? 2 : 0);
|
||||
amgCG.apply(x,b,r);
|
||||
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user