Varnish Parts of Well Implementation

In particular

  * Split some long lines
  * Reverse conditions to reduce nesting
  * Mark potentially unused arguments as [[maybe_unused]]
  * Try to remove redundant calculations
  * Mark some objets 'const' where possible

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -76,8 +76,8 @@ namespace Opm {
                                             simulator.gridView().comm())
         , simulator_(simulator)
     {
-        this->terminal_output_ = ((simulator.gridView().comm().rank() == 0) &&
+        this->terminal_output_ = (simulator.gridView().comm().rank() == 0)
-                                   Parameters::Get<Parameters::EnableTerminalOutput>());
+            && Parameters::Get<Parameters::EnableTerminalOutput>();
 
         local_num_cells_ = simulator_.gridView().size(0);
 
@@ -96,13 +96,13 @@ namespace Opm {
         this->alternative_well_rate_init_ =
             Parameters::Get<Parameters::AlternativeWellRateInit>();
 
-        using SourceDataSpan = 
+        using SourceDataSpan = typename
-            typename PAvgDynamicSourceData<Scalar>::template SourceDataSpan<Scalar>;
+            PAvgDynamicSourceData<Scalar>::template SourceDataSpan<Scalar>;
+
         this->wbpCalculationService_
             .localCellIndex([this](const std::size_t globalIndex)
             { return this->compressedIndexForInterior(globalIndex); })
-            .evalCellSource([this](const int                                                               localCell,
+            .evalCellSource([this](const int localCell, SourceDataSpan sourceTerms)
-                                   SourceDataSpan                                                          sourceTerms)
             {
                 using Item = typename SourceDataSpan::Item;
 
@@ -828,30 +828,33 @@ namespace Opm {
     BlackoilWellModel<TypeTag>::
     initializeWellState(const int timeStepIdx)
     {
-        std::vector<Scalar> cellPressures(this->local_num_cells_, 0.0);
+        const auto pressIx = []()
-        std::vector<Scalar> cellTemperatures(this->local_num_cells_, 0.0);
+        {
-        ElementContext elemCtx(simulator_);
+            if (Indices::oilEnabled)   { return FluidSystem::oilPhaseIdx;   }
+            if (Indices::waterEnabled) { return FluidSystem::waterPhaseIdx; }
 
-        const auto& gridView = simulator_.vanguard().gridView();
+            return FluidSystem::gasPhaseIdx;
+        }();
+
+        auto cellPressures = std::vector<Scalar>(this->local_num_cells_, Scalar{0});
+        auto cellTemperatures = std::vector<Scalar>(this->local_num_cells_, Scalar{0});
+
+        auto elemCtx = ElementContext { this->simulator_ };
+        const auto& gridView = this->simulator_.vanguard().gridView();
 
         OPM_BEGIN_PARALLEL_TRY_CATCH();
         for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
             elemCtx.updatePrimaryStencil(elem);
             elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
 
+            const auto ix = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
             const auto& fs = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0).fluidState();
-            // copy of get perfpressure in Standard well except for value
+
-            Scalar& perf_pressure = cellPressures[elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0)];
+            cellPressures[ix] = fs.pressure(pressIx).value();
-            if (Indices::oilEnabled) {
+            cellTemperatures[ix] = fs.temperature(0).value();
-                perf_pressure = fs.pressure(FluidSystem::oilPhaseIdx).value();
-            } else if (Indices::waterEnabled) {
-                perf_pressure = fs.pressure(FluidSystem::waterPhaseIdx).value();
-            } else {
-                perf_pressure = fs.pressure(FluidSystem::gasPhaseIdx).value();
         }
-            cellTemperatures[elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0)] = fs.temperature(0).value();
+        OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::initializeWellState() failed: ",
-        }
+                                   this->simulator_.vanguard().grid().comm());
-        OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::initializeWellState() failed: ", simulator_.vanguard().grid().comm());
 
         this->wellState().init(cellPressures, cellTemperatures, this->schedule(), this->wells_ecl_,
                                this->local_parallel_well_info_, timeStepIdx,

opm/simulators/wells/StandardWellConnections.cpp

@@ -298,9 +298,11 @@ computePropertiesForPressures(const WellState<Scalar>& well_state,
 {
     const int nperf = well_.numPerfs();
     const PhaseUsage& pu = well_.phaseUsage();
-    props.b_perf.resize(nperf * well_.numComponents());
+
-    props.surf_dens_perf.resize(nperf * well_.numComponents());
+    props.b_perf        .resize(nperf * this->well_.numComponents());
-    const auto& ws = well_state.well(well_.indexOfWell());
+    props.surf_dens_perf.resize(nperf * this->well_.numComponents());
+
+    const auto& ws = well_state.well(this->well_.indexOfWell());
 
     static constexpr int Water = BlackoilPhases::Aqua;
     static constexpr int Oil = BlackoilPhases::Liquid;
@@ -310,12 +312,13 @@ computePropertiesForPressures(const WellState<Scalar>& well_state,
     const bool oilPresent = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
     const bool gasPresent = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
 
-    //rs and rv are only used if both oil and gas is present
+    // Rs/Rv are used only if both oil and gas are present.
     if (oilPresent && gasPresent) {
         props.rsmax_perf.resize(nperf);
         props.rvmax_perf.resize(nperf);
     }
-      //rvw is only used if both water and gas is present
+
+    // Rsw/Rvw are used only if both water and gas are present.
     if (waterPresent && gasPresent) {
         props.rvwmax_perf.resize(nperf);
         props.rswmax_perf.resize(nperf);
@@ -323,9 +326,8 @@ computePropertiesForPressures(const WellState<Scalar>& well_state,
 
     // Compute the average pressure in each well block
     const auto& perf_press = ws.perf_data.pressure;
-    auto p_above =  well_.parallelWellInfo().communicateAboveValues(ws.bhp,
+    const auto  p_above    = this->well_.parallelWellInfo()
-                                                                    perf_press.data(),
+        .communicateAboveValues(ws.bhp, perf_press.data(), nperf);
-                                                                    nperf);
 
     for (int perf = 0; perf < nperf; ++perf) {
         const int cell_idx = well_.cells()[perf];
@@ -351,17 +353,23 @@ computePropertiesForPressures(const WellState<Scalar>& well_state,
                     rsw = std::min(rsw, props.rswmax_perf[perf]);
                 }
             }
-            props.b_perf[waterCompIdx + perf * well_.numComponents()] = FluidSystem::waterPvt().inverseFormationVolumeFactor(region_idx, temperature, p_avg, rsw, saltConcentration);
+
+            props.b_perf[waterCompIdx + perf * well_.numComponents()] = FluidSystem::waterPvt()
+                .inverseFormationVolumeFactor(region_idx, temperature, p_avg, rsw, saltConcentration);
         }
 
         if (gasPresent) {
             const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
             const int gaspos = gasCompIdx + perf * well_.numComponents();
+
             Scalar rvw = 0.0;
             Scalar rv = 0.0;
             if (oilPresent) {
-                const Scalar oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]); //in order to handle negative rates in producers
+                // in order to handle negative rates in producers
-                props.rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(region_idx, temperature, p_avg);
+                const Scalar oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]);
+                props.rvmax_perf[perf] = FluidSystem::gasPvt()
+                    .saturatedOilVaporizationFactor(region_idx, temperature, p_avg);
+
                 if (oilrate > 0) {
                     const Scalar gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0);
                     if (gasrate > 0) {
@@ -370,27 +378,41 @@ computePropertiesForPressures(const WellState<Scalar>& well_state,
                     rv = std::min(rv, props.rvmax_perf[perf]);
                 }
             }
+
             if (waterPresent) {
-                const Scalar waterrate = std::abs(ws.surface_rates[pu.phase_pos[Water]]); //in order to handle negative rates in producers
+                // in order to handle negative rates in producers
-                props.rvwmax_perf[perf] = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(region_idx, temperature, p_avg);
+                const Scalar waterrate = std::abs(ws.surface_rates[pu.phase_pos[Water]]);
+                props.rvwmax_perf[perf] = FluidSystem::gasPvt()
+                    .saturatedWaterVaporizationFactor(region_idx, temperature, p_avg);
+
                 if (waterrate > 0) {
-                    const Scalar gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0);
+                    const Scalar gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]])
+                        - (Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0);
+
                     if (gasrate > 0) {
                         rvw = waterrate / gasrate;
                     }
+
                     rvw = std::min(rvw, props.rvwmax_perf[perf]);
                 }
             }
-            props.b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(region_idx, temperature, p_avg, rv, rvw);
+
+            props.b_perf[gaspos] = FluidSystem::gasPvt()
+                .inverseFormationVolumeFactor(region_idx, temperature, p_avg, rv, rvw);
         }
 
         if (oilPresent) {
             const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
             const int oilpos = oilCompIdx + perf * well_.numComponents();
+
             Scalar rs = 0.0;
             if (gasPresent) {
-                props.rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(region_idx, temperature, p_avg);
+                props.rsmax_perf[perf] = FluidSystem::oilPvt()
-                const Scalar gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0);
+                    .saturatedGasDissolutionFactor(region_idx, temperature, p_avg);
+
+                const Scalar gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]])
+                    - (Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0);
+
                 if (gasrate > 0) {
                     const Scalar oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]);
                     if (oilrate > 0) {
@@ -399,7 +421,9 @@ computePropertiesForPressures(const WellState<Scalar>& well_state,
                     rs = std::min(rs, props.rsmax_perf[perf]);
                 }
             }
-            props.b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(region_idx, temperature, p_avg, rs);
+
+            props.b_perf[oilpos] = FluidSystem::oilPvt()
+                .inverseFormationVolumeFactor(region_idx, temperature, p_avg, rs);
         }
 
         // Surface density.
@@ -409,7 +433,8 @@ computePropertiesForPressures(const WellState<Scalar>& well_state,
             }
 
             const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
-            props.surf_dens_perf[well_.numComponents() * perf  + compIdx] = FluidSystem::referenceDensity( phaseIdx, region_idx );
+            props.surf_dens_perf[well_.numComponents() * perf  + compIdx] =
+                FluidSystem::referenceDensity( phaseIdx, region_idx );
         }
 
         // We use cell values for solvent injector
@@ -705,4 +730,4 @@ INSTANTIATE_TYPE(double)
 INSTANTIATE_TYPE(float)
 #endif
 
-}
+} // namespace Opm

opm/simulators/wells/StandardWell_impl.hpp

@@ -2026,18 +2026,25 @@ namespace Opm
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::
-    updateWaterThroughput(const double dt, WellState<Scalar>& well_state) const
+    updateWaterThroughput([[maybe_unused]] const double dt,
+                          WellState<Scalar>&            well_state) const
     {
         if constexpr (Base::has_polymermw) {
-            if (this->isInjector()) {
+            if (!this->isInjector()) {
-                auto& ws = well_state.well(this->index_of_well_);
+                return;
-                auto& perf_water_throughput = ws.perf_data.water_throughput;
-                for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
-                    const Scalar perf_water_vel = this->primary_variables_.value(Bhp + 1 + perf);
-                    // we do not consider the formation damage due to water flowing from reservoir into wellbore
-                    if (perf_water_vel > 0.) {
-                        perf_water_throughput[perf] += perf_water_vel * dt;
             }
+
+            auto& perf_water_throughput = well_state.well(this->index_of_well_)
+                .perf_data.water_throughput;
+
+            for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
+                const Scalar perf_water_vel =
+                    this->primary_variables_.value(Bhp + 1 + perf);
+
+                // we do not consider the formation damage due to water
+                // flowing from reservoir into wellbore
+                if (perf_water_vel > Scalar{0}) {
+                    perf_water_throughput[perf] += perf_water_vel * dt;
                 }
             }
         }

opm/simulators/wells/WellInterface_impl.hpp

@@ -1477,10 +1477,10 @@ namespace Opm
                             const WellState<Scalar>& well_state,
                             DeferredLogger& deferred_logger) const
     {
-        // Check if well is stopped or under zero rate control, either directly or from group
+        // Check if well is stopped or under zero rate control, either
-        return (this->wellIsStopped() || wellUnderZeroRateTarget(simulator,
+        // directly or from group.
-                                                                 well_state,
+        return this->wellIsStopped()
-                                                                 deferred_logger));
+            || this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger);
     }
 
     template<typename TypeTag>