Merge pull request #4 from bska/master

Document everal low-level OPM-Core interfaces related to mimetic solvers
This commit is contained in:
Alf Birger Rustad 2012-07-05 07:02:03 -07:00
commit 819061fcc2
14 changed files with 580 additions and 87 deletions

View File

@ -6,7 +6,7 @@ SUBDIRS = . tests examples tutorials
# ----------------------------------------------------------------------
# Declare products (i.e., the library)
lib_LTLIBRARIES = libopmcore.la
lib_LTLIBRARIES = lib/libopmcore.la
# ----------------------------------------------------------------------
# Build-time flags needed to build libopmcore.la
@ -18,7 +18,7 @@ $(BOOST_CPPFLAGS)
# Link-time flags needed both to successfully link the library and to
# (transitively) convey inter-library dependency information.
libopmcore_la_LDFLAGS = \
lib_libopmcore_la_LDFLAGS = \
$(BOOST_LDFLAGS) \
$(BOOST_FILESYSTEM_LIB) \
$(BOOST_SYSTEM_LIB) \
@ -31,7 +31,11 @@ $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
#
# Please try to keep the list sorted.
libopmcore_la_SOURCES = \
# List of sources that should be built but not distributed. See AGMG
# support below for additional details.
nodist_lib_libopmcore_la_SOURCES =
lib_libopmcore_la_SOURCES = \
opm/core/GridManager.cpp \
opm/core/eclipse/EclipseGridInspector.cpp \
opm/core/eclipse/EclipseGridParser.cpp \
@ -256,7 +260,7 @@ opm/core/wells/WellsManager.hpp
# Optional library constituents.
if UMFPACK
libopmcore_la_SOURCES += \
lib_libopmcore_la_SOURCES += \
opm/core/linalg/call_umfpack.c \
opm/core/linalg/LinearSolverUmfpack.cpp
@ -267,7 +271,7 @@ endif
if DUNE_ISTL
libopmcore_la_SOURCES += \
lib_libopmcore_la_SOURCES += \
opm/core/linalg/LinearSolverIstl.cpp
nobase_include_HEADERS += \
@ -276,14 +280,16 @@ endif
if BUILD_AGMG
libopmcore_la_SOURCES += \
nodist_lib_libopmcore_la_SOURCES += \
$(AGMG_SRCDIR)/dagmg.f90 \
$(AGMG_SRCDIR)/dagmg_mumps.f90 \
$(AGMG_SRCDIR)/dagmg_mumps.f90
lib_libopmcore_la_SOURCES += \
opm/core/linalg/LinearSolverAGMG.cpp
nobase_include_HEADERS += \
opm/core/linalg/LinearSolverAGMG.hpp
libopmcore_la_LDFLAGS += \
lib_libopmcore_la_LDFLAGS += \
$(FCLIBS)
endif

View File

@ -5,7 +5,7 @@ $(BOOST_CPPFLAGS)
# All targets link to the library
LDADD = \
$(top_builddir)/libopmcore.la \
$(top_builddir)/lib/libopmcore.la \
$(BOOST_FILESYSTEM_LIB) \
$(BOOST_SYSTEM_LIB)

View File

@ -41,8 +41,8 @@ AC_DEFUN([OPM_DYNLINK_BOOST_TEST],
AC_REQUIRE([OPM_BOOST_BASE])
AC_REQUIRE([AX_BOOST_UNIT_TEST_FRAMEWORK])
_opm_LIBS_SAVE="$LIBS"
_opm_CPPFLAGS_SAVE="$CPPFLAGS"
_opm_LIBS_SAVE="${LIBS}"
_opm_CPPFLAGS_SAVE="${CPPFLAGS}"
LIBS="${BOOST_LDFLAGS} ${BOOST_UNIT_TEST_FRAMEWORK_LIB} ${LIBS}"
CPPFLAGS="${BOOST_CPPFLAGS} ${CPPFLAGS}"
@ -66,12 +66,12 @@ AC_CACHE_CHECK([if the Boost.Test library can be linked dynamically],dnl
AC_LANG_POP([C++])
LIBS="$_opm_LIBS_SAVE"
CPPFLAGS="$_opm_CPPFLAGS_SAVE"
LIBS="${_opm_LIBS_SAVE}"
CPPFLAGS="${_opm_CPPFLAGS_SAVE}"
AS_IF([test x"$opm_cv_boost_link_static" = x"yes" -o \
x"$opm_cv_boost_link_dynamic" = x"yes"],
[AS_IF([test x"$opm_cv_boost_link_dynamic" = x"yes"],
AS_IF([test x"${opm_cv_boost_link_static}" = x"yes" -o \
x"${opm_cv_boost_link_dynamic}" = x"yes"],
[AS_IF([test x"${opm_cv_boost_link_dynamic}" = x"yes"],
[AC_DEFINE([HAVE_DYNAMIC_BOOST_TEST], [1],
[Define to `1' if Boost.Test should use BOOST_TEST_DYN_LINK])]
[:])[]dnl

View File

@ -780,9 +780,11 @@ struct WelspecsLine
int fluids_in_place_reg_numb_; // Fluids in place region number
WelspecsLine() :
datum_depth_BHP_(-1.0), drain_rad_(0.0), spec_inflow_("STD"),
shut_in_("SHUT"), crossflow_("YES"), pressure_table_number_(0),
density_calc_type_("SEG"), fluids_in_place_reg_numb_(0)
name_(), group_(), I_(-1), J_(-1),
datum_depth_BHP_(-1.0), pref_phase_(), drain_rad_(0.0),
spec_inflow_("STD"), shut_in_("SHUT"), crossflow_("YES"),
pressure_table_number_(0), density_calc_type_("SEG"),
fluids_in_place_reg_numb_(0)
{}
};
@ -1360,8 +1362,8 @@ struct WgrupconLine
bool available_for_group_control_;
double guide_rate_;
std::string phase_;
WgrupconLine() :
available_for_group_control_(true)
WgrupconLine()
: well_(), available_for_group_control_(true), guide_rate_(-1.0), phase_()
{
}
};
@ -1580,6 +1582,7 @@ struct WeltargLine
double new_value_; // New value of this quantity
WeltargLine()
: well_(), control_change_(), new_value_(-1.0)
{
}
};
@ -1739,6 +1742,13 @@ struct EquilLine
// initial fluids in place calculation.
EquilLine()
{
datum_depth_ = datum_depth_pressure_ = 0.0;
water_oil_contact_depth_ = oil_water_cap_pressure_ = 0.0;
gas_oil_contact_depth_ = gas_oil_cap_pressure_ = 0.0;
live_oil_table_index_ = 0;
wet_gas_table_index_ = 0;
N_ = 0;
}
};
@ -2128,7 +2138,8 @@ struct WpolymerLine
WpolymerLine()
{
well_ = polymer_group_ = salt_group_ = "";
well_ = polymer_group_ = salt_group_ = "";
polymer_concentration_ = salt_concentration_ = 0.0;
}
};

View File

@ -99,6 +99,7 @@ namespace Opm
/// Default constructor.
SaturationPropsBasic::SaturationPropsBasic()
: num_phases_(0), relperm_func_(Constant)
{
}

View File

@ -27,8 +27,9 @@
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/pressure/mimetic/mimetic.h> // for updating gpress
#include <opm/core/GridAdapter.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <stdexcept>
#include <cassert>
/// @brief

View File

@ -358,28 +358,193 @@ hybsys_schur_comp_gen(int nc, const int *pconn,
const double *Binv, const double *C2,
const double *P, struct hybsys *sys);
/**
* Compute elemental contributions to global, symmetric system of
* simultaneous linear equations from cell<->well connections.
*
* Specifically, for a well @c w intersecting a cell @c c, this function
* computes the elemental contributions
* \f[
* (F_1)_{wc} = C_{wc}^\mathsf{T} B_{wc}^{-1} D_{wc} = \mathit{WI}_{wc}
* \f]
* and
* \f[
* L_{wc} = C_{wc}^\mathsf{T} B_{wc}^{-1} C_{wc} = \mathit{WI}_{wc}
* \f]
* and incorporates the contributions into the global system quantities
* as appropriate.
*
* This function modifies <CODE>sys->L</CODE> and <CODE>wsys->F1</CODE>.
*
* @param[in] nc Total number of grid cells.
* @param[in] cwpos Indirection array that defines each cell's
* connecting wells. Values typically computed
* using function derive_cell_wells().
* @param[in] WI Peaceman well connection indices. Array of
* size <CODE>cwpos[nc]</CODE>. Must incorporate
* effects of multiple phases (i.e., total mobility)
* if applicable.
* @param[in,out] sys Hybrid system management structure allocated
* using hybsys_allocate_symm() and initialised
* using hybsys_init() and/or filled using function
* hybsys_schur_comp_symm().
* @param[in,out] wsys Hybrid well-system management structure obtained
* from function hybsys_well_allocate_symm().
*/
void
hybsys_well_schur_comp_symm(int nc, const int *cwpos,
double *WI,
struct hybsys *sys,
struct hybsys_well *wsys);
/**
* Compute final (symmetric) Schur complement contributions to
* global system of simultaneous linear equations.
*
* This function forms the coefficient matrix
* \f[
* S_c = D^\mathsf{T}B_c^{-1}D - F_c^\mathsf{T}L_c^{-1}F_c
* \f]
* and similar right-hand side \f$r_c\f$ elemental contributions.
* These values must be subsequently assembled into the global system
* using function hybsys_global_assemble_cell() after imposing any
* applicable boundary conditions.
*
* This function overwrites the fields @c S and @c r of the hybrid system
* structure.
*
* @param[in] c Cell for which to compute local contributions.
* @param[in] nconn Number of connections (faces) of cell @c c.
* @param[in] p1 Start address (into @c gpress) of the gravity
* contributions of cell @c c.
* @param[in] p2 Start address (into @c Binv) of the inverse
* inner product of cell @c c.
* @param[in] gpress Gravity contributions of all cells. Must
* include effects of multiple phases if applicable.
* @param[in] src Explicit source terms for all cells.
* @param[in] Binv Inverse inner products for all cells. Must
* include effects of multiple phases if applicable.
* @param[in,out] sys Hybrid system management structure allocated
* using hybsys_allocate_symm() and initialised
* using hybsys_init() and/or filled using function
* hybsys_schur_comp_symm() and
* hybsys_well_schur_comp_symm() if applicable.
*/
void
hybsys_cellcontrib_symm(int c, int nconn, int p1, int p2,
const double *gpress, const double *src,
const double *Binv, struct hybsys *sys);
/**
* Compute final (non-symmetric) Schur complement contributions to
* global system of simultaneous linear equations.
*
* This function forms the coefficient matrix
* \f[
* S_c = D^\mathsf{T}B_c^{-1}D - (F_1)_c^\mathsf{T}L_c^{-1}(F_2)_c
* \f]
* and similar right-hand side \f$r_c\f$ elemental contributions.
* These values must be subsequently assembled into the global system
* using function hybsys_global_assemble_cell() after imposing any
* applicable boundary conditions.
*
* This function overwrites the fields @c S and @c r of the hybrid system
* structure.
*
* @param[in] c Cell for which to compute local contributions.
* @param[in] nconn Number of connections (faces) of cell @c c.
* @param[in] p1 Start address (into @c gpress) of the gravity
* contributions of cell @c c.
* @param[in] p2 Start address (into @c Binv) of the inverse
* inner product of cell @c c.
* @param[in] gpress Gravity contributions of all cells. Must
* include effects of multiple phases if applicable.
* @param[in] src Explicit source terms for all cells.
* @param[in] Binv Inverse inner products for all cells. Must
* include effects of multiple phases if applicable.
* @param[in,out] sys Hybrid system management structure allocated
* using hybsys_allocate_symm() and initialised
* using hybsys_init() and/or filled using functions
* hybsys_schur_comp_unsymm() or hybsys_schur_comp_gen().
*/
void
hybsys_cellcontrib_unsymm(int c, int nconn, int p1, int p2,
const double *gpress, const double *src,
const double *Binv, struct hybsys *sys);
/**
* Form elemental direct contributions to global system of simultaneous linear
* equations from cell<->well interactions.
*
* Plays a role similar to function hybsys_cellcontrib_symm(), but for wells.
*
* @param[in] c Cell for which to compute cell<->well Schur complement
* @param[in] ngconn Number of inter-cell connections (faces) of cell @c c.
* @param[in] p1 Start index (into <CODE>sys->F1</CODE>) of cell @c c.
* @param[in] cwpos Indirection array that defines each cell's connecting
* wells. Must coincide with equally named parameter of
* function hybsys_well_schur_comp_symm().
* @param[in] WI Peaceman well connection indices. Array of
* size <CODE>pwconn[nc]</CODE>. Must coincide with
* equally named parameter of contribution function
* hybsys_well_schur_comp_symm().
* @param[in] wdp Well connection gravity pressure adjustments.
* One scalar for each well connection in an array of size
* <CODE>pwconn[nc]</CODE>.
* @param[in,out] sys Hybrid system management structure filled using
* functions hybsys_schur_comp_unsymm() or
* hybsys_schur_comp_gen().
* @param[in,out] wsys Hybrid well-system management structure filled using
* function hybsys_well_schur_comp_symm().
*/
void
hybsys_well_cellcontrib_symm(int c, int ngconn, int p1,
const int *cwpos,
const double *WI, const double *wdp,
struct hybsys *sys, struct hybsys_well *wsys);
/**
* Recover cell pressures and outward fluxes (with respect to cells--i.e., the
* ``half-face fluxes'') through back substitution after solving a symmetric
* (i.e., incompressible) Schur complement system of simultaneous linear
* equations.
*
* Specifically, given the solution \f$\pi\f$ to the global system of
* simultaneous linear equations, \f$A\pi=b\f$, that arises as a result of the
* Schur complement analysis, this function recovers the cell pressures \f$p\f$
* and outward fluxes \f$v\f$ defined by
* \f[
* \begin{aligned}
* Lp &= g - C_2^\mathsf{T}B^{-1}G + F_2\pi \\
* Bv &= G + C_1p - D\pi
* \end{aligned}.
* \f]
*
* @param[in] nc Total number of grid cells.
* @param[in] pconn Cell-to-face start pointers.
* @param[in] conn Cell-to-face mapping.
* @param[in] gpress Gravity contributions of all cells. Must coincide with
* equally named parameter in calls to cell contribution
* functions such as hybsys_cellcontrib_symm().
* @param[in] Binv Inverse inner products for all cells. Must coincide
* with equally named parameter in calls to contribution
* functions such as hybsys_cellcontrib_symm().
* @param[in] sys Hybrid system management structure coinciding with
* equally named parameter in contribution functions such
* as hybsys_cellcontrib_symm() or
* hybsys_cellcontrib_unsymm().
* @param[in] pi Solution (interface/contact pressure) obtained from
* solving the global system \f$A\pi = b\f$.
* @param[out] press Cell pressures, \f$p\f$. Array of size @c nc.
* @param[out] flux Outward interface fluxes, \f$v\f$. Array of size
* <CODE>pconn[nc]</CODE>.
* @param[in,out] work Scratch array for temporary results. Array of size at
* least \f$\max_c \{ \mathit{pconn}_{c + 1}
* - \mathit{pconn}_c \} \f$.
*/
void
hybsys_compute_press_flux(int nc, const int *pconn, const int *conn,
const double *gpress,
@ -387,6 +552,52 @@ hybsys_compute_press_flux(int nc, const int *pconn, const int *conn,
const double *pi, double *press, double *flux,
double *work);
/**
* Recover well pressures (i.e., bottom-hole pressure values) and well
* connection (perforation) fluxes.
*
* Specifically, this function performs the same role (i.e., back-substitution)
* for wells as function hybsys_compute_press_flux() does for grid cells and
* grid contacts (interfaces).
*
* @param[in] nc Total number of grid cells.
* @param[in] pgconn Cell-to-face start pointers.
* @param[in] nf Total number of grid faces.
* @param[in] nw Total number of wells.
* @param[in] pwconn Cell-to-well start pointers. If <CODE>nw > 0</CODE>,
* then this parameter must coincide with the @c cwpos
* array used in call to hybsys_well_schur_comp_symm().
* @param[in] wconn Cell-to-well mapping.
* @param[in] Binv Inverse inner products for all cells. Must coincide
* with equally named parameter in calls to contribution
* functions such as hybsys_well_cellcontrib_symm().
* @param[in] WI Peaceman well connection indices. Array of
* size <CODE>pwconn[nc]</CODE>. Must coincide with
* equally named parameter of contribution function
* hybsys_well_cellcontrib_symm().
* @param[in] wdp Well connection gravity pressure adjustments.
* @param[in] sys Hybrid system management structure coinciding with
* equally named parameter in contribution functions such
* as hybsys_cellcontrib_symm() and
* hybsys_well_cellcontrib_symm().
* @param[in] wsys Hybrid well-system management structure. Must coincide
* with equally named paramter of contribution function
* hybsys_well_cellcontrib_symm().
* @param[in] pi Solution (interface/contact pressure and well BHPs)
* obtained from solving the global system \f$A\pi = b\f$.
* @param[in] cpress Cell pressures, \f$p\f$, obtained from a previous call
* to function hybsys_compute_press_flux().
* @param[in] cflux Outward fluxes, \f$v\f$, obtained from a previous call
* to function hybsys_compute_press_flux().
* @param[out] wpress Well (i.e., bottom-hole) pressures. Array of size
* @c nw.
* @param[out] wflux Well connection (perforation) fluxes. Array of size
* <CODE>pwconn[nw]</CODE>.
* @param[in,out] work Scratch array for storing intermediate results. Array
* of size at least \f$\max_w \{ \mathit{pwconn}_{w + 1}
* - \mathit{pwconn}_w\}\f$.
*/
void
hybsys_compute_press_flux_well(int nc, const int *pgconn, int nf,
int nw, const int *pwconn, const int *wconn,

View File

@ -20,6 +20,16 @@
#ifndef OPM_IFS_TPFA_HEADER_INCLUDED
#define OPM_IFS_TPFA_HEADER_INCLUDED
/**
* \file
* Interfaces and data structures to assemble a system of simultaneous linear
* equations discretising a flow problem that is either incompressible or
* features rock compressibility using the two-point flux approximation method.
*
* Includes support for reconstructing the Darcy flux field as well as well
* connection fluxes.
*/
#include <opm/core/grid.h>
#ifdef __cplusplus
@ -31,37 +41,66 @@ struct CSRMatrix;
struct FlowBoundaryConditions;
struct Wells;
/**
* Main data structure presenting a view of an assembled system of simultaneous
* linear equations which may be solved using external software.
*/
struct ifs_tpfa_data {
struct CSRMatrix *A;
double *b;
double *x;
struct CSRMatrix *A; /**< Coefficient matrix */
double *b; /**< Right-hand side */
double *x; /**< Solution */
struct ifs_tpfa_impl *pimpl;
struct ifs_tpfa_impl *pimpl; /**< Internal management structure */
};
/**
* Solution variables.
*/
struct ifs_tpfa_solution {
double *cell_press;
double *face_flux ;
double *cell_press; /**< Cell pressures */
double *face_flux ; /**< Interface fluxes */
double *well_press; /* BHP */
double *well_flux ; /* Perforation (total) fluxes */
double *well_press; /**< Bottom-hole pressures for each well */
double *well_flux ; /**< Well connection total fluxes */
};
/**
* Driving forces pertaining to a particular model setup.
*/
struct ifs_tpfa_forces {
const double *src;
const struct FlowBoundaryConditions *bc ;
const double *src; /**< Explicit source terms */
const struct FlowBoundaryConditions *bc ; /**< Boundary conditions */
const struct Wells *W ;
const double *totmob;
const double *wdp ;
const struct Wells *W ; /**< Well topology */
const double *totmob; /**< Total mobility in each cell */
const double *wdp ; /**< Gravity adjustment at each perforation */
};
/**
* Allocate TPFA management structure capable of assembling a system of
* simultaneous linear equations corresponding to a particular grid and well
* configuration.
*
* @param[in] G Grid.
* @param[in] W Well topology.
* @return Fully formed TPFA management structure if successful, @c NULL in case
* of allocation failure.
*/
struct ifs_tpfa_data *
ifs_tpfa_construct(struct UnstructuredGrid *G,
struct Wells *W);
/**
*
* @param[in] G
* @param[in] F
* @param[in] trans
* @param[in] gpress
* @param[in,out] h
* @return
*/
int
ifs_tpfa_assemble(struct UnstructuredGrid *G ,
const struct ifs_tpfa_forces *F ,

View File

@ -20,23 +20,88 @@
#ifndef OPM_TRANS_TPFA_HEADER_INCLUDED
#define OPM_TRANS_TPFA_HEADER_INCLUDED
/**
* \file
* Routines to assist in the calculation of two-point transmissibilities.
*/
#include <opm/core/grid.h>
#ifdef __cplusplus
extern "C" {
#endif
/**
* Calculate static, one-sided transmissibilities for use in the two-point flux
* approximation method.
*
* The one-sided transmissibilities are defined by the formula
* \f[
* t_i = \frac{\vec{n}_f \mathsf{K}_c \vec{c}_c}{\lVert \vec{c}_c \rVert^2}
* \f]
* in which @c i is the half-face index corresponding to the cell-face index
* pair <CODE>(c,f)</CODE> and \f$\vec{c}_{cf} = \Bar{x}_f - \Bar{x}_c\f$ is the
* centroid difference vector.
*
* @param[in] G Grid.
* @param[in] perm Permeability. One symmetric, positive definite tensor
* per grid cell.
* @param[out] htrans One-sided transmissibilities. Array of size at least
* <CODE>G->cell_facepos[ G->number_of_cells ]</CODE>.
*/
void
tpfa_htrans_compute(struct UnstructuredGrid *G, const double *perm, double *htrans);
tpfa_htrans_compute(struct UnstructuredGrid *G ,
const double *perm ,
double *htrans);
/**
* Compute two-point transmissibilities from one-sided transmissibilities.
*
* The two-point transmissibilities are given by the simple, near-harmonic
* average (save a factor of two)
* \f[
* \mathsf{T}_f = \big(\frac{1}{t_1} + \frac{1}{t_2}\big)^{-1}
* = \frac{t_1t_2}{t_1 + t_2}
* \f]
* in which \f$t_1\f$ and \f$t_2\f$ are the one-sided transmissibilities that
* connect the neighbouring cells of face @c f.
*
* @param[in] G Grid.
* @param[in] htrans One-sided transmissibilities as defined by function
* tpfa_htrans_compute().
* @param[out] trans Interface, two-point transmissibilities. Array of size
* at least <CODE>G->number_of_faces</CODE>.
*/
void
tpfa_trans_compute(struct UnstructuredGrid *G, const double *htrans, double *trans);
tpfa_trans_compute(struct UnstructuredGrid *G ,
const double *htrans,
double *trans );
/**
* Calculate effective two-point transmissibilities from one-sided, total
* mobility weighted, transmissibilities.
*
* Specifically, compute the following product
* \f[
* \mathsf{T}_f = \big(\frac{1}{\lambda_1t_1} + \frac{1}{\lambda_2t_2}\big)^{-1}
* = \lambda_1\lambda_2 \frac{t_1t_2}{t_1 + t_2}
* \f]
* in which \f$t_1\f$ and \f$t_2\f$ are the one-sided, static transmissibility
* values connecting the cells of face @c f and \f$\lambda_1\f$ and
* \f$\lambda_2\f$ denote the total mobilities of the respective cells.
*
* @param[in] G Grid.
* @param[in] totmob Total mobilities. One positive scalar value for each cell.
* @param[in] htrans One-sided transmissibilities as defined by function
* tpfa_htrans_compute().
* @param[out] trans Effective, two-point transmissibilities. Array of size at
* least <CODE>G->number_of_faces</CODE>.
*/
void
tpfa_eff_trans_compute(struct UnstructuredGrid *G,
const double *totmob,
const double *htrans,
double *trans);
tpfa_eff_trans_compute(struct UnstructuredGrid *G ,
const double *totmob,
const double *htrans,
double *trans );
#ifdef __cplusplus
}

View File

@ -47,7 +47,9 @@ namespace Opm {
public:
ImplicitAssembly(Model& model)
: model_(model)
: model_(model),
nconn_(-1) ,
asm_buffer_()
{}
template <class Grid ,

View File

@ -33,6 +33,12 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
/**
* \file
* Numerical model and support classes needed to model transport of two
* incompressible fluid phases. Intended for the ImplicitTransport system.
*/
#ifndef OPM_SINGLEPOINTUPWINDTWOPHASE_HPP_HEADER
#define OPM_SINGLEPOINTUPWINDTWOPHASE_HPP_HEADER
@ -45,8 +51,43 @@
namespace Opm {
namespace spu_2p {
/**
* Internal class to manage the direct and derived quantities needed to
* formulate the fluid transport system.
*
* Note: This class elides off-diagonal elements of the phase mobility
* Jacobian, \f$(\partial_{s_\beta} \lambda_\alpha)_{\alpha\beta}\f$.
* These elements are assumed to be strictly equal to zero. In other
* words, the relative permeability of phase \f$\alpha\f$ is assumed to
* depend only on the saturation of phase \f$\alpha\f$. This convention
* allows storing only the two diagonals of the mobility Jacobian per
* grid cell.
*
* The static gravity term is the scalar value
* \f[
* \Delta G_i = \mathsf{T}_f\, \vec{g}\cdot(\Bar{x}_f - \Bar{x}_c)
* \f]
* in which @c i is the half face index corresponding to the cell-face
* pair <CODE>(f,c)</CODE> and \f$\mathsf{T}_f\f$ is the absolute
* (bacground) two-point transmissibility of face @c f.
*
* The fluid transport problem is formulated in terms of saturation
* changes, \f$\Delta s\f$, per cell. These changes are the primary
* degrees of freedom in this model.
*
* Capillary pressures are defined by the fluid model, but usually
* correspond to \f$p_w - p_n\f$ (e.g., \f$p_\mathit{oil} -
* p_\mathit{water}\f$).
*/
class ModelParameterStorage {
public:
/**
* Constructor.
*
* @param[in] nc Total number of grid cells.
* @param[in] totconn Total number of connections, accumulated per
* cell (``half faces'').
*/
ModelParameterStorage(int nc, int totconn)
: drho_(0.0), mob_(0), dmob_(0),
porevol_(0), dg_(0), ds_(0), pc_(0), dpc_(0), trans_(0),
@ -74,44 +115,160 @@ namespace Opm {
trans_ = dpc_ + (1 * nc );
}
/**
* Modifiable density difference.
* @return Reference to modifiable internal representation of fluid
* phase density difference.
*/
double& drho () { return drho_ ; }
/**
* Read-only density difference.
* @return Read-only value of current fluid phase difference value.
*/
double drho () const { return drho_ ; }
/**
* Phase mobility in cell.
* @param[in] c Cell.
* @return Read-write reference to two consecutive phase mobilities
* in cell @c c.
*/
double* mob (int c) { return mob_ + (2*c + 0); }
/**
* Phase mobility in cell.
* @param[in] c Cell.
* @return Read-only reference to two consecutive phase mobilities
* in cell @c c.
*/
const double* mob (int c) const { return mob_ + (2*c + 0); }
/**
* Diagonal elements of phase mobility derivative (Jacobian).
*
* @param[in] c Cell.
* @return Read-write reference to diagonal elements of phase
* mobility Jacobian in cell @c c.
*/
double* dmob (int c) { return dmob_ + (2*c + 0); }
/**
* Diagonal elements of phase mobility derivative (Jacobian).
*
* @param[in] c Cell.
* @return Read-only reference to two consecutive diagonal elements
* of phase mobility Jacobian in cell @c c.
*/
const double* dmob (int c) const { return dmob_ + (2*c + 0); }
/**
* Retrieve pore volumes for all cells.
* @return Modifiable vector of pore volumes for all cells.
*/
double* porevol() { return porevol_ ; }
/**
* Pore volume of single cell.
* @param[in] c Cell.
* @return Pore volume of cell @c c.
*/
double porevol(int c) const { return porevol_[c] ; }
/**
* Static gravity term associated to single half face.
*
* @param[in] i Half face index corresponding to particular
* cell-face pair.
* @return Read-write reference to static gravity term of single
* half face.
*/
double& dg(int i) { return dg_[i] ; }
/**
* Static gravity term associated to single half face.
* @param[in] i Half face index corresponding to particular
* cell-face pair.
* @return Read-only reference to static gravity term of single half
* face.
*/
double dg(int i) const { return dg_[i] ; }
/**
* Saturation change in particular cell.
*
* @param[in] c
* @return Read-write reference to saturation change (scalar) in
* cell @c c.
*/
double& ds(int c) { return ds_[c] ; }
/**
* Saturation change in particular cell.
*
* @param[in] c
* @return Read-only reference to saturation change (scalar) in cell
* @c c.
*/
double ds(int c) const { return ds_[c] ; }
/**
* Capillary pressure in particular cell.
*
* @param[in] c Cell.
* @return Read-write reference to capillary pressure in cell @c c.
*/
double& pc(int c) { return pc_[c] ; }
/**
* Capillary pressure in particular cell.
*
* @param[in] c Cell
* @return Read-only reference to capillary pressure in cell @c c.
*/
double pc(int c) const { return pc_[c] ; }
/**
* Derivative of capillary pressure with respect to saturation.
*
* @param[in] c Cell
* @return Read-write reference to capillary pressure derivative
* with respect to primary saturation in cell @c c.
*/
double& dpc(int c) { return dpc_[c] ; }
/**
* Derivative of capillary pressure with respect to saturation.
*
* @param[in] c Cell
* @return Read-only reference to capillary pressure derivative with
* respect to primary saturation in cell @c c.
*/
double dpc(int c) const { return dpc_[c] ; }
/**
* Background (absolute) face transmissibility of particular face.
*
* @param[in] f Face
* @return Read-write reference to background face transmissibility
* of face @c f.
*/
double& trans(int f) { return trans_[f] ; }
/**
* Background (absolute) face transmissibility of particular face.
*
* @param[in] f Face
* @return Read-only reference to bacground face transmissibility of
* face @c f.
*/
double trans(int f) const { return trans_[f] ; }
private:
double drho_ ;
double *mob_ ;
double *dmob_ ;
double *porevol_;
double *dg_ ;
double *ds_ ;
double *pc_ ;
double *dpc_ ;
double *trans_ ;
double drho_ ; /**< Fluid phase density difference */
double *mob_ ; /**< Fluid phase mobility in all cells */
double *dmob_ ; /**< Derivative of phase mobility in all cells */
double *porevol_; /**< Pore volume in all cells */
double *dg_ ; /**< Static gravity term on all half faces */
double *ds_ ; /**< Saturation change in all cells */
double *pc_ ; /**< Capillary pressure in all cells */
double *dpc_ ; /**< Derivative of cap. pressure in all cells */
double *trans_ ; /**< Absolute transmissibility on all faces */
/**
* Data storage from which individual quantities are managed.
*/
std::vector<double> data_;
};
}
@ -149,42 +306,42 @@ namespace Opm {
std::copy(porevol.begin(), porevol.end(), store_.porevol());
}
void makefhfQPeriodic( const std::vector<int>& p_faces,const std::vector<int>& hf_faces,
const std::vector<int>& nb_faces)
void
makefhfQPeriodic(const std::vector<int>& p_faces ,
const std::vector<int>& hf_faces,
const std::vector<int>& nb_faces)
{
assert (! p_faces.empty());
assert (p_faces.size() == hf_faces.size());
assert (hf_faces.size() == nb_faces.size());
std::vector<int> nbhf(hf_faces.size());
for(unsigned int i=0; i<p_faces.size(); ++i){
int nbf = nb_faces[i];
if(f2hf_[2*nbf] == -1){
nbhf[i] = f2hf_[2*nbf+1];
}else{
assert(f2hf_[2*nbf+1]==-1);
nbhf[i] = f2hf_[2*nbf];
}
for (std::vector<int>::size_type i = 0; i < p_faces.size(); ++i) {
const int nbf = nb_faces[i];
assert (2*std::vector<int>::size_type(nbf) + 1 < f2hf_.size());
assert ((f2hf_[2*nbf + 0] < 0) ^ (f2hf_[2*nbf + 1] < 0));
const int p = (f2hf_[2*nbf + 0] < 0) ? 1 : 0; // "Self"
nbhf[ i ] = f2hf_[2*nbf + p];
}
for(unsigned int i=0; i<p_faces.size(); ++i){
int f = p_faces[i];
int hf = hf_faces[i];
bool changed=false;
for (std::vector<int>::size_type i = 0; i < p_faces.size(); ++i) {
const int f = p_faces [i];
const int hf = hf_faces[i];
if(f2hf_[2*f] == hf){
assert(f2hf_[2*f+1]==-1);
}else{
assert(f2hf_[2*f]==-1);
f2hf_[2*f]=nbhf[i];
changed=true;
}
if(!changed){
if(f2hf_[2*f+1]== hf){
assert(f2hf_[2*f]==-1);
}else{
assert(f2hf_[2*f+1]==-1);
f2hf_[2*f+1]=nbhf[i];
changed=true;
}
}
assert(changed);
assert (0 <= f);
assert (0 <= hf);
assert (2*std::vector<int>::size_type(f) + 1 < f2hf_.size());
assert ((f2hf_[2*f + 0] < 0 ) ^ (f2hf_[2*f + 1] < 0 ));
assert ((f2hf_[2*f + 0] == hf) ^ (f2hf_[2*f + 1] == hf));
const int p = (f2hf_[2*f + 0] == hf) ? 1 : 0; // "Other"
f2hf_[2*f + p] = nbhf[ i ];
}
}

View File

@ -57,11 +57,11 @@ namespace Opm
/// restarted by a call to start().
void stop();
/// \ret the number of running seconds that have passed
/// \return the number of running seconds that have passed
/// since last call to start(), secsSinceLast() or
/// secsSinceStart()
double secsSinceLast();
/// \ret the number of running seconds that have passed
/// \return the number of running seconds that have passed
/// since last call to start().
double secsSinceStart();

View File

@ -5,7 +5,7 @@ $(BOOST_CPPFLAGS)
LDFLAGS = $(BOOST_LDFLAGS)
LDADD = $(top_builddir)/libopmcore.la
LDADD = $(top_builddir)/lib/libopmcore.la
noinst_PROGRAMS = \

View File

@ -4,7 +4,7 @@ $(BOOST_CPPFLAGS)
LDFLAGS = $(BOOST_LDFLAGS)
LDADD = $(top_builddir)/libopmcore.la
LDADD = $(top_builddir)/lib/libopmcore.la
noinst_PROGRAMS = tutorial1
tutorial1_SOURCES = tutorial1.cpp