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 #4268 from akva2/stdwell_eq

Browse Source 
added: StandardWellEquations
This commit is contained in:
Bård Skaflestad
2022-11-22 11:07:04 +01:00
committed by GitHub
parent 58d8ca144e d2c028fb74
commit c4fdbedf93
13 changed files with 609 additions and 332 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
CMakeLists_files.cmake
View File

@@ -96,6 +96,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/simulators/wells/PerfData.cpp
opm/simulators/wells/SegmentState.cpp
opm/simulators/wells/SingleWellState.cpp
opm/simulators/wells/StandardWellEquations.cpp
opm/simulators/wells/StandardWellEval.cpp
opm/simulators/wells/StandardWellGeneric.cpp
opm/simulators/wells/TargetCalculator.cpp
@@ -380,6 +381,8 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/wells/SingleWellState.hpp
opm/simulators/wells/StandardWell.hpp
opm/simulators/wells/StandardWell_impl.hpp
opm/simulators/wells/StandardWellEquations.hpp
opm/simulators/wells/StandardWellEval.hpp
opm/simulators/wells/TargetCalculator.hpp
opm/simulators/wells/VFPHelpers.hpp
opm/simulators/wells/VFPInjProperties.hpp

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

@@ -320,7 +320,7 @@ namespace Opm {
Simulator& ebosSimulator_;
// a vector of all the wells.
std::vector<WellInterfacePtr > well_container_{};
std::vector<WellInterfacePtr> well_container_{};
std::vector<bool> is_cell_perforated_{};

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

@@ -1203,7 +1203,7 @@ namespace Opm {
auto& well = well_container_[i];
std::shared_ptr<StandardWell<TypeTag> > derived = std::dynamic_pointer_cast<StandardWell<TypeTag> >(well);
if (derived) {
wellContribs.addNumBlocks(derived->getNumBlocks());
wellContribs.addNumBlocks(derived->linSys().getNumBlocks());
}
}
@@ -1213,9 +1213,9 @@ namespace Opm {
for(unsigned int i = 0; i < well_container_.size(); i++){
auto& well = well_container_[i];
// maybe WellInterface could implement addWellContribution()
auto derived_std = std::dynamic_pointer_cast<StandardWell<TypeTag> >(well);
auto derived_std = std::dynamic_pointer_cast<StandardWell<TypeTag>>(well);
if (derived_std) {
derived_std->addWellContribution(wellContribs);
derived_std->linSys().extract(derived_std->numStaticWellEq, wellContribs);
} else {
auto derived_ms = std::dynamic_pointer_cast<MultisegmentWell<TypeTag> >(well);
if (derived_ms) {

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

@@ -19,6 +19,7 @@
*/
#include <opm/simulators/linalg/SmallDenseMatrixUtils.hpp>
#include <opm/simulators/wells/MSWellHelpers.hpp>
#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>

3
opm/simulators/wells/StandardWell.hpp
View File

@@ -253,9 +253,6 @@ namespace Opm
protected:
bool regularize_;
// xw = inv(D)*(rw - C*x)
void recoverSolutionWell(const BVector& x, BVectorWell& xw) const;
// updating the well_state based on well solution dwells
void updateWellState(const BVectorWell& dwells,
WellState& well_state,

416
opm/simulators/wells/StandardWellEquations.cpp Normal file
View File

@@ -0,0 +1,416 @@
/*
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/StandardWellEquations.hpp>
#include <opm/common/Exceptions.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/WellInterfaceGeneric.hpp>
#include <algorithm>
#include <cmath>
namespace Opm
{
template<class Scalar, int numEq>
StandardWellEquations<Scalar,numEq>::
StandardWellEquations(const ParallelWellInfo& parallel_well_info)
: parallelB_(duneB_, parallel_well_info)
{
duneB_.setBuildMode(OffDiagMatWell::row_wise);
duneC_.setBuildMode(OffDiagMatWell::row_wise),
invDuneD_.setBuildMode(DiagMatWell::row_wise);
}
template<class Scalar, int numEq>
void StandardWellEquations<Scalar,numEq>::
init(const int num_cells,
const int numWellEq,
const int numPerfs,
const std::vector<int>& cells)
{
// setup sparsity pattern for the matrices
//[A C^T [x = [ res
// B D] x_well] res_well]
// set the size of the matrices
duneD_.setSize(1, 1, 1);
duneB_.setSize(1, num_cells, numPerfs);
duneC_.setSize(1, num_cells, numPerfs);
for (auto row = duneD_.createbegin(),
end = duneD_.createend(); row != end; ++row) {
// Add nonzeros for diagonal
row.insert(row.index());
}
// the block size is run-time determined now
duneD_[0][0].resize(numWellEq, numWellEq);
for (auto row = duneB_.createbegin(),
end = duneB_.createend(); row != end; ++row) {
for (int perf = 0 ; perf < numPerfs; ++perf) {
const int cell_idx = cells[perf];
row.insert(cell_idx);
}
}
for (int perf = 0 ; perf < numPerfs; ++perf) {
const int cell_idx = cells[perf];
// the block size is run-time determined now
duneB_[0][cell_idx].resize(numWellEq, numEq);
}
// make the C^T matrix
for (auto row = duneC_.createbegin(),
end = duneC_.createend(); row != end; ++row) {
for (int perf = 0; perf < numPerfs; ++perf) {
const int cell_idx = cells[perf];
row.insert(cell_idx);
}
}
for (int perf = 0; perf < numPerfs; ++perf) {
const int cell_idx = cells[perf];
duneC_[0][cell_idx].resize(numWellEq, numEq);
}
resWell_.resize(1);
// the block size of resWell_ is also run-time determined now
resWell_[0].resize(numWellEq);
// resize temporary class variables
Bx_.resize(duneB_.N());
for (unsigned i = 0; i < duneB_.N(); ++i) {
Bx_[i].resize(numWellEq);
}
invDrw_.resize(duneD_.N());
for (unsigned i = 0; i < duneD_.N(); ++i) {
invDrw_[i].resize(numWellEq);
}
}
template<class Scalar, int numEq>
void StandardWellEquations<Scalar,numEq>::clear()
{
duneB_ = 0.0;
duneC_ = 0.0;
duneD_ = 0.0;
resWell_ = 0.0;
}
template<class Scalar, int numEq>
void StandardWellEquations<Scalar,numEq>::apply(const BVector& x, BVector& Ax) const
{
assert(Bx_.size() == duneB_.N());
assert(invDrw_.size() == invDuneD_.N());
// Bx_ = duneB_ * x
parallelB_.mv(x, Bx_);
// invDBx = invDuneD_ * Bx_
// TODO: with this, we modified the content of the invDrw_.
// Is it necessary to do this to save some memory?
auto& invDBx = invDrw_;
invDuneD_.mv(Bx_, invDBx);
// Ax = Ax - duneC_^T * invDBx
duneC_.mmtv(invDBx, Ax);
}
template<class Scalar, int numEq>
void StandardWellEquations<Scalar,numEq>::apply(BVector& r) const
{
assert(invDrw_.size() == invDuneD_.N());
// invDrw_ = invDuneD_ * resWell_
invDuneD_.mv(resWell_, invDrw_);
// r = r - duneC_^T * invDrw_
duneC_.mmtv(invDrw_, r);
}
template<class Scalar, int numEq>
void StandardWellEquations<Scalar,numEq>::invert()
{
try {
invDuneD_ = duneD_; // Not strictly need if not cpr with well contributions is used
detail::invertMatrix(invDuneD_[0][0]);
} catch (NumericalProblem&) {
// for singular matrices, use identity as the inverse
invDuneD_[0][0] = 0.0;
for (size_t i = 0; i < invDuneD_[0][0].rows(); ++i) {
invDuneD_[0][0][i][i] = 1.0;
}
}
}
template<class Scalar, int numEq>
void StandardWellEquations<Scalar,numEq>::solve(BVectorWell& dx_well) const
{
invDuneD_.mv(resWell_, dx_well);
}
template<class Scalar, int numEq>
void StandardWellEquations<Scalar,numEq>::
recoverSolutionWell(const BVector& x, BVectorWell& xw) const
{
BVectorWell resWell = resWell_;
// resWell = resWell - B * x
parallelB_.mmv(x, resWell);
// xw = D^-1 * resWell
invDuneD_.mv(resWell, xw);
}
template<class Scalar, int numEq>
void StandardWellEquations<Scalar,numEq>::
extract(const int numStaticWellEq,
WellContributions& wellContribs) const
{
std::vector<int> colIndices;
std::vector<double> nnzValues;
colIndices.reserve(duneB_.nonzeroes());
nnzValues.reserve(duneB_.nonzeroes() * numStaticWellEq * numEq);
// duneC
for (auto colC = duneC_[0].begin(),
endC = duneC_[0].end(); colC != endC; ++colC )
{
colIndices.emplace_back(colC.index());
for (int i = 0; i < numStaticWellEq; ++i) {
for (int j = 0; j < numEq; ++j) {
nnzValues.emplace_back((*colC)[i][j]);
}
}
}
wellContribs.addMatrix(WellContributions::MatrixType::C,
colIndices.data(), nnzValues.data(), duneC_.nonzeroes());
// invDuneD
colIndices.clear();
nnzValues.clear();
colIndices.emplace_back(0);
for (int i = 0; i < numStaticWellEq; ++i)
{
for (int j = 0; j < numStaticWellEq; ++j) {
nnzValues.emplace_back(invDuneD_[0][0][i][j]);
}
}
wellContribs.addMatrix(WellContributions::MatrixType::D,
colIndices.data(), nnzValues.data(), 1);
// duneB
colIndices.clear();
nnzValues.clear();
for (auto colB = duneB_[0].begin(),
endB = duneB_[0].end(); colB != endB; ++colB )
{
colIndices.emplace_back(colB.index());
for (int i = 0; i < numStaticWellEq; ++i) {
for (int j = 0; j < numEq; ++j) {
nnzValues.emplace_back((*colB)[i][j]);
}
}
}
wellContribs.addMatrix(WellContributions::MatrixType::B,
colIndices.data(), nnzValues.data(), duneB_.nonzeroes());
}
template<class Scalar, int numEq>
template<class SparseMatrixAdapter>
void StandardWellEquations<Scalar,numEq>::
extract(SparseMatrixAdapter& jacobian) 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.
typename SparseMatrixAdapter::MatrixBlock tmpMat;
Dune::DynamicMatrix<Scalar> tmp;
for (auto colC = duneC_[0].begin(),
endC = duneC_[0].end(); colC != endC; ++colC)
{
const auto row_index = colC.index();
for (auto colB = duneB_[0].begin(),
endB = duneB_[0].end(); colB != endB; ++colB)
{
detail::multMatrix(invDuneD_[0][0], (*colB), tmp);
detail::negativeMultMatrixTransposed((*colC), tmp, tmpMat);
jacobian.addToBlock(row_index, colB.index(), tmpMat);
}
}
}
template<class Scalar, int numEq>
unsigned int StandardWellEquations<Scalar,numEq>::
getNumBlocks() const
{
return duneB_.nonzeroes();
}
template<class Scalar, int numEq>
template<class PressureMatrix>
void StandardWellEquations<Scalar,numEq>::
extractCPRPressureMatrix(PressureMatrix& jacobian,
const BVector& weights,
const int pressureVarIndex,
const bool use_well_weights,
const WellInterfaceGeneric& well,
const int bhp_var_index,
const WellState& well_state) const
{
// This adds pressure quation for cpr
// For use_well_weights=true
// weights lamda = inv(D)'v v = 0 v(bhpInd) = 1
// the well equations are summed i lambda' B(:,pressureVarINd) -> B lambda'*D(:,bhpInd) -> D
// For use_well_weights = false
// weights lambda = \sum_i w /n where ths sum is over weights of all perforation cells
// in the case of pressure controlled trivial equations are used and bhp C=B=0
// then the flow part of the well equations are summed lambda'*B(1:n,pressureVarInd) -> B lambda'*D(1:n,bhpInd) -> D
// For bouth
// C -> w'C(:,bhpInd) where w is weights of the perforation cell
// Add the well contributions in cpr
// use_well_weights is a quasiimpes formulation which is not implemented in multisegment
int nperf = 0;
auto cell_weights = weights[0];// not need for not(use_well_weights)
cell_weights = 0.0;
assert(duneC_.M() == weights.size());
const int welldof_ind = duneC_.M() + well.indexOfWell();
// do not assume anything about pressure controlled with use_well_weights (work fine with the assumtion also)
if (!well.isPressureControlled(well_state) || use_well_weights) {
// make coupling for reservoir to well
for (auto colC = duneC_[0].begin(),
endC = duneC_[0].end(); colC != endC; ++colC) {
const auto row_ind = colC.index();
const auto& bw = weights[row_ind];
double matel = 0;
assert((*colC).M() == bw.size());
for (size_t i = 0; i < bw.size(); ++i) {
matel += (*colC)[bhp_var_index][i] * bw[i];
}
jacobian[row_ind][welldof_ind] = matel;
cell_weights += bw;
nperf += 1;
}
}
cell_weights /= nperf;
BVectorWell bweights(1);
size_t blockSz = duneD_[0][0].N();
bweights[0].resize(blockSz);
bweights[0] = 0.0;
double diagElem = 0;
if (use_well_weights ) {
// calculate weighs and set diagonal element
//NB! use this options without treating pressure controlled separated
//NB! calculate quasiimpes well weights NB do not work well with trueimpes reservoir weights
double abs_max = 0;
BVectorWell rhs(1);
rhs[0].resize(blockSz);
rhs[0][bhp_var_index] = 1.0;
DiagMatrixBlockWellType inv_diag_block = invDuneD_[0][0];
DiagMatrixBlockWellType inv_diag_block_transpose =
Opm::wellhelpers::transposeDenseDynMatrix(inv_diag_block);
for (size_t i = 0; i < blockSz; ++i) {
bweights[0][i] = 0;
for (size_t j = 0; j < blockSz; ++j) {
bweights[0][i] += inv_diag_block_transpose[i][j] * rhs[0][j];
}
abs_max = std::max(abs_max, std::fabs(bweights[0][i]));
}
assert(abs_max > 0.0);
for (size_t i = 0; i < blockSz; ++i) {
bweights[0][i] /= abs_max;
}
diagElem = 1.0 / abs_max;
} else {
// set diagonal element
if (well.isPressureControlled(well_state)) {
bweights[0][blockSz-1] = 1.0;
diagElem = 1.0; // better scaling could have used the calculation below if weights were calculated
} else {
for (size_t i = 0; i < cell_weights.size(); ++i) {
bweights[0][i] = cell_weights[i];
}
bweights[0][blockSz-1] = 0.0;
diagElem = 0.0;
const auto& locmat = duneD_[0][0];
for (size_t i = 0; i < cell_weights.size(); ++i) {
diagElem += locmat[i][bhp_var_index] * cell_weights[i];
}
}
}
//
jacobian[welldof_ind][welldof_ind] = diagElem;
// set the matrix elements for well reservoir coupling
if (!well.isPressureControlled(well_state) || use_well_weights) {
for (auto colB = duneB_[0].begin(),
endB = duneB_[0].end(); colB != endB; ++colB) {
const auto col_index = colB.index();
const auto& bw = bweights[0];
double matel = 0;
for (size_t i = 0; i < bw.size(); ++i) {
matel += (*colB)[i][pressureVarIndex] * bw[i];
}
jacobian[welldof_ind][col_index] = matel;
}
}
}
template<class Scalar, int numEq>
void StandardWellEquations<Scalar,numEq>::
sumDistributed(Parallel::Communication comm)
{
// accumulate resWell_ and duneD_ in parallel to get effects of all perforations (might be distributed)
wellhelpers::sumDistributedWellEntries(duneD_[0][0], resWell_[0], comm);
}
#define INSTANCE(N) \
template class StandardWellEquations<double,N>; \
template void StandardWellEquations<double,N>:: \
extract(Linear::IstlSparseMatrixAdapter<MatrixBlock<double,N,N>>&) const; \
template void StandardWellEquations<double,N>:: \
extractCPRPressureMatrix(Dune::BCRSMatrix<MatrixBlock<double,1,1>>&, \
const typename StandardWellEquations<double,N>::BVector&, \
const int, \
const bool, \
const WellInterfaceGeneric&, \
const int, \
const WellState&) const;
INSTANCE(1)
INSTANCE(2)
INSTANCE(3)
INSTANCE(4)
INSTANCE(5)
INSTANCE(6)
}

140
opm/simulators/wells/StandardWellEquations.hpp Normal file
View File

@@ -0,0 +1,140 @@
/*
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_STANDARDWELL_EQUATIONS_HEADER_INCLUDED
#define OPM_STANDARDWELL_EQUATIONS_HEADER_INCLUDED
#include <opm/simulators/utils/ParallelCommunication.hpp>
#include <opm/simulators/wells/WellHelpers.hpp>
#include <dune/common/dynmatrix.hh>
#include <dune/common/dynvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
namespace Opm
{
class ParallelWellInfo;
class WellContributions;
class WellInterfaceGeneric;
class WellState;
template<class Scalar, int numEq>
class StandardWellEquations
{
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::DynamicVector<Scalar>;
using BVectorWell = Dune::BlockVector<VectorBlockWellType>;
// the matrix type for the diagonal matrix D
using DiagMatrixBlockWellType = Dune::DynamicMatrix<Scalar>;
using DiagMatWell = Dune::BCRSMatrix<DiagMatrixBlockWellType>;
// the matrix type for the non-diagonal matrix B and C^T
using OffDiagMatrixBlockWellType = Dune::DynamicMatrix<Scalar>;
using OffDiagMatWell = Dune::BCRSMatrix<OffDiagMatrixBlockWellType>;
// block vector type
using BVector = Dune::BlockVector<Dune::FieldVector<Scalar,numEq>>;
StandardWellEquations(const ParallelWellInfo& parallel_well_info);
//! \brief Setup sparsity pattern for the matrices.
//! \param num_cells Total number of cells
//! \param numWellEq Number of well equations
//! \param numPerfs Number of perforations
//! \param cells Cell indices for perforations
void init(const int num_cells,
const int numWellEq,
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 Apply inverted D matrix to residual and store in vector.
void solve(BVectorWell& dx_well) const;
//! \brief Invert D matrix.
void invert();
//! \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(const int numStaticWellEq,
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 bhp_var_index,
const WellState& well_state) const;
//! \brief Get the number of blocks of the C and B matrices.
unsigned int getNumBlocks() const;
//! \brief Sum with off-process contribution.
void sumDistributed(Parallel::Communication comm);
// two off-diagonal matrices
OffDiagMatWell duneB_;
OffDiagMatWell duneC_;
// diagonal matrix for the well
DiagMatWell invDuneD_;
DiagMatWell duneD_;
// Wrapper for the parallel application of B for distributed wells
wellhelpers::ParallelStandardWellB<Scalar> parallelB_;
// residuals of the well equations
BVectorWell resWell_;
// several vector used in the matrix calculation
mutable BVectorWell Bx_;
mutable BVectorWell invDrw_;
};
}
#endif // OPM_STANDARDWELL_EQUATIONS_HEADER_INCLUDED

114
opm/simulators/wells/StandardWellEval.cpp
View File

@@ -53,6 +53,7 @@ StandardWellEval(const WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif)
: StandardWellGeneric<Scalar>(baseif)
, baseif_(baseif)
, F0_(numWellConservationEq)
, linSys_(baseif_.parallelWellInfo())
{
}
@@ -444,9 +445,9 @@ assembleControlEq(const WellState& well_state,
// using control_eq to update the matrix and residuals
// TODO: we should use a different index system for the well equations
this->resWell_[0][Bhp] = control_eq.value();
this->linSys_.resWell_[0][Bhp] = control_eq.value();
for (int pv_idx = 0; pv_idx < numWellEq_; ++pv_idx) {
this->duneD_[0][0][Bhp][pv_idx] = control_eq.derivative(pv_idx + Indices::numEq);
this->linSys_.duneD_[0][0][Bhp][pv_idx] = control_eq.derivative(pv_idx + Indices::numEq);
}
}
@@ -773,7 +774,7 @@ getWellConvergence(const WellState& well_state,
res.resize(numWellEq_);
for (int eq_idx = 0; eq_idx < numWellEq_; ++eq_idx) {
// magnitude of the residual matters
res[eq_idx] = std::abs(this->resWell_[0][eq_idx]);
res[eq_idx] = std::abs(this->linSys_.resWell_[0][eq_idx]);
}
std::vector<double> well_flux_residual(baseif_.numComponents());
@@ -810,7 +811,7 @@ getWellConvergence(const WellState& well_state,
WellConvergence(baseif_).
checkConvergenceControlEq(well_state,
{1.e3, 1.e4, 1.e-4, 1.e-6, maxResidualAllowed},
std::abs(this->resWell_[0][Bhp]),
std::abs(this->linSys_.resWell_[0][Bhp]),
report,
deferred_logger);
@@ -1039,109 +1040,8 @@ init(std::vector<double>& perf_depth,
primary_variables_evaluation_.resize(numWellEq_, EvalWell{numWellEq_ + Indices::numEq, 0.0});
// setup sparsity pattern for the matrices
//[A C^T [x = [ res
// B D] x_well] res_well]
// set the size of the matrices
this->duneD_.setSize(1, 1, 1);
this->duneB_.setSize(1, num_cells, baseif_.numPerfs());
this->duneC_.setSize(1, num_cells, baseif_.numPerfs());
for (auto row=this->duneD_.createbegin(), end = this->duneD_.createend(); row!=end; ++row) {
// Add nonzeros for diagonal
row.insert(row.index());
}
// the block size is run-time determined now
this->duneD_[0][0].resize(numWellEq_, numWellEq_);
for (auto row = this->duneB_.createbegin(), end = this->duneB_.createend(); row!=end; ++row) {
for (int perf = 0 ; perf < baseif_.numPerfs(); ++perf) {
const int cell_idx = baseif_.cells()[perf];
row.insert(cell_idx);
}
}
for (int perf = 0 ; perf < baseif_.numPerfs(); ++perf) {
const int cell_idx = baseif_.cells()[perf];
// the block size is run-time determined now
this->duneB_[0][cell_idx].resize(numWellEq_, Indices::numEq);
}
// make the C^T matrix
for (auto row = this->duneC_.createbegin(), end = this->duneC_.createend(); row!=end; ++row) {
for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
const int cell_idx = baseif_.cells()[perf];
row.insert(cell_idx);
}
}
for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
const int cell_idx = baseif_.cells()[perf];
this->duneC_[0][cell_idx].resize(numWellEq_, Indices::numEq);
}
this->resWell_.resize(1);
// the block size of resWell_ is also run-time determined now
this->resWell_[0].resize(numWellEq_);
// resize temporary class variables
this->Bx_.resize( this->duneB_.N() );
for (unsigned i = 0; i < this->duneB_.N(); ++i) {
this->Bx_[i].resize(numWellEq_);
}
this->invDrw_.resize( this->duneD_.N() );
for (unsigned i = 0; i < this->duneD_.N(); ++i) {
this->invDrw_[i].resize(numWellEq_);
}
}
template<class FluidSystem, class Indices, class Scalar>
void
StandardWellEval<FluidSystem,Indices,Scalar>::
addWellContribution(WellContributions& wellContribs) const
{
std::vector<int> colIndices;
std::vector<double> nnzValues;
colIndices.reserve(this->duneB_.nonzeroes());
nnzValues.reserve(this->duneB_.nonzeroes()*numStaticWellEq * Indices::numEq);
// duneC
for ( auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC )
{
colIndices.emplace_back(colC.index());
for (int i = 0; i < numStaticWellEq; ++i) {
for (int j = 0; j < Indices::numEq; ++j) {
nnzValues.emplace_back((*colC)[i][j]);
}
}
}
wellContribs.addMatrix(WellContributions::MatrixType::C, colIndices.data(), nnzValues.data(), this->duneC_.nonzeroes());
// invDuneD
colIndices.clear();
nnzValues.clear();
colIndices.emplace_back(0);
for (int i = 0; i < numStaticWellEq; ++i)
{
for (int j = 0; j < numStaticWellEq; ++j) {
nnzValues.emplace_back(this->invDuneD_[0][0][i][j]);
}
}
wellContribs.addMatrix(WellContributions::MatrixType::D, colIndices.data(), nnzValues.data(), 1);
// duneB
colIndices.clear();
nnzValues.clear();
for ( auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB )
{
colIndices.emplace_back(colB.index());
for (int i = 0; i < numStaticWellEq; ++i) {
for (int j = 0; j < Indices::numEq; ++j) {
nnzValues.emplace_back((*colB)[i][j]);
}
}
}
wellContribs.addMatrix(WellContributions::MatrixType::B, colIndices.data(), nnzValues.data(), this->duneB_.nonzeroes());
this->linSys_.init(num_cells, this->numWellEq_,
baseif_.numPerfs(), baseif_.cells());
}
#define INSTANCE(...) \

12
opm/simulators/wells/StandardWellEval.hpp
View File

@@ -23,6 +23,7 @@
#ifndef OPM_STANDARDWELL_EVAL_HEADER_INCLUDED
#define OPM_STANDARDWELL_EVAL_HEADER_INCLUDED
#include <opm/simulators/wells/StandardWellEquations.hpp>
#include <opm/simulators/wells/StandardWellGeneric.hpp>
#include <opm/material/densead/DynamicEvaluation.hpp>
@@ -52,7 +53,11 @@ protected:
static constexpr int numWellControlEq = 1;
// number of the well equations that will always be used
// based on the solution strategy, there might be other well equations be introduced
public:
static constexpr int numStaticWellEq = numWellConservationEq + numWellControlEq;
protected:
// the index for Bhp in primary variables and also the index of well control equation
// they both will be the last one in their respective system.
// TODO: we should have indices for the well equations and well primary variables separately
@@ -96,8 +101,9 @@ public:
using Eval = DenseAd::Evaluation<Scalar, Indices::numEq>;
using BVectorWell = typename StandardWellGeneric<Scalar>::BVectorWell;
/// add the contribution (C, D^-1, B matrices) of this Well to the WellContributions object
void addWellContribution(WellContributions& wellContribs) const;
//! \brief Returns a const reference to equation system.
const StandardWellEquations<Scalar,Indices::numEq>& linSys() const
{ return linSys_; }
protected:
StandardWellEval(const WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif);
@@ -190,6 +196,8 @@ protected:
// the saturations in the well bore under surface conditions at the beginning of the time step
std::vector<double> F0_;
StandardWellEquations<Scalar,Indices::numEq> linSys_; //!< Linear equation system
};
}

12
opm/simulators/wells/StandardWellGeneric.cpp
View File

@@ -52,11 +52,7 @@ StandardWellGeneric(const WellInterfaceGeneric& baseif)
: baseif_(baseif)
, perf_densities_(baseif_.numPerfs())
, perf_pressure_diffs_(baseif_.numPerfs())
, parallelB_(duneB_, baseif_.parallelWellInfo())
{
duneB_.setBuildMode(OffDiagMatWell::row_wise);
duneC_.setBuildMode(OffDiagMatWell::row_wise);
invDuneD_.setBuildMode(DiagMatWell::row_wise);
}
@@ -149,14 +145,6 @@ computeConnectionPressureDelta()
baseif_.parallelWellInfo().partialSumPerfValues(beg, end);
}
template<class Scalar>
unsigned int
StandardWellGeneric<Scalar>::
getNumBlocks() const
{
return duneB_.nonzeroes();
}
template class StandardWellGeneric<double>;
}

22
opm/simulators/wells/StandardWellGeneric.hpp
View File

@@ -65,11 +65,6 @@ protected:
using OffDiagMatrixBlockWellType = Dune::DynamicMatrix<Scalar>;
using OffDiagMatWell = Dune::BCRSMatrix<OffDiagMatrixBlockWellType>;
public:
/// get the number of blocks of the C and B matrices, used to allocate memory in a WellContributions object
unsigned int getNumBlocks() const;
protected:
StandardWellGeneric(const WellInterfaceGeneric& baseif);
// calculate a relaxation factor to avoid overshoot of total rates
@@ -85,28 +80,11 @@ protected:
// Base interface reference
const WellInterfaceGeneric& baseif_;
// residuals of the well equations
BVectorWell resWell_;
// densities of the fluid in each perforation
std::vector<double> perf_densities_;
// pressure drop between different perforations
std::vector<double> perf_pressure_diffs_;
// two off-diagonal matrices
OffDiagMatWell duneB_;
OffDiagMatWell duneC_;
// diagonal matrix for the well
DiagMatWell invDuneD_;
DiagMatWell duneD_;
// Wrapper for the parallel application of B for distributed wells
wellhelpers::ParallelStandardWellB<Scalar> parallelB_;
// several vector used in the matrix calculation
mutable BVectorWell Bx_;
mutable BVectorWell invDrw_;
double getRho() const
{
return this->perf_densities_.empty() ? 0.0 : perf_densities_[0];

209
opm/simulators/wells/StandardWell_impl.hpp
View File

@@ -19,9 +19,7 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/common/utility/numeric/RootFinders.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/linalg/SmallDenseMatrixUtils.hpp>
#include <opm/simulators/wells/VFPHelpers.hpp>
#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
#include <opm/simulators/wells/WellConvergence.hpp>
@@ -432,10 +430,7 @@ namespace Opm
if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
// clear all entries
this->duneB_ = 0.0;
this->duneC_ = 0.0;
this->duneD_ = 0.0;
this->resWell_ = 0.0;
this->linSys_.clear();
assembleWellEqWithoutIterationImpl(ebosSimulator, dt, well_state, group_state, deferred_logger);
}
@@ -488,17 +483,17 @@ namespace Opm
connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective);
// subtract sum of phase fluxes in the well equations.
this->resWell_[0][componentIdx] += cq_s_effective.value();
this->linSys_.resWell_[0][componentIdx] += cq_s_effective.value();
// assemble the jacobians
for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
// also need to consider the efficiency factor when manipulating the jacobians.
this->duneC_[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+Indices::numEq); // intput in transformed matrix
this->duneD_[0][0][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx+Indices::numEq);
this->linSys_.duneC_[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+Indices::numEq); // intput in transformed matrix
this->linSys_.duneD_[0][0][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx+Indices::numEq);
}
for (int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
this->duneB_[0][cell_idx][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx);
this->linSys_.duneB_[0][cell_idx][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx);
}
// Store the perforation phase flux for later usage.
@@ -512,7 +507,7 @@ namespace Opm
if constexpr (has_zFraction) {
for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
this->duneC_[0][cell_idx][pvIdx][Indices::contiZfracEqIdx] -= cq_s_zfrac_effective.derivative(pvIdx+Indices::numEq);
this->linSys_.duneC_[0][cell_idx][pvIdx][Indices::contiZfracEqIdx] -= cq_s_zfrac_effective.derivative(pvIdx+Indices::numEq);
}
}
}
@@ -529,8 +524,8 @@ namespace Opm
}
// accumulate resWell_ and duneD_ in parallel to get effects of all perforations (might be distributed)
wellhelpers::sumDistributedWellEntries(this->duneD_[0][0], this->resWell_[0],
this->parallel_well_info_.communication());
this->linSys_.sumDistributed(this->parallel_well_info_.communication());
// add vol * dF/dt + Q to the well equations;
for (int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
// TODO: following the development in MSW, we need to convert the volume of the wellbore to be surface volume
@@ -542,9 +537,9 @@ namespace Opm
}
resWell_loc -= this->getQs(componentIdx) * this->well_efficiency_factor_;
for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
this->duneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+Indices::numEq);
this->linSys_.duneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+Indices::numEq);
}
this->resWell_[0][componentIdx] += resWell_loc.value();
this->linSys_.resWell_[0][componentIdx] += resWell_loc.value();
}
const auto& summaryState = ebosSimulator.vanguard().summaryState();
@@ -554,14 +549,7 @@ namespace Opm
// do the local inversion of D.
try {
this->invDuneD_ = this->duneD_; // Not strictly need if not cpr with well contributions is used
detail::invertMatrix(this->invDuneD_[0][0]);
} catch (NumericalProblem&) {
// for singular matrices, use identity as the inverse
this->invDuneD_[0][0] = 0.0;
for (size_t i = 0; i < this->invDuneD_[0][0].rows(); ++i) {
this->invDuneD_[0][0][i][i] = 1.0;
}
this->linSys_.invert();
} catch( ... ) {
OPM_DEFLOG_THROW(NumericalIssue,"Error when inverting local well equations for well " + name(), deferred_logger);
}
@@ -1689,7 +1677,7 @@ namespace Opm
// which is why we do not put the assembleWellEq here.
BVectorWell dx_well(1);
dx_well[0].resize(this->numWellEq_);
this->invDuneD_.mv(this->resWell_, dx_well);
this->linSys_.solve( dx_well);
updateWellState(dx_well, well_state, deferred_logger);
}
@@ -1725,20 +1713,8 @@ namespace Opm
// Contributions are already in the matrix itself
return;
}
assert( this->Bx_.size() == this->duneB_.N() );
assert( this->invDrw_.size() == this->invDuneD_.N() );
// Bx_ = duneB_ * x
this->parallelB_.mv(x, this->Bx_);
// invDBx = invDuneD_ * Bx_
// TODO: with this, we modified the content of the invDrw_.
// Is it necessary to do this to save some memory?
BVectorWell& invDBx = this->invDrw_;
this->invDuneD_.mv(this->Bx_, invDBx);
// Ax = Ax - duneC_^T * invDBx
this->duneC_.mmtv(invDBx,Ax);
this->linSys_.apply(x, Ax);
}
@@ -1751,29 +1727,9 @@ namespace Opm
{
if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
assert( this->invDrw_.size() == this->invDuneD_.N() );
// invDrw_ = invDuneD_ * resWell_
this->invDuneD_.mv(this->resWell_, this->invDrw_);
// r = r - duneC_^T * invDrw_
this->duneC_.mmtv(this->invDrw_, r);
this->linSys_.apply(r);
}
template<typename TypeTag>
void
StandardWell<TypeTag>::
recoverSolutionWell(const BVector& x, BVectorWell& xw) const
{
if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
BVectorWell resWell = this->resWell_;
// resWell = resWell - B * x
this->parallelB_.mmv(x, resWell);
// xw = D^-1 * resWell
this->invDuneD_.mv(resWell, xw);
}
@@ -1789,7 +1745,7 @@ namespace Opm
BVectorWell xw(1);
xw[0].resize(this->numWellEq_);
recoverSolutionWell(x, xw);
this->linSys_.recoverSolutionWell(x, xw);
updateWellState(xw, well_state, deferred_logger);
}
@@ -2171,24 +2127,7 @@ namespace Opm
void
StandardWell<TypeTag>::addWellContributions(SparseMatrixAdapter& jacobian) 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.
typename SparseMatrixAdapter::MatrixBlock tmpMat;
Dune::DynamicMatrix<Scalar> tmp;
for ( auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC )
{
const auto row_index = colC.index();
for ( auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB )
{
detail::multMatrix(this->invDuneD_[0][0], (*colB), tmp);
detail::negativeMultMatrixTransposed((*colC), tmp, tmpMat);
jacobian.addToBlock( row_index, colB.index(), tmpMat );
}
}
this->linSys_.extract(jacobian);
}
@@ -2200,105 +2139,13 @@ namespace Opm
const bool use_well_weights,
const WellState& well_state) const
{
// This adds pressure quation for cpr
// For use_well_weights=true
// weights lamda = inv(D)'v v = 0 v(bhpInd) = 1
// the well equations are summed i lambda' B(:,pressureVarINd) -> B lambda'*D(:,bhpInd) -> D
// For use_well_weights = false
// weights lambda = \sum_i w /n where ths sum is over weights of all perforation cells
// in the case of pressure controlled trivial equations are used and bhp C=B=0
// then the flow part of the well equations are summed lambda'*B(1:n,pressureVarInd) -> B lambda'*D(1:n,bhpInd) -> D
// For bouth
// C -> w'C(:,bhpInd) where w is weights of the perforation cell
// Add the well contributions in cpr
// use_well_weights is a quasiimpes formulation which is not implemented in multisegment
int bhp_var_index = Bhp;
int nperf = 0;
auto cell_weights = weights[0];// not need for not(use_well_weights)
cell_weights = 0.0;
assert(this->duneC_.M() == weights.size());
const int welldof_ind = this->duneC_.M() + this->index_of_well_;
// do not assume anything about pressure controlled with use_well_weights (work fine with the assumtion also)
if( not( this->isPressureControlled(well_state) ) || use_well_weights ){
// make coupling for reservoir to well
for (auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC) {
const auto row_ind = colC.index();
const auto& bw = weights[row_ind];
double matel = 0;
assert((*colC).M() == bw.size());
for (size_t i = 0; i < bw.size(); ++i) {
matel += (*colC)[bhp_var_index][i] * bw[i];
}
jacobian[row_ind][welldof_ind] = matel;
cell_weights += bw;
nperf += 1;
}
}
cell_weights /= nperf;
BVectorWell bweights(1);
size_t blockSz = this->numWellEq_;
bweights[0].resize(blockSz);
bweights[0] = 0.0;
double diagElem = 0;
{
if ( use_well_weights ){
// calculate weighs and set diagonal element
//NB! use this options without treating pressure controlled separated
//NB! calculate quasiimpes well weights NB do not work well with trueimpes reservoir weights
double abs_max = 0;
BVectorWell rhs(1);
rhs[0].resize(blockSz);
rhs[0][bhp_var_index] = 1.0;
DiagMatrixBlockWellType inv_diag_block = this->invDuneD_[0][0];
DiagMatrixBlockWellType inv_diag_block_transpose = Opm::wellhelpers::transposeDenseDynMatrix(inv_diag_block);
for (size_t i = 0; i < blockSz; ++i) {
bweights[0][i] = 0;
for (size_t j = 0; j < blockSz; ++j) {
bweights[0][i] += inv_diag_block_transpose[i][j]*rhs[0][j];
}
abs_max = std::max(abs_max, std::fabs(bweights[0][i]));
}
assert( abs_max > 0.0 );
for (size_t i = 0; i < blockSz; ++i) {
bweights[0][i] /= abs_max;
}
diagElem = 1.0/abs_max;
}else{
// set diagonal element
if( this->isPressureControlled(well_state) ){
bweights[0][blockSz-1] = 1.0;
diagElem = 1.0;// better scaling could have used the calculation below if weights were calculated
}else{
for (size_t i = 0; i < cell_weights.size(); ++i) {
bweights[0][i] = cell_weights[i];
}
bweights[0][blockSz-1] = 0.0;
diagElem = 0.0;
const auto& locmat = this->duneD_[0][0];
for (size_t i = 0; i < cell_weights.size(); ++i) {
diagElem += locmat[i][bhp_var_index]*cell_weights[i];
}
}
}
}
//
jacobian[welldof_ind][welldof_ind] = diagElem;
// set the matrix elements for well reservoir coupling
if( not( this->isPressureControlled(well_state) ) || use_well_weights ){
for (auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB) {
const auto col_index = colB.index();
const auto& bw = bweights[0];
double matel = 0;
for (size_t i = 0; i < bw.size(); ++i) {
matel += (*colB)[i][pressureVarIndex] * bw[i];
}
jacobian[welldof_ind][col_index] = matel;
}
}
this->linSys_.extractCPRPressureMatrix(jacobian,
weights,
pressureVarIndex,
use_well_weights,
*this,
Bhp,
well_state);
}
@@ -2466,7 +2313,7 @@ namespace Opm
// equation for the water velocity
const EvalWell eq_wat_vel = this->primary_variables_evaluation_[wat_vel_index] - water_velocity;
this->resWell_[0][wat_vel_index] = eq_wat_vel.value();
this->linSys_.resWell_[0][wat_vel_index] = eq_wat_vel.value();
const auto& ws = well_state.well(this->index_of_well_);
const auto& perf_data = ws.perf_data;
@@ -2481,15 +2328,15 @@ namespace Opm
const EvalWell eq_pskin = this->primary_variables_evaluation_[pskin_index]
- pskin(throughput, this->primary_variables_evaluation_[wat_vel_index], poly_conc, deferred_logger);
this->resWell_[0][pskin_index] = eq_pskin.value();
this->linSys_.resWell_[0][pskin_index] = eq_pskin.value();
for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
this->duneD_[0][0][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx+Indices::numEq);
this->duneD_[0][0][pskin_index][pvIdx] = eq_pskin.derivative(pvIdx+Indices::numEq);
this->linSys_.duneD_[0][0][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx+Indices::numEq);
this->linSys_.duneD_[0][0][pskin_index][pvIdx] = eq_pskin.derivative(pvIdx+Indices::numEq);
}
// the water velocity is impacted by the reservoir primary varaibles. It needs to enter matrix B
for (int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
this->duneB_[0][cell_idx][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx);
this->linSys_.duneB_[0][cell_idx][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx);
}
}

1
opm/simulators/wells/WellInterface.hpp
View File

@@ -71,7 +71,6 @@ class WellInterface : public WellInterfaceIndices<GetPropType<TypeTag, Propertie
GetPropType<TypeTag, Properties::Scalar>>
{
public:
using ModelParameters = BlackoilModelParametersEbos<TypeTag>;
using Grid = GetPropType<TypeTag, Properties::Grid>;