Added preliminary sketch of compressible tpfa solver.

This solver is:
 - using the residual based assembler,
 - aiming to include the nonlinear Newton iterations
   (therefore also the (re-)evaluation of fluid data).

opm/core/pressure/CompressibleTpfa.cpp

@@ -0,0 +1,189 @@
+/*
+  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/CompressibleTpfa.hpp>
+#include <opm/core/pressure/tpfa/cfs_tpfa_residual.h>
+#include <opm/core/pressure/tpfa/compr_quant_general.h>
+#include <opm/core/pressure/tpfa/compr_source.h>
+#include <opm/core/pressure/tpfa/trans_tpfa.h>
+#include <opm/core/linalg/LinearSolverInterface.hpp>
+#include <opm/core/linalg/sparse_sys.h>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/newwells.h>
+
+#include <algorithm>
+
+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.
+    /// \param[in] wells         The wells argument. Will be used in solution, 
+    ///                          is ignored if NULL
+    CompressibleTpfa::CompressibleTpfa(const UnstructuredGrid& g,
+                                       const double* permeability,
+                                       const double* /* gravity */, // ???
+                                       const LinearSolverInterface& linsolver,
+                                       const struct Wells* wells,
+                                       const int num_phases)
+	: grid_(g),
+          linsolver_(linsolver),
+	  htrans_(g.cell_facepos[ g.number_of_cells ]),
+	  trans_ (g.number_of_faces),
+          wells_(wells)
+    {
+        if (wells_ && (wells_->number_of_phases != num_phases)) {
+            THROW("Inconsistent number of phases specified: "
+                  << wells_->number_of_phases << " != " << num_phases);
+        }
+        const int num_dofs = g.number_of_cells + wells ? wells->number_of_wells : 0;
+        pressure_increment_.resize(num_dofs);
+	UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
+	tpfa_htrans_compute(gg, permeability, &htrans_[0]);
+	tpfa_trans_compute(gg, &htrans_[0], &trans_[0]);
+        cfs_tpfa_res_wells w;
+        w.W = const_cast<struct Wells*>(wells_);
+        w.data = NULL;
+	h_ = cfs_tpfa_res_construct(gg, &w, num_phases);
+    }
+
+
+
+
+    /// Destructor.
+    CompressibleTpfa::~CompressibleTpfa()
+    {
+	cfs_tpfa_res_destroy(h_);
+    }
+
+
+
+
+    /// Solve pressure equation, by Newton iterations.
+    void CompressibleTpfa::solve()
+    {
+        // Set up dynamic data.
+        computeDynamicData();
+
+        // Assemble J and F.
+        assemble();
+
+        bool residual_ok = false; // Replace with tolerance check.
+        while (!residual_ok) {
+            // Solve for increment in Newton method:
+            //   incr = x_{n+1} - x_{n} = -J^{-1}F
+            // (J is Jacobian matrix, F is residual)
+            solveIncrement();
+
+            // Update pressure vars with increment.
+
+            // Set up dynamic data.
+            computeDynamicData();
+
+            // Assemble J and F.
+            assemble();
+
+            // Check for convergence.
+            // Include both tolerance check for residual
+            // and solution change.
+        }
+
+        // Write to output parameters.
+        // computeResults(...);
+    }
+
+
+
+
+    /// Solve pressure equation, by Newton iterations.
+    void CompressibleTpfa::computeDynamicData()
+    {
+    }
+
+
+
+
+    /// Compute the residual and Jacobian.
+    void CompressibleTpfa::assemble()
+    {
+        // Arguments or members?
+        const double dt = 0.0;
+        const double* z = NULL;
+        const double* cell_press = NULL;
+        const double* well_bhp = NULL;
+        const double* porevol = NULL;
+
+	UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
+        CompletionData completion_data;
+        completion_data.gpot = 0; // TODO
+        completion_data.A = 0; // TODO
+        completion_data.phasemob = 0; // TODO
+        cfs_tpfa_res_wells wells_tmp;
+        wells_tmp.W = const_cast<Wells*>(wells_);
+        wells_tmp.data = &completion_data;
+        cfs_tpfa_res_forces forces;
+        forces.wells = &wells_tmp;
+        forces.src = NULL; // Check if it is legal to leave it as NULL.
+        compr_quantities_gen cq;
+        cq.Ac = 0; // TODO
+        cq.dAc = 0; // TODO
+        cq.Af = 0; // TODO
+        cq.phasemobf = 0; // TODO
+        cq.voldiscr = 0; // TODO
+        // TODO: gravcapf_ must be set already.
+        cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0],
+                              &gravcapf_[0], &cell_press[0], &well_bhp[0],
+                              &porevol[0], h_);
+    }
+
+
+
+
+    /// Computes pressure_increment_.
+    void CompressibleTpfa::solveIncrement()
+    {
+        // Increment is equal to -J^{-1}F
+	linsolver_.solve(h_->J, h_->F, &pressure_increment_[0]);        
+        std::transform(pressure_increment_.begin(), pressure_increment_.end(),
+                       pressure_increment_.begin(), std::negate<double>());
+    }
+
+
+
+
+    /// Compute the output.
+    void CompressibleTpfa::computeResults(std::vector<double>& // pressure
+                                          ,
+                                          std::vector<double>& // faceflux
+                                          ,
+                                          std::vector<double>& // well_bhp
+                                          ,
+                                          std::vector<double>& // well_rate
+                                          )
+    {
+    }
+
+} // namespace Opm

opm/core/pressure/CompressibleTpfa.hpp

@@ -0,0 +1,102 @@
+/*
+  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_COMPRESSIBLETPFA_HEADER_INCLUDED
+#define OPM_COMPRESSIBLETPFA_HEADER_INCLUDED
+
+
+#include <vector>
+
+struct UnstructuredGrid;
+struct cfs_tpfa_res_data;
+struct Wells;
+struct FlowBoundaryConditions;
+
+namespace Opm
+{
+
+    class LinearSolverInterface;
+
+    /// Encapsulating a tpfa pressure solver for the compressible-fluid case.
+    /// Supports gravity, wells 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.
+    class CompressibleTpfa
+    {
+    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.
+        /// \param[in] wells         The wells argument. Will be used in solution, 
+        ///                          is ignored if NULL
+        /// \param[in] num_phases    Must be 2 or 3.
+	CompressibleTpfa(const UnstructuredGrid& g,
+                         const double* permeability,
+                         const double* gravity,
+                         const LinearSolverInterface& linsolver,
+                         const struct Wells* wells,
+                         const int num_phases);
+
+	/// Destructor.
+	~CompressibleTpfa();
+
+        void solve();
+
+    private:
+        void computeDynamicData();
+        void assemble();
+        void solveIncrement();
+
+	void computeResults(std::vector<double>& pressure,
+                            std::vector<double>& faceflux,
+                            std::vector<double>& well_bhp,
+                            std::vector<double>& well_rate);
+
+        // ------ Data that will remain unmodified after construction. ------
+	const UnstructuredGrid& grid_;
+        const LinearSolverInterface& linsolver_;
+	std::vector<double> htrans_;
+	std::vector<double> trans_ ;
+        const Wells* wells_;   // Outside may modify controls (only) between calls to solve().
+
+        // ------ Internal data for the cfs_tpfa_res solver. ------
+	struct cfs_tpfa_res_data* h_;
+
+        // ------ Data that will be modified for every solve. ------
+
+        // ------ Data that will be modified for every solver iteration. ------
+        // The update to be applied to the pressures (cell and bhp).
+        std::vector<double> pressure_increment_;
+
+
+
+        // Gravity and capillary contributions (per face).
+        std::vector<double> gravcapf_;
+
+
+    };
+
+} // namespace Opm
+
+
+#endif // OPM_COMPRESSIBLETPFA_HEADER_INCLUDED