mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Merge pull request #1397 from blattms/system-amg-rebased
Add support for CPR (and variants) for flow_ebos
This commit is contained in:
Atgeirr Flø Rasmussen
GitHub
commit
8070e9d229
@ -138,6 +138,7 @@ endif()
|
||||
list (APPEND TEST_SOURCE_FILES
|
||||
tests/test_autodiffhelpers.cpp
|
||||
tests/test_autodiffmatrix.cpp
|
||||
tests/test_blackoil_amg.cpp
|
||||
tests/test_block.cpp
|
||||
tests/test_boprops_ad.cpp
|
||||
tests/test_rateconverter.cpp
|
||||
@ -259,11 +260,11 @@ list (APPEND PROGRAM_SOURCE_FILES
|
||||
# originally generated with the command:
|
||||
# find opm -name '*.h*' -a ! -name '*-pch.hpp' -printf '\t%p\n' | sort
|
||||
list (APPEND PUBLIC_HEADER_FILES
|
||||
opm/autodiff/AdditionalObjectDeleter.hpp
|
||||
opm/autodiff/AutoDiffBlock.hpp
|
||||
opm/autodiff/AutoDiffHelpers.hpp
|
||||
opm/autodiff/AutoDiffMatrix.hpp
|
||||
opm/autodiff/AutoDiff.hpp
|
||||
opm/autodiff/BlackoilAmg.hpp
|
||||
opm/autodiff/BlackoilDetails.hpp
|
||||
opm/autodiff/BlackoilLegacyDetails.hpp
|
||||
opm/autodiff/BlackoilModel.hpp
|
||||
|
||||
@ -1,58 +0,0 @@
|
||||
/*
|
||||
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
|
||||
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_ADDITIONALOBJECTDELETER_HEADER_INCLUDED
|
||||
#define OPM_ADDITIONALOBJECTDELETER_HEADER_INCLUDED
|
||||
namespace Opm
|
||||
{
|
||||
//! \brief A custom deleter that will delete an additional object, too.
|
||||
//!
|
||||
//! In dune-istl most parallel preconditioners hold a reference to
|
||||
//! a sequential preconditioner.
|
||||
//! In CPRPreconditioner and NewtonIterationBlackoilInterleaved we use unique_ptr
|
||||
//! for the memory management.
|
||||
//! Ergo we need to construct the sequential preconditioner with new and
|
||||
//! make sure that it gets deleted together with the enclosing parallel
|
||||
//! preconditioner. Therefore this deleter stores a pointer to it and deletes
|
||||
//! it during destruction.
|
||||
//! \tparam The type of the additional object to be deleted.
|
||||
template<class T>
|
||||
class AdditionalObjectDeleter
|
||||
{
|
||||
public:
|
||||
//! \brief empty constructor.
|
||||
AdditionalObjectDeleter()
|
||||
: additional_object_()
|
||||
{}
|
||||
//! \brief Constructor taking the object that needs to deleted.
|
||||
AdditionalObjectDeleter(T& additional_object)
|
||||
: additional_object_(&additional_object){}
|
||||
//! \brief Delete an object and the additional one.
|
||||
template<class T1>
|
||||
void operator()(T1* pt)
|
||||
{
|
||||
delete pt;
|
||||
delete additional_object_;
|
||||
}
|
||||
private:
|
||||
T* additional_object_;
|
||||
};
|
||||
|
||||
}
|
||||
#endif
|
||||
1005
opm/autodiff/BlackoilAmg.hpp
Normal file
1005
opm/autodiff/BlackoilAmg.hpp
Normal file
File diff suppressed because it is too large
Load Diff
@ -47,12 +47,60 @@
|
||||
#include <opm/common/ErrorMacros.hpp>
|
||||
#include <opm/common/Exceptions.hpp>
|
||||
|
||||
#include <opm/autodiff/AdditionalObjectDeleter.hpp>
|
||||
#include <opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp>
|
||||
#include <opm/autodiff/ParallelOverlappingILU0.hpp>
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
template<typename O, typename S, typename C,
|
||||
typename P, std::size_t COMPONENT_INDEX>
|
||||
class BlackoilAmg;
|
||||
|
||||
/**
|
||||
* \brief Parameters used to configure the CPRPreconditioner.
|
||||
*/
|
||||
struct CPRParameter
|
||||
{
|
||||
double cpr_relax_;
|
||||
double cpr_solver_tol_;
|
||||
int cpr_ilu_n_;
|
||||
int cpr_max_ell_iter_;
|
||||
bool cpr_use_amg_;
|
||||
bool cpr_use_bicgstab_;
|
||||
bool cpr_solver_verbose_;
|
||||
bool cpr_pressure_aggregation_;
|
||||
|
||||
CPRParameter() { reset(); }
|
||||
|
||||
CPRParameter( const ParameterGroup& param)
|
||||
{
|
||||
// reset values to default
|
||||
reset();
|
||||
|
||||
cpr_relax_ = param.getDefault("cpr_relax", cpr_relax_);
|
||||
cpr_solver_tol_ = param.getDefault("cpr_solver_tol", cpr_solver_tol_);
|
||||
cpr_ilu_n_ = param.getDefault("cpr_ilu_n", cpr_ilu_n_);
|
||||
cpr_max_ell_iter_ = param.getDefault("cpr_max_elliptic_iter",cpr_max_ell_iter_);
|
||||
cpr_use_amg_ = param.getDefault("cpr_use_amg", cpr_use_amg_);
|
||||
cpr_use_bicgstab_ = param.getDefault("cpr_use_bicgstab", cpr_use_bicgstab_);
|
||||
cpr_solver_verbose_ = param.getDefault("cpr_solver_verbose", cpr_solver_verbose_);
|
||||
cpr_pressure_aggregation_ = param.getDefault("cpr_pressure_aggregation", cpr_pressure_aggregation_);
|
||||
}
|
||||
|
||||
void reset()
|
||||
{
|
||||
cpr_relax_ = 1.0;
|
||||
cpr_solver_tol_ = 1e-2;
|
||||
cpr_ilu_n_ = 0;
|
||||
cpr_max_ell_iter_ = 25;
|
||||
cpr_use_amg_ = true;
|
||||
cpr_use_bicgstab_ = true;
|
||||
cpr_solver_verbose_ = false;
|
||||
cpr_pressure_aggregation_ = false;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
namespace ISTLUtility
|
||||
{
|
||||
///
|
||||
@ -67,47 +115,13 @@ template<class M, class X, class Y, class P>
|
||||
struct CPRSelector
|
||||
{
|
||||
/// \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 Dune::SeqILU0<M,X,X> 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 P&)
|
||||
{
|
||||
return new Operator(m);
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
#if HAVE_MPI
|
||||
/// \copydoc CPRSelector<M,X,X,Y,P>
|
||||
template<class M, class X, class Y, class I1, class I2>
|
||||
struct CPRSelector<M,X,Y,Dune::OwnerOverlapCopyCommunication<I1,I2> >
|
||||
{
|
||||
/// \brief The information about the parallelization and communication
|
||||
typedef Dune::OwnerOverlapCopyCommunication<I1,I2> ParallelInformation;
|
||||
using ParallelInformation = P;
|
||||
/// \brief The operator type;
|
||||
typedef Dune::OverlappingSchwarzOperator<M,X,X,ParallelInformation> Operator;
|
||||
/// \brief The type of the preconditioner used for the elliptic part.
|
||||
//typedef Dune::BlockPreconditioner<X, X, ParallelInformation, Dune::SeqILU0<M,X,X> >
|
||||
//EllipticPreconditioner;
|
||||
//typedef ParallelRestrictedOverlappingSchwarz<X, X, ParallelInformation, Dune::SeqILU0<M,X,X> >
|
||||
//EllipticPreconditioner;
|
||||
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,
|
||||
AdditionalObjectDeleter<Dune::SeqILU0<M,X,X> > >
|
||||
typedef std::unique_ptr<EllipticPreconditioner>
|
||||
EllipticPreconditionerPointer;
|
||||
|
||||
typedef EllipticPreconditioner Smoother;
|
||||
@ -122,154 +136,116 @@ struct CPRSelector<M,X,Y,Dune::OwnerOverlapCopyCommunication<I1,I2> >
|
||||
}
|
||||
};
|
||||
|
||||
//! \brief Creates the deleter needed for the parallel ILU preconditioners.
|
||||
//! \tparam ILU The type of the underlying sequential ILU preconditioner.
|
||||
//! \tparam I1 The global index type.
|
||||
//! \tparam I2 The local index type.
|
||||
//! \param ilu A reference to the wrapped preconditioner
|
||||
//! \param p The parallel information for template parameter deduction.
|
||||
template<class ILU, class I1, class I2>
|
||||
AdditionalObjectDeleter<ILU>
|
||||
createParallelDeleter(ILU& ilu, const Dune::OwnerOverlapCopyCommunication<I1,I2>& p)
|
||||
{
|
||||
(void) p;
|
||||
return AdditionalObjectDeleter<ILU>(ilu);
|
||||
}
|
||||
#endif
|
||||
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>
|
||||
std::shared_ptr<Dune::SeqILU0<M,X,X> >
|
||||
createILU0Ptr(const M& A, double relax, const Dune::Amg::SequentialInformation&)
|
||||
template<class M, class X, class C>
|
||||
std::shared_ptr<ParallelOverlappingILU0<M,X,X,C> >
|
||||
createILU0Ptr(const M& A, const C& comm, double relax)
|
||||
{
|
||||
return std::shared_ptr<Dune::SeqILU0<M,X,X> >(new Dune::SeqILU0<M,X,X>( A, relax) );
|
||||
return std::make_shared<ParallelOverlappingILU0<M,X,X,C> >(A, comm, relax);
|
||||
}
|
||||
//! \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>
|
||||
std::shared_ptr<Dune::SeqILUn<M,X,X> >
|
||||
createILUnPtr(const M& A, int ilu_n, double relax, const Dune::Amg::SequentialInformation&)
|
||||
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)
|
||||
{
|
||||
return std::shared_ptr<Dune::SeqILUn<M,X,X> >(new Dune::SeqILUn<M,X,X>( A, ilu_n, relax) );
|
||||
return std::make_shared<ParallelOverlappingILU0<M,X,X,C> >( A, comm, ilu_n, relax);
|
||||
}
|
||||
|
||||
#if HAVE_MPI
|
||||
template<class ILU, class I1, class I2>
|
||||
struct SelectParallelILUSharedPtr
|
||||
{
|
||||
typedef std::shared_ptr<
|
||||
Dune::BlockPreconditioner<
|
||||
typename ILU::range_type,
|
||||
typename ILU::domain_type,
|
||||
Dune::OwnerOverlapCopyCommunication<I1,I2>,
|
||||
ILU
|
||||
>
|
||||
> type;
|
||||
};
|
||||
|
||||
//! \brief Creates and initializes a shared pointer to an ILUn preconditioner.
|
||||
//! \param A The matrix of the linear system to solve.
|
||||
//! \param relax The relaxation factor to use.
|
||||
/// \param comm The object describing the parallelization information and communication.
|
||||
template<class M, class X, class I1, class I2>
|
||||
typename SelectParallelILUSharedPtr<Dune::SeqILU0<M,X,X>, I1, I2>::type
|
||||
createILU0Ptr(const M& A, double relax,
|
||||
const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
|
||||
{
|
||||
typedef Dune::BlockPreconditioner<
|
||||
X,
|
||||
X,
|
||||
Dune::OwnerOverlapCopyCommunication<I1,I2>,
|
||||
Dune::SeqILU0<M,X,X>
|
||||
> PointerType;
|
||||
Dune::SeqILU0<M,X,X>* ilu = nullptr;
|
||||
int ilu_setup_successful = 1;
|
||||
std::string message;
|
||||
try {
|
||||
ilu = new Dune::SeqILU0<M,X,X>(A, 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 ( ilu ) // not null if constructor succeeded
|
||||
{
|
||||
// prevent memory leak
|
||||
delete ilu;
|
||||
}
|
||||
throw Dune::MatrixBlockError();
|
||||
}
|
||||
return typename SelectParallelILUSharedPtr<Dune::SeqILU0<M,X,X>, I1, I2>
|
||||
::type ( new PointerType(*ilu, comm), createParallelDeleter(*ilu, comm));
|
||||
}
|
||||
|
||||
//! \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.
|
||||
/// \param comm The object describing the parallelization information and communication.
|
||||
template<class M, class X, class I1, class I2>
|
||||
typename SelectParallelILUSharedPtr<Dune::SeqILUn<M,X,X>, I1, I2>::type
|
||||
createILUnPtr(const M& A, int ilu_n, double relax,
|
||||
const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
|
||||
{
|
||||
typedef Dune::BlockPreconditioner<
|
||||
X,
|
||||
X,
|
||||
Dune::OwnerOverlapCopyCommunication<I1,I2>,
|
||||
Dune::SeqILUn<M,X,X>
|
||||
> PointerType;
|
||||
Dune::SeqILUn<M,X,X>* ilu = new Dune::SeqILUn<M,X,X>( A, ilu_n, relax);
|
||||
|
||||
return typename SelectParallelILUSharedPtr<Dune::SeqILUn<M,X,X>, I1, I2>::type
|
||||
(new PointerType(*ilu, comm),createParallelDeleter(*ilu, comm));
|
||||
}
|
||||
#endif
|
||||
|
||||
/// \brief Creates the elliptic preconditioner (ILU0)
|
||||
/// \param Ae The matrix of the elliptic system.
|
||||
/// \param relax The relaxation parameter for ILU0
|
||||
template<class M, class X=typename M::range_type>
|
||||
std::unique_ptr<Dune::SeqILU0<M,X,X> >
|
||||
createEllipticPreconditionerPointer(const M& Ae, double relax,
|
||||
const Dune::Amg::SequentialInformation&)
|
||||
{
|
||||
return std::unique_ptr<Dune::SeqILU0<M,X,X> >(new Dune::SeqILU0<M,X,X>(Ae, relax));
|
||||
}
|
||||
|
||||
#if HAVE_MPI
|
||||
/// \brief Creates the elliptic preconditioner (ILU0)
|
||||
/// \param Ae The matrix of the elliptic system.
|
||||
/// \param relax The relaxation parameter for ILU0.
|
||||
/// \param comm The object describing the parallelization information and communication.
|
||||
template<class M, class X=typename M::range_type, class I1, class I2>
|
||||
typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
|
||||
::EllipticPreconditionerPointer
|
||||
template<class M, class X=typename M::range_type, class P>
|
||||
typename CPRSelector<M,X,X,P>::EllipticPreconditionerPointer
|
||||
createEllipticPreconditionerPointer(const M& Ae, double relax,
|
||||
const Dune::OwnerOverlapCopyCommunication<I1,I2>& comm)
|
||||
const P& comm)
|
||||
{
|
||||
typedef typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
|
||||
typedef typename CPRSelector<M,X,X,P >
|
||||
::EllipticPreconditioner ParallelPreconditioner;
|
||||
|
||||
//Dune::SeqILU0<M,X,X>* ilu=new Dune::SeqILU0<M,X,X>(Ae, relax);
|
||||
typedef typename CPRSelector<M,X,X,Dune::OwnerOverlapCopyCommunication<I1,I2> >
|
||||
typedef typename CPRSelector<M,X,X,P>
|
||||
::EllipticPreconditionerPointer EllipticPreconditionerPointer;
|
||||
return EllipticPreconditionerPointer(new ParallelPreconditioner(Ae, comm, relax));
|
||||
}
|
||||
|
||||
#endif // #if HAVE_MPI
|
||||
template < class C, class Op, class P, class S, std::size_t index >
|
||||
inline void
|
||||
createAMGPreconditionerPointer(Op& opA, const double relax, const P& comm,
|
||||
std::unique_ptr< BlackoilAmg<Op,S,C,P,index> >& amgPtr,
|
||||
const CPRParameter& params)
|
||||
{
|
||||
using AMG = BlackoilAmg<Op,S,C,P,index>;
|
||||
// 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 );
|
||||
|
||||
// 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;
|
||||
|
||||
amgPtr.reset( new AMG( params, opA, criterion, smootherArgs, comm ) );
|
||||
}
|
||||
|
||||
template < class C, class Op, class P, class AMG >
|
||||
inline void
|
||||
createAMGPreconditionerPointer(Op& opA, const double relax, 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;
|
||||
|
||||
amgPtr.reset( new AMG(opA, criterion, smootherArgs, comm ) );
|
||||
}
|
||||
|
||||
/// \brief Creates the elliptic preconditioner (ILU0)
|
||||
/// \param opA The operator representing the matrix of the system.
|
||||
@ -292,65 +268,11 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std:
|
||||
// The coarsening criterion used in the AMG
|
||||
typedef Dune::Amg::CoarsenCriterion<CritBase> Criterion;
|
||||
|
||||
// TODO: revise choice of parameters
|
||||
int coarsenTarget=1200;
|
||||
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;
|
||||
|
||||
amgPtr.reset( new AMG(opA, criterion, smootherArgs, comm ) );
|
||||
createAMGPreconditionerPointer<Criterion>(opA, relax, comm, amgPtr);
|
||||
}
|
||||
|
||||
} // end namespace ISTLUtility
|
||||
|
||||
struct CPRParameter
|
||||
{
|
||||
double cpr_relax_;
|
||||
double cpr_solver_tol_;
|
||||
int cpr_ilu_n_;
|
||||
int cpr_max_ell_iter_;
|
||||
bool cpr_use_amg_;
|
||||
bool cpr_use_bicgstab_;
|
||||
bool cpr_solver_verbose_;
|
||||
|
||||
CPRParameter() { reset(); }
|
||||
|
||||
CPRParameter( const ParameterGroup& param)
|
||||
{
|
||||
// reset values to default
|
||||
reset();
|
||||
|
||||
cpr_relax_ = param.getDefault("cpr_relax", cpr_relax_);
|
||||
cpr_solver_tol_ = param.getDefault("cpr_solver_tol", cpr_solver_tol_);
|
||||
cpr_ilu_n_ = param.getDefault("cpr_ilu_n", cpr_ilu_n_);
|
||||
cpr_max_ell_iter_ = param.getDefault("cpr_max_elliptic_iter",cpr_max_ell_iter_);
|
||||
cpr_use_amg_ = param.getDefault("cpr_use_amg", cpr_use_amg_);
|
||||
cpr_use_bicgstab_ = param.getDefault("cpr_use_bicgstab", cpr_use_bicgstab_);
|
||||
cpr_solver_verbose_ = param.getDefault("cpr_solver_verbose", cpr_solver_verbose_);
|
||||
}
|
||||
|
||||
void reset()
|
||||
{
|
||||
cpr_relax_ = 1.0;
|
||||
cpr_solver_tol_ = 1e-2;
|
||||
cpr_ilu_n_ = 0;
|
||||
cpr_max_ell_iter_ = 25;
|
||||
cpr_use_amg_ = true;
|
||||
cpr_use_bicgstab_ = true;
|
||||
cpr_solver_verbose_ = false;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
/*!
|
||||
\brief CPR preconditioner.
|
||||
|
||||
@ -448,10 +370,10 @@ createAMGPreconditionerPointer( Op& opA, const double relax, const P& comm, std:
|
||||
|
||||
// create the preconditioner for the whole system.
|
||||
if( param_.cpr_ilu_n_ == 0 ) {
|
||||
pre_ = ISTLUtility::createILU0Ptr<M,X>( A_, param_.cpr_relax_, comm_ );
|
||||
pre_ = ISTLUtility::createILU0Ptr<M,X>( A_, comm, param_.cpr_relax_ );
|
||||
}
|
||||
else {
|
||||
pre_ = ISTLUtility::createILUnPtr<M,X>( A_, param_.cpr_ilu_n_, param_.cpr_relax_, comm_ );
|
||||
pre_ = ISTLUtility::createILUnPtr<M,X>( A_, comm, param_.cpr_ilu_n_, param_.cpr_relax_);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -577,7 +577,14 @@ namespace Opm
|
||||
void setupLinearSolver()
|
||||
{
|
||||
typedef typename BlackoilModelEbos<TypeTag> :: ISTLSolverType ISTLSolverType;
|
||||
|
||||
const std::string cprSolver = "cpr";
|
||||
if (!param_.has("solver_approach") )
|
||||
{
|
||||
if ( eclState().getSimulationConfig().useCPR() )
|
||||
{
|
||||
param_.insertParameter("solver_approach", cprSolver);
|
||||
}
|
||||
}
|
||||
extractParallelGridInformationToISTL(grid(), parallel_information_);
|
||||
fis_solver_.reset( new ISTLSolverType( param_, parallel_information_ ) );
|
||||
}
|
||||
|
||||
@ -20,7 +20,7 @@
|
||||
#ifndef OPM_ISTLSOLVER_HEADER_INCLUDED
|
||||
#define OPM_ISTLSOLVER_HEADER_INCLUDED
|
||||
|
||||
#include <opm/autodiff/AdditionalObjectDeleter.hpp>
|
||||
#include <opm/autodiff/BlackoilAmg.hpp>
|
||||
#include <opm/autodiff/CPRPreconditioner.hpp>
|
||||
#include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
|
||||
#include <opm/autodiff/NewtonIterationUtilities.hpp>
|
||||
@ -425,13 +425,11 @@ namespace Detail
|
||||
parallelInformation_arg.copyOwnerToAll(istlb, istlb);
|
||||
|
||||
#if FLOW_SUPPORT_AMG // activate AMG if either flow_ebos is used or UMFPack is not available
|
||||
if( parameters_.linear_solver_use_amg_ )
|
||||
if( parameters_.linear_solver_use_amg_ || parameters_.use_cpr_)
|
||||
{
|
||||
typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm> CPRSelectorType;
|
||||
typedef typename CPRSelectorType::AMG AMG;
|
||||
typedef typename CPRSelectorType::Operator MatrixOperator;
|
||||
|
||||
std::unique_ptr< AMG > amg;
|
||||
std::unique_ptr< MatrixOperator > opA;
|
||||
|
||||
if( ! std::is_same< LinearOperator, MatrixOperator > :: value )
|
||||
@ -440,13 +438,35 @@ namespace Detail
|
||||
opA.reset( CPRSelectorType::makeOperator( linearOperator.getmat(), parallelInformation_arg ) );
|
||||
}
|
||||
|
||||
const double relax = 1.0;
|
||||
const double relax = parameters_.ilu_relaxation_;
|
||||
if ( parameters_.use_cpr_ )
|
||||
{
|
||||
using Matrix = typename MatrixOperator::matrix_type;
|
||||
using CouplingMetric = Dune::Amg::Diagonal<pressureIndex>;
|
||||
using CritBase = Dune::Amg::SymmetricCriterion<Matrix, CouplingMetric>;
|
||||
using Criterion = Dune::Amg::CoarsenCriterion<CritBase>;
|
||||
using AMG = typename ISTLUtility
|
||||
::BlackoilAmgSelector< Matrix, Vector, Vector,POrComm, Criterion, pressureIndex >::AMG;
|
||||
|
||||
// Construct preconditioner.
|
||||
constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax );
|
||||
std::unique_ptr< AMG > amg;
|
||||
// Construct preconditioner.
|
||||
Criterion crit(15, 2000);
|
||||
constructAMGPrecond<Criterion>( linearOperator, parallelInformation_arg, amg, opA, relax );
|
||||
|
||||
// Solve.
|
||||
solve(linearOperator, x, istlb, *sp, *amg, result);
|
||||
// Solve.
|
||||
solve(linearOperator, x, istlb, *sp, *amg, result);
|
||||
}
|
||||
else
|
||||
{
|
||||
typedef typename CPRSelectorType::AMG AMG;
|
||||
std::unique_ptr< AMG > amg;
|
||||
|
||||
// Construct preconditioner.
|
||||
constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax );
|
||||
|
||||
// Solve.
|
||||
solve(linearOperator, x, istlb, *sp, *amg, result);
|
||||
}
|
||||
}
|
||||
else
|
||||
#endif
|
||||
@ -516,6 +536,20 @@ namespace Detail
|
||||
ISTLUtility::template createAMGPreconditionerPointer<pressureIndex>( opA, relax, comm, amg );
|
||||
}
|
||||
|
||||
template <class C, class LinearOperator, class MatrixOperator, class POrComm, class AMG >
|
||||
void
|
||||
constructAMGPrecond(LinearOperator& /* linearOperator */, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax ) const
|
||||
{
|
||||
ISTLUtility::template createAMGPreconditionerPointer<C>( *opA, relax, comm, amg, parameters_ );
|
||||
}
|
||||
|
||||
|
||||
template <class C, class MatrixOperator, class POrComm, class AMG >
|
||||
void
|
||||
constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax ) const
|
||||
{
|
||||
ISTLUtility::template createAMGPreconditionerPointer<C>( opA, relax, comm, amg, parameters_ );
|
||||
}
|
||||
/// \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
|
||||
|
||||
@ -27,7 +27,6 @@
|
||||
#define FLOW_SUPPORT_AMG !defined(HAVE_UMFPACK)
|
||||
|
||||
#include <opm/autodiff/DuneMatrix.hpp>
|
||||
#include <opm/autodiff/AdditionalObjectDeleter.hpp>
|
||||
#include <opm/autodiff/CPRPreconditioner.hpp>
|
||||
#include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
|
||||
#include <opm/autodiff/NewtonIterationUtilities.hpp>
|
||||
|
||||
@ -25,6 +25,7 @@
|
||||
#ifndef OPM_NEWTONITERATIONBLACKOILINTERLEAVED_HEADER_INCLUDED
|
||||
#define OPM_NEWTONITERATIONBLACKOILINTERLEAVED_HEADER_INCLUDED
|
||||
|
||||
#include <opm/autodiff/CPRPreconditioner.hpp>
|
||||
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
|
||||
#include <opm/common/utility/parameters/ParameterGroup.hpp>
|
||||
|
||||
@ -35,6 +36,7 @@ namespace Opm
|
||||
{
|
||||
/// This class carries all parameters for the NewtonIterationBlackoilInterleaved class
|
||||
struct NewtonIterationBlackoilInterleavedParameters
|
||||
: public CPRParameter
|
||||
{
|
||||
double linear_solver_reduction_;
|
||||
double ilu_relaxation_;
|
||||
@ -46,10 +48,12 @@ namespace Opm
|
||||
bool require_full_sparsity_pattern_;
|
||||
bool ignoreConvergenceFailure_;
|
||||
bool linear_solver_use_amg_;
|
||||
bool use_cpr_;
|
||||
|
||||
NewtonIterationBlackoilInterleavedParameters() { reset(); }
|
||||
// read values from parameter class
|
||||
NewtonIterationBlackoilInterleavedParameters( const ParameterGroup& param )
|
||||
: CPRParameter(param)
|
||||
{
|
||||
// set default parameters
|
||||
reset();
|
||||
@ -65,11 +69,16 @@ namespace Opm
|
||||
linear_solver_use_amg_ = param.getDefault("linear_solver_use_amg", linear_solver_use_amg_ );
|
||||
ilu_relaxation_ = param.getDefault("ilu_relaxation", ilu_relaxation_ );
|
||||
ilu_fillin_level_ = param.getDefault("ilu_fillin_level", ilu_fillin_level_ );
|
||||
|
||||
// Check whether to use cpr approach
|
||||
const std::string cprSolver = "cpr";
|
||||
use_cpr_ = ( param.getDefault("solver_approach", std::string()) == cprSolver );
|
||||
}
|
||||
|
||||
// set default values
|
||||
void reset()
|
||||
{
|
||||
use_cpr_ = false;
|
||||
newton_use_gmres_ = false;
|
||||
linear_solver_reduction_ = 1e-2;
|
||||
linear_solver_maxiter_ = 150;
|
||||
|
||||
@ -21,7 +21,7 @@
|
||||
#define OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
|
||||
|
||||
#include <opm/common/Exceptions.hpp>
|
||||
|
||||
#include <opm/common/ErrorMacros.hpp>
|
||||
#include <dune/common/version.hh>
|
||||
#include <dune/istl/preconditioner.hh>
|
||||
#include <dune/istl/paamg/smoother.hh>
|
||||
@ -85,6 +85,13 @@ namespace Opm
|
||||
template<class M, class CRS, class InvVector>
|
||||
void convertToCRS(const M& A, CRS& lower, CRS& upper, InvVector& inv )
|
||||
{
|
||||
// No need to do anything for 0 rows. Return to prevent indexing a
|
||||
// a zero sized array.
|
||||
if ( A.N() == 0 )
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
typedef typename M :: size_type size_type;
|
||||
|
||||
lower.resize( A.N() );
|
||||
@ -265,6 +272,26 @@ public:
|
||||
init( reinterpret_cast<const Matrix&>(A), n );
|
||||
}
|
||||
|
||||
/*! \brief Constructor gets all parameters to operate the prec.
|
||||
\param A The matrix to operate on.
|
||||
\param comm communication object, e.g. Dune::OwnerOverlapCopyCommunication
|
||||
\param n ILU fill in level (for testing). This does not work in parallel.
|
||||
\param w The relaxation factor.
|
||||
*/
|
||||
template<class BlockType, class Alloc>
|
||||
ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
|
||||
const ParallelInfo& comm, const int n, const field_type w )
|
||||
: lower_(),
|
||||
upper_(),
|
||||
inv_(),
|
||||
comm_(&comm), w_(w),
|
||||
relaxation_( std::abs( w - 1.0 ) > 1e-15 )
|
||||
{
|
||||
// BlockMatrix is a Subclass of FieldMatrix that just adds
|
||||
// methods. Therefore this cast should be safe.
|
||||
init( reinterpret_cast<const Matrix&>(A), n );
|
||||
}
|
||||
|
||||
/*! \brief Constructor.
|
||||
|
||||
Constructor gets all parameters to operate the prec.
|
||||
@ -328,8 +355,7 @@ public:
|
||||
const size_type lastRow = iEnd - 1;
|
||||
if( iEnd != upper_.rows() )
|
||||
{
|
||||
std::abort();
|
||||
// OPM_THROW(std::logic_error,"ILU: lower and upper rows must be the same");
|
||||
OPM_THROW(std::logic_error,"ILU: number of lower and upper rows must be the same");
|
||||
}
|
||||
|
||||
// lower triangular solve
|
||||
@ -393,6 +419,22 @@ public:
|
||||
protected:
|
||||
void init( const Matrix& A, const int iluIteration )
|
||||
{
|
||||
// (For older DUNE versions the communicator might be
|
||||
// invalid if redistribution in AMG happened on the coarset level.
|
||||
// Therefore we check for nonzero size
|
||||
if ( comm_ && comm_->communicator().size()<=0 )
|
||||
{
|
||||
if ( A.N() > 0 )
|
||||
{
|
||||
OPM_THROW(std::logic_error, "Expected a matrix with zero rows for an invalid communicator.");
|
||||
}
|
||||
else
|
||||
{
|
||||
// simply set the communicator to null
|
||||
comm_ = nullptr;
|
||||
}
|
||||
}
|
||||
|
||||
int ilu_setup_successful = 1;
|
||||
std::string message;
|
||||
const int rank = ( comm_ ) ? comm_->communicator().rank() : 0;
|
||||
|
||||
330
tests/test_blackoil_amg.cpp
Normal file
330
tests/test_blackoil_amg.cpp
Normal file
@ -0,0 +1,330 @@
|
||||
/*
|
||||
Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services
|
||||
|
||||
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>
|
||||
#if HAVE_DYNAMIC_BOOST_TEST && HAVE_MPI
|
||||
#define BOOST_TEST_DYN_LINK
|
||||
#define BOOST_TEST_MODULE BlackoilAmgTest
|
||||
#define BOOST_TEST_NO_MAIN
|
||||
#include <boost/test/unit_test.hpp>
|
||||
#include <opm/autodiff/BlackoilAmg.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/parallel/collectivecommunication.hh>
|
||||
#include <dune/istl/bcrsmatrix.hh>
|
||||
#include <dune/istl/owneroverlapcopy.hh>
|
||||
|
||||
class MPIError {
|
||||
public:
|
||||
/** @brief Constructor. */
|
||||
MPIError(std::string s, int e) : errorstring(s), errorcode(e){}
|
||||
/** @brief The error string. */
|
||||
std::string errorstring;
|
||||
/** @brief The mpi error code. */
|
||||
int errorcode;
|
||||
};
|
||||
|
||||
void MPI_err_handler(MPI_Comm *, int *err_code, ...){
|
||||
char *err_string=new char[MPI_MAX_ERROR_STRING];
|
||||
int err_length;
|
||||
MPI_Error_string(*err_code, err_string, &err_length);
|
||||
std::string s(err_string, err_length);
|
||||
std::cerr << "An MPI Error ocurred:"<<std::endl<<s<<std::endl;
|
||||
delete[] err_string;
|
||||
throw MPIError(s, *err_code);
|
||||
}
|
||||
|
||||
typedef Dune::OwnerOverlapCopyAttributeSet GridAttributes;
|
||||
typedef GridAttributes::AttributeSet GridFlag;
|
||||
typedef Dune::ParallelLocalIndex<GridFlag> LocalIndex;
|
||||
|
||||
template<class M, class G, class L, int n>
|
||||
void setupPattern(int N, M& mat, Dune::ParallelIndexSet<G,L,n>& indices, int overlapStart, int overlapEnd,
|
||||
int start, int end);
|
||||
|
||||
template<class M>
|
||||
void fillValues(int N, M& mat, int overlapStart, int overlapEnd, int start, int end);
|
||||
|
||||
|
||||
template<class M, class G, class L, class C, int n>
|
||||
M setupAnisotropic2d(int N, Dune::ParallelIndexSet<G,L,n>& indices,
|
||||
const Dune::CollectiveCommunication<C>& p, int *nout, typename M::block_type::value_type eps=1.0);
|
||||
|
||||
|
||||
template<class M, class G, class L, int s>
|
||||
void setupPattern(int N, M& mat, Dune::ParallelIndexSet<G,L,s>& indices, int overlapStart, int overlapEnd,
|
||||
int start, int end)
|
||||
{
|
||||
int n = overlapEnd - overlapStart;
|
||||
|
||||
typename M::CreateIterator iter = mat.createbegin();
|
||||
indices.beginResize();
|
||||
|
||||
for(int j=0; j < N; j++)
|
||||
for(int i=overlapStart; i < overlapEnd; i++, ++iter) {
|
||||
int global = j*N+i;
|
||||
GridFlag flag = GridAttributes::owner;
|
||||
bool isPublic = false;
|
||||
|
||||
if((i<start && i > 0) || (i>= end && i < N-1))
|
||||
flag=GridAttributes::copy;
|
||||
|
||||
if(i<start+1 || i>= end-1) {
|
||||
isPublic = true;
|
||||
indices.add(global, LocalIndex(iter.index(), flag, isPublic));
|
||||
}
|
||||
|
||||
|
||||
iter.insert(iter.index());
|
||||
|
||||
// i direction
|
||||
if(i > overlapStart )
|
||||
// We have a left neighbour
|
||||
iter.insert(iter.index()-1);
|
||||
|
||||
if(i < overlapEnd-1)
|
||||
// We have a rigt neighbour
|
||||
iter.insert(iter.index()+1);
|
||||
|
||||
// j direction
|
||||
// Overlap is a dirichlet border, discard neighbours
|
||||
if(flag != GridAttributes::copy) {
|
||||
if(j>0)
|
||||
// lower neighbour
|
||||
iter.insert(iter.index()-n);
|
||||
if(j<N-1)
|
||||
// upper neighbour
|
||||
iter.insert(iter.index()+n);
|
||||
}
|
||||
}
|
||||
indices.endResize();
|
||||
}
|
||||
|
||||
template<class M, class T>
|
||||
void fillValues(int N, M& mat, int overlapStart, int overlapEnd, int start, int end, T eps)
|
||||
{
|
||||
DUNE_UNUSED_PARAMETER(N);
|
||||
typedef typename M::block_type Block;
|
||||
Block dval(0), bone(0), bmone(0), beps(0);
|
||||
|
||||
for(typename Block::RowIterator b = dval.begin(); b != dval.end(); ++b)
|
||||
b->operator[](b.index())=2.0+2.0*eps;
|
||||
|
||||
for(typename Block::RowIterator b=bone.begin(); b != bone.end(); ++b)
|
||||
b->operator[](b.index())=1.0;
|
||||
|
||||
for(typename Block::RowIterator b=bmone.begin(); b != bmone.end(); ++b)
|
||||
b->operator[](b.index())=-1.0;
|
||||
|
||||
for(typename Block::RowIterator b=beps.begin(); b != beps.end(); ++b)
|
||||
b->operator[](b.index())=-eps;
|
||||
|
||||
int n = overlapEnd-overlapStart;
|
||||
typedef typename M::ColIterator ColIterator;
|
||||
typedef typename M::RowIterator RowIterator;
|
||||
|
||||
for (RowIterator i = mat.begin(); i != mat.end(); ++i) {
|
||||
// calculate coordinate
|
||||
int y = i.index() / n;
|
||||
int x = overlapStart + i.index() - y * n;
|
||||
|
||||
ColIterator endi = (*i).end();
|
||||
|
||||
if(x<start || x >= end) { // || x==0 || x==N-1 || y==0 || y==N-1){
|
||||
// overlap node is dirichlet border
|
||||
ColIterator j = (*i).begin();
|
||||
|
||||
for(; j.index() < i.index(); ++j)
|
||||
*j=0;
|
||||
|
||||
*j = bone;
|
||||
|
||||
for(++j; j != endi; ++j)
|
||||
*j=0;
|
||||
}else{
|
||||
for(ColIterator j = (*i).begin(); j != endi; ++j)
|
||||
if(j.index() == i.index())
|
||||
*j=dval;
|
||||
else if(j.index()+1==i.index() || j.index()-1==i.index())
|
||||
*j=beps;
|
||||
else
|
||||
*j=bmone;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<class V, class G, class L, int s>
|
||||
void setBoundary(V& lhs, V& rhs, const G& n, Dune::ParallelIndexSet<G,L,s>& indices)
|
||||
{
|
||||
typedef typename Dune::ParallelIndexSet<G,L,s>::const_iterator Iter;
|
||||
for(Iter i=indices.begin(); i != indices.end(); ++i) {
|
||||
G x = i->global()/n;
|
||||
G y = i->global()%n;
|
||||
|
||||
if(x==0 || y ==0 || x==n-1 || y==n-1) {
|
||||
//double h = 1.0 / ((double) (n-1));
|
||||
lhs[i->local()]=rhs[i->local()]=0; //((double)x)*((double)y)*h*h;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<class V, class G>
|
||||
void setBoundary(V& lhs, V& rhs, const G& N)
|
||||
{
|
||||
typedef typename V::block_type Block;
|
||||
typedef typename Block::value_type T;
|
||||
for(int j=0; j < N; ++j)
|
||||
for(int i=0; i < N; i++)
|
||||
if(i==0 || j ==0 || i==N-1 || j==N-1) {
|
||||
T h = 1.0 / ((double) (N-1));
|
||||
T x, y;
|
||||
if(i==N-1)
|
||||
x=1;
|
||||
else
|
||||
x = ((T) i)*h;
|
||||
if(j==N-1)
|
||||
y = 1;
|
||||
else
|
||||
y = ((T) j)*h;
|
||||
|
||||
lhs[j*N+i]=rhs[j*N+i]=0; //x*y;
|
||||
}
|
||||
}
|
||||
|
||||
template<class M, class G, class L, class C, int s>
|
||||
M setupAnisotropic2d(int N, Dune::ParallelIndexSet<G,L,s>& indices, const Dune::CollectiveCommunication<C>& p, int *nout, typename M::block_type::value_type eps)
|
||||
{
|
||||
int procs=p.size(), rank=p.rank();
|
||||
|
||||
typedef M BCRSMat;
|
||||
|
||||
// calculate size of local matrix in the distributed direction
|
||||
int start, end, overlapStart, overlapEnd;
|
||||
|
||||
int n = N/procs; // number of unknowns per process
|
||||
int bigger = N%procs; // number of process with n+1 unknows
|
||||
|
||||
// Compute owner region
|
||||
if(rank<bigger) {
|
||||
start = rank*(n+1);
|
||||
end = start+(n+1);
|
||||
}else{
|
||||
start = bigger*(n+1) + (rank-bigger) * n;
|
||||
end = start+n;
|
||||
}
|
||||
|
||||
// Compute overlap region
|
||||
if(start>0)
|
||||
overlapStart = start - 1;
|
||||
else
|
||||
overlapStart = start;
|
||||
|
||||
if(end<N)
|
||||
overlapEnd = end + 1;
|
||||
else
|
||||
overlapEnd = end;
|
||||
|
||||
int noKnown = overlapEnd-overlapStart;
|
||||
|
||||
*nout = noKnown;
|
||||
|
||||
BCRSMat mat(noKnown*N, noKnown*N, noKnown*N*5, BCRSMat::row_wise);
|
||||
|
||||
setupPattern(N, mat, indices, overlapStart, overlapEnd, start, end);
|
||||
fillValues(N, mat, overlapStart, overlapEnd, start, end, eps);
|
||||
|
||||
// Dune::printmatrix(std::cout,mat,"aniso 2d","row",9,1);
|
||||
|
||||
return mat;
|
||||
}
|
||||
|
||||
//BOOST_AUTO_TEST_CASE(runBlackoilAmgLaplace)
|
||||
void runBlackoilAmgLaplace()
|
||||
{
|
||||
const int BS=2, N=100;
|
||||
typedef Dune::FieldMatrix<double,BS,BS> MatrixBlock;
|
||||
typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
|
||||
typedef Dune::FieldVector<double,BS> VectorBlock;
|
||||
typedef Dune::BlockVector<VectorBlock> Vector;
|
||||
typedef int GlobalId;
|
||||
typedef Dune::OwnerOverlapCopyCommunication<GlobalId> Communication;
|
||||
typedef Dune::OverlappingSchwarzOperator<BCRSMat,Vector,Vector,Communication> Operator;
|
||||
int argc;
|
||||
char** argv;
|
||||
const auto& ccomm = Dune::MPIHelper::instance(argc, argv).getCollectiveCommunication();
|
||||
|
||||
Communication comm(ccomm);
|
||||
int n=0;
|
||||
BCRSMat mat = setupAnisotropic2d<BCRSMat>(N, comm.indexSet(), comm.communicator(), &n, 1);
|
||||
|
||||
|
||||
comm.remoteIndices().template rebuild<false>();
|
||||
|
||||
Vector b(mat.N()), x(mat.M());
|
||||
|
||||
b=0;
|
||||
x=100;
|
||||
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;
|
||||
|
||||
Opm::BlackoilAmg<Operator,ParSmoother,Criterion,Communication,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);
|
||||
|
||||
}
|
||||
|
||||
bool init_unit_test_func()
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
int main(int argc, char** argv)
|
||||
{
|
||||
const auto& helper = Dune::MPIHelper::instance(argc, argv);
|
||||
DUNE_UNUSED_PARAMETER(helper);
|
||||
boost::unit_test::unit_test_main(&init_unit_test_func,
|
||||
argc, argv);
|
||||
}
|
||||
#else
|
||||
int main () { return 0; }
|
||||
#endif // #if HAVE_DYNAMIC_BOOST_TEST
|
||||
Loading…
Reference in New Issue
Block a user