OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merged branch from bska.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2010-10-26 11:16:43 +02:00
commit 8f4bc05767
32 changed files with 1649 additions and 475 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
Makefile.am
View File

@ -12,17 +12,22 @@ libopmpressure_HEADERS = \
blas_lapack.h \
coarse_conn.h \
coarse_sys.h \
compr_quant.h \
dfs.h \
flow_bc.h \
fsh.h \
fsh_common.h \
grid.h \
hash_set.h \
hybsys_global.h \
hybsys.h \
ifs_tpfa.h \
ifsh.h \
ifsh_ms.h \
mimetic.h \
partition.h \
sparse_sys.h \
trans_tpfa.h \
well.h \
GridAdapter.hpp \
HybridPressureSolver.hpp
@ -31,16 +36,21 @@ HybridPressureSolver.hpp
libopmpressure_la_SOURCES = \
coarse_conn.c \
coarse_sys.c \
compr_quant.c \
dfs.c \
flow_bc.c \
fsh.c \
fsh_common.c \
hash_set.c \
hybsys.c \
hybsys_global.c \
ifs_tpfa.c \
ifsh.c \
ifsh_ms.c \
mimetic.c \
partition.c \
sparse_sys.c \
trans_tpfa.c \
well.c

6
blas_lapack.h
View File

@ -17,8 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef BLAS_LAPACK_H_INCLUDED
#define BLAS_LAPACK_H_INCLUDED
#ifndef OPM_BLAS_LAPACK_HEADER_INCLUDED
#define OPM_BLAS_LAPACK_HEADER_INCLUDED
#ifdef __cplusplus
extern "C" {
@ -101,4 +101,4 @@ double ddot_(const MAT_SIZE_T *n, const double *x, const MAT_SIZE_T *incx,
}
#endif
#endif /* BLAS_LAPACK_H_INCLUDED */
#endif /* OPM_BLAS_LAPACK_HEADER_INCLUDED */

6
coarse_conn.h
View File

@ -17,8 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef COARSE_CONN_H_INCLUDED
#define COARSE_CONN_H_INCLUDED
#ifndef OPM_COARSE_CONN_HEADER_INCLUDED
#define OPM_COARSE_CONN_HEADER_INCLUDED
#ifdef __cplusplus
extern "C" {
@ -51,4 +51,4 @@ coarse_topology_destroy(struct coarse_topology *t);
}
#endif
#endif /* COARSE_CONN_H_INCLUDED */
#endif /* OPM_COARSE_CONN_HEADER_INCLUDED */

14
coarse_sys.c
View File

@ -311,6 +311,7 @@ enumerate_active_bf(struct coarse_topology *ct,
b_out = (b1 == b_in) ? b2 : b1;
if (b_out >= 0) {
/* Restricted to internal cf's for now. */
m->bfno[cf] = act++;
}
}
@ -1053,19 +1054,6 @@ define_csr_sparsity(size_t nc, size_t m, struct bf_asm_data *bf_asm)
}
/* ---------------------------------------------------------------------- */
/* v = zeros([n, 1]) */
/* ---------------------------------------------------------------------- */
static void
vector_zero(size_t n, double *v)
/* ---------------------------------------------------------------------- */
{
size_t i;
for (i = 0; i < n; i++) { v[i] = 0.0; }
}
/* ---------------------------------------------------------------------- */
/* Assemble system of linear equations corresponding to local
* discretisation of flow problem on domain connected to coarse face

6
coarse_sys.h
View File

@ -17,8 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef COARSE_SYS_H_INCLUDED
#define COARSE_SYS_H_INCLUDED
#ifndef OPM_COARSE_SYS_HEADER_INCLUDED
#define OPM_COARSE_SYS_HEADER_INCLUDED
#ifdef __cplusplus
extern "C" {
@ -85,4 +85,4 @@ coarse_sys_compute_Binv(int nb,
}
#endif
#endif /* COARSE_SYS_H_INCLUDED */
#endif /* OPM_COARSE_SYS_HEADER_INCLUDED */

81
compr_quant.c Normal file
View File

@ -0,0 +1,81 @@
/*
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/>.
*/
#include "compr_quant.h"
/* ---------------------------------------------------------------------- */
/* Compute B \ (V') == zeta(cellNo) .* faceFlux2CellFlux(fflux) */
/* ---------------------------------------------------------------------- */
void
compr_flux_term(grid_t *G,
const double *fflux,
const double *zeta,
double *Biv)
/* ---------------------------------------------------------------------- */
{
int c, i, f;
double s;
for (c = i = 0; c < G->number_of_cells; c++) {
for (; i < G->cell_facepos[c + 1]; i++) {
f = G->cell_faces[i];
s = 2.0*(c == G->face_cells[2*f + 0]) - 1.0;
Biv[i] = zeta[c] * s * fflux[f];
}
}
}
/* ---------------------------------------------------------------------- */
/* Compute P == ct .* pv ./ dt */
/* ---------------------------------------------------------------------- */
void
compr_accum_term(size_t nc,
double dt,
const double *porevol,
const double *totcompr,
double *P)
/* ---------------------------------------------------------------------- */
{
size_t c;
for (c = 0; c < nc; c++) {
P[c] = -totcompr[c] * porevol[c] / dt;
}
}
/* ---------------------------------------------------------------------- */
/* Add pressure accumulation term (P*p_0) to cell sources */
/* ---------------------------------------------------------------------- */
void
compr_src_add_press_accum(size_t nc,
const double *p0,
const double *P,
double *src)
/* ---------------------------------------------------------------------- */
{
size_t c;
for (c = 0; c < nc; c++) {
src[c] -= P[c] * p0[c];
}
}

46
compr_quant.h Normal file
View File

@ -0,0 +1,46 @@
/*
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_COMPR_QUANT_HEADER_INCLUDED
#define OPM_COMPR_QUANT_HEADER_INCLUDED
#include <stddef.h>
#include "grid.h"
void
compr_flux_term(grid_t *G,
const double *fflux,
const double *zeta,
double *Biv);
void
compr_accum_term(size_t nc,
double dt,
const double *porevol,
const double *totcompr,
double *P);
void
compr_src_add_press_accum(size_t nc,
const double *p0,
const double *P,
double *src);
#endif /* OPM_COMPR_QUANT_HEADER_INCLUDED */

4
dfs.h
View File

@ -17,8 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef DFS_H_INCLUDED
#define DFS_H_INCLUDED
#ifndef OPM_DFS_HEADER_INCLUDED
#define OPM_DFS_HEADER_INCLUDED
#ifdef __cplusplus
extern "C" {

6
flow_bc.h
View File

@ -17,8 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef FLOW_BC_H_INCLUDED
#define FLOW_BC_H_INCLUDED
#ifndef OPM_FLOW_BC_HEADER_INCLUDED
#define OPM_FLOW_BC_HEADER_INCLUDED
#include <stddef.h>
@ -46,4 +46,4 @@ deallocate_flowbc(flowbc_t *fbc);
}
#endif
#endif /* FLOW_BC_H_INCLUDED */
#endif /* OPM_FLOW_BC_HEADER_INCLUDED */

304
fsh.c Normal file
View File

@ -0,0 +1,304 @@
/*
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/>.
*/
#include <assert.h>
#include <limits.h>
#include <math.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include "fsh_common.h"
#include "fsh.h"
#include "fsh_common_impl.h"
#include "hybsys.h"
#include "hybsys_global.h"
#if defined MAX
#undef MAX
#endif
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
/* ---------------------------------------------------------------------- */
static void
fsh_compute_table_sz(grid_t *G, well_t *W, int max_ngconn,
size_t *nnu, size_t *idata_sz, size_t *ddata_sz)
/* ---------------------------------------------------------------------- */
{
int nc, ngconn_tot;
*nnu = G->number_of_faces;
nc = G->number_of_cells;
ngconn_tot = G->cell_facepos[nc];
*idata_sz = nc + 1; /* gdof_pos */
*idata_sz += ngconn_tot; /* gdof */
*idata_sz += max_ngconn; /* iwork */
*ddata_sz = 2 * (*nnu); /* rhs + soln */
*ddata_sz += ngconn_tot; /* cflux */
*ddata_sz += max_ngconn; /* work */
if (W != NULL) {
*nnu += W->number_of_wells;
/* cwell_pos */
*idata_sz += nc + 1;
/* cwells */
*idata_sz += 2 * W->well_connpos[ W->number_of_wells ];
/* rhs + soln */
*ddata_sz += 2 * W->number_of_wells;
/* WI, wdp */
*ddata_sz += 2 * W->well_connpos[ W->number_of_wells ];
}
}
/* ---------------------------------------------------------------------- */
static int
fsh_assemble_grid(flowbc_t *bc,
const double *Binv,
const double *gpress,
const double *src,
struct fsh_data *h)
/* ---------------------------------------------------------------------- */
{
int c, n, nc, p1, p2;
int npp;
int *pgconn, *gconn;
nc = h->pimpl->nc;
pgconn = h->pimpl->gdof_pos;
gconn = h->pimpl->gdof;
p1 = p2 = npp = 0;
for (c = 0; c < nc; c++) {
n = pgconn[c + 1] - pgconn[c];
hybsys_cellcontrib_unsymm(c, n, p1, p2, gpress, src, Binv,
h->pimpl->sys);
npp += fsh_impose_bc(n, gconn + p1, bc, h->pimpl);
hybsys_global_assemble_cell(n, gconn + p1,
h->pimpl->sys->S,
h->pimpl->sys->r, h->A, h->b);
p1 += n;
p2 += n * n;
}
return npp;
}
/* ======================================================================
* Public routines follow.
* ====================================================================== */
/* ---------------------------------------------------------------------- */
/* Allocate and define supporting structures for assembling the global
* system of linear equations to couple the grid (reservoir)
* connections represented by 'G' and, if present (i.e., non-NULL),
* the well connections represented by 'W'. */
/* ---------------------------------------------------------------------- */
struct fsh_data *
fsh_construct(grid_t *G, well_t *W)
/* ---------------------------------------------------------------------- */
{
int nc, ngconn_tot;
size_t idata_sz, ddata_sz, nnu;
struct fsh_data *new;
assert (G != NULL);
/* Allocate master structure, define system matrix sparsity */
new = malloc(1 * sizeof *new);
if (new != NULL) {
new->A = hybsys_define_globconn(G, W);
new->pimpl = NULL;
if (new->A == NULL) {
fsh_destroy(new);
new = NULL;
}
}
/* Allocate implementation structure */
if (new != NULL) {
fsh_count_grid_dof(G, &new->max_ngconn, &new->sum_ngconn2);
fsh_compute_table_sz(G, W, new->max_ngconn,
&nnu, &idata_sz, &ddata_sz);
new->pimpl = fsh_impl_allocate_basic(idata_sz, ddata_sz);
if (new->pimpl == NULL) {
fsh_destroy(new);
new = NULL;
}
}
/* Allocate Schur complement contributions. Unsymmetric system. */
if (new != NULL) {
nc = G->number_of_cells;
ngconn_tot = G->cell_facepos[nc];
fsh_define_linsys_arrays(new);
fsh_define_impl_arrays(nc, nnu, ngconn_tot, new->max_ngconn,
W, new->pimpl);
new->pimpl->sys = hybsys_allocate_unsymm(new->max_ngconn,
nc, ngconn_tot);
if (W != NULL) {
fsh_define_cell_wells(nc, W, new->pimpl);
new->pimpl->wsys =
hybsys_well_allocate_unsymm(new->max_ngconn, nc,
new->pimpl->cwell_pos);
}
if ((new->pimpl->sys == NULL) ||
((W != NULL) && (new->pimpl->wsys == NULL))) {
/* Failed to allocate ->sys or ->wsys (if W != NULL) */
fsh_destroy(new);
new = NULL;
}
}
if (new != NULL) {
/* All allocations succeded. Fill metadata and return. */
new->pimpl->nc = nc;
new->pimpl->nf = G->number_of_faces;
new->pimpl->nw = (W != NULL) ? W->number_of_wells : 0;
memcpy(new->pimpl->gdof_pos,
G->cell_facepos ,
(nc + 1) * sizeof *new->pimpl->gdof_pos);
memcpy(new->pimpl->gdof ,
G->cell_faces ,
ngconn_tot * sizeof *new->pimpl->gdof);
hybsys_init(new->max_ngconn, new->pimpl->sys);
}
return new;
}
/* ---------------------------------------------------------------------- */
/* Assemble global system of linear equations
*
* fsh->A * fsh->x = fsh->b
*/
/* ---------------------------------------------------------------------- */
void
fsh_assemble(flowbc_t *bc,
const double *src,
const double *Binv,
const double *Biv,
const double *P,
const double *gpress,
well_control_t *wctrl,
const double *WI,
const double *BivW,
const double *wdp,
struct fsh_data *h)
/* ---------------------------------------------------------------------- */
{
int npp; /* Number of prescribed pressure values */
hybsys_schur_comp_unsymm(h->pimpl->nc,
h->pimpl->gdof_pos,
Binv, Biv, P, h->pimpl->sys);
npp = fsh_assemble_grid(bc, Binv, gpress, src, h);
if (npp == 0) {
h->A->sa[0] *= 2; /* Remove zero eigenvalue */
}
}
/* ---------------------------------------------------------------------- */
/* Compute cell pressures (cpress) and interface fluxes (fflux) from
* current solution of system of linear equations, h->x. Back
* substitution process, projected half-contact fluxes. */
/* ---------------------------------------------------------------------- */
void
fsh_press_flux(grid_t *G,
const double *Binv, const double *gpress,
struct fsh_data *h,
double *cpress, double *fflux,
double *wpress, double *wflux)
/* ---------------------------------------------------------------------- */
{
int c, f, i;
double s;
hybsys_compute_press_flux(G->number_of_cells,
G->cell_facepos,
G->cell_faces,
gpress, Binv,
h->pimpl->sys,
h->x, cpress, h->pimpl->cflux,
h->pimpl->work);
if (h->pimpl->nw > 0) {
assert ((wpress != NULL) && (wflux != NULL));
hybsys_compute_press_flux_well(G->number_of_cells, G->cell_facepos,
G->number_of_faces, h->pimpl->nw,
h->pimpl->cwell_pos, h->pimpl->cwells,
Binv, h->pimpl->WI,
h->pimpl->wdp, h->pimpl->sys,
h->pimpl->wsys, h->x, cpress,
h->pimpl->cflux, wpress, wflux,
h->pimpl->work);
}
for (f = 0; f < G->number_of_faces; f++) { fflux[f] = 0.0; }
i = 0;
for (c = 0; c < G->number_of_cells; c++) {
for (; i < G->cell_facepos[c + 1]; i++) {
f = G->cell_faces[i];
s = 2.0*(G->face_cells[2*f + 0] == c) - 1.0;
fflux[f] += s * h->pimpl->cflux[i];
}
}
for (f = 0; f < G->number_of_faces; f++) {
i = (G->face_cells[2*f + 0] >= 0) +
(G->face_cells[2*f + 1] >= 0);
fflux[f] /= i;
}
}

71
fsh.h Normal file
View File

@ -0,0 +1,71 @@
/*
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_FSH_HEADER_INCLUDED
#define OPM_FHS_HEADER_INCLUDED
#include "grid.h"
#include "well.h"
#include "flow_bc.h"
#ifdef __cplusplus
extern "C" {
#endif
struct fsh_data;
/** Constructs incompressible hybrid flow-solver data object for a
* given grid and well pattern.
*/
struct fsh_data *
fsh_construct(grid_t *G, well_t *W);
/** Assembles the hybridized linear system for face pressures.
*/
void
fsh_assemble(flowbc_t *bc,
const double *src,
const double *Binv,
const double *Biv,
const double *P,
const double *gpress,
well_control_t *wctrl,
const double *WI,
const double *BivW,
const double *wdp,
struct fsh_data *h);
/** Computes cell pressures, face fluxes, well pressures and well
* fluxes from face pressures.
*/
void
fsh_press_flux(grid_t *G,
const double *Binv, const double *gpress,
struct fsh_data *h,
double *cpress, double *fflux,
double *wpress, double *wflux);
#ifdef __cplusplus
}
#endif
#endif /* OPM_IFSH_HEADER_INCLUDED */

237
fsh_common.c Normal file
View File

@ -0,0 +1,237 @@
/*
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/>.
*/
#include <assert.h>
#include <limits.h>
#include <math.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include "grid.h"
#include "flow_bc.h"
#include "well.h"
#include "fsh_common.h"
#include "fsh_common_impl.h"
#include "hybsys.h"
#include "hybsys_global.h"
#if defined MAX
#undef MAX
#endif
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
/* ---------------------------------------------------------------------- */
/* Release dynamic memory resources associated to internal data of a
* particular (incompressible) flow solver instance. */
/* ---------------------------------------------------------------------- */
static void
fsh_destroy_impl(struct fsh_impl *pimpl)
/* ---------------------------------------------------------------------- */
{
if (pimpl != NULL) {
hybsys_well_free(pimpl->wsys );
hybsys_free (pimpl->sys );
free (pimpl->ddata);
free (pimpl->idata);
}
free(pimpl);
}
/* ---------------------------------------------------------------------- */
/* Eliminate 'npp' prescribed (Dirichlet) conditions (locdof,dofval)
* from global system of linear equations. Move known values to RHS
* whilst zeroing coefficient matrix contributions. Local system of
* dimension 'n'. */
/* ---------------------------------------------------------------------- */
static void
fsh_explicit_elimination(int n, int npp,
int *locdof, double *dofval,
double *S, double *r)
/* ---------------------------------------------------------------------- */
{
int i, p;
for (i = 0; i < npp; i++) {
for (p = (locdof[i] + 1) % n; p != locdof[i]; p = (p + 1) % n) {
/* Subtract known values from RHS. */
r[p] -= S[p + locdof[i]*n] * dofval[i];
/* Eliminate DOF (zero row/col "locdof[i]"; leave diagonal). */
S[p + locdof[i]*n] = 0.0;
S[locdof[i] + p*n] = 0.0;
}
/* Produce trivial equation S(i,i)*x(i) = S(i,i)*x0(i). */
r[p] = S[p * (n + 1)] * dofval[i];
}
}
/* ---------------------------------------------------------------------- */
/* Release memory resources for ifsh data-handle 'h' */
/* ---------------------------------------------------------------------- */
void
fsh_destroy(struct fsh_data *h)
/* ---------------------------------------------------------------------- */
{
if (h != NULL) {
fsh_destroy_impl(h->pimpl);
csrmatrix_delete(h->A );
}
free(h);
}
/* ---------------------------------------------------------------------- */
struct fsh_impl *
fsh_impl_allocate_basic(size_t idata_sz, size_t ddata_sz)
/* ---------------------------------------------------------------------- */
{
struct fsh_impl *new;
new = malloc(1 * sizeof *new);
if (new != NULL) {
new->idata = malloc(idata_sz * sizeof *new->idata);
new->ddata = malloc(ddata_sz * sizeof *new->ddata);
new->sys = NULL;
new->wsys = NULL;
if ((new->idata == NULL) || (new->ddata == NULL)) {
fsh_destroy_impl(new);
new = NULL;
}
}
return new;
}
/* ---------------------------------------------------------------------- */
/* Determine nnz (=sum(diff(facePos)^2)) and max(diff(facePos) for grid */
/* ---------------------------------------------------------------------- */
void
fsh_count_grid_dof(grid_t *G, int *max_ngdof, size_t *sum_ngdof2)
/* ---------------------------------------------------------------------- */
{
int c, n;
*max_ngdof = INT_MIN;
*sum_ngdof2 = 0;
for (c = 0; c < G->number_of_cells; c++) {
n = G->cell_facepos[c + 1] - G->cell_facepos[c];
*max_ngdof = MAX(*max_ngdof, n);
*sum_ngdof2 += n * n;
}
}
/* ---------------------------------------------------------------------- */
/* Impose boundary conditions on local contribution to global system. */
/* ---------------------------------------------------------------------- */
int
fsh_impose_bc(int nconn, int *conn, flowbc_t *bc,
struct fsh_impl *pimpl)
/* ---------------------------------------------------------------------- */
{
int i, npp, f;
npp = 0;
for (i = 0; i < nconn; i++) {
f = conn[i];
if (bc->type[f] == PRESSURE) {
pimpl->work [npp] = bc->bcval[f];
pimpl->iwork[npp] = i;
npp += 1;
} else if (bc->type[f] == FLUX) {
pimpl->sys->r[i] -= bc->bcval[f];
}
}
if (npp > 0) {
fsh_explicit_elimination(nconn, npp,
pimpl->iwork,
pimpl->work,
pimpl->sys->S,
pimpl->sys->r);
}
return npp;
}
/* ---------------------------------------------------------------------- */
void
fsh_define_impl_arrays(size_t nc,
size_t nnu,
size_t nhf,
size_t max_ncf,
well_t *W,
struct fsh_impl *pimpl)
/* ---------------------------------------------------------------------- */
{
pimpl->cflux = pimpl->ddata + 2 * nnu;
pimpl->work = pimpl->cflux + nhf;
pimpl->gdof_pos = pimpl->idata;
pimpl->gdof = pimpl->gdof_pos + (nc + 1);
pimpl->iwork = pimpl->gdof + nhf;
if (W != NULL) {
pimpl->cwell_pos = pimpl->iwork + max_ncf;
pimpl->cwells = pimpl->cwell_pos + nc + 1;
pimpl->WI = pimpl->work + max_ncf;
pimpl->wdp = pimpl->WI + W->well_connpos[ W->number_of_wells ];
}
}
/* ---------------------------------------------------------------------- */
void
fsh_define_cell_wells(size_t nc, well_t *W, struct fsh_impl *pimpl)
/* ---------------------------------------------------------------------- */
{
memset(pimpl->cwell_pos, 0, (nc + 1) * sizeof *pimpl->cwell_pos);
derive_cell_wells(nc, W, pimpl->cwell_pos, pimpl->cwells);
}
/* ---------------------------------------------------------------------- */
void
fsh_define_linsys_arrays(struct fsh_data *h)
/* ---------------------------------------------------------------------- */
{
h->b = h->pimpl->ddata;
h->x = h->b + h->A->m;
}

61
fsh_common.h Normal file
View File

@ -0,0 +1,61 @@
/*
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_IFSH_HEADER_INCLUDED
#define OPM_IFHS_HEADER_INCLUDED
#ifdef __cplusplus
extern "C" {
#endif
struct CSRMatrix;
struct fsh_impl;
/** Contains the linear system for assembly, as well as internal data
* for the assembly routines.
*/
struct fsh_data {
/* Let \f$n_i\f$ be the number of connections/faces of grid cell
* number \f$i\f$. Then max_ngconn = \f$\max_i n_i\f$
*/
int max_ngconn;
/* With n_i as above, sum_ngconn2 = \f$\sum_i n_i^2\f$ */
size_t sum_ngconn2;
/* Linear system */
struct CSRMatrix *A; /* Coefficient matrix */
double *b; /* System RHS */
double *x; /* Solution */
/* Private implementational details. */
struct fsh_impl *pimpl;
};
/** Destroys the fsh data object */
void
fsh_destroy(struct fsh_data *h);
#ifdef __cplusplus
}
#endif
#endif /* OPM_IFSH_HEADER_INCLUDED */

79
fsh_common_impl.h Normal file
View File

@ -0,0 +1,79 @@
/*
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_FSH_COMMON_IMPL_HEADER_INCLUDED
#define OPM_FSH_COMMON_IMPL_HEADER_INCLUDED
/* Internal header. Don't install. */
struct fsh_impl {
int nc, nf, nw; /* Number of cells, faces, wells */
/* Topology */
int *gdof_pos; /* Pointers, grid DOFs (== cell_facepos) */
int *gdof; /* Grid DOFs (== cell_faces) */
int *cwell_pos; /* Start pointers, well DOFs (c->w) */
int *cwells; /* Well DOFs (c->w) */
/* Discretisation data */
double *WI; /* Permuted well production indices */
double *wdp; /* Permuted well gravity pressures */
double *cflux; /* Cell (half-contact) fluxes */
struct hybsys *sys; /* Hybrid cell contribs */
struct hybsys_well *wsys; /* Hybrid cell contribs from wells */
double *work; /* Scratch array, floating point */
int *iwork; /* Scratch array, integers */
/* Linear storage goes here... */
int *idata; /* Actual storage array, integers */
double *ddata; /* Actual storage array, floating point */
};
struct fsh_impl *
fsh_impl_allocate_basic(size_t idata_sz, size_t ddata_sz);
void
fsh_count_grid_dof(grid_t *G, int *max_ngdof, size_t *sum_ngdof2);
int
fsh_impose_bc(int ndof,
int *dof,
flowbc_t *bc,
struct fsh_impl *pimpl);
void
fsh_define_impl_arrays(size_t nc,
size_t nnu,
size_t nhf,
size_t max_ncf,
well_t *W,
struct fsh_impl *pimpl);
void
fsh_define_cell_wells(size_t nc, well_t *W, struct fsh_impl *pimpl);
void
fsh_define_sys_arrays(struct fsh_data *h);
#endif /* OPM_FSH_COMMON_IMPL_HEADER_INCLUDED */

6
grid.h
View File

@ -17,8 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_GRID_H_INCLUDED
#define OPM_GRID_H_INCLUDED
#ifndef OPM_GRID_HEADER_INCLUDED
#define OPM_GRID_HEADER_INCLUDED
#ifdef __cplusplus
@ -66,4 +66,4 @@ typedef struct {
}
#endif
#endif /* OPM_GRID_H_INCLUDED */
#endif /* OPM_GRID_HEADER_INCLUDED */

6
hash_set.h
View File

@ -17,8 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef HASH_SET_H_INCLUDED
#define HASH_SET_H_INCLUDED
#ifndef OPM_HASH_SET_HEADER_INCLUDED
#define OPM_HASH_SET_HEADER_INCLUDED
#ifdef __cplusplus
extern "C" {
@ -67,4 +67,4 @@ hash_set_count_elms(const struct hash_set *set);
}
#endif
#endif /* HASH_SET_H_INCLUDED */
#endif /* OPM_HASH_SET_HEADER_INCLUDED */

6
hybsys.h
View File

@ -17,8 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef HYBSYS_H_INCLUDED
#define HYBSYS_H_INCLUDED
#ifndef OPM_HYBSYS_HEADER_INCLUDED
#define OPM_HYBSYS_HEADER_INCLUDED
#ifdef __cplusplus
extern "C" {
@ -154,4 +154,4 @@ hybsys_compute_press_flux_well(int nc, const int *pgconn, int nf,
}
#endif
#endif /* HYBSYS_H_INCLUDED */
#endif /* OPM_HYBSYS_HEADER_INCLUDED */

6
hybsys_global.h
View File

@ -17,8 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef HYBSYS_GLOBAL_H_INCLUDED
#define HYBSYS_GLOBAL_H_INCLUDED
#ifndef OPM_HYBSYS_GLOBAL_HEADER_INCLUDED
#define OPM_HYBSYS_GLOBAL_HEADER_INCLUDED
#ifdef __cplusplus
extern "C" {
@ -55,4 +55,4 @@ hybsys_global_assemble_well_sym(int ngconn_tot,
}
#endif
#endif /* HYBSYS_GLOBAL_H_INCLUDED */
#endif /* OPM_HYBSYS_GLOBAL_HEADER_INCLUDED */

235
ifs_tpfa.c Normal file
View File

@ -0,0 +1,235 @@
#include <assert.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include "ifs_tpfa.h"
#include "sparse_sys.h"
struct ifs_tpfa_impl {
double *fgrav; /* Accumulated grav contrib/face */
/* Linear storage */
double *ddata;
};
/* ---------------------------------------------------------------------- */
static void
impl_deallocate(struct ifs_tpfa_impl *pimpl)
/* ---------------------------------------------------------------------- */
{
if (pimpl != NULL) {
free(pimpl->ddata);
}
free(pimpl);
}
/* ---------------------------------------------------------------------- */
static struct ifs_tpfa_impl *
impl_allocate(grid_t *G)
/* ---------------------------------------------------------------------- */
{
struct ifs_tpfa_impl *new;
size_t ddata_sz;
ddata_sz = 2 * G->number_of_cells; /* b, x */
ddata_sz += 1 * G->number_of_faces; /* fgrav */
new = malloc(1 * sizeof *new);
if (new != NULL) {
new->ddata = malloc(ddata_sz * sizeof *new->ddata);
if (new->ddata == NULL) {
impl_deallocate(new);
new = NULL;
}
}
return new;
}
/* ---------------------------------------------------------------------- */
static struct CSRMatrix *
ifs_tpfa_construct_matrix(grid_t *G)
/* ---------------------------------------------------------------------- */
{
int f, c1, c2;
size_t nnz;
struct CSRMatrix *A;
A = csrmatrix_new_count_nnz(G->number_of_cells);
if (A != NULL) {
/* Self connections */
for (c1 = 0; c1 < G->number_of_cells; c1++) {
A->ia[ c1 + 1 ] = 1;
}
/* Other connections */
for (f = 0; f < G->number_of_faces; f++) {
c1 = G->face_cells[2*f + 0];
c2 = G->face_cells[2*f + 1];
if ((c1 >= 0) && (c2 >= 0)) {
A->ia[ c1 + 1 ] += 1;
A->ia[ c2 + 1 ] += 1;
}
}
nnz = csrmatrix_new_elms_pushback(A);
if (nnz == 0) {
csrmatrix_delete(A);
A = NULL;
}
}
if (A != NULL) {
/* Fill self connections */
for (c1 = 0; c1 < G->number_of_cells; c1++) {
A->ja[ A->ia[ c1 + 1 ] ++ ] = c1;
}
/* Fill other connections */
for (f = 0; f < G->number_of_faces; f++) {
c1 = G->face_cells[2*f + 0];
c2 = G->face_cells[2*f + 1];
if ((c1 >= 0) && (c2 >= 0)) {
A->ja[ A->ia[ c1 + 1 ] ++ ] = c2;
A->ja[ A->ia[ c2 + 1 ] ++ ] = c1;
}
}
assert ((size_t) A->ia[ G->number_of_cells ] == nnz);
/* Guarantee sorted rows */
csrmatrix_sortrows(A);
}
return A;
}
/* ======================================================================
* Public interface below separator.
* ====================================================================== */
/* ---------------------------------------------------------------------- */
struct ifs_tpfa_data *
ifs_tpfa_construct(grid_t *G)
/* ---------------------------------------------------------------------- */
{
struct ifs_tpfa_data *new;
new = malloc(1 * sizeof *new);
if (new != NULL) {
new->pimpl = impl_allocate(G);
new->A = ifs_tpfa_construct_matrix(G);
if ((new->pimpl == NULL) || (new->A == NULL)) {
ifs_tpfa_destroy(new);
new = NULL;
}
}
if (new != NULL) {
new->b = new->pimpl->ddata;
new->x = new->b + new->A->m;
new->pimpl->fgrav = new->x + new->A->m;
}
return new;
}
/* ---------------------------------------------------------------------- */
void
ifs_tpfa_assemble(grid_t *G,
const double *trans,
const double *src,
const double *gpress,
struct ifs_tpfa_data *h)
/* ---------------------------------------------------------------------- */
{
int c1, c2, c, i, f, j1, j2;
csrmatrix_zero( h->A);
vector_zero (h->A->m, h->b);
vector_zero (G->number_of_faces, h->pimpl->fgrav);
for (c = i = 0; c < G->number_of_cells; c++) {
j1 = csrmatrix_elm_index(c, c, h->A);
for (; i < G->cell_facepos[c + 1]; i++) {
f = G->cell_faces[i];
c1 = G->face_cells[2*f + 0];
c2 = G->face_cells[2*f + 1];
c2 = (c1 == c) ? c2 : c1;
if (c2 >= 0) {
j2 = csrmatrix_elm_index(c, c2, h->A);
h->A->sa[j1] += trans[f];
h->A->sa[j2] -= trans[f];
}
}
h->b[c] += src[c];
}
h->A->sa[0] *= 2;
}
/* ---------------------------------------------------------------------- */
void
ifs_tpfa_press_flux(grid_t *G,
const double *trans,
struct ifs_tpfa_data *h,
double *cpress,
double *fflux)
/* ---------------------------------------------------------------------- */
{
int c1, c2, f;
/* Assign cell pressure directly from solution vector */
memcpy(cpress, h->x, G->number_of_cells * sizeof *cpress);
for (f = 0; f < G->number_of_faces; f++) {
c1 = G->face_cells[2*f + 0];
c2 = G->face_cells[2*f + 1];
if ((c1 >= 0) && (c2 >= 0)) {
fflux[f] = trans[f] * (cpress[c1] - cpress[c2]);
} else {
fflux[f] = 0.0;
}
}
}
/* ---------------------------------------------------------------------- */
void
ifs_tpfa_destroy(struct ifs_tpfa_data *h)
/* ---------------------------------------------------------------------- */
{
if (h != NULL) {
csrmatrix_delete(h->A);
impl_deallocate (h->pimpl);
}
free(h);
}

57
ifs_tpfa.h Normal file
View File

@ -0,0 +1,57 @@
/*
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_IFS_TPFA_HEADER_INCLUDED
#define OPM_IFS_TPFA_HEADER_INCLUDED
#include "grid.h"
struct ifs_tpfa_impl;
struct CSRMatrix;
struct ifs_tpfa_data {
struct CSRMatrix *A;
double *b;
double *x;
struct ifs_tpfa_impl *pimpl;
};
struct ifs_tpfa_data *
ifs_tpfa_construct(grid_t *G);
void
ifs_tpfa_assemble(grid_t *G,
const double *trans,
const double *src,
const double *gpress,
struct ifs_tpfa_data *h);
void
ifs_tpfa_press_flux(grid_t *G,
const double *trans,
struct ifs_tpfa_data *h,
double *cpress,
double *fflux);
void
ifs_tpfa_destroy(struct ifs_tpfa_data *h);
#endif /* OPM_IFS_TPFA_HEADER_INCLUDED */

522
ifsh.c
View File

@ -24,7 +24,9 @@
#include <stdlib.h>
#include <string.h>
#include "fsh_common.h"
#include "ifsh.h"
#include "fsh_common_impl.h"
#include "hybsys.h"
#include "hybsys_global.h"
@ -33,59 +35,10 @@
#endif
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
struct ifsh_impl {
int nc, nf, nw; /* Number of cells, faces, wells */
double *Binv; /* Effective inverse IP */
double *gpress; /* Effective gravity pressure */
double *cflux; /* Cell (half-contact) fluxes */
double *WI; /* Effective inverse IP, wells */
double *wdp; /* Effective grav. press, wells */
int *pgconn; /* Start pointers, grid connections */
int *gconn; /* Grid connections */
int *cwpos; /* Start pointers, c->w mapping */
int *cwells; /* c->w mapping */
struct hybsys *sys; /* Hybrid cell contribs */
struct hybsys_well *wsys; /* Hybrid cell contribs from wells */
double *work; /* Scratch array, floating point */
int *iwork; /* Scratch array, integers */
/* Linear storage goes here... */
int *idata; /* Actual storage array, integers */
double *ddata; /* Actual storage array, floating point */
};
/* ---------------------------------------------------------------------- */
/* Determine nnz (=sum(diff(facePos)^2)) and max(diff(facePos) for grid */
/* ---------------------------------------------------------------------- */
static void
count_grid_connections(grid_t *G, int *max_ngconn, size_t *sum_ngconn2)
/* ---------------------------------------------------------------------- */
{
int c, n;
*max_ngconn = INT_MIN;
*sum_ngconn2 = 0;
for (c = 0; c < G->number_of_cells; c++) {
n = G->cell_facepos[c + 1] - G->cell_facepos[c];
*max_ngconn = MAX(*max_ngconn, n);
*sum_ngconn2 += n * n;
}
}
/* ---------------------------------------------------------------------- */
static void
ifsh_compute_table_sz(grid_t *G, well_t *W,
int max_ngconn, size_t sum_ngconn2,
ifsh_compute_table_sz(grid_t *G, well_t *W, int max_ngconn,
size_t *nnu, size_t *idata_sz, size_t *ddata_sz)
/* ---------------------------------------------------------------------- */
{
@ -96,18 +49,18 @@ ifsh_compute_table_sz(grid_t *G, well_t *W,
nc = G->number_of_cells;
ngconn_tot = G->cell_facepos[nc];
*idata_sz = max_ngconn; /* iwork */
*idata_sz = nc + 1; /* gdof_pos */
*idata_sz += ngconn_tot; /* gdof */
*idata_sz += max_ngconn; /* iwork */
*ddata_sz = 2 * (*nnu); /* rhs + soln */
*ddata_sz += sum_ngconn2; /* Binv */
*ddata_sz += ngconn_tot; /* gpress */
*ddata_sz += ngconn_tot; /* cflux */
*ddata_sz += max_ngconn; /* work */
if (W != NULL) {
*nnu += W->number_of_wells;
/* cwpos */
/* cwell_pos */
*idata_sz += nc + 1;
/* cwells */
@ -122,277 +75,59 @@ ifsh_compute_table_sz(grid_t *G, well_t *W,
}
/* ---------------------------------------------------------------------- */
/* Allocate and define supporting structures for assembling the global
* system of linear equations to couple the grid (reservoir)
* connections represented by 'G' and, if present (i.e., non-NULL),
* the well connections represented by 'W'. */
/* ---------------------------------------------------------------------- */
struct ifsh_data *
ifsh_construct(grid_t *G, well_t *W)
/* ---------------------------------------------------------------------- */
{
int nc, ngconn_tot;
size_t idata_sz, ddata_sz, nnu;
struct ifsh_data *new;
assert (G != NULL);
new = malloc(1 * sizeof *new);
if (new != NULL) {
new->A = hybsys_define_globconn(G, W);
new->pimpl = NULL;
if (new->A == NULL) {
ifsh_destroy(new);
new = NULL;
}
}
if (new != NULL) {
new->pimpl = malloc(1 * sizeof *new->pimpl);
if (new->pimpl == NULL) {
ifsh_destroy(new);
new = NULL;
} else {
new->pimpl->sys = NULL; new->pimpl->wsys = NULL;
new->pimpl->idata = NULL; new->pimpl->ddata = NULL;
}
}
if (new != NULL) {
count_grid_connections(G, &new->max_ngconn, &new->sum_ngconn2);
ifsh_compute_table_sz(G, W, new->max_ngconn, new->sum_ngconn2,
&nnu, &idata_sz, &ddata_sz);
nc = G->number_of_cells;
ngconn_tot = G->cell_facepos[nc];
new->pimpl->idata = malloc(idata_sz * sizeof *new->pimpl->idata);
new->pimpl->ddata = malloc(ddata_sz * sizeof *new->pimpl->ddata);
new->pimpl->sys = hybsys_allocate_symm(new->max_ngconn,
nc, ngconn_tot);
if ((new->pimpl->idata == NULL) ||
(new->pimpl->ddata == NULL) || (new->pimpl->sys == NULL)) {
ifsh_destroy(new);
new = NULL;
} else {
new->b = new->pimpl->ddata;
new->x = new->b + nnu;
new->pimpl->Binv = new->x + nnu;
new->pimpl->gpress = new->pimpl->Binv + new->sum_ngconn2;
new->pimpl->cflux = new->pimpl->gpress + ngconn_tot;
new->pimpl->work = new->pimpl->cflux + ngconn_tot;
new->pimpl->iwork = new->pimpl->idata;
hybsys_init(new->max_ngconn, new->pimpl->sys);
}
}
if ((new != NULL) && (W != NULL)) {
new->pimpl->cwpos = new->pimpl->iwork + new->max_ngconn;
new->pimpl->cwells = new->pimpl->cwpos + nc + 1;
memset(new->pimpl->cwpos, 0, (nc + 1) * sizeof *new->pimpl->cwpos);
derive_cell_wells(nc, W, new->pimpl->cwpos, new->pimpl->cwells);
new->pimpl->wsys = hybsys_well_allocate_symm(new->max_ngconn, nc,
new->pimpl->cwpos);
if (new->pimpl->wsys == NULL) {
ifsh_destroy(new);
new = NULL;
} else {
new->pimpl->WI = new->pimpl->work + new->max_ngconn;
new->pimpl->wdp = new->pimpl->WI + W->well_connpos[ W->number_of_wells ];
}
}
if (new != NULL) {
new->pimpl->nc = nc;
new->pimpl->nf = G->number_of_faces;
new->pimpl->nw = (W != NULL) ? W->number_of_wells : 0;
/* Contentious. Imposes severe restrictions on lifetime(G). */
new->pimpl->pgconn = G->cell_facepos;
new->pimpl->gconn = G->cell_faces;
}
return new;
}
/* ---------------------------------------------------------------------- */
/* Release dynamic memory resources associated to internal data of a
* particular (incompressible) flow solver instance. */
/* ---------------------------------------------------------------------- */
static void
ifsh_destroy_impl(struct ifsh_impl *pimpl)
ifsh_set_effective_well_params(const double *WI,
const double *wdp,
struct fsh_data *ifsh)
/* ---------------------------------------------------------------------- */
{
if (pimpl != NULL) {
hybsys_well_free(pimpl->wsys );
hybsys_free (pimpl->sys );
free (pimpl->ddata);
free (pimpl->idata);
}
int c, nc, i, perf;
int *cwpos, *cwells;
double *wsys_WI, *wsys_wdp;
free(pimpl);
}
nc = ifsh->pimpl->nc;
cwpos = ifsh->pimpl->cwell_pos;
cwells = ifsh->pimpl->cwells;
wsys_WI = ifsh->pimpl->WI;
wsys_wdp = ifsh->pimpl->wdp;
for (c = i = 0; c < nc; c++) {
for (; i < cwpos[c + 1]; i++) {
perf = cwells[2*i + 1];
/* ---------------------------------------------------------------------- */
/* Release memory resources for ifsh data-handle 'h' */
/* ---------------------------------------------------------------------- */
void
ifsh_destroy(struct ifsh_data *h)
/* ---------------------------------------------------------------------- */
{
if (h != NULL) {
ifsh_destroy_impl(h->pimpl);
csrmatrix_delete (h->A );
}
free(h);
}
/* ---------------------------------------------------------------------- */
/* Eliminate 'npp' prescribed (Dirichlet) conditions (locdof,dofval)
* from global system of linear equations. Move known values to RHS
* whilst zeroing coefficient matrix contributions. Local system of
* dimension 'n'. */
/* ---------------------------------------------------------------------- */
static void
ifsh_explicit_elimination(int n, int npp,
int *locdof, double *dofval,
double *S, double *r)
/* ---------------------------------------------------------------------- */
{
int i, p;
for (i = 0; i < npp; i++) {
for (p = (locdof[i] + 1) % n; p != locdof[i]; p = (p + 1) % n) {
/* Subtract known values from RHS. */
r[p] -= S[p + locdof[i]*n] * dofval[i];
/* Eliminate DOF (zero row/col "locdof[i]"; leave diagonal). */
S[p + locdof[i]*n] = 0.0;
S[locdof[i] + p*n] = 0.0;
}
/* Produce trivial equation S(i,i)*x(i) = S(i,i)*x0(i). */
r[p] = S[p * (n + 1)] * dofval[i];
}
}
/* ---------------------------------------------------------------------- */
/* Impose boundary conditions on local contribution to global system. */
/* ---------------------------------------------------------------------- */
static int
ifsh_impose_bc(int nconn, int *conn, flowbc_t *bc,
struct ifsh_impl *pimpl)
/* ---------------------------------------------------------------------- */
{
int i, npp, f;
npp = 0;
for (i = 0; i < nconn; i++) {
f = conn[i];
if (bc->type[f] == PRESSURE) {
pimpl->work [npp] = bc->bcval[f];
pimpl->iwork[npp] = i;
npp += 1;
} else if (bc->type[f] == FLUX) {
pimpl->sys->r[i] -= bc->bcval[f];
wsys_WI [i] = WI [perf];
wsys_wdp[i] = wdp[perf];
}
}
if (npp > 0) {
ifsh_explicit_elimination(nconn, npp,
pimpl->iwork,
pimpl->work,
pimpl->sys->S,
pimpl->sys->r);
}
return npp;
}
/* ---------------------------------------------------------------------- */
static void
ifsh_set_effective_params(const double *Binv,
const double *gpress,
const double *totmob,
const double *omega,
struct ifsh_data *ifsh)
/* ---------------------------------------------------------------------- */
{
int c, n, nc, p1, p2, i;
int *pgconn, *gconn;
double *BI, *gp;
nc = ifsh->pimpl->nc;
pgconn = ifsh->pimpl->pgconn;
gconn = ifsh->pimpl->gconn;
BI = ifsh->pimpl->Binv;
gp = ifsh->pimpl->gpress;
p1 = p2 = 0;
for (c = 0; c < nc; c++) {
n = pgconn[c + 1] - pgconn[c];
for (i = 0; i < n ; i++) {
gp[p1 + i] = omega[c] * gpress[p1 + i];
}
for (i = 0; i < n * n; i++) {
BI[p2 + i] = totmob[c] * Binv[p2 + i];
}
p1 += n;
p2 += n * n;
}
}
/* ---------------------------------------------------------------------- */
static int
ifsh_assemble_grid(flowbc_t *bc,
const double *src,
struct ifsh_data *ifsh)
ifsh_assemble_grid(flowbc_t *bc,
const double *Binv,
const double *gpress,
const double *src,
struct fsh_data *ifsh)
/* ---------------------------------------------------------------------- */
{
int c, n, nc, p1, p2;
int npp;
int *pgconn, *gconn;
double *BI, *gp;
nc = ifsh->pimpl->nc;
pgconn = ifsh->pimpl->pgconn;
gconn = ifsh->pimpl->gconn;
BI = ifsh->pimpl->Binv;
gp = ifsh->pimpl->gpress;
pgconn = ifsh->pimpl->gdof_pos;
gconn = ifsh->pimpl->gdof;
p1 = p2 = npp = 0;
for (c = 0; c < nc; c++) {
n = pgconn[c + 1] - pgconn[c];
hybsys_cellcontrib_symm(c, n, p1, p2, gp, src, BI,
hybsys_cellcontrib_symm(c, n, p1, p2, gpress, src, Binv,
ifsh->pimpl->sys);
npp += ifsh_impose_bc(n, gconn + p1, bc, ifsh->pimpl);
npp += fsh_impose_bc(n, gconn + p1, bc, ifsh->pimpl);
hybsys_global_assemble_cell(n, gconn + p1,
ifsh->pimpl->sys->S,
@ -409,40 +144,10 @@ ifsh_assemble_grid(flowbc_t *bc,
/* ---------------------------------------------------------------------- */
static void
ifsh_set_effective_well_params(const double *WI,
const double *wdp,
const double *totmob,
const double *omega,
struct ifsh_data *ifsh)
/* ---------------------------------------------------------------------- */
{
int c, nc, i, perf;
int *cwpos, *cwells;
double *wsys_WI, *wsys_wdp;
nc = ifsh->pimpl->nc;
cwpos = ifsh->pimpl->cwpos;
cwells = ifsh->pimpl->cwells;
wsys_WI = ifsh->pimpl->WI;
wsys_wdp = ifsh->pimpl->wdp;
for (c = i = 0; c < nc; c++) {
for (; i < cwpos[c + 1]; i++) {
perf = cwells[2*i + 1];
wsys_WI [i] = totmob[c] * WI [perf];
wsys_wdp[i] = omega [c] * wdp[perf];
}
}
}
/* ---------------------------------------------------------------------- */
static void
ifsh_impose_well_control(int c,
flowbc_t *bc,
well_control_t *wctrl,
struct ifsh_data *ifsh)
ifsh_impose_well_control(int c,
flowbc_t *bc,
well_control_t *wctrl,
struct fsh_data *ifsh)
/* ---------------------------------------------------------------------- */
{
int ngconn, nwconn, i, w1, w2, wg, f;
@ -454,10 +159,10 @@ ifsh_impose_well_control(int c,
/* Enforce symmetric system */
assert (ifsh->pimpl->wsys->r2w == ifsh->pimpl->wsys->w2r);
pgconn = ifsh->pimpl->pgconn;
pwconn = ifsh->pimpl->cwpos;
pgconn = ifsh->pimpl->gdof_pos;
pwconn = ifsh->pimpl->cwell_pos;
gconn = ifsh->pimpl->gconn + pgconn[c];
gconn = ifsh->pimpl->gdof + pgconn[c];
wconn = ifsh->pimpl->cwells + 2*pwconn[c];
ngconn = pgconn[c + 1] - pgconn[c];
@ -517,9 +222,9 @@ ifsh_impose_well_control(int c,
/* ---------------------------------------------------------------------- */
static int
ifsh_assemble_well(flowbc_t *bc,
well_control_t *wctrl,
struct ifsh_data *ifsh)
ifsh_assemble_well(flowbc_t *bc,
well_control_t *wctrl,
struct fsh_data *ifsh)
/* ---------------------------------------------------------------------- */
{
int npp;
@ -529,9 +234,9 @@ ifsh_assemble_well(flowbc_t *bc,
nc = ifsh->pimpl->nc;
pgconn = ifsh->pimpl->pgconn;
gconn = ifsh->pimpl->gconn;
pwconn = ifsh->pimpl->cwpos;
pgconn = ifsh->pimpl->gdof_pos;
gconn = ifsh->pimpl->gdof;
pwconn = ifsh->pimpl->cwell_pos;
wconn = ifsh->pimpl->cwells;
for (c = 0; c < nc; c++) {
@ -575,6 +280,105 @@ ifsh_assemble_well(flowbc_t *bc,
}
/* ======================================================================
* Public routines follow.
* ====================================================================== */
/* ---------------------------------------------------------------------- */
/* Allocate and define supporting structures for assembling the global
* system of linear equations to couple the grid (reservoir)
* connections represented by 'G' and, if present (i.e., non-NULL),
* the well connections represented by 'W'. */
/* ---------------------------------------------------------------------- */
struct fsh_data *
ifsh_construct(grid_t *G, well_t *W)
/* ---------------------------------------------------------------------- */
{
int nc, ngconn_tot;
size_t idata_sz, ddata_sz, nnu;
struct fsh_data *new;
assert (G != NULL);
/* Allocate master structure, define system matrix sparsity */
new = malloc(1 * sizeof *new);
if (new != NULL) {
new->A = hybsys_define_globconn(G, W);
new->pimpl = NULL;
if (new->A == NULL) {
fsh_destroy(new);
new = NULL;
}
}
/* Allocate implementation structure */
if (new != NULL) {
fsh_count_grid_dof(G, &new->max_ngconn, &new->sum_ngconn2);
ifsh_compute_table_sz(G, W, new->max_ngconn,
&nnu, &idata_sz, &ddata_sz);
new->pimpl = fsh_impl_allocate_basic(idata_sz, ddata_sz);
if (new->pimpl == NULL) {
fsh_destroy(new);
new = NULL;
}
}
/* Allocate Schur complement contributions. Symmetric system. */
if (new != NULL) {
nc = G->number_of_cells;
ngconn_tot = G->cell_facepos[nc];
fsh_define_linsys_arrays(new);
fsh_define_impl_arrays(nc, nnu, ngconn_tot, new->max_ngconn,
W, new->pimpl);
new->pimpl->sys = hybsys_allocate_symm(new->max_ngconn,
nc, ngconn_tot);
if (W != NULL) {
fsh_define_cell_wells(nc, W, new->pimpl);
new->pimpl->wsys = hybsys_well_allocate_symm(new->max_ngconn, nc,
new->pimpl->cwell_pos);
}
if ((new->pimpl->sys == NULL) ||
((W != NULL) && (new->pimpl->wsys == NULL))) {
/* Failed to allocate ->sys or ->wsys (if W != NULL) */
fsh_destroy(new);
new = NULL;
}
}
if (new != NULL) {
/* All allocations succeded. Fill metadata and return. */
new->pimpl->nc = nc;
new->pimpl->nf = G->number_of_faces;
new->pimpl->nw = (W != NULL) ? W->number_of_wells : 0;
memcpy(new->pimpl->gdof_pos,
G->cell_facepos ,
(nc + 1) * sizeof *new->pimpl->gdof_pos);
memcpy(new->pimpl->gdof ,
G->cell_faces ,
ngconn_tot * sizeof *new->pimpl->gdof);
hybsys_init(new->max_ngconn, new->pimpl->sys);
}
return new;
}
/* ---------------------------------------------------------------------- */
/* Assemble global system of linear equations
*
@ -586,37 +390,32 @@ ifsh_assemble_well(flowbc_t *bc,
/* ---------------------------------------------------------------------- */
void
ifsh_assemble(flowbc_t *bc,
double *src,
double *Binv,
double *gpress,
const double *src,
const double *Binv,
const double *gpress,
well_control_t *wctrl,
double *WI,
double *wdp,
double *totmob, /* \sum_i \lambda_i */
double *omega, /* \sum_i \rho_i f_i */
struct ifsh_data *ifsh)
const double *WI,
const double *wdp,
struct fsh_data *ifsh)
/* ---------------------------------------------------------------------- */
{
int npp; /* Number of prescribed pressure values */
ifsh_set_effective_params(Binv, gpress, totmob, omega, ifsh);
hybsys_schur_comp_symm(ifsh->pimpl->nc,
ifsh->pimpl->pgconn,
ifsh->pimpl->Binv,
ifsh->pimpl->sys);
ifsh->pimpl->gdof_pos,
Binv, ifsh->pimpl->sys);
if (ifsh->pimpl->nw > 0) {
ifsh_set_effective_well_params(WI, wdp, totmob, omega, ifsh);
ifsh_set_effective_well_params(WI, wdp, ifsh);
hybsys_well_schur_comp_symm(ifsh->pimpl->nc,
ifsh->pimpl->cwpos,
ifsh->pimpl->cwell_pos,
ifsh->pimpl->WI,
ifsh->pimpl->sys,
ifsh->pimpl->wsys);
}
npp = ifsh_assemble_grid(bc, src, ifsh);
npp = ifsh_assemble_grid(bc, Binv, gpress, src, ifsh);
if (ifsh->pimpl->nw > 0) {
npp += ifsh_assemble_well(bc, wctrl, ifsh);
@ -634,7 +433,9 @@ ifsh_assemble(flowbc_t *bc,
* substitution process, projected half-contact fluxes. */
/* ---------------------------------------------------------------------- */
void
ifsh_press_flux(grid_t *G, struct ifsh_data *h,
ifsh_press_flux(grid_t *G,
const double *Binv, const double *gpress,
struct fsh_data *h,
double *cpress, double *fflux,
double *wpress, double *wflux)
/* ---------------------------------------------------------------------- */
@ -645,8 +446,7 @@ ifsh_press_flux(grid_t *G, struct ifsh_data *h,
hybsys_compute_press_flux(G->number_of_cells,
G->cell_facepos,
G->cell_faces,
h->pimpl->gpress,
h->pimpl->Binv,
gpress, Binv,
h->pimpl->sys,
h->x, cpress, h->pimpl->cflux,
h->pimpl->work);
@ -655,8 +455,8 @@ ifsh_press_flux(grid_t *G, struct ifsh_data *h,
assert ((wpress != NULL) && (wflux != NULL));
hybsys_compute_press_flux_well(G->number_of_cells, G->cell_facepos,
G->number_of_faces, h->pimpl->nw,
h->pimpl->cwpos, h->pimpl->cwells,
h->pimpl->Binv, h->pimpl->WI,
h->pimpl->cwell_pos, h->pimpl->cwells,
Binv, h->pimpl->WI,
h->pimpl->wdp, h->pimpl->sys,
h->pimpl->wsys, h->x, cpress,
h->pimpl->cflux, wpress, wflux,

69
ifsh.h
View File

@ -17,14 +17,12 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef IFSH_H_INCLUDED
#define IFHS_H_INCLUDED
#ifndef OPM_IFSH_HEADER_INCLUDED
#define OPM_IFHS_HEADER_INCLUDED
#include "grid.h"
#include "well.h"
#include "flow_bc.h"
#include "sparse_sys.h"
#ifdef __cplusplus
extern "C" {
@ -38,48 +36,19 @@ extern "C" {
* and face fluxes is also included.
*/
struct fsh_data;
struct ifsh_impl;
/** Contains the linear system for assembly, as well as internal data
* for the assembly routines.
*/
struct ifsh_data {
/* Let \f$n_i\f$ be the number of connections/faces of grid cell
* number \f$i\f$. Then max_ngconn = \f$\max_i n_i\f$
*/
int max_ngconn;
/* With n_i as above, sum_ngconn2 = \f$\sum_i n_i^2\f$ */
size_t sum_ngconn2;
/* Linear system */
struct CSRMatrix *A; /* Coefficient matrix */
double *b; /* System RHS */
double *x; /* Solution */
/* Private implementational details. */
struct ifsh_impl *pimpl;
};
/** Constructs the ifsh_data object for a given grid and well
* pattern.
/** Constructs incompressible hybrid flow-solver data object for a
* given grid and well pattern.
*
* @param G The grid
* @param W The wells
*/
struct ifsh_data *
struct fsh_data *
ifsh_construct(grid_t *G, well_t *W);
/** Destroys the ifsh_data object */
void
ifsh_destroy(struct ifsh_data *h);
/** Assembles the hybridized linear system for face pressures.
*
* This routine produces no output, other than changing the linear
@ -107,16 +76,14 @@ ifsh_destroy(struct ifsh_data *h);
* be modified). Must already be constructed.
*/
void
ifsh_assemble(flowbc_t *bc,
double *src,
double *Binv,
double *gpress,
well_control_t *wctrl,
double *WI,
double *wdp,
double *totmob, /* \sum_i \lambda_i */
double *omega, /* \sum_i \rho_i f_i */
struct ifsh_data *h);
ifsh_assemble(flowbc_t *bc,
const double *src,
const double *Binv,
const double *gpress,
well_control_t *wctrl,
const double *WI,
const double *wdp,
struct fsh_data *h);
/** Computes cell pressures, face fluxes, well pressures and well
* fluxes from face pressures.
@ -131,7 +98,9 @@ ifsh_assemble(flowbc_t *bc,
* @param wflux[out] \TODO
*/
void
ifsh_press_flux(grid_t *G, struct ifsh_data *h,
ifsh_press_flux(grid_t *G,
const double *Binv, const double *gpress,
struct fsh_data *h,
double *cpress, double *fflux,
double *wpress, double *wflux);
@ -141,4 +110,4 @@ ifsh_press_flux(grid_t *G, struct ifsh_data *h,
#endif
#endif /* IFSH_H_INCLUDED */
#endif /* OPM_IFSH_HEADER_INCLUDED */

11
ifsh_ms.c
View File

@ -348,17 +348,6 @@ ifsh_ms_impl_construct(grid_t *G ,
}
/* ---------------------------------------------------------------------- */
static void
vector_zero(size_t m, double *v)
/* ---------------------------------------------------------------------- */
{
size_t i;
for (i = 0; i < m; i++) { v[i] = 0.0; }
}
/* ---------------------------------------------------------------------- */
static void
average_flux(size_t nf, const int *N, double *flux) /* Poor name */

6
ifsh_ms.h
View File

@ -17,8 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef IFSH_MS_H_INCLUDED
#define IFSH_MS_H_INCLUDED
#ifndef OPM_IFSH_MS_HEADER_INCLUDED
#define OPM_IFSH_MS_HEADER_INCLUDED
#ifdef __cplusplus
extern "C" {
@ -68,4 +68,4 @@ ifsh_ms_press_flux(grid_t *G, struct ifsh_ms_data *h,
}
#endif
#endif /* IFSH_MS_H_INCLUDED */
#endif /* OPM_IFSH_MS_HEADER_INCLUDED */

65
mimetic.c
View File

@ -183,3 +183,68 @@ mim_ip_linpress_exact(int nf, int nconn, int d,
dgemm_("No Transpose", "Transpose", &m, &m, &n,
&a1, work, &ld1, N, &ld1, &a2, Binv, &ldBinv);
}
/* ---------------------------------------------------------------------- */
void
mim_ip_compute_gpress(int nc, int d, const double *grav,
const int *pconn, const int *conn,
const double *fcentroid, const double *ccentroid,
double *gpress)
/* ---------------------------------------------------------------------- */
{
int c, i, j;
const double *cc, *fc;
for (c = i = 0; c < nc; c++) {
cc = ccentroid + (c * d);
for (; i < pconn[c + 1]; i++) {
fc = fcentroid + (conn[i] * d);
gpress[i] = 0.0;
for (j = 0; j < d; j++) {
gpress[i] += grav[j] * (fc[j] - cc[j]);
}
}
}
}
/* inv(B) <- \lambda_t(s)*inv(B) */
/* ---------------------------------------------------------------------- */
void
mim_ip_mobility_update(int nc, const int *pconn, const double *totmob,
double *Binv)
/* ---------------------------------------------------------------------- */
{
int c, i, n, p2;
for (c = p2 = 0; c < nc; c++) {
n = pconn[c + 1] - pconn[c];
for (i = 0; i < n * n; i++) {
Binv[p2 + i] *= totmob[c];
}
p2 += n * n;
}
}
/* G <- \sum_i \rho_i f_i(s) * G */
/* ---------------------------------------------------------------------- */
void
mim_ip_density_update(int nc, const int *pconn, const double *omega,
double *gpress)
/* ---------------------------------------------------------------------- */
{
int c, i;
for (c = i = 0; c < nc; c++) {
for (; i < pconn[c + 1]; i++) {
gpress[i] *= omega[c];
}
}
}

22
mimetic.h
View File

@ -17,8 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef MIMETIC_H_INCLUDED
#define MIMETIC_H_INCLUDED
#ifndef OPM_MIMETIC_HEADER_INCLUDED
#define OPM_MIMETIC_HEADER_INCLUDED
#ifdef __cplusplus
extern "C" {
@ -81,8 +81,24 @@ void mim_ip_simple_all(int ncells, int d, int max_ncf, int *ncf,
double *farea, double *ccentroid, double *cvol,
double *perm, double *Binv);
void
mim_ip_compute_gpress(int nc, int d, const double *grav,
const int *pconn, const int *conn,
const double *fcentroid, const double *ccentroid,
double *gpress);
/* inv(B) <- \lambda_t(s)*inv(B) */
void
mim_ip_mobility_update(int nc, const int *pconn, const double *totmob,
double *Binv);
/* G <- \sum_i \rho_i f_i(s) * G */
void
mim_ip_density_update(int nc, const int *pconn, const double *omega,
double *gpress);
#ifdef __cplusplus
}
#endif
#endif /* MIMETIC_H_INCLUDED */
#endif /* OPM_MIMETIC_HEADER_INCLUDED */

4
partition.h
View File

@ -17,8 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef PARTITION_H_INCLUDED
#define PARTITION_H_INCLUDED
#ifndef OPM_PARTITION_HEADER_INCLUDED
#define OPM_PARTITION_HEADER_INCLUDED
#ifdef __cplusplus
extern "C" {

13
sparse_sys.c
View File

@ -192,3 +192,16 @@ csrmatrix_zero(struct CSRMatrix *A)
for (i = 0; i < A->nnz; i++) { A->sa[i] = 0.0; }
}
/* ---------------------------------------------------------------------- */
/* v = zeros([n, 1]) */
/* ---------------------------------------------------------------------- */
void
vector_zero(size_t n, double *v)
/* ---------------------------------------------------------------------- */
{
size_t i;
for (i = 0; i < n; i++) { v[i] = 0.0; }
}

12
sparse_sys.h
View File

@ -17,8 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef SPARSE_SYS_H_INCLUDED
#define SPARSE_SYS_H_INCLUDED
#ifndef OPM_SPARSE_SYS_HEADER_INCLUDED
#define OPM_SPARSE_SYS_HEADER_INCLUDED
#include <stddef.h>
@ -72,8 +72,14 @@ csrmatrix_delete(struct CSRMatrix *A);
void
csrmatrix_zero(struct CSRMatrix *A);
/* ---------------------------------------------------------------------- */
/* v = zeros([n, 1]) */
/* ---------------------------------------------------------------------- */
void
vector_zero(size_t n, double *v);
#ifdef __cplusplus
}
#endif
#endif /* SPARSE_SYS_H_INCLUDED */
#endif /* OPM_SPARSE_SYS_HEADER_INCLUDED */

110
trans_tpfa.c Normal file
View File

@ -0,0 +1,110 @@
#include <assert.h>
#include <math.h>
#include "blas_lapack.h"
#include "trans_tpfa.h"
/* ---------------------------------------------------------------------- */
/* htrans <- sum(C(:,i) .* K(cellNo,:) .* N(:,j), 2) ./ sum(C.*C, 2) */
/* ---------------------------------------------------------------------- */
void
tpfa_htrans_compute(grid_t *G, const double *perm, double *htrans)
/* ---------------------------------------------------------------------- */
{
int c, d, f, i, j;
double s, dist, denom;
double Kn[3];
double *cc, *fc, *n;
double *K;
MAT_SIZE_T nrows, ncols, ldA, incx, incy;
double a1, a2;
d = G->dimensions;
nrows = ncols = ldA = d;
incx = incy = 1 ;
a1 = 1.0; a2 = 0.0 ;
for (c = i = 0; c < G->number_of_cells; c++) {
K = perm + (c * d * d);
cc = G->cell_centroids + (c * d);
for (; i < G->cell_facepos[c + 1]; i++) {
f = G->cell_faces[i];
s = 2.0*(G->face_cells[2*f + 0] == c) - 1.0;
n = G->face_normals + (f * d);
fc = G->face_centroids + (f * d);
dgemv_("No Transpose", &nrows, &ncols,
&a1, K, &ldA, n, &incx, &a2, &Kn[0], &incy);
htrans[i] = denom = 0.0;
for (j = 0; j < d; j++) {
dist = fc[j] - cc[j];
htrans[i] += s * dist * Kn[j];
denom += dist * dist;
}
assert (denom > 0);
htrans[i] /= denom;
htrans[i] = fabs(htrans[i]);
}
}
}
/* ---------------------------------------------------------------------- */
void
tpfa_trans_compute(grid_t *G, const double *htrans, double *trans)
/* ---------------------------------------------------------------------- */
{
int c, i, f;
for (f = 0; f < G->number_of_faces; f++) {
trans[f] = 0.0;
}
for (c = i = 0; c < G->number_of_cells; c++) {
for (; i < G->cell_facepos[c + 1]; i++) {
f = G->cell_faces[i];
trans[f] += 1.0 / htrans[i];
}
}
for (f = 0; f < G->number_of_faces; f++) {
trans[f] = 1.0 / trans[f];
}
}
/* ---------------------------------------------------------------------- */
void
tpfa_eff_trans_compute(grid_t *G,
const double *totmob,
const double *htrans,
double *trans)
/* ---------------------------------------------------------------------- */
{
int c, i, f;
for (f = 0; f < G->number_of_faces; f++) {
trans[f] = 0.0;
}
for (c = i = 0; c < G->number_of_cells; c++) {
for (; i < G->cell_facepos[c + 1]; i++) {
f = G->cell_faces[i];
trans[f] += 1.0 / (totmob[c] * htrans[i]);
}
}
for (f = 0; f < G->number_of_faces; f++) {
trans[f] = 1.0 / trans[f];
}
}

37
trans_tpfa.h Normal file
View File

@ -0,0 +1,37 @@
/*
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_TRANS_TPFA_HEADER_INCLUDED
#define OPM_TRANS_TPFA_HEADER_INCLUDED
#include "grid.h"
void
tpfa_htrans_compute(grid_t *G, const double *perm, double *htrans);
void
tpfa_trans_compute(grid_t *G, const double *htrans, double *trans);
void
tpfa_eff_trans_compute(grid_t *G,
const double *totmob,
const double *htrans,
double *trans);
#endif /* OPM_TRANS_TPFA_HEADER_INCLUDED */

6
well.h
View File

@ -17,8 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef WELL_H_INCLUDED
#define WELL_H_INCLUDED
#ifndef OPM_WELL_HEADER_INCLUDED
#define OPM_WELL_HEADER_INCLUDED
#ifdef __cplusplus
extern "C" {
@ -54,4 +54,4 @@ derive_cell_wells(int nc, well_t *W, int *cwpos, int *cwells);
}
#endif
#endif /* WELL_H_INCLUDED */
#endif /* OPM_WELL_HEADER_INCLUDED */