Added sketch of interface to compressible tpfa solver.

.hgsubstate

@@ -1 +1 @@
-00c15b54de2e838a2cb85ff20cb9d4b8edde1f34 dune/porsol/opmpressure
+335bb8e27ad82ccf4b60c2a0042041a704d50691 dune/porsol/opmpressure

dune/porsol/mimetic/TpfaCompressible.hpp

@@ -0,0 +1,477 @@
+/*
+  Copyright 2010 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_TPFACOMPRESSIBLE_HEADER_INCLUDED
+#define OPM_TPFACOMPRESSIBLE_HEADER_INCLUDED
+
+
+#include "../opmpressure/src/TPFACompressiblePressureSolver.hpp"
+
+#include <dune/common/ErrorMacros.hpp>
+#include <dune/common/SparseTable.hpp>
+#include <dune/porsol/common/LinearSolverISTL.hpp>
+
+
+namespace Dune
+{
+
+
+    template <class GridInterface,
+              class RockInterface,
+              class BCInterface>
+    class TpfaCompressible
+    {
+    public:
+        typedef TPFACompressiblePressureSolver PressureSolver;
+
+        /// @brief
+        ///    Default constructor. Does nothing.
+        TpfaCompressible()
+            : pgrid_(0)
+        {
+        }
+
+
+        /// @brief
+        ///    All-in-one initialization routine.  Enumerates all grid
+        ///    connections, allocates sufficient space, defines the
+        ///    structure of the global system of linear equations for
+        ///    the contact pressures, and computes the permeability
+        ///    dependent inner products for all of the grid's cells.
+        ///
+        /// @param [in] g
+        ///    The grid.
+        ///
+        /// @param [in] r
+        ///    The reservoir properties of each grid cell.
+        ///
+        /// @param [in] grav
+        ///    Gravity vector.  Its Euclidian two-norm value
+        ///    represents the strength of the gravity field (in units
+        ///    of m/s^2) while its direction is the direction of
+        ///    gravity in the current model.
+        ///
+        /// @param [in] bc
+        ///    The boundary conditions describing how the current flow
+        ///    problem interacts with the outside world.  This is used
+        ///    only for the purpose of introducing additional
+        ///    couplings in the case of periodic boundary conditions.
+        ///    The specific values of the boundary conditions are not
+        ///    inspected in @code init() @endcode.
+        template<class Point>
+        void init(const GridInterface&      g,
+                  const RockInterface&      r,
+                  const Point&              grav,
+                  const BCInterface&        bc)
+        {
+            pgrid_ = &g;
+            // Extract perm tensors.
+            const double* perm = &(r.permeability(0)(0,0));
+            std::vector<double> poro(g.numberOfCells(), 1.0);
+            for (int i = 0; i < g.numberOfCells(); ++i) {
+                poro[i] = r.porosity(i);
+            }
+            // Check that we only have noflow boundary conditions.
+            for (int i = 0; i < bc.size(); ++i) {
+                if (bc.flowCond(i).isPeriodic()) {
+                    THROW("Periodic BCs are not handled yet by ifsh.");
+                }
+            }
+            // Initialize 
+            psolver_.init(g.grid(), perm, &poro[0]);
+        }
+
+
+
+
+
+        /// @brief
+        ///    Construct and solve system of linear equations for the
+        ///    pressure values on each interface/contact between
+        ///    neighbouring grid cells.  Recover cell pressure and
+        ///    interface fluxes.  Following a call to @code solve()
+        ///    @encode, you may recover the flow solution from the
+        ///    @code getSolution() @endcode method.
+        ///
+        /// @tparam FluidInterface
+        ///    Type presenting an interface to fluid properties such
+        ///    as density, mobility &c.  The type is expected to
+        ///    provide methods @code phaseMobilities() @endcode and
+        ///    @code phaseDensities() @endcode for phase mobility and
+        ///    density in a single cell, respectively.
+        ///
+        /// @param [in] r
+        ///    The fluid properties of each grid cell.  In method
+        ///    @code solve() @endcode we query this object for the
+        ///    phase mobilities (i.e., @code r.phaseMobilities()
+        ///    @endcode) and the phase densities (i.e., @code
+        ///    phaseDensities() @encode) of each phase.
+        ///
+        /// @param [in] sat
+        ///    Saturation of primary phase.  One scalar value for each
+        ///    grid cell.  This parameter currently limits @code
+        ///    IncompFlowSolverHybrid @endcode to two-phase flow
+        ///    problems.
+        ///
+        /// @param [in] bc
+        ///    The boundary conditions describing how the current flow
+        ///    problem interacts with the outside world.  Method @code
+        ///    solve() @endcode inspects the actual values of the
+        ///    boundary conditions whilst forming the system of linear
+        ///    equations.
+        ///
+        ///    Specifically, the @code bc.flowCond(bid) @endcode
+        ///    method is expected to yield a valid @code FlowBC
+        ///    @endcode object for which the methods @code pressure()
+        ///    @endcode, @code pressureDifference() @endcode, and
+        ///    @code outflux() @endcode yield valid responses
+        ///    depending on the type of the object.
+        ///
+        /// @param [in] src
+        ///    Explicit source terms.  One scalar value for each grid
+        ///    cell representing the rate (in units of m^3/s) of fluid
+        ///    being injected into (>0) or extracted from (<0) a given
+        ///    grid cell.
+        ///
+        /// @param [in] residual_tolerance
+        ///    Control parameter for iterative linear solver software.
+        ///    The iteration process is terminated when the norm of
+        ///    the linear system residual is less than @code
+        ///    residual_tolerance @endcode times the initial residual.
+        ///
+        /// @param [in] linsolver_verbosity
+        ///    Control parameter for iterative linear solver software.
+        ///    Verbosity level 0 prints nothing, level 1 prints
+        ///    summary information, level 2 prints data for each
+        ///    iteration.
+        ///
+        /// @param [in] linsolver_type
+        ///    Control parameter for iterative linear solver software.
+        ///    Type 0 selects a BiCGStab solver, type 1 selects AMG/CG.
+        ///
+        template<class FluidInterface>
+        void solve(const FluidInterface& fl,
+                   const std::vector<double>& cell_pressure,
+                   // const std::vector<typename FluidInterface::ComponentVec>& z,
+                   const std::vector<double>& sat,
+                   const BCInterface& bc,
+                   const std::vector<double>& src,
+                   const double dt,
+                   const double residual_tolerance = 1e-8,
+                   const int linsolver_verbosity = 1,
+                   const int linsolver_type = 1)
+        {
+            // Build bctypes and bcvalues.
+            int num_faces = pgrid_->numberOfFaces();
+            std::vector<PressureSolver::FlowBCTypes> bctypes(num_faces, PressureSolver::FBC_UNSET);
+            std::vector<double> bcvalues(num_faces, 0.0);
+            for (int face = 0; face < num_faces; ++face) {
+                int bid = pgrid_->grid().boundaryId(face);
+                FlowBC face_bc = bc.flowCond(bid);
+                if (face_bc.isDirichlet()) {
+                    bctypes[face] = PressureSolver::FBC_PRESSURE;
+                    bcvalues[face] = face_bc.pressure();
+                } else if (face_bc.isNeumann()) {
+                    bctypes[face] = PressureSolver::FBC_FLUX;
+                    bcvalues[face] = face_bc.outflux(); // TODO: may have to switch sign here depending on orientation.
+                    if (bcvalues[face] != 0.0) {
+                        THROW("Nonzero Neumann conditions not yet properly implemented (signs must be fixed)");
+                    }
+                } else {
+                    THROW("Unhandled boundary condition type.");
+                }
+            }
+
+            // Compute fluid properties.
+//             int num_cells = z.size();
+//             const int np = FluidInterface::numPhases;
+//             const int nc = FluidInterface::numComponents;
+            int num_cells = sat.size();
+            const int np = 3;
+            const int nc = 3;
+            BOOST_STATIC_ASSERT(np == nc);
+            std::vector<double> totcompr(num_cells, 0.0);
+            std::vector<double> voldiscr(num_cells, 0.0);
+            std::vector<double> cellA(num_cells*nc*np, 0.0);
+            std::vector<double> faceA(num_faces*nc*np, 0.0);
+            std::vector<double> phasemobf(num_faces*np, 0.0);
+            std::vector<double> phasemobc(num_cells*np, 0.0); // Just a helper
+            boost::array<double, np> mob;
+            BOOST_STATIC_ASSERT(np == 3);
+            mob[2] = 0.0; // Dummy third phase.
+            // Set cellA to unity for now.
+            for (int cell = 0; cell < num_cells; ++cell) {
+                fl.phaseMobilities(cell, sat[cell], mob);
+                std::copy(mob.begin(), mob.end(), phasemobc.begin() + cell*np);
+                for (int row = 0; row < nc; ++row) {
+                    for (int col = 0; col < np; ++col) {
+                        if (row == col) {
+                            cellA[cell*nc*np + col*nc + row] = 1.0; // Column-wise storage in cellA.
+                        }
+                    }
+                }
+            }
+            // Set cellB to unity for now. Set phasemobf from 'old' interface, to average of cells'.
+            for (int face = 0; face < num_faces; ++face) {
+                int c[2] = { pgrid_->grid().faceCell(face, 0), pgrid_->grid().faceCell(face, 1) };
+                for (int phase = 0; phase < np; ++phase) {
+                    double aver = 0.0;
+                    int num = 0;
+                    for (int j = 0; j < 2; ++j) {
+                        if (c[j] >= 0) {
+                            aver += phasemobc[np*c[j] + phase];
+                            ++num;
+                        }
+                    }
+                    aver /= double(num);
+                    phasemobf[np*face + phase] = aver;
+                }
+                for (int row = 0; row < nc; ++row) {
+                    for (int col = 0; col < np; ++col) {
+                        if (row == col) {
+                            faceA[face*nc*np + col*nc + row] = 1.0; // Column-wise storage in faceA.
+                        }
+                    }
+                }
+            }
+
+            // Assemble system matrix and rhs.
+            psolver_.assemble(src, bctypes, bcvalues, dt, totcompr, voldiscr, cellA, faceA, phasemobf, cell_pressure);
+
+            // Solve system.
+            PressureSolver::LinearSystem s;
+            psolver_.linearSystem(s);
+            parameter::ParameterGroup params;
+            params.insertParameter("linsolver_tolerance", boost::lexical_cast<std::string>(residual_tolerance));
+            params.insertParameter("linsolver_verbosity", boost::lexical_cast<std::string>(linsolver_verbosity));
+            params.insertParameter("linsolver_type", boost::lexical_cast<std::string>(linsolver_type));
+            linsolver_.init(params);
+            LinearSolverISTL::LinearSolverResults res = linsolver_.solve(s.n, s.nnz, s.ia, s.ja, s.sa, s.b, s.x);
+            if (!res.converged) {
+                THROW("Linear solver failed to converge in " << res.iterations << " iterations.\n"
+                      << "Residual reduction achieved is " << res.reduction << '\n');
+            }
+            flow_solution_.clear();
+            std::vector<double> face_flux;
+            psolver_.computePressuresAndFluxes(flow_solution_.pressure_, face_flux);
+            std::vector<double> cell_flux;
+            psolver_.faceFluxToCellFlux(face_flux, cell_flux);
+            const std::vector<int>& ncf = psolver_.numCellFaces();
+            flow_solution_.outflux_.assign(cell_flux.begin(), cell_flux.end(),
+                                           ncf.begin(), ncf.end());
+        }
+
+    private:
+        /// A helper class for postProcessFluxes.
+        class FaceFluxes
+        {
+        public:
+            FaceFluxes(int sz)
+                : fluxes_(sz, 0.0), visited_(sz, 0), max_modification_(0.0)
+            {
+            }
+            void put(double flux, int f_ix) {
+                ASSERT(visited_[f_ix] == 0 || visited_[f_ix] == 1);
+                double sign = visited_[f_ix] ? -1.0 : 1.0;
+                fluxes_[f_ix] += sign*flux;
+                ++visited_[f_ix];
+            }
+            void get(double& flux, int f_ix) {
+                ASSERT(visited_[f_ix] == 0 || visited_[f_ix] == 1);
+                double sign = visited_[f_ix] ? -1.0 : 1.0;
+                double new_flux = 0.5*sign*fluxes_[f_ix];
+                double diff = std::fabs(flux - new_flux);
+                max_modification_ = std::max(max_modification_, diff);
+                flux = new_flux;
+                ++visited_[f_ix];
+            }
+            void resetVisited()
+            {
+                std::fill(visited_.begin(), visited_.end(), 0);
+            }
+
+            double maxMod() const
+            {
+                return max_modification_;
+            }
+        private:
+            std::vector<double> fluxes_;
+            std::vector<int> visited_;
+            double max_modification_;
+
+        };
+
+    public:
+        /// @brief
+        ///    Postprocess the solution fluxes.
+        ///    This method modifies the solution object so that
+        ///    out-fluxes of twin faces (that is, the two faces on a
+        ///    cell-cell intersection) will be made antisymmetric.
+        ///
+        /// @return
+        ///    The maximum modification made to the fluxes.
+        double postProcessFluxes()
+        {
+            typedef typename GridInterface::CellIterator CI;
+            typedef typename CI           ::FaceIterator FI;
+            SparseTable<double>& cflux = flow_solution_.outflux_;
+
+            FaceFluxes face_fluxes(pgrid_->numberOfFaces());
+            // First pass: compute projected fluxes.
+            for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
+                const int cell_index = c->index();
+                for (FI f = c->facebegin(); f != c->faceend(); ++f) {
+                    int f_ix = f->index();
+                    double flux = cflux[cell_index][f->localIndex()];
+                    if (f->boundary()) {
+                        continue;
+                    } else {
+                        face_fluxes.put(flux, f_ix);
+                    }
+                }
+            }
+            face_fluxes.resetVisited();
+            // Second pass: set all fluxes to the projected ones.
+            for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
+                const int cell_index = c->index();
+                for (FI f = c->facebegin(); f != c->faceend(); ++f) {
+                    int f_ix = f->index();
+                    double& flux = cflux[cell_index][f->localIndex()];
+                    if (f->boundary()) {
+                        continue;
+                    } else {
+                        face_fluxes.get(flux, f_ix);
+                    }
+                }
+            }
+            return face_fluxes.maxMod();
+        }
+
+
+
+
+        /// @brief
+        ///    Type representing the solution to a given flow problem.
+        class FlowSolution {
+        public:
+            friend class TpfaCompressible;
+
+            /// @brief
+            ///    The element type of the matrix representation of
+            ///    the mimetic inner product.  Assumed to be a
+            ///    floating point type, and usually, @code Scalar
+            ///    @endcode is an alias for @code double @endcode.
+            typedef typename GridInterface::Scalar       Scalar;
+
+            /// @brief
+            ///    Convenience alias for the grid interface's cell
+            ///    iterator.
+            typedef typename GridInterface::CellIterator CI;
+
+            /// @brief
+            ///    Convenience alias for the cell's face iterator.
+            typedef typename CI           ::FaceIterator FI;
+
+            /// @brief
+            ///    Retrieve the current cell pressure in a given cell.
+            ///
+            /// @param [in] c
+            ///    Cell for which to retrieve the current cell
+            ///    pressure.
+            ///
+            /// @return
+            ///    Current cell pressure in cell @code *c @endcode.
+            Scalar pressure(const CI& c) const
+            {
+                return pressure_[c->index()];
+            }
+
+            /// @brief
+            ///    Retrieve current flux across given face in
+            ///    direction of outward normal vector.
+            ///
+            /// @param [in] f
+            ///    Face across which to retrieve the current outward
+            ///    flux.
+            ///
+            /// @return
+            ///    Current outward flux across face @code *f @endcode.
+            Scalar outflux (const FI& f) const
+            {
+                return outflux_[f->cellIndex()][f->localIndex()];
+            }
+        private:
+            std::vector<Scalar> pressure_;
+            SparseTable<Scalar> outflux_;
+
+            void clear()
+            {
+                pressure_.clear();
+                outflux_.clear();
+            }
+        };
+
+
+
+
+
+        /// @brief
+        ///    Type representing the solution to the problem defined
+        ///    by the parameters to @code solve() @endcode.  Always a
+        ///    reference-to-const.  The @code SolutionType @endcode
+        ///    exposes methods @code pressure(c) @endcode and @code
+        ///    outflux(f) @endcode from which the cell pressure in
+        ///    cell @code *c @endcode and outward-pointing flux across
+        ///    interface @code *f @endcode may be recovered.
+        typedef const FlowSolution& SolutionType;
+
+
+
+
+
+        /// @brief
+        ///    Recover the solution to the problem defined by the
+        ///    parameters to method @code solve() @endcode.  This
+        ///    solution is meaningless without a previous call to
+        ///    method @code solve() @endcode.
+        ///
+        /// @return
+        ///    The current solution.
+        SolutionType getSolution()
+        {
+            return flow_solution_;
+        }
+
+
+
+
+
+    private:
+        const GridInterface* pgrid_;
+        PressureSolver psolver_;
+        LinearSolverISTL linsolver_;
+        FlowSolution flow_solution_;
+    };
+
+
+} // namespace Dune
+
+
+
+#endif // OPM_TPFACOMPRESSIBLE_HEADER_INCLUDED

examples/Makefile.am

@@ -10,6 +10,7 @@ noinst_PROGRAMS = grdecl_to_legacy_test \
                   mimetic_periodic_test \
 		  spe10_test \
 		  tpfa_solver_test \
+		  tpfa_compressible_solver_test \
 		  known_answer_test \
                   aniso_implicitcap_test \
                   implicitcap_test \
@@ -36,6 +37,9 @@ spe10_test_SOURCES = spe10_test.cpp
 tpfa_solver_test_SOURCES = tpfa_solver_test.cpp
 tpfa_solver_test_LDADD = $(LDADD) ../lib/libduneporsol.la
 
+tpfa_compressible_solver_test_SOURCES = tpfa_compressible_solver_test.cpp
+tpfa_compressible_solver_test_LDADD = $(LDADD) ../lib/libduneporsol.la
+
 known_answer_test_SOURCES = known_answer_test.cpp
 
 simulator_test_SOURCES = simulator_test.cpp

examples/tpfa_compressible_solver_test.cpp

@@ -0,0 +1,128 @@
+/*
+  Copyright 2010 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/>.
+*/
+
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <algorithm>
+#include <iostream>
+#include <iomanip>
+
+#include <boost/static_assert.hpp>
+
+#include <dune/common/array.hh>
+#include <dune/common/mpihelper.hh>
+#include <dune/common/Units.hpp>
+
+#include <dune/porsol/common/SimulatorUtilities.hpp>
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/CpGrid.hpp>
+#include <dune/common/EclipseGridParser.hpp>
+#include <dune/common/EclipseGridInspector.hpp>
+
+#include <dune/porsol/common/fortran.hpp>
+#include <dune/porsol/common/blas_lapack.hpp>
+#include <dune/porsol/common/Matrix.hpp>
+#include <dune/porsol/common/GridInterfaceEuler.hpp>
+#include <dune/porsol/common/ReservoirPropertyCapillary.hpp>
+#include <dune/porsol/common/BoundaryConditions.hpp>
+
+#include <dune/porsol/mimetic/TpfaCompressible.hpp>
+#include <dune/common/param/ParameterGroup.hpp>
+#include <dune/porsol/common/setupGridAndProps.hpp>
+
+
+template<int dim, class GI, class RI>
+void test_flowsolver(const GI& g, const RI& r)
+{
+    typedef typename GI::CellIterator                   CI;
+    typedef typename CI::FaceIterator                   FI;
+    typedef Dune::BasicBoundaryConditions<true, false>  FBC;
+    typedef Dune::TpfaCompressible<GI, RI, FBC> FlowSolver;
+    // typedef Dune::TpfaInterface<GI, RI, FBC> FlowSolver;
+    // typedef Dune::IfshInterface<GI, RI, FBC> FlowSolver;
+    // typedef Dune::IncompFlowSolverHybrid<GI, RI, FBC,
+    //                                   Dune::MimeticIPEvaluator> FlowSolver;
+
+    FlowSolver solver;
+
+    typedef Dune::FlowBC BC;
+    FBC flow_bc(7);
+
+//     flow_bc.flowCond(5) = BC(BC::Dirichlet, 100.0*Dune::unit::barsa);
+//     flow_bc.flowCond(6) = BC(BC::Dirichlet, 0.0*Dune::unit::barsa);
+
+    typename CI::Vector gravity(0.0);
+//     gravity[2] = Dune::unit::gravity;
+
+    solver.init(g, r, gravity, flow_bc);
+
+    std::vector<double> src(g.numberOfCells(), 0.0);
+    std::vector<double> sat(g.numberOfCells(), 0.0);
+    if (g.numberOfCells() > 1) {
+        src[0]     = 1.0;
+        src.back() = -1.0;
+    }
+    std::vector<double> cell_pressure(g.numberOfCells(), 0.0);
+    double dt = 1;
+
+    solver.solve(r, cell_pressure, sat, flow_bc, src, dt, 1e-8, 3, 1);
+
+
+    typedef typename FlowSolver::SolutionType FlowSolution;
+    FlowSolution soln = solver.getSolution();
+
+    std::vector<typename GI::Vector> cell_velocity;
+    estimateCellVelocity(cell_velocity, g, soln);
+    // Dune's vtk writer wants multi-component data to be flattened.
+    std::vector<double> cell_velocity_flat(&*cell_velocity.front().begin(),
+                                           &*cell_velocity.back().end());
+    getCellPressure(cell_pressure, g, soln);
+
+    Dune::VTKWriter<typename GI::GridType::LeafGridView> vtkwriter(g.grid().leafView());
+    vtkwriter.addCellData(cell_velocity_flat, "velocity", dim);
+    vtkwriter.addCellData(cell_pressure, "pressure");
+    vtkwriter.write("testsolution-" + boost::lexical_cast<std::string>(0),
+                    Dune::VTKOptions::ascii);
+}
+
+
+using namespace Dune;
+
+int main(int argc, char** argv)
+{
+    Dune::parameter::ParameterGroup param(argc, argv);
+    Dune::MPIHelper::instance(argc,argv);
+
+    // Make a grid and props.
+    Dune::CpGrid grid;
+    ReservoirPropertyCapillary<3> res_prop;
+    setupGridAndProps(param, grid, res_prop);
+
+    // Make the grid interface.
+    Dune::GridInterfaceEuler<Dune::CpGrid> g(grid);
+
+    // Run test.
+    test_flowsolver<3>(g, res_prop);
+}
+