Merged changes

This commit is contained in:
Ingeborg Ligaarden 2011-11-25 14:47:01 +01:00
commit 8595ea52b0
14 changed files with 1980 additions and 17 deletions

View File

@ -7,8 +7,12 @@ libopmpressure_HEADERS = \
blas_lapack.h \
coarse_conn.h \
coarse_sys.h \
compr_bc.h \
compr_quant.h \
compr_quant_general.h \
compr_source.h \
cfs_tpfa.h \
cfs_tpfa_residual.h \
dfs.h \
flow_bc.h \
fsh.h \
@ -32,8 +36,12 @@ libopmpressure_la_SOURCES = \
cfsh.c \
coarse_conn.c \
coarse_sys.c \
compr_bc.c \
compr_quant.c \
compr_quant_general.c \
compr_source.c \
cfs_tpfa.c \
cfs_tpfa_residual.c \
dfs.c \
flow_bc.c \
fsh.c \

View File

@ -70,7 +70,7 @@ void dgetrs_(const char *trans, const MAT_SIZE_T *n,
const double *A , const MAT_SIZE_T *lda,
const MAT_SIZE_T *ipiv , double *B,
const MAT_SIZE_T *ldb , MAT_SIZE_T *info);
/* A <- chol(A) */
void dpotrf_(const char *uplo, const MAT_SIZE_T *n,
double *A , const MAT_SIZE_T *lda,

1190
src/cfs_tpfa_residual.c Normal file

File diff suppressed because it is too large Load Diff

140
src/cfs_tpfa_residual.h Normal file
View File

@ -0,0 +1,140 @@
/*
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_CFS_TPFA_HEADER_INCLUDED
#define OPM_CFS_TPFA_HEADER_INCLUDED
#include "grid.h"
#ifdef __cplusplus
extern "C" {
#endif
struct cfs_tpfa_res_impl;
struct CSRMatrix;
struct compr_quantities_gen;
struct compr_src;
struct compr_bc;
struct WellCompletions;
struct WellControls;
struct completion_data;
struct cfs_tpfa_res_wells {
struct WellCompletions *conn;
struct WellControls *ctrl;
struct completion_data *data;
};
struct cfs_tpfa_res_forces {
struct cfs_tpfa_res_wells *W ;
struct compr_bc *bc ;
struct compr_src *src;
};
struct cfs_tpfa_res_data {
struct CSRMatrix *J;
double *F;
struct cfs_tpfa_res_impl *pimpl;
};
struct cfs_tpfa_res_data *
cfs_tpfa_res_construct(grid_t *G ,
struct WellCompletions *wconn ,
int nphases);
void
cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h);
void
cfs_tpfa_res_assemble(grid_t *G,
double dt,
struct cfs_tpfa_res_forces *forces,
const double *zc,
struct compr_quantities_gen *cq,
const double *trans,
const double *gravcap_f,
const double *cpress,
const double *wpress,
const double *porevol,
struct cfs_tpfa_res_data *h);
void
cfs_tpfa_res_flux(grid_t *G,
flowbc_t *bc,
int np,
const double *trans,
const double *pmobf,
const double *gravcap_f,
const double *cpress,
double *fflux);
void
cfs_tpfa_res_fpress(grid_t *G,
flowbc_t *bc,
int np,
const double *htrans,
const double *pmobf,
const double *gravcap_f,
struct cfs_tpfa_res_data *h,
const double *cpress,
const double *fflux,
double *fpress);
#if 0
void
cfs_tpfa_retrieve_masstrans(grid_t *G,
int np,
struct cfs_tpfa_data *h,
double *masstrans_f);
void
cfs_tpfa_retrieve_gravtrans(grid_t *G,
int np,
struct cfs_tpfa_data *h,
double *gravtrans_f);
double
cfs_tpfa_impes_maxtime(grid_t *G,
struct compr_quantities *cq,
const double *trans,
const double *porevol,
struct cfs_tpfa_data *h,
const double *dpmobf,
const double *surf_dens,
const double *gravity);
void
cfs_tpfa_expl_mass_transport(grid_t *G,
well_t *W,
struct completion_data *wdata,
int np,
double dt,
const double *porevol,
struct cfs_tpfa_data *h,
double *surf_vol);
#endif
#ifdef __cplusplus
}
#endif
#endif /* OPM_CFS_TPFA_HEADER_INCLUDED */

View File

@ -975,13 +975,6 @@ unenumerate_local_dofs(size_t cf,
{
int *b, *c, i;
/* Note that memset()-ing all of m->loc_fno (g->number_of_faces
* entries) may or may not be faster than this loop. The loop is
* entirely random access to m->loc_fno which ruins locality, but
* typically visits only a small subset of the actual elements.
*
* Profiling/measurement needed. */
for (b = ct->neighbours + 2*(cf + 0);
b != ct->neighbours + 2*(cf + 1); b++) {
if (*b >= 0) {

182
src/compr_bc.c Normal file
View File

@ -0,0 +1,182 @@
/*===========================================================================
//
// File: compr_bc.c
//
// Created: 2011-10-24 16:07:17+0200
//
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
Copyright 2011 SINTEF ICT, Applied Mathematics.
Copyright 2011 Statoil ASA.
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 <stdlib.h>
#include <string.h>
#include "compr_bc.h"
static int
expand_source_tables(int alloc, struct compr_bc *bc)
{
enum compr_bc_type *t;
int *f;
double *p, *v, *s;
t = realloc(bc->type , alloc * 1 * sizeof *t);
f = realloc(bc->face , alloc * 1 * sizeof *f);
p = realloc(bc->press , alloc * 1 * sizeof *p);
v = realloc(bc->flux , alloc * 1 * sizeof *v);
s = realloc(bc->saturation, alloc * bc->nphases * sizeof *s);
if ((t == NULL) || (f == NULL) ||
(p == NULL) || (v == NULL) || (s == NULL)) {
free(s); free(v); free(p); free(f); free(t);
alloc = 0;
} else {
bc->type = t; bc->face = f;
bc->press = p; bc->flux = v;
bc->saturation = s;
}
return alloc;
}
/* ======================================================================
* Public methods below separator.
* ====================================================================== */
/* ---------------------------------------------------------------------- */
struct compr_bc *
compr_bc_allocate(int np, int nbc)
/* ---------------------------------------------------------------------- */
{
int status;
struct compr_bc *bc;
if (np <= 0) {
bc = NULL;
} else {
bc = malloc(1 * sizeof *bc);
if (bc != NULL) {
bc->nbc = 0 ;
bc->cpty = 0 ;
bc->nphases = np;
bc->type = NULL;
bc->face = NULL;
bc->press = NULL;
bc->flux = NULL;
bc->saturation = NULL;
if (nbc > 0) {
status = expand_source_tables(nbc, bc);
} else {
status = 1;
}
if (status <= 0) {
compr_bc_deallocate(bc);
bc = NULL;
}
}
}
return bc;
}
/* ---------------------------------------------------------------------- */
void
compr_bc_deallocate(struct compr_bc *bc)
/* ---------------------------------------------------------------------- */
{
if (bc != NULL) {
free(bc->saturation);
free(bc->flux );
free(bc->press );
free(bc->face );
free(bc->type );
}
free(bc);
}
/* ---------------------------------------------------------------------- */
int
compr_bc_append(enum compr_bc_type t ,
int f ,
int np ,
double p ,
double v ,
const double *sat,
struct compr_bc *bc )
/* ---------------------------------------------------------------------- */
{
int alloc, status;
if (np != bc->nphases) {
return -1;
}
if (bc->nbc == bc->cpty) {
alloc = (bc->cpty > 0) ? 2 * bc->cpty : 1;
status = expand_source_tables(alloc, bc);
} else {
assert (bc->nbc < bc->cpty);
status = 1;
}
if (status > 0) {
bc->type [ bc->nbc ] = t;
bc->face [ bc->nbc ] = f;
bc->press[ bc->nbc ] = p;
bc->flux [ bc->nbc ] = v;
memcpy(bc->saturation + (bc->nbc * np), sat,
np * sizeof *bc->saturation);
bc->nbc += 1;
}
return status > 0;
}
/* ---------------------------------------------------------------------- */
void
compr_bc_clear(struct compr_bc *bc)
/* ---------------------------------------------------------------------- */
{
bc->nbc = 0;
}

83
src/compr_bc.h Normal file
View File

@ -0,0 +1,83 @@
/*===========================================================================
//
// File: compr_bc.h
//
// Created: 2011-10-24 15:58:16+0200
//
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
Copyright 2011 SINTEF ICT, Applied Mathematics.
Copyright 2011 Statoil ASA.
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_BC_H_HEADER
#define OPM_COMPR_BC_H_HEADER
#ifdef __cplusplus
extern "C" {
#endif
enum compr_bc_type { BC_PRESSURE, BC_RESV };
struct compr_bc {
int nbc;
int cpty;
int nphases;
enum compr_bc_type *type;
int *face;
double *press;
double *flux;
double *saturation;
};
struct compr_bc *
compr_bc_allocate(int np, int nbc);
void
compr_bc_deallocate(struct compr_bc *bc);
int
compr_bc_append(enum compr_bc_type type,
int f ,
int np ,
double p ,
double v ,
const double *sat ,
struct compr_bc *bc );
void
compr_bc_clear(struct compr_bc *bc);
#ifdef __cplusplus
}
#endif
#endif /* OPM_COMPR_BC_H_HEADER */

72
src/compr_quant_general.c Normal file
View File

@ -0,0 +1,72 @@
/*
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 <stdlib.h>
#include "sparse_sys.h"
#include "compr_quant_general.h"
void
compr_quantities_gen_deallocate(struct compr_quantities_gen *cq)
{
if (cq != NULL) {
free(cq->Ac);
}
free(cq);
}
struct compr_quantities_gen *
compr_quantities_gen_allocate(size_t nc, size_t nf, int np)
{
size_t alloc_sz, np2;
struct compr_quantities_gen *cq;
cq = malloc(1 * sizeof *cq);
if (cq != NULL) {
np2 = np * np;
alloc_sz = np2 * nc; /* Ac */
alloc_sz += np2 * nc; /* dAc */
alloc_sz += np2 * nf; /* Af */
alloc_sz += np * nf; /* phasemobf */
alloc_sz += 1 * nc; /* voldiscr */
cq->Ac = malloc(alloc_sz * sizeof *cq->Ac);
if (cq->Ac == NULL) {
compr_quantities_gen_deallocate(cq);
cq = NULL;
} else {
cq->dAc = cq->Ac + (np2 * nc);
cq->Af = cq->dAc + (np2 * nc);
cq->phasemobf = cq->Af + (np2 * nf);
cq->voldiscr = cq->phasemobf + (np * nf);
cq->nphases = np;
vector_zero(alloc_sz, cq->Ac);
}
}
return cq;
}

49
src/compr_quant_general.h Normal file
View File

@ -0,0 +1,49 @@
/*
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>
#ifdef __cplusplus
extern "C" {
#endif
struct compr_quantities_gen {
int nphases; /* Number of phases/components */
double *Ac; /* RB^{-1} per cell */
double *dAc; /* d/dp (RB^{-1}) per cell */
double *Af; /* RB^{-1} per face */
double *phasemobf; /* Phase mobility per face */
double *voldiscr; /* Volume discrepancy per cell */
};
struct compr_quantities_gen *
compr_quantities_gen_allocate(size_t nc, size_t nf, int np);
void
compr_quantities_gen_deallocate(struct compr_quantities_gen *cq);
#ifdef __cplusplus
}
#endif
#endif /* OPM_COMPR_QUANT_HEADER_INCLUDED */

169
src/compr_source.c Normal file
View File

@ -0,0 +1,169 @@
/*===========================================================================
//
// File: compr_source.c
//
// Created: 2011-10-19 19:22:24+0200
//
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
Copyright 2011 SINTEF ICT, Applied Mathematics.
Copyright 2011 Statoil ASA.
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 <stdlib.h>
#include <string.h>
#include "compr_source.h"
static int
expand_source_tables(int alloc, struct compr_src *src)
{
int *c;
double *v, *s;
c = realloc(src->cell , alloc * 1 * sizeof *c);
v = realloc(src->flux , alloc * 1 * sizeof *v);
s = realloc(src->saturation, alloc * src->nphases * sizeof *s);
if ((c == NULL) || (v == NULL) || (s == NULL)) {
free(s); free(v); free(c);
alloc = 0;
} else {
src->cell = c; src->cpty = alloc;
src->flux = v; src->saturation = s;
}
return alloc;
}
/* ======================================================================
* Public methods below separator.
* ====================================================================== */
/* ---------------------------------------------------------------------- */
struct compr_src *
compr_src_allocate(int np, int nsrc)
/* ---------------------------------------------------------------------- */
{
int status;
struct compr_src *src;
if (np <= 0) {
src = NULL;
} else {
src = malloc(1 * sizeof *src);
if (src != NULL) {
src->nsrc = 0 ;
src->cpty = 0 ;
src->nphases = np;
src->cell = NULL;
src->flux = NULL;
src->saturation = NULL;
if (nsrc > 0) {
status = expand_source_tables(nsrc, src);
} else {
status = 1;
}
if (status <= 0) {
compr_src_deallocate(src);
src = NULL;
}
}
}
return src;
}
/* ---------------------------------------------------------------------- */
void
compr_src_deallocate(struct compr_src *src)
/* ---------------------------------------------------------------------- */
{
if (src != NULL) {
free(src->saturation);
free(src->flux );
free(src->cell );
}
free(src);
}
/* ---------------------------------------------------------------------- */
int
append_compr_source_term(int c ,
int np ,
double v ,
const double *sat,
struct compr_src *src)
/* ---------------------------------------------------------------------- */
{
int alloc, status;
if (np != src->nphases) {
return -1;
}
if (src->nsrc == src->cpty) {
alloc = (src->cpty > 0) ? 2 * src->cpty : 1;
status = expand_source_tables(alloc, src);
} else {
assert (src->nsrc < src->cpty);
status = 1;
}
if (status > 0) {
src->cell[ src->nsrc ] = c;
src->flux[ src->nsrc ] = v;
memcpy(src->saturation + (src->nsrc * np), sat,
np * sizeof *src->saturation);
src->nsrc += 1;
}
return status > 0;
}
/* ---------------------------------------------------------------------- */
void
clear_compr_source_term(struct compr_src *src)
/* ---------------------------------------------------------------------- */
{
src->nsrc = 0;
}

76
src/compr_source.h Normal file
View File

@ -0,0 +1,76 @@
/*===========================================================================
//
// File: compr_source.h
//
// Created: 2011-10-19 19:14:30+0200
//
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
Copyright 2011 SINTEF ICT, Applied Mathematics.
Copyright 2011 Statoil ASA.
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_SOURCE_H_HEADER
#define OPM_COMPR_SOURCE_H_HEADER
#ifdef __cplusplus
extern "C" {
#endif
struct compr_src {
int nsrc;
int cpty;
int nphases;
int *cell;
double *flux;
double *saturation;
};
struct compr_src *
compr_src_allocate(int np, int nsrc);
void
compr_src_deallocate(struct compr_src *src);
int
append_compr_source_term(int c ,
int np ,
double v ,
const double *sat,
struct compr_src *src);
void
clear_compr_source_term(struct compr_src *src);
#ifdef __cplusplus
}
#endif
#endif /* OPM_COMPR_SOURCE_H_HEADER */

View File

@ -182,9 +182,7 @@ void
csrmatrix_zero(struct CSRMatrix *A)
/* ---------------------------------------------------------------------- */
{
size_t i;
for (i = 0; i < A->nnz; i++) { A->sa[i] = 0.0; }
vector_zero(A->nnz, A->sa);
}
@ -247,7 +245,7 @@ vector_write(size_t n, const double *v, const char *fn)
if (fp != NULL) {
vector_write_stream(n, v, fp);
}
fclose(fp);
}

View File

@ -674,7 +674,7 @@ main(void)
for (i = 0; i < G->number_of_cells; i++) {
fprintf(stderr, "press(%02d) = %g;\n", i + 1, cpress[i]);
}
cfs_tpfa_destroy(h);
deallocate_cq(cq);
deallocate_flowbc(bc);

View File

@ -28,17 +28,17 @@ extern "C" {
enum well_type { INJECTOR, PRODUCER };
enum well_control { BHP , RATE };
typedef struct {
struct WellCompletions {
int number_of_wells;
int *well_connpos;
int *well_cells;
} well_t;
};
typedef struct {
struct WellControls {
enum well_type *type;
enum well_control *ctrl;
double *target;
} well_control_t;
};
struct completion_data {
double *WI; /* Productivity index */
@ -47,6 +47,9 @@ struct completion_data {
double *phasemob; /* Phase mobility */
};
typedef struct WellCompletions well_t;
typedef struct WellControls well_control_t;
int
allocate_cell_wells(int nc, well_t *W, int **cwpos, int **cwells);