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 #4302 from akva2/msw_equations

Browse Source 
Added: MultisegmentWellEquations
This commit is contained in:
Bård Skaflestad 2022-12-05 11:05:31 +01:00 committed by GitHub
commit 4f38217fc8
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
12 changed files with 673 additions and 395 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
CMakeLists_files.cmake
View File

@ -90,6 +90,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/simulators/wells/GlobalWellInfo.cpp
opm/simulators/wells/GroupState.cpp
opm/simulators/wells/MSWellHelpers.cpp
opm/simulators/wells/MultisegmentWellEquations.cpp
opm/simulators/wells/MultisegmentWellEval.cpp
opm/simulators/wells/MultisegmentWellGeneric.cpp
opm/simulators/wells/ParallelWellInfo.cpp
@ -371,6 +372,9 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/wells/MSWellHelpers.hpp
opm/simulators/wells/MultisegmentWell.hpp
opm/simulators/wells/MultisegmentWell_impl.hpp
opm/simulators/wells/MultisegmentWellEquations.hpp
opm/simulators/wells/MultisegmentWellEval.hpp
opm/simulators/wells/MultisegmentWellGeneric.hpp
opm/simulators/wells/ParallelWellInfo.hpp
opm/simulators/wells/PerfData.hpp
opm/simulators/wells/PerforationData.hpp

2
opm/simulators/wells/BlackoilWellModel_impl.hpp
View File

@ -1219,7 +1219,7 @@ namespace Opm {
} else {
auto derived_ms = std::dynamic_pointer_cast<MultisegmentWell<TypeTag> >(well);
if (derived_ms) {
derived_ms->addWellContribution(wellContribs);
derived_ms->linSys().extract(wellContribs);
} else {
OpmLog::warning("Warning unknown type of well");
}

36
opm/simulators/wells/MSWellHelpers.cpp
View File

@ -100,16 +100,10 @@ namespace mswellhelpers
/// Applies umfpack and checks for singularity
template <typename MatrixType, typename VectorType>
VectorType
applyUMFPack(const MatrixType& D,
std::shared_ptr<Dune::UMFPack<MatrixType>>& linsolver,
applyUMFPack(Dune::UMFPack<MatrixType>& linsolver,
VectorType x)
{
#if HAVE_UMFPACK
if (!linsolver)
{
linsolver = std::make_shared<Dune::UMFPack<MatrixType>>(D, 0);
}
// The copy of x seems mandatory for calling UMFPack!
VectorType y(x.size());
y = 0.;
@ -118,7 +112,7 @@ applyUMFPack(const MatrixType& D,
Dune::InverseOperatorResult res;
// Solve
linsolver->apply(y, x, res);
linsolver.apply(y, x, res);
// Checking if there is any inf or nan in y
// it will be the solution before we find a way to catch the singularity of the matrix
@ -141,24 +135,24 @@ applyUMFPack(const MatrixType& D,
template <typename VectorType, typename MatrixType>
Dune::Matrix<typename MatrixType::block_type>
invertWithUMFPack(const MatrixType& D, std::shared_ptr<Dune::UMFPack<MatrixType> >& linsolver)
invertWithUMFPack(const int size,
const int bsize,
Dune::UMFPack<MatrixType>& linsolver)
{
#if HAVE_UMFPACK
const int sz = D.M();
const int bsz = D[0][0].M();
VectorType e(sz);
VectorType e(size);
e = 0.0;
// Make a full block matrix.
Dune::Matrix<typename MatrixType::block_type> inv(sz, sz);
Dune::Matrix<typename MatrixType::block_type> inv(size, size);
// Create inverse by passing basis vectors to the solver.
for (int ii = 0; ii < sz; ++ii) {
for (int jj = 0; jj < bsz; ++jj) {
for (int ii = 0; ii < size; ++ii) {
for (int jj = 0; jj < bsize; ++jj) {
e[ii][jj] = 1.0;
auto col = applyUMFPack(D, linsolver, e);
for (int cc = 0; cc < sz; ++cc) {
for (int dd = 0; dd < bsz; ++dd) {
auto col = applyUMFPack(linsolver, e);
for (int cc = 0; cc < size; ++cc) {
for (int dd = 0; dd < bsize; ++dd) {
inv[cc][ii][dd][jj] = col[cc][dd];
}
}
@ -328,12 +322,10 @@ template<int Dim>
using Mat = Dune::BCRSMatrix<Dune::FieldMatrix<double,Dim,Dim>>;
#define INSTANCE_UMF(Dim) \
template Vec<Dim> applyUMFPack<Mat<Dim>,Vec<Dim>>(const Mat<Dim>&, \
std::shared_ptr<Dune::UMFPack<Mat<Dim>>>&, \
template Vec<Dim> applyUMFPack<Mat<Dim>,Vec<Dim>>(Dune::UMFPack<Mat<Dim>>&, \
Vec<Dim>); \
template Dune::Matrix<typename Mat<Dim>::block_type> \
invertWithUMFPack<Vec<Dim>,Mat<Dim>>(const Mat<Dim>& D, \
std::shared_ptr<Dune::UMFPack<Mat<Dim>>>&);
invertWithUMFPack<Vec<Dim>,Mat<Dim>>(const int, const int, Dune::UMFPack<Mat<Dim>>&);
INSTANCE_UMF(2)
INSTANCE_UMF(3)

8
opm/simulators/wells/MSWellHelpers.hpp
View File

@ -42,8 +42,7 @@ namespace mswellhelpers
/// Applies umfpack and checks for singularity
template <typename MatrixType, typename VectorType>
VectorType
applyUMFPack(const MatrixType& D,
std::shared_ptr<Dune::UMFPack<MatrixType>>& linsolver,
applyUMFPack(Dune::UMFPack<MatrixType>& linsolver,
VectorType x);
@ -51,8 +50,9 @@ namespace mswellhelpers
/// Applies umfpack and checks for singularity
template <typename VectorType, typename MatrixType>
Dune::Matrix<typename MatrixType::block_type>
invertWithUMFPack(const MatrixType& D,
std::shared_ptr<Dune::UMFPack<MatrixType> >& linsolver);
invertWithUMFPack(const int size,
const int bsize,
Dune::UMFPack<MatrixType>& linsolver);

15
opm/simulators/wells/MultisegmentWell.hpp
View File

@ -65,10 +65,9 @@ namespace Opm
using typename Base::BVector;
using typename Base::Eval;
using typename MSWEval::Equations;
using typename MSWEval::EvalWell;
using typename MSWEval::BVectorWell;
using typename MSWEval::DiagMatWell;
using typename MSWEval::OffDiagMatrixBlockWellType;
using MSWEval::GFrac;
using MSWEval::WFrac;
using MSWEval::WQTotal;
@ -138,13 +137,13 @@ namespace Opm
WellState& well_state,
DeferredLogger& deferred_logger) const override;
virtual void addWellContributions(SparseMatrixAdapter& jacobian) const override;
void addWellContributions(SparseMatrixAdapter& jacobian) const override;
virtual void addWellPressureEquations(PressureMatrix& mat,
const BVector& x,
const int pressureVarIndex,
const bool use_well_weights,
const WellState& well_state) const override;
void addWellPressureEquations(PressureMatrix& mat,
const BVector& x,
const int pressureVarIndex,
const bool use_well_weights,
const WellState& well_state) const override;
virtual std::vector<double> computeCurrentWellRates(const Simulator& ebosSimulator,
DeferredLogger& deferred_logger) const override;

397
opm/simulators/wells/MultisegmentWellEquations.cpp Normal file
View File

@ -0,0 +1,397 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
Copyright 2016 - 2017 IRIS 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/>.
*/
#include <config.h>
#include <opm/simulators/wells/MultisegmentWellEquations.hpp>
#include <dune/istl/umfpack.hh>
#include <opm/common/ErrorMacros.hpp>
#include <opm/simulators/linalg/bda/WellContributions.hpp>
#include <opm/simulators/linalg/istlsparsematrixadapter.hh>
#include <opm/simulators/linalg/matrixblock.hh>
#include <opm/simulators/linalg/SmallDenseMatrixUtils.hpp>
#include <opm/simulators/wells/MSWellHelpers.hpp>
#include <opm/simulators/wells/MultisegmentWellGeneric.hpp>
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
#include <stdexcept>
namespace Opm {
template<class Scalar, int numWellEq, int numEq>
MultisegmentWellEquations<Scalar,numWellEq,numEq>::
MultisegmentWellEquations(const MultisegmentWellGeneric<Scalar>& well)
: well_(well)
{
}
template<class Scalar, int numWellEq, int numEq>
void MultisegmentWellEquations<Scalar,numWellEq,numEq>::
init(const int num_cells,
const int numPerfs,
const std::vector<int>& cells)
{
duneB_.setBuildMode(OffDiagMatWell::row_wise);
duneC_.setBuildMode(OffDiagMatWell::row_wise);
duneD_.setBuildMode(DiagMatWell::row_wise);
// set the size and patterns for all the matrices and vectors
// [A C^T [x = [ res
// B D] x_well] res_well]
// calculating the NNZ for duneD_
// NNZ = number_of_segments + 2 * (number_of_inlets / number_of_outlets)
{
int nnz_d = well_.numberOfSegments();
for (const std::vector<int>& inlets : well_.segmentInlets()) {
nnz_d += 2 * inlets.size();
}
duneD_.setSize(well_.numberOfSegments(), well_.numberOfSegments(), nnz_d);
}
duneB_.setSize(well_.numberOfSegments(), num_cells, numPerfs);
duneC_.setSize(well_.numberOfSegments(), num_cells, numPerfs);
// we need to add the off diagonal ones
for (auto row = duneD_.createbegin(),
end = duneD_.createend(); row != end; ++row) {
// the number of the row corrspnds to the segment now
const int seg = row.index();
// adding the item related to outlet relation
const Segment& segment = well_.segmentSet()[seg];
const int outlet_segment_number = segment.outletSegment();
if (outlet_segment_number > 0) { // if there is a outlet_segment
const int outlet_segment_index = well_.segmentNumberToIndex(outlet_segment_number);
row.insert(outlet_segment_index);
}
// Add nonzeros for diagonal
row.insert(seg);
// insert the item related to its inlets
for (const int& inlet : well_.segmentInlets()[seg]) {
row.insert(inlet);
}
}
// make the C matrix
for (auto row = duneC_.createbegin(),
end = duneC_.createend(); row != end; ++row) {
// the number of the row corresponds to the segment number now.
for (const int& perf : well_.segmentPerforations()[row.index()]) {
const int cell_idx = cells[perf];
row.insert(cell_idx);
}
}
// make the B^T matrix
for (auto row = duneB_.createbegin(),
end = duneB_.createend(); row != end; ++row) {
// the number of the row corresponds to the segment number now.
for (const int& perf : well_.segmentPerforations()[row.index()]) {
const int cell_idx = cells[perf];
row.insert(cell_idx);
}
}
resWell_.resize(well_.numberOfSegments());
}
template<class Scalar, int numWellEq, int numEq>
void MultisegmentWellEquations<Scalar,numWellEq,numEq>::clear()
{
duneB_ = 0.0;
duneC_ = 0.0;
duneD_ = 0.0;
resWell_ = 0.0;
duneDSolver_.reset();
}
template<class Scalar, int numWellEq, int numEq>
void MultisegmentWellEquations<Scalar,numWellEq,numEq>::
apply(const BVector& x, BVector& Ax) const
{
BVectorWell Bx(duneB_.N());
duneB_.mv(x, Bx);
// invDBx = duneD^-1 * Bx_
const BVectorWell invDBx = mswellhelpers::applyUMFPack(*duneDSolver_, Bx);
// Ax = Ax - duneC_^T * invDBx
duneC_.mmtv(invDBx,Ax);
}
template<class Scalar, int numWellEq, int numEq>
void MultisegmentWellEquations<Scalar,numWellEq,numEq>::
apply(BVector& r) const
{
// invDrw_ = duneD^-1 * resWell_
const BVectorWell invDrw = mswellhelpers::applyUMFPack(*duneDSolver_, resWell_);
// r = r - duneC_^T * invDrw
duneC_.mmtv(invDrw, r);
}
template<class Scalar, int numWellEq, int numEq>
void MultisegmentWellEquations<Scalar,numWellEq,numEq>::createSolver()
{
#if HAVE_UMFPACK
if (duneDSolver_) {
return;
}
duneDSolver_ = std::make_shared<Dune::UMFPack<DiagMatWell>>(duneD_, 0);
#else
OPM_THROW(std::runtime_error, "MultisegmentWell support requires UMFPACK. "
"Reconfigure opm-simulators with SuiteSparse/UMFPACK support and recompile.");
#endif
}
template<class Scalar, int numWellEq, int numEq>
typename MultisegmentWellEquations<Scalar,numWellEq,numEq>::BVectorWell
MultisegmentWellEquations<Scalar,numWellEq,numEq>::solve() const
{
return mswellhelpers::applyUMFPack(*duneDSolver_, resWell_);
}
template<class Scalar, int numWellEq, int numEq>
void MultisegmentWellEquations<Scalar,numWellEq,numEq>::
recoverSolutionWell(const BVector& x, BVectorWell& xw) const
{
BVectorWell resWell = resWell_;
// resWell = resWell - B * x
duneB_.mmv(x, resWell);
// xw = D^-1 * resWell
xw = mswellhelpers::applyUMFPack(*duneDSolver_, resWell);
}
template<class Scalar, int numWellEq, int numEq>
void MultisegmentWellEquations<Scalar,numWellEq,numEq>::
extract(WellContributions& wellContribs) const
{
unsigned int Mb = duneB_.N(); // number of blockrows in duneB_, duneC_ and duneD_
unsigned int BnumBlocks = duneB_.nonzeroes();
unsigned int DnumBlocks = duneD_.nonzeroes();
// duneC
std::vector<unsigned int> Ccols;
std::vector<double> Cvals;
Ccols.reserve(BnumBlocks);
Cvals.reserve(BnumBlocks * numEq * numWellEq);
for (auto rowC = duneC_.begin(); rowC != duneC_.end(); ++rowC) {
for (auto colC = rowC->begin(), endC = rowC->end(); colC != endC; ++colC) {
Ccols.emplace_back(colC.index());
for (int i = 0; i < numWellEq; ++i) {
for (int j = 0; j < numEq; ++j) {
Cvals.emplace_back((*colC)[i][j]);
}
}
}
}
// duneD
Dune::UMFPack<DiagMatWell> umfpackMatrix(duneD_, 0);
double* Dvals = umfpackMatrix.getInternalMatrix().getValues();
auto* Dcols = umfpackMatrix.getInternalMatrix().getColStart();
auto* Drows = umfpackMatrix.getInternalMatrix().getRowIndex();
// duneB
std::vector<unsigned int> Bcols;
std::vector<unsigned int> Brows;
std::vector<double> Bvals;
Bcols.reserve(BnumBlocks);
Brows.reserve(Mb+1);
Bvals.reserve(BnumBlocks * numEq * numWellEq);
Brows.emplace_back(0);
unsigned int sumBlocks = 0;
for (auto rowB = duneB_.begin(); rowB != duneB_.end(); ++rowB) {
int sizeRow = 0;
for (auto colB = rowB->begin(), endB = rowB->end(); colB != endB; ++colB) {
Bcols.emplace_back(colB.index());
for (int i = 0; i < numWellEq; ++i) {
for (int j = 0; j < numEq; ++j) {
Bvals.emplace_back((*colB)[i][j]);
}
}
sizeRow++;
}
sumBlocks += sizeRow;
Brows.emplace_back(sumBlocks);
}
wellContribs.addMultisegmentWellContribution(numEq,
numWellEq,
Mb,
Bvals,
Bcols,
Brows,
DnumBlocks,
Dvals,
Dcols,
Drows,
Cvals);
}
template<class Scalar, int numWellEq, int numEq>
template<class SparseMatrixAdapter>
void MultisegmentWellEquations<Scalar,numWellEq,numEq>::
extract(SparseMatrixAdapter& jacobian) const
{
const auto invDuneD = mswellhelpers::invertWithUMFPack<BVectorWell>(numWellEq,
numEq,
*duneDSolver_);
// We need to change matrix A as follows
// A -= C^T D^-1 B
// D is a (nseg x nseg) block matrix with (4 x 4) blocks.
// B and C are (nseg x ncells) block matrices with (4 x 4 blocks).
// They have nonzeros at (i, j) only if this well has a
// perforation at cell j connected to segment i. The code
// assumes that no cell is connected to more than one segment,
// i.e. the columns of B/C have no more than one nonzero.
for (size_t rowC = 0; rowC < duneC_.N(); ++rowC) {
for (auto colC = duneC_[rowC].begin(),
endC = duneC_[rowC].end(); colC != endC; ++colC) {
const auto row_index = colC.index();
for (size_t rowB = 0; rowB < duneB_.N(); ++rowB) {
for (auto colB = duneB_[rowB].begin(),
endB = duneB_[rowB].end(); colB != endB; ++colB) {
const auto col_index = colB.index();
OffDiagMatrixBlockWellType tmp1;
detail::multMatrixImpl(invDuneD[rowC][rowB], (*colB), tmp1, std::true_type());
typename SparseMatrixAdapter::MatrixBlock tmp2;
detail::multMatrixTransposedImpl((*colC), tmp1, tmp2, std::false_type());
jacobian.addToBlock(row_index, col_index, tmp2);
}
}
}
}
}
template<class Scalar, int numWellEq, int numEq>
template<class PressureMatrix>
void MultisegmentWellEquations<Scalar,numWellEq,numEq>::
extractCPRPressureMatrix(PressureMatrix& jacobian,
const BVector& weights,
const int pressureVarIndex,
const bool /*use_well_weights*/,
const WellInterfaceGeneric& well,
const int seg_pressure_var_ind,
const WellState& well_state) const
{
// Add the pressure contribution to the cpr system for the well
// Add for coupling from well to reservoir
const int welldof_ind = duneC_.M() + well.indexOfWell();
if (!well.isPressureControlled(well_state)) {
for (size_t rowC = 0; rowC < duneC_.N(); ++rowC) {
for (auto colC = duneC_[rowC].begin(),
endC = duneC_[rowC].end(); colC != endC; ++colC) {
const auto row_index = colC.index();
const auto& bw = weights[row_index];
double matel = 0.0;
for(size_t i = 0; i< bw.size(); ++i){
matel += bw[i]*(*colC)[seg_pressure_var_ind][i];
}
jacobian[row_index][welldof_ind] += matel;
}
}
}
// make cpr weights for well by pure avarage of reservoir weights of the perforations
if (!well.isPressureControlled(well_state)) {
auto well_weight = weights[0];
well_weight = 0.0;
int num_perfs = 0;
for (size_t rowB = 0; rowB < duneB_.N(); ++rowB) {
for (auto colB = duneB_[rowB].begin(),
endB = duneB_[rowB].end(); colB != endB; ++colB) {
const auto col_index = colB.index();
const auto& bw = weights[col_index];
well_weight += bw;
num_perfs += 1;
}
}
well_weight /= num_perfs;
assert(num_perfs > 0);
// Add for coupling from reservoir to well and caclulate diag elelement corresping to incompressible standard well
double diag_ell = 0.0;
for (size_t rowB = 0; rowB < duneB_.N(); ++rowB) {
const auto& bw = well_weight;
for (auto colB = duneB_[rowB].begin(),
endB = duneB_[rowB].end(); colB != endB; ++colB) {
const auto col_index = colB.index();
double matel = 0.0;
for(size_t i = 0; i< bw.size(); ++i){
matel += bw[i] *(*colB)[i][pressureVarIndex];
}
jacobian[welldof_ind][col_index] += matel;
diag_ell -= matel;
}
}
#define EXTRA_DEBUG_MSW 0
#if EXTRA_DEBUG_MSW
if (diag_ell <= 0.0) {
std::stringstream msg;
msg << "Diagonal element for cprw on "
<< this->name()
<< " is " << diag_ell;
OpmLog::warning(msg.str());
}
#endif
#undef EXTRA_DEBUG_MSW
jacobian[welldof_ind][welldof_ind] = diag_ell;
} else {
jacobian[welldof_ind][welldof_ind] = 1.0; // maybe we could have used diag_ell if calculated
}
}
#define INSTANCE(numWellEq, numEq) \
template class MultisegmentWellEquations<double,numWellEq,numEq>; \
template void MultisegmentWellEquations<double,numWellEq,numEq>:: \
extract(Linear::IstlSparseMatrixAdapter<MatrixBlock<double,numEq,numEq>>&) const; \
template void MultisegmentWellEquations<double,numWellEq,numEq>:: \
extractCPRPressureMatrix(Dune::BCRSMatrix<MatrixBlock<double,1,1>>&, \
const MultisegmentWellEquations<double,numWellEq,numEq>::BVector&, \
const int, \
const bool, \
const WellInterfaceGeneric&, \
const int, \
const WellState&) const;
INSTANCE(2,1)
INSTANCE(2,2)
INSTANCE(2,6)
INSTANCE(3,2)
INSTANCE(3,3)
INSTANCE(3,4)
INSTANCE(4,3)
INSTANCE(4,4)
INSTANCE(4,5)
}

133
opm/simulators/wells/MultisegmentWellEquations.hpp Normal file
View File

@ -0,0 +1,133 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
Copyright 2016 - 2017 IRIS 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_MULTISEGMENTWELL_EQUATIONS_HEADER_INCLUDED
#define OPM_MULTISEGMENTWELL_EQUATIONS_HEADER_INCLUDED
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <memory>
namespace Dune {
template<class M> class UMFPack;
}
namespace Opm
{
template<class Scalar> class MultisegmentWellGeneric;
class WellContributions;
class WellInterfaceGeneric;
class WellState;
template<class Scalar, int numWellEq, int numEq>
class MultisegmentWellEquations
{
public:
// sparsity pattern for the matrices
// [A C^T [x = [ res
// B D ] x_well] res_well]
// the vector type for the res_well and x_well
using VectorBlockWellType = Dune::FieldVector<Scalar,numWellEq>;
using BVectorWell = Dune::BlockVector<VectorBlockWellType>;
using VectorBlockType = Dune::FieldVector<Scalar,numEq>;
using BVector = Dune::BlockVector<VectorBlockType>;
// the matrix type for the diagonal matrix D
using DiagMatrixBlockWellType = Dune::FieldMatrix<Scalar,numWellEq,numWellEq>;
using DiagMatWell = Dune::BCRSMatrix<DiagMatrixBlockWellType>;
// the matrix type for the non-diagonal matrix B and C^T
using OffDiagMatrixBlockWellType = Dune::FieldMatrix<Scalar,numWellEq,numEq>;
using OffDiagMatWell = Dune::BCRSMatrix<OffDiagMatrixBlockWellType>;
MultisegmentWellEquations(const MultisegmentWellGeneric<Scalar>& well);
//! \brief Setup sparsity pattern for the matrices.
//! \param num_cells Total number of cells
//! \param numPerfs Number of perforations
//! \param cells Cell indices for perforations
void init(const int num_cells,
const int numPerfs,
const std::vector<int>& cells);
//! \brief Set all coefficients to 0.
void clear();
//! \brief Apply linear operator to vector.
void apply(const BVector& x, BVector& Ax) const;
//! \brief Apply linear operator to vector.
void apply(BVector& r) const;
//! \brief Compute the LU-decomposition of D matrix.
void createSolver();
//! \brief Apply inverted D matrix to residual and return result.
BVectorWell solve() const;
//! \brief Recover well solution.
//! \details xw = inv(D)*(rw - C*x)
void recoverSolutionWell(const BVector& x, BVectorWell& xw) const;
//! \brief Add the matrices of this well to the WellContributions object.
void extract(WellContributions& wellContribs) const;
//! \brief Add the matrices of this well to the sparse matrix adapter.
template<class SparseMatrixAdapter>
void extract(SparseMatrixAdapter& jacobian) const;
//! \brief Extract CPR pressure matrix.
template<class PressureMatrix>
void extractCPRPressureMatrix(PressureMatrix& jacobian,
const BVector& weights,
const int pressureVarIndex,
const bool /*use_well_weights*/,
const WellInterfaceGeneric& well,
const int seg_pressure_var_ind,
const WellState& well_state) const;
// two off-diagonal matrices
OffDiagMatWell duneB_;
OffDiagMatWell duneC_;
// "diagonal" matrix for the well. It has offdiagonal entries for inlets and outlets.
DiagMatWell duneD_;
/// \brief solver for diagonal matrix
///
/// This is a shared_ptr as MultisegmentWell is copied in computeWellPotentials...
mutable std::shared_ptr<Dune::UMFPack<DiagMatWell> > duneDSolver_;
// residuals of the well equations
BVectorWell resWell_;
private:
const MultisegmentWellGeneric<Scalar>& well_; //!< Reference to well
};
}
#endif // OPM_MULTISEGMENTWELLWELL_EQUATIONS_HEADER_INCLUDED

202
opm/simulators/wells/MultisegmentWellEval.cpp
View File

@ -30,7 +30,6 @@
#include <opm/models/blackoil/blackoilonephaseindices.hh>
#include <opm/models/blackoil/blackoiltwophaseindices.hh>
#include <opm/simulators/linalg/bda/WellContributions.hpp>
#include <opm/simulators/timestepping/ConvergenceReport.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/wells/MSWellHelpers.hpp>
@ -57,6 +56,7 @@ MultisegmentWellEval<FluidSystem,Indices,Scalar>::
MultisegmentWellEval(WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif)
: MultisegmentWellGeneric<Scalar>(baseif)
, baseif_(baseif)
, linSys_(*this)
, upwinding_segments_(this->numberOfSegments(), 0)
, segment_densities_(this->numberOfSegments(), 0.0)
, segment_mass_rates_(this->numberOfSegments(), 0.0)
@ -74,67 +74,7 @@ void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
initMatrixAndVectors(const int num_cells)
{
duneB_.setBuildMode(OffDiagMatWell::row_wise);
duneC_.setBuildMode(OffDiagMatWell::row_wise);
duneD_.setBuildMode(DiagMatWell::row_wise);
// set the size and patterns for all the matrices and vectors
// [A C^T [x = [ res
// B D] x_well] res_well]
// calculatiing the NNZ for duneD_
// NNZ = number_of_segments + 2 * (number_of_inlets / number_of_outlets)
{
int nnz_d = this->numberOfSegments();
for (const std::vector<int>& inlets : this->segment_inlets_) {
nnz_d += 2 * inlets.size();
}
duneD_.setSize(this->numberOfSegments(), this->numberOfSegments(), nnz_d);
}
duneB_.setSize(this->numberOfSegments(), num_cells, baseif_.numPerfs());
duneC_.setSize(this->numberOfSegments(), num_cells, baseif_.numPerfs());
// we need to add the off diagonal ones
for (auto row = duneD_.createbegin(), end = duneD_.createend(); row != end; ++row) {
// the number of the row corrspnds to the segment now
const int seg = row.index();
// adding the item related to outlet relation
const Segment& segment = this->segmentSet()[seg];
const int outlet_segment_number = segment.outletSegment();
if (outlet_segment_number > 0) { // if there is a outlet_segment
const int outlet_segment_index = this->segmentNumberToIndex(outlet_segment_number);
row.insert(outlet_segment_index);
}
// Add nonzeros for diagonal
row.insert(seg);
// insert the item related to its inlets
for (const int& inlet : this->segment_inlets_[seg]) {
row.insert(inlet);
}
}
// make the C matrix
for (auto row = duneC_.createbegin(), end = duneC_.createend(); row != end; ++row) {
// the number of the row corresponds to the segment number now.
for (const int& perf : this->segment_perforations_[row.index()]) {
const int cell_idx = baseif_.cells()[perf];
row.insert(cell_idx);
}
}
// make the B^T matrix
for (auto row = duneB_.createbegin(), end = duneB_.createend(); row != end; ++row) {
// the number of the row corresponds to the segment number now.
for (const int& perf : this->segment_perforations_[row.index()]) {
const int cell_idx = baseif_.cells()[perf];
row.insert(cell_idx);
}
}
resWell_.resize(this->numberOfSegments());
linSys_.init(num_cells, baseif_.numPerfs(), baseif_.cells());
primary_variables_.resize(this->numberOfSegments());
primary_variables_evaluation_.resize(this->numberOfSegments());
}
@ -172,7 +112,7 @@ getWellConvergence(const WellState& well_state,
std::vector<std::vector<double>> abs_residual(this->numberOfSegments(), std::vector<double>(numWellEq, 0.0));
for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
abs_residual[seg][eq_idx] = std::abs(resWell_[seg][eq_idx]);
abs_residual[seg][eq_idx] = std::abs(linSys_.resWell_[seg][eq_idx]);
}
}
@ -237,7 +177,7 @@ getWellConvergence(const WellState& well_state,
tolerance_wells,
tolerance_wells,
max_residual_allowed},
std::abs(resWell_[0][SPres]),
std::abs(linSys_.resWell_[0][SPres]),
report,
deferred_logger);
@ -448,20 +388,6 @@ updatePrimaryVariables(const WellState& well_state) const
}
}
template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
recoverSolutionWell(const BVector& x, BVectorWell& xw) const
{
if (!baseif_.isOperableAndSolvable() && !baseif_.wellIsStopped()) return;
BVectorWell resWell = resWell_;
// resWell = resWell - B * x
duneB_.mmv(x, resWell);
// xw = D^-1 * resWell
xw = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell);
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
@ -1289,9 +1215,9 @@ assembleControlEq(const WellState& well_state,
}
// using control_eq to update the matrix and residuals
resWell_[0][SPres] = control_eq.value();
linSys_.resWell_[0][SPres] = control_eq.value();
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
duneD_[0][0][SPres][pv_idx] = control_eq.derivative(pv_idx + Indices::numEq);
linSys_.duneD_[0][0][SPres][pv_idx] = control_eq.derivative(pv_idx + Indices::numEq);
}
}
@ -1336,14 +1262,14 @@ handleAccelerationPressureLoss(const int seg,
auto& segments = well_state.well(baseif_.indexOfWell()).segments;
segments.pressure_drop_accel[seg] = accelerationPressureLoss.value();
resWell_[seg][SPres] -= accelerationPressureLoss.value();
duneD_[seg][seg][SPres][SPres] -= accelerationPressureLoss.derivative(SPres + Indices::numEq);
duneD_[seg][seg][SPres][WQTotal] -= accelerationPressureLoss.derivative(WQTotal + Indices::numEq);
linSys_.resWell_[seg][SPres] -= accelerationPressureLoss.value();
linSys_.duneD_[seg][seg][SPres][SPres] -= accelerationPressureLoss.derivative(SPres + Indices::numEq);
linSys_.duneD_[seg][seg][SPres][WQTotal] -= accelerationPressureLoss.derivative(WQTotal + Indices::numEq);
if (has_wfrac_variable) {
duneD_[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + Indices::numEq);
linSys_.duneD_[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + Indices::numEq);
}
if (has_gfrac_variable) {
duneD_[seg][seg_upwind][SPres][GFrac] -= accelerationPressureLoss.derivative(GFrac + Indices::numEq);
linSys_.duneD_[seg][seg_upwind][SPres][GFrac] -= accelerationPressureLoss.derivative(GFrac + Indices::numEq);
}
}
@ -1374,24 +1300,24 @@ assembleDefaultPressureEq(const int seg,
segments.pressure_drop_friction[seg] = friction_pressure_drop.value();
}
resWell_[seg][SPres] = pressure_equation.value();
linSys_.resWell_[seg][SPres] = pressure_equation.value();
const int seg_upwind = upwinding_segments_[seg];
duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq);
duneD_[seg][seg][SPres][WQTotal] += pressure_equation.derivative(WQTotal + Indices::numEq);
linSys_.duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq);
linSys_.duneD_[seg][seg][SPres][WQTotal] += pressure_equation.derivative(WQTotal + Indices::numEq);
if (has_wfrac_variable) {
duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + Indices::numEq);
linSys_.duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + Indices::numEq);
}
if (has_gfrac_variable) {
duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + Indices::numEq);
linSys_.duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + Indices::numEq);
}
// contribution from the outlet segment
const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());
const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index);
resWell_[seg][SPres] -= outlet_pressure.value();
linSys_.resWell_[seg][SPres] -= outlet_pressure.value();
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + Indices::numEq);
linSys_.duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + Indices::numEq);
}
if (this->accelerationalPressureLossConsidered()) {
@ -1609,8 +1535,8 @@ assembleICDPressureEq(const int seg,
if (const auto& segment = this->segmentSet()[seg];
(segment.segmentType() == Segment::SegmentType::VALVE) &&
(segment.valve().status() == Opm::ICDStatus::SHUT) ) { // we use a zero rate equation to handle SHUT valve
resWell_[seg][SPres] = this->primary_variables_evaluation_[seg][WQTotal].value();
duneD_[seg][seg][SPres][WQTotal] = 1.;
linSys_.resWell_[seg][SPres] = this->primary_variables_evaluation_[seg][WQTotal].value();
linSys_.duneD_[seg][seg][SPres][WQTotal] = 1.;
auto& ws = well_state.well(baseif_.indexOfWell());
ws.segments.pressure_drop_friction[seg] = 0.;
@ -1644,23 +1570,23 @@ assembleICDPressureEq(const int seg,
ws.segments.pressure_drop_friction[seg] = icd_pressure_drop.value();
const int seg_upwind = upwinding_segments_[seg];
resWell_[seg][SPres] = pressure_equation.value();
duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq);
duneD_[seg][seg][SPres][WQTotal] += pressure_equation.derivative(WQTotal + Indices::numEq);
linSys_.resWell_[seg][SPres] = pressure_equation.value();
linSys_.duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq);
linSys_.duneD_[seg][seg][SPres][WQTotal] += pressure_equation.derivative(WQTotal + Indices::numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + Indices::numEq);
linSys_.duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + Indices::numEq);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + Indices::numEq);
linSys_.duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + Indices::numEq);
}
// contribution from the outlet segment
const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());
const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index);
resWell_[seg][SPres] -= outlet_pressure.value();
linSys_.resWell_[seg][SPres] -= outlet_pressure.value();
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + Indices::numEq);
linSys_.duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + Indices::numEq);
}
}
@ -1697,10 +1623,10 @@ getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
double residual = 0.;
if (eq_idx < baseif_.numComponents()) {
residual = std::abs(resWell_[seg][eq_idx]) * B_avg[eq_idx];
residual = std::abs(linSys_.resWell_[seg][eq_idx]) * B_avg[eq_idx];
} else {
if (seg > 0) {
residual = std::abs(resWell_[seg][eq_idx]);
residual = std::abs(linSys_.resWell_[seg][eq_idx]);
}
}
if (std::isnan(residual) || std::isinf(residual)) {
@ -1717,7 +1643,7 @@ getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
// handling the control equation residual
{
const double control_residual = std::abs(resWell_[0][numWellEq - 1]);
const double control_residual = std::abs(linSys_.resWell_[0][numWellEq - 1]);
if (std::isnan(control_residual) || std::isinf(control_residual)) {
deferred_logger.debug("nan or inf value for control residal get for well " + baseif_.name());
return {false, residuals};
@ -1858,74 +1784,6 @@ updateUpwindingSegments()
}
}
template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
addWellContribution(WellContributions& wellContribs) const
{
unsigned int Mb = duneB_.N(); // number of blockrows in duneB_, duneC_ and duneD_
unsigned int BnumBlocks = duneB_.nonzeroes();
unsigned int DnumBlocks = duneD_.nonzeroes();
// duneC
std::vector<unsigned int> Ccols;
std::vector<double> Cvals;
Ccols.reserve(BnumBlocks);
Cvals.reserve(BnumBlocks * Indices::numEq * numWellEq);
for (auto rowC = duneC_.begin(); rowC != duneC_.end(); ++rowC) {
for (auto colC = rowC->begin(), endC = rowC->end(); colC != endC; ++colC) {
Ccols.emplace_back(colC.index());
for (int i = 0; i < numWellEq; ++i) {
for (int j = 0; j < Indices::numEq; ++j) {
Cvals.emplace_back((*colC)[i][j]);
}
}
}
}
// duneD
Dune::UMFPack<DiagMatWell> umfpackMatrix(duneD_, 0);
double *Dvals = umfpackMatrix.getInternalMatrix().getValues();
auto *Dcols = umfpackMatrix.getInternalMatrix().getColStart();
auto *Drows = umfpackMatrix.getInternalMatrix().getRowIndex();
// duneB
std::vector<unsigned int> Bcols;
std::vector<unsigned int> Brows;
std::vector<double> Bvals;
Bcols.reserve(BnumBlocks);
Brows.reserve(Mb+1);
Bvals.reserve(BnumBlocks * Indices::numEq * numWellEq);
Brows.emplace_back(0);
unsigned int sumBlocks = 0;
for (auto rowB = duneB_.begin(); rowB != duneB_.end(); ++rowB) {
int sizeRow = 0;
for (auto colB = rowB->begin(), endB = rowB->end(); colB != endB; ++colB) {
Bcols.emplace_back(colB.index());
for (int i = 0; i < numWellEq; ++i) {
for (int j = 0; j < Indices::numEq; ++j) {
Bvals.emplace_back((*colB)[i][j]);
}
}
sizeRow++;
}
sumBlocks += sizeRow;
Brows.emplace_back(sumBlocks);
}
wellContribs.addMultisegmentWellContribution(Indices::numEq,
numWellEq,
Mb,
Bvals,
Bcols,
Brows,
DnumBlocks,
Dvals,
Dcols,
Drows,
Cvals);
}
#define INSTANCE(...) \
template class MultisegmentWellEval<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>;

55
opm/simulators/wells/MultisegmentWellEval.hpp
View File

@ -22,17 +22,13 @@
#ifndef OPM_MULTISEGMENTWELL_EVAL_HEADER_INCLUDED
#define OPM_MULTISEGMENTWELL_EVAL_HEADER_INCLUDED
#include <opm/simulators/wells/MultisegmentWellEquations.hpp>
#include <opm/simulators/wells/MultisegmentWellGeneric.hpp>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/input/eclipse/Schedule/Well/Well.hpp>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <array>
#include <memory>
#include <utility>
@ -55,10 +51,6 @@ class WellState;
template<typename FluidSystem, typename Indices, typename Scalar>
class MultisegmentWellEval : public MultisegmentWellGeneric<Scalar>
{
public:
/// add the contribution (C, D, B matrices) of this Well to the WellContributions object
void addWellContribution(WellContributions& wellContribs) const;
protected:
// TODO: for now, not considering the polymer, solvent and so on to simplify the development process.
@ -91,24 +83,10 @@ protected:
// the number of well equations TODO: it should have a more general strategy for it
static constexpr int numWellEq = Indices::numPhases + 1;
// sparsity pattern for the matrices
// [A C^T [x = [ res
// B D ] x_well] res_well]
using Equations = MultisegmentWellEquations<Scalar,numWellEq,Indices::numEq>;
// the vector type for the res_well and x_well
using VectorBlockWellType = Dune::FieldVector<Scalar, numWellEq>;
using BVectorWell = Dune::BlockVector<VectorBlockWellType>;
using VectorBlockType = Dune::FieldVector<Scalar, Indices::numEq>;
using BVector = Dune::BlockVector<VectorBlockType>;
// the matrix type for the diagonal matrix D
using DiagMatrixBlockWellType = Dune::FieldMatrix<Scalar, numWellEq, numWellEq>;
using DiagMatWell = Dune::BCRSMatrix<DiagMatrixBlockWellType>;
// the matrix type for the non-diagonal matrix B and C^T
using OffDiagMatrixBlockWellType = Dune::FieldMatrix<Scalar, numWellEq, Indices::numEq>;
using OffDiagMatWell = Dune::BCRSMatrix<OffDiagMatrixBlockWellType>;
using BVector = typename Equations::BVector;
using BVectorWell = typename Equations::BVectorWell;
// TODO: for more efficient implementation, we should have EvalReservoir, EvalWell, and EvalRerservoirAndWell
// EvalR (Eval), EvalW, EvalRW
@ -116,6 +94,12 @@ protected:
using EvalWell = DenseAd::Evaluation<double, /*size=*/Indices::numEq + numWellEq>;
using Eval = DenseAd::Evaluation<Scalar, /*size=*/Indices::numEq>;
public:
//! \brief Returns a const reference to equation system.
const Equations& linSys() const
{ return linSys_; }
protected:
MultisegmentWellEval(WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif);
void initMatrixAndVectors(const int num_cells);
@ -160,10 +144,6 @@ protected:
// handling the overshooting and undershooting of the fractions
void processFractions(const int seg) const;
// xw = inv(D)*(rw - C*x)
void recoverSolutionWell(const BVector& x,
BVectorWell& xw) const;
void updatePrimaryVariables(const WellState& well_state) const;
void updateUpwindingSegments();
@ -245,20 +225,7 @@ protected:
const WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif_;
// TODO, the following should go to a class for computing purpose
// two off-diagonal matrices
OffDiagMatWell duneB_;
OffDiagMatWell duneC_;
// "diagonal" matrix for the well. It has offdiagonal entries for inlets and outlets.
DiagMatWell duneD_;
/// \brief solver for diagonal matrix
///
/// This is a shared_ptr as MultisegmentWell is copied in computeWellPotentials...
mutable std::shared_ptr<Dune::UMFPack<DiagMatWell> > duneDSolver_;
// residuals of the well equations
BVectorWell resWell_;
Equations linSys_; //!< The equation system
// the values for the primary varibles
// based on different solutioin strategies, the wells can have different primary variables

15
opm/simulators/wells/MultisegmentWellGeneric.cpp
View File

@ -170,6 +170,21 @@ compPressureDrop() const
return segmentSet().compPressureDrop();
}
template<typename Scalar>
const std::vector<std::vector<int>>&
MultisegmentWellGeneric<Scalar>::
segmentInlets() const
{
return segment_inlets_;
}
template<typename Scalar>
const std::vector<std::vector<int>>&
MultisegmentWellGeneric<Scalar>::
segmentPerforations() const
{
return segment_perforations_;
}
template<typename Scalar>
int

25
opm/simulators/wells/MultisegmentWellGeneric.hpp
View File

@ -40,19 +40,16 @@ class WellState;
template <typename Scalar>
class MultisegmentWellGeneric
{
protected:
MultisegmentWellGeneric(WellInterfaceGeneric& baseif);
public:
//! \brief Returns the inlet segments for each segment.
const std::vector<std::vector<int>>& segmentInlets() const;
// scale the segment rates and pressure based on well rates and bhp
void scaleSegmentRatesWithWellRates(WellState& well_state) const;
void scaleSegmentPressuresWithBhp(WellState& well_state) const;
//! \brief Returns the perforation index for each segment.
const std::vector<std::vector<int>>& segmentPerforations() const;
// get the WellSegments from the well_ecl_
const WellSegments& segmentSet() const;
// components of the pressure drop to be included
WellSegments::CompPressureDrop compPressureDrop() const;
// segment number is an ID of the segment, it is specified in the deck
// get the loation of the segment with a segment number in the segmentSet
int segmentNumberToIndex(const int segment_number) const;
@ -60,6 +57,16 @@ protected:
/// number of segments for this well
int numberOfSegments() const;
protected:
MultisegmentWellGeneric(WellInterfaceGeneric& baseif);
// scale the segment rates and pressure based on well rates and bhp
void scaleSegmentRatesWithWellRates(WellState& well_state) const;
void scaleSegmentPressuresWithBhp(WellState& well_state) const;
// components of the pressure drop to be included
WellSegments::CompPressureDrop compPressureDrop() const;
/// Detect oscillation or stagnation based on the residual measure history
void detectOscillations(const std::vector<double>& measure_history,
const int it,
@ -82,7 +89,7 @@ protected:
// belonging to this segment
std::vector<std::vector<int>> segment_perforations_;
// the inlet segments for each segment. It is for convinience and efficiency reason
// the inlet segments for each segment. It is for convenience and efficiency reason
std::vector<std::vector<int>> segment_inlets_;
std::vector<double> segment_depth_diffs_;

176
opm/simulators/wells/MultisegmentWell_impl.hpp
View File

@ -19,8 +19,6 @@
*/
#include <opm/simulators/linalg/SmallDenseMatrixUtils.hpp>
#include <opm/simulators/wells/MSWellHelpers.hpp>
#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/input/eclipse/Schedule/MSW/Valve.hpp>
@ -196,22 +194,16 @@ namespace Opm
MultisegmentWell<TypeTag>::
apply(const BVector& x, BVector& Ax) const
{
if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
return;
}
if ( this->param_.matrix_add_well_contributions_ )
{
if (this->param_.matrix_add_well_contributions_) {
// Contributions are already in the matrix itself
return;
}
BVectorWell Bx(this->duneB_.N());
this->duneB_.mv(x, Bx);
// invDBx = duneD^-1 * Bx_
const BVectorWell invDBx = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, Bx);
// Ax = Ax - duneC_^T * invDBx
this->duneC_.mmtv(invDBx,Ax);
this->linSys_.apply(x, Ax);
}
@ -223,12 +215,11 @@ namespace Opm
MultisegmentWell<TypeTag>::
apply(BVector& r) const
{
if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
return;
}
// invDrw_ = duneD^-1 * resWell_
const BVectorWell invDrw = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
// r = r - duneC_^T * invDrw
this->duneC_.mmtv(invDrw, r);
this->linSys_.apply(r);
}
@ -240,10 +231,12 @@ namespace Opm
WellState& well_state,
DeferredLogger& deferred_logger)
{
if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
return;
}
BVectorWell xw(1);
this->recoverSolutionWell(x, xw);
this->linSys_.recoverSolutionWell(x, xw);
updateWellState(xw, well_state, deferred_logger);
}
@ -526,7 +519,7 @@ namespace Opm
// We assemble the well equations, then we check the convergence,
// which is why we do not put the assembleWellEq here.
const BVectorWell dx_well = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
const BVectorWell dx_well = this->linSys_.solve();
updateWellState(dx_well, well_state, deferred_logger);
}
@ -727,31 +720,7 @@ namespace Opm
MultisegmentWell<TypeTag>::
addWellContributions(SparseMatrixAdapter& jacobian) const
{
const auto invDuneD = mswellhelpers::invertWithUMFPack<BVectorWell>(this->duneD_, this->duneDSolver_);
// We need to change matrix A as follows
// A -= C^T D^-1 B
// D is a (nseg x nseg) block matrix with (4 x 4) blocks.
// B and C are (nseg x ncells) block matrices with (4 x 4 blocks).
// They have nonzeros at (i, j) only if this well has a
// perforation at cell j connected to segment i. The code
// assumes that no cell is connected to more than one segment,
// i.e. the columns of B/C have no more than one nonzero.
for (size_t rowC = 0; rowC < this->duneC_.N(); ++rowC) {
for (auto colC = this->duneC_[rowC].begin(), endC = this->duneC_[rowC].end(); colC != endC; ++colC) {
const auto row_index = colC.index();
for (size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) {
for (auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) {
const auto col_index = colB.index();
OffDiagMatrixBlockWellType tmp1;
detail::multMatrixImpl(invDuneD[rowC][rowB], (*colB), tmp1, std::true_type());
typename SparseMatrixAdapter::MatrixBlock tmp2;
detail::multMatrixTransposedImpl((*colC), tmp1, tmp2, std::false_type());
jacobian.addToBlock(row_index, col_index, tmp2);
}
}
}
}
this->linSys_.extract(jacobian);
}
@ -761,76 +730,17 @@ namespace Opm
addWellPressureEquations(PressureMatrix& jacobian,
const BVector& weights,
const int pressureVarIndex,
const bool /*use_well_weights*/,
const bool use_well_weights,
const WellState& well_state) const
{
// Add the pressure contribution to the cpr system for the well
// Add for coupling from well to reservoir
const auto seg_pressure_var_ind = this->SPres;
const int welldof_ind = this->duneC_.M() + this->index_of_well_;
if(not(this->isPressureControlled(well_state))){
for (size_t rowC = 0; rowC < this->duneC_.N(); ++rowC) {
for (auto colC = this->duneC_[rowC].begin(), endC = this->duneC_[rowC].end(); colC != endC; ++colC) {
const auto row_index = colC.index();
const auto& bw = weights[row_index];
double matel = 0.0;
for(size_t i = 0; i< bw.size(); ++i){
matel += bw[i]*(*colC)[seg_pressure_var_ind][i];
}
jacobian[row_index][welldof_ind] += matel;
}
}
}
// make cpr weights for well by pure avarage of reservoir weights of the perforations
if(not(this->isPressureControlled(well_state))){
auto well_weight = weights[0];
well_weight = 0.0;
int num_perfs = 0;
for (size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) {
for (auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) {
const auto col_index = colB.index();
const auto& bw = weights[col_index];
well_weight += bw;
num_perfs += 1;
}
}
well_weight /= num_perfs;
assert(num_perfs>0);
// Add for coupling from reservoir to well and caclulate diag elelement corresping to incompressible standard well
double diag_ell = 0.0;
for (size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) {
const auto& bw = well_weight;
for (auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) {
const auto col_index = colB.index();
double matel = 0.0;
for(size_t i = 0; i< bw.size(); ++i){
matel += bw[i] *(*colB)[i][pressureVarIndex];
}
jacobian[welldof_ind][col_index] += matel;
diag_ell -= matel;
}
}
#define EXTRA_DEBUG_MSW 0
#if EXTRA_DEBUG_MSW
if(not(diag_ell > 0.0)){
std::stringstream msg;
msg << "Diagonal element for cprw on "
<< this->name()
<< " is " << diag_ell;
OpmLog::warning(msg.str());
}
#endif
#undef EXTRA_DEBUG_MSW
jacobian[welldof_ind][welldof_ind] = diag_ell;
}else{
jacobian[welldof_ind][welldof_ind] = 1.0; // maybe we could have used diag_ell if calculated
}
this->linSys_.extractCPRPressureMatrix(jacobian,
weights,
pressureVarIndex,
use_well_weights,
*this,
this->SPres,
well_state);
}
@ -1494,7 +1404,7 @@ namespace Opm
assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
const BVectorWell dx_well = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
const BVectorWell dx_well = this->linSys_.solve();
if (it > this->param_.strict_inner_iter_wells_) {
relax_convergence = true;
@ -1622,13 +1532,7 @@ namespace Opm
computeSegmentFluidProperties(ebosSimulator, deferred_logger);
// clear all entries
this->duneB_ = 0.0;
this->duneC_ = 0.0;
this->duneD_ = 0.0;
this->resWell_ = 0.0;
this->duneDSolver_.reset();
this->linSys_.clear();
auto& ws = well_state.well(this->index_of_well_);
ws.dissolved_gas_rate = 0;
@ -1659,9 +1563,9 @@ namespace Opm
const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->surfaceVolumeFraction(seg, comp_idx)
- segment_fluid_initial_[seg][comp_idx]) / dt;
this->resWell_[seg][comp_idx] += accumulation_term.value();
this->linSys_.resWell_[seg][comp_idx] += accumulation_term.value();
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
this->duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + Indices::numEq);
this->linSys_.duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + Indices::numEq);
}
}
}
@ -1673,13 +1577,13 @@ namespace Opm
const int seg_upwind = this->upwinding_segments_[seg];
// segment_rate contains the derivatives with respect to WQTotal in seg,
// and WFrac and GFrac in seg_upwind
this->resWell_[seg][comp_idx] -= segment_rate.value();
this->duneD_[seg][seg][comp_idx][WQTotal] -= segment_rate.derivative(WQTotal + Indices::numEq);
this->linSys_.resWell_[seg][comp_idx] -= segment_rate.value();
this->linSys_.duneD_[seg][seg][comp_idx][WQTotal] -= segment_rate.derivative(WQTotal + Indices::numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
this->duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + Indices::numEq);
this->linSys_.duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + Indices::numEq);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
this->duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + Indices::numEq);
this->linSys_.duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + Indices::numEq);
}
// pressure derivative should be zero
}
@ -1694,13 +1598,13 @@ namespace Opm
const int inlet_upwind = this->upwinding_segments_[inlet];
// inlet_rate contains the derivatives with respect to WQTotal in inlet,
// and WFrac and GFrac in inlet_upwind
this->resWell_[seg][comp_idx] += inlet_rate.value();
this->duneD_[seg][inlet][comp_idx][WQTotal] += inlet_rate.derivative(WQTotal + Indices::numEq);
this->linSys_.resWell_[seg][comp_idx] += inlet_rate.value();
this->linSys_.duneD_[seg][inlet][comp_idx][WQTotal] += inlet_rate.derivative(WQTotal + Indices::numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
this->duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + Indices::numEq);
this->linSys_.duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + Indices::numEq);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
this->duneD_[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + Indices::numEq);
this->linSys_.duneD_[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + Indices::numEq);
}
// pressure derivative should be zero
}
@ -1744,21 +1648,21 @@ namespace Opm
this->connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective);
// subtract sum of phase fluxes in the well equations.
this->resWell_[seg][comp_idx] += cq_s_effective.value();
this->linSys_.resWell_[seg][comp_idx] += cq_s_effective.value();
// assemble the jacobians
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
// also need to consider the efficiency factor when manipulating the jacobians.
this->duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + Indices::numEq); // intput in transformed matrix
this->linSys_.duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + Indices::numEq); // intput in transformed matrix
// the index name for the D should be eq_idx / pv_idx
this->duneD_[seg][seg][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx + Indices::numEq);
this->linSys_.duneD_[seg][seg][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx + Indices::numEq);
}
for (int pv_idx = 0; pv_idx < Indices::numEq; ++pv_idx) {
// also need to consider the efficiency factor when manipulating the jacobians.
this->duneB_[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx);
this->linSys_.duneB_[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx);
}
}
}
@ -1780,6 +1684,8 @@ namespace Opm
this->assemblePressureEq(seg, unit_system, well_state, deferred_logger);
}
}
this->linSys_.createSolver();
}