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

add: MultisegmentWellEquations

Browse Source 
this is a container for the multisegment well equation system
This commit is contained in:
Arne Morten Kvarving
2022-12-01 09:23:52 +01:00
parent e23da0b6b7
commit 57f09050fc
7 changed files with 221 additions and 125 deletions
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/GlobalWellInfo.cpp
opm/simulators/wells/GroupState.cpp opm/simulators/wells/GroupState.cpp
opm/simulators/wells/MSWellHelpers.cpp opm/simulators/wells/MSWellHelpers.cpp
opm/simulators/wells/MultisegmentWellEquations.cpp
opm/simulators/wells/MultisegmentWellEval.cpp opm/simulators/wells/MultisegmentWellEval.cpp
opm/simulators/wells/MultisegmentWellGeneric.cpp opm/simulators/wells/MultisegmentWellGeneric.cpp
opm/simulators/wells/ParallelWellInfo.cpp opm/simulators/wells/ParallelWellInfo.cpp
@@ -371,6 +372,9 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/wells/MSWellHelpers.hpp opm/simulators/wells/MSWellHelpers.hpp
opm/simulators/wells/MultisegmentWell.hpp opm/simulators/wells/MultisegmentWell.hpp
opm/simulators/wells/MultisegmentWell_impl.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/ParallelWellInfo.hpp
opm/simulators/wells/PerfData.hpp opm/simulators/wells/PerfData.hpp
opm/simulators/wells/PerforationData.hpp opm/simulators/wells/PerforationData.hpp

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

@@ -65,10 +65,9 @@ namespace Opm
using typename Base::BVector; using typename Base::BVector;
using typename Base::Eval; using typename Base::Eval;
using typename MSWEval::Equations;
using typename MSWEval::EvalWell; using typename MSWEval::EvalWell;
using typename MSWEval::BVectorWell; using typename MSWEval::BVectorWell;
using typename MSWEval::DiagMatWell;
using typename MSWEval::OffDiagMatrixBlockWellType;
using MSWEval::GFrac; using MSWEval::GFrac;
using MSWEval::WFrac; using MSWEval::WFrac;
using MSWEval::WQTotal; using MSWEval::WQTotal;

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

@@ -0,0 +1,29 @@
/*
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>
namespace Opm
{
}

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

@@ -0,0 +1,80 @@
/*
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, 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>;
// 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_;
};
}
#endif // OPM_MULTISEGMENTWELLWELL_EQUATIONS_HEADER_INCLUDED

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

@@ -74,9 +74,9 @@ void
MultisegmentWellEval<FluidSystem,Indices,Scalar>:: MultisegmentWellEval<FluidSystem,Indices,Scalar>::
initMatrixAndVectors(const int num_cells) initMatrixAndVectors(const int num_cells)
{ {
duneB_.setBuildMode(OffDiagMatWell::row_wise); linSys_.duneB_.setBuildMode(Equations::OffDiagMatWell::row_wise);
duneC_.setBuildMode(OffDiagMatWell::row_wise); linSys_.duneC_.setBuildMode(Equations::OffDiagMatWell::row_wise);
duneD_.setBuildMode(DiagMatWell::row_wise); linSys_.duneD_.setBuildMode(Equations::DiagMatWell::row_wise);
// set the size and patterns for all the matrices and vectors // set the size and patterns for all the matrices and vectors
// [A C^T [x = [ res // [A C^T [x = [ res
@@ -89,13 +89,14 @@ initMatrixAndVectors(const int num_cells)
for (const std::vector<int>& inlets : this->segment_inlets_) { for (const std::vector<int>& inlets : this->segment_inlets_) {
nnz_d += 2 * inlets.size(); nnz_d += 2 * inlets.size();
} }
duneD_.setSize(this->numberOfSegments(), this->numberOfSegments(), nnz_d); linSys_.duneD_.setSize(this->numberOfSegments(), this->numberOfSegments(), nnz_d);
} }
duneB_.setSize(this->numberOfSegments(), num_cells, baseif_.numPerfs()); linSys_.duneB_.setSize(this->numberOfSegments(), num_cells, baseif_.numPerfs());
duneC_.setSize(this->numberOfSegments(), num_cells, baseif_.numPerfs()); linSys_.duneC_.setSize(this->numberOfSegments(), num_cells, baseif_.numPerfs());
// we need to add the off diagonal ones // we need to add the off diagonal ones
for (auto row = duneD_.createbegin(), end = duneD_.createend(); row != end; ++row) { for (auto row = linSys_.duneD_.createbegin(),
end = linSys_.duneD_.createend(); row != end; ++row) {
// the number of the row corrspnds to the segment now // the number of the row corrspnds to the segment now
const int seg = row.index(); const int seg = row.index();
// adding the item related to outlet relation // adding the item related to outlet relation
@@ -116,7 +117,8 @@ initMatrixAndVectors(const int num_cells)
} }
// make the C matrix // make the C matrix
for (auto row = duneC_.createbegin(), end = duneC_.createend(); row != end; ++row) { for (auto row = linSys_.duneC_.createbegin(),
end = linSys_.duneC_.createend(); row != end; ++row) {
// the number of the row corresponds to the segment number now. // the number of the row corresponds to the segment number now.
for (const int& perf : this->segment_perforations_[row.index()]) { for (const int& perf : this->segment_perforations_[row.index()]) {
const int cell_idx = baseif_.cells()[perf]; const int cell_idx = baseif_.cells()[perf];
@@ -125,7 +127,8 @@ initMatrixAndVectors(const int num_cells)
} }
// make the B^T matrix // make the B^T matrix
for (auto row = duneB_.createbegin(), end = duneB_.createend(); row != end; ++row) { for (auto row = linSys_.duneB_.createbegin(),
end = linSys_.duneB_.createend(); row != end; ++row) {
// the number of the row corresponds to the segment number now. // the number of the row corresponds to the segment number now.
for (const int& perf : this->segment_perforations_[row.index()]) { for (const int& perf : this->segment_perforations_[row.index()]) {
const int cell_idx = baseif_.cells()[perf]; const int cell_idx = baseif_.cells()[perf];
@@ -133,7 +136,7 @@ initMatrixAndVectors(const int num_cells)
} }
} }
resWell_.resize(this->numberOfSegments()); linSys_.resWell_.resize(this->numberOfSegments());
primary_variables_.resize(this->numberOfSegments()); primary_variables_.resize(this->numberOfSegments());
primary_variables_evaluation_.resize(this->numberOfSegments()); primary_variables_evaluation_.resize(this->numberOfSegments());
@@ -172,7 +175,7 @@ getWellConvergence(const WellState& well_state,
std::vector<std::vector<double>> abs_residual(this->numberOfSegments(), std::vector<double>(numWellEq, 0.0)); 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 seg = 0; seg < this->numberOfSegments(); ++seg) {
for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { 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 +240,7 @@ getWellConvergence(const WellState& well_state,
tolerance_wells, tolerance_wells,
tolerance_wells, tolerance_wells,
max_residual_allowed}, max_residual_allowed},
std::abs(resWell_[0][SPres]), std::abs(linSys_.resWell_[0][SPres]),
report, report,
deferred_logger); deferred_logger);
@@ -455,11 +458,11 @@ recoverSolutionWell(const BVector& x, BVectorWell& xw) const
{ {
if (!baseif_.isOperableAndSolvable() && !baseif_.wellIsStopped()) return; if (!baseif_.isOperableAndSolvable() && !baseif_.wellIsStopped()) return;
BVectorWell resWell = resWell_; BVectorWell resWell = linSys_.resWell_;
// resWell = resWell - B * x // resWell = resWell - B * x
duneB_.mmv(x, resWell); linSys_.duneB_.mmv(x, resWell);
// xw = D^-1 * resWell // xw = D^-1 * resWell
xw = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell); xw = mswellhelpers::applyUMFPack(linSys_.duneD_, linSys_.duneDSolver_, resWell);
} }
template<typename FluidSystem, typename Indices, typename Scalar> template<typename FluidSystem, typename Indices, typename Scalar>
@@ -1289,9 +1292,9 @@ assembleControlEq(const WellState& well_state,
} }
// using control_eq to update the matrix and residuals // 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) { 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 +1339,14 @@ handleAccelerationPressureLoss(const int seg,
auto& segments = well_state.well(baseif_.indexOfWell()).segments; auto& segments = well_state.well(baseif_.indexOfWell()).segments;
segments.pressure_drop_accel[seg] = accelerationPressureLoss.value(); segments.pressure_drop_accel[seg] = accelerationPressureLoss.value();
resWell_[seg][SPres] -= accelerationPressureLoss.value(); linSys_.resWell_[seg][SPres] -= accelerationPressureLoss.value();
duneD_[seg][seg][SPres][SPres] -= accelerationPressureLoss.derivative(SPres + Indices::numEq); linSys_.duneD_[seg][seg][SPres][SPres] -= accelerationPressureLoss.derivative(SPres + Indices::numEq);
duneD_[seg][seg][SPres][WQTotal] -= accelerationPressureLoss.derivative(WQTotal + Indices::numEq); linSys_.duneD_[seg][seg][SPres][WQTotal] -= accelerationPressureLoss.derivative(WQTotal + Indices::numEq);
if (has_wfrac_variable) { 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) { 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 +1377,24 @@ assembleDefaultPressureEq(const int seg,
segments.pressure_drop_friction[seg] = friction_pressure_drop.value(); 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]; const int seg_upwind = upwinding_segments_[seg];
duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq); linSys_.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][WQTotal] += pressure_equation.derivative(WQTotal + Indices::numEq);
if (has_wfrac_variable) { 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) { 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 // contribution from the outlet segment
const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment()); const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());
const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index); 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) { 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()) { if (this->accelerationalPressureLossConsidered()) {
@@ -1609,8 +1612,8 @@ assembleICDPressureEq(const int seg,
if (const auto& segment = this->segmentSet()[seg]; if (const auto& segment = this->segmentSet()[seg];
(segment.segmentType() == Segment::SegmentType::VALVE) && (segment.segmentType() == Segment::SegmentType::VALVE) &&
(segment.valve().status() == Opm::ICDStatus::SHUT) ) { // we use a zero rate equation to handle SHUT 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(); linSys_.resWell_[seg][SPres] = this->primary_variables_evaluation_[seg][WQTotal].value();
duneD_[seg][seg][SPres][WQTotal] = 1.; linSys_.duneD_[seg][seg][SPres][WQTotal] = 1.;
auto& ws = well_state.well(baseif_.indexOfWell()); auto& ws = well_state.well(baseif_.indexOfWell());
ws.segments.pressure_drop_friction[seg] = 0.; ws.segments.pressure_drop_friction[seg] = 0.;
@@ -1644,23 +1647,23 @@ assembleICDPressureEq(const int seg,
ws.segments.pressure_drop_friction[seg] = icd_pressure_drop.value(); ws.segments.pressure_drop_friction[seg] = icd_pressure_drop.value();
const int seg_upwind = upwinding_segments_[seg]; const int seg_upwind = upwinding_segments_[seg];
resWell_[seg][SPres] = pressure_equation.value(); linSys_.resWell_[seg][SPres] = pressure_equation.value();
duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + Indices::numEq); linSys_.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][WQTotal] += pressure_equation.derivative(WQTotal + Indices::numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { 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)) { 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 // contribution from the outlet segment
const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment()); const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());
const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index); 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) { 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 +1700,10 @@ getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
double residual = 0.; double residual = 0.;
if (eq_idx < baseif_.numComponents()) { 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 { } else {
if (seg > 0) { 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)) { if (std::isnan(residual) || std::isinf(residual)) {
@@ -1717,7 +1720,7 @@ getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
// handling the control equation residual // 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)) { if (std::isnan(control_residual) || std::isinf(control_residual)) {
deferred_logger.debug("nan or inf value for control residal get for well " + baseif_.name()); deferred_logger.debug("nan or inf value for control residal get for well " + baseif_.name());
return {false, residuals}; return {false, residuals};
@@ -1863,16 +1866,16 @@ void
MultisegmentWellEval<FluidSystem,Indices,Scalar>:: MultisegmentWellEval<FluidSystem,Indices,Scalar>::
addWellContribution(WellContributions& wellContribs) const addWellContribution(WellContributions& wellContribs) const
{ {
unsigned int Mb = duneB_.N(); // number of blockrows in duneB_, duneC_ and duneD_ unsigned int Mb = linSys_.duneB_.N(); // number of blockrows in duneB_, duneC_ and duneD_
unsigned int BnumBlocks = duneB_.nonzeroes(); unsigned int BnumBlocks = linSys_.duneB_.nonzeroes();
unsigned int DnumBlocks = duneD_.nonzeroes(); unsigned int DnumBlocks = linSys_.duneD_.nonzeroes();
// duneC // duneC
std::vector<unsigned int> Ccols; std::vector<unsigned int> Ccols;
std::vector<double> Cvals; std::vector<double> Cvals;
Ccols.reserve(BnumBlocks); Ccols.reserve(BnumBlocks);
Cvals.reserve(BnumBlocks * Indices::numEq * numWellEq); Cvals.reserve(BnumBlocks * Indices::numEq * numWellEq);
for (auto rowC = duneC_.begin(); rowC != duneC_.end(); ++rowC) { for (auto rowC = linSys_.duneC_.begin(); rowC != linSys_.duneC_.end(); ++rowC) {
for (auto colC = rowC->begin(), endC = rowC->end(); colC != endC; ++colC) { for (auto colC = rowC->begin(), endC = rowC->end(); colC != endC; ++colC) {
Ccols.emplace_back(colC.index()); Ccols.emplace_back(colC.index());
for (int i = 0; i < numWellEq; ++i) { for (int i = 0; i < numWellEq; ++i) {
@@ -1884,7 +1887,7 @@ addWellContribution(WellContributions& wellContribs) const
} }
// duneD // duneD
Dune::UMFPack<DiagMatWell> umfpackMatrix(duneD_, 0); Dune::UMFPack<typename Equations::DiagMatWell> umfpackMatrix(linSys_.duneD_, 0);
double *Dvals = umfpackMatrix.getInternalMatrix().getValues(); double *Dvals = umfpackMatrix.getInternalMatrix().getValues();
auto *Dcols = umfpackMatrix.getInternalMatrix().getColStart(); auto *Dcols = umfpackMatrix.getInternalMatrix().getColStart();
auto *Drows = umfpackMatrix.getInternalMatrix().getRowIndex(); auto *Drows = umfpackMatrix.getInternalMatrix().getRowIndex();
@@ -1898,7 +1901,7 @@ addWellContribution(WellContributions& wellContribs) const
Bvals.reserve(BnumBlocks * Indices::numEq * numWellEq); Bvals.reserve(BnumBlocks * Indices::numEq * numWellEq);
Brows.emplace_back(0); Brows.emplace_back(0);
unsigned int sumBlocks = 0; unsigned int sumBlocks = 0;
for (auto rowB = duneB_.begin(); rowB != duneB_.end(); ++rowB) { for (auto rowB = linSys_.duneB_.begin(); rowB != linSys_.duneB_.end(); ++rowB) {
int sizeRow = 0; int sizeRow = 0;
for (auto colB = rowB->begin(), endB = rowB->end(); colB != endB; ++colB) { for (auto colB = rowB->begin(), endB = rowB->end(); colB != endB; ++colB) {
Bcols.emplace_back(colB.index()); Bcols.emplace_back(colB.index());

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

@@ -22,17 +22,13 @@
#ifndef OPM_MULTISEGMENTWELL_EVAL_HEADER_INCLUDED #ifndef OPM_MULTISEGMENTWELL_EVAL_HEADER_INCLUDED
#define 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/simulators/wells/MultisegmentWellGeneric.hpp>
#include <opm/material/densead/Evaluation.hpp> #include <opm/material/densead/Evaluation.hpp>
#include <opm/input/eclipse/Schedule/Well/Well.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 <array>
#include <memory> #include <memory>
#include <utility> #include <utility>
@@ -91,24 +87,10 @@ protected:
// the number of well equations TODO: it should have a more general strategy for it // the number of well equations TODO: it should have a more general strategy for it
static constexpr int numWellEq = Indices::numPhases + 1; static constexpr int numWellEq = Indices::numPhases + 1;
// sparsity pattern for the matrices using Equations = MultisegmentWellEquations<Scalar,numWellEq,Indices::numEq>;
// [A C^T [x = [ res
// B D ] x_well] res_well]
// the vector type for the res_well and x_well using BVector = typename Equations::BVector;
using VectorBlockWellType = Dune::FieldVector<Scalar, numWellEq>; using BVectorWell = typename Equations::BVectorWell;
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>;
// TODO: for more efficient implementation, we should have EvalReservoir, EvalWell, and EvalRerservoirAndWell // TODO: for more efficient implementation, we should have EvalReservoir, EvalWell, and EvalRerservoirAndWell
// EvalR (Eval), EvalW, EvalRW // EvalR (Eval), EvalW, EvalRW
@@ -245,20 +227,7 @@ protected:
const WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif_; const WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif_;
// TODO, the following should go to a class for computing purpose Equations linSys_; //!< The equation system
// 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_;
// the values for the primary varibles // the values for the primary varibles
// based on different solutioin strategies, the wells can have different primary variables // based on different solutioin strategies, the wells can have different primary variables

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

@@ -203,15 +203,16 @@ namespace Opm
// Contributions are already in the matrix itself // Contributions are already in the matrix itself
return; return;
} }
BVectorWell Bx(this->duneB_.N()); BVectorWell Bx(this->linSys_.duneB_.N());
this->duneB_.mv(x, Bx); this->linSys_.duneB_.mv(x, Bx);
// invDBx = duneD^-1 * Bx_ // invDBx = duneD^-1 * Bx_
const BVectorWell invDBx = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, Bx); const BVectorWell invDBx = mswellhelpers::applyUMFPack(this->linSys_.duneD_,
this->linSys_.duneDSolver_, Bx);
// Ax = Ax - duneC_^T * invDBx // Ax = Ax - duneC_^T * invDBx
this->duneC_.mmtv(invDBx,Ax); this->linSys_.duneC_.mmtv(invDBx,Ax);
} }
@@ -226,9 +227,11 @@ namespace Opm
if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return; if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
// invDrw_ = duneD^-1 * resWell_ // invDrw_ = duneD^-1 * resWell_
const BVectorWell invDrw = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_); const BVectorWell invDrw = mswellhelpers::applyUMFPack(this->linSys_.duneD_,
this->linSys_.duneDSolver_,
this->linSys_.resWell_);
// r = r - duneC_^T * invDrw // r = r - duneC_^T * invDrw
this->duneC_.mmtv(invDrw, r); this->linSys_.duneC_.mmtv(invDrw, r);
} }
@@ -526,7 +529,9 @@ namespace Opm
// We assemble the well equations, then we check the convergence, // We assemble the well equations, then we check the convergence,
// which is why we do not put the assembleWellEq here. // 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 = mswellhelpers::applyUMFPack(this->linSys_.duneD_,
this->linSys_.duneDSolver_,
this->linSys_.resWell_);
updateWellState(dx_well, well_state, deferred_logger); updateWellState(dx_well, well_state, deferred_logger);
} }
@@ -727,7 +732,7 @@ namespace Opm
MultisegmentWell<TypeTag>:: MultisegmentWell<TypeTag>::
addWellContributions(SparseMatrixAdapter& jacobian) const addWellContributions(SparseMatrixAdapter& jacobian) const
{ {
const auto invDuneD = mswellhelpers::invertWithUMFPack<BVectorWell>(this->duneD_, this->duneDSolver_); const auto invDuneD = mswellhelpers::invertWithUMFPack<BVectorWell>(this->linSys_.duneD_, this->linSys_.duneDSolver_);
// We need to change matrix A as follows // We need to change matrix A as follows
// A -= C^T D^-1 B // A -= C^T D^-1 B
@@ -737,13 +742,15 @@ namespace Opm
// perforation at cell j connected to segment i. The code // perforation at cell j connected to segment i. The code
// assumes that no cell is connected to more than one segment, // assumes that no cell is connected to more than one segment,
// i.e. the columns of B/C have no more than one nonzero. // i.e. the columns of B/C have no more than one nonzero.
for (size_t rowC = 0; rowC < this->duneC_.N(); ++rowC) { for (size_t rowC = 0; rowC < this->linSys_.duneC_.N(); ++rowC) {
for (auto colC = this->duneC_[rowC].begin(), endC = this->duneC_[rowC].end(); colC != endC; ++colC) { for (auto colC = this->linSys_.duneC_[rowC].begin(),
endC = this->linSys_.duneC_[rowC].end(); colC != endC; ++colC) {
const auto row_index = colC.index(); const auto row_index = colC.index();
for (size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) { for (size_t rowB = 0; rowB < this->linSys_.duneB_.N(); ++rowB) {
for (auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) { for (auto colB = this->linSys_.duneB_[rowB].begin(),
endB = this->linSys_.duneB_[rowB].end(); colB != endB; ++colB) {
const auto col_index = colB.index(); const auto col_index = colB.index();
OffDiagMatrixBlockWellType tmp1; typename Equations::OffDiagMatrixBlockWellType tmp1;
detail::multMatrixImpl(invDuneD[rowC][rowB], (*colB), tmp1, std::true_type()); detail::multMatrixImpl(invDuneD[rowC][rowB], (*colB), tmp1, std::true_type());
typename SparseMatrixAdapter::MatrixBlock tmp2; typename SparseMatrixAdapter::MatrixBlock tmp2;
detail::multMatrixTransposedImpl((*colC), tmp1, tmp2, std::false_type()); detail::multMatrixTransposedImpl((*colC), tmp1, tmp2, std::false_type());
@@ -768,10 +775,11 @@ namespace Opm
// Add for coupling from well to reservoir // Add for coupling from well to reservoir
const auto seg_pressure_var_ind = this->SPres; const auto seg_pressure_var_ind = this->SPres;
const int welldof_ind = this->duneC_.M() + this->index_of_well_; const int welldof_ind = this->linSys_.duneC_.M() + this->index_of_well_;
if(not(this->isPressureControlled(well_state))){ if(not(this->isPressureControlled(well_state))){
for (size_t rowC = 0; rowC < this->duneC_.N(); ++rowC) { for (size_t rowC = 0; rowC < this->linSys_.duneC_.N(); ++rowC) {
for (auto colC = this->duneC_[rowC].begin(), endC = this->duneC_[rowC].end(); colC != endC; ++colC) { for (auto colC = this->linSys_.duneC_[rowC].begin(),
endC = this->linSys_.duneC_[rowC].end(); colC != endC; ++colC) {
const auto row_index = colC.index(); const auto row_index = colC.index();
const auto& bw = weights[row_index]; const auto& bw = weights[row_index];
double matel = 0.0; double matel = 0.0;
@@ -789,8 +797,9 @@ namespace Opm
auto well_weight = weights[0]; auto well_weight = weights[0];
well_weight = 0.0; well_weight = 0.0;
int num_perfs = 0; int num_perfs = 0;
for (size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) { for (size_t rowB = 0; rowB < this->linSys_.duneB_.N(); ++rowB) {
for (auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) { for (auto colB = this->linSys_.duneB_[rowB].begin(),
endB = this->linSys_.duneB_[rowB].end(); colB != endB; ++colB) {
const auto col_index = colB.index(); const auto col_index = colB.index();
const auto& bw = weights[col_index]; const auto& bw = weights[col_index];
well_weight += bw; well_weight += bw;
@@ -803,9 +812,10 @@ namespace Opm
// Add for coupling from reservoir to well and caclulate diag elelement corresping to incompressible standard well // Add for coupling from reservoir to well and caclulate diag elelement corresping to incompressible standard well
double diag_ell = 0.0; double diag_ell = 0.0;
for (size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) { for (size_t rowB = 0; rowB < this->linSys_.duneB_.N(); ++rowB) {
const auto& bw = well_weight; const auto& bw = well_weight;
for (auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) { for (auto colB = this->linSys_.duneB_[rowB].begin(),
endB = this->linSys_.duneB_[rowB].end(); colB != endB; ++colB) {
const auto col_index = colB.index(); const auto col_index = colB.index();
double matel = 0.0; double matel = 0.0;
for(size_t i = 0; i< bw.size(); ++i){ for(size_t i = 0; i< bw.size(); ++i){
@@ -1494,7 +1504,9 @@ namespace Opm
assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); 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 = mswellhelpers::applyUMFPack(this->linSys_.duneD_,
this->linSys_.duneDSolver_,
this->linSys_.resWell_);
if (it > this->param_.strict_inner_iter_wells_) { if (it > this->param_.strict_inner_iter_wells_) {
relax_convergence = true; relax_convergence = true;
@@ -1622,13 +1634,13 @@ namespace Opm
computeSegmentFluidProperties(ebosSimulator, deferred_logger); computeSegmentFluidProperties(ebosSimulator, deferred_logger);
// clear all entries // clear all entries
this->duneB_ = 0.0; this->linSys_.duneB_ = 0.0;
this->duneC_ = 0.0; this->linSys_.duneC_ = 0.0;
this->duneD_ = 0.0; this->linSys_.duneD_ = 0.0;
this->resWell_ = 0.0; this->linSys_.resWell_ = 0.0;
this->duneDSolver_.reset(); this->linSys_.duneDSolver_.reset();
auto& ws = well_state.well(this->index_of_well_); auto& ws = well_state.well(this->index_of_well_);
ws.dissolved_gas_rate = 0; ws.dissolved_gas_rate = 0;
@@ -1659,9 +1671,9 @@ namespace Opm
const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->surfaceVolumeFraction(seg, comp_idx) const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->surfaceVolumeFraction(seg, comp_idx)
- segment_fluid_initial_[seg][comp_idx]) / dt; - 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) { 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 +1685,13 @@ namespace Opm
const int seg_upwind = this->upwinding_segments_[seg]; const int seg_upwind = this->upwinding_segments_[seg];
// segment_rate contains the derivatives with respect to WQTotal in seg, // segment_rate contains the derivatives with respect to WQTotal in seg,
// and WFrac and GFrac in seg_upwind // and WFrac and GFrac in seg_upwind
this->resWell_[seg][comp_idx] -= segment_rate.value(); this->linSys_.resWell_[seg][comp_idx] -= segment_rate.value();
this->duneD_[seg][seg][comp_idx][WQTotal] -= segment_rate.derivative(WQTotal + Indices::numEq); this->linSys_.duneD_[seg][seg][comp_idx][WQTotal] -= segment_rate.derivative(WQTotal + Indices::numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { 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)) { 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 // pressure derivative should be zero
} }
@@ -1694,13 +1706,13 @@ namespace Opm
const int inlet_upwind = this->upwinding_segments_[inlet]; const int inlet_upwind = this->upwinding_segments_[inlet];
// inlet_rate contains the derivatives with respect to WQTotal in inlet, // inlet_rate contains the derivatives with respect to WQTotal in inlet,
// and WFrac and GFrac in inlet_upwind // and WFrac and GFrac in inlet_upwind
this->resWell_[seg][comp_idx] += inlet_rate.value(); this->linSys_.resWell_[seg][comp_idx] += inlet_rate.value();
this->duneD_[seg][inlet][comp_idx][WQTotal] += inlet_rate.derivative(WQTotal + Indices::numEq); this->linSys_.duneD_[seg][inlet][comp_idx][WQTotal] += inlet_rate.derivative(WQTotal + Indices::numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { 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)) { 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 // pressure derivative should be zero
} }
@@ -1744,21 +1756,21 @@ namespace Opm
this->connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective); this->connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective);
// subtract sum of phase fluxes in the well equations. // 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 // assemble the jacobians
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
// also need to consider the efficiency factor when manipulating the jacobians. // 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 // 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) { for (int pv_idx = 0; pv_idx < Indices::numEq; ++pv_idx) {
// also need to consider the efficiency factor when manipulating the jacobians. // 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);
} }
} }
} }