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

Merge.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2012-03-06 22:55:59 +01:00
commit 046e9b1f2f
16 changed files with 517 additions and 180 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

53
opm/core/pressure/HybridPressureSolver.hpp
View File

@ -28,7 +28,6 @@
#include <stdexcept>
/// @brief
/// Encapsulates the ifsh (= incompressible flow solver hybrid) solver modules.
class HybridPressureSolver
@ -137,17 +136,11 @@ public:
}
// Boundary conditions.
int num_faces = grid_.c_grid()->number_of_faces;
assert(num_faces == int(bctypes.size()));
std::vector<flowbc_type> bctypes2(num_faces, UNSET);
for (int face = 0; face < num_faces; ++face) {
if (bctypes[face] == FBC_PRESSURE) {
bctypes2[face] = PRESSURE;
} else if (bctypes[face] == FBC_FLUX) {
bctypes2[face] = FLUX;
}
}
flowbc_t bc = { &bctypes2[0], const_cast<double*>(&bcvalues[0]) };
assert (bctypes.size() ==
static_cast<std::vector<FlowBCTypes>::size_type>(grid_.numFaces()));
FlowBoundaryConditions *bc = gather_boundary_conditions(bctypes, bcvalues);
// Source terms from user.
double* src = const_cast<double*>(&sources[0]); // Ugly? Yes. Safe? I think so.
@ -174,9 +167,11 @@ public:
}
// Assemble the embedded linear system.
ifsh_assemble(&bc, src, &Binv_mobilityweighted_[0], &gpress_omegaweighted_[0],
ifsh_assemble(bc, src, &Binv_mobilityweighted_[0], &gpress_omegaweighted_[0],
wctrl, WI, wdp, data_);
state_ = Assembled;
flow_conditions_destroy(bc);
}
/// Encapsulate a sparse linear system in CSR format.
@ -297,7 +292,37 @@ private:
// Gravity contributions.
std::vector<double> gpress_;
std::vector<double> gpress_omegaweighted_;
};
FlowBoundaryConditions*
gather_boundary_conditions(const std::vector<FlowBCTypes>& bctypes ,
const std::vector<double>& bcvalues)
{
FlowBoundaryConditions* fbc = flow_conditions_construct(0);
int ok = fbc != 0;
std::vector<FlowBCTypes>::size_type i;
for (i = 0; ok && (i < bctypes.size()); ++i) {
if (bctypes[ i ] == FBC_PRESSURE) {
ok = flow_conditions_append(BC_PRESSURE,
static_cast<int>(i),
bcvalues[ i ],
fbc);
}
else if (bctypes[ i ] == FBC_FLUX) {
ok = flow_conditions_append(BC_FLUX_TOTVOL,
static_cast<int>(i),
bcvalues[ i ],
fbc);
}
}
return fbc;
}
}; // class HybridPressureSolver
#endif // OPM_HYBRIDPRESSURESOLVER_HEADER_INCLUDED

5
opm/core/pressure/IncompTpfa.cpp
View File

@ -105,7 +105,10 @@ namespace Opm
}
}
ifs_tpfa_assemble(gg, &trans_[0], &src[0], &gpress_omegaweighted_[0], h_);
ifs_tpfa_forces F;
F.src = &src[0];
ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_);
linsolver_.solve(h_->A, h_->b, h_->x);

73
opm/core/pressure/TPFACompressiblePressureSolver.hpp
View File

@ -41,7 +41,7 @@ public:
/// @brief
/// Default constructor, does nothing.
TPFACompressiblePressureSolver()
: state_(Uninitialized), data_(0)
: state_(Uninitialized), data_(0), bc_(0)
{
wells_.number_of_wells = 0;
}
@ -53,6 +53,7 @@ public:
/// Destructor.
~TPFACompressiblePressureSolver()
{
flow_conditions_destroy(bc_);
cfs_tpfa_destroy(data_);
}
@ -132,7 +133,9 @@ public:
/// Boundary condition types.
enum FlowBCTypes { FBC_UNSET = UNSET, FBC_PRESSURE = PRESSURE, FBC_FLUX = FLUX};
enum FlowBCTypes { FBC_UNSET = BC_NOFLOW ,
FBC_PRESSURE = BC_PRESSURE ,
FBC_FLUX = BC_FLUX_TOTVOL };
/// @brief
/// Assemble the sparse system.
@ -167,20 +170,7 @@ public:
UnstructuredGrid* g = grid_.c_grid();
// Boundary conditions.
int num_faces = g->number_of_faces;
// assert(num_faces == int(bctypes.size()));
bctypes_.clear();
bctypes_.resize(num_faces, UNSET);
for (int face = 0; face < num_faces; ++face) {
if (bctypes[face] == FBC_PRESSURE) {
bctypes_[face] = PRESSURE;
} else if (bctypes[face] == FBC_FLUX) {
bctypes_[face] = FLUX;
}
}
bcvalues_.resize(num_faces);
std::copy(bcvalues, bcvalues + num_faces, bcvalues_.begin());
flowbc_t bc = { &bctypes_[0], &bcvalues_[0] };
gather_boundary_conditions(bctypes, bcvalues);
// Source terms from user.
double* src = const_cast<double*>(sources); // Ugly? Yes. Safe? I think so.
@ -208,7 +198,7 @@ public:
const_cast<double *>(phasemobf) };
// Call the assembly routine. After this, linearSystem() may be called.
cfs_tpfa_assemble(g, dt, wells, &bc, src,
cfs_tpfa_assemble(g, dt, wells, bc_, src,
&cq, &trans_[0], gravcapf,
wctrl, wcompl,
cell_pressure, &porevol_[0],
@ -284,9 +274,7 @@ public:
face_pressures.resize(num_faces, 0.0);
face_fluxes.clear();
face_fluxes.resize(num_faces, 0.0);
// ifs_tpfa_press_flux(grid_.c_grid(), &eff_trans_[0],
// data_, &cell_pressures[0], &face_fluxes[0]);
flowbc_t bc = { &bctypes_[0], const_cast<double*>(&bcvalues_[0]) };
int np = 3; // Number of phases.
// Wells.
@ -304,12 +292,12 @@ public:
}
cfs_tpfa_press_flux(grid_.c_grid(),
&bc, wells,
bc_, wells,
np, &trans_[0], &phasemobf_[0], &gravcapf_[0],
wcompl,
data_, &cell_pressures[0], &face_fluxes[0],
wpress, wflux);
cfs_tpfa_fpress(grid_.c_grid(), &bc, np, &htrans_[0],
cfs_tpfa_fpress(grid_.c_grid(), bc_, np, &htrans_[0],
&phasemobf_[0], &gravcapf_[0], data_, &cell_pressures[0],
&face_fluxes[0], &face_pressures[0]);
}
@ -430,8 +418,7 @@ private:
double gravity_[3];
// Boundary conditions.
std::vector<flowbc_type> bctypes_;
std::vector<double> bcvalues_;
FlowBoundaryConditions *bc_;
// Well data
well_t wells_;
@ -496,7 +483,43 @@ private:
wcompl_.phasemob = &well_phasemob_storage_[0];
}
};
void
gather_boundary_conditions(const FlowBCTypes* bctypes ,
const double* bcvalues)
{
if (bc_ == 0) {
bc_ = flow_conditions_construct(0);
}
else {
flow_conditions_clear(bc_);
}
int ok = bc_ != 0;
for (std::size_t i = 0, nf = grid_.numFaces(); ok && (i < nf); ++i) {
if (bctypes[ i ] == FBC_PRESSURE) {
ok = flow_conditions_append(BC_PRESSURE,
static_cast<int>(i),
bcvalues[ i ],
bc_);
}
else if (bctypes[ i ] == FBC_FLUX) {
ok = flow_conditions_append(BC_FLUX_TOTVOL,
static_cast<int>(i),
bcvalues[ i ],
bc_);
}
}
if (! ok) {
flow_conditions_destroy(bc_);
bc_ = 0;
}
}
}; // class TPFACompressiblePressureSolver

23
opm/core/pressure/TPFAPressureSolver.hpp
View File

@ -122,22 +122,6 @@ public:
}
UnstructuredGrid* g = grid_.c_grid();
// Boundary conditions.
int num_faces = g->number_of_faces;
assert(num_faces == int(bctypes.size()));
std::vector<flowbc_type> bctypes2(num_faces, UNSET);
for (int face = 0; face < num_faces; ++face) {
if (bctypes[face] != FBC_FLUX || bcvalues[face] != 0.0) {
throw std::logic_error("TPFAPressureSolver currently only supports noflow bcs.");
}
if (bctypes[face] == FBC_PRESSURE) {
bctypes2[face] = PRESSURE;
} else if (bctypes[face] == FBC_FLUX) {
bctypes2[face] = FLUX;
}
}
// flowbc_t bc = { &bctypes2[0], const_cast<double*>(&bcvalues[0]) };
// Source terms from user.
double* src = const_cast<double*>(&sources[0]); // Ugly? Yes. Safe? I think so.
@ -147,7 +131,7 @@ public:
// double* wdp = 0;
// Compute effective transmissibilities.
eff_trans_.resize(num_faces);
eff_trans_.resize(grid_.numFaces());
tpfa_eff_trans_compute(g, &total_mobilities[0], &htrans_[0], &eff_trans_[0]);
// Update gravity term.
@ -162,8 +146,11 @@ public:
data_->b[i] = 0.0;
}
ifs_tpfa_forces f;
f.src = &src[0];
// Assemble the embedded linear system.
ifs_tpfa_assemble(g, &eff_trans_[0], src, &gpress_[0], data_);
ifs_tpfa_assemble(g, &f, &eff_trans_[0], &gpress_[0], data_);
state_ = Assembled;
}

8
opm/core/pressure/cfsh.c
View File

@ -34,7 +34,7 @@
/* ---------------------------------------------------------------------- */
static int
cfsh_assemble_grid(flowbc_t *bc,
cfsh_assemble_grid(struct FlowBoundaryConditions *bc,
const double *Binv,
const double *gpress,
const double *src,
@ -126,8 +126,8 @@ cfsh_construct(struct UnstructuredGrid *G, well_t *W)
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);
fsh_define_impl_arrays(nc, G->number_of_faces, nnu,
ngconn_tot, new->max_ngconn, W, new->pimpl);
new->pimpl->sys = hybsys_allocate_unsymm(new->max_ngconn,
nc, ngconn_tot);
@ -177,7 +177,7 @@ cfsh_construct(struct UnstructuredGrid *G, well_t *W)
*/
/* ---------------------------------------------------------------------- */
void
cfsh_assemble(flowbc_t *bc,
cfsh_assemble(struct FlowBoundaryConditions *bc,
const double *src,
const double *Binv,
const double *Biv,

174
opm/core/pressure/flow_bc.c
View File

@ -1,5 +1,5 @@
/*
Copyright 2010 SINTEF ICT, Applied Mathematics.
Copyright 2010, 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
@ -17,54 +17,172 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <assert.h>
#include <stdlib.h>
#include "flow_bc.h"
#include <opm/core/pressure/flow_bc.h>
/* Create structure to hold flow boundary conditions for 'nf' faces.
*
* Return fully allocated structure or NULL in case of allocation
* failure. */
/* ---------------------------------------------------------------------- */
flowbc_t *
allocate_flowbc(size_t nf)
/* Compute an appropriate array dimension to minimise total number of
* (re-)allocations. */
/* ---------------------------------------------------------------------- */
static size_t
alloc_size(size_t n, size_t c)
/* ---------------------------------------------------------------------- */
{
size_t i;
flowbc_t *new;
if (c < n) {
c *= 2; /* log_2(n) allocations */
new = malloc(1 * sizeof *new);
if (new != NULL) {
new->type = malloc(nf * sizeof *new->type);
new->bcval = malloc(nf * sizeof *new->bcval);
if ((new->type == NULL) || (new->bcval == NULL)) {
deallocate_flowbc(new);
new = NULL;
} else {
for (i = 0; i < nf; i++) {
new->type [i] = UNSET;
new->bcval[i] = 0.0;
}
if (c < n) {
c = n; /* Typically for the first few allocs */
}
}
return new;
return c;
}
/* Release memory resources for dynamically allocated flow boundary
* condition structure. */
/* ---------------------------------------------------------------------- */
/* Put structure in a well-defined, initial state */
/* ---------------------------------------------------------------------- */
static void
initialise_structure(struct FlowBoundaryConditions *fbc)
/* ---------------------------------------------------------------------- */
{
fbc->nbc = 0;
fbc->cpty = 0;
fbc->type = NULL;
fbc->value = NULL;
fbc->face = NULL;
}
/* ---------------------------------------------------------------------- */
static int
expand_tables(size_t nbc,
struct FlowBoundaryConditions *fbc)
/* ---------------------------------------------------------------------- */
{
int ok;
size_t alloc_sz;
void *p1, *p2, *p3;
ok = nbc <= fbc->cpty;
if (! ok) {
alloc_sz = alloc_size(nbc, fbc->cpty);
p1 = realloc(fbc->type , alloc_sz * sizeof *fbc->type );
p2 = realloc(fbc->value, alloc_sz * sizeof *fbc->value);
p3 = realloc(fbc->face , alloc_sz * sizeof *fbc->face );
ok = (p1 != NULL) && (p2 != NULL) && (p3 != NULL);
if (ok) {
fbc->type = p1;
fbc->value = p2;
fbc->face = p3;
fbc->cpty = alloc_sz;
} else {
free(p3); free(p2); free(p1);
}
}
return ok;
}
/* ======================================================================
* Public interface below separator
* ====================================================================== */
/* ---------------------------------------------------------------------- */
/* Allocate a 'FlowBoundaryConditions' structure, initially capable of
* managing 'nbc' individual boundary conditions. */
/* ---------------------------------------------------------------------- */
struct FlowBoundaryConditions *
flow_conditions_construct(size_t nbc)
/* ---------------------------------------------------------------------- */
{
int ok;
struct FlowBoundaryConditions *fbc;
fbc = malloc(1 * sizeof *fbc);
if (fbc != NULL) {
initialise_structure(fbc);
ok = expand_tables(nbc, fbc);
if (! ok) {
flow_conditions_destroy(fbc);
fbc = NULL;
}
}
return fbc;
}
/* ---------------------------------------------------------------------- */
/* Release memory resources managed by 'fbc', including the containing
* 'struct' pointer, 'fbc'. */
/* ---------------------------------------------------------------------- */
void
deallocate_flowbc(flowbc_t *fbc)
flow_conditions_destroy(struct FlowBoundaryConditions *fbc)
/* ---------------------------------------------------------------------- */
{
if (fbc != NULL) {
free(fbc->bcval);
free(fbc->face );
free(fbc->value);
free(fbc->type );
}
free(fbc);
}
/* ---------------------------------------------------------------------- */
/* Append a new boundary condition to existing set.
*
* Return one (1) if successful, and zero (0) otherwise. */
/* ---------------------------------------------------------------------- */
int
flow_conditions_append(enum FlowBCType type ,
int face ,
double value,
struct FlowBoundaryConditions *fbc )
/* ---------------------------------------------------------------------- */
{
int ok;
ok = expand_tables(fbc->nbc + 1, fbc);
if (ok) {
fbc->type [ fbc->nbc ] = type ;
fbc->value[ fbc->nbc ] = value;
fbc->face [ fbc->nbc ] = face ;
fbc->nbc += 1;
}
return ok;
}
/* ---------------------------------------------------------------------- */
/* Clear existing set of boundary conditions */
/* ---------------------------------------------------------------------- */
void
flow_conditions_clear(struct FlowBoundaryConditions *fbc)
/* ---------------------------------------------------------------------- */
{
assert (fbc != NULL);
fbc->nbc = 0;
}

48
opm/core/pressure/flow_bc.h
View File

@ -22,26 +22,56 @@
#include <stddef.h>
#ifdef __cplusplus
extern "C" {
#endif
enum flowbc_type { UNSET, PRESSURE, FLUX };
enum FlowBCType { BC_NOFLOW ,
BC_PRESSURE ,
BC_FLUX_TOTVOL };
typedef struct {
enum flowbc_type *type;
double *bcval;
} flowbc_t;
/* Boundary condition structure.
*
* Condition i (in [0 .. nbc-1]) affects (outer) interface face[i], is
* of type type[i], and specifies a target value of value[i].
*
* The field 'cpty' is for internal use by the implementation. */
struct FlowBoundaryConditions {
size_t nbc; /* Current number of bdry. conditions */
size_t cpty; /* Internal management. Do not touch */
enum FlowBCType *type; /* Condition type */
double *value; /* Condition value (target) */
int *face; /* Outer faces affected by ind. target */
};
flowbc_t *
allocate_flowbc(size_t nf);
/* Allocate a 'FlowBoundaryConditions' structure, initially capable of
* managing 'nbc' individual boundary conditions. */
struct FlowBoundaryConditions *
flow_conditions_construct(size_t nbc);
/* Release memory resources managed by 'fbc', including the containing
* 'struct' pointer, 'fbc'. */
void
deallocate_flowbc(flowbc_t *fbc);
flow_conditions_destroy(struct FlowBoundaryConditions *fbc);
/* Append a new boundary condition to existing set.
*
* Return one (1) if successful, and zero (0) otherwise. */
int
flow_conditions_append(enum FlowBCType type ,
int face ,
double value,
struct FlowBoundaryConditions *fbc );
/* Clear existing set of boundary conditions */
void
flow_conditions_clear(struct FlowBoundaryConditions *fbc);
#ifdef __cplusplus
}
#endif

4
opm/core/pressure/fsh.h
View File

@ -79,7 +79,7 @@ cfsh_construct(struct UnstructuredGrid *G, well_t *W);
/** Assembles the hybridized linear system for face pressures.
*/
void
cfsh_assemble(flowbc_t *bc,
cfsh_assemble(struct FlowBoundaryConditions *bc,
const double *src,
const double *Binv,
const double *Biv,
@ -132,7 +132,7 @@ ifsh_construct(struct UnstructuredGrid *G, well_t *W);
* be modified). Must already be constructed.
*/
void
ifsh_assemble(flowbc_t *bc,
ifsh_assemble(struct FlowBoundaryConditions *bc,
const double *src,
const double *Binv,
const double *gpress,

49
opm/core/pressure/fsh_common_impl.c
View File

@ -157,23 +157,28 @@ fsh_count_grid_dof(struct UnstructuredGrid *G, int *max_ngdof, size_t *sum_ngdof
/* Impose boundary conditions on local contribution to global system. */
/* ---------------------------------------------------------------------- */
int
fsh_impose_bc(int nconn, int *conn, flowbc_t *bc,
fsh_impose_bc(int nconn, int *conn, struct FlowBoundaryConditions *bc,
struct fsh_impl *pimpl)
/* ---------------------------------------------------------------------- */
{
int i, npp, f;
int i, j, 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;
j = pimpl->bdry_condition[ f ];
npp += 1;
} else if (bc->type[f] == FLUX) {
pimpl->sys->r[i] -= bc->bcval[f];
if (j != -1) {
if (bc->type[j] == BC_PRESSURE) {
pimpl->work [npp] = bc->value[j];
pimpl->iwork[npp] = i;
npp += 1;
} else if (bc->type[j] == BC_FLUX_TOTVOL) {
pimpl->sys->r[i] -= bc->value[j];
}
}
}
@ -189,9 +194,32 @@ fsh_impose_bc(int nconn, int *conn, flowbc_t *bc,
}
/* ---------------------------------------------------------------------- */
void
fsh_map_bdry_condition(struct FlowBoundaryConditions *fbc ,
struct fsh_impl *pimpl)
/* ---------------------------------------------------------------------- */
{
int f, i;
for (i = 0; i < pimpl->nf; i++) { pimpl->bdry_condition[ i ] = -1; }
if (fbc != NULL) {
for (i = 0; ((size_t) i) < fbc->nbc; i++) {
f = fbc->face[ i ];
assert ((0 <= f) && (f < pimpl->nf));
pimpl->bdry_condition[ f ] = i;
}
}
}
/* ---------------------------------------------------------------------- */
void
fsh_define_impl_arrays(size_t nc,
size_t nf,
size_t nnu,
size_t nhf,
size_t max_ncf,
@ -206,8 +234,10 @@ fsh_define_impl_arrays(size_t nc,
pimpl->gdof = pimpl->gdof_pos + (nc + 1);
pimpl->iwork = pimpl->gdof + nhf;
pimpl->bdry_condition = pimpl->iwork + max_ncf;
if (W != NULL) {
pimpl->cwell_pos = pimpl->iwork + max_ncf;
pimpl->cwell_pos = pimpl->bdry_condition + nf;
pimpl->cwells = pimpl->cwell_pos + nc + 1;
pimpl->WI = pimpl->work + max_ncf;
@ -256,6 +286,7 @@ fsh_compute_table_sz(struct UnstructuredGrid *G, well_t *W, int max_ngconn,
*idata_sz = nc + 1; /* gdof_pos */
*idata_sz += ngconn_tot; /* gdof */
*idata_sz += G->number_of_faces; /* bdry_condition */
*idata_sz += max_ngconn; /* iwork */
*ddata_sz = 2 * (*nnu); /* rhs + soln */

9
opm/core/pressure/fsh_common_impl.h
View File

@ -44,6 +44,8 @@ struct fsh_impl {
double *work; /* Scratch array, floating point */
int *iwork; /* Scratch array, integers */
int *bdry_condition; /* Map face->boundary condition ID */
/* Linear storage goes here... */
int *idata; /* Actual storage array, integers */
double *ddata; /* Actual storage array, floating point */
@ -56,14 +58,19 @@ fsh_impl_allocate_basic(size_t idata_sz, size_t ddata_sz);
void
fsh_count_grid_dof(struct UnstructuredGrid *G, int *max_ngdof, size_t *sum_ngdof2);
void
fsh_map_bdry_condition(struct FlowBoundaryConditions *fbc ,
struct fsh_impl *pimpl);
int
fsh_impose_bc(int ndof,
int *dof,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
struct fsh_impl *pimpl);
void
fsh_define_impl_arrays(size_t nc,
size_t nf,
size_t nnu,
size_t nhf,
size_t max_ncf,

28
opm/core/pressure/ifsh.c
View File

@ -60,7 +60,7 @@ ifsh_set_effective_well_params(const double *WI,
/* ---------------------------------------------------------------------- */
static int
ifsh_assemble_grid(flowbc_t *bc,
ifsh_assemble_grid(struct FlowBoundaryConditions *bc,
const double *Binv,
const double *gpress,
const double *src,
@ -100,12 +100,12 @@ ifsh_assemble_grid(flowbc_t *bc,
/* ---------------------------------------------------------------------- */
static void
ifsh_impose_well_control(int c,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
well_control_t *wctrl,
struct fsh_data *ifsh)
/* ---------------------------------------------------------------------- */
{
int ngconn, nwconn, i, w1, w2, wg, f;
int ngconn, nwconn, i, j, w1, w2, wg, f;
int *pgconn, *gconn, *pwconn, *wconn;
double bhp;
@ -130,11 +130,12 @@ ifsh_impose_well_control(int c,
/* Adapt local system to prescribed boundary pressures (r->w) */
for (i = 0; i < ngconn; i++) {
f = gconn[i];
j = ifsh->pimpl->bdry_condition[ f ];
if (bc->type[f] == PRESSURE) {
if (j != -1) {
for (w1 = 0; w1 < nwconn; w1++) {
/* Eliminate prescribed (boundary) pressure value */
r [ngconn + w1] -= r2w[i + w1*ngconn] * bc->bcval[f];
r [ngconn + w1] -= r2w[i + w1*ngconn] * bc->value[j];
r2w[i + w1*ngconn] = 0.0;
}
@ -152,8 +153,13 @@ ifsh_impose_well_control(int c,
/* Well->reservoir */
for (i = 0; i < ngconn; i++) {
assert ((bc->type[gconn[i]] != PRESSURE) ||
#ifndef NDEBUG
j = ifsh->pimpl->bdry_condition[ gconn[i] ];
assert ((j == -1) ||
(bc->type[j] != BC_PRESSURE) ||
!(fabs(r2w[i + w1*ngconn]) > 0.0));
#endif
r [i] -= r2w[i + w1*ngconn] * bhp;
r2w[i + w1*ngconn] = 0.0;
@ -177,7 +183,7 @@ ifsh_impose_well_control(int c,
/* ---------------------------------------------------------------------- */
static int
ifsh_assemble_well(flowbc_t *bc,
ifsh_assemble_well(struct FlowBoundaryConditions *bc,
well_control_t *wctrl,
struct fsh_data *ifsh)
/* ---------------------------------------------------------------------- */
@ -291,8 +297,8 @@ ifsh_construct(struct UnstructuredGrid *G, well_t *W)
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);
fsh_define_impl_arrays(nc, G->number_of_faces, nnu,
ngconn_tot, new->max_ngconn, W, new->pimpl);
new->pimpl->sys = hybsys_allocate_symm(new->max_ngconn,
nc, ngconn_tot);
@ -343,7 +349,7 @@ ifsh_construct(struct UnstructuredGrid *G, well_t *W)
* pressure gpress, boundary conditions bc, and source terms src. */
/* ---------------------------------------------------------------------- */
void
ifsh_assemble(flowbc_t *bc,
ifsh_assemble(struct FlowBoundaryConditions *bc,
const double *src,
const double *Binv,
const double *gpress,
@ -355,6 +361,8 @@ ifsh_assemble(flowbc_t *bc,
{
int npp; /* Number of prescribed pressure values */
fsh_map_bdry_condition(bc, ifsh->pimpl);
hybsys_schur_comp_symm(ifsh->pimpl->nc,
ifsh->pimpl->gdof_pos,
Binv, ifsh->pimpl->sys);

199
opm/core/pressure/tpfa/cfs_tpfa.c
View File

@ -36,6 +36,9 @@ struct densrat_util {
struct cfs_tpfa_impl {
/* Mapping (cell,face) -> half-face */
int *f2hf;
/* Reservoir flow */
double *ctrans;
double *accum;
@ -140,7 +143,7 @@ static struct cfs_tpfa_impl *
impl_allocate(struct UnstructuredGrid *G, well_t *W, int np)
/* ---------------------------------------------------------------------- */
{
size_t nnu, ngconn, nwperf;
size_t nnu, nf, ngconn, nwperf;
struct cfs_tpfa_impl *new;
size_t ddata_sz;
@ -158,10 +161,11 @@ impl_allocate(struct UnstructuredGrid *G, well_t *W, int np)
ddata_sz = 2 * nnu; /* b, x */
/* Reservoir */
nf = G->number_of_faces;
ddata_sz += 1 * ngconn; /* ctrans */
ddata_sz += 1 * G->number_of_cells; /* accum */
ddata_sz += np * G->number_of_faces; /* masstrans_f */
ddata_sz += np * G->number_of_faces; /* gravtrans_f */
ddata_sz += np * nf; /* masstrans_f */
ddata_sz += np * nf; /* gravtrans_f */
/* Wells */
ddata_sz += 1 * nwperf; /* wtrans */
@ -169,15 +173,17 @@ impl_allocate(struct UnstructuredGrid *G, well_t *W, int np)
ddata_sz += np * nwperf; /* masstrans_p */
ddata_sz += np * nwperf; /* gravtrans_p */
ddata_sz += 1 * G->number_of_faces; /* scratch_f */
ddata_sz += 1 * nf; /* scratch_f */
new = malloc(1 * sizeof *new);
if (new != NULL) {
new->f2hf = malloc(2 * nf * sizeof *new->f2hf );
new->ddata = malloc(ddata_sz * sizeof *new->ddata);
new->ratio = allocate_densrat(G, W, np);
if (new->ddata == NULL || new->ratio == NULL) {
if (new->f2hf == NULL ||
new->ddata == NULL || new->ratio == NULL) {
impl_deallocate(new);
new = NULL;
}
@ -392,7 +398,7 @@ set_dynamic_trans(struct UnstructuredGrid *G ,
/* ---------------------------------------------------------------------- */
static void
set_dynamic_grav(struct UnstructuredGrid *G ,
flowbc_t *bc ,
struct FlowBoundaryConditions *bc ,
const double *trans ,
const double *gravcap_f,
struct compr_quantities *cq ,
@ -405,7 +411,7 @@ set_dynamic_grav(struct UnstructuredGrid *G ,
c1 = G->face_cells[2*f + 0];
c2 = G->face_cells[2*f + 1];
if (((c1 >= 0) && (c2 >= 0)) || (bc->type[f] == PRESSURE)) {
if ((c1 >= 0) && (c2 >= 0)) {
for (p = 0; p < cq->nphases; p++, i++) {
ratio->x[i] = trans[f] * gravcap_f[i] * cq->phasemobf[i];
}
@ -415,6 +421,18 @@ set_dynamic_grav(struct UnstructuredGrid *G ,
}
}
}
if (bc != NULL) {
for (i = 0; ((size_t) i) < bc->nbc; i++) {
if (bc->type[ i ] == BC_PRESSURE) {
f = bc->face[ i ];
for (p = cq->nphases * f; p < cq->nphases * (f + 1); p++) {
ratio->x[p] = trans[f] * gravcap_f[p] * cq->phasemobf[p];
}
}
}
}
}
@ -561,7 +579,7 @@ static void
compute_psys_contrib(struct UnstructuredGrid *G,
well_t *W,
struct completion_data *wdata,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
struct compr_quantities *cq,
double dt,
const double *trans,
@ -626,20 +644,16 @@ compute_psys_contrib(struct UnstructuredGrid *G,
/* ---------------------------------------------------------------------- */
static int
assemble_cell_contrib(struct UnstructuredGrid *G,
flowbc_t *bc,
const double *src,
struct cfs_tpfa_data *h)
static void
assemble_cell_contrib(struct UnstructuredGrid *G ,
const double *src,
struct cfs_tpfa_data *h )
/* ---------------------------------------------------------------------- */
{
int c1, c2, c, i, f, j1, j2;
int is_neumann;
const double *ctrans = h->pimpl->ctrans;
is_neumann = 1;
for (c = i = 0; c < G->number_of_cells; c++) {
j1 = csrmatrix_elm_index(c, c, h->A);
@ -656,11 +670,6 @@ assemble_cell_contrib(struct UnstructuredGrid *G,
h->A->sa[j1] += ctrans[i];
h->A->sa[j2] -= ctrans[i];
} else if (bc->type[f] == PRESSURE) {
is_neumann = 0;
h->A->sa[j1] += ctrans[i];
h->b [c ] += ctrans[i] * bc->bcval[f];
}
}
@ -669,6 +678,47 @@ assemble_cell_contrib(struct UnstructuredGrid *G,
/* Compressible accumulation term */
h->A->sa[j1] += h->pimpl->accum[c];
}
}
/* ---------------------------------------------------------------------- */
static int
assemble_bc_contrib(struct UnstructuredGrid *G ,
struct FlowBoundaryConditions *fbc,
struct cfs_tpfa_data *h )
/* ---------------------------------------------------------------------- */
{
int c1, c2, c, f, hf;
int is_neumann;
size_t p, j;
const double *ctrans = h->pimpl->ctrans;
is_neumann = 1;
for (p = 0; p < fbc->nbc; p++) {
f = fbc->face[ p ];
c1 = G->face_cells[2*f + 0];
c2 = G->face_cells[2*f + 1];
assert ((c1 < 0) ^ (c2 < 0));
c = (c1 >= 0) ? c1 : c2;
if (fbc->type[ p ] == BC_PRESSURE) {
is_neumann = 0;
hf = h->pimpl->f2hf[2*f + (c1 < 0)];
j = csrmatrix_elm_index(c, c, h->A);
h->A->sa[ j ] += ctrans[ hf ];
h->b [ c ] += ctrans[ hf ] * fbc->value[ p ];
}
/* Other boundary condition types not handled. */
}
return is_neumann;
}
@ -735,7 +785,7 @@ assemble_well_contrib(size_t nc,
/* ---------------------------------------------------------------------- */
static void
compute_fpress(struct UnstructuredGrid *G,
flowbc_t *bc,
struct FlowBoundaryConditions *fbc,
int np,
const double *htrans,
const double *pmobf,
@ -746,7 +796,7 @@ compute_fpress(struct UnstructuredGrid *G,
double *scratch_f)
/* ---------------------------------------------------------------------- */
{
int c, i, f, c1, c2;
int c, i, f;
/* Suppress warning about unused parameters. */
(void) np; (void) pmobf; (void) gravcap_f; (void) fflux;
@ -777,12 +827,13 @@ compute_fpress(struct UnstructuredGrid *G,
for (f = 0; f < G->number_of_faces; f++) {
fpress[f] /= scratch_f[f];
}
c1 = G->face_cells[2*f + 0];
c2 = G->face_cells[2*f + 1];
if (((c1 < 0) || (c2 < 0)) && (bc->type[f] == PRESSURE)) {
fpress[f] = bc->bcval[f];
if (fbc != NULL) {
for (i = 0; ((size_t) i) < fbc->nbc; i++) {
if (fbc->type[ i ] == BC_PRESSURE) {
fpress[ fbc->face[ i ] ] = fbc->value[ i ];
}
}
}
}
@ -791,7 +842,7 @@ compute_fpress(struct UnstructuredGrid *G,
/* ---------------------------------------------------------------------- */
static void
compute_flux(struct UnstructuredGrid *G,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
int np,
const double *trans,
const double *pmobf,
@ -800,15 +851,14 @@ compute_flux(struct UnstructuredGrid *G,
double *fflux)
/* ---------------------------------------------------------------------- */
{
int f, c1, c2, p;
int f, c1, c2, p, i;
double t, dp, g;
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)) && (bc->type[f] == FLUX)) {
fflux[f] = bc->bcval[f];
if ((c1 < 0) || (c2 < 0)) {
continue;
}
@ -817,24 +867,42 @@ compute_flux(struct UnstructuredGrid *G,
t += pmobf[f*np + p];
g += pmobf[f*np + p] * gravcap_f[f*np + p];
}
/* t *= trans[f]; */
if ((c1 >= 0) && (c2 >= 0)) {
dp = cpress[c1] - cpress[c2];
} else if (bc->type[f] == PRESSURE) {
if (c1 < 0) {
dp = bc->bcval[f] - cpress[c2];
/* g = -g; */
} else {
dp = cpress[c1] - bc->bcval[f];
}
} else {
/* No BC -> no-flow (== pressure drop offsets gravity) */
dp = -g / t;
}
dp = cpress[c1] - cpress[c2];
fflux[f] = trans[f] * (t*dp + g);
}
if (bc != NULL) {
for (i = 0; ((size_t) i) < bc->nbc; i++) {
f = bc->face[ i ];
if (bc->type[ i ] == BC_FLUX_TOTVOL) {
t = 2*(G->face_cells[2*f + 0] < 0) - 1.0;
fflux[ f ] = t * bc->value[ i ];
}
else if (bc->type[ i ] == BC_PRESSURE) {
c1 = G->face_cells[2*f + 0];
c2 = G->face_cells[2*f + 1];
t = g = 0.0;
for (p = 0; p < np; p++) {
t += pmobf[f*np + p];
g += pmobf[f*np + p] * gravcap_f[f*np + p];
}
if (c1 < 0) {
dp = bc->value[ i ] - cpress[c1];
}
else {
dp = cpress[c2] - bc->value[ i ];
}
fflux[f] = trans[f] * (t*dp + g);
}
}
}
}
@ -889,6 +957,29 @@ is_incompr(int nc, struct compr_quantities *cq)
}
/* ---------------------------------------------------------------------- */
static void
compute_f2hf_mapping(struct UnstructuredGrid *G ,
struct cfs_tpfa_data *data)
/* ---------------------------------------------------------------------- */
{
int c, nc, i, f, p;
for (i = 0; i < 2 * G->number_of_faces; i++) {
data->pimpl->f2hf[ i ] = -1;
}
for (c = i = 0, nc = G->number_of_cells; c < nc; c++) {
for (; i < G->cell_facepos[c + 1]; i++) {
f = G->cell_faces[ i ];
p = G->cell_faces[2*f + 0] != c;
data->pimpl->f2hf[2*f + p] = i;
}
}
}
/* ======================================================================
* Public interface below separator.
@ -924,6 +1015,8 @@ cfs_tpfa_construct(struct UnstructuredGrid *G, well_t *W, int nphases)
nwconn = W->well_connpos[ W->number_of_wells ];
}
compute_f2hf_mapping(G, new);
/* Allocate linear system components */
new->b = new->pimpl->ddata + 0;
new->x = new->b + new->A->m;
@ -953,7 +1046,7 @@ void
cfs_tpfa_assemble(struct UnstructuredGrid *G,
double dt,
well_t *W,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
const double *src,
struct compr_quantities *cq,
const double *trans,
@ -973,7 +1066,13 @@ cfs_tpfa_assemble(struct UnstructuredGrid *G,
compute_psys_contrib(G, W, wdata, bc, cq, dt,
trans, gravcap_f, cpress0, porevol, h);
res_is_neumann = assemble_cell_contrib(G, bc, src, h);
res_is_neumann = 1;
assemble_cell_contrib(G, src, h);
if (bc != NULL) {
res_is_neumann = assemble_bc_contrib(G, bc, h);
}
if ((W != NULL) && (wctrl != NULL)) {
assert (wdata != NULL);
@ -993,7 +1092,7 @@ cfs_tpfa_assemble(struct UnstructuredGrid *G,
/* ---------------------------------------------------------------------- */
void
cfs_tpfa_press_flux(struct UnstructuredGrid *G,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
well_t *W,
int np,
const double *trans,
@ -1028,7 +1127,7 @@ cfs_tpfa_press_flux(struct UnstructuredGrid *G,
/* ---------------------------------------------------------------------- */
void
cfs_tpfa_fpress(struct UnstructuredGrid *G,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
int np,
const double *htrans,
const double *pmobf,

6
opm/core/pressure/tpfa/cfs_tpfa.h
View File

@ -48,7 +48,7 @@ void
cfs_tpfa_assemble(struct UnstructuredGrid *G,
double dt,
well_t *W,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
const double *src,
struct compr_quantities *cq,
const double *trans,
@ -61,7 +61,7 @@ cfs_tpfa_assemble(struct UnstructuredGrid *G,
void
cfs_tpfa_press_flux(struct UnstructuredGrid *G,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
well_t *W,
int np,
const double *trans,
@ -76,7 +76,7 @@ cfs_tpfa_press_flux(struct UnstructuredGrid *G,
void
cfs_tpfa_fpress(struct UnstructuredGrid *G,
flowbc_t *bc,
struct FlowBoundaryConditions *bc,
int np,
const double *htrans,
const double *pmobf,

6
opm/core/pressure/tpfa/ifs_tpfa.c
View File

@ -186,8 +186,8 @@ ifs_tpfa_construct(struct UnstructuredGrid *G)
/* ---------------------------------------------------------------------- */
void
ifs_tpfa_assemble(struct UnstructuredGrid *G,
const struct ifs_tpfa_forces *F,
const double *trans,
const double *src,
const double *gpress,
struct ifs_tpfa_data *h)
/* ---------------------------------------------------------------------- */
@ -223,7 +223,9 @@ ifs_tpfa_assemble(struct UnstructuredGrid *G,
}
}
h->b[c] += src[c];
if ((F != NULL) && (F->src != NULL)) {
h->b[c] += F->src[c];
}
}
h->A->sa[0] *= 2;

7
opm/core/pressure/tpfa/ifs_tpfa.h
View File

@ -38,13 +38,18 @@ struct ifs_tpfa_data {
};
struct ifs_tpfa_forces {
const double *src;
};
struct ifs_tpfa_data *
ifs_tpfa_construct(struct UnstructuredGrid *G);
void
ifs_tpfa_assemble(struct UnstructuredGrid *G,
const struct ifs_tpfa_forces *F,
const double *trans,
const double *src,
const double *gpress,
struct ifs_tpfa_data *h);

5
tests/test_cfs_tpfa.c
View File

@ -615,7 +615,7 @@ int
main(void)
{
int i, nphases, dim;
flowbc_t *bc;
struct FlowBoundaryConditions *bc;
double dt;
double *htrans, *trans, *perm, *totmob, *src;
double *gravcap_f, *porevol, *cpress0;
@ -655,7 +655,7 @@ main(void)
wpress = malloc(W->number_of_wells * sizeof *wpress);
wflux = malloc(W->well_connpos[ W->number_of_wells ] * sizeof *wflux);
bc = allocate_flowbc(G->number_of_faces);
bc = NULL;
set_homoperm(G->number_of_cells, G->dimensions, perm);
@ -690,7 +690,6 @@ main(void)
cfs_tpfa_destroy(h);
deallocate_cq(cq);
deallocate_flowbc(bc);
free(wflux); free(wpress);
free(fflux); free(fpress); free(cpress);