Added class IncompTpfa.

Makefile.am

@@ -74,6 +74,7 @@ opm/core/pressure/msmfem/ifsh_ms.c                      \
 opm/core/pressure/msmfem/dfs.c                          \
 opm/core/pressure/msmfem/coarse_sys.c                   \
 opm/core/pressure/ifsh.c                                \
+opm/core/pressure/IncompTpfa.cpp                        \
 opm/core/pressure/mimetic/mimetic.c                     \
 opm/core/pressure/mimetic/hybsys_global.c               \
 opm/core/pressure/mimetic/hybsys.c                      \
@@ -157,6 +158,7 @@ opm/core/GridAdapter.hpp                                \
 opm/core/GridManager.hpp                                \
 opm/core/pressure/fsh.h                                 \
 opm/core/pressure/HybridPressureSolver.hpp              \
+opm/core/pressure/IncompTpfa.hpp                        \
 opm/core/pressure/TPFAPressureSolver.hpp                \
 opm/core/pressure/tpfa/compr_quant.h                    \
 opm/core/pressure/tpfa/compr_source.h                   \

opm/core/pressure/IncompTpfa.cpp

@@ -0,0 +1,112 @@
+/*
+  Copyright 2012 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <opm/core/pressure/IncompTpfa.hpp>
+#include <opm/core/pressure/tpfa/ifs_tpfa.h>
+#include <opm/core/pressure/tpfa/trans_tpfa.h>
+#include <opm/core/pressure/mimetic/mimetic.h>
+#include <opm/core/transport/CSRMatrixUmfpackSolver.hpp>
+
+
+namespace Opm
+{
+
+
+
+    /// Construct solver.
+    /// \param[in] g             A 2d or 3d grid.
+    /// \param[in] permeability  Array of permeability tensors, the array
+    ///                          should have size N*D^2, if D == g.dimensions
+    ///                          and N == g.number_of_cells.
+    /// \param[in] gravity       Gravity vector. If nonzero, the array should
+    ///                          have D elements.
+    IncompTpfa::IncompTpfa(const UnstructuredGrid& g,
+			   const double* permeability,
+			   const double* gravity)
+	: grid_(g),
+	  htrans_(g.cell_facepos[ g.number_of_cells ]),
+	  trans_ (g.number_of_faces),
+	  gpress_(g.cell_facepos[ g.number_of_cells ], 0.0),
+	  gpress_omegaweighted_(g.cell_facepos[ g.number_of_cells ], 0.0)
+    {
+	UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
+	tpfa_htrans_compute(gg, permeability, &htrans_[0]);
+	if (gravity) {
+	    mim_ip_compute_gpress(gg->number_of_cells, gg->dimensions, gravity,
+				  gg->cell_facepos, gg->cell_faces,
+				  gg->face_centroids, gg->cell_centroids,
+				  &gpress_[0]);
+	}
+	h_ = ifs_tpfa_construct(gg);
+    }
+
+
+
+
+    /// Destructor.
+    IncompTpfa::~IncompTpfa()
+    {
+	ifs_tpfa_destroy(h_);
+    }
+
+
+
+    /// Assemble and solve pressure system.
+    /// \param[in]  totmob     Must contain N total mobility values (one per cell).
+    ///                        totmob = \sum_{p} kr_p/mu_p.
+    /// \param[in]  omega      Must contain N mobility-weighted density values (one per cell).
+    ///                        omega = \frac{\sum_{p} mob_p rho_p}{\sum_p rho_p}.
+    /// \param[in]  src        Must contain N souce rates (one per cell).
+    ///                        Positive values represent total inflow rates,
+    ///                        negative values represent total outflow rates.
+    /// \param[out] pressure   Will contain N cell-pressure values.
+    /// \param[out] faceflux   Will contain F signed face flux values.
+    void IncompTpfa::solve(const std::vector<double>& totmob,
+		      const std::vector<double>& omega,
+		      const std::vector<double>& src,
+		      std::vector<double>& pressure,
+		      std::vector<double>& faceflux)
+    {
+	UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
+	tpfa_eff_trans_compute(gg, &totmob[0], &htrans_[0], &trans_[0]);
+
+	if (!omega.empty()) {
+	    mim_ip_density_update(gg->number_of_cells, gg->cell_facepos,
+				  &omega[0],
+				  &gpress_[0], &gpress_omegaweighted_[0]);
+	}
+
+	ifs_tpfa_assemble(gg, &trans_[0], &src[0], &gpress_omegaweighted_[0], h_);
+
+	using Opm::ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver;
+	CSRMatrixUmfpackSolver linsolve;
+	linsolve.solve(h_->A, h_->b, h_->x);
+
+	pressure.resize(grid_.number_of_cells);
+	faceflux.resize(grid_.number_of_faces);
+
+	ifs_tpfa_press_flux(gg, &trans_[0], h_,
+			    &pressure[0],
+			    &faceflux[0]);
+    }
+
+
+
+
+} // namespace Opm

opm/core/pressure/IncompTpfa.hpp

@@ -0,0 +1,83 @@
+/*
+  Copyright 2012 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_INCOMPTPFA_HEADER_INCLUDED
+#define OPM_INCOMPTPFA_HEADER_INCLUDED
+
+
+#include <vector>
+
+struct UnstructuredGrid;
+struct ifs_tpfa_data;
+
+
+namespace Opm
+{
+
+    /// Encapsulating a tpfa pressure solver for the incompressible case.
+    /// Supports gravity and simple sources as driving forces.
+    /// Below we use the shortcuts D for the number of dimensions, N
+    /// for the number of cells and F for the number of faces.
+    /// Note: we intend to add wells in the future.
+    class IncompTpfa
+    {
+    public:
+	/// Construct solver.
+	/// \param[in] g             A 2d or 3d grid.
+	/// \param[in] permeability  Array of permeability tensors, the array
+	///                          should have size N*D^2, if D == g.dimensions
+	///                          and N == g.number_of_cells.
+	/// \param[in] gravity       Gravity vector. If nonzero, the array should
+	///                          have D elements.
+	IncompTpfa(const UnstructuredGrid& g,
+		   const double* permeability,
+		   const double* gravity);
+
+	/// Destructor.
+	~IncompTpfa();
+
+	/// Assemble and solve pressure system.
+	/// \param[in]  totmob     Must contain N total mobility values (one per cell).
+	///                        totmob = \sum_{p} kr_p/mu_p.
+	/// \param[in]  omega      Must contain N mobility-weighted density values (one per cell).
+	///                        omega = \frac{\sum_{p} mob_p rho_p}{\sum_p rho_p}.
+	/// \param[in]  src        Must contain N souce rates (one per cell).
+	///                        Positive values represent total inflow rates,
+	///                        negative values represent total outflow rates.
+	/// \param[out] pressure   Will contain N cell-pressure values.
+	/// \param[out] faceflux   Will contain F signed face flux values.
+	void solve(const std::vector<double>& totmob,
+		   const std::vector<double>& omega,
+		   const std::vector<double>& src,
+		   std::vector<double>& pressure,
+		   std::vector<double>& faceflux);
+
+    private:
+	const UnstructuredGrid& grid_;
+	::std::vector<double> htrans_;
+	::std::vector<double> trans_ ;
+	::std::vector<double> gpress_;
+	::std::vector<double> gpress_omegaweighted_;
+
+	struct ifs_tpfa_data* h_;
+    };
+
+} // namespace Opm
+
+#endif // OPM_INCOMPTPFA_HEADER_INCLUDED