Merge pull request #1422 from blattms/resort-to-split-operator

If requested use matrix with well contributions only for the preconditioner.
This commit is contained in:
Atgeirr Flø Rasmussen 2018-03-07 14:33:09 +01:00 committed by GitHub
commit 717a9180e5
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
8 changed files with 45 additions and 21 deletions

View File

@ -374,6 +374,16 @@ namespace Opm {
OPM_THROW(Opm::NumericalIssue,"Error encounted when solving well equations");
}
auto& ebosJac = ebosSimulator_.model().linearizer().matrix();
if (param_.matrix_add_well_contributions_) {
wellModel().addWellContributions(ebosJac);
}
if ( param_.preconditioner_add_well_contributions_ &&
! param_.matrix_add_well_contributions_ ) {
matrix_for_preconditioner_ .reset(new Mat(ebosJac));
wellModel().addWellContributions(*matrix_for_preconditioner_);
}
return wellModel().lastReport();
}
@ -484,12 +494,12 @@ namespace Opm {
// set initial guess
x = 0.0;
const Mat& actual_mat_for_prec = matrix_for_preconditioner_ ? *matrix_for_preconditioner_.get() : ebosJac;
// Solve system.
if( isParallel() )
{
typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, true > Operator;
Operator opA(ebosJac, wellModel(),
param_.matrix_add_well_contributions_,
Operator opA(ebosJac, actual_mat_for_prec, wellModel(),
istlSolver().parallelInformation() );
assert( opA.comm() );
istlSolver().solve( opA, x, ebosResid, *(opA.comm()) );
@ -497,8 +507,7 @@ namespace Opm {
else
{
typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, false > Operator;
Operator opA(ebosJac, wellModel(),
param_.matrix_add_well_contributions_ );
Operator opA(ebosJac, actual_mat_for_prec, wellModel());
istlSolver().solve( opA, x, ebosResid );
}
}
@ -545,11 +554,11 @@ namespace Opm {
#endif
//! constructor: just store a reference to a matrix
WellModelMatrixAdapter (const M& A, const WellModel& wellMod,
bool matrix_add_well_contributions,
WellModelMatrixAdapter (const M& A,
const M& A_for_precond,
const WellModel& wellMod,
const boost::any& parallelInformation = boost::any() )
: A_( A ), wellMod_( wellMod ), comm_(),
matrix_add_well_contributions_(matrix_add_well_contributions)
: A_( A ), A_for_precond_(A_for_precond), wellMod_( wellMod ), comm_()
{
#if HAVE_MPI
if( parallelInformation.type() == typeid(ParallelISTLInformation) )
@ -588,7 +597,7 @@ namespace Opm {
#endif
}
virtual const matrix_type& getmat() const { return A_; }
virtual const matrix_type& getmat() const { return A_for_precond_; }
communication_type* comm()
{
@ -597,9 +606,9 @@ namespace Opm {
protected:
const matrix_type& A_ ;
const matrix_type& A_for_precond_ ;
const WellModel& wellMod_;
std::unique_ptr< communication_type > comm_;
bool matrix_add_well_contributions_;
};
/// Apply an update to the primary variables, chopped if appropriate.
@ -1057,6 +1066,8 @@ namespace Opm {
double current_relaxation_;
BVector dx_old_;
std::unique_ptr<Mat> matrix_for_preconditioner_;
public:
/// return the StandardWells object
BlackoilWellModel<TypeTag>&

View File

@ -66,6 +66,7 @@ namespace Opm
use_update_stabilization_ = param.getDefault("use_update_stabilization", use_update_stabilization_);
deck_file_name_ = param.template get<std::string>("deck_filename");
matrix_add_well_contributions_ = param.getDefault("matrix_add_well_contributions", matrix_add_well_contributions_);
preconditioner_add_well_contributions_ = param.getDefault("preconditioner_add_well_contributions", preconditioner_add_well_contributions_);
}
@ -95,6 +96,7 @@ namespace Opm
use_update_stabilization_ = true;
use_multisegment_well_ = false;
matrix_add_well_contributions_ = false;
preconditioner_add_well_contributions_ = false;
}

View File

@ -91,9 +91,12 @@ namespace Opm
/// The file name of the deck
std::string deck_file_name_;
// Whether to add influences of wells between cells to the matrix
// Whether to add influences of wells between cells to the matrix and preconditioner matrix
bool matrix_add_well_contributions_;
// Whether to add influences of wells between cells to the preconditioner matrix only
bool preconditioner_add_well_contributions_;
/// Construct from user parameters or defaults.
explicit BlackoilModelParameters( const ParameterGroup& param );

View File

@ -153,6 +153,13 @@ namespace Opm {
const SimulatorReport& lastReport() const;
void addWellContributions(Mat& mat)
{
for ( const auto& well: well_container_ ) {
well->addWellContributions(mat);
}
}
protected:
Simulator& ebosSimulator_;

View File

@ -326,10 +326,7 @@ namespace Opm {
}
for (auto& well : well_container_) {
if ( ! well->jacobianContainsWellContributions() )
{
well->apply(x, Ax);
}
well->apply(x, Ax);
}
}

View File

@ -179,7 +179,8 @@ public:
WellState wellStateDummy; //not used. Only passed to make the old interfaces happy
if ( model_param_.matrix_add_well_contributions_ )
if ( model_param_.matrix_add_well_contributions_ ||
model_param_.preconditioner_add_well_contributions_ )
{
ebosSimulator_.model().clearAuxiliaryModules();
auto auxMod = std::make_shared<WellConnectionAuxiliaryModule<TypeTag> >(schedule(), grid());

View File

@ -671,11 +671,6 @@ namespace Opm
// specializations in for 3x3 and 4x4 matrices.
auto original = invDuneD_[0][0];
Dune::FMatrixHelp::invertMatrix(original, invDuneD_[0][0]);
if ( param_.matrix_add_well_contributions_ )
{
addWellContributions( ebosJac );
}
}
@ -1604,6 +1599,11 @@ namespace Opm
StandardWell<TypeTag>::
apply(const BVector& x, BVector& Ax) const
{
if ( param_.matrix_add_well_contributions_ )
{
// Contributions are already in the matrix itself
return;
}
assert( Bx_.size() == duneB_.N() );
assert( invDrw_.size() == invDuneD_.N() );

View File

@ -212,6 +212,9 @@ namespace Opm
// updating the voidage rates in well_state when requested
void calculateReservoirRates(WellState& well_state) const;
// Add well contributions to matrix
virtual void addWellContributions(Mat&) const
{}
protected:
// to indicate a invalid connection