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

using different size of block size of well and reservoir

Browse Source
This commit is contained in:
Kai Bao 2017-08-03 15:14:36 +02:00
parent 894529be57
commit 3e26a8b467
3 changed files with 47 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
opm/autodiff/StandardWell.hpp
View File

@ -44,6 +44,7 @@ namespace Opm
using typename Base::MaterialLaw;
using typename Base::ModelParameters;
using typename Base::BlackoilIndices;
using typename Base::PolymerModule;
// the positions of the primary variables for StandardWell
// there are three primary variables, the second and the third ones are F_w and F_g
@ -56,22 +57,37 @@ namespace Opm
};
using typename Base::Scalar;
using typename Base::VectorBlockType;
using typename Base::MatrixBlockType;
using Base::numEq;
// TODO: with flow_ebos，for a 2P deck, // TODO: for the 2p deck, numEq will be 3, a dummy phase is already added from the reservoir side.
// it will cause problem here without processing the dummy phase.
static const int numWellEq = GET_PROP_VALUE(TypeTag, EnablePolymer)? numEq-1 : numEq; // number of wellEq is only numEq - 1 for polymer
using typename Base::Mat;
using typename Base::BVector;
using typename Base::Eval;
using typename Base::PolymerModule;
using Base::numEq;
static const int numWellEq = GET_PROP_VALUE(TypeTag, EnablePolymer)? numEq-1 : numEq; // //numEq; //number of wellEq is only numEq for polymer
typedef Dune::FieldVector<Scalar, numWellEq> VectorBlockWellType;
typedef Dune::BlockVector<VectorBlockWellType> BVectorWell;
// sparsity pattern for the matrices
//[A C^T [x = [ res
// B D ] x_well] res_well]
// the matrix type for the diagonal matrix D
typedef Dune::FieldMatrix<Scalar, numWellEq, numWellEq > DiagMatrixBlockWellType;
typedef Dune::BCRSMatrix <DiagMatrixBlockWellType> DiagMatWell;
// the matrix type for the non-diagonal matrix B and C^T
typedef Dune::FieldMatrix<Scalar, numWellEq, numEq> OffDiagMatrixBlockWellType;
typedef Dune::BCRSMatrix<OffDiagMatrixBlockWellType> OffDiagMatWell;
typedef DenseAd::Evaluation<double, /*size=*/numEq + numWellEq> EvalWell;
// TODO: should these go to WellInterface?
static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx;
static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx;
static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx;
static const int polymerConcentrationIdx = BlackoilIndices::polymerConcentrationIdx;
typedef DenseAd::Evaluation<double, /*size=*/numEq + numWellEq> EvalWell;
StandardWell(const Well* well, const int time_step, const Wells* wells);
@ -119,7 +135,7 @@ namespace Opm
const int num_cells);
// Update the well_state based on solution
void updateWellState(const BVector& dwells,
void updateWellState(const BVectorWell& dwells,
const BlackoilModelParameters& param,
WellState& well_state) const;
@ -181,10 +197,10 @@ namespace Opm
protected:
// TODO: maybe this function can go to some helper file.
void localInvert(Mat& istlA) const;
void localInvert(DiagMatWell& istlA) const;
// xw = inv(D)*(rw - C*x)
void recoverSolutionWell(const BVector& x, BVector& xw) const;
void recoverSolutionWell(const BVector& x, BVectorWell& xw) const;
// TODO: decide wether to use member function to refer to private member later
using Base::vfp_properties_;
@ -208,17 +224,17 @@ namespace Opm
// TODO: probably, they should be moved to the WellInterface, when
// we decide the template paramters.
// two off-diagonal matrices
Mat duneB_;
Mat duneC_;
OffDiagMatWell duneB_;
OffDiagMatWell duneC_;
// diagonal matrix for the well
Mat invDuneD_;
DiagMatWell invDuneD_;
// several vector used in the matrix calculation
mutable BVector Bx_;
mutable BVector invDrw_;
mutable BVectorWell Bx_;
mutable BVectorWell invDrw_;
mutable BVector scaleAddRes_;
BVector resWell_;
BVectorWell resWell_;
std::vector<EvalWell> well_variables_;
std::vector<double> F0_;

32
opm/autodiff/StandardWell_impl.hpp
View File

@ -30,9 +30,9 @@ namespace Opm
, well_variables_(numWellEq) // the number of the primary variables
, F0_(numWellEq)
{
duneB_.setBuildMode( Mat::row_wise );
duneC_.setBuildMode( Mat::row_wise );
invDuneD_.setBuildMode( Mat::row_wise );
duneB_.setBuildMode( OffDiagMatWell::row_wise );
duneC_.setBuildMode( OffDiagMatWell::row_wise );
invDuneD_.setBuildMode( DiagMatWell::row_wise );
}
@ -152,9 +152,7 @@ namespace Opm
setWellVariables(const WellState& well_state)
{
const int nw = well_state.bhp().size();
// for two-phase numComp < numWellEq
const int numComp = numComponents();
for (int eqIdx = 0; eqIdx < numComp; ++eqIdx) {
for (int eqIdx = 0; eqIdx < numWellEq; ++eqIdx) {
const unsigned int idx = nw * eqIdx + indexOfWell();
assert( eqIdx < well_variables_.size() );
assert( idx < well_state.wellSolutions().size() );
@ -665,10 +663,10 @@ namespace Opm
}
// add a trivial equation for the dummy phase for 2p cases (Only support water + oil)
if ( numComp < numWellEq ) {
/* if ( numComp < numWellEq ) {
assert(!active()[ Gas ]);
invDuneD_[0][0][Gas][Gas] = 1.0;
}
} */
// Store the perforation phase flux for later usage.
if (has_solvent && componentIdx == contiSolventEqIdx) {// if (flowPhaseToEbosCompIdx(componentIdx) == Solvent)
@ -707,9 +705,9 @@ namespace Opm
resWell_[0][componentIdx] += resWell_loc.value();
// add trivial equation for polymer
if (has_polymer) {
/* if (has_polymer) {
invDuneD_[0][0][contiPolymerEqIdx][polymerConcentrationIdx] = 1.0;
}
} */
}
// do the local inversion of D.
@ -859,7 +857,7 @@ namespace Opm
template<typename TypeTag>
void
StandardWell<TypeTag>::
updateWellState(const BVector& dwells,
updateWellState(const BVectorWell& dwells,
const BlackoilModelParameters& param,
WellState& well_state) const
{
@ -1173,7 +1171,7 @@ namespace Opm
template<typename TypeTag>
void
StandardWell<TypeTag>::
localInvert(Mat& istlA) const
localInvert(DiagMatWell& istlA) const
{
for (auto row = istlA.begin(), rowend = istlA.end(); row != rowend; ++row ) {
for (auto col = row->begin(), colend = row->end(); col != colend; ++col ) {
@ -1890,7 +1888,7 @@ namespace Opm
{
// We assemble the well equations, then we check the convergence,
// which is why we do not put the assembleWellEq here.
BVector dx_well(1);
BVectorWell dx_well(1);
invDuneD_.mv(resWell_, dx_well);
updateWellState(dx_well, param, well_state);
@ -1927,7 +1925,7 @@ namespace Opm
// invDBx = invDuneD_ * Bx_
// TODO: with this, we modified the content of the invDrw_.
// Is it necessary to do this to save some memory?
BVector& invDBx = invDrw_;
BVectorWell& invDBx = invDrw_;
invDuneD_.mv(Bx_, invDBx);
// Ax = Ax - duneC_^T * invDBx
@ -1957,9 +1955,9 @@ namespace Opm
template<typename TypeTag>
void
StandardWell<TypeTag>::
recoverSolutionWell(const BVector& x, BVector& xw) const
recoverSolutionWell(const BVector& x, BVectorWell& xw) const
{
BVector resWell = resWell_;
BVectorWell resWell = resWell_;
// resWell = resWell - B * x
duneB_.mmv(x, resWell);
// xw = D^-1 * resWell
@ -1977,7 +1975,7 @@ namespace Opm
const ModelParameters& param,
WellState& well_state) const
{
BVector xw(1);
BVectorWell xw(1);
recoverSolutionWell(x, xw);
updateWellState(xw, param, well_state);
}

2
opm/autodiff/WellInterface_impl.hpp
View File

@ -364,7 +364,7 @@ namespace Opm
{
// TODO: how about two phase polymer
if (numPhases() == 2) {
return 2;
return 2;
}
int numComp = FluidSystem::numComponents;