diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index b2f0a51e5..0ccd767c1 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -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 diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index e5c8f398b..e7138ff87 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -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; diff --git a/opm/simulators/wells/MultisegmentWellEquations.cpp b/opm/simulators/wells/MultisegmentWellEquations.cpp new file mode 100644 index 000000000..954644706 --- /dev/null +++ b/opm/simulators/wells/MultisegmentWellEquations.cpp @@ -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 . +*/ + +#include +#include + + +namespace Opm +{ + +} diff --git a/opm/simulators/wells/MultisegmentWellEquations.hpp b/opm/simulators/wells/MultisegmentWellEquations.hpp new file mode 100644 index 000000000..cfa62aee4 --- /dev/null +++ b/opm/simulators/wells/MultisegmentWellEquations.hpp @@ -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 . +*/ + +#ifndef OPM_MULTISEGMENTWELL_EQUATIONS_HEADER_INCLUDED +#define OPM_MULTISEGMENTWELL_EQUATIONS_HEADER_INCLUDED + +#include +#include +#include +#include + +#include + +namespace Dune { +template class UMFPack; +} + +namespace Opm +{ + +template +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; + using BVectorWell = Dune::BlockVector; + + using VectorBlockType = Dune::FieldVector; + using BVector = Dune::BlockVector; + + // the matrix type for the diagonal matrix D + using DiagMatrixBlockWellType = Dune::FieldMatrix; + using DiagMatWell = Dune::BCRSMatrix; + + // the matrix type for the non-diagonal matrix B and C^T + using OffDiagMatrixBlockWellType = Dune::FieldMatrix; + using OffDiagMatWell = Dune::BCRSMatrix; + + // 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 > duneDSolver_; + + // residuals of the well equations + BVectorWell resWell_; +}; + +} + +#endif // OPM_MULTISEGMENTWELLWELL_EQUATIONS_HEADER_INCLUDED diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index e1c3f53e8..8d473f718 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -74,9 +74,9 @@ void MultisegmentWellEval:: initMatrixAndVectors(const int num_cells) { - duneB_.setBuildMode(OffDiagMatWell::row_wise); - duneC_.setBuildMode(OffDiagMatWell::row_wise); - duneD_.setBuildMode(DiagMatWell::row_wise); + linSys_.duneB_.setBuildMode(Equations::OffDiagMatWell::row_wise); + linSys_.duneC_.setBuildMode(Equations::OffDiagMatWell::row_wise); + linSys_.duneD_.setBuildMode(Equations::DiagMatWell::row_wise); // set the size and patterns for all the matrices and vectors // [A C^T [x = [ res @@ -89,13 +89,14 @@ initMatrixAndVectors(const int num_cells) for (const std::vector& inlets : this->segment_inlets_) { 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()); - duneC_.setSize(this->numberOfSegments(), num_cells, baseif_.numPerfs()); + linSys_.duneB_.setSize(this->numberOfSegments(), num_cells, baseif_.numPerfs()); + linSys_.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) { + for (auto row = linSys_.duneD_.createbegin(), + end = linSys_.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 @@ -116,7 +117,8 @@ initMatrixAndVectors(const int num_cells) } // 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. for (const int& perf : this->segment_perforations_[row.index()]) { const int cell_idx = baseif_.cells()[perf]; @@ -125,7 +127,8 @@ initMatrixAndVectors(const int num_cells) } // 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. for (const int& perf : this->segment_perforations_[row.index()]) { 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_evaluation_.resize(this->numberOfSegments()); @@ -172,7 +175,7 @@ getWellConvergence(const WellState& well_state, std::vector> abs_residual(this->numberOfSegments(), std::vector(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 +240,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); @@ -455,11 +458,11 @@ recoverSolutionWell(const BVector& x, BVectorWell& xw) const { if (!baseif_.isOperableAndSolvable() && !baseif_.wellIsStopped()) return; - BVectorWell resWell = resWell_; + BVectorWell resWell = linSys_.resWell_; // resWell = resWell - B * x - duneB_.mmv(x, resWell); + linSys_.duneB_.mmv(x, resWell); // xw = D^-1 * resWell - xw = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell); + xw = mswellhelpers::applyUMFPack(linSys_.duneD_, linSys_.duneDSolver_, resWell); } template @@ -1289,9 +1292,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 +1339,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 +1377,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 +1612,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 +1647,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 +1700,10 @@ getFiniteWellResiduals(const std::vector& 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 +1720,7 @@ getFiniteWellResiduals(const std::vector& 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}; @@ -1863,16 +1866,16 @@ void MultisegmentWellEval:: 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(); + unsigned int Mb = linSys_.duneB_.N(); // number of blockrows in duneB_, duneC_ and duneD_ + unsigned int BnumBlocks = linSys_.duneB_.nonzeroes(); + unsigned int DnumBlocks = linSys_.duneD_.nonzeroes(); // duneC std::vector Ccols; std::vector Cvals; Ccols.reserve(BnumBlocks); 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) { Ccols.emplace_back(colC.index()); for (int i = 0; i < numWellEq; ++i) { @@ -1884,7 +1887,7 @@ addWellContribution(WellContributions& wellContribs) const } // duneD - Dune::UMFPack umfpackMatrix(duneD_, 0); + Dune::UMFPack umfpackMatrix(linSys_.duneD_, 0); double *Dvals = umfpackMatrix.getInternalMatrix().getValues(); auto *Dcols = umfpackMatrix.getInternalMatrix().getColStart(); auto *Drows = umfpackMatrix.getInternalMatrix().getRowIndex(); @@ -1898,7 +1901,7 @@ addWellContribution(WellContributions& wellContribs) const Bvals.reserve(BnumBlocks * Indices::numEq * numWellEq); Brows.emplace_back(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; for (auto colB = rowB->begin(), endB = rowB->end(); colB != endB; ++colB) { Bcols.emplace_back(colB.index()); diff --git a/opm/simulators/wells/MultisegmentWellEval.hpp b/opm/simulators/wells/MultisegmentWellEval.hpp index d53c7d478..1b8ae3250 100644 --- a/opm/simulators/wells/MultisegmentWellEval.hpp +++ b/opm/simulators/wells/MultisegmentWellEval.hpp @@ -22,17 +22,13 @@ #ifndef OPM_MULTISEGMENTWELL_EVAL_HEADER_INCLUDED #define OPM_MULTISEGMENTWELL_EVAL_HEADER_INCLUDED +#include #include #include #include -#include -#include -#include -#include - #include #include #include @@ -91,24 +87,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; - // the vector type for the res_well and x_well - using VectorBlockWellType = Dune::FieldVector; - using BVectorWell = Dune::BlockVector; - - using VectorBlockType = Dune::FieldVector; - using BVector = Dune::BlockVector; - - // the matrix type for the diagonal matrix D - using DiagMatrixBlockWellType = Dune::FieldMatrix; - using DiagMatWell = Dune::BCRSMatrix; - - // the matrix type for the non-diagonal matrix B and C^T - using OffDiagMatrixBlockWellType = Dune::FieldMatrix; - using OffDiagMatWell = Dune::BCRSMatrix; + 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 @@ -245,20 +227,7 @@ protected: const WellInterfaceIndices& 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 > 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 diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 194e81fd1..87dd6cacb 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -203,15 +203,16 @@ namespace Opm // Contributions are already in the matrix itself 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_ - 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 - this->duneC_.mmtv(invDBx,Ax); + this->linSys_.duneC_.mmtv(invDBx,Ax); } @@ -226,9 +227,11 @@ namespace Opm if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return; // 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 - 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, // 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); } @@ -727,7 +732,7 @@ namespace Opm MultisegmentWell:: addWellContributions(SparseMatrixAdapter& jacobian) const { - const auto invDuneD = mswellhelpers::invertWithUMFPack(this->duneD_, this->duneDSolver_); + const auto invDuneD = mswellhelpers::invertWithUMFPack(this->linSys_.duneD_, this->linSys_.duneDSolver_); // We need to change matrix A as follows // A -= C^T D^-1 B @@ -737,13 +742,15 @@ namespace Opm // 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) { + for (size_t rowC = 0; rowC < this->linSys_.duneC_.N(); ++rowC) { + for (auto colC = this->linSys_.duneC_[rowC].begin(), + endC = this->linSys_.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) { + for (size_t rowB = 0; rowB < this->linSys_.duneB_.N(); ++rowB) { + for (auto colB = this->linSys_.duneB_[rowB].begin(), + endB = this->linSys_.duneB_[rowB].end(); colB != endB; ++colB) { const auto col_index = colB.index(); - OffDiagMatrixBlockWellType tmp1; + typename Equations::OffDiagMatrixBlockWellType tmp1; detail::multMatrixImpl(invDuneD[rowC][rowB], (*colB), tmp1, std::true_type()); typename SparseMatrixAdapter::MatrixBlock tmp2; detail::multMatrixTransposedImpl((*colC), tmp1, tmp2, std::false_type()); @@ -768,10 +775,11 @@ namespace Opm // 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_; + const int welldof_ind = this->linSys_.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) { + for (size_t rowC = 0; rowC < this->linSys_.duneC_.N(); ++rowC) { + 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& bw = weights[row_index]; double matel = 0.0; @@ -789,8 +797,9 @@ namespace Opm 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) { + for (size_t rowB = 0; rowB < this->linSys_.duneB_.N(); ++rowB) { + 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& bw = weights[col_index]; 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 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; - 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(); double matel = 0.0; 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); - 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_) { relax_convergence = true; @@ -1622,13 +1634,13 @@ namespace Opm computeSegmentFluidProperties(ebosSimulator, deferred_logger); // clear all entries - this->duneB_ = 0.0; - this->duneC_ = 0.0; + this->linSys_.duneB_ = 0.0; + this->linSys_.duneC_ = 0.0; - this->duneD_ = 0.0; - this->resWell_ = 0.0; + this->linSys_.duneD_ = 0.0; + this->linSys_.resWell_ = 0.0; - this->duneDSolver_.reset(); + this->linSys_.duneDSolver_.reset(); auto& ws = well_state.well(this->index_of_well_); 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) - 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 +1685,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 +1706,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 +1756,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); } } }