Merge pull request #622 from totto82/refact_sovlent

Refactor the solvent model

opm/autodiff/BlackoilModelBase.hpp

@@ -387,11 +387,19 @@ namespace Opm {
         void computeWellConnectionPressures(const SolutionState& state,
                                             const WellState& xw);
 
+        void computePropertiesForWellConnectionPressures(const SolutionState& state,
+                                                         const WellState& xw,
+                                                         std::vector<double>& b_perf,
+                                                         std::vector<double>& rsmax_perf,
+                                                         std::vector<double>& rvmax_perf,
+                                                         std::vector<double>& surf_dens_perf);
+
         void
         assembleMassBalanceEq(const SolutionState& state);
 
         void
-        extractWellPerfProperties(std::vector<ADB>& mob_perfcells,
+        extractWellPerfProperties(const SolutionState& state,
+                                  std::vector<ADB>& mob_perfcells,
                                   std::vector<ADB>& b_perfcells) const;
 
         bool

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -795,18 +795,14 @@ namespace detail {
         }
     }
 
-
-
     template <class Grid, class Implementation>
-    void BlackoilModelBase<Grid, Implementation>::computeWellConnectionPressures(const SolutionState& state,
+    void BlackoilModelBase<Grid, Implementation>::computePropertiesForWellConnectionPressures(const SolutionState& state,
-                                                                                 const WellState& xw)
+                                                                                              const WellState& xw,
+                                                                                              std::vector<double>& b_perf,
+                                                                                              std::vector<double>& rsmax_perf,
+                                                                                              std::vector<double>& rvmax_perf,
+                                                                                              std::vector<double>& surf_dens_perf)
     {
-        if( ! localWellsActive() ) return ;
-
-        using namespace Opm::AutoDiffGrid;
-        // 1. Compute properties required by computeConnectionPressureDelta().
-        //    Note that some of the complexity of this part is due to the function
-        //    taking std::vector<double> arguments, and not Eigen objects.
         const int nperf = wells().well_connpos[wells().number_of_wells];
         const int nw = wells().number_of_wells;
 
@@ -837,8 +833,6 @@ namespace detail {
         }
         const PhaseUsage& pu = fluid_.phaseUsage();
         DataBlock b(nperf, pu.num_phases);
-        std::vector<double> rsmax_perf(nperf, 0.0);
-        std::vector<double> rvmax_perf(nperf, 0.0);
         if (pu.phase_used[BlackoilPhases::Aqua]) {
             const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value();
             b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
@@ -851,6 +845,8 @@ namespace detail {
             b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
             const V rssat = fluidRsSat(avg_press, perf_so, well_cells);
             rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
+        } else {
+            rsmax_perf.assign(nperf, 0.0);
         }
         if (pu.phase_used[BlackoilPhases::Vapour]) {
             const ADB perf_rv = subset(state.rv, well_cells);
@@ -858,13 +854,12 @@ namespace detail {
             b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
             const V rvsat = fluidRvSat(avg_press, perf_so, well_cells);
             rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
+        } else {
+            rvmax_perf.assign(nperf, 0.0);
         }
         // b is row major, so can just copy data.
-        std::vector<double> b_perf(b.data(), b.data() + nperf * pu.num_phases);
+        b_perf.assign(b.data(), b.data() + nperf * pu.num_phases);
-        // Extract well connection depths.
+
-        const V depth = cellCentroidsZToEigen(grid_);
-        const V pdepth = subset(depth, well_cells);
-        std::vector<double> perf_depth(pdepth.data(), pdepth.data() + nperf);
 
         // Surface density.
         // The compute density segment wants the surface densities as
@@ -873,10 +868,25 @@ namespace detail {
         for (int phase = 1; phase < pu.num_phases; ++phase) {
             rho += superset(fluid_.surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases);
         }
-        std::vector<double> surf_dens_perf(rho.data(), rho.data() + nperf * pu.num_phases);
+        surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases);
+    }
 
-        // Gravity
+
-        double grav = detail::getGravity(geo_.gravity(), dimensions(grid_));
+    template <class Grid, class Implementation>
+    void BlackoilModelBase<Grid, Implementation>::computeWellConnectionPressures(const SolutionState& state,
+                                                                                 const WellState& xw)
+    {
+        if( ! localWellsActive() ) return ;
+
+        using namespace Opm::AutoDiffGrid;
+        // 1. Compute properties required by computeConnectionPressureDelta().
+        //    Note that some of the complexity of this part is due to the function
+        //    taking std::vector<double> arguments, and not Eigen objects.
+        std::vector<double> b_perf;
+        std::vector<double> rsmax_perf;
+        std::vector<double> rvmax_perf;
+        std::vector<double> surf_dens_perf;
+        asImpl().computePropertiesForWellConnectionPressures(state, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
                 
         // 2. Compute densities
         std::vector<double> cd =
@@ -884,6 +894,17 @@ namespace detail {
                         wells(), xw, fluid_.phaseUsage(),
                         b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
 
+        const int nperf = wells().well_connpos[wells().number_of_wells];
+        const std::vector<int>& well_cells = wops_.well_cells;
+
+        // Extract well connection depths.
+        const V depth = cellCentroidsZToEigen(grid_);
+        const V pdepth = subset(depth, well_cells);
+        std::vector<double> perf_depth(pdepth.data(), pdepth.data() + nperf);
+
+        // Gravity
+        double grav = detail::getGravity(geo_.gravity(), dimensions(grid_));
+
         // 3. Compute pressure deltas
         std::vector<double> cdp =
                 WellDensitySegmented::computeConnectionPressureDelta(
@@ -954,7 +975,7 @@ namespace detail {
 
         std::vector<ADB> mob_perfcells;
         std::vector<ADB> b_perfcells;
-        asImpl().extractWellPerfProperties(mob_perfcells, b_perfcells);
+        asImpl().extractWellPerfProperties(state, mob_perfcells, b_perfcells);
         if (param_.solve_welleq_initially_ && initial_assembly) {
             // solve the well equations as a pre-processing step
             asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state);
@@ -1101,7 +1122,8 @@ namespace detail {
 
     template <class Grid, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::extractWellPerfProperties(std::vector<ADB>& mob_perfcells,
+    BlackoilModelBase<Grid, Implementation>::extractWellPerfProperties(const SolutionState&,
+                                                                       std::vector<ADB>& mob_perfcells,
                                                                        std::vector<ADB>& b_perfcells) const
     {
         // If we have wells, extract the mobilities and b-factors for

opm/autodiff/BlackoilMultiSegmentModel_impl.hpp

@@ -630,7 +630,7 @@ namespace Opm {
 
         std::vector<ADB> mob_perfcells;
         std::vector<ADB> b_perfcells;
-        asImpl().extractWellPerfProperties(mob_perfcells, b_perfcells);
+        asImpl().extractWellPerfProperties(state, mob_perfcells, b_perfcells);
         if (param_.solve_welleq_initially_ && initial_assembly) {
             // solve the well equations as a pre-processing step
             asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state);

opm/autodiff/BlackoilSolventModel.hpp

@@ -87,14 +87,6 @@ namespace Opm {
                          ReservoirState& reservoir_state,
                          WellState& well_state);
 
-        /// Assemble the residual and Jacobian of the nonlinear system.
-        /// \param[in]      reservoir_state   reservoir state variables
-        /// \param[in, out] well_state        well state variables
-        /// \param[in]      initial_assembly  pass true if this is the first call to assemble() in this timestep
-        void assemble(const ReservoirState& reservoir_state,
-                      WellState& well_state,
-                      const bool initial_assembly);
-
 
     protected:
 
@@ -159,6 +151,7 @@ namespace Opm {
         using Base::updateWellControls;
         using Base::computeWellConnectionPressures;
         using Base::addWellControlEq;
+        using Base::computePropertiesForWellConnectionPressures;
 
         std::vector<ADB>
         computeRelPerm(const SolutionState& state) const;
@@ -212,8 +205,12 @@ namespace Opm {
                                            const SolutionState& state,
                                            WellState& xw);
 
-        void computeWellConnectionPressures(const SolutionState& state,
+        void computePropertiesForWellConnectionPressures(const SolutionState& state,
-                                            const WellState& xw);
+                                                         const WellState& xw,
+                                                         std::vector<double>& b_perf,
+                                                         std::vector<double>& rsmax_perf,
+                                                         std::vector<double>& rvmax_perf,
+                                                         std::vector<double>& surf_dens_perf);
 
         void updateEquationsScaling();
 
@@ -229,6 +226,11 @@ namespace Opm {
         const std::vector<PhasePresence>
         phaseCondition() const {return this->phaseCondition_;}
 
+        void extractWellPerfProperties(const SolutionState& state,
+                                       std::vector<ADB>& mob_perfcells,
+                                       std::vector<ADB>& b_perfcells);
+
+
         // compute effective viscosities (mu_eff_) and effective b factors (b_eff_)  using the ToddLongstaff model
         void computeEffectiveProperties(const SolutionState&  state);
 

opm/autodiff/BlackoilSolventModel_impl.hpp

@@ -135,7 +135,7 @@ namespace Opm {
         std::vector<V> vars0 = Base::variableStateInitials(x, xw);
         assert(int(vars0.size()) == fluid_.numPhases() + 2);
 
-        // Initial polymer concentration.
+        // Initial solvent concentration.
         if (has_solvent_) {
             const auto& solvent_saturation = x.getCellData( BlackoilSolventState::SSOL );
             const int nc = solvent_saturation.size();
@@ -257,6 +257,10 @@ namespace Opm {
     BlackoilSolventModel<Grid>::computeAccum(const SolutionState& state,
                                              const int            aix  )
     {
+
+        if (is_miscible_) {
+            computeEffectiveProperties(state);
+        }
         Base::computeAccum(state, aix);
 
         // Compute accumulation of the solvent
@@ -378,10 +382,13 @@ namespace Opm {
     }
 
     template <class Grid>
-    void BlackoilSolventModel<Grid>::computeWellConnectionPressures(const SolutionState& state,
+    void BlackoilSolventModel<Grid>::computePropertiesForWellConnectionPressures(const SolutionState& state,
-                                                                        const WellState& xw)
+                                                                                 const WellState& xw,
+                                                                                 std::vector<double>& b_perf,
+                                                                                 std::vector<double>& rsmax_perf,
+                                                                                 std::vector<double>& rvmax_perf,
+                                                                                 std::vector<double>& surf_dens_perf)
     {
-        if( ! Base::localWellsActive() ) return ;
 
         using namespace Opm::AutoDiffGrid;
         // 1. Compute properties required by computeConnectionPressureDelta().
@@ -417,8 +424,7 @@ namespace Opm {
 
         const PhaseUsage& pu = fluid_.phaseUsage();
         DataBlock b(nperf, pu.num_phases);
-        std::vector<double> rsmax_perf(nperf, 0.0);
+
-        std::vector<double> rvmax_perf(nperf, 0.0);
         const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value();
         if (pu.phase_used[BlackoilPhases::Aqua]) {
             b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
@@ -435,6 +441,8 @@ namespace Opm {
             b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
             const V rssat = fluidRsSat(avg_press, perf_so, well_cells);
             rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
+        } else {
+            rsmax_perf.assign(0.0, nperf);
         }
         V surf_dens_copy = superset(fluid_.surfaceDensity(0, well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases);
         for (int phase = 1; phase < pu.num_phases; ++phase) {
@@ -492,41 +500,18 @@ namespace Opm {
 
             const V rvsat = fluidRvSat(avg_press, perf_so, well_cells);
             rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
+        } else {
+            rvmax_perf.assign(0.0, nperf);
         }
 
         // b and surf_dens_perf is row major, so can just copy data.
-        std::vector<double> b_perf(b.data(), b.data() + nperf * pu.num_phases);        
+        b_perf.assign(b.data(), b.data() + nperf * pu.num_phases);
-        std::vector<double> surf_dens_perf(surf_dens_copy.data(), surf_dens_copy.data() + nperf * pu.num_phases);
+        surf_dens_perf.assign(surf_dens_copy.data(), surf_dens_copy.data() + nperf * pu.num_phases);
-
-        // Extract well connection depths.
-        const V depth = cellCentroidsZToEigen(grid_);
-        const V pdepth = subset(depth, well_cells);
-        std::vector<double> perf_depth(pdepth.data(), pdepth.data() + nperf);
-
-        // Gravity
-        double grav = detail::getGravity(geo_.gravity(), dimensions(grid_));
-
-        // 2. Compute densities
-        std::vector<double> cd =
-                WellDensitySegmented::computeConnectionDensities(
-                        wells(), xw, fluid_.phaseUsage(),
-                        b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
-
-        // 3. Compute pressure deltas
-        std::vector<double> cdp =
-                WellDensitySegmented::computeConnectionPressureDelta(
-                        wells(), perf_depth, cd, grav);
-
-        // 4. Store the results
-        Base::well_perforation_densities_ = Eigen::Map<const V>(cd.data(), nperf);
-        Base::well_perforation_pressure_diffs_ = Eigen::Map<const V>(cdp.data(), nperf);
     }
 
 
 
 
-
-
     template <class Grid>
     void BlackoilSolventModel<Grid>::updateState(const V& dx,
                                                   ReservoirState& reservoir_state,
@@ -771,6 +756,44 @@ namespace Opm {
 
     }
 
+    template <class Grid>
+    void
+    BlackoilSolventModel<Grid>::extractWellPerfProperties(const SolutionState& state,
+                                                          std::vector<ADB>& mob_perfcells,
+                                                          std::vector<ADB>& b_perfcells)
+    {
+        Base::extractWellPerfProperties(state, mob_perfcells, b_perfcells);
+        if (has_solvent_) {
+            int gas_pos = fluid_.phaseUsage().phase_pos[Gas];
+            const std::vector<int>& well_cells = wops_.well_cells;
+            const int nperf = well_cells.size();
+            // Gas and solvent is combinded and solved together
+            // The input in the well equation is then the
+            // total gas phase = hydro carbon gas + solvent gas
+
+            // The total mobility is the sum of the solvent and gas mobiliy
+            mob_perfcells[gas_pos] += subset(rq_[solvent_pos_].mob, well_cells);
+
+            // A weighted sum of the b-factors of gas and solvent are used.
+            const int nc = Opm::AutoDiffGrid::numCells(grid_);
+
+            const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+            const ADB zero = ADB::constant(V::Zero(nc));
+            const ADB& ss = state.solvent_saturation;
+            const ADB& sg = (active_[ Gas ]
+                             ? state.saturation[ pu.phase_pos[ Gas ] ]
+                             : zero);
+
+            Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
+            ADB F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
+            V ones = V::Constant(nperf,1.0);
+
+            b_perfcells[gas_pos] = (ones - F_solvent) * b_perfcells[gas_pos];
+            b_perfcells[gas_pos] += (F_solvent * subset(rq_[solvent_pos_].b, well_cells));
+
+        }
+    }
+
     template <class Grid>
     void
     BlackoilSolventModel<Grid>::computeEffectiveProperties(const SolutionState&    state)
@@ -965,114 +988,20 @@ namespace Opm {
                      const ADB& ss) const
     {
         std::vector<ADB> pressures = Base::computePressures(po, sw, so, sg);
-        // The imiscible capillary pressure is evaluated using the total gas saturation (sg + ss)
-        std::vector<ADB> pressures_imisc = Base::computePressures(po, sw, so, sg + ss);
-
-        // Pressure effects on capillary pressure miscibility
-        const ADB pmisc = solvent_props_.pressureMiscibilityFunction(po, cells_);
-        // Only the pcog is effected by the miscibility. Since pg = po + pcog, changing pg is eqvivalent
-        // to changing the gas pressure directly.
-        const int  nc = cells_.size();
-        const V ones = V::Constant(nc, 1.0);
-        pressures[Gas] = ( pmisc * pressures[Gas] + ((ones - pmisc) * pressures_imisc[Gas]));
-        return pressures;
-    }
-
-
-    template <class Grid>
-    void
-    BlackoilSolventModel<Grid>::assemble(const ReservoirState& reservoir_state,
-                                         WellState& well_state,
-                                         const bool initial_assembly)
-    {
-
-        using namespace Opm::AutoDiffGrid;
-
-        // Possibly switch well controls and updating well state to
-        // get reasonable initial conditions for the wells
-        updateWellControls(well_state);
-
-        // Create the primary variables.
-        SolutionState state = variableState(reservoir_state, well_state);
-
-        if (initial_assembly) {
-            // Create the (constant, derivativeless) initial state.
-            SolutionState state0 = state;
-            makeConstantState(state0);
-            // Compute initial accumulation contributions
-            // and well connection pressures.
-            if (is_miscible_) {
-                computeEffectiveProperties(state0);
-            }
-
-            computeAccum(state0, 0);
-            computeWellConnectionPressures(state0, well_state);
-        }
-        if (is_miscible_) {
-            computeEffectiveProperties(state);
-        }
-
-        // -------- Mass balance equations --------
-        assembleMassBalanceEq(state);
-
-        // -------- Well equations ----------
-        if ( ! wellsActive() ) {
-            return;
-        }
-
-        V aliveWells;
-
-        const int np = wells().number_of_phases;
-        std::vector<ADB> cq_s(np, ADB::null());
-
-        const int nw = wells().number_of_wells;
-        const int nperf = wells().well_connpos[nw];
-        const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
-
-        std::vector<ADB> mob_perfcells(np, ADB::null());
-        std::vector<ADB> b_perfcells(np, ADB::null());
-        for (int phase = 0; phase < np; ++phase) {
-            mob_perfcells[phase] = subset(rq_[phase].mob, well_cells);
-            b_perfcells[phase] = subset(rq_[phase].b, well_cells);
-        }
 
         if (has_solvent_) {
-            int gas_pos = fluid_.phaseUsage().phase_pos[Gas];
+            // The imiscible capillary pressure is evaluated using the total gas saturation (sg + ss)
-            // Gas and solvent is combinded and solved together
+            std::vector<ADB> pressures_imisc = Base::computePressures(po, sw, so, sg + ss);
-            // The input in the well equation is then the
-            // total gas phase = hydro carbon gas + solvent gas
-
-            // The total mobility is the sum of the solvent and gas mobiliy
-            mob_perfcells[gas_pos] += subset(rq_[solvent_pos_].mob, well_cells);
-
-            // A weighted sum of the b-factors of gas and solvent are used.
-            const int nc = Opm::AutoDiffGrid::numCells(grid_);
-
-            const Opm::PhaseUsage& pu = fluid_.phaseUsage();
-            const ADB zero = ADB::constant(V::Zero(nc));
-            const ADB& ss = state.solvent_saturation;
-            const ADB& sg = (active_[ Gas ]
-                             ? state.saturation[ pu.phase_pos[ Gas ] ]
-                             : zero);
-
-            Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
-            ADB F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
-            V ones = V::Constant(nperf,1.0);
-
-            b_perfcells[gas_pos] = (ones - F_solvent) * b_perfcells[gas_pos];
-            b_perfcells[gas_pos] += (F_solvent * subset(rq_[solvent_pos_].b, well_cells));
 
+            // Pressure effects on capillary pressure miscibility
+            const ADB pmisc = solvent_props_.pressureMiscibilityFunction(po, cells_);
+            // Only the pcog is effected by the miscibility. Since pg = po + pcog, changing pg is eqvivalent
+            // to changing the gas pressure directly.
+            const int  nc = cells_.size();
+            const V ones = V::Constant(nc, 1.0);
+            pressures[Gas] = ( pmisc * pressures[Gas] + ((ones - pmisc) * pressures_imisc[Gas]));
         }
-        if (param_.solve_welleq_initially_ && initial_assembly) {
+        return pressures;
-            // solve the well equations as a pre-processing step
-            Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state);
-        }
-        Base::computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
-        Base::updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
-        Base::addWellFluxEq(cq_s, state);
-        addWellContributionToMassBalanceEq(cq_s, state, well_state);
-        Base::addWellControlEq(state, well_state, aliveWells);
-
     }