Merged.

examples/wells_example.cpp

@@ -4,6 +4,7 @@
 
 
 #include "opm/core/utility/initState.hpp"
+#include "opm/core/utility/SimulatorTimer.hpp"
 #include <opm/core/WellsManager.hpp>
 #include <opm/core/GridManager.hpp>
 #include <opm/core/pressure/IncompTpfa.hpp>
@@ -14,43 +15,51 @@
 #include <opm/core/TwophaseState.hpp>
 #include <opm/core/pressure/FlowBCManager.hpp>
 #include <opm/core/linalg/LinearSolverFactory.hpp>
-int main(int argc, char** argv) {
+#include <opm/core/fluid/RockCompressibility.hpp>
+
+int main(int argc, char** argv)
+{
 
     using namespace Opm::parameter;
     using namespace Opm;
-    ParameterGroup parameters( argc, argv, false );
-    std::string file_name = parameters.getDefault<std::string>("inputdeck", "data.data");
+    ParameterGroup parameters(argc, argv, false);
+    std::string file_name = parameters.getDefault<std::string > ("inputdeck", "data.data");
+
+    SimulatorTimer simtimer;
+    simtimer.init(parameters);
 
     // Read input file
     EclipseGridParser parser(file_name);
     std::cout << "Done!" << std::endl;
     // Setup grid
     GridManager grid(parser);
-    
+
     // Finally handle the wells
     WellsManager wells(parser, *grid.c_grid(), NULL);
-    
+
     std::vector<int> global_cells(grid.c_grid()->global_cell, grid.c_grid()->global_cell + grid.c_grid()->number_of_cells);
 
     double gravity[3] = {0.0, 0.0, parameters.getDefault<double>("gravity", 0.0)};
     IncompPropertiesFromDeck incomp_properties(parser, global_cells);
 
+    RockCompressibility rock_comp(parser);
+
     Opm::LinearSolverFactory linsolver(parameters);
 
     // EXPERIMENT_ISTL
-    IncompTpfa pressure_solver(*grid.c_grid(), incomp_properties.permeability(), 
-                               gravity, linsolver,  wells.c_wells());
-    
-    
+    IncompTpfa pressure_solver(*grid.c_grid(), incomp_properties.permeability(),
+            gravity, linsolver, wells.c_wells());
+
+
     std::vector<int> all_cells;
-    for(int i = 0; i < grid.c_grid()->number_of_cells; i++) {
+    for (int i = 0; i < grid.c_grid()->number_of_cells; i++) {
         all_cells.push_back(i);
     }
-            
+
     Opm::TwophaseState state;
-    
+
     initStateTwophaseFromDeck(*grid.c_grid(), incomp_properties, parser, gravity[2], state);
-    
+
     // Compute phase mobilities
     std::vector<double> phase_mob;
     computePhaseMobilities(incomp_properties, all_cells, state.saturation(), phase_mob);
@@ -58,10 +67,10 @@ int main(int argc, char** argv) {
     std::vector<double> totmob;
     std::vector<double> omega;
     computeTotalMobilityOmega(incomp_properties, all_cells, state.saturation(), totmob, omega);
-    
+
     std::vector<double> wdp;
     computeWDP(*wells.c_wells(), *grid.c_grid(), state.saturation(), incomp_properties.density(), gravity[2], true, wdp);
-    
+
     std::vector<double> src;
     Opm::FlowBCManager bcs;
 
@@ -70,87 +79,114 @@ int main(int argc, char** argv) {
 
     std::vector<double> well_bhp;
     std::vector<double> well_rate_per_cell;
-    pressure_solver.solve(totmob, omega, src, wdp, bcs.c_bcs(), pressure, face_flux, well_bhp, well_rate_per_cell);
-    std::cout << "Solved" << std::endl;
+    std::vector<double> rc;
+    rc.resize(grid.c_grid()->number_of_cells);
     
-    for(size_t i = 0; i < well_rate_per_cell.size(); i++) {
-        std::cout << well_rate_per_cell[i] << std::endl;
+    int nl_pressure_maxiter = 100;
+    double nl_pressure_tolerance = 0.0;
+    if (rock_comp.isActive()) {
+        nl_pressure_maxiter = parameters.getDefault("nl_pressure_maxiter", 10);
+        nl_pressure_tolerance = parameters.getDefault("nl_pressure_tolerance", 1.0); // in Pascal
     }
-    std::vector<double> well_rate;
-    
-    // This will be refactored into a separate function once done.
+
+    const int num_cells = grid.c_grid()->number_of_cells;
+    std::vector<double> porevol;
+    if (rock_comp.isActive()) {
+        computePorevolume(*grid.c_grid(), incomp_properties, rock_comp, state.pressure(), porevol);
+    } else {
+        computePorevolume(*grid.c_grid(), incomp_properties, porevol);
+    }
+    if (rock_comp.isActive()) {
+        std::vector<double> initial_pressure = state.pressure();
+        std::vector<double> prev_pressure;
+        for (int iter = 0; iter < nl_pressure_maxiter; ++iter) {
+            prev_pressure = state.pressure();
+            for (int cell = 0; cell < num_cells; ++cell) {
+                rc[cell] = rock_comp.rockComp(state.pressure()[cell]);
+            }
+            state.pressure() = initial_pressure;
+            pressure_solver.solve(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc, simtimer.currentStepLength(),
+                    state.pressure(), state.faceflux(), well_bhp, well_rate_per_cell);
+            double max_change = 0.0;
+            for (int cell = 0; cell < num_cells; ++cell) {
+                max_change = std::max(max_change, std::fabs(state.pressure()[cell] - prev_pressure[cell]));
+            }
+            std::cout << "Pressure iter " << iter << "   max change = " << max_change << std::endl;
+            if (max_change < nl_pressure_tolerance) {
+                break;
+            }
+        }
+        computePorevolume(*grid.c_grid(), incomp_properties, rock_comp, state.pressure(), porevol);
+    } else {
+        pressure_solver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(),
+                well_bhp, well_rate_per_cell);
+    }
+
     const int np = incomp_properties.numPhases();
     std::vector<double> fractional_flows(grid.c_grid()->number_of_cells*np, 0.0);
-    for (int cell = 0; cell < grid.c_grid()->number_of_cells; ++cell) {
-        double phase_sum = 0.0;
-        for (int phase = 0; phase < np; ++phase) {
-            phase_sum += phase_mob[cell*np + phase];
-        }
-        for (int phase = 0; phase < np; ++phase) {
-            fractional_flows[cell*np + phase] = phase_mob[cell*np + phase] / phase_sum;
-        }
-    }
-    // End stuff that needs to be refactored into a seperated function
-    
-    computeFlowRatePerWell(*wells.c_wells(), well_rate_per_cell, well_rate);
-    
+    computeFractionalFlow(incomp_properties, all_cells, state.saturation(), fractional_flows);
+
     // This will be refactored into a separate function once done
     std::vector<double> well_resflows(wells.c_wells()->number_of_wells*np, 0.0);
-    for ( int wix = 0; wix < wells.c_wells()->number_of_wells; ++wix) {
-        for (int i = wells.c_wells()->well_connpos[wix]; i < wells.c_wells()->well_connpos[wix+1]; ++i) {
-            const int cell = wells.c_wells()->well_cells[i];
-            for (int phase = 0; phase < np; ++phase) {
-                well_resflows[wix*np + phase] += well_rate_per_cell[i]*fractional_flows[cell*np + phase];
-            }
-        }
-    }
-
+    computePhaseFlowRatesPerWell(*wells.c_wells(), well_rate_per_cell, fractional_flows, well_resflows);
     // We approximate (for _testing_ that resflows = surfaceflows)
-    for (int iter = 0; iter < 10 && !wells.conditionsMet(well_bhp, well_resflows, well_resflows); ++iter) {
+    for (int wc_iter = 0; wc_iter < 10 && !wells.conditionsMet(well_bhp, well_resflows, well_resflows); ++wc_iter) {
         std::cout << "Conditions not met for well, trying again" << std::endl;
-        pressure_solver.solve(totmob, omega, src, wdp, bcs.c_bcs(), pressure, face_flux, well_bhp, well_rate_per_cell);
-        std::cout << "Solved" << std::endl;
-
-
-        for (int wix = 0; wix < wells.c_wells()->number_of_wells; ++wix) {
-            for (int phase = 0; phase < np; ++phase) {
-                // Reset
-                well_resflows[wix * np + phase] = 0.0;
-            }
-            for (int i = wells.c_wells()->well_connpos[wix]; i < wells.c_wells()->well_connpos[wix + 1]; ++i) {
-                const int cell = wells.c_wells()->well_cells[i];
-                for (int phase = 0; phase < np; ++phase) {
-                    well_resflows[wix * np + phase] += well_rate_per_cell[i] * fractional_flows[cell * np + phase];
+        if (rock_comp.isActive()) {
+            std::vector<double> initial_pressure = state.pressure();
+            std::vector<double> prev_pressure;
+            for (int iter = 0; iter < nl_pressure_maxiter; ++iter) {
+                prev_pressure = state.pressure();
+                for (int cell = 0; cell < num_cells; ++cell) {
+                    rc[cell] = rock_comp.rockComp(state.pressure()[cell]);
+                }
+                state.pressure() = initial_pressure;
+                pressure_solver.solve(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc, simtimer.currentStepLength(),
+                        state.pressure(), state.faceflux(), well_bhp, well_rate_per_cell);
+                double max_change = 0.0;
+                for (int cell = 0; cell < num_cells; ++cell) {
+                    max_change = std::max(max_change, std::fabs(state.pressure()[cell] - prev_pressure[cell]));
+                }
+                std::cout << "Pressure iter " << iter << "   max change = " << max_change << std::endl;
+                if (max_change < nl_pressure_tolerance) {
+                    break;
                 }
             }
+            computePorevolume(*grid.c_grid(), incomp_properties, rock_comp, state.pressure(), porevol);
+        } else {
+            pressure_solver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(),
+                    well_bhp, well_rate_per_cell);
         }
+        std::cout << "Solved" << std::endl;
+
+        computePhaseFlowRatesPerWell(*wells.c_wells(), well_rate_per_cell, fractional_flows, well_resflows);
     }
 
 #if 0
     std::vector<double> porevol;
     computePorevolume(*grid->c_grid(), incomp_properties, porevol);
-    
-    
+
+
 
     TwophaseFluid fluid(incomp_properties);
-    TransportModel  model  (fluid, *grid->c_grid(), porevol, gravity[2], true);
+    TransportModel model(fluid, *grid->c_grid(), porevol, gravity[2], true);
 
     TransportSolver tsolver(model);
 
     TransportSource* tsrc = create_transport_source(2, 2);
-    double ssrc[]   = { 1.0, 0.0 };
-    double ssink[]  = { 0.0, 1.0 };
-    double zdummy[] = { 0.0, 0.0 };
-    
+    double ssrc[] = {1.0, 0.0};
+    double ssink[] = {0.0, 1.0};
+    double zdummy[] = {0.0, 0.0};
+
     {
         int well_cell_index = 0;
         for (int well = 0; well < wells.c_wells()->number_of_wells; ++well) {
-            for( int cell = wells.c_wells()->well_connpos[well]; cell < wells.c_wells()->well_connpos[well + 1]; ++cell) {
+            for (int cell = wells.c_wells()->well_connpos[well]; cell < wells.c_wells()->well_connpos[well + 1]; ++cell) {
                 if (well_rate_per_cell[well_cell_index] > 0.0) {
-                    append_transport_source(well_cell_index, 2, 0, 
+                    append_transport_source(well_cell_index, 2, 0,
                             well_rate_per_cell[well_cell_index], ssrc, zdummy, tsrc);
                 } else if (well_rate_per_cell[well_cell_index] < 0.0) {
-                    append_transport_source(well_cell_index, 2, 0, 
+                    append_transport_source(well_cell_index, 2, 0,
                             well_rate_per_cell[well_cell_index], ssink, zdummy, tsrc);
                 }
             }
@@ -158,7 +194,7 @@ int main(int argc, char** argv) {
     }
 
     tsolver.solve(*grid->c_grid(), tsrc, stepsize, ctrl, state, linsolve, rpt);
-     
+
     Opm::computeInjectedProduced(*props, state.saturation(), src, stepsize, injected, produced);
 #endif
     return 0;

opm/core/WellsGroup.cpp

@@ -1000,6 +1000,7 @@ namespace Opm
                         production_specification.gas_max_rate_ = line.gas_max_rate_;
                         production_specification.liquid_max_rate_ = line.liquid_max_rate_;
                         production_specification.procedure_ = toProductionProcedure(line.procedure_);
+                        production_specification.reservoir_flow_max_rate_ = line.resv_max_rate_;
                     }
                 }
             }

opm/core/eclipse/SpecialEclipseFields.hpp

@@ -1211,11 +1211,13 @@ struct GconprodLine
     double water_max_rate_;        // Water rate target or upper limit
     double gas_max_rate_;          // Gas rate target or upper limit
     double liquid_max_rate_;       // Liquid rate target or upper limit
+    double resv_max_rate_;         // Reservoir liquid rate target or upper limit
     std::string procedure_;        // Procedure on exceeding a maximum rate limit
     // Default values
     GconprodLine() :
 	oil_max_rate_(-1.0E20), water_max_rate_(-1.0E20),
-	gas_max_rate_(-1.0E20), liquid_max_rate_(-1.0E20)
+	gas_max_rate_(-1.0E20), liquid_max_rate_(-1.0E20),
+        resv_max_rate_(-1.0E20)
     {
     }
 };
@@ -1269,11 +1271,42 @@ struct GCONPROD : public SpecialBase
                     is.putback('/');
                     procedure = "NONE";
                 }
+                int read_ignored = 0;
+                int to_read = 6;
+                for (; read_ignored < to_read; ++read_ignored) {
+                    std::string ignored_value = readString(is);
+                    if (ignored_value[ignored_value.size()-1]=='*') {
+                        // we've got defaulted argument, increment
+                        
+                        int num_defaulted;
+                        std::string num_defaulted_str = ignored_value.substr(0, ignored_value.size()-1);
+                        std::istringstream(num_defaulted_str) >> num_defaulted;
+                        read_ignored += num_defaulted;
+                        if (read_ignored >= to_read) {
+                            break;
+                        }
+                    } else if(ignored_value[0] == '/') {
+                        is.putback('/');
+                        break;
+                    }
+                }
+                if ( read_ignored <= to_read) {
+                    // Not defaulted, we can read
+                    std::string reservoir_volume = readString(is);
+                    if (reservoir_volume[reservoir_volume.size()-1] == '*') {
+                        // Defaulted, we're not doing anything
+                    } else if(reservoir_volume[0] == '/') {
+                        is.putback('/');
+                    } else {
+                        std::istringstream(reservoir_volume) >> gconprod_line.resv_max_rate_;
+                    }
+                }
             } else {
                 procedure = "NONE";
             }
             gconprod_line.procedure_ = procedure;
 
+            
             gconprod.push_back(gconprod_line);
             // HACK! Ignore any further items
             if (num_read == num_to_read) {
@@ -1294,7 +1327,8 @@ struct GCONPROD : public SpecialBase
 	       << gconprod[i].water_max_rate_ << "  " 
 	       << gconprod[i].gas_max_rate_ << "  " 
 	       << gconprod[i].liquid_max_rate_ << "  " 
-	       << gconprod[i].procedure_
+	       << gconprod[i].procedure_ << " " 
+               << gconprod[i].resv_max_rate_
 	       << std::endl; 
 	}
 	os << std::endl;
@@ -1304,11 +1338,14 @@ struct GCONPROD : public SpecialBase
     {
 	double lrat = units.liqvol_s / units.time;
 	double grat = units.gasvol_s / units.time;
+      	double resv = units.liqvol_r / units.time;
+
 	for (int i=0; i<(int) gconprod.size(); ++i) {
 	    gconprod[i].oil_max_rate_ *= lrat;
 	    gconprod[i].water_max_rate_ *= lrat;
 	    gconprod[i].gas_max_rate_ *= grat;
 	    gconprod[i].liquid_max_rate_ *= lrat;
+            gconprod[i].resv_max_rate_ *= resv;
 	}
     }
 };

opm/core/utility/miscUtilities.cpp

@@ -271,6 +271,32 @@ namespace Opm
         }
     }
 
+    /// Computes the fractional flow for each cell in the cells argument
+    /// @param[in] props                rock and fluid properties
+    /// @param[in] cells                cells with which the saturation values are associated
+    /// @param[in] saturations          saturation values (for all phases)
+    /// @param[out] fractional_flow     the fractional flow for each phase for each cell.
+
+    void computeFractionalFlow(const Opm::IncompPropertiesInterface& props,
+                               const std::vector<int>& cells,
+                               const std::vector<double>& saturations,
+                               std::vector<double>& fractional_flows)
+    {
+        const int num_phases = props.numPhases();
+        std::vector<double> pc_mobs(cells.size() * num_phases);
+        computePhaseMobilities(props, cells, saturations, pc_mobs);
+        fractional_flows.resize(cells.size() * num_phases);
+        for (size_t i = 0; i < cells.size(); ++i) {
+            double phase_sum = 0.0;
+            for (int phase = 0; phase < num_phases; ++phase) {
+                phase_sum += pc_mobs[i * num_phases + phase];
+            }
+            for (int phase = 0; phase < num_phases; ++phase) {
+                fractional_flows[i * num_phases + phase] = pc_mobs[i * num_phases + phase] / phase_sum;
+            }
+        }
+    }
+
     /// Compute two-phase transport source terms from face fluxes,
     /// and pressure equation source terms. This puts boundary flows
     /// into the source terms for the transport equation.
@@ -495,6 +521,34 @@ namespace Opm
     }
 
 
+    /// Computes the phase flow rate per well
+    /// \param[in] wells                        The wells for which the flow rate should be computed
+    /// \param[in] flow_rates_per_well_cell     The total flow rate for each cell (ordered the same
+    ///                                         way as the wells struct
+    /// \param[in] fractional_flows             the fractional flow for each cell in each well
+    /// \param[out] phase_flow_per_well         Will contain the phase flow per well
+    void computePhaseFlowRatesPerWell(const Wells& wells,
+                                      const std::vector<double>& flow_rates_per_well_cell,
+                                      const std::vector<double>& fractional_flows,
+                                      std::vector<double>& phase_flow_per_well)
+    {
+        const int np = wells.number_of_phases;
+        const int nw = wells.number_of_wells;
+        phase_flow_per_well.resize(nw * np);
+        for (int wix = 0; wix < nw; ++wix) {
+            for (int phase = 0; phase < np; ++phase) {
+                // Reset vector
+                phase_flow_per_well[wix*np + phase] = 0.0;
+            }
+            for (int i = wells.well_connpos[wix]; i < wells.well_connpos[wix + 1]; ++i) {
+                const int cell = wells.well_cells[i];
+                for (int phase = 0; phase < np; ++phase) {
+                    phase_flow_per_well[wix * np + phase] += flow_rates_per_well_cell[i] * fractional_flows[cell * np + phase];
+                }
+            }
+        }
+    }
+
 
     void Watercut::push(double time, double fraction, double produced)
     {

opm/core/utility/miscUtilities.hpp

@@ -129,6 +129,17 @@ namespace Opm
                                 const std::vector<int>&               cells,
                                 const std::vector<double>&            s    ,
                                 std::vector<double>&                  pmobc);
+    
+
+    /// Computes the fractional flow for each cell in the cells argument
+    /// @param[in] props                rock and fluid properties
+    /// @param[in] cells                cells with which the saturation values are associated
+    /// @param[in] saturations          saturation values (for all phases)
+    /// @param[out] fractional_flow     the fractional flow for each phase for each cell.
+    void computeFractionalFlow(const Opm::IncompPropertiesInterface& props,
+                               const std::vector<int>& cells,
+                               const std::vector<double>& saturations,
+                               std::vector<double>& fractional_flows);
 
 
     /// Compute two-phase transport source terms from face fluxes,
@@ -203,6 +214,18 @@ namespace Opm
     void computeFlowRatePerWell(const Wells& wells, const std::vector<double>& flow_rates_per_cell,
                                 std::vector<double>& flow_rates_per_well);
     
+    /// Computes the phase flow rate per well
+    /// \param[in] wells The wells for which the flow rate should be computed
+    /// \param[in] flow_rates_per_well_cell The total flow rate for each cell (ordered the same
+    ///                                 way as the wells struct
+    /// \param[in] fractional_flows    the fractional flow for each cell in each well
+    /// \param[out] phase_flow_per_well Will contain the phase flow per well
+    void computePhaseFlowRatesPerWell(const Wells& wells, 
+                                      const std::vector<double>& flow_rates_per_well_cell,
+                                      const std::vector<double>& fractional_flows,
+                                      std::vector<double>& phase_flow_per_well);
+ 
+    
     /// Encapsulates the watercut curves.
     class Watercut
     {