Fix pressure renormalization conditions. Some formatting fixes.

Checking for both rock compressibility and pressure conditions
is necessary before we renormalize pressure.

examples/sim_2p_incomp_reorder.cpp

@@ -277,7 +277,6 @@ main(int argc, char** argv)
         }
     }
 
-    epoch_os.close();
     std::cout << "\n\n================    End of simulation     ===============\n\n";
     rep.report(std::cout);
 
@@ -285,7 +284,6 @@ main(int argc, char** argv)
       std::string filename = output_dir + "/walltime.param";
       std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
       rep.reportParam(tot_os);
-      tot_os.close();
     }
 
 }

opm/core/simulator/SimulatorIncompTwophase.cpp

@@ -96,6 +96,7 @@ namespace Opm
         WellsManager& wells_manager_;
         const Wells* wells_;
         const std::vector<double>& src_;
+        const FlowBoundaryConditions* bcs_;
         // Solvers
         IncompTpfa psolver_;
         TransportModelTwophase tsolver_;
@@ -224,6 +225,35 @@ namespace Opm
     }
 
 
+    static bool allNeumannBCs(const FlowBoundaryConditions* bcs)
+    {
+        if (bcs == NULL) {
+            return true;
+        } else {
+            return std::find(bcs->type, bcs->type + bcs->nbc, BC_PRESSURE)
+                == bcs->type + bcs->nbc;
+        }
+    }
+
+
+    static bool allRateWells(const Wells* wells)
+    {
+        if (wells == NULL) {
+            return true;
+        }
+        const int nw = wells->number_of_wells;
+        for (int w = 0; w < nw; ++w) {
+            const WellControls* wc = wells->ctrls[w];
+            if (wc->current >= 0) {
+                if (wc->type[wc->current] == BHP) {
+                    return false;
+                }
+            }
+        }
+        return true;
+    }
+
+
 
 
 
@@ -242,6 +272,7 @@ namespace Opm
           wells_manager_(wells_manager),
           wells_(wells_manager.c_wells()),
           src_(src),
+          bcs_(bcs),
           psolver_(grid, props, rock_comp, linsolver,
                    param.getDefault("nl_pressure_residual_tolerance", 0.0),
                    param.getDefault("nl_pressure_change_tolerance", 1.0),
@@ -363,8 +394,12 @@ namespace Opm
                 std::vector<double> initial_pressure = state.pressure();
                 psolver_.solve(timer.currentStepLength(), state, well_state);
 
-                // Renormalize if incompressible.
-                if (!rock_comp_->isActive()) {
+                // Renormalize pressure if rock is incompressible, and
+                // there are no pressure conditions (bcs or wells).
+                // It is deemed sufficient for now to renormalize
+                // using geometric volume instead of pore volume.
+                if ((rock_comp_ == NULL || !rock_comp_->isActive())
+                    && allNeumannBCs(bcs_) && allRateWells(wells_)) {
                     // Compute average pressures of previous and last
                     // step, and total volume.
                     double av_prev_press = 0.0;
@@ -381,7 +416,8 @@ namespace Opm
                     for (int cell = 0; cell < num_cells; ++cell) {
                         state.pressure()[cell] += ren_const;
                     }
-                    for (int well = 0; well < wells_->number_of_wells; ++well) {
+                    const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells;
+                    for (int well = 0; well < num_wells; ++well) {
                         well_state.bhp()[well] += ren_const;
                     }
                 }
@@ -488,8 +524,6 @@ namespace Opm
             if (output_) {
                 sreport.reportParam(tstep_os);
             }
-
-
         }
 
         if (output_) {