From 6ee73daad8e081412612553820d5f008713c3811 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 18 Nov 2010 16:02:48 +0100 Subject: [PATCH] Initial attempt at C++ interface for compressible TPFA solver. --- src/TPFACompressiblePressureSolver.hpp | 301 +++++++++++++++++++++++++ 1 file changed, 301 insertions(+) create mode 100644 src/TPFACompressiblePressureSolver.hpp diff --git a/src/TPFACompressiblePressureSolver.hpp b/src/TPFACompressiblePressureSolver.hpp new file mode 100644 index 00000000..75b75805 --- /dev/null +++ b/src/TPFACompressiblePressureSolver.hpp @@ -0,0 +1,301 @@ +/* + 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_TPFACOMPRESSIBLEPRESSURESOLVER_HEADER_INCLUDED +#define OPM_TPFACOMPRESSIBLEPRESSURESOLVER_HEADER_INCLUDED + + +#include "cfs_tpfa.h" +#include "trans_tpfa.h" +#include "sparse_sys.h" +#include "flow_bc.h" +#include "compr_quant.h" +#include "GridAdapter.hpp" +#include + + + +/// @brief +/// Encapsulates the cfs_tpfa (= compressible flow solver +/// two-point flux approximation) solver modules. +class TPFACompressiblePressureSolver +{ +public: + /// @brief + /// Default constructor, does nothing. + TPFACompressiblePressureSolver() + : state_(Uninitialized), data_(0) + { + } + + /// @brief + /// Destructor. + ~TPFACompressiblePressureSolver() + { + cfs_tpfa_destroy(data_); + } + + /// @brief + /// Initialize the solver's structures for a given grid (at some point also well pattern). + /// @tparam Grid This must conform to the SimpleGrid concept. + /// @param grid The grid object. + /// @param perm Permeability. It should contain dim*dim entries (a full tensor) for each cell. + /// @param gravity Array containing gravity acceleration vector. It should contain dim entries. + template + void init(const Grid& grid, const double* perm, const double* porosity) + { + // Build C grid structure. + grid_.init(grid); + + // Build (empty for now) C well structure. + // well_t* w = 0; + + // Initialize data. + int num_phases = 3; + data_ = cfs_tpfa_construct(grid_.c_grid(), num_phases); + if (!data_) { + throw std::runtime_error("Failed to initialize cfs_tpfa solver."); + } + + // Compute half-transmissibilities + int num_cells = grid.numCells(); + int ngconn = grid_.c_grid()->cell_facepos[num_cells]; + ncf_.resize(num_cells); + int count = 0; + for (int cell = 0; cell < num_cells; ++cell) { + int num_local_faces = grid.numCellFaces(cell); + ncf_[cell] = num_local_faces; + } + assert(count == ngconn); + htrans_.resize(ngconn); + tpfa_htrans_compute(grid_.c_grid(), perm, &htrans_[0]); + + // Compute transmissibilities. + trans_.resize(grid_.numFaces()); + tpfa_trans_compute(grid_.c_grid(), &htrans_[0], &trans_[0]); + + // Compute pore volumes. + porevol_.resize(num_cells); + for (int i = 0; i < num_cells; ++i) { + porevol_[i] = porosity[i]*grid.cellVolume(i); + } + + state_ = Initialized; + } + + + enum FlowBCTypes { FBC_UNSET, FBC_PRESSURE, FBC_FLUX }; + + /// @brief + /// Assemble the sparse system. + /// You must call init() prior to calling assemble(). + /// @param sources Source terms, one per cell. Positive numbers + /// are sources, negative are sinks. + /// @param total_mobilities Scalar total mobilities, one per cell. + /// @param omegas Gravity term, one per cell. In a multi-phase + /// flow setting this is equal to + /// \f[ \omega = \sum_{p} \frac{\lambda_p}{\lambda_t} \rho_p \f] + /// where \f$\lambda_p\f$ is a phase mobility, \f$\rho_p\f$ is a + /// phase density and \f$\lambda_t\f$ is the total mobility. + void assemble(const std::vector& sources, + const std::vector& bctypes, + const std::vector& bcvalues, + const double dt, + const std::vector& totcompr, + const std::vector& voldiscr, + const std::vector& cellA, // num phases^2 * num cells, fortran ordering! + const std::vector& faceA, // num phases^2 * num faces, fortran ordering! + const std::vector& phasemobf, + const std::vector& cell_pressure) + { + if (state_ == Uninitialized) { + throw std::runtime_error("Error in TPFACompressiblePressureSolver::assemble(): You must call init() prior to calling assemble()."); + } + grid_t* g = grid_.c_grid(); + + // Boundary conditions. + int num_faces = g->number_of_faces; + assert(num_faces == int(bctypes.size())); + bctypes_.clear(); + bctypes_.resize(num_faces, UNSET); + for (int face = 0; face < num_faces; ++face) { + if (bctypes[face] != FBC_FLUX || bcvalues[face] != 0.0) { + throw std::logic_error("TPFACompressiblePressureSolver currently only supports noflow bcs."); + } + if (bctypes[face] == FBC_PRESSURE) { + bctypes_[face] = PRESSURE; + } else if (bctypes[face] == FBC_FLUX) { + bctypes_[face] = FLUX; + } + } + bcvalues_ = bcvalues; + flowbc_t bc = { &bctypes_[0], const_cast(&bcvalues_[0]) }; + + // Source terms from user. + double* src = const_cast(&sources[0]); // Ugly? Yes. Safe? I think so. + + // All well related things are zero. +// well_control_t* wctrl = 0; +// double* WI = 0; +// double* wdp = 0; + + // Assemble the embedded linear system. + compr_quantities cq = { 3, &totcompr[0], &voldiscr[0], &cellA[0], &faceA[0], &phasemobf[0] }; + std::vector gravcap_f(3*num_faces, 0.0); + cfs_tpfa_assemble(g, dt, &bc, src, + &cq, &trans_[0], &gravcap_f[0], &cell_pressure[0], &porevol_[0], + data_); + phasemobf_ = phasemobf; + state_ = Assembled; + } + + /// Encapsulate a sparse linear system in CSR format. + struct LinearSystem + { + int n; + int nnz; + int* ia; + int* ja; + double* sa; + double* b; + double* x; + }; + + /// @brief + /// Access the linear system assembled. + /// You must call assemble() prior to calling linearSystem(). + /// @param[out] s The linear system encapsulation to modify. + /// After this call, s will point to linear system structures + /// that are owned and allocated internally. + void linearSystem(LinearSystem& s) + + { + if (state_ != Assembled) { + throw std::runtime_error("Error in TPFACompressiblePressureSolver::linearSystem(): " + "You must call assemble() prior to calling linearSystem()."); + } + s.n = data_->A->n; + s.nnz = data_->A->nnz; + s.ia = data_->A->ia; + s.ja = data_->A->ja; + s.sa = data_->A->sa; + s.b = data_->b; + s.x = data_->x; + } + + /// @brief + /// Compute cell pressures and face fluxes. + /// You must call assemble() (and solve the linear system accessed + /// by calling linearSystem()) prior to calling + /// computePressuresAndFluxes(). + /// @param[out] cell_pressures Cell pressure values. + /// @param[out] face_areas Face flux values. + void computePressuresAndFluxes(std::vector& cell_pressures, + std::vector& face_fluxes) + { + if (state_ != Assembled) { + throw std::runtime_error("Error in TPFACompressiblePressureSolver::computePressuresAndFluxes(): " + "You must call assemble() (and solve the linear system) " + "prior to calling computePressuresAndFluxes()."); + } + int num_cells = grid_.c_grid()->number_of_cells; + int num_faces = grid_.c_grid()->number_of_faces; + cell_pressures.clear(); + cell_pressures.resize(num_cells, 0.0); + face_fluxes.clear(); + face_fluxes.resize(num_faces, 0.0); +// ifs_tpfa_press_flux(grid_.c_grid(), &eff_trans_[0], +// data_, &cell_pressures[0], &face_fluxes[0]); + flowbc_t bc = { &bctypes_[0], const_cast(&bcvalues_[0]) }; + int np = 3; // Number of phases. + cfs_tpfa_press_flux(grid_.c_grid(), + &bc, np, &trans_[0], &htrans_[0], &phasemobf_[0], + data_, &cell_pressures[0], &face_fluxes[0]); + } + + /// @brief + /// Compute cell fluxes from face fluxes. + /// You must call assemble() (and solve the linear system accessed + /// by calling linearSystem()) prior to calling + /// faceFluxToCellFlux(). + /// @param face_fluxes + /// @param face_areas Face flux values (usually output from computePressuresAndFluxes()). + /// @param[out] cell_fluxes Cell-wise flux values. + /// They are given in cell order, and for each cell there is + /// one value for each adjacent face (in the same order as the + /// cell-face topology of the grid). Positive values represent + /// fluxes out of the cell. + void faceFluxToCellFlux(const std::vector& face_fluxes, + std::vector& cell_fluxes) + { + if (state_ != Assembled) { + throw std::runtime_error("Error in TPFACompressiblePressureSolver::faceFluxToCellFlux(): " + "You must call assemble() (and solve the linear system) " + "prior to calling faceFluxToCellFlux()."); + } + const grid_t& g = *(grid_.c_grid()); + int num_cells = g.number_of_cells; + cell_fluxes.resize(g.cell_facepos[num_cells]); + for (int cell = 0; cell < num_cells; ++cell) { + for (int hface = g.cell_facepos[cell]; hface < g.cell_facepos[cell + 1]; ++hface) { + int face = g.cell_faces[hface]; + bool pos = (g.face_cells[2*face] == cell); + cell_fluxes[hface] = pos ? face_fluxes[face] : -face_fluxes[face]; + } + } + } + + /// @brief + /// Access the number of connections (faces) per cell. Deprecated, will be removed. + const std::vector& numCellFaces() + { + return ncf_; + } + +private: + // Disabling copy and assigment for now. + TPFACompressiblePressureSolver(const TPFACompressiblePressureSolver&); + TPFACompressiblePressureSolver& operator=(const TPFACompressiblePressureSolver&); + + enum State { Uninitialized, Initialized, Assembled }; + State state_; + + // Solver data. + cfs_tpfa_data* data_; + // Grid. + GridAdapter grid_; + // Number of faces per cell. + std::vector ncf_; + // Transmissibility storage. + std::vector htrans_; + std::vector trans_; + // Pore volumes. + std::vector porevol_; + // Phase mobilities per face. + std::vector phasemobf_; + + // Boundary conditions. + std::vector bctypes_; + std::vector bcvalues_; +}; + + + + +#endif // OPM_TPFACOMPRESSIBLEPRESSURESOLVER_HEADER_INCLUDED