OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1183 from blattms/matrix-with-well-connections

Browse Source 
Add possibility to use matrix with well connections for preconditioner.
This commit is contained in:
Atgeirr Flø Rasmussen
2018-02-27 16:25:06 +01:00
committed by GitHub
parent f0b8e586a6 44d31f9dcc
commit 10887ddaca
14 changed files with 307 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists_files.cmake
View File

@@ -147,6 +147,7 @@ list (APPEND TEST_SOURCE_FILES
tests/test_welldensitysegmented.cpp
tests/test_vfpproperties.cpp
tests/test_singlecellsolves.cpp
tests/test_multmatrixtransposed.cpp
tests/test_multiphaseupwind.cpp
tests/test_wellmodel.cpp
# tests/test_thresholdpressure.cpp
@@ -315,6 +316,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/SimulatorIncompTwophaseAd.hpp
opm/autodiff/SimulatorSequentialBlackoil.hpp
opm/autodiff/TransportSolverTwophaseAd.hpp
opm/autodiff/WellConnectionAuxiliaryModule.hpp
opm/autodiff/WellDensitySegmented.hpp
opm/autodiff/WellStateFullyImplicitBlackoil.hpp
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp

19
opm/autodiff/BlackoilModelEbos.hpp
View File

@@ -29,6 +29,7 @@
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/BlackoilWellModel.hpp>
#include <opm/autodiff/WellConnectionAuxiliaryModule.hpp>
#include <opm/autodiff/BlackoilDetails.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
@@ -487,14 +488,17 @@ namespace Opm {
if( isParallel() )
{
typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, true > Operator;
Operator opA(ebosJac, wellModel(), istlSolver().parallelInformation() );
Operator opA(ebosJac, wellModel(),
param_.matrix_add_well_contributions_,
istlSolver().parallelInformation() );
assert( opA.comm() );
istlSolver().solve( opA, x, ebosResid, *(opA.comm()) );
}
else
{
typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, false > Operator;
Operator opA(ebosJac, wellModel());
Operator opA(ebosJac, wellModel(),
param_.matrix_add_well_contributions_ );
istlSolver().solve( opA, x, ebosResid );
}
}
@@ -541,8 +545,11 @@ namespace Opm {
#endif
//! constructor: just store a reference to a matrix
WellModelMatrixAdapter (const M& A, const WellModel& wellMod, const boost::any& parallelInformation = boost::any() )
: A_( A ), wellMod_( wellMod ), comm_()
WellModelMatrixAdapter (const M& A, const WellModel& wellMod,
bool matrix_add_well_contributions,
const boost::any& parallelInformation = boost::any() )
: A_( A ), wellMod_( wellMod ), comm_(),
matrix_add_well_contributions_(matrix_add_well_contributions)
{
#if HAVE_MPI
if( parallelInformation.type() == typeid(ParallelISTLInformation) )
@@ -557,6 +564,7 @@ namespace Opm {
virtual void apply( const X& x, Y& y ) const
{
A_.mv( x, y );
// add well model modification to y
wellMod_.apply(x, y );
@@ -570,6 +578,7 @@ namespace Opm {
virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
{
A_.usmv(alpha,x,y);
// add scaled well model modification to y
wellMod_.applyScaleAdd( alpha, x, y );
@@ -590,6 +599,7 @@ namespace Opm {
const matrix_type& A_ ;
const WellModel& wellMod_;
std::unique_ptr< communication_type > comm_;
bool matrix_add_well_contributions_;
};
/// Apply an update to the primary variables, chopped if appropriate.
@@ -1067,6 +1077,7 @@ namespace Opm {
private:
double dpMaxRel() const { return param_.dp_max_rel_; }
double dsMax() const { return param_.ds_max_; }
double drMaxRel() const { return param_.dr_max_rel_; }

2
opm/autodiff/BlackoilModelParameters.cpp
View File

@@ -65,6 +65,7 @@ namespace Opm
update_equations_scaling_ = param.getDefault("update_equations_scaling", update_equations_scaling_);
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_);
}
@@ -93,6 +94,7 @@ namespace Opm
update_equations_scaling_ = false;
use_update_stabilization_ = true;
use_multisegment_well_ = false;
matrix_add_well_contributions_ = false;
}

3
opm/autodiff/BlackoilModelParameters.hpp
View File

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

10
opm/autodiff/BlackoilWellModel.hpp
View File

@@ -84,6 +84,15 @@ namespace Opm {
typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
typedef Dune::BlockVector<VectorBlockType> BVector;
#if DUNE_VERSION_NEWER_REV(DUNE_ISTL, 2 , 5, 1)
// 3x3 matrix block inversion was unstable from at least 2.3 until and
// including 2.5.0
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
#else
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
#endif
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
// For the conversion between the surface volume rate and resrevoir voidage rate
@@ -143,6 +152,7 @@ namespace Opm {
const SimulatorReport& lastReport() const;
protected:
Simulator& ebosSimulator_;

9
opm/autodiff/BlackoilWellModel_impl.hpp
View File

@@ -294,10 +294,6 @@ namespace Opm {
}
}
// applying the well residual to reservoir residuals
// r = r - duneC_^T * invDuneD_ * resWell_
template<typename TypeTag>
@@ -330,7 +326,10 @@ namespace Opm {
}
for (auto& well : well_container_) {
well->apply(x, Ax);
if ( ! well->jacobianContainsWellContributions() )
{
well->apply(x, Ax);
}
}
}

21
opm/autodiff/ISTLSolver.hpp
View File

@@ -310,6 +310,27 @@ public:
namespace Opm
{
namespace Detail
{
//! calculates ret = A^T * B
template< class K, int m, int n, int p >
static inline void multMatrixTransposed ( const Dune::FieldMatrix< K, n, m > &A,
const Dune::FieldMatrix< K, n, p > &B,
Dune::FieldMatrix< K, m, p > &ret )
{
typedef typename Dune::FieldMatrix< K, m, p > :: size_type size_type;
for( size_type i = 0; i < m; ++i )
{
for( size_type j = 0; j < p; ++j )
{
ret[ i ][ j ] = K( 0 );
for( size_type k = 0; k < n; ++k )
ret[ i ][ j ] += A[ k ][ i ] * B[ k ][ j ];
}
}
}
}
/// This class solves the fully implicit black-oil system by
/// solving the reduced system (after eliminating well variables)
/// as a block-structured matrix (one block for all cell variables) for a fixed

7
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@@ -179,6 +179,13 @@ public:
WellState wellStateDummy; //not used. Only passed to make the old interfaces happy
if ( model_param_.matrix_add_well_contributions_ )
{
ebosSimulator_.model().clearAuxiliaryModules();
auto auxMod = std::make_shared<WellConnectionAuxiliaryModule<TypeTag> >(schedule(), grid());
ebosSimulator_.model().addAuxiliaryModule(auxMod);
}
// Main simulation loop.
while (!timer.done()) {
// Report timestep.

16
opm/autodiff/StandardWell.hpp
View File

@@ -27,6 +27,7 @@
#include <opm/autodiff/WellInterface.hpp>
#include <opm/autodiff/ISTLSolver.hpp>
#include <opm/autodiff/RateConverter.hpp>
#include <opm/autodiff/ISTLSolver.hpp>
namespace Opm
{
@@ -88,15 +89,8 @@ namespace Opm
typedef Dune::FieldVector<Scalar, numWellEq> VectorBlockWellType;
typedef Dune::BlockVector<VectorBlockWellType> BVectorWell;
#if DUNE_VERSION_NEWER_REV(DUNE_ISTL, 2 , 5, 1)
// 3x3 matrix block inversion was unstable from at least 2.3 until and
// including 2.5.0
// the matrix type for the diagonal matrix D
typedef Dune::FieldMatrix<Scalar, numWellEq, numWellEq > DiagMatrixBlockWellType;
#else
// the matrix type for the diagonal matrix D
typedef Dune::MatrixBlock<Scalar, numWellEq, numWellEq > DiagMatrixBlockWellType;
#endif
typedef Dune::BCRSMatrix <DiagMatrixBlockWellType> DiagMatWell;
@@ -157,6 +151,14 @@ namespace Opm
virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
const WellState& well_state); // should be const?
virtual void addWellContributions(Mat& mat) const;
/// \brief Wether the Jacobian will also have well contributions in it.
virtual bool jacobianContainsWellContributions() const
{
return param_.matrix_add_well_contributions_;
}
protected:
// protected functions from the Base class

43
opm/autodiff/StandardWell_impl.hpp
View File

@@ -22,6 +22,7 @@
namespace Opm
{
template<typename TypeTag>
StandardWell<TypeTag>::
StandardWell(const Well* well, const int time_step, const Wells* wells,
@@ -666,7 +667,15 @@ namespace Opm
}
// do the local inversion of D.
invDuneD_[0][0].invert();
// we do this manually with invertMatrix to always get our
// 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 );
}
}
@@ -2073,4 +2082,36 @@ namespace Opm
mob[waterCompIdx] /= shear_factor;
}
}
template<typename TypeTag>
void
StandardWell<TypeTag>::addWellContributions(Mat& mat) const
{
// We need to change matrx A as follows
// A -= C^T D^-1 B
// D is diagonal
// B and C have 1 row, nc colums and nonzero
// at (0,j) only if this well has a perforation at cell j.
for ( auto colC = duneC_[0].begin(), endC = duneC_[0].end(); colC != endC; ++colC )
{
const auto row_index = colC.index();
auto& row = mat[row_index];
auto col = row.begin();
for ( auto colB = duneB_[0].begin(), endB = duneB_[0].end(); colB != endB; ++colB )
{
const auto col_index = colB.index();
// Move col to index col_index
while ( col != row.end() && col.index() < col_index ) ++col;
assert(col != row.end() && col.index() == col_index);
Dune::FieldMatrix<Scalar, numWellEq, numEq> tmp;
typename Mat::block_type tmp1;
Dune::FMatrixHelp::multMatrix(invDuneD_[0][0], (*colB), tmp);
Detail::multMatrixTransposed((*colC), tmp, tmp1);
(*col) -= tmp1;
}
}
}
}

130
opm/autodiff/WellConnectionAuxiliaryModule.hpp Normal file
View File

@@ -0,0 +1,130 @@
/*
Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2017 Statoil ASA.
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_WELLCONNECTIONAUXILIARYMODULE_HEADER_INCLUDED
#define OPM_WELLCONNECTIONAUXILIARYMODULE_HEADER_INCLUDED
#include <ewoms/aux/baseauxiliarymodule.hh>
#include <opm/grid/CpGrid.hpp>
#include <opm/core/wells.h>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
namespace Opm
{
template<class TypeTag>
class WellConnectionAuxiliaryModule
: public Ewoms::BaseAuxiliaryModule<TypeTag>
{
typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
public:
using NeighborSet = typename
Ewoms::BaseAuxiliaryModule<TypeTag>::NeighborSet;
WellConnectionAuxiliaryModule(const Schedule& schedule,
const Dune::CpGrid& grid)
{
// Create cartesian to compressed mapping
const auto& globalCell = grid.globalCell();
const auto& cartesianSize = grid.logicalCartesianSize();
auto size = cartesianSize[0]*cartesianSize[1]*cartesianSize[2];
std::vector<int> cartesianToCompressed(size, -1);
auto begin = globalCell.begin();
for ( auto cell = begin, end= globalCell.end(); cell != end; ++cell )
{
cartesianToCompressed[ *cell ] = cell - begin;
}
int last_time_step = schedule.getTimeMap().size() - 1;
const auto& schedule_wells = schedule.getWells();
wells_.reserve(schedule_wells.size());
// initialize the additional cell connections introduced by wells.
for ( const auto well : schedule_wells )
{
std::vector<int> compressed_well_perforations;
// All possible completions of the well
const auto& completionSet = well->getCompletions(last_time_step);
compressed_well_perforations.reserve(completionSet.size());
for ( size_t c=0; c < completionSet.size(); c++ )
{
const auto& completion = completionSet.get(c);
int i = completion.getI();
int j = completion.getJ();
int k = completion.getK();
int cart_grid_idx = i + cartesianSize[0]*(j + cartesianSize[1]*k);
int compressed_idx = cartesianToCompressed[cart_grid_idx];
if ( compressed_idx >= 0 ) // Ignore completions in inactive/remote cells.
{
compressed_well_perforations.push_back(compressed_idx);
}
}
if ( ! compressed_well_perforations.empty() )
{
std::sort(compressed_well_perforations.begin(),
compressed_well_perforations.end());
wells_.push_back(compressed_well_perforations);
}
}
}
unsigned numDofs() const
{
// No extra dofs are inserted for wells.
return 0;
}
void addNeighbors(std::vector<NeighborSet>& neighbors) const
{
for(const auto well_perforations : wells_)
{
for(const auto& perforation : well_perforations)
neighbors[perforation].insert(well_perforations.begin(),
well_perforations.end());
}
}
void applyInitial()
{}
void linearize(JacobianMatrix& , GlobalEqVector&)
{
// Linearization is done in StandardDenseWells
}
private:
std::vector<std::vector<int> > wells_;
};
} // end namespace OPM
#endif

9
opm/autodiff/WellInterface.hpp
View File

@@ -24,7 +24,8 @@
#define OPM_WELLINTERFACE_HEADER_INCLUDED
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/Exceptions.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/core/wells.h>
@@ -202,6 +203,12 @@ namespace Opm
virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
const WellState& well_state) = 0; // should be const?
/// \brief Wether the Jacobian will also have well contributions in it.
virtual bool jacobianContainsWellContributions() const
{
return false;
}
// updating the voidage rates in well_state when requested
void calculateReservoirRates(WellState& well_state) const;

1
opm/autodiff/WellInterface_impl.hpp
View File

@@ -847,7 +847,6 @@ namespace Opm
template<typename TypeTag>
void
WellInterface<TypeTag>::calculateReservoirRates(WellState& well_state) const

54
tests/test_multmatrixtransposed.cpp Normal file
View File

@@ -0,0 +1,54 @@
/*
Copyright 2018 Statoil ASA
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>
#define BOOST_TEST_MODULE MultMatrixTransposed
#include <boost/test/unit_test.hpp>
#include <opm/autodiff/ISTLSolver.hpp>
using namespace Dune;
using namespace Opm::Detail;
BOOST_AUTO_TEST_CASE(testmultmatrixtrans)
{
using Scalar = FieldMatrix<double,1,1>;
Scalar a1 = {{ 4 }}, b1 = {{ 2 }}, res1 = {{}}, resExpect1 = {{ 8 }};
multMatrixTransposed(a1, b1, res1);
BOOST_CHECK_EQUAL(res1, resExpect1);
using Mat2 = FieldMatrix<double,2,2>;
Mat2 a2 = {{ 1, 2 }, { 3, 4} }, b2 = {{ 3, 4 }, { 5, 6} };
Mat2 res2 = {{}, {}}, resExpect2 = {{ 18, 22 }, {26, 32} };
multMatrixTransposed(a2, b2, res2);
BOOST_CHECK_EQUAL(res2, resExpect2);
using Mat3 = FieldMatrix<double,3,3>;
Mat3 a3 = {{ 1, 2, 3 }, { 3, 4, 5}, {6, 7, 8} };
Mat3 b3 = {{ 3, 4, 5 }, { 5, 6, 7}, {7, 8, 9} };
Mat3 res3 = {{}, {}, {}};
Mat3 resExpect3 = {{ 60, 70, 80 }, {75, 88, 101}, { 90, 106, 122 } };
multMatrixTransposed(a3, b3, res3);
BOOST_CHECK_EQUAL(res3, resExpect3);
}