diff --git a/.hgsubstate b/.hgsubstate index 3411215..2ec64ad 100644 --- a/.hgsubstate +++ b/.hgsubstate @@ -1 +1 @@ -00c15b54de2e838a2cb85ff20cb9d4b8edde1f34 dune/porsol/opmpressure +335bb8e27ad82ccf4b60c2a0042041a704d50691 dune/porsol/opmpressure diff --git a/dune/porsol/mimetic/TpfaCompressible.hpp b/dune/porsol/mimetic/TpfaCompressible.hpp new file mode 100644 index 0000000..5d1590e --- /dev/null +++ b/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 . +*/ + +#ifndef OPM_TPFACOMPRESSIBLE_HEADER_INCLUDED +#define OPM_TPFACOMPRESSIBLE_HEADER_INCLUDED + + +#include "../opmpressure/src/TPFACompressiblePressureSolver.hpp" + +#include +#include +#include + + +namespace Dune +{ + + + template + 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 + 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 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 + void solve(const FluidInterface& fl, + const std::vector& cell_pressure, + // const std::vector& z, + const std::vector& sat, + const BCInterface& bc, + const std::vector& 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 bctypes(num_faces, PressureSolver::FBC_UNSET); + std::vector 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 totcompr(num_cells, 0.0); + std::vector voldiscr(num_cells, 0.0); + std::vector cellA(num_cells*nc*np, 0.0); + std::vector faceA(num_faces*nc*np, 0.0); + std::vector phasemobf(num_faces*np, 0.0); + std::vector phasemobc(num_cells*np, 0.0); // Just a helper + boost::array 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(residual_tolerance)); + params.insertParameter("linsolver_verbosity", boost::lexical_cast(linsolver_verbosity)); + params.insertParameter("linsolver_type", boost::lexical_cast(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 face_flux; + psolver_.computePressuresAndFluxes(flow_solution_.pressure_, face_flux); + std::vector cell_flux; + psolver_.faceFluxToCellFlux(face_flux, cell_flux); + const std::vector& 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 fluxes_; + std::vector 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& 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 pressure_; + SparseTable 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 diff --git a/examples/Makefile.am b/examples/Makefile.am index 19fb4f0..b75b317 100644 --- a/examples/Makefile.am +++ b/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 diff --git a/examples/tpfa_compressible_solver_test.cpp b/examples/tpfa_compressible_solver_test.cpp new file mode 100644 index 0000000..a79df1b --- /dev/null +++ b/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 . +*/ + + +#if HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + +template +void test_flowsolver(const GI& g, const RI& r) +{ + typedef typename GI::CellIterator CI; + typedef typename CI::FaceIterator FI; + typedef Dune::BasicBoundaryConditions FBC; + typedef Dune::TpfaCompressible FlowSolver; + // typedef Dune::TpfaInterface FlowSolver; + // typedef Dune::IfshInterface FlowSolver; + // typedef Dune::IncompFlowSolverHybrid 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 src(g.numberOfCells(), 0.0); + std::vector sat(g.numberOfCells(), 0.0); + if (g.numberOfCells() > 1) { + src[0] = 1.0; + src.back() = -1.0; + } + std::vector 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 cell_velocity; + estimateCellVelocity(cell_velocity, g, soln); + // Dune's vtk writer wants multi-component data to be flattened. + std::vector cell_velocity_flat(&*cell_velocity.front().begin(), + &*cell_velocity.back().end()); + getCellPressure(cell_pressure, g, soln); + + Dune::VTKWriter vtkwriter(g.grid().leafView()); + vtkwriter.addCellData(cell_velocity_flat, "velocity", dim); + vtkwriter.addCellData(cell_pressure, "pressure"); + vtkwriter.write("testsolution-" + boost::lexical_cast(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 g(grid); + + // Run test. + test_flowsolver<3>(g, res_prop); +} +