Merge from upstream.

opm/core/utility/miscUtilities.cpp

@@ -45,8 +45,33 @@ namespace Opm
     }
 
 
+    /// @brief Computes total saturated volumes over all grid cells.
+    /// @param[in]  pv        the pore volume by cell.
+    /// @param[in]  s         saturation values (for all P phases)
+    /// @param[out] sat_vol   must point to a valid array with P elements,
+    ///                       where P = s.size()/pv.size().
+    ///                       For each phase p, we compute
+    ///                       sat_vol_p = sum_i s_p_i pv_i
+    void computeSaturatedVol(const std::vector<double>& pv,
+			     const std::vector<double>& s,
+			     double* sat_vol)
+    {
+	const int num_cells = pv.size();
+	const int np = s.size()/pv.size();
+	if (int(s.size()) != num_cells*np) {
+	    THROW("Sizes of s and pv vectors do not match.");
+	}
+	std::fill(sat_vol, sat_vol + np, 0.0);
+	for (int c = 0; c < num_cells; ++c) {
+	    for (int p = 0; p < np; ++p) {
+		sat_vol[p] += pv[c]*s[np*c + p];
+	    }
+	}
+    }
+
+
     /// @brief Computes average saturations over all grid cells.
-    /// @param[out] pv        the pore volume by cell.
+    /// @param[in]  pv        the pore volume by cell.
     /// @param[in]  s         saturation values (for all P phases)
     /// @param[out] aver_sat  must point to a valid array with P elements,
     ///                       where P = s.size()/pv.size().
@@ -78,6 +103,54 @@ namespace Opm
     }
 
 
+    /// @brief Computes injected and produced volumes of all phases.
+    /// Note 1: assumes that only the first phase is injected.
+    /// Note 2: assumes that transport has been done with an
+    ///         implicit method, i.e. that the current state
+    ///         gives the mobilities used for the preceding timestep.
+    /// @param[in]  props     fluid and rock properties.
+    /// @param[in]  s         saturation values (for all P phases)
+    /// @param[in]  src       if < 0: total outflow, if > 0: first phase inflow.
+    /// @param[in]  dt        timestep used
+    /// @param[out] injected  must point to a valid array with P elements,
+    ///                       where P = s.size()/src.size().
+    /// @param[out] produced  must also point to a valid array with P elements.
+    void computeInjectedProduced(const IncompPropertiesInterface& props,
+				 const std::vector<double>& s,
+				 const std::vector<double>& src,
+				 const double dt,
+				 double* injected,
+				 double* produced)
+    {
+	const int num_cells = src.size();
+	const int np = s.size()/src.size();
+	if (int(s.size()) != num_cells*np) {
+	    THROW("Sizes of s and src vectors do not match.");
+	}
+	std::fill(injected, injected + np, 0.0);
+	std::fill(produced, produced + np, 0.0);
+	const double* visc = props.viscosity();
+	std::vector<double> mob(np);
+	for (int c = 0; c < num_cells; ++c) {
+	    if (src[c] > 0.0) {
+		injected[0] += src[c]*dt;
+	    } else if (src[c] < 0.0) {
+		const double flux = -src[c]*dt;
+		const double* sat = &s[np*c];
+		props.relperm(1, sat, &c, &mob[0], 0);
+		double totmob = 0.0;
+		for (int p = 0; p < np; ++p) {
+		    mob[p] /= visc[p];
+		    totmob += mob[p];
+		}
+		for (int p = 0; p < np; ++p) {
+		    produced[p] += (mob[p]/totmob)*flux;
+		}
+	    }
+	}
+    }
+
+
 
     /// @brief Computes total mobility for a set of saturation values.
     /// @param[in]  props     rock and fluid properties
@@ -164,6 +237,48 @@ namespace Opm
         }
     }
 
+    /// 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.
+    /// \param[in]  grid          The grid used.
+    /// \param[in]  src           Pressure eq. source terms. The sign convention is:
+    ///                           (+) positive  total inflow (positive velocity divergence)
+    ///                           (-) negative  total outflow
+    /// \param[in]  faceflux      Signed face fluxes, typically the result from a flow solver.
+    /// \param[in]  inflow_frac   Fraction of inflow that consists of first phase.
+    ///                           Example: if only water is injected, inflow_frac == 1.0.
+    ///                           Note: it is not possible (with this method) to use different fractions
+    ///                           for different inflow sources, be they source terms of boundary flows.
+    /// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign:
+    ///                           (+) positive  inflow of first phase (water)
+    ///                           (-) negative  total outflow of both phases
+    void computeTransportSource(const UnstructuredGrid& grid,
+				const std::vector<double>& src,
+				const std::vector<double>& faceflux,
+				const double inflow_frac,
+				std::vector<double>& transport_src)
+    {
+	int nc = grid.number_of_cells;
+	transport_src.resize(nc);
+	for (int c = 0; c < nc; ++c) {
+	    transport_src[c] = 0.0;
+	    transport_src[c] += src[c] > 0.0 ? inflow_frac*src[c] : src[c];
+	    for (int hf = grid.cell_facepos[c]; hf < grid.cell_facepos[c + 1]; ++hf) {
+		int f = grid.cell_faces[hf];
+		const int* f2c = &grid.face_cells[2*f];
+		double bdy_influx = 0.0;
+		if (f2c[0] == c && f2c[1] == -1) {
+		    bdy_influx = -faceflux[f];
+		} else if (f2c[0] == -1 && f2c[1] == c) {
+		    bdy_influx = faceflux[f];
+		}
+		if (bdy_influx != 0.0) {
+		    transport_src[c] += bdy_influx > 0.0 ? inflow_frac*bdy_influx : bdy_influx;
+		}
+	    }
+	}
+    }
+
     /// @brief Estimates a scalar cell velocity from face fluxes.
     /// @param[in]  grid            a grid
     /// @param[in]  face_flux       signed per-face fluxes

opm/core/utility/miscUtilities.hpp

@@ -38,6 +38,18 @@ namespace Opm
 			   std::vector<double>& porevol);
 
 
+    /// @brief Computes total saturated volumes over all grid cells.
+    /// @param[out] pv        the pore volume by cell.
+    /// @param[in]  s         saturation values (for all P phases)
+    /// @param[out] sat_vol   must point to a valid array with P elements,
+    ///                       where P = s.size()/pv.size().
+    ///                       For each phase p, we compute
+    ///                       sat_vol_p = sum_i s_p_i pv_i
+    void computeSaturatedVol(const std::vector<double>& pv,
+			     const std::vector<double>& s,
+			     double* sat_vol);
+
+
     /// @brief Computes average saturations over all grid cells.
     /// @param[out] pv        the pore volume by cell.
     /// @param[in]  s         saturation values (for all P phases)
@@ -50,6 +62,25 @@ namespace Opm
 			   double* aver_sat);
 
 
+    /// @brief Computes injected and produced volumes of all phases.
+    /// Note 1: assumes that only the first phase is injected.
+    /// Note 2: assumes that transport has been done with an
+    ///         implicit method, i.e. that the current state
+    ///         gives the mobilities used for the preceding timestep.
+    /// @param[in]  props     fluid and rock properties.
+    /// @param[in]  s         saturation values (for all P phases)
+    /// @param[in]  src       if < 0: total outflow, if > 0: first phase inflow.
+    /// @param[in]  dt        timestep used
+    /// @param[out] injected  must point to a valid array with P elements,
+    ///                       where P = s.size()/src.size().
+    /// @param[out] produced  must also point to a valid array with P elements.
+    void computeInjectedProduced(const IncompPropertiesInterface& props,
+				 const std::vector<double>& s,
+				 const std::vector<double>& src,
+				 const double dt,
+				 double* injected,
+				 double* produced);
+
     /// @brief Computes total mobility for a set of saturation values.
     /// @param[in]  props     rock and fluid properties
     /// @param[in]  cells     cells with which the saturation values are associated
@@ -73,11 +104,35 @@ namespace Opm
 				   std::vector<double>& totmob,
 				   std::vector<double>& omega);
 
+
     void computePhaseMobilities(const Opm::IncompPropertiesInterface& props,
                                 const std::vector<int>&               cells,
                                 const std::vector<double>&            s    ,
                                 std::vector<double>&                  pmobc);
 
+
+    /// 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.
+    /// \param[in]  grid          The grid used.
+    /// \param[in]  src           Pressure eq. source terms. The sign convention is:
+    ///                           (+) positive  total inflow (positive velocity divergence)
+    ///                           (-) negative  total outflow
+    /// \param[in]  faceflux      Signed face fluxes, typically the result from a flow solver.
+    /// \param[in]  inflow_frac   Fraction of inflow that consists of first phase.
+    ///                           Example: if only water is injected, inflow_frac == 1.0.
+    ///                           Note: it is not possible (with this method) to use different fractions
+    ///                           for different inflow sources, be they source terms of boundary flows.
+    /// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign:
+    ///                           (+) positive  inflow of first phase (water)
+    ///                           (-) negative  total outflow of both phases
+    void computeTransportSource(const UnstructuredGrid& grid,
+				const std::vector<double>& src,
+				const std::vector<double>& faceflux,
+				const double inflow_frac,
+				std::vector<double>& transport_src);
+
+
     /// @brief Estimates a scalar cell velocity from face fluxes.
     /// @param[in]  grid            a grid
     /// @param[in]  face_flux       signed per-face fluxes
@@ -91,7 +146,7 @@ namespace Opm
     void toWaterSat(const std::vector<double>& sboth,
 		    std::vector<double>& sw);
 
-    /// Make a a vector of interleaved water and oil saturations from
+    /// Make a vector of interleaved water and oil saturations from
     /// a vector of water saturations.
     void toBothSat(const std::vector<double>& sw,
 		   std::vector<double>& sboth);

opm/core/utility/writeVtkData.hpp

@@ -42,7 +42,7 @@ namespace Opm
 		      std::ostream& os);
 
     /// Vtk output for general grids.
-    void writeVtkData(const UnstructuredGrid* grid,
+    void writeVtkData(const UnstructuredGrid& grid,
 		      const DataMap& data,
 		      std::ostream& os);
 } // namespace Opm