Added Newton compressible fluid solver.

opm/core/pressure/IncompTpfa.cpp

@@ -26,6 +26,7 @@
 #include <opm/core/linalg/sparse_sys.h>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/newwells.h>
+#include <string.h>
 
 namespace Opm
 {
@@ -233,9 +234,108 @@ namespace Opm
             soln.well_flux = &well_rate[0];
             soln.well_press = &well_bhp[0];
         }
+
         ifs_tpfa_press_flux(gg, &F, &trans_[0], h_, &soln);
     }
 
     
+    void IncompTpfa::solveIncrement(const std::vector<double>& totmob,
+				    const std::vector<double>& omega,
+				    const std::vector<double>& src,
+				    const std::vector<double>& wdp,
+				    const FlowBoundaryConditions* bcs,
+				    const std::vector<double>& porevol,
+				    const std::vector<double>& rock_comp,
+				    const std::vector<double>& prev_pressure,
+				    const std::vector<double>& initial_porevol,
+				    const double dt,
+				    std::vector<double>& pressure_increment)
+
+    {
+	UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
+	tpfa_eff_trans_compute(gg, &totmob[0], &htrans_[0], &trans_[0]);
+
+	if (!omega.empty()) {
+	    if (gpress_.empty()) {
+		THROW("Nozero omega argument given, but gravity was null in constructor.");
+	    }
+	    mim_ip_density_update(gg->number_of_cells, gg->cell_facepos,
+				  &omega[0],
+				  &gpress_[0], &gpress_omegaweighted_[0]);
+	} else {
+	    if (!gpress_.empty()) {
+		THROW("Empty omega argument given, but gravity was non-null in constructor.");
+	    }
+	}
+
+        ifs_tpfa_forces F = { NULL, NULL, wells_, NULL, NULL };
+        if (! src.empty()) { F.src = &src[0]; }
+        F.bc = bcs;
+        F.totmob = &totmob[0];
+        F.wdp = &wdp[0];
+        
+        if (rock_comp.empty()) {
+            ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_);
+        } else {
+            ifs_tpfa_assemble_comprock_increment(gg, &F, &trans_[0], &gpress_omegaweighted_[0],
+                                       &porevol[0], &rock_comp[0], dt, &prev_pressure[0], 
+				       &initial_porevol[0], h_);
+        }
+
+	linsolver_.solve(h_->A, h_->b, &pressure_increment[0]);
+
+    }
+
+    void IncompTpfa::computeFaceFlux(const std::vector<double>& totmob,
+				     const std::vector<double>& omega,
+				     const std::vector<double>& src,
+				     const std::vector<double>& wdp,
+				     const FlowBoundaryConditions* bcs,
+				     std::vector<double>& pressure,
+				     std::vector<double>& faceflux,
+				     std::vector<double>& well_bhp,
+				     std::vector<double>& well_rate) {
+
+	UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
+	tpfa_eff_trans_compute(gg, &totmob[0], &htrans_[0], &trans_[0]);
+
+	if (!omega.empty()) {
+	    if (gpress_.empty()) {
+		THROW("Nozero omega argument given, but gravity was null in constructor.");
+	    }
+	    mim_ip_density_update(gg->number_of_cells, gg->cell_facepos,
+				  &omega[0],
+				  &gpress_[0], &gpress_omegaweighted_[0]);
+	} else {
+	    if (!gpress_.empty()) {
+		THROW("Empty omega argument given, but gravity was non-null in constructor.");
+	    }
+	}
+
+        ifs_tpfa_forces F = { NULL, NULL, wells_, NULL, NULL };
+        if (! src.empty()) { F.src = &src[0]; }
+        F.bc = bcs;
+        F.totmob = &totmob[0];
+        F.wdp = &wdp[0];
+
+	ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_);
+	
+	faceflux.resize(grid_.number_of_faces);
+
+        ifs_tpfa_solution soln = { NULL, NULL, NULL, NULL };
+        soln.cell_press = &pressure[0];
+        soln.face_flux  = &faceflux[0];
+
+        if(wells_ != NULL) {
+            well_bhp.resize(wells_->number_of_wells);
+            well_rate.resize(wells_->well_connpos[ wells_->number_of_wells ]);
+            soln.well_flux = &well_rate[0];
+            soln.well_press = &well_bhp[0];
+        }
+
+	memcpy(h_->x, &pressure[0], grid_.number_of_cells * sizeof *(h_->x));
+
+        ifs_tpfa_press_flux(gg, &F, &trans_[0], h_, &soln);
+    }
 
 } // namespace Opm

opm/core/pressure/IncompTpfa.hpp

@@ -121,6 +121,28 @@ namespace Opm
                    std::vector<double>& well_bhp,
                    std::vector<double>& well_rate);
 
+        void solveIncrement(const std::vector<double>& totmob,
+			    const std::vector<double>& omega,
+			    const std::vector<double>& src,
+			    const std::vector<double>& wdp,
+			    const FlowBoundaryConditions* bcs,
+			    const std::vector<double>& porevol,
+			    const std::vector<double>& rock_comp,
+			    const std::vector<double>& prev_pressure,
+			    const std::vector<double>& initial_porevol,
+			    const double dt,
+			    std::vector<double>& pressure_increment);
+
+	void computeFaceFlux(const std::vector<double>& totmob,
+					 const std::vector<double>& omega,
+					 const std::vector<double>& src,
+					 const std::vector<double>& wdp,
+					 const FlowBoundaryConditions* bcs,
+					 std::vector<double>& pressure,
+					 std::vector<double>& faceflux,
+					 std::vector<double>& well_bhp,
+					 std::vector<double>& well_rate);
+
         /// Expose read-only reference to internal half-transmissibility.
         const ::std::vector<double>& getHalfTrans() const { return htrans_; }
 

opm/core/pressure/tpfa/ifs_tpfa.c

@@ -11,6 +11,16 @@
 #include <opm/core/pressure/tpfa/ifs_tpfa.h>
 
 
+void mult_csr_matrix(const struct CSRMatrix* A, const double* u, double* v) {
+    size_t i,j;
+    for (j = 0; j < A->m; ++j) {
+	v[j] = 0;
+	for (i = (size_t) (A->ia[j]); i < (size_t) (A->ia[j+1]); ++i) {
+	    v[j] += A->sa[i]*u[A->ja[i]];
+	}
+    }
+}
+
 struct ifs_tpfa_impl {
     double *fgrav;              /* Accumulated grav contrib/face */
 
@@ -725,6 +735,49 @@ ifs_tpfa_assemble_comprock(struct UnstructuredGrid      *G        ,
     return ok;
 }
 
+/* ---------------------------------------------------------------------- */
+int
+ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid      *G        ,
+				     const struct ifs_tpfa_forces *F        ,
+				     const double                 *trans    ,
+				     const double                 *gpress   ,
+				     const double                 *porevol  ,
+				     const double                 *rock_comp,
+				     const double                  dt       ,
+				     const double                 *prev_pressure,
+				     const double                 *initial_porevolume,
+				     struct ifs_tpfa_data         *h        )
+/* ---------------------------------------------------------------------- */
+{
+    int c;
+    size_t j;
+    int system_singular, ok;
+    double *v;
+
+    assemble_incompressible(G, F, trans, gpress, h, &system_singular, &ok);
+
+    v = malloc(h->A->m * sizeof *v);
+
+    mult_csr_matrix(h->A, prev_pressure, v);
+
+    /* We want to solve a Newton step for the residual
+     * (porevol(pressure)-porevol(initial_pressure))/dt + residual_for_imcompressible 
+     *
+     */
+
+    if (ok) {
+	for (c = 0; c < G->number_of_cells; ++c) {
+            j = csrmatrix_elm_index(c, c, h->A);
+	    h->A->sa[j] += porevol[c] * rock_comp[c] / dt;
+            h->b[c]     += -(porevol[c] - initial_porevolume[c])/dt - v[c];
+        }
+    }
+
+    free(v);
+
+    return ok;
+}
+
 
 /* ---------------------------------------------------------------------- */
 void

opm/core/pressure/tpfa/ifs_tpfa.h

@@ -61,6 +61,10 @@ struct ifs_tpfa_data *
 ifs_tpfa_construct(struct UnstructuredGrid *G,
                    struct Wells            *W);
 
+
+void mult_csr_matrix(const struct CSRMatrix* A, const double* u, double* v);
+
+
 int
 ifs_tpfa_assemble(struct UnstructuredGrid      *G     ,
                   const struct ifs_tpfa_forces *F     ,
@@ -78,6 +82,19 @@ ifs_tpfa_assemble_comprock(struct UnstructuredGrid      *G        ,
                            const double                  dt       ,
                            const double                 *pressure ,
                            struct ifs_tpfa_data         *h        );
+int
+ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid      *G        ,
+				     const struct ifs_tpfa_forces *F        ,
+				     const double                 *trans    ,
+				     const double                 *gpress   ,
+				     const double                 *porevol  ,
+				     const double                 *rock_comp,
+				     const double                  dt       ,
+				     const double                 *prev_pressure ,
+				     const double                 *initial_porevolume,
+				     struct ifs_tpfa_data         *h        );
+
+
 void
 ifs_tpfa_press_flux(struct UnstructuredGrid      *G    ,
                     const struct ifs_tpfa_forces *F    ,