From bc45c1e80a068a5bf20bcc2b6b59bca722910d9f Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 20 Dec 2024 12:31:29 +0100 Subject: [PATCH 1/3] changed: remove unused code in WellConnectionAuxilliaryModule --- CMakeLists_files.cmake | 1 - .../wells/WellConnectionAuxiliaryModule.cpp | 83 ------------------- .../wells/WellConnectionAuxiliaryModule.hpp | 17 ---- 3 files changed, 101 deletions(-) delete mode 100644 opm/simulators/wells/WellConnectionAuxiliaryModule.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index ba25c05fb..e91d5729a 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -201,7 +201,6 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/wells/VFPProdProperties.cpp opm/simulators/wells/WellAssemble.cpp opm/simulators/wells/WellBhpThpCalculator.cpp - opm/simulators/wells/WellConnectionAuxiliaryModule.cpp opm/simulators/wells/WellConstraints.cpp opm/simulators/wells/WellConvergence.cpp opm/simulators/wells/WellFilterCake.cpp diff --git a/opm/simulators/wells/WellConnectionAuxiliaryModule.cpp b/opm/simulators/wells/WellConnectionAuxiliaryModule.cpp deleted file mode 100644 index dd65295ea..000000000 --- a/opm/simulators/wells/WellConnectionAuxiliaryModule.cpp +++ /dev/null @@ -1,83 +0,0 @@ -/* - Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services - Copyright 2017 Statoil ASA. - - 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 - -#include - -#include -#include -#include - -#include - -namespace Opm -{ - -WellConnectionAuxiliaryModuleGeneric:: -WellConnectionAuxiliaryModuleGeneric(const Schedule& schedule, - const Dune::CpGrid& grid) -{ - // Create cartesian to compressed mapping - const auto& globalCell = grid.globalCell(); - const auto& cartesianSize = grid.logicalCartesianSize(); - - auto size = cartesianSize[0] * cartesianSize[1] * cartesianSize[2]; - - std::vector cartesianToCompressed(size, -1); - auto begin = globalCell.begin(); - - for (auto cell = begin, end = globalCell.end(); cell != end; ++cell) - { - cartesianToCompressed[ *cell ] = cell - begin; - } - - const auto& schedule_wells = schedule.getWellsatEnd(); - wells_.reserve(schedule_wells.size()); - - // initialize the additional cell connections introduced by wells. - for (const auto& well : schedule_wells) - { - std::vector compressed_well_perforations; - // All possible completions of the well - const auto& completionSet = well.getConnections(); - compressed_well_perforations.reserve(completionSet.size()); - - for (const auto& completion : completionSet) - { - int compressed_idx = cartesianToCompressed[completion.global_index()]; - if (compressed_idx >= 0) // Ignore completions in inactive/remote cells. - { - compressed_well_perforations.push_back(compressed_idx); - } - } - - if (!compressed_well_perforations.empty()) - { - std::sort(compressed_well_perforations.begin(), - compressed_well_perforations.end()); - - wells_.push_back(compressed_well_perforations); - } - } -} - -} // end namespace Opm diff --git a/opm/simulators/wells/WellConnectionAuxiliaryModule.hpp b/opm/simulators/wells/WellConnectionAuxiliaryModule.hpp index f07b252eb..89fea0072 100644 --- a/opm/simulators/wells/WellConnectionAuxiliaryModule.hpp +++ b/opm/simulators/wells/WellConnectionAuxiliaryModule.hpp @@ -32,19 +32,9 @@ namespace Opm class Schedule; -class WellConnectionAuxiliaryModuleGeneric -{ -protected: - WellConnectionAuxiliaryModuleGeneric(const Schedule& schedule, - const Dune::CpGrid& grid); - - std::vector > wells_; -}; - template class WellConnectionAuxiliaryModule : public BaseAuxiliaryModule - , private WellConnectionAuxiliaryModuleGeneric { using GlobalEqVector = GetPropType; using SparseMatrixAdapter = GetPropType; @@ -56,7 +46,6 @@ public: WellConnectionAuxiliaryModule(const Schedule& schedule, const Dune::CpGrid& grid) - : WellConnectionAuxiliaryModuleGeneric(schedule, grid) { } @@ -68,12 +57,6 @@ public: void addNeighbors(std::vector& neighbors) const { - for (const auto& well_perforations : wells_) - { - for (const auto& perforation : well_perforations) - neighbors[perforation].insert(well_perforations.begin(), - well_perforations.end()); - } } void applyInitial() From 2800f35b73417bcb7a6f760ccc671272f255044b Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 20 Dec 2024 12:40:00 +0100 Subject: [PATCH 2/3] changed: move the auxillary well module from BlackoilModule to WellConnectionAuxillary this is an adapter for the non-linear solver. it is still model using is-a but this can now be changed if desired. --- opm/simulators/wells/BlackoilWellModel.hpp | 115 ++++--------- .../wells/BlackoilWellModelGeneric.hpp | 3 +- .../wells/BlackoilWellModel_impl.hpp | 134 ++------------- .../wells/WellConnectionAuxiliaryModule.hpp | 159 +++++++++++++++--- 4 files changed, 192 insertions(+), 219 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index c0c3e5ade..83d0cad25 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -24,13 +24,11 @@ #ifndef OPM_BLACKOILWELLMODEL_HEADER_INCLUDED #define OPM_BLACKOILWELLMODEL_HEADER_INCLUDED -#include +#include +#include +#include -#include -#include -#include -#include -#include +#include #include @@ -39,7 +37,7 @@ #include #include -#include +#include #include #include @@ -47,6 +45,11 @@ #include +#include +#include + +#include + #include #include #include @@ -63,22 +66,18 @@ #include #include #include -#include +#include #include #include #include #include +#include -#include -#include - -#include -#include -#include - -#include - -#include +#include +#include +#include +#include +#include namespace Opm { @@ -88,7 +87,7 @@ template class WellContributions; /// Class for handling the blackoil well model. template - class BlackoilWellModel : public BaseAuxiliaryModule + class BlackoilWellModel : public WellConnectionAuxiliaryModule> , public BlackoilWellModelGeneric> { @@ -106,8 +105,9 @@ template class WellContributions; using SparseMatrixAdapter = GetPropType; using ModelParameters = BlackoilModelParameters; + using WellConnectionModule = WellConnectionAuxiliaryModule>; + constexpr static std::size_t pressureVarIndex = GetPropType::pressureSwitchIdx; - typedef typename BaseAuxiliaryModule::NeighborSet NeighborSet; static const int numEq = Indices::numEq; static const int solventSaturationIdx = Indices::solventSaturationIdx; @@ -139,51 +139,6 @@ template class WellContributions; void init(); void initWellContainer(const int reportStepIdx) override; - ///////////// - // - ///////////// - unsigned numDofs() const override - // No extra dofs are inserted for wells. (we use a Schur complement.) - { return 0; } - - void addNeighbors(std::vector& neighbors) const override; - - void applyInitial() override - {} - - void linearize(SparseMatrixAdapter& jacobian, GlobalEqVector& res) override; - void linearizeDomain(const Domain& domain, SparseMatrixAdapter& jacobian, GlobalEqVector& res); - - void postSolve(GlobalEqVector& deltaX) override - { - recoverWellSolutionAndUpdateWellState(deltaX); - } - - void postSolveDomain(GlobalEqVector& deltaX, const Domain& domain) - { - recoverWellSolutionAndUpdateWellStateDomain(deltaX, domain); - } - - ///////////// - // - ///////////// - - template - void deserialize(Restarter& /* res */) - { - // TODO (?) - } - - /*! - * \brief This method writes the complete state of the well - * to the harddisk. - */ - template - void serialize(Restarter& /* res*/) - { - // TODO (?) - } - void beginEpisode() { OPM_TIMEBLOCK(beginEpsiode); @@ -371,6 +326,22 @@ template class WellContributions; auto end() const { return well_container_.end(); } bool empty() const { return well_container_.empty(); } + bool addMatrixContributions() const { return param_.matrix_add_well_contributions_; } + + int compressedIndexForInterior(int cartesian_cell_idx) const override + { + return simulator_.vanguard().compressedIndexForInterior(cartesian_cell_idx); + } + + // using the solution x to recover the solution xw for wells and applying + // xw to update Well State + void recoverWellSolutionAndUpdateWellState(const BVector& x); + + // using the solution x to recover the solution xw for wells and applying + // xw to update Well State + void recoverWellSolutionAndUpdateWellStateDomain(const BVector& x, + const Domain& domain); + protected: Simulator& simulator_; @@ -470,14 +441,6 @@ template class WellContributions; // called at the end of a report step void endReportStep(); - // using the solution x to recover the solution xw for wells and applying - // xw to update Well State - void recoverWellSolutionAndUpdateWellState(const BVector& x); - - // using the solution x to recover the solution xw for wells and applying - // xw to update Well State - void recoverWellSolutionAndUpdateWellStateDomain(const BVector& x, const Domain& domain); - // setting the well_solutions_ based on well_state. void updatePrimaryVariables(DeferredLogger& deferred_logger); @@ -529,10 +492,6 @@ template class WellContributions; void computeWellTemperature(); - int compressedIndexForInterior(int cartesian_cell_idx) const override { - return simulator_.vanguard().compressedIndexForInterior(cartesian_cell_idx); - } - private: BlackoilWellModel(Simulator& simulator, const PhaseUsage& pu); @@ -543,8 +502,6 @@ template class WellContributions; // Their state is not relevant between function calls, so they can // (and must) be mutable, as the functions using them are const. mutable BVector x_local_; - mutable BVector res_local_; - mutable GlobalEqVector linearize_res_local_; }; diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.hpp b/opm/simulators/wells/BlackoilWellModelGeneric.hpp index 80f083c25..6d2ff6f09 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.hpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.hpp @@ -214,6 +214,8 @@ public: return well_domain_; } + std::vector getCellsForConnections(const Well& well) const; + bool reportStepStarts() const { return report_step_starts_; } bool shouldBalanceNetwork(const int reportStepIndex, @@ -435,7 +437,6 @@ protected: /// \brief get compressed index for interior cells (-1, otherwise virtual int compressedIndexForInterior(int cartesian_cell_idx) const = 0; - std::vector getCellsForConnections(const Well& well) const; std::vector> getMaxWellConnections() const; std::vector getWellsForTesting(const int timeStepIdx, diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 34606f29a..c45e9658a 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -57,10 +57,6 @@ #include #endif -#if HAVE_MPI -#include -#endif - #include #include #include @@ -73,7 +69,8 @@ namespace Opm { template BlackoilWellModel:: BlackoilWellModel(Simulator& simulator, const PhaseUsage& phase_usage) - : BlackoilWellModelGeneric(simulator.vanguard().schedule(), + : WellConnectionModule(*this, simulator.gridView().comm()) + , BlackoilWellModelGeneric(simulator.vanguard().schedule(), simulator.vanguard().summaryState(), simulator.vanguard().eclState(), phase_usage, @@ -190,116 +187,6 @@ namespace Opm { } } - - template - void - BlackoilWellModel:: - addNeighbors(std::vector& neighbors) const - { - if (!param_.matrix_add_well_contributions_) { - return; - } - - // Create cartesian to compressed mapping - const auto& schedule_wells = this->schedule().getWellsatEnd(); - auto possibleFutureConnections = this->schedule().getPossibleFutureConnections(); - -#if HAVE_MPI - // Communicate Map to other processes, since it is only available on rank 0 - const auto& comm = this->simulator_.vanguard().grid().comm(); - Parallel::MpiSerializer ser(comm); - ser.broadcast(possibleFutureConnections); -#endif - // initialize the additional cell connections introduced by wells. - for (const auto& well : schedule_wells) - { - std::vector wellCells = this->getCellsForConnections(well); - // Now add the cells of the possible future connections - const auto possibleFutureConnectionSetIt = possibleFutureConnections.find(well.name()); - if (possibleFutureConnectionSetIt != possibleFutureConnections.end()) { - for (auto& global_index : possibleFutureConnectionSetIt->second) { - int compressed_idx = compressedIndexForInterior(global_index); - if (compressed_idx >= 0) { // Ignore connections in inactive/remote cells. - wellCells.push_back(compressed_idx); - } - } - } - for (int cellIdx : wellCells) { - neighbors[cellIdx].insert(wellCells.begin(), - wellCells.end()); - } - } - } - - - template - void - BlackoilWellModel:: - linearize(SparseMatrixAdapter& jacobian, GlobalEqVector& res) - { - OPM_BEGIN_PARALLEL_TRY_CATCH(); - for (const auto& well: well_container_) { - // Modifiy the Jacobian with explicit Schur complement - // contributions if requested. - if (param_.matrix_add_well_contributions_) { - well->addWellContributions(jacobian); - } - // Apply as Schur complement the well residual to reservoir residuals: - // r = r - duneC_^T * invDuneD_ * resWell_ - // Well equations B and C uses only the perforated cells, so need to apply on local residual - const auto& cells = well->cells(); - linearize_res_local_.resize(cells.size()); - - for (size_t i = 0; i < cells.size(); ++i) { - linearize_res_local_[i] = res[cells[i]]; - } - - well->apply(linearize_res_local_); - - for (size_t i = 0; i < cells.size(); ++i) { - res[cells[i]] = linearize_res_local_[i]; - } - } - OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::linearize failed: ", - simulator_.gridView().comm()); - } - - - template - void - BlackoilWellModel:: - linearizeDomain(const Domain& domain, SparseMatrixAdapter& jacobian, GlobalEqVector& res) - { - // Note: no point in trying to do a parallel gathering - // try/catch here, as this function is not called in - // parallel but for each individual domain of each rank. - for (const auto& well: well_container_) { - if (this->well_domain_.at(well->name()) == domain.index) { - // Modifiy the Jacobian with explicit Schur complement - // contributions if requested. - if (param_.matrix_add_well_contributions_) { - well->addWellContributions(jacobian); - } - // Apply as Schur complement the well residual to reservoir residuals: - // r = r - duneC_^T * invDuneD_ * resWell_ - // Well equations B and C uses only the perforated cells, so need to apply on local residual - const auto& cells = well->cells(); - linearize_res_local_.resize(cells.size()); - - for (size_t i = 0; i < cells.size(); ++i) { - linearize_res_local_[i] = res[cells[i]]; - } - - well->apply(linearize_res_local_); - - for (size_t i = 0; i < cells.size(); ++i) { - res[cells[i]] = linearize_res_local_[i]; - } - } - } - } - - template void BlackoilWellModel:: @@ -1711,7 +1598,9 @@ namespace Opm { template void BlackoilWellModel:: - addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights, const bool use_well_weights) const + addWellPressureEquations(PressureMatrix& jacobian, + const BVector& weights, + const bool use_well_weights) const { int nw = this->numLocalWellsEnd(); int rdofs = local_num_cells_; @@ -1720,8 +1609,12 @@ namespace Opm { jacobian[wdof][wdof] = 1.0;// better scaling ? } - for ( const auto& well : well_container_ ) { - well->addWellPressureEquations(jacobian, weights, pressureVarIndex, use_well_weights, this->wellState()); + for (const auto& well : well_container_) { + well->addWellPressureEquations(jacobian, + weights, + pressureVarIndex, + use_well_weights, + this->wellState()); } } @@ -1810,14 +1703,15 @@ namespace Opm { DeferredLogger local_deferredLogger; OPM_BEGIN_PARALLEL_TRY_CATCH(); { - for (auto& well : well_container_) { + for (const auto& well : well_container_) { const auto& cells = well->cells(); x_local_.resize(cells.size()); for (size_t i = 0; i < cells.size(); ++i) { x_local_[i] = x[cells[i]]; } - well->recoverWellSolutionAndUpdateWellState(simulator_, x_local_, this->wellState(), local_deferredLogger); + well->recoverWellSolutionAndUpdateWellState(simulator_, x_local_, + this->wellState(), local_deferredLogger); } } OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, diff --git a/opm/simulators/wells/WellConnectionAuxiliaryModule.hpp b/opm/simulators/wells/WellConnectionAuxiliaryModule.hpp index 89fea0072..e3e0ad3c0 100644 --- a/opm/simulators/wells/WellConnectionAuxiliaryModule.hpp +++ b/opm/simulators/wells/WellConnectionAuxiliaryModule.hpp @@ -23,51 +23,172 @@ #include -#include +#include -namespace Dune { class CpGrid; } +#include +#include -namespace Opm -{ - -class Schedule; - -template -class WellConnectionAuxiliaryModule - : public BaseAuxiliaryModule +#if HAVE_MPI +#include +#endif + +namespace Opm { + +template +class WellConnectionAuxiliaryModule : public BaseAuxiliaryModule { + using Grid = GetPropType; using GlobalEqVector = GetPropType; using SparseMatrixAdapter = GetPropType; public: - using NeighborSet = typename ::Opm::BaseAuxiliaryModule::NeighborSet; - WellConnectionAuxiliaryModule(const Schedule& schedule, - const Dune::CpGrid& grid) + using Domain = SubDomain; + + WellConnectionAuxiliaryModule(Model& model, Parallel::Communication comm) + : model_(model) + , lin_comm_(std::move(comm)) { } - unsigned numDofs() const + unsigned numDofs() const override { // No extra dofs are inserted for wells. return 0; } - void addNeighbors(std::vector& neighbors) const + void addNeighbors(std::vector& neighbors) const override { + if (!model_.addMatrixContributions()) { + return; + } + + // Create cartesian to compressed mapping + const auto& schedule_wells = model_.schedule().getWellsatEnd(); + auto possibleFutureConnections = model_.schedule().getPossibleFutureConnections(); + +#if HAVE_MPI + // Communicate Map to other processes, since it is only available on rank 0 + Parallel::MpiSerializer ser(lin_comm_); + ser.broadcast(possibleFutureConnections); +#endif + // initialize the additional cell connections introduced by wells. + for (const auto& well : schedule_wells) + { + std::vector wellCells = model_.getCellsForConnections(well); + // Now add the cells of the possible future connections + const auto possibleFutureConnectionSetIt = possibleFutureConnections.find(well.name()); + if (possibleFutureConnectionSetIt != possibleFutureConnections.end()) { + for (auto& global_index : possibleFutureConnectionSetIt->second) { + int compressed_idx = model_.compressedIndexForInterior(global_index); + if (compressed_idx >= 0) { // Ignore connections in inactive/remote cells. + wellCells.push_back(compressed_idx); + } + } + } + for (int cellIdx : wellCells) { + neighbors[cellIdx].insert(wellCells.begin(), + wellCells.end()); + } + } } - void applyInitial() + void applyInitial() override {} - void linearize(SparseMatrixAdapter& , GlobalEqVector&) + void linearize(SparseMatrixAdapter& jacobian, GlobalEqVector& res) override { - // Linearization is done in StandardDenseWells + OPM_BEGIN_PARALLEL_TRY_CATCH(); + for (const auto& well : model_) { + // Modifiy the Jacobian with explicit Schur complement + // contributions if requested. + if (model_.addMatrixContributions()) { + well->addWellContributions(jacobian); + } + // Apply as Schur complement the well residual to reservoir residuals: + // r = r - duneC_^T * invDuneD_ * resWell_ + // Well equations B and C uses only the perforated cells, so need to apply on local residual + const auto& cells = well->cells(); + linearize_res_local_.resize(cells.size()); + + for (size_t i = 0; i < cells.size(); ++i) { + linearize_res_local_[i] = res[cells[i]]; + } + + well->apply(linearize_res_local_); + + for (size_t i = 0; i < cells.size(); ++i) { + res[cells[i]] = linearize_res_local_[i]; + } + } + OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::linearize failed: ", lin_comm_); } - private: + void postSolve(GlobalEqVector& deltaX) override + { + model_.recoverWellSolutionAndUpdateWellState(deltaX); + } + + void linearizeDomain(const Domain& domain, + SparseMatrixAdapter& jacobian, + GlobalEqVector& res) + { + // Note: no point in trying to do a parallel gathering + // try/catch here, as this function is not called in + // parallel but for each individual domain of each rank. + for (const auto& well : model_) { + if (model_.well_domain().at(well->name()) == domain.index) { + // Modifiy the Jacobian with explicit Schur complement + // contributions if requested. + if (model_.addMatrixContributions()) { + well->addWellContributions(jacobian); + } + // Apply as Schur complement the well residual to reservoir residuals: + // r = r - duneC_^T * invDuneD_ * resWell_ + // Well equations B and C uses only the perforated cells, so need to apply on local residual + const auto& cells = well->cells(); + linearize_res_local_.resize(cells.size()); + + for (size_t i = 0; i < cells.size(); ++i) { + linearize_res_local_[i] = res[cells[i]]; + } + + well->apply(linearize_res_local_); + + for (size_t i = 0; i < cells.size(); ++i) { + res[cells[i]] = linearize_res_local_[i]; + } + } + } + } + + void postSolveDomain(GlobalEqVector& deltaX, const Domain& domain) + { + model_.recoverWellSolutionAndUpdateWellStateDomain(deltaX, domain); + } + + template + void deserialize(Restarter& /* res */) + { + // TODO (?) + } + + /*! + * \brief This method writes the complete state of the well + * to the harddisk. + */ + template + void serialize(Restarter& /* res*/) + { + // TODO (?) + } + +private: + Model& model_; + GlobalEqVector linearize_res_local_{}; + Parallel::Communication lin_comm_; }; } // end namespace OPM From 9e87c50b0b7bdd7a1467af02587da545e11f3001 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 20 Dec 2024 12:46:06 +0100 Subject: [PATCH 3/3] split out linearization for a single well use this to unify code between regular and domain --- .../wells/WellConnectionAuxiliaryModule.hpp | 65 +++++++------------ 1 file changed, 25 insertions(+), 40 deletions(-) diff --git a/opm/simulators/wells/WellConnectionAuxiliaryModule.hpp b/opm/simulators/wells/WellConnectionAuxiliaryModule.hpp index e3e0ad3c0..72ff79ed6 100644 --- a/opm/simulators/wells/WellConnectionAuxiliaryModule.hpp +++ b/opm/simulators/wells/WellConnectionAuxiliaryModule.hpp @@ -102,26 +102,7 @@ public: { OPM_BEGIN_PARALLEL_TRY_CATCH(); for (const auto& well : model_) { - // Modifiy the Jacobian with explicit Schur complement - // contributions if requested. - if (model_.addMatrixContributions()) { - well->addWellContributions(jacobian); - } - // Apply as Schur complement the well residual to reservoir residuals: - // r = r - duneC_^T * invDuneD_ * resWell_ - // Well equations B and C uses only the perforated cells, so need to apply on local residual - const auto& cells = well->cells(); - linearize_res_local_.resize(cells.size()); - - for (size_t i = 0; i < cells.size(); ++i) { - linearize_res_local_[i] = res[cells[i]]; - } - - well->apply(linearize_res_local_); - - for (size_t i = 0; i < cells.size(); ++i) { - res[cells[i]] = linearize_res_local_[i]; - } + this->linearizeSingleWell(jacobian, res, well); } OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::linearize failed: ", lin_comm_); } @@ -140,26 +121,7 @@ public: // parallel but for each individual domain of each rank. for (const auto& well : model_) { if (model_.well_domain().at(well->name()) == domain.index) { - // Modifiy the Jacobian with explicit Schur complement - // contributions if requested. - if (model_.addMatrixContributions()) { - well->addWellContributions(jacobian); - } - // Apply as Schur complement the well residual to reservoir residuals: - // r = r - duneC_^T * invDuneD_ * resWell_ - // Well equations B and C uses only the perforated cells, so need to apply on local residual - const auto& cells = well->cells(); - linearize_res_local_.resize(cells.size()); - - for (size_t i = 0; i < cells.size(); ++i) { - linearize_res_local_[i] = res[cells[i]]; - } - - well->apply(linearize_res_local_); - - for (size_t i = 0; i < cells.size(); ++i) { - res[cells[i]] = linearize_res_local_[i]; - } + this->linearizeSingleWell(jacobian, res, well); } } } @@ -186,6 +148,29 @@ public: } private: + template + void linearizeSingleWell(SparseMatrixAdapter& jacobian, + GlobalEqVector& res, + const WellType& well) + { + if (model_.addMatrixContributions()) { + well->addWellContributions(jacobian); + } + + const auto& cells = well->cells(); + linearize_res_local_.resize(cells.size()); + + for (size_t i = 0; i < cells.size(); ++i) { + linearize_res_local_[i] = res[cells[i]]; + } + + well->apply(linearize_res_local_); + + for (size_t i = 0; i < cells.size(); ++i) { + res[cells[i]] = linearize_res_local_[i]; + } + } + Model& model_; GlobalEqVector linearize_res_local_{}; Parallel::Communication lin_comm_;