Merge pull request #1184 from akva2/remove_unused_sources
remove unused code
This commit is contained in:
Atgeirr Flø Rasmussen
GitHub
commit
c3121f0a1a
@ -45,28 +45,13 @@ list (APPEND MAIN_SOURCE_FILES
|
||||
opm/core/pressure/FlowBCManager.cpp
|
||||
opm/core/pressure/IncompTpfa.cpp
|
||||
opm/core/pressure/IncompTpfaSinglePhase.cpp
|
||||
opm/core/pressure/cfsh.c
|
||||
opm/core/pressure/flow_bc.c
|
||||
opm/core/pressure/fsh.c
|
||||
opm/core/pressure/fsh_common_impl.c
|
||||
opm/core/pressure/ifsh.c
|
||||
opm/core/pressure/legacy_well.c
|
||||
opm/core/pressure/mimetic/hybsys.c
|
||||
opm/core/pressure/mimetic/hybsys_global.c
|
||||
opm/core/pressure/mimetic/mimetic.c
|
||||
opm/core/pressure/msmfem/coarse_conn.c
|
||||
opm/core/pressure/msmfem/coarse_sys.c
|
||||
opm/core/pressure/msmfem/dfs.c
|
||||
opm/core/pressure/msmfem/hash_set.c
|
||||
opm/core/pressure/msmfem/ifsh_ms.c
|
||||
opm/core/pressure/msmfem/partition.c
|
||||
opm/core/pressure/tpfa/TransTpfa.cpp
|
||||
opm/core/pressure/tpfa/cfs_tpfa.c
|
||||
opm/core/pressure/tpfa/cfs_tpfa_residual.c
|
||||
opm/core/pressure/tpfa/compr_bc.c
|
||||
opm/core/pressure/tpfa/compr_quant.c
|
||||
opm/core/pressure/tpfa/compr_quant_general.c
|
||||
opm/core/pressure/tpfa/compr_source.c
|
||||
opm/core/pressure/tpfa/ifs_tpfa.c
|
||||
opm/core/pressure/tpfa/trans_tpfa.c
|
||||
opm/core/props/BlackoilPropertiesBasic.cpp
|
||||
@ -88,8 +73,6 @@ list (APPEND MAIN_SOURCE_FILES
|
||||
opm/core/transport/TransportSolverTwophaseInterface.cpp
|
||||
opm/core/transport/implicit/TransportSolverTwophaseImplicit.cpp
|
||||
opm/core/transport/implicit/transport_source.c
|
||||
opm/core/transport/minimal/spu_explicit.c
|
||||
opm/core/transport/minimal/spu_implicit.c
|
||||
opm/core/transport/reorder/ReorderSolverInterface.cpp
|
||||
opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp
|
||||
opm/core/transport/reorder/TransportSolverTwophaseReorder.cpp
|
||||
@ -233,25 +216,15 @@ list (APPEND PUBLIC_HEADER_FILES
|
||||
opm/core/pressure/CompressibleTpfa.hpp
|
||||
opm/core/pressure/FlowBCManager.hpp
|
||||
opm/core/pressure/IncompTpfa.hpp
|
||||
opm/core/pressure/IncompTpfaSinglePhase.hpp
|
||||
opm/core/pressure/flow_bc.h
|
||||
opm/core/pressure/fsh.h
|
||||
opm/core/pressure/fsh_common_impl.h
|
||||
opm/core/pressure/legacy_well.h
|
||||
opm/core/pressure/mimetic/hybsys.h
|
||||
opm/core/pressure/mimetic/hybsys_global.h
|
||||
opm/core/pressure/mimetic/mimetic.h
|
||||
opm/core/pressure/msmfem/coarse_conn.h
|
||||
opm/core/pressure/msmfem/coarse_sys.h
|
||||
opm/core/pressure/msmfem/dfs.h
|
||||
opm/core/pressure/msmfem/hash_set.h
|
||||
opm/core/pressure/msmfem/ifsh_ms.h
|
||||
opm/core/pressure/msmfem/partition.h
|
||||
opm/core/pressure/tpfa/TransTpfa.hpp
|
||||
opm/core/pressure/tpfa/TransTpfa_impl.hpp
|
||||
opm/core/pressure/tpfa/cfs_tpfa.h
|
||||
opm/core/pressure/tpfa/cfs_tpfa_residual.h
|
||||
opm/core/pressure/tpfa/compr_bc.h
|
||||
opm/core/pressure/tpfa/compr_quant.h
|
||||
opm/core/pressure/tpfa/compr_quant_general.h
|
||||
opm/core/pressure/tpfa/compr_source.h
|
||||
@ -305,8 +278,6 @@ list (APPEND PUBLIC_HEADER_FILES
|
||||
opm/core/transport/implicit/SinglePointUpwindTwoPhase.hpp
|
||||
opm/core/transport/implicit/TransportSolverTwophaseImplicit.hpp
|
||||
opm/core/transport/implicit/transport_source.h
|
||||
opm/core/transport/minimal/spu_explicit.h
|
||||
opm/core/transport/minimal/spu_implicit.h
|
||||
opm/core/transport/reorder/ReorderSolverInterface.hpp
|
||||
opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp
|
||||
opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp
|
||||
|
||||
@ -1,210 +0,0 @@
|
||||
/*
|
||||
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 "config.h"
|
||||
#include <assert.h>
|
||||
#include <limits.h>
|
||||
#include <math.h>
|
||||
#include <stddef.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include <opm/core/pressure/fsh.h>
|
||||
#include <opm/core/pressure/fsh_common_impl.h>
|
||||
#include <opm/core/pressure/mimetic/hybsys.h>
|
||||
#include <opm/core/pressure/mimetic/hybsys_global.h>
|
||||
|
||||
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static int
|
||||
cfsh_assemble_grid(struct FlowBoundaryConditions *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 *
|
||||
cfsh_construct(struct UnstructuredGrid *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, 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);
|
||||
|
||||
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
|
||||
cfsh_assemble(struct FlowBoundaryConditions *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 */
|
||||
|
||||
/* Suppress warnings about unused parameters. */
|
||||
(void) wctrl; (void) WI; (void) BivW; (void) wdp;
|
||||
|
||||
hybsys_schur_comp_unsymm(h->pimpl->nc,
|
||||
h->pimpl->gdof_pos,
|
||||
Binv, Biv, P, h->pimpl->sys);
|
||||
|
||||
fsh_map_bdry_condition(bc, h->pimpl);
|
||||
|
||||
npp = cfsh_assemble_grid(bc, Binv, gpress, src, h);
|
||||
|
||||
if (npp == 0) {
|
||||
h->A->sa[0] *= 2; /* Remove zero eigenvalue */
|
||||
}
|
||||
}
|
||||
@ -1,90 +0,0 @@
|
||||
/*
|
||||
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 "config.h"
|
||||
#include <assert.h>
|
||||
#include <limits.h>
|
||||
#include <math.h>
|
||||
#include <stddef.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include <opm/core/pressure/fsh.h>
|
||||
#include <opm/core/pressure/fsh_common_impl.h>
|
||||
#include <opm/core/pressure/mimetic/hybsys.h>
|
||||
#include <opm/core/pressure/mimetic/hybsys_global.h>
|
||||
|
||||
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
/* 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(struct UnstructuredGrid *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;
|
||||
}
|
||||
}
|
||||
@ -1,243 +0,0 @@
|
||||
/*
|
||||
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_FSH_HEADER_INCLUDED
|
||||
|
||||
/**
|
||||
* \file
|
||||
* Routines and data structures to support the construction and
|
||||
* formation of hybridized pressure solvers based on Schur
|
||||
* complement reductions.
|
||||
*
|
||||
* Pressure solvers based on this strategy will be structured
|
||||
* according to the following scheme
|
||||
* -# Construct @c FSH data object suitable for the particular
|
||||
* problem using either of the functions cfsh_construct() or
|
||||
* ifsh_construct() for compressible or incompressible flow,
|
||||
* respectively
|
||||
* -# Compute static discretisation quantities, for instance
|
||||
* using functions mim_ip_simple_all() and mim_ip_compute_gpress()
|
||||
* -# For each time step or non-linear iteration:
|
||||
* -# Compute dynamic discretisation quantities incorporating
|
||||
* effects of multiple phases and non-linear feedback.
|
||||
* -# Assemble system of simultaneous linear equations using
|
||||
* functions cfsh_assemble() or ifsh_assemble()
|
||||
* -# Solve the resulting system of linear equations, available
|
||||
* in the @c A and @c b objects of the @c FSH data object,
|
||||
* using some linear solver software. The solution should
|
||||
* be placed in the @c x object of the @c FSH object.
|
||||
* -# Reconstruct derived quantities such as cell pressures and
|
||||
* interface fluxes using function fsh_press_flux().
|
||||
* Function fsh_press_flux() relies on the solution to the
|
||||
* linear system being stored in the @c x object.
|
||||
* -# Release resources using function fsh_destroy() at end of
|
||||
* simulation.
|
||||
*/
|
||||
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/pressure/legacy_well.h>
|
||||
#include <opm/core/pressure/flow_bc.h>
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
|
||||
struct CSRMatrix;
|
||||
struct fsh_impl;
|
||||
|
||||
/**
|
||||
* Main data structure of hybridized pressure solvers based on Schur
|
||||
* complement reductions. Mainly intended to present a common view
|
||||
* of a Schur complement system of simultaneous linear equations
|
||||
* and to hold the solution of said system.
|
||||
*/
|
||||
struct fsh_data {
|
||||
/**
|
||||
* Maximum number of connections in any grid cell,
|
||||
* \f[
|
||||
* \mathit{max\_ngconn} = \max_c \{ n_c \}
|
||||
* \f]
|
||||
* in which \f$n_c\f$ denotes the number connections
|
||||
* (i.e., faces) of cell \f$c\f$.
|
||||
*/
|
||||
int max_ngconn;
|
||||
|
||||
/**
|
||||
* Sum of squared number of connections in all grid cells,
|
||||
* \f[
|
||||
* \mathit{sum\_ngconn2} = \sum_c n_c^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;
|
||||
};
|
||||
|
||||
|
||||
|
||||
/**
|
||||
* Dispose of all memory associated to <CODE>FSH</CODE> object.
|
||||
*
|
||||
* @param[in,out] h <CODE>FSH</CODE> object. Following a call
|
||||
* to function fsh_destroy(), the pointer is
|
||||
* invalid.
|
||||
*/
|
||||
void
|
||||
fsh_destroy(struct fsh_data *h);
|
||||
|
||||
|
||||
/**
|
||||
* Construct compressible hybrid flow-solver data object for a
|
||||
* given grid and well pattern.
|
||||
*
|
||||
* @param[in] G Grid.
|
||||
* @param[in] W Well topology. Ignored.
|
||||
* @return Fully formed data object suitable for use in a
|
||||
* compressible pressure solver. @c NULL in case of construction
|
||||
* failure.
|
||||
*/
|
||||
struct fsh_data *
|
||||
cfsh_construct(struct UnstructuredGrid *G, well_t *W);
|
||||
|
||||
|
||||
|
||||
/**
|
||||
* Form Schur-complement system of simultaneous linear equations
|
||||
* arising in compressible flow using a hybridized formulation.
|
||||
*
|
||||
* Upon returning from function cfsh_assemble(), the resulting
|
||||
* system of simultaneous linear equations is stored in
|
||||
* <CODE>h->A</CODE> and <CODE>h->b</CODE>.
|
||||
*
|
||||
* @param[in] bc Boundary conditions.
|
||||
* @param[in] src Explicit source terms.
|
||||
* @param[in] Binv Inverse of block-diagonal matrix \f$B\f$
|
||||
* Typically computed using functions
|
||||
* mim_ip_simple_all() and
|
||||
* mim_ip_mobility_update().
|
||||
* @param[in] Biv \f$B^{-1}v\f$.
|
||||
* @param[in] P Compressible accumulation term.
|
||||
* @param[in] gpress Gravity pressure.
|
||||
* @param[in] wctrl Well controls. Ignored.
|
||||
* @param[in] WI Well indices. Ignored.
|
||||
* @param[in] BivW \f$B^{-1}v\f$ for wells. Ignored.
|
||||
* @param[in] wdp Gravity pressure along well track. Ignored.
|
||||
* @param[in,out] h Data object.
|
||||
*/
|
||||
void
|
||||
cfsh_assemble(struct FlowBoundaryConditions *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);
|
||||
|
||||
|
||||
/**
|
||||
* Construct incompressible hybrid flow-solver data object for a
|
||||
* given grid and well pattern.
|
||||
*
|
||||
* @param G Grid.
|
||||
* @param W Well topology.
|
||||
* @return Fully formed data object suitable for use in an
|
||||
* incompressible pressure solver. @c NULL in case of construction
|
||||
* failure.
|
||||
*/
|
||||
struct fsh_data *
|
||||
ifsh_construct(struct UnstructuredGrid *G, well_t *W);
|
||||
|
||||
|
||||
/**
|
||||
* Form Schur-complement system of simultaneous linear equations
|
||||
* arising in compressible flow using a hybridized formulation.
|
||||
*
|
||||
* Upon returning from function cfsh_assemble(), the resulting
|
||||
* system of simultaneous linear equations is stored in
|
||||
* <CODE>h->A</CODE> and <CODE>h->b</CODE>.
|
||||
*
|
||||
* @param[in] bc Boundary conditions.
|
||||
* @param[in] src Explicit source terms.
|
||||
* @param[in] Binv Inverse of block-diagonal matrix \f$B\f$
|
||||
* Typically computed using functions
|
||||
* mim_ip_simple_all() and
|
||||
* mim_ip_mobility_update().
|
||||
* @param[in] gpress Gravity pressure.
|
||||
* @param[in] wctrl Well controls.
|
||||
* @param[in] WI Well indices.
|
||||
* @param[in] wdp Gravity pressure along well track.
|
||||
* @param[in,out] h Data object.
|
||||
*/
|
||||
void
|
||||
ifsh_assemble(struct FlowBoundaryConditions *bc,
|
||||
const double *src,
|
||||
const double *Binv,
|
||||
const double *gpress,
|
||||
well_control_t *wctrl,
|
||||
const double *WI,
|
||||
const double *wdp,
|
||||
struct fsh_data *h);
|
||||
|
||||
|
||||
|
||||
|
||||
/**
|
||||
* Compute cell pressures (cpress) and interface fluxes (fflux) from
|
||||
* current solution of system of linear equations, <CODE>h->x</CODE>.
|
||||
* Back substitution process, projected half-contact fluxes.
|
||||
*
|
||||
* @param[in] G Grid.
|
||||
* @param[in] Binv Inverse of block-diagonal matrix \f$B\f$
|
||||
* Must coincide with the matrix used to
|
||||
* form the system of linear equations
|
||||
* currently stored in the data object.
|
||||
* @param[in] gpress Gravity pressure. Must coincide with
|
||||
* the array used to form the system of
|
||||
* linear equations.
|
||||
* @param[in] h Data object.
|
||||
* @param[out] cpress Cell pressure.
|
||||
* @param[out] fflux Interface fluxes.
|
||||
* @param[out] wpress Well pressure.
|
||||
* @param[out] wflux Well perforation fluxes.
|
||||
*/
|
||||
void
|
||||
fsh_press_flux(struct UnstructuredGrid *G,
|
||||
const double *Binv, const double *gpress,
|
||||
struct fsh_data *h,
|
||||
double *cpress, double *fflux,
|
||||
double *wpress, double *wflux);
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
#endif /* OPM_FSH_HEADER_INCLUDED */
|
||||
@ -1,317 +0,0 @@
|
||||
/*
|
||||
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 "config.h"
|
||||
#include <assert.h>
|
||||
#include <limits.h>
|
||||
#include <math.h>
|
||||
#include <stddef.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/pressure/legacy_well.h>
|
||||
#include <opm/core/pressure/flow_bc.h>
|
||||
|
||||
#include <opm/core/pressure/fsh.h>
|
||||
#include <opm/core/pressure/fsh_common_impl.h>
|
||||
|
||||
#include <opm/core/pressure/mimetic/hybsys.h>
|
||||
#include <opm/core/pressure/mimetic/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(struct UnstructuredGrid *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, struct FlowBoundaryConditions *bc,
|
||||
struct fsh_impl *pimpl)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int i, j, npp, f;
|
||||
|
||||
npp = 0;
|
||||
for (i = 0; i < nconn; i++) {
|
||||
f = conn[i];
|
||||
|
||||
j = pimpl->bdry_condition[ 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];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (npp > 0) {
|
||||
fsh_explicit_elimination(nconn, npp,
|
||||
pimpl->iwork,
|
||||
pimpl->work,
|
||||
pimpl->sys->S,
|
||||
pimpl->sys->r);
|
||||
}
|
||||
|
||||
return npp;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
fsh_map_bdry_condition(struct FlowBoundaryConditions *fbc ,
|
||||
struct fsh_impl *pimpl)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int f, i;
|
||||
size_t j;
|
||||
|
||||
for (i = 0; i < pimpl->nf; i++) { pimpl->bdry_condition[ i ] = -1; }
|
||||
|
||||
if (fbc != NULL) {
|
||||
assert (fbc->cond_pos[0] == 0);
|
||||
|
||||
for (i = 0, j = 0; ((size_t) i) < fbc->nbc; i++) {
|
||||
for (; j < fbc->cond_pos[i + 1]; j++) {
|
||||
f = fbc->face[ j ];
|
||||
|
||||
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,
|
||||
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;
|
||||
|
||||
pimpl->bdry_condition = pimpl->iwork + max_ncf;
|
||||
|
||||
if (W != NULL) {
|
||||
pimpl->cwell_pos = pimpl->bdry_condition + nf;
|
||||
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)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
size_t i;
|
||||
|
||||
for (i = 0; i < nc + 1; i++) {
|
||||
pimpl->cwell_pos[i] = 0;
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
fsh_compute_table_sz(struct UnstructuredGrid *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 += G->number_of_faces; /* bdry_condition */
|
||||
*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 ];
|
||||
}
|
||||
}
|
||||
@ -1,90 +0,0 @@
|
||||
/*
|
||||
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 */
|
||||
|
||||
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 */
|
||||
};
|
||||
|
||||
|
||||
struct fsh_impl *
|
||||
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,
|
||||
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,
|
||||
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_linsys_arrays(struct fsh_data *h);
|
||||
|
||||
void
|
||||
fsh_compute_table_sz(struct UnstructuredGrid *G, well_t *W, int max_ngconn,
|
||||
size_t *nnu, size_t *idata_sz, size_t *ddata_sz);
|
||||
|
||||
#endif /* OPM_FSH_COMMON_IMPL_HEADER_INCLUDED */
|
||||
@ -1,390 +0,0 @@
|
||||
/*
|
||||
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 "config.h"
|
||||
#include <assert.h>
|
||||
#include <limits.h>
|
||||
#include <math.h>
|
||||
#include <stddef.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include <opm/core/pressure/fsh.h>
|
||||
#include <opm/core/pressure/fsh_common_impl.h>
|
||||
#include <opm/core/pressure/mimetic/hybsys.h>
|
||||
#include <opm/core/pressure/mimetic/hybsys_global.h>
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
ifsh_set_effective_well_params(const double *WI,
|
||||
const double *wdp,
|
||||
struct fsh_data *ifsh)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int c, nc, i, perf;
|
||||
int *cwpos, *cwells;
|
||||
double *wsys_WI, *wsys_wdp;
|
||||
|
||||
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];
|
||||
|
||||
wsys_WI [i] = WI [perf];
|
||||
wsys_wdp[i] = wdp[perf];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static int
|
||||
ifsh_assemble_grid(struct FlowBoundaryConditions *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;
|
||||
|
||||
nc = ifsh->pimpl->nc;
|
||||
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, gpress, src, Binv,
|
||||
ifsh->pimpl->sys);
|
||||
|
||||
npp += fsh_impose_bc(n, gconn + p1, bc, ifsh->pimpl);
|
||||
|
||||
hybsys_global_assemble_cell(n, gconn + p1,
|
||||
ifsh->pimpl->sys->S,
|
||||
ifsh->pimpl->sys->r,
|
||||
ifsh->A, ifsh->b);
|
||||
|
||||
p1 += n;
|
||||
p2 += n * n;
|
||||
}
|
||||
|
||||
return npp;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
ifsh_impose_well_control(int c,
|
||||
struct FlowBoundaryConditions *bc,
|
||||
well_control_t *wctrl,
|
||||
struct fsh_data *ifsh)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int ngconn, nwconn, i, j, w1, w2, wg, f;
|
||||
int *pgconn, *gconn, *pwconn, *wconn;
|
||||
|
||||
double bhp;
|
||||
double *r, *r2w, *w2w;
|
||||
|
||||
/* Enforce symmetric system */
|
||||
assert (ifsh->pimpl->wsys->r2w == ifsh->pimpl->wsys->w2r);
|
||||
|
||||
pgconn = ifsh->pimpl->gdof_pos;
|
||||
pwconn = ifsh->pimpl->cwell_pos;
|
||||
|
||||
gconn = ifsh->pimpl->gdof + pgconn[c];
|
||||
wconn = ifsh->pimpl->cwells + 2*pwconn[c];
|
||||
|
||||
ngconn = pgconn[c + 1] - pgconn[c];
|
||||
nwconn = pwconn[c + 1] - pwconn[c];
|
||||
|
||||
r2w = ifsh->pimpl->wsys->r2w;
|
||||
w2w = ifsh->pimpl->wsys->w2w;
|
||||
r = ifsh->pimpl->wsys->r ;
|
||||
|
||||
/* 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 (j != -1) {
|
||||
for (w1 = 0; w1 < nwconn; w1++) {
|
||||
/* Eliminate prescribed (boundary) pressure value */
|
||||
r [ngconn + w1] -= r2w[i + w1*ngconn] * bc->value[j];
|
||||
r2w[i + w1*ngconn] = 0.0;
|
||||
}
|
||||
|
||||
r[i] = 0.0; /* RHS value handled in *reservoir* asm */
|
||||
}
|
||||
}
|
||||
|
||||
/* Adapt local system to prescribed well (bottom-hole) pressures;
|
||||
* w->r and w->w. */
|
||||
for (w1 = 0; w1 < nwconn; w1++) {
|
||||
wg = wconn[2*w1 + 0];
|
||||
|
||||
if (wctrl->ctrl[wg] == BHP) {
|
||||
bhp = wctrl->target[wg];
|
||||
|
||||
/* Well->reservoir */
|
||||
for (i = 0; i < ngconn; i++) {
|
||||
#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;
|
||||
}
|
||||
|
||||
/* Well->well */
|
||||
for (w2 = (w1 + 1) % nwconn; w2 != w1; w2 = (w2 + 1) % nwconn) {
|
||||
r [ngconn + w2] -= w2w[w2 + w1*nwconn] * bhp;
|
||||
w2w[w2 + w1*ngconn] = 0.0;
|
||||
w2w[w1 + w2*ngconn] = 0.0;
|
||||
}
|
||||
|
||||
/* Assemble final well equation of the form S*p_bh = S*p_bh^0 */
|
||||
assert (fabs(w2w[w1 * (nwconn + 1)]) > 0.0);
|
||||
|
||||
r[ngconn + w1] = w2w[w1 * (nwconn + 1)] * bhp;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static int
|
||||
ifsh_assemble_well(struct FlowBoundaryConditions *bc,
|
||||
well_control_t *wctrl,
|
||||
struct fsh_data *ifsh)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int npp;
|
||||
int ngconn, nwconn, c, nc, w;
|
||||
|
||||
int *pgconn, *gconn, *pwconn, *wconn;
|
||||
|
||||
nc = ifsh->pimpl->nc;
|
||||
|
||||
pgconn = ifsh->pimpl->gdof_pos;
|
||||
gconn = ifsh->pimpl->gdof;
|
||||
pwconn = ifsh->pimpl->cwell_pos;
|
||||
wconn = ifsh->pimpl->cwells;
|
||||
|
||||
for (c = 0; c < nc; c++) {
|
||||
ngconn = pgconn[c + 1] - pgconn[c];
|
||||
nwconn = pwconn[c + 1] - pwconn[c];
|
||||
|
||||
if (nwconn > 0) {
|
||||
hybsys_well_cellcontrib_symm(c, ngconn, pgconn[c],
|
||||
pwconn,
|
||||
ifsh->pimpl->WI,
|
||||
ifsh->pimpl->wdp,
|
||||
ifsh->pimpl->sys,
|
||||
ifsh->pimpl->wsys);
|
||||
|
||||
ifsh_impose_well_control(c, bc, wctrl, ifsh);
|
||||
|
||||
hybsys_global_assemble_well_sym(ifsh->pimpl->nf,
|
||||
ngconn, gconn + pgconn[c],
|
||||
nwconn, wconn + 2*pwconn[c] + 0,
|
||||
ifsh->pimpl->wsys->r2w,
|
||||
ifsh->pimpl->wsys->w2w,
|
||||
ifsh->pimpl->wsys->r,
|
||||
ifsh->A, ifsh->b);
|
||||
}
|
||||
}
|
||||
|
||||
npp = 0;
|
||||
for (w = 0; w < ifsh->pimpl->nw; w++) {
|
||||
if (wctrl->ctrl[w] == BHP) {
|
||||
npp += 1;
|
||||
} else if (wctrl->ctrl[w] == RATE) {
|
||||
/* Impose total rate constraint.
|
||||
*
|
||||
* Note sign resulting from ->target[w] denoting
|
||||
* *injection* flux. */
|
||||
ifsh->b[ifsh->pimpl->nf + w] -= - wctrl->target[w];
|
||||
}
|
||||
}
|
||||
|
||||
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 *
|
||||
ifsh_construct(struct UnstructuredGrid *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. 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, 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);
|
||||
|
||||
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
|
||||
*
|
||||
* ifsh->A * ifsh->x = ifsh->b
|
||||
*
|
||||
* from effective local inner product matrices Binv, effective gravity
|
||||
* pressure gpress, boundary conditions bc, and source terms src. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
ifsh_assemble(struct FlowBoundaryConditions *bc,
|
||||
const double *src,
|
||||
const double *Binv,
|
||||
const double *gpress,
|
||||
well_control_t *wctrl,
|
||||
const double *WI,
|
||||
const double *wdp,
|
||||
struct fsh_data *ifsh)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
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);
|
||||
|
||||
if (ifsh->pimpl->nw > 0) {
|
||||
ifsh_set_effective_well_params(WI, wdp, ifsh);
|
||||
|
||||
hybsys_well_schur_comp_symm(ifsh->pimpl->nc,
|
||||
ifsh->pimpl->cwell_pos,
|
||||
ifsh->pimpl->WI,
|
||||
ifsh->pimpl->sys,
|
||||
ifsh->pimpl->wsys);
|
||||
}
|
||||
|
||||
npp = ifsh_assemble_grid(bc, Binv, gpress, src, ifsh);
|
||||
|
||||
if (ifsh->pimpl->nw > 0) {
|
||||
npp += ifsh_assemble_well(bc, wctrl, ifsh);
|
||||
}
|
||||
|
||||
if (npp == 0) {
|
||||
ifsh->A->sa[0] *= 2; /* Remove zero eigenvalue */
|
||||
}
|
||||
}
|
||||
@ -1,102 +0,0 @@
|
||||
/*
|
||||
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 "config.h"
|
||||
#include <stdlib.h>
|
||||
|
||||
#include <opm/core/pressure/legacy_well.h>
|
||||
|
||||
|
||||
/* Release memory resources for cell->well mapping. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
deallocate_cell_wells(int *cwpos, int *cwells)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
free(cwells);
|
||||
free(cwpos);
|
||||
}
|
||||
|
||||
|
||||
/* Allocate memory resources for cell->well mapping.
|
||||
*
|
||||
* Returns 1 if successful and 0 if not. CSR array pair set to NULL
|
||||
* unless allocation succeeds. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
int
|
||||
allocate_cell_wells(int nc, well_t *W, int **cwpos, int **cwells)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int i, totwconn;
|
||||
|
||||
totwconn = W->well_connpos[W->number_of_wells];
|
||||
|
||||
*cwpos = malloc((nc + 1) * sizeof **cwpos );
|
||||
*cwells = malloc(2 * totwconn * sizeof **cwells);
|
||||
|
||||
if ((*cwpos == NULL) || (*cwells == NULL)) {
|
||||
deallocate_cell_wells(*cwpos, *cwells);
|
||||
|
||||
*cwpos = NULL;
|
||||
*cwells = NULL;
|
||||
|
||||
totwconn = 0;
|
||||
} else {
|
||||
for (i = 0; i < nc + 1; i++) {
|
||||
(*cwpos)[i] = 0;
|
||||
}
|
||||
}
|
||||
|
||||
return totwconn;
|
||||
}
|
||||
|
||||
|
||||
/* Derive cell->well mapping from well->cell (connection) mapping. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
derive_cell_wells(int nc, well_t *W, int *cwpos, int *cwells)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int i, w, *c, *connpos;
|
||||
|
||||
connpos = W->well_connpos;
|
||||
|
||||
c = W->well_cells;
|
||||
for (w = 0; w < W->number_of_wells; w++) {
|
||||
for (; c != W->well_cells + connpos[w + 1]; c++) {
|
||||
cwpos[*c + 1] += 1;
|
||||
}
|
||||
}
|
||||
|
||||
for (i = 1; i <= nc; i++) {
|
||||
cwpos[0] += cwpos[i];
|
||||
cwpos[i] = cwpos[0] - cwpos[i];
|
||||
}
|
||||
|
||||
cwpos[0] = 0;
|
||||
c = W->well_cells;
|
||||
for (w = 0; w < W->number_of_wells; w++) {
|
||||
for (; c != W->well_cells + connpos[w + 1]; c++) {
|
||||
cwells[ 2*cwpos[*c + 1] + 0 ] = w;
|
||||
cwells[ 2*cwpos[*c + 1] + 1 ] = c - W->well_cells;
|
||||
|
||||
cwpos[*c + 1] += 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -87,48 +87,6 @@ typedef struct LegacyWellCompletions well_t;
|
||||
*/
|
||||
typedef struct LegacyWellControls well_control_t;
|
||||
|
||||
/**
|
||||
* Allocate cell-to-well mapping (as a sparse array).
|
||||
*
|
||||
* @param[in] nc Total number of cells.
|
||||
* @param[in] W Well topology (well-to-cell mapping).
|
||||
* @param[out] cwpos Indirection array. Points to array of size
|
||||
* <CODE>nc + 1</CODE> if successful.
|
||||
* @param[out] cwells Cell-to-well mapping. Points to array
|
||||
* of size <CODE>W->well_connpos[
|
||||
* W->number_of_wells]</CODE> if successful.
|
||||
* @return Positive number (size of <CODE>*cwells</CODE>)
|
||||
* if successful. Zero in case of allocation failure.
|
||||
*/
|
||||
int
|
||||
allocate_cell_wells(int nc, well_t *W, int **cwpos, int **cwells);
|
||||
|
||||
/**
|
||||
* Dispose of memory resources allocated using function
|
||||
* allocate_cell_wells().
|
||||
*
|
||||
* Following a call to deallocate_cell_wells(), the input pointers
|
||||
* are no longer valid.
|
||||
*
|
||||
* @param[in,out] cvpos Cell-to-well start pointers.
|
||||
* @param[in,out] cwells Cell-to-well mapping.
|
||||
*/
|
||||
void
|
||||
deallocate_cell_wells(int *cvpos, int *cwells);
|
||||
|
||||
/**
|
||||
* Construct cell-to-well mapping (i.e., transpose the
|
||||
* well-to-cell mapping represented by <CODE>W->well_cells</CODE>).
|
||||
*
|
||||
* @param[in] nc Total number of cells.
|
||||
* @param[in] W Well topology (well-to-cell mapping).
|
||||
* @param[out] cwpos Cell-to-well start pointers.
|
||||
* @param[out] cwells Cell-to-well mapping.
|
||||
*/
|
||||
void
|
||||
derive_cell_wells(int nc, well_t *W, int *cwpos, int *cwells);
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
@ -1,694 +0,0 @@
|
||||
/*
|
||||
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 "config.h"
|
||||
#include <assert.h>
|
||||
#include <stddef.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include <opm/core/linalg/blas_lapack.h>
|
||||
#include <opm/core/pressure/mimetic/hybsys.h>
|
||||
|
||||
|
||||
#if defined(MAX)
|
||||
#undef MAX
|
||||
#endif
|
||||
|
||||
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
struct hybsys *
|
||||
hybsys_allocate_symm(int max_nconn, int nc, int nconn_tot)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
struct hybsys *new;
|
||||
|
||||
new = malloc(1 * sizeof *new);
|
||||
if (new != NULL) {
|
||||
new->one = malloc(max_nconn * sizeof *new->one);
|
||||
new->r = malloc(max_nconn * sizeof *new->r );
|
||||
new->S = malloc(max_nconn * max_nconn * sizeof *new->S );
|
||||
new->L = malloc(nc * sizeof *new->L );
|
||||
new->q = malloc(nc * sizeof *new->q );
|
||||
new->F1 = malloc(nconn_tot * sizeof *new->F1 );
|
||||
|
||||
if ((new->one == NULL) || (new->r == NULL) || (new->S == NULL) ||
|
||||
(new->L == NULL) || (new->q == NULL) || (new->F1 == NULL)) {
|
||||
hybsys_free(new);
|
||||
|
||||
new = NULL;
|
||||
} else {
|
||||
new->F2 = new->F1;
|
||||
}
|
||||
}
|
||||
|
||||
return new;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
struct hybsys *
|
||||
hybsys_allocate_unsymm(int max_nconn, int nc, int nconn_tot)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
struct hybsys *new;
|
||||
|
||||
new = hybsys_allocate_symm(max_nconn, nc, nconn_tot);
|
||||
|
||||
if (new != NULL) {
|
||||
new->F2 = malloc(nconn_tot * sizeof *new->F2);
|
||||
|
||||
if (new->F2 == NULL) {
|
||||
hybsys_free(new);
|
||||
new = NULL;
|
||||
}
|
||||
}
|
||||
|
||||
return new;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
hybsys_well_count_conn(int nc, const int *cwpos,
|
||||
int *max_nw, size_t *sum_nwc)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int c, nw;
|
||||
|
||||
*max_nw = 0;
|
||||
*sum_nwc = 0;
|
||||
|
||||
for (c = 0; c < nc; c++) {
|
||||
nw = cwpos[c + 1] - cwpos[c];
|
||||
|
||||
assert (nw >= 0);
|
||||
|
||||
*max_nw = MAX(*max_nw, nw);
|
||||
*sum_nwc += nw;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
struct hybsys_well *
|
||||
hybsys_well_allocate_symm(int max_nconn, int nc, int *cwpos)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int max_nw;
|
||||
size_t sum_nwc, alloc_sz;
|
||||
|
||||
struct hybsys_well *new;
|
||||
|
||||
assert (cwpos[nc] > cwpos[0]); /* Else no wells. */
|
||||
|
||||
new = malloc(1 * sizeof *new);
|
||||
|
||||
if (new != NULL) {
|
||||
hybsys_well_count_conn(nc, cwpos, &max_nw, &sum_nwc);
|
||||
|
||||
alloc_sz = sum_nwc; /* F1 */
|
||||
alloc_sz += max_nconn + max_nw; /* r */
|
||||
alloc_sz += max_nw * max_nconn; /* w2r */
|
||||
alloc_sz += max_nw * max_nw; /* w2w */
|
||||
|
||||
new->data = malloc(alloc_sz * sizeof *new->data);
|
||||
|
||||
if (new->data != NULL) {
|
||||
new->F1 = new->data;
|
||||
new->F2 = new->F1;
|
||||
|
||||
new->r = new->F2 + sum_nwc;
|
||||
|
||||
new->w2r = new->r + max_nconn + max_nw;
|
||||
new->r2w = new->w2r;
|
||||
|
||||
new->w2w = new->r2w + (max_nw * max_nconn);
|
||||
} else {
|
||||
hybsys_well_free(new);
|
||||
new = NULL;
|
||||
}
|
||||
}
|
||||
|
||||
return new;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
struct hybsys_well *
|
||||
hybsys_well_allocate_unsymm(int max_nconn, int nc, int *cwpos)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int max_nw;
|
||||
size_t sum_nwc, alloc_sz;
|
||||
|
||||
struct hybsys_well *new;
|
||||
|
||||
assert (cwpos[nc] > cwpos[0]); /* Else no wells. */
|
||||
|
||||
new = malloc(1 * sizeof *new);
|
||||
|
||||
if (new != NULL) {
|
||||
hybsys_well_count_conn(nc, cwpos, &max_nw, &sum_nwc);
|
||||
|
||||
alloc_sz = 2 * sum_nwc; /* F1, F2 */
|
||||
alloc_sz += max_nconn + max_nw; /* r */
|
||||
alloc_sz += 2 * max_nw * max_nconn; /* w2r, r2w */
|
||||
alloc_sz += max_nw * max_nw; /* w2w */
|
||||
|
||||
new->data = malloc(alloc_sz * sizeof *new->data);
|
||||
|
||||
if (new->data != NULL) {
|
||||
new->F1 = new->data;
|
||||
new->F2 = new->F1 + sum_nwc;
|
||||
|
||||
new->r = new->F2 + sum_nwc;
|
||||
|
||||
new->w2r = new->r + max_nconn + max_nw;
|
||||
new->r2w = new->w2r + (max_nw * max_nconn);
|
||||
|
||||
new->w2w = new->r2w + (max_nw * max_nconn);
|
||||
} else {
|
||||
hybsys_well_free(new);
|
||||
new = NULL;
|
||||
}
|
||||
}
|
||||
|
||||
return new;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
hybsys_free(struct hybsys *sys)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
if (sys != NULL) {
|
||||
if (sys->F2 != sys->F1) { free(sys->F2); } /* unsymmetric system */
|
||||
|
||||
free(sys->F1 );
|
||||
free(sys->q );
|
||||
free(sys->L );
|
||||
free(sys->S );
|
||||
free(sys->r );
|
||||
free(sys->one);
|
||||
}
|
||||
|
||||
free(sys);
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
hybsys_well_free(struct hybsys_well *wsys)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
if (wsys != NULL) {
|
||||
free(wsys->data);
|
||||
}
|
||||
|
||||
free(wsys);
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
hybsys_init(int max_nconn, struct hybsys *sys)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int i;
|
||||
|
||||
for (i = 0; i < max_nconn; i++) {
|
||||
sys->one[i] = 1.0;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
hybsys_schur_comp_symm(int nc, const int *pconn,
|
||||
const double *Binv, struct hybsys *sys)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int c, p1, p2, nconn;
|
||||
double a1, a2;
|
||||
|
||||
MAT_SIZE_T incx, incy;
|
||||
MAT_SIZE_T nrows, ncols, lda;
|
||||
|
||||
incx = incy = 1;
|
||||
p1 = p2 = 0;
|
||||
|
||||
for (c = 0; c < nc; c++) {
|
||||
p1 = pconn[c + 0];
|
||||
nconn = pconn[c + 1] - pconn[c];
|
||||
nrows = ncols = lda = nconn;
|
||||
|
||||
/* F <- C' * inv(B) == (inv(B) * ones(n,1))' in single cell */
|
||||
a1 = 1.0; a2 = 0.0;
|
||||
dgemv_("No Transpose" , &nrows, &ncols,
|
||||
&a1, &Binv[p2] , &lda, sys->one, &incx,
|
||||
&a2, &sys->F1[p1], &incy);
|
||||
|
||||
/* L <- C' * inv(B) * C == SUM(F) == ones(n,1)' * F */
|
||||
sys->L[c] = ddot_(&nrows, sys->one, &incx, &sys->F1[p1], &incy);
|
||||
|
||||
p2 += nconn * nconn;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
hybsys_schur_comp_unsymm(int nc, const int *pconn,
|
||||
const double *Binv, const double *BIV,
|
||||
const double *P, struct hybsys *sys)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int c, p1, p2, nconn;
|
||||
double a1, a2;
|
||||
|
||||
MAT_SIZE_T incx, incy;
|
||||
MAT_SIZE_T nrows, ncols, lda;
|
||||
|
||||
assert ((sys->F2 != sys->F1) &&
|
||||
(sys->F2 != NULL));
|
||||
|
||||
incx = incy = 1;
|
||||
p2 = 0;
|
||||
|
||||
for (c = 0; c < nc; c++) {
|
||||
p1 = pconn[c + 0];
|
||||
nconn = pconn[c + 1] - pconn[c];
|
||||
|
||||
nrows = ncols = lda = nconn;
|
||||
|
||||
/* F1 <- C' * inv(B) */
|
||||
a1 = 1.0; a2 = 0.0;
|
||||
dgemv_("No Transpose" , &nrows, &ncols,
|
||||
&a1, &Binv[p2] , &lda, sys->one, &incx,
|
||||
&a2, &sys->F1[p1], &incy);
|
||||
|
||||
/* F2 <- (C - V)' * inv(B) == F1 - V'*inv(B) */
|
||||
a1 = -1.0;
|
||||
memcpy(&sys->F2[p1], &sys->F1[p1], nconn * sizeof sys->F2[p1]);
|
||||
daxpy_(&nrows, &a1, &BIV[p1], &incx, &sys->F2[p1], &incy);
|
||||
|
||||
/* L <- (C - V)' * inv(B) * C - P */
|
||||
sys->L[c] = ddot_(&nrows, sys->one, &incx, &sys->F1[p1], &incy);
|
||||
sys->L[c] -= ddot_(&nrows, sys->one, &incx, &BIV[p1] , &incy);
|
||||
sys->L[c] -= P[c];
|
||||
|
||||
p2 += nconn * nconn;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
hybsys_schur_comp_gen(int nc, const int *pconn,
|
||||
const double *Binv, const double *C2,
|
||||
const double *P, struct hybsys *sys)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int c, p1, p2, nconn;
|
||||
double a1, a2;
|
||||
|
||||
MAT_SIZE_T incx, incy;
|
||||
MAT_SIZE_T nrows, ncols, lda;
|
||||
|
||||
assert ((sys->F2 != sys->F1) &&
|
||||
(sys->F2 != NULL));
|
||||
|
||||
incx = incy = 1;
|
||||
p2 = 0;
|
||||
|
||||
for (c = 0; c < nc; c++) {
|
||||
p1 = pconn[c + 0];
|
||||
nconn = pconn[c + 1] - pconn[c];
|
||||
|
||||
nrows = ncols = lda = nconn;
|
||||
|
||||
/* F1 <- C' * inv(B) */
|
||||
a1 = 1.0; a2 = 0.0;
|
||||
dgemv_("No Transpose" , &nrows, &ncols,
|
||||
&a1, &Binv[p2] , &lda, sys->one, &incx,
|
||||
&a2, &sys->F1[p1], &incy);
|
||||
|
||||
/* F2 <- C2' * inv(B) */
|
||||
dgemv_("No Transpose" , &nrows, &ncols,
|
||||
&a1, &Binv[p2] , &lda, &C2[p1], &incx,
|
||||
&a2, &sys->F2[p1], &incy);
|
||||
|
||||
/* L <- C2' * inv(B) * C - P == F2'*ones(n,1) - P */
|
||||
sys->L[c] = ddot_(&nrows, sys->one, &incx, &sys->F2[p1], &incy);
|
||||
sys->L[c] -= P[c];
|
||||
|
||||
p2 += nconn * nconn;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
hybsys_well_schur_comp_symm(int nc, const int *cwpos,
|
||||
double *WI,
|
||||
struct hybsys *sys,
|
||||
struct hybsys_well *wsys)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int c, i;
|
||||
|
||||
for (c = i = 0; c < nc; c++) {
|
||||
for (; i < cwpos[c + 1]; i++) {
|
||||
wsys->F1[i] = WI[i];
|
||||
sys->L [c] += WI[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
hybsys_cellmat_symm_core(int nconn, const double *Binv, double L,
|
||||
const double *F, double *S)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int i, j;
|
||||
MAT_SIZE_T n, k, ldA, ldC;
|
||||
double a1, a2;
|
||||
|
||||
/* S <- D' * inv(B) * D == inv(B) in single cell */
|
||||
memcpy(S, Binv, nconn * nconn * sizeof *S);
|
||||
|
||||
/* S <- S - F'*inv(L)*F */
|
||||
n = ldA = ldC = nconn;
|
||||
k = 1;
|
||||
a1 = -1.0 / L;
|
||||
a2 = 1.0;
|
||||
dsyrk_("Upper Triangular", "No Transpose", &n, &k,
|
||||
&a1, F, &ldA, &a2, S, &ldC);
|
||||
|
||||
/* Account for DSYRK only updating the upper triangular part of S */
|
||||
for (j = 0; j < nconn; j++) {
|
||||
for (i = j + 1; i < nconn; i++) {
|
||||
S[i + j*nconn] = S[j + i*nconn];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
hybsys_cellmat_unsymm_core(int nconn, const double *Binv, double L,
|
||||
const double *F1, const double *F2,
|
||||
double *S)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
MAT_SIZE_T m, n, k, ldF1, ldF2, ldS;
|
||||
double a1, a2;
|
||||
|
||||
/* S <- D' * inv(B) * D == inv(B) in single cell */
|
||||
memcpy(S, Binv, nconn * nconn * sizeof *S);
|
||||
|
||||
/* S <- S - F1'*inv(L)*F2 */
|
||||
a1 = -1.0 / L;
|
||||
a2 = 1.0;
|
||||
|
||||
m = n = nconn;
|
||||
k = 1;
|
||||
ldF1 = ldF2 = 1;
|
||||
ldS = nconn;
|
||||
|
||||
dgemm_("Transpose", "No Transpose", &m, &n, &k,
|
||||
&a1, F1, &ldF1, F2, &ldF2, &a2, S, &ldS);
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static double
|
||||
hybsys_cellrhs_core(int nconn, const double *gpress, double src,
|
||||
const double *Binv, double L, const double *F1,
|
||||
const double *F2, double *R)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
MAT_SIZE_T n, k, ldA, incx, incy;
|
||||
double a1, a2;
|
||||
|
||||
/* r <- inv(B)*gpress + F1'*inv(L)*(src - F2*gpress)
|
||||
* == inv(B)*gpress + F1'*inv(L)*(src - C2'*inv(B)*gpress) */
|
||||
k = 1;
|
||||
a1 = 1.0; a2 = 0.0;
|
||||
incx = incy = 1;
|
||||
|
||||
n = k = ldA = nconn;
|
||||
|
||||
dgemv_("No Transpose", &n, &k,
|
||||
&a1, Binv, &ldA, gpress, &incx,
|
||||
&a2, R , &incy);
|
||||
|
||||
src -= ddot_(&n, F2, &incx, gpress, &incy);
|
||||
|
||||
a1 = src / L;
|
||||
daxpy_(&n, &a1, F1, &incx, R, &incy);
|
||||
|
||||
return src;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
hybsys_cellcontrib_symm(int c, int nconn, int p1, int p2,
|
||||
const double *gpress, const double *src,
|
||||
const double *Binv, struct hybsys *sys)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
hybsys_cellmat_symm_core(nconn, &Binv[p2],
|
||||
sys->L[c], &sys->F1[p1],
|
||||
sys->S);
|
||||
|
||||
sys->q[c] = hybsys_cellrhs_core(nconn, &gpress[p1], src[c], &Binv[p2],
|
||||
sys->L[c], &sys->F1[p1], &sys->F1[p1],
|
||||
sys->r);
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
hybsys_cellcontrib_unsymm(int c, int nconn, int p1, int p2,
|
||||
const double *gpress, const double *src,
|
||||
const double *Binv, struct hybsys *sys)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
assert ((sys->F2 != sys->F1) &&
|
||||
(sys->F2 != NULL));
|
||||
|
||||
hybsys_cellmat_unsymm_core(nconn, &Binv[p2],
|
||||
sys->L[c], &sys->F1[p1], &sys->F2[p1],
|
||||
sys->S);
|
||||
|
||||
sys->q[c] = hybsys_cellrhs_core(nconn, &gpress[p1], src[c], &Binv[p2],
|
||||
sys->L[c], &sys->F1[p1], &sys->F2[p1],
|
||||
sys->r);
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
hybsys_well_cellcontrib_symm(int c, int ngconn, int p1,
|
||||
const int *cwpos,
|
||||
const double *WI, const double *wdp,
|
||||
struct hybsys *sys, struct hybsys_well *wsys)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int i, w, nw, wp1;
|
||||
MAT_SIZE_T mm, nn, kk, ld1, ld2, ld3, incx, incy;
|
||||
double a1, a2, q;
|
||||
|
||||
nw = cwpos[c + 1] - cwpos[c];
|
||||
wp1 = cwpos[c];
|
||||
|
||||
/* -------------------------------------------------------------- */
|
||||
/* w2r = - F1(r)'*F2(w)/L, r2w = w2r' */
|
||||
mm = ngconn; ld1 = 1;
|
||||
nn = nw; ld2 = 1;
|
||||
kk = 1; ld3 = ngconn;
|
||||
|
||||
a1 = -1.0 / sys->L[c];
|
||||
a2 = 0.0;
|
||||
|
||||
dgemm_("Transpose", "No Transpose", &mm, &nn, &kk,
|
||||
&a1, &sys->F1[p1], &ld1, &wsys->F2[wp1], &ld2,
|
||||
&a2, wsys->w2r, &ld3);
|
||||
|
||||
/* -------------------------------------------------------------- */
|
||||
/* w2w = BI - F1(w)'*F2(w)/L */
|
||||
mm = nw; ld1 = 1;
|
||||
nn = nw; ld2 = 1;
|
||||
kk = 1; ld3 = nw;
|
||||
|
||||
a1 = -1.0 / sys->L[c];
|
||||
a2 = 0.0;
|
||||
|
||||
dgemm_("Transpose", "No Transpose", &mm, &nn, &kk,
|
||||
&a1, &wsys->F1[wp1], &ld1, &wsys->F2[wp1], &ld2,
|
||||
&a2, wsys->w2w, &ld3);
|
||||
|
||||
for (w = 0; w < nw; w++) {
|
||||
wsys->w2w[w * (nw + 1)] += WI[wp1 + w];
|
||||
}
|
||||
|
||||
/* -------------------------------------------------------------- */
|
||||
/* Global RHS contributions */
|
||||
mm = nw;
|
||||
incx = incy = 1;
|
||||
q = ddot_(&mm, &wsys->F2[wp1], &incx, &wdp[wp1], &incy);
|
||||
|
||||
a1 = -q / sys->L[c];
|
||||
for (i = 0; i < ngconn; i++) {
|
||||
wsys->r[i] = a1 * sys->F1[p1 + i];
|
||||
}
|
||||
|
||||
sys->q[c] -= q;
|
||||
a1 = sys->q[c] / sys->L[c];
|
||||
for (w = 0; w < nw; w++) {
|
||||
wsys->r[ngconn + w] = a1*wsys->F1[wp1 + w] +
|
||||
WI[wp1 + w] * wdp[wp1 + w];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
hybsys_compute_press_flux(int nc, const int *pconn, const int *conn,
|
||||
const double *gpress,
|
||||
const double *Binv, const struct hybsys *sys,
|
||||
const double *pi, double *press, double *flux,
|
||||
double *work)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int c, i, nconn, p1, p2;
|
||||
double a1, a2;
|
||||
|
||||
MAT_SIZE_T incx, incy, nrows, ncols, lda;
|
||||
|
||||
incx = incy = 1;
|
||||
|
||||
p2 = 0;
|
||||
a1 = 1.0;
|
||||
a2 = 0.0;
|
||||
for (c = 0; c < nc; c++) {
|
||||
p1 = pconn[c + 0];
|
||||
nconn = pconn[c + 1] - p1;
|
||||
|
||||
/* Serialise interface pressures for cell */
|
||||
for (i = 0; i < nconn; i++) {
|
||||
/* work[i] = pi[conn[p1 + i]] - gpress[p1 + i]; */
|
||||
work[i] = pi[conn[p1 + i]];
|
||||
}
|
||||
|
||||
nrows = ncols = lda = nconn;
|
||||
|
||||
/* Solve Lp = g - F2*f + F2*pi (for cell pressure) */
|
||||
press[c] = sys->q[c]; /* src[c]; */
|
||||
press[c] += ddot_(&nrows, &sys->F2[p1], &incx, work, &incy);
|
||||
press[c] /= sys->L[c];
|
||||
|
||||
/* Form rhs of system B*v = f + C*p - D*pi */
|
||||
for (i = 0; i < nconn; i++) {
|
||||
work[i] = gpress[p1 + i] + press[c] - work[i];
|
||||
}
|
||||
|
||||
/* Solve resulting system (-> half face fluxes) */
|
||||
dgemv_("No Transpose", &nrows, &ncols,
|
||||
&a1, &Binv[p2], &lda, work, &incx,
|
||||
&a2, &flux[p1], &incy);
|
||||
|
||||
p2 += nconn * nconn;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
hybsys_compute_press_flux_well(int nc, const int *pgconn, int nf,
|
||||
int nw, const int *pwconn, const int *wconn,
|
||||
const double *Binv,
|
||||
const double *WI,
|
||||
const double *wdp,
|
||||
const struct hybsys *sys,
|
||||
const struct hybsys_well *wsys,
|
||||
const double *pi,
|
||||
double *cpress, double *cflux,
|
||||
double *wpress, double *wflux,
|
||||
double *work)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int c, w, wg, perf;
|
||||
int ngconn, nwconn;
|
||||
size_t gp1, gp2, wp1;
|
||||
|
||||
MAT_SIZE_T mm, nn, incx, incy, ld;
|
||||
|
||||
double dcp, one;
|
||||
|
||||
gp2 = 0;
|
||||
for (c = 0; c < nc; c++) {
|
||||
ngconn = pgconn[c + 1] - pgconn[c];
|
||||
nwconn = pwconn[c + 1] - pwconn[c];
|
||||
|
||||
if (nwconn > 0) {
|
||||
dcp = 0.0;
|
||||
|
||||
gp1 = pgconn[c];
|
||||
wp1 = pwconn[c];
|
||||
|
||||
for (w = 0; w < nwconn; w++) {
|
||||
wg = wconn[2*(wp1 + w) + 0];
|
||||
work[w] = pi[nf + wg];
|
||||
}
|
||||
|
||||
mm = nwconn; incx = incy = 1;
|
||||
dcp = ddot_(&mm, &wsys->F2[wp1], &incx, work, &incy);
|
||||
dcp /= sys->L[c];
|
||||
|
||||
cpress[c] += dcp;
|
||||
|
||||
mm = nn = ld = ngconn;
|
||||
one = 1.0;
|
||||
dgemv_("No Transpose", &mm, &nn,
|
||||
&dcp, &Binv[gp2], &ld, sys->one , &incx,
|
||||
&one, &cflux[gp1], &incy);
|
||||
|
||||
for (w = 0; w < nwconn; w++) {
|
||||
perf = wconn[2*(wp1 + w) + 1];
|
||||
|
||||
wflux[perf] = wdp[wp1 + w] + cpress[c] - work[w];
|
||||
wflux[perf] *= - WI [wp1 + w]; /* Sign => positive inj. */
|
||||
}
|
||||
}
|
||||
|
||||
gp2 += ngconn + ngconn;
|
||||
}
|
||||
|
||||
/* Assign well BHP from linsolve output */
|
||||
memcpy(wpress, pi + nf, nw * sizeof *wpress);
|
||||
}
|
||||
@ -1,619 +0,0 @@
|
||||
/*
|
||||
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_HYBSYS_HEADER_INCLUDED
|
||||
#define OPM_HYBSYS_HEADER_INCLUDED
|
||||
|
||||
/**
|
||||
* \file
|
||||
* Routines and data structures to manage local contributions to a
|
||||
* global system of simultaneous linear equations arising from a
|
||||
* Schur complement reduction of an original block system.
|
||||
*
|
||||
* Specifically, these data structures and related routines compute
|
||||
* and store the elemental (cell-based) contributions of the Schur
|
||||
* complement reduction of the block system of simultaneous linear
|
||||
* equations
|
||||
* \f[
|
||||
* \begin{pmatrix}
|
||||
* B & C_1 & D \\
|
||||
* C_2^\mathsf{T} & P & 0 \\
|
||||
* D^\mathsf{T} & 0 & 0
|
||||
* \end{pmatrix}
|
||||
* \begin{pmatrix}
|
||||
* v \\ -p \\ \pi
|
||||
* \end{pmatrix} = \begin{pmatrix}
|
||||
* G \\ g \\ h
|
||||
* \end{pmatrix}
|
||||
* \f]
|
||||
* in which \f$G\f$ accounts for effects of gravity. The traditional
|
||||
* Schurcomplement reduction (block Gaussian elimination) then produces
|
||||
* the equivalent system of simultaneous linear equations
|
||||
* \f[
|
||||
* \begin{pmatrix}
|
||||
* B & C_1 & D \\
|
||||
* 0 & -L & -F_2 \\
|
||||
* 0 & 0 & A
|
||||
* \end{pmatrix}
|
||||
* \begin{pmatrix}
|
||||
* v \\ -p \\ \pi
|
||||
* \end{pmatrix} = \begin{pmatrix}
|
||||
* G \\ g - C_2^\mathsf{T}B^{-1}G \\ b
|
||||
* \end{pmatrix}.
|
||||
* \f]
|
||||
* Here, the matrix \f$A\f$ and the right hand side vector \f$b\f$ are given
|
||||
* by
|
||||
* \f[
|
||||
* \begin{aligned}
|
||||
* A &= D^\mathsf{T}B^{-1}D - F_1^\mathsf{T}L^{-1}F_2 \\
|
||||
* b &= D^\mathsf{T}B^{-1}G +
|
||||
* F_1^\mathsf{T}L^{-1}(g - C_2^\mathsf{T}B^{-1}G) - h,
|
||||
* \end{aligned}
|
||||
* \f]
|
||||
* and the component matrices \f$F_1\f$, \f$F_2\f$, and \f$L\f$ are given
|
||||
* by
|
||||
* \f[
|
||||
* F_1 = C_1^\mathsf{T}B^{-1}D, \quad
|
||||
* F_2 = C_2^\mathsf{T}B^{-1}D, \quad
|
||||
* L = C_2^\mathsf{T}B^{-1}C_1 - P.
|
||||
* \f]
|
||||
* In the case of incompressible flow, the matrix \f$C_2\f$ is the same
|
||||
* as \f$C_1\f$ and \f$P=0\f$ whence the coefficient matrix \f$A\f$ of
|
||||
* the Schur complement system \f$A\pi=b\f$ is symmetric.
|
||||
*
|
||||
* A great deal of simplification arises from the simple characterisation
|
||||
* of the \f$C_1\f$ and \f$D\f$ matrices. Specifically,
|
||||
* \f[
|
||||
* (C_1)_{ij} = \begin{cases}
|
||||
* 1, &\quad i\in\{\mathit{pconn}_j, \dots, \mathit{pconn}_{j+1}-1\}, \\
|
||||
* 0, &\quad \text{otherwise},
|
||||
* \end{cases}
|
||||
* \f]
|
||||
* and
|
||||
* \f[
|
||||
* (D)_{ij} = \begin{cases}
|
||||
* 1, &\quad \mathit{conn}_i = j, \\
|
||||
* 0, &\quad \text{otherwise}.
|
||||
* \end{cases}
|
||||
* \f]
|
||||
* When viewed in the context of a single cell, then the \f$D\f$ matrix
|
||||
* is, effectively, the identity with the \f$\mathit{conn}\f$ array
|
||||
* simply affecting a finite-element style redistribution (assembly)
|
||||
* of the local contributions. This module leverages that property
|
||||
* extensively.
|
||||
*/
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
/**
|
||||
* Elemental contributions (from cells) to block system of simultaneous
|
||||
* linear equations. Mixes quantities of single cells (@c r,
|
||||
* @c S and @c one) with those that cater to all cells (@c L,
|
||||
* @c q, @c F1, and--possibly--@c F2).
|
||||
*/
|
||||
struct hybsys {
|
||||
double *L; /**< \f$C_2^\mathsf{T}B^{-1}C - P\f$, all cells */
|
||||
double *q; /**< \f$g - F_2 G\f$, all cells */
|
||||
double *F1; /**< \f$C_1^\mathsf{T}B^{-1}\f$, all cells */
|
||||
double *F2; /**< \f$C_2^\mathsf{T}B^{-1}\f$, all cells*/
|
||||
double *r; /**< Data buffer for system right-hand side, single cell */
|
||||
double *S; /**< Data buffer system matrix, single cell */
|
||||
double *one; /**< \f$(1,1,\dots,1)^\mathsf{T}\f$, single cell */
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
* Elemental contributions (from wells) to block system of simultaneous
|
||||
* linear equations. Mixes quantities of single cell connections (@c r,
|
||||
* @c w2r, @c r2w, and @c w2w) and those that pertain to all well
|
||||
* connections (perforations) in concert (@c F1 and @c F2).
|
||||
*/
|
||||
struct hybsys_well {
|
||||
double *F1; /**< \f$C_1^\mathsf{T}B^{-1}\f$, all connections. */
|
||||
double *F2; /**< \f$C_2^\mathsf{T}B^{-1}\f$, all connections. */
|
||||
double *r; /**< Data buffer for system right-hand side, single cell. */
|
||||
|
||||
double *w2r; /**< Well-to-reservoir connection strength, single cell. */
|
||||
double *r2w; /**< Reservoir-to-well connection strength, single cell. */
|
||||
double *w2w; /**< Aggregate well-to-well connection strength. */
|
||||
|
||||
double *data; /**< Linear storage array. Structure undisclosed. */
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
* Allocate a hybrid system management structure suitable for discretising
|
||||
* a symmetric (i.e., incompressible) flow problem on a grid model of
|
||||
* given size.
|
||||
*
|
||||
* @param[in] max_nconn Maximum number of single cell faces.
|
||||
* @param[in] nc Total number of grid cells.
|
||||
* @param[in] nconn_tot Aggregate number of cell faces for all cells.
|
||||
* @return Fully formed hybrid system management structure if successful or
|
||||
* @c NULL in case of allocation failure.
|
||||
*/
|
||||
struct hybsys *
|
||||
hybsys_allocate_symm(int max_nconn, int nc, int nconn_tot);
|
||||
|
||||
|
||||
/**
|
||||
* Allocate a hybrid system management structure suitable for discretising
|
||||
* an unsymmetric (i.e., compressible) flow problem on a grid model of
|
||||
* given size.
|
||||
*
|
||||
* @param[in] max_nconn Maximum number of single cell faces.
|
||||
* @param[in] nc Total number of grid cells.
|
||||
* @param[in] nconn_tot Aggregate number of cell faces for all cells.
|
||||
* @return Fully formed hybrid system management structure if successful or
|
||||
* @c NULL in case of allocation failure.
|
||||
*/
|
||||
struct hybsys *
|
||||
hybsys_allocate_unsymm(int max_nconn, int nc, int nconn_tot);
|
||||
|
||||
|
||||
/**
|
||||
* Allocate a hybrid system management structure suitable for discretising
|
||||
* an incompressible (i.e., symmetric) well flow problem on a grid model
|
||||
* of given size.
|
||||
*
|
||||
* @param[in] max_nconn Maximum number of single cell faces.
|
||||
* @param[in] nc Total number of grid cells.
|
||||
* @param[in] cwpos Indirection array that defines each cell's
|
||||
* connecting wells. Values typically computed
|
||||
* using function derive_cell_wells().
|
||||
* @return Fully formed hybrid system management structure if successful or
|
||||
* @c NULL in case of allocation failure.
|
||||
*/
|
||||
struct hybsys_well *
|
||||
hybsys_well_allocate_symm(int max_nconn, int nc, int *cwpos);
|
||||
|
||||
|
||||
/**
|
||||
* Allocate a hybrid system management structure suitable for discretising
|
||||
* a compressible (i.e., unsymmetric) well flow problem on a grid model
|
||||
* of given size.
|
||||
*
|
||||
* @param[in] max_nconn Maximum number of single cell faces.
|
||||
* @param[in] nc Total number of grid cells.
|
||||
* @param[in] cwpos Indirection array that defines each cell's
|
||||
* connecting wells. Values typically computed
|
||||
* using function derive_cell_wells().
|
||||
* @return Fully formed hybrid system management structure if successful
|
||||
* or @c NULL in case of allocation failure.
|
||||
*/
|
||||
struct hybsys_well *
|
||||
hybsys_well_allocate_unsymm(int max_nconn, int nc, int *cwpos);
|
||||
|
||||
|
||||
/**
|
||||
* Dispose of memory resources previously obtained through one of the
|
||||
* allocation functions, hybsys_allocate_symm() or
|
||||
* hybsys_allocate_unsymm().
|
||||
*
|
||||
* Following a call to hybsys_free(), the input pointer is no longer
|
||||
* valid. <CODE>hybsys_free(NULL)</CODE> does nothing.
|
||||
*
|
||||
* @param[in,out] sys Previously allocated hybrid system management
|
||||
* structure (or @c NULL).
|
||||
*/
|
||||
void
|
||||
hybsys_free(struct hybsys *sys);
|
||||
|
||||
/**
|
||||
* Dispose of memory resources previously obtained through one of the
|
||||
* allocation functions, hybsys_well_allocate_symm() or
|
||||
* hybsys_well_allocate_unsymm().
|
||||
*
|
||||
* Following a call to hybsys_well_free(), the input pointer is
|
||||
* no longer valid. <CODE>hybsys_well_free(NULL)</CODE> does nothing.
|
||||
*
|
||||
* @param[in,out] wsys Previously allocated hybrid system management
|
||||
* structure (or @c NULL).
|
||||
*/
|
||||
void
|
||||
hybsys_well_free(struct hybsys_well *wsys);
|
||||
|
||||
|
||||
/**
|
||||
* Perform post-construction dynamic initialisation of system
|
||||
* structure obtained from function hybsys_allocate_symm() or
|
||||
* hybsys_allocate_unsymm().
|
||||
*
|
||||
* @param[in] max_nconn Maximum number of single cell faces.
|
||||
* Must coincide with the equally named
|
||||
* parameter of functions hybsys_allocate_symm()
|
||||
* or hybsys_allocate_unsymm().
|
||||
* @param[in,out] sys Previously allocated hybrid system management
|
||||
* structure.
|
||||
*/
|
||||
void
|
||||
hybsys_init(int max_nconn, struct hybsys *sys);
|
||||
|
||||
/**
|
||||
* Compute elemental (per-cell) contributions to symmetric Schur
|
||||
* system of simultaneous linear equations.
|
||||
*
|
||||
* This function assumes that the coefficient matrix of the hybrid
|
||||
* system of linear equations is that of the introduction with the
|
||||
* additional provision that \f$C_1=C_2=C\f$ and that \f$P=0\f$.
|
||||
* In other words, this function assumes that the coefficient matrix
|
||||
* is of the form
|
||||
* \f[
|
||||
* \begin{pmatrix}
|
||||
* B & C & D \\
|
||||
* C^\mathsf{T} & 0 & 0 \\
|
||||
* D^\mathsf{T} & 0 & 0
|
||||
* \end{pmatrix}.
|
||||
* \f]
|
||||
* This function fills the @c F1 and @c L fields of the management
|
||||
* structure.
|
||||
*
|
||||
* @param[in] nc Total number of grid cells.
|
||||
* @param[in] pconn Cell-to-face start pointers.
|
||||
* @param[in] Binv Inverse inner product results, usually
|
||||
* computed using mim_ip_simple_all() and
|
||||
* mim_ip_mobility_update().
|
||||
* @param[in,out] sys Hybrid system management structure allocated
|
||||
* using hybsys_allocate_symm() and initialised
|
||||
* using hybsys_init().
|
||||
*/
|
||||
void
|
||||
hybsys_schur_comp_symm(int nc, const int *pconn,
|
||||
const double *Binv, struct hybsys *sys);
|
||||
|
||||
|
||||
/**
|
||||
* Compute elemental (per-cell) contributions to unsymmetric Schur
|
||||
* system of simultaneous linear equations.
|
||||
*
|
||||
* This function assumes that the coefficient matrix of the hybrid
|
||||
* system of linear equations is that of the introduction with the
|
||||
* additional provision that \f$C_2=C_1-V\f$. In other words, this
|
||||
* function assumes that the coefficient matrix is of the form
|
||||
* \f[
|
||||
* \begin{pmatrix}
|
||||
* B & C & D \\
|
||||
* (C-V)^\mathsf{T} & P & 0 \\
|
||||
* D^\mathsf{T} & 0 & 0
|
||||
* \end{pmatrix}.
|
||||
* \f]
|
||||
* This matrix arises in the ``\f$v^2\f$'' phase compressibility
|
||||
* formulation of the compressible black-oil model. This function
|
||||
* fills the @c F1, @c F2 and @c L fields of the management structure.
|
||||
*
|
||||
* @param[in] nc Total number of grid cells.
|
||||
* @param[in] pconn Cell-to-face start pointers.
|
||||
* @param[in] Binv Inverse inner product results, usually
|
||||
* computed using mim_ip_simple_all() and
|
||||
* mim_ip_mobility_update().
|
||||
* @param[in] BIV \f$B^{-1}v\f$ in which \f$v\f$ is the flux
|
||||
* field of a previous time step or non-linear
|
||||
* iteration.
|
||||
* @param[in] P Per cell compressible accumulation term. One
|
||||
* scalar per cell.
|
||||
* @param[in,out] sys Hybrid system management structure allocated
|
||||
* using hybsys_allocate_symm() and initialised
|
||||
* using hybsys_init().
|
||||
*/
|
||||
void
|
||||
hybsys_schur_comp_unsymm(int nc, const int *pconn,
|
||||
const double *Binv, const double *BIV,
|
||||
const double *P, struct hybsys *sys);
|
||||
|
||||
|
||||
/**
|
||||
* Compute elemental (per-cell) contributions to unsymmetric Schur
|
||||
* system of simultaneous linear equations.
|
||||
*
|
||||
* This function assumes that the coefficient matrix of the hybrid
|
||||
* system of linear equations is that of the introduction with no
|
||||
* additional provisions. In other words, this
|
||||
* function assumes that the coefficient matrix is of the form
|
||||
* \f[
|
||||
* \begin{pmatrix}
|
||||
* B & C_1 & D \\
|
||||
* C_2^\mathsf{T} & P & 0 \\
|
||||
* D^\mathsf{T} & 0 & 0
|
||||
* \end{pmatrix}.
|
||||
* \f]
|
||||
* This function fills the @c F1, @c F2 and @c L fields of
|
||||
* the management structure.
|
||||
*
|
||||
* @param[in] nc Total number of grid cells.
|
||||
* @param[in] pconn Cell-to-face start pointers.
|
||||
* @param[in] Binv Inverse inner product results, usually
|
||||
* computed using mim_ip_simple_all() and
|
||||
* mim_ip_mobility_update().
|
||||
* @param[in] C2 Explicit representation of the \f$C_2\f$
|
||||
* matrix as a linear array. Assumed to only
|
||||
* contain the (structurally) non-zero matrix
|
||||
* elements (that correspond to the non-zero
|
||||
* structure of \f$C_1\f$).
|
||||
* @param[in] P Per cell compressible accumulation term. One
|
||||
* scalar per cell.
|
||||
* @param[in,out] sys Hybrid system management structure allocated
|
||||
* using hybsys_allocate_symm() and initialised
|
||||
* using hybsys_init().
|
||||
*/
|
||||
void
|
||||
hybsys_schur_comp_gen(int nc, const int *pconn,
|
||||
const double *Binv, const double *C2,
|
||||
const double *P, struct hybsys *sys);
|
||||
|
||||
/**
|
||||
* Compute elemental contributions to global, symmetric system of
|
||||
* simultaneous linear equations from cell<->well connections.
|
||||
*
|
||||
* Specifically, for a well @c w intersecting a cell @c c, this function
|
||||
* computes the elemental contributions
|
||||
* \f[
|
||||
* (F_1)_{wc} = C_{wc}^\mathsf{T} B_{wc}^{-1} D_{wc} = \mathit{WI}_{wc}
|
||||
* \f]
|
||||
* and
|
||||
* \f[
|
||||
* L_{wc} = C_{wc}^\mathsf{T} B_{wc}^{-1} C_{wc} = \mathit{WI}_{wc}
|
||||
* \f]
|
||||
* and incorporates the contributions into the global system quantities
|
||||
* as appropriate.
|
||||
*
|
||||
* This function modifies <CODE>sys->L</CODE> and <CODE>wsys->F1</CODE>.
|
||||
*
|
||||
* @param[in] nc Total number of grid cells.
|
||||
* @param[in] cwpos Indirection array that defines each cell's
|
||||
* connecting wells. Values typically computed
|
||||
* using function derive_cell_wells().
|
||||
* @param[in] WI Peaceman well connection indices. Array of
|
||||
* size <CODE>cwpos[nc]</CODE>. Must incorporate
|
||||
* effects of multiple phases (i.e., total mobility)
|
||||
* if applicable.
|
||||
* @param[in,out] sys Hybrid system management structure allocated
|
||||
* using hybsys_allocate_symm() and initialised
|
||||
* using hybsys_init() and/or filled using function
|
||||
* hybsys_schur_comp_symm().
|
||||
* @param[in,out] wsys Hybrid well-system management structure obtained
|
||||
* from function hybsys_well_allocate_symm().
|
||||
*/
|
||||
void
|
||||
hybsys_well_schur_comp_symm(int nc, const int *cwpos,
|
||||
double *WI,
|
||||
struct hybsys *sys,
|
||||
struct hybsys_well *wsys);
|
||||
|
||||
/**
|
||||
* Compute final (symmetric) Schur complement contributions to
|
||||
* global system of simultaneous linear equations.
|
||||
*
|
||||
* This function forms the coefficient matrix
|
||||
* \f[
|
||||
* S_c = D^\mathsf{T}B_c^{-1}D - F_c^\mathsf{T}L_c^{-1}F_c
|
||||
* \f]
|
||||
* and similar right-hand side \f$r_c\f$ elemental contributions.
|
||||
* These values must be subsequently assembled into the global system
|
||||
* using function hybsys_global_assemble_cell() after imposing any
|
||||
* applicable boundary conditions.
|
||||
*
|
||||
* This function overwrites the fields @c S and @c r of the hybrid system
|
||||
* structure.
|
||||
*
|
||||
* @param[in] c Cell for which to compute local contributions.
|
||||
* @param[in] nconn Number of connections (faces) of cell @c c.
|
||||
* @param[in] p1 Start address (into @c gpress) of the gravity
|
||||
* contributions of cell @c c.
|
||||
* @param[in] p2 Start address (into @c Binv) of the inverse
|
||||
* inner product of cell @c c.
|
||||
* @param[in] gpress Gravity contributions of all cells. Must
|
||||
* include effects of multiple phases if applicable.
|
||||
* @param[in] src Explicit source terms for all cells.
|
||||
* @param[in] Binv Inverse inner products for all cells. Must
|
||||
* include effects of multiple phases if applicable.
|
||||
* @param[in,out] sys Hybrid system management structure allocated
|
||||
* using hybsys_allocate_symm() and initialised
|
||||
* using hybsys_init() and/or filled using function
|
||||
* hybsys_schur_comp_symm() and
|
||||
* hybsys_well_schur_comp_symm() if applicable.
|
||||
*/
|
||||
void
|
||||
hybsys_cellcontrib_symm(int c, int nconn, int p1, int p2,
|
||||
const double *gpress, const double *src,
|
||||
const double *Binv, struct hybsys *sys);
|
||||
|
||||
|
||||
/**
|
||||
* Compute final (non-symmetric) Schur complement contributions to
|
||||
* global system of simultaneous linear equations.
|
||||
*
|
||||
* This function forms the coefficient matrix
|
||||
* \f[
|
||||
* S_c = D^\mathsf{T}B_c^{-1}D - (F_1)_c^\mathsf{T}L_c^{-1}(F_2)_c
|
||||
* \f]
|
||||
* and similar right-hand side \f$r_c\f$ elemental contributions.
|
||||
* These values must be subsequently assembled into the global system
|
||||
* using function hybsys_global_assemble_cell() after imposing any
|
||||
* applicable boundary conditions.
|
||||
*
|
||||
* This function overwrites the fields @c S and @c r of the hybrid system
|
||||
* structure.
|
||||
*
|
||||
* @param[in] c Cell for which to compute local contributions.
|
||||
* @param[in] nconn Number of connections (faces) of cell @c c.
|
||||
* @param[in] p1 Start address (into @c gpress) of the gravity
|
||||
* contributions of cell @c c.
|
||||
* @param[in] p2 Start address (into @c Binv) of the inverse
|
||||
* inner product of cell @c c.
|
||||
* @param[in] gpress Gravity contributions of all cells. Must
|
||||
* include effects of multiple phases if applicable.
|
||||
* @param[in] src Explicit source terms for all cells.
|
||||
* @param[in] Binv Inverse inner products for all cells. Must
|
||||
* include effects of multiple phases if applicable.
|
||||
* @param[in,out] sys Hybrid system management structure allocated
|
||||
* using hybsys_allocate_symm() and initialised
|
||||
* using hybsys_init() and/or filled using functions
|
||||
* hybsys_schur_comp_unsymm() or hybsys_schur_comp_gen().
|
||||
*/
|
||||
void
|
||||
hybsys_cellcontrib_unsymm(int c, int nconn, int p1, int p2,
|
||||
const double *gpress, const double *src,
|
||||
const double *Binv, struct hybsys *sys);
|
||||
|
||||
|
||||
/**
|
||||
* Form elemental direct contributions to global system of simultaneous linear
|
||||
* equations from cell<->well interactions.
|
||||
*
|
||||
* Plays a role similar to function hybsys_cellcontrib_symm(), but for wells.
|
||||
*
|
||||
* @param[in] c Cell for which to compute cell<->well Schur complement
|
||||
* @param[in] ngconn Number of inter-cell connections (faces) of cell @c c.
|
||||
* @param[in] p1 Start index (into <CODE>sys->F1</CODE>) of cell @c c.
|
||||
* @param[in] cwpos Indirection array that defines each cell's connecting
|
||||
* wells. Must coincide with equally named parameter of
|
||||
* function hybsys_well_schur_comp_symm().
|
||||
* @param[in] WI Peaceman well connection indices. Array of
|
||||
* size <CODE>pwconn[nc]</CODE>. Must coincide with
|
||||
* equally named parameter of contribution function
|
||||
* hybsys_well_schur_comp_symm().
|
||||
* @param[in] wdp Well connection gravity pressure adjustments.
|
||||
* One scalar for each well connection in an array of size
|
||||
* <CODE>pwconn[nc]</CODE>.
|
||||
* @param[in,out] sys Hybrid system management structure filled using
|
||||
* functions hybsys_schur_comp_unsymm() or
|
||||
* hybsys_schur_comp_gen().
|
||||
* @param[in,out] wsys Hybrid well-system management structure filled using
|
||||
* function hybsys_well_schur_comp_symm().
|
||||
*/
|
||||
void
|
||||
hybsys_well_cellcontrib_symm(int c, int ngconn, int p1,
|
||||
const int *cwpos,
|
||||
const double *WI, const double *wdp,
|
||||
struct hybsys *sys, struct hybsys_well *wsys);
|
||||
|
||||
|
||||
/**
|
||||
* Recover cell pressures and outward fluxes (with respect to cells--i.e., the
|
||||
* ``half-face fluxes'') through back substitution after solving a symmetric
|
||||
* (i.e., incompressible) Schur complement system of simultaneous linear
|
||||
* equations.
|
||||
*
|
||||
* Specifically, given the solution \f$\pi\f$ to the global system of
|
||||
* simultaneous linear equations, \f$A\pi=b\f$, that arises as a result of the
|
||||
* Schur complement analysis, this function recovers the cell pressures \f$p\f$
|
||||
* and outward fluxes \f$v\f$ defined by
|
||||
* \f[
|
||||
* \begin{aligned}
|
||||
* Lp &= g - C_2^\mathsf{T}B^{-1}G + F_2\pi \\
|
||||
* Bv &= G + C_1p - D\pi
|
||||
* \end{aligned}.
|
||||
* \f]
|
||||
*
|
||||
* @param[in] nc Total number of grid cells.
|
||||
* @param[in] pconn Cell-to-face start pointers.
|
||||
* @param[in] conn Cell-to-face mapping.
|
||||
* @param[in] gpress Gravity contributions of all cells. Must coincide with
|
||||
* equally named parameter in calls to cell contribution
|
||||
* functions such as hybsys_cellcontrib_symm().
|
||||
* @param[in] Binv Inverse inner products for all cells. Must coincide
|
||||
* with equally named parameter in calls to contribution
|
||||
* functions such as hybsys_cellcontrib_symm().
|
||||
* @param[in] sys Hybrid system management structure coinciding with
|
||||
* equally named parameter in contribution functions such
|
||||
* as hybsys_cellcontrib_symm() or
|
||||
* hybsys_cellcontrib_unsymm().
|
||||
* @param[in] pi Solution (interface/contact pressure) obtained from
|
||||
* solving the global system \f$A\pi = b\f$.
|
||||
* @param[out] press Cell pressures, \f$p\f$. Array of size @c nc.
|
||||
* @param[out] flux Outward interface fluxes, \f$v\f$. Array of size
|
||||
* <CODE>pconn[nc]</CODE>.
|
||||
* @param[in,out] work Scratch array for temporary results. Array of size at
|
||||
* least \f$\max_c \{ \mathit{pconn}_{c + 1}
|
||||
* - \mathit{pconn}_c \} \f$.
|
||||
*/
|
||||
void
|
||||
hybsys_compute_press_flux(int nc, const int *pconn, const int *conn,
|
||||
const double *gpress,
|
||||
const double *Binv, const struct hybsys *sys,
|
||||
const double *pi, double *press, double *flux,
|
||||
double *work);
|
||||
|
||||
|
||||
/**
|
||||
* Recover well pressures (i.e., bottom-hole pressure values) and well
|
||||
* connection (perforation) fluxes.
|
||||
*
|
||||
* Specifically, this function performs the same role (i.e., back-substitution)
|
||||
* for wells as function hybsys_compute_press_flux() does for grid cells and
|
||||
* grid contacts (interfaces).
|
||||
*
|
||||
* @param[in] nc Total number of grid cells.
|
||||
* @param[in] pgconn Cell-to-face start pointers.
|
||||
* @param[in] nf Total number of grid faces.
|
||||
* @param[in] nw Total number of wells.
|
||||
* @param[in] pwconn Cell-to-well start pointers. If <CODE>nw > 0</CODE>,
|
||||
* then this parameter must coincide with the @c cwpos
|
||||
* array used in call to hybsys_well_schur_comp_symm().
|
||||
* @param[in] wconn Cell-to-well mapping.
|
||||
* @param[in] Binv Inverse inner products for all cells. Must coincide
|
||||
* with equally named parameter in calls to contribution
|
||||
* functions such as hybsys_well_cellcontrib_symm().
|
||||
* @param[in] WI Peaceman well connection indices. Array of
|
||||
* size <CODE>pwconn[nc]</CODE>. Must coincide with
|
||||
* equally named parameter of contribution function
|
||||
* hybsys_well_cellcontrib_symm().
|
||||
* @param[in] wdp Well connection gravity pressure adjustments.
|
||||
* @param[in] sys Hybrid system management structure coinciding with
|
||||
* equally named parameter in contribution functions such
|
||||
* as hybsys_cellcontrib_symm() and
|
||||
* hybsys_well_cellcontrib_symm().
|
||||
* @param[in] wsys Hybrid well-system management structure. Must coincide
|
||||
* with equally named paramter of contribution function
|
||||
* hybsys_well_cellcontrib_symm().
|
||||
* @param[in] pi Solution (interface/contact pressure and well BHPs)
|
||||
* obtained from solving the global system \f$A\pi = b\f$.
|
||||
* @param[in] cpress Cell pressures, \f$p\f$, obtained from a previous call
|
||||
* to function hybsys_compute_press_flux().
|
||||
* @param[in] cflux Outward fluxes, \f$v\f$, obtained from a previous call
|
||||
* to function hybsys_compute_press_flux().
|
||||
* @param[out] wpress Well (i.e., bottom-hole) pressures. Array of size
|
||||
* @c nw.
|
||||
* @param[out] wflux Well connection (perforation) fluxes. Array of size
|
||||
* <CODE>pwconn[nw]</CODE>.
|
||||
* @param[in,out] work Scratch array for storing intermediate results. Array
|
||||
* of size at least \f$\max_w \{ \mathit{pwconn}_{w + 1}
|
||||
* - \mathit{pwconn}_w\}\f$.
|
||||
*/
|
||||
void
|
||||
hybsys_compute_press_flux_well(int nc, const int *pgconn, int nf,
|
||||
int nw, const int *pwconn, const int *wconn,
|
||||
const double *Binv,
|
||||
const double *WI,
|
||||
const double *wdp,
|
||||
const struct hybsys *sys,
|
||||
const struct hybsys_well *wsys,
|
||||
const double *pi,
|
||||
double *cpress, double *cflux,
|
||||
double *wpress, double *wflux,
|
||||
double *work);
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
#endif /* OPM_HYBSYS_HEADER_INCLUDED */
|
||||
@ -1,418 +0,0 @@
|
||||
/*
|
||||
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 "config.h"
|
||||
#include <assert.h>
|
||||
#include <stddef.h>
|
||||
#include <stdlib.h>
|
||||
|
||||
/* ---> Layering violation! */
|
||||
#include <opm/core/pressure/msmfem/hash_set.h>
|
||||
/* <--- Layering violation! */
|
||||
|
||||
#include <opm/core/pressure/mimetic/hybsys_global.h>
|
||||
|
||||
|
||||
#if defined MAX
|
||||
#undef MAX
|
||||
#endif
|
||||
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
|
||||
|
||||
|
||||
/* Release memory resources for each individual well's DOF set (and
|
||||
* the aggregate well DOF set structure). */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
deallocate_well_dofset(size_t nw, struct hash_set **wia)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
size_t w;
|
||||
|
||||
if (wia != NULL) {
|
||||
for (w = 0; w < nw; w++) {
|
||||
hash_set_deallocate(wia[w]);
|
||||
}
|
||||
}
|
||||
|
||||
free(wia);
|
||||
}
|
||||
|
||||
|
||||
/* Allocate DOF set for each well.
|
||||
*
|
||||
* Returns fully created vector of W->number_of_wells DOF sets if
|
||||
* successful and NULL otherwise. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static struct hash_set **
|
||||
allocate_well_dofset(struct UnstructuredGrid *G, well_t *W)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int w, i, c, ok;
|
||||
int set_size;
|
||||
struct hash_set **wia;
|
||||
|
||||
wia = malloc(W->number_of_wells * sizeof *wia);
|
||||
|
||||
if (wia != NULL) {
|
||||
ok = 1;
|
||||
for (w = 0; ok && (w < W->number_of_wells); w++) {
|
||||
set_size = 1; /* Approximate expected set size */
|
||||
|
||||
for (i = W->well_connpos[w]; i < W->well_connpos[w + 1]; i++) {
|
||||
c = W->well_cells[i];
|
||||
|
||||
set_size += G->cell_facepos[c + 1] - G->cell_facepos[c];
|
||||
}
|
||||
|
||||
wia[w] = hash_set_allocate(set_size);
|
||||
|
||||
ok = wia[w] != NULL;
|
||||
}
|
||||
|
||||
if (!ok) {
|
||||
deallocate_well_dofset(w, wia);
|
||||
wia = NULL;
|
||||
}
|
||||
}
|
||||
|
||||
return wia;
|
||||
}
|
||||
|
||||
|
||||
/* Count number of coefficient matrix connections (non-zeros) per row
|
||||
* from grid contributions. Add self connections for all rows. Well
|
||||
* connection counts will be overwritten in count_conn_per_row_well()
|
||||
* if applicable. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
count_conn_per_row_grid(struct UnstructuredGrid *G, struct CSRMatrix *A)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int c, nc, *ia, *ja;
|
||||
size_t m;
|
||||
|
||||
nc = G->number_of_cells;
|
||||
ia = G->cell_facepos;
|
||||
ja = G->cell_faces;
|
||||
|
||||
A->ia[0] = 0;
|
||||
for (m = 0; m < A->m; m++) {
|
||||
A->ia[m + 1] = 1;
|
||||
}
|
||||
|
||||
for (c = 0; c < nc; c++) {
|
||||
for (; ja != G->cell_faces + ia[c + 1]; ja++) {
|
||||
A->ia[ *ja + 1 ] += ia[c + 1] - ia[c] - 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* Count number of coefficient matrix connections per row from well
|
||||
* contributions. Increment connection count for other-to-well,
|
||||
* define count for well-to-other.
|
||||
*
|
||||
* Fails if unable to insert a DOF into a DOF set.
|
||||
*
|
||||
* Returns 1 if successful, and zero otherwise. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static int
|
||||
count_conn_per_row_well(struct UnstructuredGrid *G, well_t *W,
|
||||
int *cwpos,
|
||||
int *cwells,
|
||||
struct hash_set **wia,
|
||||
struct CSRMatrix *A)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int c, w, nc, nf, nwdof, ok, i, dof, *ia, *ja, *wja;
|
||||
|
||||
size_t m;
|
||||
|
||||
nc = G->number_of_cells;
|
||||
nf = G->number_of_faces;
|
||||
ia = G->cell_facepos;
|
||||
wja = cwells;
|
||||
|
||||
ok = 1;
|
||||
for (c = 0; c < nc; c++) {
|
||||
for (; ok && (wja != cwells + 2*cwpos[c + 1]); wja += 2) {
|
||||
for ( ja = G->cell_faces + ia[c + 0];
|
||||
ok && (ja != G->cell_faces + ia[c + 1]); ja++) {
|
||||
dof = *ja;
|
||||
ok = hash_set_insert(dof, wia[*wja]) == dof;
|
||||
}
|
||||
|
||||
for (i = cwpos[c]; ok && (i < cwpos[c + 1]); i++) {
|
||||
dof = nf + cwells[2*i + 0];
|
||||
ok = hash_set_insert(dof, wia[*wja]) == dof;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (ok) {
|
||||
for (w = 0; w < W->number_of_wells; w++) {
|
||||
nwdof = 0;
|
||||
|
||||
for (m = 0; m < wia[w]->m; m++) {
|
||||
if ((dof = wia[w]->s[m]) >= 0) {
|
||||
A->ia[ dof + 1 ] += 1; /* face to well */
|
||||
nwdof += 1; /* well to face */
|
||||
}
|
||||
}
|
||||
|
||||
A->ia[ nf + w + 1 ] = nwdof;
|
||||
}
|
||||
}
|
||||
|
||||
return ok;
|
||||
}
|
||||
|
||||
|
||||
/* Fill self-connections (i.e., diagonal of coefficient matrix) for
|
||||
* all DOFs. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
fill_self_connections(struct CSRMatrix *A)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
size_t r;
|
||||
|
||||
for (r = 0; r < A->m; r++) {
|
||||
A->ja[ A->ia[r + 1] ++ ] = r;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* Fill self-to-other DOF connections (i.e., define 'ja') for grid. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
fill_grid_connections(struct UnstructuredGrid *G, struct CSRMatrix *A)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int c, i, j, n;
|
||||
int dof1, dof2;
|
||||
|
||||
int *ia, *ja;
|
||||
|
||||
ia = G->cell_facepos;
|
||||
ja = G->cell_faces;
|
||||
|
||||
for (c = 0; c < G->number_of_cells; c++) {
|
||||
n = ia[c + 1] - ia[c];
|
||||
|
||||
for (i = 0; i < n; i++) {
|
||||
dof1 = ja[ ia[c] + i ];
|
||||
|
||||
if (dof1 >= 0) {
|
||||
for (j = (i + 1) % n; j != i; j = (j + 1) % n) {
|
||||
dof2 = ja[ ia[c] + j ];
|
||||
|
||||
if (dof2 >= 0) {
|
||||
A->ja[ A->ia[dof1 + 1] ++ ] = dof2;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* Fill self-to-other and other-to-self DOF connections ('ja') for wells. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
fill_well_connections(int nf, int nw,
|
||||
struct hash_set **wia,
|
||||
struct CSRMatrix *A)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int w, dof;
|
||||
size_t i;
|
||||
|
||||
for (w = 0; w < nw; w++) {
|
||||
for (i = 0; i < wia[w]->m; i++) {
|
||||
dof = wia[w]->s[i];
|
||||
|
||||
if (dof >= 0) {
|
||||
if (dof < nf) { /* Connect face to well */
|
||||
A->ja[ A->ia[ dof + 1 ] ++ ] = nf + w;
|
||||
}
|
||||
|
||||
if (dof != nf + w) {/* Connect well to dof (avoid self) */
|
||||
A->ja[ A->ia[ nf + w + 1 ] ++ ] = dof;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* Define pressure system coefficient matrix sparsity structure
|
||||
* (A->ia, A->ja) from grid and well contributions. Allocate
|
||||
* coefficient matrix elements (A->sa).
|
||||
*
|
||||
* Returns fully defined CSR matrix structure if successful or NULL
|
||||
* otherwise. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
struct CSRMatrix *
|
||||
hybsys_define_globconn(struct UnstructuredGrid *G, well_t *W)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int nw, ok;
|
||||
int *cwell_pos = NULL, *cwells = NULL;
|
||||
|
||||
struct hash_set **wia;
|
||||
struct CSRMatrix *A;
|
||||
|
||||
assert (G != NULL);
|
||||
|
||||
ok = 1;
|
||||
nw = 0;
|
||||
wia = NULL;
|
||||
if (W != NULL) {
|
||||
ok = allocate_cell_wells(G->number_of_cells, W, &cwell_pos, &cwells);
|
||||
if (ok) {
|
||||
derive_cell_wells(G->number_of_cells, W, cwell_pos, cwells);
|
||||
}
|
||||
nw = ok ? W->number_of_wells : 0;
|
||||
|
||||
if (nw > 0) {
|
||||
wia = allocate_well_dofset(G, W);
|
||||
if (wia == NULL) { nw = 0; }
|
||||
}
|
||||
}
|
||||
|
||||
A = csrmatrix_new_count_nnz(G->number_of_faces + nw);
|
||||
|
||||
if (A != NULL) {
|
||||
count_conn_per_row_grid(G, A);
|
||||
|
||||
if (nw > 0) {
|
||||
assert (nw == W->number_of_wells);
|
||||
ok = count_conn_per_row_well(G, W, cwell_pos,
|
||||
cwells, wia, A) > 0;
|
||||
} else {
|
||||
ok = 1;
|
||||
}
|
||||
|
||||
ok = ok && (csrmatrix_new_elms_pushback(A) > 0);
|
||||
|
||||
if (ok) {
|
||||
fill_self_connections(A);
|
||||
fill_grid_connections(G, A);
|
||||
fill_well_connections(G->number_of_faces, nw, wia, A);
|
||||
|
||||
csrmatrix_sortrows(A);
|
||||
} else {
|
||||
csrmatrix_delete(A);
|
||||
A = NULL;
|
||||
}
|
||||
}
|
||||
|
||||
deallocate_well_dofset(nw, wia);
|
||||
deallocate_cell_wells(cwell_pos, cwells);
|
||||
|
||||
return A;
|
||||
}
|
||||
|
||||
|
||||
/* Assemble (hybrid) cell contributions into global system coefficient
|
||||
* matrix and right hand side. Traditional FEM assembly process.
|
||||
* Boundary conditions assumed enforced outside.
|
||||
*
|
||||
* Local coefficient matrix contributions assumed organised in row
|
||||
* major format (row index cycling the most rapidly--Fortran
|
||||
* conventions). Convention immaterial if matrix is symmetric. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
hybsys_global_assemble_cell(int nconn, int *conn,
|
||||
const double *S,
|
||||
const double *r,
|
||||
struct CSRMatrix *A,
|
||||
double *b)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int il, jl; /* local */
|
||||
size_t ig, jg; /* global */
|
||||
|
||||
for (il = 0; il < nconn; il++) {
|
||||
assert ((0 <= conn[il]) && ((size_t) conn[il] < A->m));
|
||||
|
||||
ig = conn[il];
|
||||
|
||||
for (jl = 0; jl < nconn; jl++) {
|
||||
jg = csrmatrix_elm_index(ig, conn[jl], A);
|
||||
|
||||
assert ((((size_t)(A->ia[ig])) <= jg) &&
|
||||
(jg < ((size_t)(A->ia[ig + 1]))));
|
||||
|
||||
A->sa[jg] += S[il + jl*nconn]; /* Row major per cell */
|
||||
}
|
||||
|
||||
b[ig] += r[il];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
hybsys_global_assemble_well_sym(int ngconn_tot,
|
||||
int ngconn, const int *gconn,
|
||||
int nwconn, const int *wconn,
|
||||
const double *r2w,
|
||||
const double *w2w,
|
||||
const double *r,
|
||||
struct CSRMatrix *A,
|
||||
double *b)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int il, wl1, wl2;
|
||||
size_t ig, wg1, wg2, jg, jw;
|
||||
|
||||
/* Global matrix contributions for this cell's wells */
|
||||
for (wl1 = 0; wl1 < nwconn; wl1++) {
|
||||
wg1 = ngconn_tot + wconn[2*wl1 + 0];
|
||||
|
||||
for (il = 0; il < ngconn; il++) {
|
||||
ig = gconn[il];
|
||||
|
||||
jw = csrmatrix_elm_index(ig, wg1, A);
|
||||
jg = csrmatrix_elm_index(wg1, ig, A);
|
||||
|
||||
A->sa[jw] += r2w[il + wl1*ngconn]; /* Row major per cell */
|
||||
A->sa[jg] += r2w[il + wl1*ngconn]; /* Symmetry */
|
||||
}
|
||||
|
||||
for (wl2 = 0; wl2 < nwconn; wl2++) {
|
||||
wg2 = ngconn_tot + wconn[2*wl2 + 0];
|
||||
|
||||
jw = csrmatrix_elm_index(wg1, wg2, A);
|
||||
|
||||
A->sa[jw] += w2w[wl1 + wl2*nwconn]; /* Row major per well */
|
||||
}
|
||||
}
|
||||
|
||||
/* Global right-hand side contributions */
|
||||
for (il = 0; il < ngconn; il++) {
|
||||
b[gconn[il]] += r[il];
|
||||
}
|
||||
|
||||
for (wl1 = 0; wl1 < nwconn; wl1++) {
|
||||
b[ngconn_tot + wconn[2*wl1 + 0]] += r[ngconn + wl1];
|
||||
}
|
||||
}
|
||||
@ -1,212 +0,0 @@
|
||||
/*
|
||||
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_HYBSYS_GLOBAL_HEADER_INCLUDED
|
||||
#define OPM_HYBSYS_GLOBAL_HEADER_INCLUDED
|
||||
|
||||
/**
|
||||
* \file
|
||||
* Routines to assist in the formation and assembly of a global
|
||||
* system of simultaneous linear equations derived from a Schur
|
||||
* complement reduction of an original hybrid block system.
|
||||
*
|
||||
* We assume that the original block system of linear equations
|
||||
* is given by
|
||||
* \f[
|
||||
* \begin{pmatrix}
|
||||
* B & C & D \\
|
||||
* C^\mathsf{T} & 0 & 0 \\
|
||||
* D^\mathsf{T} & 0 & 0
|
||||
* \end{pmatrix}
|
||||
* \begin{pmatrix}
|
||||
* v \\ -p \\ \pi
|
||||
* \end{pmatrix} = \begin{pmatrix}
|
||||
* f \\ g \\ h
|
||||
* \end{pmatrix}
|
||||
* \f]
|
||||
* in which the block matrices \f$C\f$ and \f$D\f$ are assumed to
|
||||
* have a particularly simple structure and the matrix \f$B\f$ is
|
||||
* block diagonal.
|
||||
*
|
||||
* The Schur complement reduction process (a block Gaussian elimination)
|
||||
* then produces the following block system of simultaneous linear
|
||||
* equations
|
||||
* \f[
|
||||
* \begin{pmatrix}
|
||||
* B & C & D \\
|
||||
* 0 & -L & -F \\
|
||||
* 0 & 0 & A
|
||||
* \end{pmatrix}
|
||||
* \begin{pmatrix}
|
||||
* v \\ -p \\ \pi
|
||||
* \end{pmatrix} = \begin{pmatrix}
|
||||
* f \\ g - C^\mathsf{T}B^{-1} f \\ b
|
||||
* \end{pmatrix}
|
||||
* \f]
|
||||
* in which
|
||||
* \f[
|
||||
* \begin{aligned}
|
||||
* A &= D^\mathsf{T}B^{-1}D - F^\mathsf{T}L^{-1}F \text{ and} \\
|
||||
* b &= D^\mathsf{T}B^{-1}f + F^\mathsf{T}L^{-1} (g - C^\mathsf{T}B^{-1}f) - h.
|
||||
* \end{aligned}
|
||||
* \f] The component matrices \f$F\f$
|
||||
* and \f$L\f$ are given by
|
||||
* \f[
|
||||
* \begin{aligned}
|
||||
* L &= C^\mathsf{T} B^{-1} C \\
|
||||
* F &= C^\mathsf{T} B^{-1} D.
|
||||
* \end{aligned}
|
||||
* \f]
|
||||
* The primary degrees of freedom, \f$\pi\f$, may then be recovered
|
||||
* by solving the Schur complement system
|
||||
* \f[A\pi = b\f]
|
||||
* from which the derived quantities \f$p\f$ and \f$v\f$ may be
|
||||
* computed through a back-substitution process.
|
||||
*
|
||||
* The functions in this module assist in the creation of the sparsity
|
||||
* pattern of matrix \f$A\f$ and in the global assembling of values
|
||||
* into the matrix \f$A\f$ and the right-hand side vector \f$b\f$.
|
||||
* Specifically, function hybsys_define_globconn() builds the
|
||||
* sparsity pattern of \f$A\f$ while functions hybsys_global_assemble_cell()
|
||||
* and hybsys_global_assemble_well_sym() perform the task of
|
||||
* adding matrix and right-hand side values from local contributions.
|
||||
*/
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/pressure/legacy_well.h>
|
||||
#include <opm/core/linalg/sparse_sys.h>
|
||||
|
||||
|
||||
/**
|
||||
* Construct sparse matrix capable of managing a (reduced) hybrid
|
||||
* system of simultaneous linear equations with one degree of freedom
|
||||
* for each grid interface and each well.
|
||||
*
|
||||
* The return value is suitable for use in pressure solvers based on
|
||||
* Schur-complement reductions such as the mimetic finite-difference
|
||||
* or multiscale mixed finite-element classes of discretisations. In
|
||||
* typical applications, the matrix will be cleared using function
|
||||
* csrmatrix_zero() and then filled using an assembly process similar
|
||||
* to the traditional finite-element algorithm.
|
||||
*
|
||||
* Functions hybsys_global_assemble_cell() and
|
||||
* hybsys_global_assemble_well_sys() may be used to assist such an
|
||||
* assembly process.
|
||||
*
|
||||
* @param[in] G Grid.
|
||||
* @param[in] W Well topology. @c NULL in a model without wells.
|
||||
* @return Fully formed and structurally consistent sparse matrix
|
||||
* whose number of rows (i.e., degrees-of-freedom) equals
|
||||
* the number of grid faces (<CODE>G->number_of_faces</CODE>)
|
||||
* plus the number of wells (<CODE>W->number_of_wells</CODE>).
|
||||
*/
|
||||
struct CSRMatrix *
|
||||
hybsys_define_globconn(struct UnstructuredGrid *G, well_t *W);
|
||||
|
||||
|
||||
/**
|
||||
* Assemble local contributions into global system of simultaneous
|
||||
* linear equations.
|
||||
*
|
||||
* The contributions will typically have been computed using
|
||||
* function hybsys_schur_comp_symm() and function
|
||||
* hybsys_cellcontrib_symm() or some of the related functions.
|
||||
*
|
||||
* @param[in] nconn Number of cell faces.
|
||||
* @param[in] l2g Local-to-global mapping of cell's primary
|
||||
* degrees of freedom (i.e., the faces).
|
||||
* @param[in] S Single cell local contribution to global
|
||||
* coefficient matrix. An
|
||||
* \f$\mathit{nconn}\times\mathit{nconn}\f$
|
||||
* dense matrix in column major (Fortran) order.
|
||||
* @param[in] r Single cell local contribution to global
|
||||
* system right-hand side. An
|
||||
* \f$\mathit{nconn}\times 1\f$ dense vector.
|
||||
* @param[in,out] A Global coefficient matrix (of Schur
|
||||
* complement system).
|
||||
* @param[in,out] b Global system right-hand side (of Schur
|
||||
* complement system).
|
||||
*/
|
||||
void
|
||||
hybsys_global_assemble_cell(int nconn, int *l2g,
|
||||
const double *S,
|
||||
const double *r,
|
||||
struct CSRMatrix *A,
|
||||
double *b);
|
||||
|
||||
|
||||
/**
|
||||
* Assemble local contributions from single cell's well connections
|
||||
* (perforations) into global system of simultaneous linear equations.
|
||||
*
|
||||
* This function assumes that the connection strength from cell to well
|
||||
* equals the connection strength from well to cell. In other words,
|
||||
* that the numerical values of the well contributions are symmetric.
|
||||
*
|
||||
* The contributions are typically computed using functions
|
||||
* hybsys_well_schur_comp_symm() and hybsys_well_cellcontrib_symm().
|
||||
*
|
||||
* @param[in] ngconn_tot Total number of grid connections.
|
||||
* Expected to equal
|
||||
* <CODE>G->number_of_faces</CODE>
|
||||
* when @c G is the grid used to form the
|
||||
* original matrix in
|
||||
* hybsys_define_globconn().
|
||||
* @param[in] ngconn Number of grid connections referenced by
|
||||
* given cell.
|
||||
* @param[in] gconn Actual grid connections (DOFs) referenced
|
||||
* by given cell. Pointer to @c ngconn
|
||||
* consecutive DOF indices.
|
||||
* @param[in] nwconn Number of well connections intersecting
|
||||
* given cell. Typically \f$\mathit{ngconn} = 1\f$.
|
||||
* @param[in] wconn Actual well connections (DOFs) intersecting
|
||||
* given cell. Pointer to @c nwconn consecutive
|
||||
* DOF indices.
|
||||
* @param[in] r2w Reservoir-to-well connection strengths.
|
||||
* @param[in] w2w Well-to-well-connection strenghts.
|
||||
* @param[in] r Single cell local contribution to global
|
||||
* system right-hand side. An
|
||||
* \f$(\mathit{ngconn} + \mathit{nwconn})
|
||||
* \times 1\f$ dense vector.
|
||||
* @param[in,out] A Global coefficient matrix (of Schur
|
||||
* complement system).
|
||||
* @param[in,out] b Global system right-hand side (of Schur
|
||||
* complement system).
|
||||
*/
|
||||
void
|
||||
hybsys_global_assemble_well_sym(int ngconn_tot,
|
||||
int ngconn, const int *gconn,
|
||||
int nwconn, const int *wconn,
|
||||
const double *r2w,
|
||||
const double *w2w,
|
||||
const double *r,
|
||||
struct CSRMatrix *A,
|
||||
double *b);
|
||||
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
#endif /* OPM_HYBSYS_GLOBAL_HEADER_INCLUDED */
|
||||
@ -1,626 +0,0 @@
|
||||
/*
|
||||
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 "config.h"
|
||||
#include <assert.h>
|
||||
#include <limits.h>
|
||||
#include <stddef.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include <opm/core/pressure/msmfem/hash_set.h>
|
||||
#include <opm/core/pressure/msmfem/coarse_conn.h>
|
||||
|
||||
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
|
||||
#define MIN(a,b) (-MAX(-(a), -(b)))
|
||||
|
||||
|
||||
/* ======================================================================
|
||||
* Data structures
|
||||
* ====================================================================== */
|
||||
|
||||
/* Individual block connection. */
|
||||
struct block_neighbour {
|
||||
int b; /* Neighbouring block */
|
||||
struct hash_set *fconns; /* Constituent connections */
|
||||
};
|
||||
|
||||
|
||||
/* Adjacency list of single block (directed graph) */
|
||||
struct block_neighbours {
|
||||
int nneigh; /* Number of neighbours. */
|
||||
int cpty; /* Neighbour capacity. */
|
||||
struct block_neighbour **neigh; /* Actual neighbours (sorted on neigh[i]->b) */
|
||||
};
|
||||
|
||||
|
||||
/* ======================================================================
|
||||
* Operations
|
||||
* ====================================================================== */
|
||||
|
||||
|
||||
/* Relase dynamic memory resources for single block neighbour 'bn'. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
block_neighbour_deallocate(struct block_neighbour *bn)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
if (bn != NULL) {
|
||||
hash_set_deallocate(bn->fconns);
|
||||
}
|
||||
|
||||
free(bn);
|
||||
}
|
||||
|
||||
|
||||
/* Construct empty block neighbour connection capable of holding
|
||||
* 'nconn' fine-scale connections (e.g., fine-scale interfaces).
|
||||
* The fine-scale table is not allocated unless nconn > 0. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static struct block_neighbour *
|
||||
block_neighbour_allocate(int nconn)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
struct block_neighbour *new;
|
||||
|
||||
new = malloc(1 * sizeof *new);
|
||||
if (new != NULL) {
|
||||
if (nconn > 0) {
|
||||
new->fconns = hash_set_allocate(nconn);
|
||||
|
||||
if (new->fconns != NULL) {
|
||||
new->b = INT_MIN;
|
||||
} else {
|
||||
block_neighbour_deallocate(new);
|
||||
new = NULL;
|
||||
}
|
||||
} else {
|
||||
new->b = INT_MIN;
|
||||
new->fconns = NULL;
|
||||
}
|
||||
}
|
||||
|
||||
return new;
|
||||
}
|
||||
|
||||
|
||||
/* Insert fine-scale connection 'fconn' into block neighbour
|
||||
* connection 'bn', but only if the bn->fconns table has been allocated. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static int
|
||||
block_neighbour_insert_fconn(int fconn, struct block_neighbour *bn)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int ret;
|
||||
|
||||
assert (bn != NULL);
|
||||
|
||||
ret = 0;
|
||||
|
||||
if (bn->fconns != NULL) {
|
||||
ret = hash_set_insert(fconn, bn->fconns);
|
||||
}
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
/* Relase dynamic memory resources for single-block adjacency list 'bns'. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
block_neighbours_deallocate(struct block_neighbours *bns)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int i;
|
||||
|
||||
if (bns != NULL) {
|
||||
if (bns->neigh != NULL) {
|
||||
for (i = bns->nneigh - 1; i >= 0; i--) {
|
||||
block_neighbour_deallocate(bns->neigh[i]);
|
||||
}
|
||||
}
|
||||
|
||||
free(bns->neigh);
|
||||
}
|
||||
|
||||
free(bns);
|
||||
}
|
||||
|
||||
|
||||
/* Allocate a single-block adjacency list capable of holding 'nneigh'
|
||||
* connections. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static struct block_neighbours *
|
||||
block_neighbours_allocate(int nneigh)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int i;
|
||||
struct block_neighbours *new;
|
||||
|
||||
new = malloc(1 * sizeof *new);
|
||||
|
||||
if (new != NULL) {
|
||||
if (nneigh > 0) {
|
||||
new->neigh = malloc(nneigh * sizeof *new->neigh);
|
||||
|
||||
if (new->neigh != NULL) {
|
||||
for (i = 0; i < nneigh; i++) { new->neigh[i] = NULL; }
|
||||
new->nneigh = 0;
|
||||
new->cpty = nneigh;
|
||||
} else {
|
||||
block_neighbours_deallocate(new);
|
||||
new = NULL;
|
||||
}
|
||||
} else {
|
||||
new->nneigh = 0;
|
||||
new->cpty = 0;
|
||||
new->neigh = NULL;
|
||||
}
|
||||
}
|
||||
|
||||
return new;
|
||||
}
|
||||
|
||||
|
||||
/* Increase size of single-block adjacency list 'bns' to hold 'nneigh'
|
||||
* coarse-scale connections. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static int
|
||||
block_neighbours_expand(int nneigh, struct block_neighbours *bns)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int ret;
|
||||
struct block_neighbour **neigh;
|
||||
|
||||
assert (bns != NULL);
|
||||
|
||||
neigh = realloc(bns->neigh, nneigh * sizeof *neigh);
|
||||
|
||||
if (neigh != NULL) {
|
||||
bns->neigh = neigh;
|
||||
bns->cpty = nneigh;
|
||||
|
||||
for (ret = bns->nneigh; ret < bns->cpty; ret++) {
|
||||
bns->neigh[ret] = NULL;
|
||||
}
|
||||
} else {
|
||||
ret = -1;
|
||||
}
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
/* Insert fine-scale connection 'fconn' into single-block adjacency
|
||||
* list 'bns' in slot corresponding to connection 'b'.
|
||||
*
|
||||
* New coarse-scale connections are assumed to hold 'expct_nconn'
|
||||
* fine-scale connections.*/
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static int
|
||||
block_neighbours_insert_neighbour(int b, int fconn, int expct_nconn,
|
||||
struct block_neighbours *bns)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int i, j, p, t, nmove, ret;
|
||||
|
||||
assert (bns != NULL);
|
||||
|
||||
ret = 1;
|
||||
if ((bns->neigh == NULL) || (bns->cpty == 0)) {
|
||||
ret = block_neighbours_expand(1, bns);
|
||||
}
|
||||
|
||||
if (ret == 1) {
|
||||
/* bns->neigh points to table containing at least one slot. */
|
||||
i = 0;
|
||||
j = bns->nneigh;
|
||||
|
||||
while (i < j) {
|
||||
p = (i + j) / 2;
|
||||
|
||||
assert (bns->neigh[p] != NULL);
|
||||
|
||||
t = bns->neigh[p]->b;
|
||||
|
||||
if (t < b) { i = p + 1; }
|
||||
else if (t > b) { j = p + 0; }
|
||||
else { i = j = p; }
|
||||
}
|
||||
|
||||
if ((i < bns->nneigh) &&
|
||||
(bns->neigh[i] != NULL) && (bns->neigh[i]->b == b)) {
|
||||
ret = block_neighbour_insert_fconn(fconn, bns->neigh[i]);
|
||||
} else {
|
||||
if (bns->nneigh == bns->cpty) {
|
||||
assert (bns->cpty >= 1);
|
||||
ret = block_neighbours_expand(2 * bns->cpty, bns);
|
||||
}
|
||||
|
||||
if (ret >= 0) {
|
||||
if (i < bns->nneigh) {
|
||||
nmove = bns->nneigh - i;
|
||||
|
||||
memmove(bns->neigh + i + 1, bns->neigh + i + 0,
|
||||
nmove * sizeof *bns->neigh);
|
||||
}
|
||||
|
||||
bns->neigh[i] = block_neighbour_allocate(expct_nconn);
|
||||
|
||||
if (bns->neigh[i] != NULL) {
|
||||
ret = block_neighbour_insert_fconn(fconn, bns->neigh[i]);
|
||||
|
||||
bns->neigh[i]->b = b;
|
||||
bns->nneigh += 1;
|
||||
} else {
|
||||
ret = -1;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
/* Count number of (presumably) contiguously numbered blocks
|
||||
* represented by partition vector 'p'. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static int
|
||||
count_blocks(int nc, const int *p)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int i, max_blk;
|
||||
|
||||
max_blk = -1;
|
||||
for (i = 0; i < nc; i++) {
|
||||
max_blk = MAX(max_blk, p[i]);
|
||||
}
|
||||
|
||||
return max_blk + 1;
|
||||
}
|
||||
|
||||
|
||||
/* Derive coarse-scale block faces from fine-scale neighbour-ship
|
||||
* definition 'neighbours' ('nfinef' connections) and partition vector
|
||||
* 'p' (representing 'nblk' coarse blocks). Inter-block faces keyed
|
||||
* off minimum block number if internal and valid block number if
|
||||
* external.
|
||||
*
|
||||
* Fine-scale constituents of each coarse face are computed if
|
||||
* 'expct_nconn' is positive, in which case 'expct_nconn' is
|
||||
* interpreted as the expected number of constituents in each coarse
|
||||
* face and used as an initial size of a hash_set.
|
||||
*
|
||||
* Return number of coarse faces if successful and -1 otherwise. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static int
|
||||
derive_block_faces(int nfinef, int nblk, int expct_nconn,
|
||||
const int *p, const int *neighbours,
|
||||
struct block_neighbours **bns)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int f, c1, b1, c2, b2, b_in, b_out;
|
||||
int ret;
|
||||
|
||||
ret = 0;
|
||||
for (f = 0; (f < nfinef) && (0 <= ret); f++) {
|
||||
c1 = neighbours[2*f + 0]; b1 = (c1 >= 0) ? p[c1] : -1;
|
||||
c2 = neighbours[2*f + 1]; b2 = (c2 >= 0) ? p[c2] : -1;
|
||||
|
||||
assert ((b1 >= 0) || (b2 >= 0));
|
||||
|
||||
if ((b1 >= 0) && (b2 >= 0)) {
|
||||
b_in = MIN(b1, b2);
|
||||
b_out = MAX(b1, b2);
|
||||
} else if (b1 >= 0) { /* (b2 == -1) */
|
||||
b_in = b1;
|
||||
b_out = b2;
|
||||
} else {/*(b2 >= 0) *//* (b1 == -1) */
|
||||
b_in = b2;
|
||||
b_out = b1;
|
||||
}
|
||||
|
||||
if (b_in != b_out) {
|
||||
/* Block boundary */
|
||||
if (bns[b_in] == NULL) {
|
||||
bns[b_in] = block_neighbours_allocate(1);
|
||||
}
|
||||
|
||||
if (bns[b_in] != NULL) {
|
||||
ret = block_neighbours_insert_neighbour(b_out, f,
|
||||
expct_nconn,
|
||||
bns[b_in]);
|
||||
} else {
|
||||
ret = -1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (ret >= 0) {
|
||||
ret = 0;
|
||||
|
||||
for (b1 = 0; b1 < nblk; b1++) {
|
||||
if (bns[b1] != NULL) {
|
||||
ret += bns[b1]->nneigh;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
/* Create coarse-scale neighbour-ship definition from block-to-block
|
||||
* connectivity information ('bns') keyed off block numbers. Set
|
||||
* start pointers for CSR push-back build mode.
|
||||
*
|
||||
* Cannot fail. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
coarse_topology_build_coarsef(int nblk, struct block_neighbours **bns,
|
||||
int *neighbours, int *blkfacepos,
|
||||
size_t *nblkf, size_t *nsubf)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int b, n, coarse_f;
|
||||
|
||||
coarse_f = 0;
|
||||
*nsubf = 0;
|
||||
|
||||
for (b = 0; b < nblk; b++) {
|
||||
if (bns[b] != NULL) {
|
||||
for (n = 0; n < bns[b]->nneigh; n++) {
|
||||
neighbours[2*coarse_f + 0] = b;
|
||||
neighbours[2*coarse_f + 1] = bns[b]->neigh[n]->b;
|
||||
|
||||
coarse_f += 1;
|
||||
blkfacepos[b + 1] += 1;
|
||||
|
||||
if (bns[b]->neigh[n]->b >= 0) {
|
||||
blkfacepos[bns[b]->neigh[n]->b + 1] += 1;
|
||||
}
|
||||
|
||||
if (bns[b]->neigh[n]->fconns != NULL) {
|
||||
*nsubf += hash_set_count_elms(bns[b]->neigh[n]->fconns);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* Derive start pointers */
|
||||
for (b = 1; b <= nblk; b++) {
|
||||
blkfacepos[0] += blkfacepos[b];
|
||||
blkfacepos[b] = blkfacepos[0] - blkfacepos[b];
|
||||
}
|
||||
*nblkf = blkfacepos[0];
|
||||
blkfacepos[0] = 0;
|
||||
}
|
||||
|
||||
|
||||
/* Create coarse-scale block-to-face mapping and, if requested,
|
||||
* coarse-scale constituent faces for each coarse face.
|
||||
*
|
||||
* Constituent faces requested if subfacepos and subfaces non-NULL.
|
||||
* In this case, all coarse faces must carry sub-face information.
|
||||
*
|
||||
* Returns 1 if successful (i.e., no sub-face information requested or
|
||||
* sub-face information requested and available for all coarse faces)
|
||||
* and zero otherwise. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static int
|
||||
coarse_topology_build_final(int ncoarse_f, int nblk,
|
||||
const int *neighbours,
|
||||
int *blkfacepos, int *blkfaces,
|
||||
struct block_neighbours **bns,
|
||||
int *subfacepos, int *subfaces)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int coarse_f, b1, b2, n, subpos, subface_valid = 1;
|
||||
size_t i;
|
||||
struct hash_set *set;
|
||||
|
||||
assert ((subfacepos == NULL) == (subfaces == NULL));
|
||||
|
||||
for (coarse_f = 0; coarse_f < ncoarse_f; coarse_f++) {
|
||||
b1 = neighbours[2*coarse_f + 0];
|
||||
b2 = neighbours[2*coarse_f + 1];
|
||||
|
||||
assert (b1 != b2);
|
||||
|
||||
if (b1 >= 0) { blkfaces[blkfacepos[b1 + 1] ++] = coarse_f; }
|
||||
if (b2 >= 0) { blkfaces[blkfacepos[b2 + 1] ++] = coarse_f; }
|
||||
}
|
||||
|
||||
if (subfacepos != NULL) {
|
||||
coarse_f = 0;
|
||||
subpos = 0;
|
||||
|
||||
for (b1 = 0; (b1 < nblk) && subface_valid; b1++) {
|
||||
if (bns[b1] != NULL) {
|
||||
for (n = 0; n < bns[b1]->nneigh; n++) {
|
||||
set = bns[b1]->neigh[n]->fconns;
|
||||
subface_valid = set != NULL;
|
||||
|
||||
if (subface_valid) {
|
||||
for (i = 0; i < set->m; i++) {
|
||||
if (set->s[i] != -1) {
|
||||
subfaces[subpos ++] = set->s[i];
|
||||
}
|
||||
}
|
||||
} else {
|
||||
break;
|
||||
}
|
||||
|
||||
subfacepos[++ coarse_f] = subpos;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return (subfacepos == NULL) || subface_valid;
|
||||
}
|
||||
|
||||
|
||||
/* Allocate and assemble coarse-grid structure from non-linear
|
||||
* block-to-block connection information keyed off block numbers. The
|
||||
* final coarse grid consists of 'ncoarse_f' coarse faces numbered
|
||||
* 0..ncoarse_f-1 and 'nblk' coarse blocks numbered 0..nblk-1.
|
||||
*
|
||||
* Returns fully assembled coarse-grid structure if successful or NULL
|
||||
* otherwise. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static struct coarse_topology *
|
||||
coarse_topology_build(int ncoarse_f, int nblk,
|
||||
struct block_neighbours **bns)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int i;
|
||||
int subface_valid;
|
||||
size_t nblkf, nsubf;
|
||||
struct coarse_topology *new;
|
||||
|
||||
new = malloc(1 * sizeof *new);
|
||||
if (new != NULL) {
|
||||
new->neighbours = malloc(2 * ncoarse_f * sizeof *new->neighbours);
|
||||
new->blkfacepos = malloc((nblk + 1) * sizeof *new->blkfacepos);
|
||||
|
||||
new->blkfaces = NULL;
|
||||
new->subfacepos = NULL;
|
||||
new->subfaces = NULL;
|
||||
|
||||
if ((new->neighbours == NULL) ||
|
||||
(new->blkfacepos == NULL)) {
|
||||
coarse_topology_destroy(new);
|
||||
new = NULL;
|
||||
} else {
|
||||
for (i = 0; i < 2 * ncoarse_f; i++) {
|
||||
new->neighbours[i] = INT_MIN;
|
||||
}
|
||||
for (i = 0; i < nblk + 1; i++) {
|
||||
new->blkfacepos[i] = 0;
|
||||
}
|
||||
|
||||
coarse_topology_build_coarsef(nblk, bns, new->neighbours,
|
||||
new->blkfacepos, &nblkf, &nsubf);
|
||||
|
||||
if (nsubf > 0) {
|
||||
new->subfacepos = malloc((ncoarse_f + 1) * sizeof *new->subfacepos);
|
||||
new->subfaces = malloc(nsubf * sizeof *new->subfaces);
|
||||
|
||||
if ((new->subfacepos == NULL) || (new->subfaces == NULL)) {
|
||||
free(new->subfaces); new->subfaces = NULL;
|
||||
free(new->subfacepos); new->subfacepos = NULL;
|
||||
} else {
|
||||
for (i = 0; i < ncoarse_f + 1; i++) {
|
||||
new->subfacepos[i] = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
new->blkfaces = malloc(nblkf * sizeof *new->blkfaces);
|
||||
|
||||
if (new->blkfaces == NULL) {
|
||||
coarse_topology_destroy(new);
|
||||
new = NULL;
|
||||
} else {
|
||||
subface_valid = coarse_topology_build_final(ncoarse_f, nblk,
|
||||
new->neighbours,
|
||||
new->blkfacepos,
|
||||
new->blkfaces,
|
||||
bns,
|
||||
new->subfacepos,
|
||||
new->subfaces);
|
||||
|
||||
if (!subface_valid) {
|
||||
free(new->subfaces); new->subfaces = NULL;
|
||||
free(new->subfacepos); new->subfacepos = NULL;
|
||||
} else {
|
||||
new->nblocks = nblk;
|
||||
new->nfaces = ncoarse_f;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return new;
|
||||
}
|
||||
|
||||
|
||||
/* Create coarse-grid topology structure from fine-scale
|
||||
* neighbour-ship definition 'neighbours' and partition vector 'p'.
|
||||
*
|
||||
* Returns fully allocated and assembled coarse-grid structure if
|
||||
* successful and NULL otherwise. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
struct coarse_topology *
|
||||
coarse_topology_create(int nc, int nf, int expct_nconn,
|
||||
const int *p, const int *neighbours)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int b, nblocks, ncoarse_f;
|
||||
|
||||
struct block_neighbours **bns;
|
||||
struct coarse_topology *topo;
|
||||
|
||||
nblocks = count_blocks(nc, p);
|
||||
|
||||
bns = malloc(nblocks * sizeof *bns);
|
||||
if (bns != NULL) {
|
||||
for (b = 0; b < nblocks; b++) {
|
||||
bns[b] = NULL;
|
||||
}
|
||||
|
||||
ncoarse_f = derive_block_faces(nf, nblocks, expct_nconn,
|
||||
p, neighbours, bns);
|
||||
|
||||
topo = coarse_topology_build(ncoarse_f, nblocks, bns);
|
||||
|
||||
for (b = 0; b < nblocks; b++) {
|
||||
block_neighbours_deallocate(bns[b]);
|
||||
}
|
||||
|
||||
free(bns);
|
||||
} else {
|
||||
topo = NULL;
|
||||
}
|
||||
|
||||
return topo;
|
||||
}
|
||||
|
||||
|
||||
/* Release memory resources for dynamically allocated coarse-grid
|
||||
* topology structure 't'. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
coarse_topology_destroy(struct coarse_topology *t)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
if (t != NULL) {
|
||||
free(t->subfaces);
|
||||
free(t->subfacepos);
|
||||
|
||||
free(t->blkfaces);
|
||||
free(t->blkfacepos);
|
||||
|
||||
free(t->neighbours);
|
||||
}
|
||||
|
||||
free(t);
|
||||
}
|
||||
@ -1,54 +0,0 @@
|
||||
/*
|
||||
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_COARSE_CONN_HEADER_INCLUDED
|
||||
#define OPM_COARSE_CONN_HEADER_INCLUDED
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
struct coarse_topology {
|
||||
int nblocks; /* Number of blocks */
|
||||
int nfaces; /* Number of faces */
|
||||
|
||||
int *neighbours; /* Neighbourship definition */
|
||||
|
||||
int *blkfacepos; /* Index into blkfaces */
|
||||
int *blkfaces; /* Faces per block */
|
||||
|
||||
int *subfacepos; /* Index into subfaces */
|
||||
int *subfaces; /* FS faces per coarse face */
|
||||
};
|
||||
|
||||
|
||||
struct coarse_topology *
|
||||
coarse_topology_create(int nc, int nf, int expct_nconn,
|
||||
const int *p, const int *neighbours);
|
||||
|
||||
|
||||
void
|
||||
coarse_topology_destroy(struct coarse_topology *t);
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
#endif /* OPM_COARSE_CONN_HEADER_INCLUDED */
|
||||
File diff suppressed because it is too large
Load Diff
@ -1,98 +0,0 @@
|
||||
/*
|
||||
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_COARSE_SYS_HEADER_INCLUDED
|
||||
#define OPM_COARSE_SYS_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/grid.h>
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
struct coarse_sys {
|
||||
int *dof2conn; /* Map dof->connection (coarse interface) */
|
||||
int *blkdof_pos; /* Start pointers to each block's dofs */
|
||||
int *basis_pos; /* Start pointers to each block's bf's */
|
||||
int *cell_ip_pos; /* Start pointers to each block's IP */
|
||||
|
||||
int *blkdof; /* Each block's dofs */
|
||||
double *basis; /* All basis functions */
|
||||
double *cell_ip; /* Fine-scale IP contributions */
|
||||
double *Binv; /* Coarse-scale inverse IP per block */
|
||||
};
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
struct coarse_topology;
|
||||
struct CSRMatrix;
|
||||
|
||||
typedef void (*LocalSolver)(struct CSRMatrix *A,
|
||||
double *b,
|
||||
double *x);
|
||||
|
||||
struct coarse_sys *
|
||||
coarse_sys_construct(struct UnstructuredGrid *g, const int *p,
|
||||
struct coarse_topology *ct,
|
||||
const double *perm,
|
||||
const double *src,
|
||||
const double *totmob,
|
||||
LocalSolver linsolve);
|
||||
|
||||
void
|
||||
coarse_sys_destroy(struct coarse_sys *sys);
|
||||
|
||||
void
|
||||
coarse_sys_compute_cell_ip(int nc,
|
||||
int max_nconn,
|
||||
int nb,
|
||||
const int *pconn,
|
||||
const double *Binv,
|
||||
const int *b2c_pos,
|
||||
const int *b2c,
|
||||
struct coarse_sys *sys);
|
||||
|
||||
void
|
||||
coarse_sys_compute_Binv(int nb,
|
||||
int max_bcells,
|
||||
const double *totmob,
|
||||
const int *b2c_pos,
|
||||
const int *b2c,
|
||||
struct coarse_sys *sys,
|
||||
double *work);
|
||||
|
||||
void
|
||||
coarse_sys_compute_fs_flux(struct UnstructuredGrid *g,
|
||||
struct coarse_topology *ct,
|
||||
struct coarse_sys *sys,
|
||||
const int *b2c_pos,
|
||||
const int *b2c,
|
||||
const double *v_c,
|
||||
double *flux,
|
||||
double *work);
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
#endif /* OPM_COARSE_SYS_HEADER_INCLUDED */
|
||||
@ -1,241 +0,0 @@
|
||||
/*
|
||||
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 "config.h"
|
||||
#include <assert.h>
|
||||
#include <limits.h>
|
||||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include <opm/core/pressure/msmfem/hash_set.h>
|
||||
|
||||
/* ======================================================================
|
||||
* Macros
|
||||
* ====================================================================== */
|
||||
#define GOLDEN_RAT (0.6180339887498949) /* (sqrt(5) - 1) / 2 */
|
||||
#define IS_POW2(x) (((x) & ((x) - 1)) == 0)
|
||||
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
|
||||
|
||||
|
||||
/* Define a hash array size (1<<p) capable of holding a set of size 'm' */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static size_t
|
||||
hash_set_size(size_t m)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
size_t i;
|
||||
|
||||
if (m == 0) {
|
||||
return 1;
|
||||
}
|
||||
|
||||
if (IS_POW2(m)) {
|
||||
return m;
|
||||
}
|
||||
|
||||
/* General case. Use next power of two. */
|
||||
/* Algorithm due to
|
||||
*
|
||||
* Warren Jr., Henry S. (2002). Hacker's Delight.
|
||||
* Addison Wesley. pp. 48. ISBN 978-0201914658
|
||||
*
|
||||
* by way of Wikipedia. */
|
||||
m -= 1;
|
||||
for (i = 1; i < CHAR_BIT * sizeof m; i <<= 1) {
|
||||
m = m | (m >> i);
|
||||
}
|
||||
|
||||
return m + 1;
|
||||
}
|
||||
|
||||
|
||||
/* Hash element 'k' into table of size 'm' (multiplication method) */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static size_t
|
||||
hash_set_idx(int k, size_t m)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
double x = fmod(k * GOLDEN_RAT, 1.0);
|
||||
double y = floor(m * x);
|
||||
return y;
|
||||
}
|
||||
|
||||
|
||||
/* Insert element 'k' into set 's' of size 'm'
|
||||
* (open addressing, double probing). */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static size_t
|
||||
hash_set_insert_core(int k, size_t m, int *s)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
size_t h1, h2, i, j;
|
||||
|
||||
assert ((0 < m) && (m < (size_t)(-1)));
|
||||
assert (IS_POW2(m));
|
||||
|
||||
j = h1 = hash_set_idx(k, m);
|
||||
assert (h1 < m);
|
||||
|
||||
if (s[j] == -1) { s[j] = k; }
|
||||
if (s[j] == k) { return j; }
|
||||
|
||||
/* Double hash probing. h2 relatively prime to 'm' */
|
||||
h2 = 2 * hash_set_idx(k, MAX(m >> 1, 1)) + 1;
|
||||
|
||||
for (i = 1; (s[j] != -1) && (s[j] != k) && (i < m); i++) {
|
||||
j += h2;
|
||||
j &= m - 1; /* Modulo m since IS_POW2(m). */
|
||||
}
|
||||
|
||||
if ((s[j] == -1) || (s[j] == k)) {
|
||||
s[j] = k; /* Possibly no-op. */
|
||||
} else {
|
||||
j = m + 1; /* Invalid. Caveat emptor. */
|
||||
}
|
||||
|
||||
return j;
|
||||
}
|
||||
|
||||
|
||||
/* Increase size of hash set 't' to hold 'm' elements whilst copying
|
||||
* existing elements. This is typically fairly expensive. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static int
|
||||
hash_set_expand(size_t m, struct hash_set *t)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int ret, *s, *p;
|
||||
size_t i;
|
||||
|
||||
assert (m > t->m);
|
||||
|
||||
s = malloc(m * sizeof *s);
|
||||
if (s != NULL) {
|
||||
for (i = 0; i < m; i++) { s[i] = -1; }
|
||||
|
||||
for (i = 0; i < t->m; i++) {
|
||||
ret = hash_set_insert_core(t->s[i], m, s);
|
||||
assert ((size_t) ret < m);
|
||||
}
|
||||
|
||||
p = t->s;
|
||||
t->s = s;
|
||||
t->m = m;
|
||||
|
||||
free(p);
|
||||
|
||||
ret = m;
|
||||
} else {
|
||||
ret = -1;
|
||||
}
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
/* Release dynamic memory resources for hash set 't'. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
hash_set_deallocate(struct hash_set *t)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
if (t != NULL) {
|
||||
free(t->s);
|
||||
}
|
||||
|
||||
free(t);
|
||||
}
|
||||
|
||||
|
||||
/* Construct an emtpy hash set capable of holding 'm' elements */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
struct hash_set *
|
||||
hash_set_allocate(int m)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
size_t i, sz;
|
||||
struct hash_set *new;
|
||||
|
||||
new = malloc(1 * sizeof *new);
|
||||
if (new != NULL) {
|
||||
sz = hash_set_size(m);
|
||||
new->s = malloc(sz * sizeof *new->s);
|
||||
|
||||
if (new->s == NULL) {
|
||||
hash_set_deallocate(new);
|
||||
|
||||
new = NULL;
|
||||
} else {
|
||||
for (i = 0; i < sz; i++) { new->s[i] = -1; }
|
||||
new->m = sz;
|
||||
}
|
||||
}
|
||||
|
||||
return new;
|
||||
}
|
||||
|
||||
|
||||
/* Insert element 'k' into hash set 't'. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
int
|
||||
hash_set_insert(int k, struct hash_set *t)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int ret;
|
||||
size_t i;
|
||||
|
||||
assert (k >= 0);
|
||||
assert (t != NULL);
|
||||
assert (IS_POW2(t->m));
|
||||
|
||||
i = hash_set_insert_core(k, t->m, t->s);
|
||||
if (i == t->m + 1) {
|
||||
/* Table full. Preferable an infrequent occurrence. Expand
|
||||
* table and re-insert key (if possible). */
|
||||
ret = hash_set_expand(t->m << 1, t);
|
||||
|
||||
if (ret > 0) {
|
||||
i = hash_set_insert_core(k, t->m, t->s);
|
||||
assert (i < t->m);
|
||||
|
||||
ret = k;
|
||||
}
|
||||
} else {
|
||||
ret = k;
|
||||
}
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
size_t
|
||||
hash_set_count_elms(const struct hash_set *set)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
size_t i, n;
|
||||
|
||||
n = 0;
|
||||
for (i = 0; i < set->m; i++) {
|
||||
n += set->s[i] != -1;
|
||||
}
|
||||
|
||||
return n;
|
||||
}
|
||||
@ -1,70 +0,0 @@
|
||||
/*
|
||||
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_HASH_SET_HEADER_INCLUDED
|
||||
#define OPM_HASH_SET_HEADER_INCLUDED
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
#include <stddef.h>
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
/* Poor-man's unordered set (ind. key insert/all key extract only). */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
struct hash_set {
|
||||
size_t m; /* Table/set capacity (1<<p for some p) */
|
||||
int *s; /* Set representation */
|
||||
};
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
/* Release dynamic memory resources for hash set 's'. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
hash_set_deallocate(struct hash_set *s);
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
/* Construct an emtpy hash set capable of holding 'm' elements */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
struct hash_set *
|
||||
hash_set_allocate(int m);
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
/* Insert element 'k' into hash set 's'. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
int
|
||||
hash_set_insert(int k, struct hash_set *s);
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
/* Count number of valid keys in a hash set. */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
size_t
|
||||
hash_set_count_elms(const struct hash_set *set);
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
#endif /* OPM_HASH_SET_HEADER_INCLUDED */
|
||||
@ -1,553 +0,0 @@
|
||||
/*
|
||||
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 "config.h"
|
||||
#include <assert.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include <opm/core/linalg/blas_lapack.h>
|
||||
#include <opm/core/linalg/sparse_sys.h>
|
||||
|
||||
#include <opm/core/pressure/mimetic/hybsys.h>
|
||||
#include <opm/core/pressure/mimetic/hybsys_global.h>
|
||||
#include <opm/core/pressure/mimetic/mimetic.h>
|
||||
|
||||
#include <opm/core/pressure/msmfem/coarse_conn.h>
|
||||
#include <opm/core/pressure/msmfem/coarse_sys.h>
|
||||
#include <opm/core/pressure/msmfem/partition.h>
|
||||
#include <opm/core/pressure/msmfem/ifsh_ms.h>
|
||||
|
||||
|
||||
#if defined(MAX)
|
||||
#undef MAX
|
||||
#endif
|
||||
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
|
||||
|
||||
|
||||
struct ifsh_ms_impl {
|
||||
int max_bcells; /* Maximum number of cells per block */
|
||||
size_t ntotdof; /* Total number of degrees of freedom */
|
||||
size_t max_nblkdof; /* Max_i ndof(i) */
|
||||
size_t sum_nblkdof2; /* Sum_i ndof(i)^2 */
|
||||
|
||||
double *gpress; /* Coarse gravity contributions */
|
||||
double *bsrc; /* Coarse-scale source terms */
|
||||
|
||||
double *bpress; /* Block pressure */
|
||||
double *hflux; /* Coarse-scale half-contact fluxes */
|
||||
double *flux; /* Coarse-scale interface fluxes */
|
||||
|
||||
double *work; /* Assembly and back subst. work array */
|
||||
|
||||
double *fs_hflux; /* Fine-scale half-contact fluxes */
|
||||
|
||||
int *p; /* Copy of partition vector */
|
||||
int *pb2c, *b2c; /* Bloc->cell mapping (inverse partition) */
|
||||
|
||||
struct hybsys *hsys;
|
||||
struct coarse_topology *ct;
|
||||
struct coarse_sys *sys;
|
||||
|
||||
/* Linear storage */
|
||||
double *ddata;
|
||||
int *idata;
|
||||
};
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
ifsh_ms_impl_destroy(struct ifsh_ms_impl *pimpl)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
if (pimpl != NULL) {
|
||||
free (pimpl->idata);
|
||||
free (pimpl->ddata);
|
||||
hybsys_free (pimpl->hsys );
|
||||
coarse_sys_destroy (pimpl->sys );
|
||||
coarse_topology_destroy(pimpl->ct );
|
||||
}
|
||||
|
||||
free(pimpl);
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static int
|
||||
max_block_cells(size_t nc, size_t nb, const int *p)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int ret, *cnt;
|
||||
size_t b, c;
|
||||
|
||||
ret = -1;
|
||||
|
||||
cnt = malloc(nb * sizeof *cnt);
|
||||
|
||||
if (cnt != NULL) {
|
||||
ret = 0;
|
||||
|
||||
for (b = 0; b < nb; b++) { cnt[b] = 0; }
|
||||
|
||||
for (c = 0; c < nc; c++) {
|
||||
cnt[p[c]] += 1;
|
||||
ret = MAX(ret, cnt[p[c]]);
|
||||
}
|
||||
}
|
||||
|
||||
free(cnt);
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
count_blkdof(struct ifsh_ms_impl *pimpl)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int b, ndof, p;
|
||||
|
||||
pimpl->ntotdof = 0;
|
||||
pimpl->max_nblkdof = 0;
|
||||
pimpl->sum_nblkdof2 = 0;
|
||||
|
||||
for (b = p = 0; b < pimpl->ct->nblocks; b++) {
|
||||
ndof = pimpl->sys->blkdof_pos[b + 1] - pimpl->sys->blkdof_pos[b];
|
||||
|
||||
for (; p < pimpl->sys->blkdof_pos[b + 1]; p++) {
|
||||
pimpl->ntotdof = MAX(pimpl->ntotdof,
|
||||
(size_t) pimpl->sys->blkdof[p]);
|
||||
}
|
||||
|
||||
pimpl->max_nblkdof = MAX(pimpl->max_nblkdof, (size_t) ndof);
|
||||
pimpl->sum_nblkdof2 += ndof * ndof;
|
||||
}
|
||||
|
||||
pimpl->ntotdof += 1; /* Account for zero-based numbering */
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static int
|
||||
ifsh_ms_vectors_construct(struct UnstructuredGrid *G, struct ifsh_ms_impl *pimpl)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
size_t nb, n_fs_gconn;
|
||||
size_t nblkdof_tot, max_pairs, work_sz;
|
||||
size_t alloc_sz;
|
||||
|
||||
nb = pimpl->ct->nblocks;
|
||||
nblkdof_tot = pimpl->sys->blkdof_pos[nb];
|
||||
max_pairs = pimpl->max_nblkdof * (pimpl->max_nblkdof + 1) / 2;
|
||||
n_fs_gconn = G->cell_facepos[ G->number_of_cells ];
|
||||
|
||||
work_sz = pimpl->max_bcells + max_pairs;
|
||||
|
||||
alloc_sz = pimpl->ntotdof; /* b */
|
||||
alloc_sz += pimpl->ntotdof; /* x */
|
||||
|
||||
alloc_sz += nblkdof_tot; /* gpress */
|
||||
alloc_sz += nb; /* bsrc */
|
||||
alloc_sz += nb; /* bpress */
|
||||
alloc_sz += nblkdof_tot; /* hflux */
|
||||
alloc_sz += pimpl->ntotdof; /* flux */
|
||||
|
||||
alloc_sz += work_sz; /* work */
|
||||
|
||||
alloc_sz += n_fs_gconn; /* fs_hflux */
|
||||
|
||||
pimpl->ddata = malloc(alloc_sz * sizeof *pimpl->ddata);
|
||||
|
||||
alloc_sz = G->number_of_cells; /* p */
|
||||
alloc_sz += nb + 1; /* pb2c */
|
||||
alloc_sz += G->number_of_cells; /* b2c */
|
||||
|
||||
pimpl->idata = malloc(alloc_sz * sizeof *pimpl->idata);
|
||||
|
||||
return (pimpl->ddata != NULL) && (pimpl->idata != NULL);
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
set_impl_pointers(struct UnstructuredGrid *G, struct ifsh_ms_data *h)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
size_t nb;
|
||||
size_t nblkdof_tot, max_pairs, work_sz;
|
||||
|
||||
nb = h->pimpl->ct->nblocks;
|
||||
nblkdof_tot = h->pimpl->sys->blkdof_pos[nb];
|
||||
max_pairs = h->pimpl->max_nblkdof * (h->pimpl->max_nblkdof + 1) / 2;
|
||||
|
||||
work_sz = h->pimpl->max_bcells + max_pairs;
|
||||
|
||||
/* Linear system */
|
||||
h->b = h->pimpl->ddata;
|
||||
h->x = h->b + h->pimpl->ntotdof;
|
||||
|
||||
/* Coarse-scale back-substitution results */
|
||||
h->pimpl->gpress = h->x + h->pimpl->ntotdof;
|
||||
h->pimpl->bsrc = h->pimpl->gpress + nblkdof_tot;
|
||||
h->pimpl->bpress = h->pimpl->bsrc + nb;
|
||||
h->pimpl->hflux = h->pimpl->bpress + nb;
|
||||
h->pimpl->flux = h->pimpl->hflux + nblkdof_tot;
|
||||
|
||||
/* Back-substitution work array */
|
||||
h->pimpl->work = h->pimpl->flux + h->pimpl->ntotdof;
|
||||
|
||||
/* Fine-scale hflux accumulation array */
|
||||
h->pimpl->fs_hflux = h->pimpl->work + work_sz;
|
||||
|
||||
/* Partition vector */
|
||||
h->pimpl->p = h->pimpl->idata;
|
||||
|
||||
/* Block->cell mapping (inverse partition vector) */
|
||||
h->pimpl->pb2c = h->pimpl->p + G->number_of_cells;
|
||||
h->pimpl->b2c = h->pimpl->pb2c + nb + 1;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static struct CSRMatrix *
|
||||
ifsh_ms_matrix_construct(size_t m, size_t nnz, size_t nb,
|
||||
const int *pdof, const int *dof)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int p, n, i1, dof1, i2, dof2;
|
||||
size_t i, b;
|
||||
struct CSRMatrix *new;
|
||||
|
||||
new = csrmatrix_new_known_nnz(m, nnz);
|
||||
|
||||
if (new != NULL) {
|
||||
for (i = 0; i < m; i++) { new->ia[ i + 1 ] = 1; } /* Count self */
|
||||
|
||||
/* Count others */
|
||||
for (b = 0, p = 0; b < nb; b++) {
|
||||
n = pdof[b + 1] - pdof[b];
|
||||
|
||||
for (; p < pdof[b + 1]; p++) {
|
||||
new->ia[ dof[p] + 1 ] += n - 1;
|
||||
}
|
||||
}
|
||||
|
||||
/* Set start pointers */
|
||||
new->ia[0] = 0;
|
||||
for (i = 1; i <= m; i++) {
|
||||
new->ia[0] += new->ia[i];
|
||||
new->ia[i] = new->ia[0] - new->ia[i];
|
||||
}
|
||||
|
||||
/* Insert self */
|
||||
for (i = 0; i < m; i++) {
|
||||
new->ja[ new->ia[ i + 1 ] ++ ] = (MAT_SIZE_T) i;
|
||||
}
|
||||
|
||||
/* Insert others */
|
||||
for (b = 0; b < nb; b++) {
|
||||
p = pdof[b];
|
||||
n = pdof[b + 1] - p;
|
||||
|
||||
for (i1 = 0; i1 < n; i1++) {
|
||||
dof1 = dof[p + i1];
|
||||
|
||||
for (i2 = (i1 + 1) % n; i2 != i1; i2 = (i2 + 1) % n) {
|
||||
dof2 = dof[p + i2];
|
||||
|
||||
new->ja[ new->ia[ dof1 + 1 ] ++ ] = dof2;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
assert (new->ia[m] <= (MAT_SIZE_T) nnz);
|
||||
|
||||
new->ia[0] = 0;
|
||||
csrmatrix_sortrows(new);
|
||||
}
|
||||
|
||||
return new;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static struct ifsh_ms_impl *
|
||||
ifsh_ms_impl_construct(struct UnstructuredGrid *G ,
|
||||
const int *p ,
|
||||
const double *perm ,
|
||||
const double *src ,
|
||||
const double *totmob,
|
||||
LocalSolver linsolve)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int max_nconn = -1, nb, nconn_tot;
|
||||
int expected_nconn, alloc_ok;
|
||||
|
||||
struct ifsh_ms_impl *new;
|
||||
|
||||
new = malloc(1 * sizeof *new);
|
||||
|
||||
if (new != NULL) {
|
||||
new->hsys = NULL; /* Crucial NULL-initialisation */
|
||||
new->sys = NULL;
|
||||
new->ddata = NULL;
|
||||
new->idata = NULL;
|
||||
|
||||
expected_nconn = 256; /* CHANGE THIS! */
|
||||
alloc_ok = 0;
|
||||
|
||||
new->ct = coarse_topology_create(G->number_of_cells,
|
||||
G->number_of_faces,
|
||||
expected_nconn, p,
|
||||
G->face_cells);
|
||||
|
||||
|
||||
if (new->ct != NULL) {
|
||||
new->max_bcells = max_block_cells(G->number_of_cells,
|
||||
new->ct->nblocks, p);
|
||||
|
||||
new->sys = coarse_sys_construct(G, p, new->ct, perm,
|
||||
src, totmob, linsolve);
|
||||
}
|
||||
|
||||
if ((new->sys != NULL) && (new->max_bcells > 0)) {
|
||||
count_blkdof(new);
|
||||
|
||||
max_nconn = new->max_nblkdof;
|
||||
nb = new->ct->nblocks;
|
||||
nconn_tot = new->sys->blkdof_pos[ nb ];
|
||||
|
||||
new->hsys = hybsys_allocate_symm(max_nconn, nb, nconn_tot);
|
||||
}
|
||||
|
||||
if (new->hsys != NULL) {
|
||||
alloc_ok = ifsh_ms_vectors_construct(G, new);
|
||||
}
|
||||
|
||||
if (! alloc_ok) {
|
||||
ifsh_ms_impl_destroy(new);
|
||||
new = NULL;
|
||||
} else {
|
||||
hybsys_init(max_nconn, new->hsys);
|
||||
}
|
||||
}
|
||||
|
||||
return new;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
static void
|
||||
average_flux(size_t nf, const int *N, double *flux) /* Poor name */
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
size_t f;
|
||||
|
||||
for (f = 0; f < nf; f++) {
|
||||
flux[f] /= ((N[2*f + 0] >= 0) + (N[2*f + 1] >= 0));
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ====================================================================== *
|
||||
* Public routines below. *
|
||||
* ====================================================================== */
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
struct ifsh_ms_data *
|
||||
ifsh_ms_construct(struct UnstructuredGrid *G ,
|
||||
const int *p ,
|
||||
const double *perm ,
|
||||
const double *src ,
|
||||
const double *totmob,
|
||||
LocalSolver linsolve)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int i;
|
||||
struct ifsh_ms_data *new;
|
||||
|
||||
new = malloc(1 * sizeof *new);
|
||||
|
||||
if (new != NULL) {
|
||||
new->pimpl = ifsh_ms_impl_construct(G, p, perm, src,
|
||||
totmob, linsolve);
|
||||
new->A = NULL;
|
||||
|
||||
if (new->pimpl != NULL) {
|
||||
new->A = ifsh_ms_matrix_construct(new->pimpl->ntotdof,
|
||||
new->pimpl->sum_nblkdof2,
|
||||
new->pimpl->ct->nblocks,
|
||||
new->pimpl->sys->blkdof_pos,
|
||||
new->pimpl->sys->blkdof);
|
||||
}
|
||||
|
||||
if (new->A == NULL) {
|
||||
ifsh_ms_destroy(new);
|
||||
new = NULL;
|
||||
} else {
|
||||
set_impl_pointers(G, new);
|
||||
|
||||
memcpy(new->pimpl->p, p, G->number_of_cells * sizeof *p);
|
||||
|
||||
for (i = 0; i < new->pimpl->ct->nblocks; i++) {
|
||||
new->pimpl->pb2c[i] = 0;
|
||||
}
|
||||
|
||||
partition_invert(G->number_of_cells, new->pimpl->p,
|
||||
new->pimpl->pb2c, new->pimpl->b2c);
|
||||
}
|
||||
}
|
||||
|
||||
return new;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
ifsh_ms_destroy(struct ifsh_ms_data *h)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
if (h != NULL) {
|
||||
ifsh_ms_impl_destroy(h->pimpl);
|
||||
csrmatrix_delete (h->A );
|
||||
}
|
||||
|
||||
free(h);
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
ifsh_ms_assemble(const double *src ,
|
||||
const double *totmob,
|
||||
struct ifsh_ms_data *h)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int b, i, nconn, p1, p2;
|
||||
|
||||
int *pb2c, *b2c;
|
||||
double *bsrc;
|
||||
|
||||
pb2c = h->pimpl->pb2c;
|
||||
b2c = h->pimpl->b2c;
|
||||
bsrc = h->pimpl->bsrc;
|
||||
|
||||
for (b = i = 0; b < h->pimpl->ct->nblocks; b++) {
|
||||
bsrc[b] = 0.0;
|
||||
|
||||
for (; i < pb2c[b + 1]; i++) { bsrc[b] += src[ b2c[i] ]; }
|
||||
}
|
||||
|
||||
coarse_sys_compute_Binv(h->pimpl->ct->nblocks,
|
||||
h->pimpl->max_bcells, totmob,
|
||||
pb2c, b2c, h->pimpl->sys, h->pimpl->work);
|
||||
|
||||
hybsys_schur_comp_symm(h->pimpl->ct->nblocks,
|
||||
h->pimpl->sys->blkdof_pos,
|
||||
h->pimpl->sys->Binv, h->pimpl->hsys);
|
||||
|
||||
csrmatrix_zero( h->A);
|
||||
vector_zero (h->A->m, h->b);
|
||||
|
||||
/* Exclude gravity */
|
||||
vector_zero(h->pimpl->sys->blkdof_pos[ h->pimpl->ct->nblocks ],
|
||||
h->pimpl->gpress);
|
||||
|
||||
p2 = 0;
|
||||
for (b = 0; b < h->pimpl->ct->nblocks; b++) {
|
||||
p1 = h->pimpl->sys->blkdof_pos[b + 0] ;
|
||||
nconn = h->pimpl->sys->blkdof_pos[b + 1] - p1;
|
||||
|
||||
hybsys_cellcontrib_symm(b, nconn, p1, p2, h->pimpl->gpress,
|
||||
bsrc, h->pimpl->sys->Binv,
|
||||
h->pimpl->hsys);
|
||||
|
||||
hybsys_global_assemble_cell(nconn, h->pimpl->sys->blkdof + p1,
|
||||
h->pimpl->hsys->S, h->pimpl->hsys->r,
|
||||
h->A, h->b);
|
||||
|
||||
p2 += nconn * nconn;
|
||||
}
|
||||
|
||||
/* Remove zero eigenvalue */
|
||||
h->A->sa[0] *= 2;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
ifsh_ms_press_flux(struct UnstructuredGrid *G, struct ifsh_ms_data *h,
|
||||
double *cpress, double *fflux)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int b, f, i, dof, *c;
|
||||
double s;
|
||||
|
||||
hybsys_compute_press_flux(h->pimpl->ct->nblocks,
|
||||
h->pimpl->sys->blkdof_pos,
|
||||
h->pimpl->sys->blkdof,
|
||||
h->pimpl->gpress, h->pimpl->sys->Binv,
|
||||
h->pimpl->hsys, h->x, h->pimpl->bpress,
|
||||
h->pimpl->hflux, h->pimpl->work);
|
||||
|
||||
vector_zero(h->pimpl->ct->nfaces, h->pimpl->flux);
|
||||
|
||||
for (b = i = 0; b < h->pimpl->ct->nblocks; b++) {
|
||||
for (; i < h->pimpl->sys->blkdof_pos[b + 1]; i++) {
|
||||
dof = h->pimpl->sys->blkdof[ i ];
|
||||
|
||||
f = h->pimpl->sys->dof2conn[ dof ];
|
||||
s = 2.0*(h->pimpl->ct->neighbours[2*f + 0] == b) - 1.0;
|
||||
|
||||
assert (f < h->pimpl->ct->nfaces);
|
||||
|
||||
h->pimpl->flux[ f ] += s * h->pimpl->hflux[ i ];
|
||||
}
|
||||
}
|
||||
|
||||
average_flux(h->pimpl->ct->nfaces,
|
||||
h->pimpl->ct->neighbours, h->pimpl->flux);
|
||||
|
||||
for (b = i = 0; b < h->pimpl->ct->nblocks; b++) {
|
||||
/* Derive coarse-scale (projected) hc fluxes for block */
|
||||
for (; i < h->pimpl->sys->blkdof_pos[b + 1]; i++) {
|
||||
dof = h->pimpl->sys->blkdof[ i ];
|
||||
|
||||
f = h->pimpl->sys->dof2conn[ dof ];
|
||||
s = 2.0*(h->pimpl->ct->neighbours[2*f + 0] == b) - 1.0;
|
||||
|
||||
h->pimpl->hflux[ i ] = s * h->pimpl->flux[ f ];
|
||||
}
|
||||
|
||||
/* Derive fine-scale pressure from block pressure */
|
||||
for (c = h->pimpl->b2c + h->pimpl->pb2c[b + 0];
|
||||
c != h->pimpl->b2c + h->pimpl->pb2c[b + 1]; c++) {
|
||||
|
||||
cpress[*c] = h->pimpl->bpress[b];
|
||||
}
|
||||
}
|
||||
|
||||
coarse_sys_compute_fs_flux(G, h->pimpl->ct, h->pimpl->sys,
|
||||
h->pimpl->pb2c, h->pimpl->b2c,
|
||||
h->pimpl->hflux,
|
||||
fflux, /* Output */
|
||||
h->pimpl->fs_hflux); /* Scratch array */
|
||||
}
|
||||
@ -1,71 +0,0 @@
|
||||
/*
|
||||
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_MS_HEADER_INCLUDED
|
||||
#define OPM_IFSH_MS_HEADER_INCLUDED
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
#include <stddef.h>
|
||||
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/pressure/msmfem/coarse_sys.h>
|
||||
|
||||
struct CSRMatrix;
|
||||
struct ifsh_ms_impl;
|
||||
|
||||
struct ifsh_ms_data {
|
||||
/* Linear system */
|
||||
struct CSRMatrix *A; /* Coefficient matrix */
|
||||
double *b; /* System RHS */
|
||||
double *x; /* Solution */
|
||||
|
||||
/* Private implementational details. */
|
||||
struct ifsh_ms_impl *pimpl;
|
||||
};
|
||||
|
||||
|
||||
struct ifsh_ms_data *
|
||||
ifsh_ms_construct(struct UnstructuredGrid *G,
|
||||
const int *p,
|
||||
const double *perm,
|
||||
const double *src,
|
||||
const double *totmob,
|
||||
LocalSolver linsolve);
|
||||
|
||||
void
|
||||
ifsh_ms_destroy(struct ifsh_ms_data *h);
|
||||
|
||||
void
|
||||
ifsh_ms_assemble(const double *src,
|
||||
const double *totmob,
|
||||
struct ifsh_ms_data *h);
|
||||
|
||||
void
|
||||
ifsh_ms_press_flux(struct UnstructuredGrid *G, struct ifsh_ms_data *h,
|
||||
double *cpress, double *fflux);
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
#endif /* OPM_IFSH_MS_HEADER_INCLUDED */
|
||||
@ -1,10 +0,0 @@
|
||||
#include "TransTpfa.hpp"
|
||||
#include <opm/core/grid.h>
|
||||
|
||||
/* Ecplicitly initialize UnstructuredGrid versions */
|
||||
|
||||
template void tpfa_htrans_compute(const UnstructuredGrid*, const double*, double*);
|
||||
|
||||
template void tpfa_trans_compute(const UnstructuredGrid*, const double*, double*);
|
||||
|
||||
template void tpfa_eff_trans_compute(const UnstructuredGrid*, const double*, const double*, double*);
|
||||
@ -1,183 +0,0 @@
|
||||
/*===========================================================================
|
||||
//
|
||||
// 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 "config.h"
|
||||
#include <assert.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include <opm/core/pressure/tpfa/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;
|
||||
}
|
||||
@ -1,83 +0,0 @@
|
||||
/*===========================================================================
|
||||
//
|
||||
// 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 */
|
||||
@ -1,73 +0,0 @@
|
||||
/*
|
||||
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 "config.h"
|
||||
#include <stdlib.h>
|
||||
|
||||
#include <opm/core/linalg/sparse_sys.h>
|
||||
#include <opm/core/pressure/tpfa/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;
|
||||
}
|
||||
@ -91,43 +91,6 @@ struct compr_quantities_gen {
|
||||
double *voldiscr;
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
* Create management structure capable of storing derived fluid quantities
|
||||
* pertaining to a reservoir model of a particular size.
|
||||
*
|
||||
* The resources acquired using function compr_quantities_gen_allocate() must
|
||||
* be released using the destructor function compr_quantities_gen_deallocate().
|
||||
*
|
||||
* @param[in] nc Number of grid cells
|
||||
* @param[in] nf Number of grid interfaces
|
||||
* @param[in] np Number of fluid phases (and components)
|
||||
*
|
||||
* @return Fully formed management structure. Returns @c NULL in case of
|
||||
* allocation failure.
|
||||
*/
|
||||
struct compr_quantities_gen *
|
||||
compr_quantities_gen_allocate(size_t nc, size_t nf, int np);
|
||||
|
||||
|
||||
/**
|
||||
* Release resources acquired in a previous call to constructor function
|
||||
* compr_quantities_gen_allocate().
|
||||
*
|
||||
* Note that
|
||||
* <CODE>
|
||||
* compr_quantities_gen_deallocate(NULL)
|
||||
* </CODE>
|
||||
* is supported and benign (i.e., this statement behaves as
|
||||
* <CODE>free(NULL)</CODE>).
|
||||
*
|
||||
* @param[in,out] cq On input - compressible quantity management structure
|
||||
* obtained through a previous call to construction function
|
||||
* compr_quantities_gen_allocate(). On output - invalid pointer.
|
||||
*/
|
||||
void
|
||||
compr_quantities_gen_deallocate(struct compr_quantities_gen *cq);
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
@ -1,170 +0,0 @@
|
||||
/*===========================================================================
|
||||
//
|
||||
// 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 "config.h"
|
||||
#include <assert.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include <opm/core/pressure/tpfa/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;
|
||||
}
|
||||
@ -90,78 +90,6 @@ struct compr_src {
|
||||
double *saturation;
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
* Create a management structure that is, initially, capable of storing a
|
||||
* specified number of sources defined by a particular number a fluid phases.
|
||||
*
|
||||
* @param[in] np Number of fluid phases. Must be positive to actually
|
||||
* allocate any sources.
|
||||
* @param[in] nsrc Initial management capacity. If positive, attempt to
|
||||
* allocate that number of source terms. Otherwise, the
|
||||
* initial capacity is treated as (and set to) zero.
|
||||
* @return Fully formed management structure if <CODE>np > 0</CODE> and
|
||||
* allocation success. @c NULL otherwise. The resources must be released using
|
||||
* destructor function compr_src_deallocate().
|
||||
*/
|
||||
struct compr_src *
|
||||
compr_src_allocate(int np, int nsrc);
|
||||
|
||||
|
||||
/**
|
||||
* Release memory resources acquired in a previous call to constructor function
|
||||
* compr_src_allocate() and, possibly, source term insertion function
|
||||
* append_compr_source_term().
|
||||
*
|
||||
* @param[in,out] src On input - source term management structure obtained
|
||||
* through a previous call to construction function compr_src_allocate().
|
||||
* On output - invalid pointer.
|
||||
*/
|
||||
void
|
||||
compr_src_deallocate(struct compr_src *src);
|
||||
|
||||
|
||||
/**
|
||||
* Insert a new explicit source term into an existing collection.
|
||||
*
|
||||
* @param[in] c Cell influenced by this source term.
|
||||
* @param[in] np Number of fluid phases. Used for consistency checking
|
||||
* only.
|
||||
* @param[in] v Source term total reservoir volume Darcy flux. Positive
|
||||
* if the source term is to be interpreted as an injection
|
||||
* source and negative otherwise.
|
||||
* @param[in] sat Injection composition for this source term. Array of size
|
||||
* @c np. Copied to internal storage, but the actual numbers
|
||||
* are not inspected unless <CODE>v > 0.0</CODE>.
|
||||
* @param[in,out] src On input - source term management structure obtained
|
||||
* through a previous call to construction function compr_src_allocate() and,
|
||||
* possibly, another call to function append_compr_source_term(). On output -
|
||||
* source term collection that includes the new source term if successful and
|
||||
* unchanged if not.
|
||||
*
|
||||
* @return One (1, true) if successful (i.e., if the new source term could be
|
||||
* included in the existing collection) and zero (0, false) if not.
|
||||
*/
|
||||
int
|
||||
append_compr_source_term(int c ,
|
||||
int np ,
|
||||
double v ,
|
||||
const double *sat,
|
||||
struct compr_src *src);
|
||||
|
||||
|
||||
/**
|
||||
* Empty source term collection while maintaining existing capacity.
|
||||
*
|
||||
* @param[in,out] src On input - an existing source term collection with a
|
||||
* given number of sources and source capacity. On output - an empty source
|
||||
* term collection (i.e., <CODE>src->nsrc == 0</CODE>) with an unchanged
|
||||
* capacity.
|
||||
*/
|
||||
void
|
||||
clear_compr_source_term(struct compr_src *src);
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
@ -1,93 +0,0 @@
|
||||
/*
|
||||
* Copyright 2010 (c) SINTEF ICT, Applied Mathematics.
|
||||
* Jostein R. Natvig <Jostein.R.Natvig at sintef.no>
|
||||
*/
|
||||
|
||||
#include "config.h"
|
||||
#include <assert.h>
|
||||
#include <stdio.h>
|
||||
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/transport/minimal/spu_explicit.h>
|
||||
|
||||
|
||||
/* Twophase mobility-weighted upwind */
|
||||
void
|
||||
spu_explicit(struct UnstructuredGrid *g, double *s0, double *s, double *mob,
|
||||
double *dflux, double *gflux, double *src,
|
||||
double dt)
|
||||
{
|
||||
int i, f, c1, c2;
|
||||
int nc = g->number_of_cells;
|
||||
int nf = g->number_of_faces;
|
||||
|
||||
double m1, m2, flux;
|
||||
|
||||
|
||||
/* Contribution form sources */
|
||||
for (i=0; i<nc; ++i) {
|
||||
s[i] = s0[i];
|
||||
|
||||
/* Injection */
|
||||
if (src[i] > 0.0) {
|
||||
/* Assume sat==1.0 in source, and f(1.0)=1.0; */
|
||||
s[i] += dt*src[i];
|
||||
}
|
||||
|
||||
/* Production */
|
||||
else {
|
||||
m1 = mob[2*i];
|
||||
m2 = mob[2*i+1];
|
||||
s[i] += dt*src[i]* m1/(m1 + m2);
|
||||
}
|
||||
}
|
||||
|
||||
/* Contribution from internal faces */
|
||||
for (f=0; f<nf; ++f) {
|
||||
c1 = g->face_cells[2*f+0];
|
||||
c2 = g->face_cells[2*f+1];
|
||||
if ((c1 !=-1) && (c2 !=-1)) {
|
||||
|
||||
if ((dflux[f]>0.0 && gflux[f]>0.0) ||
|
||||
(dflux[f]<0.0 && gflux[f]<0.0) ) {
|
||||
/* Water mobility */
|
||||
if (dflux[f]>0) {
|
||||
m1 = mob[2*c1];
|
||||
}
|
||||
else {
|
||||
m1 = mob[2*c2];
|
||||
}
|
||||
/* Oil mobility */
|
||||
if (dflux[f] - m1*gflux[f]>0) {
|
||||
m2 = mob[2*c1+1];
|
||||
}
|
||||
else {
|
||||
m2 = mob[2*c2+1];
|
||||
}
|
||||
}
|
||||
|
||||
else {
|
||||
/* Oil mobility */
|
||||
if (dflux[f]>0) {
|
||||
m2 = mob[2*c1+1];
|
||||
}
|
||||
else {
|
||||
m2 = mob[2*c2+1];
|
||||
}
|
||||
/* Water mobility */
|
||||
if (dflux[f]+m2*gflux[f]>0) {
|
||||
m1 = mob[2*c1];
|
||||
}
|
||||
else {
|
||||
m1 = mob[2*c2];
|
||||
}
|
||||
}
|
||||
|
||||
/* Water flux */
|
||||
assert(m1+m2>0.0);
|
||||
flux = m1/(m1+m2)*(dflux[f] + m2*gflux[f]);
|
||||
s[c1] -= flux*dt;
|
||||
s[c2] += flux*dt;
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -1,19 +0,0 @@
|
||||
/*
|
||||
* Copyright 2010 (c) SINTEF ICT, Applied Mathematics.
|
||||
* Jostein R. Natvig <Jostein.R.Natvig at sintef.no>
|
||||
*/
|
||||
|
||||
#ifndef SPU_EXPLICIT_H_INCLUDED
|
||||
#define SPU_EXPLICIT_H_INCLUDED
|
||||
void
|
||||
spu_explicit(struct UnstructuredGrid *g,
|
||||
double *s0,
|
||||
double *s,
|
||||
double *mob,
|
||||
double *dflux,
|
||||
double *gflux,
|
||||
double *src,
|
||||
double dt);
|
||||
|
||||
#endif /* SPU_EXPLICIT_H_INCLUDED */
|
||||
|
||||
@ -1,372 +0,0 @@
|
||||
/*
|
||||
* Copyright 2010 (c) SINTEF ICT, Applied Mathematics.
|
||||
* Jostein R. Natvig <Jostein.R.Natvig at sintef.no>
|
||||
*/
|
||||
|
||||
|
||||
#include "config.h"
|
||||
#include <assert.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <limits.h>
|
||||
#include <math.h>
|
||||
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/transport/minimal/spu_implicit.h>
|
||||
|
||||
|
||||
|
||||
|
||||
/* Assume uniformly spaced table. */
|
||||
static double
|
||||
interpolate(int n, double h, double x0, double *tab, double x)
|
||||
{
|
||||
int i;
|
||||
double a;
|
||||
|
||||
assert(h > 0);
|
||||
assert((x-x0) < h*INT_MAX);
|
||||
assert((x-x0) > h*INT_MIN);
|
||||
|
||||
if ( x < x0 ) {
|
||||
return tab[0];
|
||||
}
|
||||
|
||||
i = ((x-x0)/h);
|
||||
|
||||
assert(i>=0);
|
||||
|
||||
if (i+1 > n-1) {
|
||||
return tab[n-1];
|
||||
}
|
||||
|
||||
a = (x-x0 - i*h) / h;
|
||||
|
||||
return (1-a) * tab[i] + a * tab[i+1];
|
||||
}
|
||||
|
||||
/* Assume uniformly spaced table. */
|
||||
static double
|
||||
differentiate(int n, double h, double x0, double *tab, double x)
|
||||
{
|
||||
int i;
|
||||
|
||||
assert(h > 0);
|
||||
assert((x-x0) < h*INT_MAX);
|
||||
assert((x-x0) > h*INT_MIN);
|
||||
assert(n>1);
|
||||
|
||||
if ( x < x0 ) {
|
||||
return (tab[1]-tab[0])/h;
|
||||
}
|
||||
|
||||
i = ((x-x0)/h);
|
||||
|
||||
assert(i>=0);
|
||||
|
||||
if (i+1 > n-1) {
|
||||
return (tab[n-1]-tab[n-2])/h;
|
||||
}
|
||||
|
||||
return (tab[i+1]-tab[i])/h;
|
||||
}
|
||||
|
||||
static void
|
||||
compute_mobilities(int n, double *s, double *mob, double *dmob, int ntab, double h, double x0, double *tab)
|
||||
{
|
||||
double *tabw=tab;
|
||||
double *tabo=tab+ntab;
|
||||
|
||||
int i;
|
||||
for (i=0; i<n; ++i) {
|
||||
*mob++ = interpolate (ntab, h, x0, tabw, *s);
|
||||
*mob++ = interpolate (ntab, h, x0, tabo, *s);
|
||||
*dmob++ = differentiate(ntab, h, x0, tabw, *s);
|
||||
*dmob++ = differentiate(ntab, h, x0, tabo, *s++);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
/*
|
||||
* To compute the derivative of the (mobility-weighted) upwind flux in
|
||||
* terms of the mobilities and the first derivatives of mobilities, we
|
||||
* use
|
||||
*
|
||||
* mw
|
||||
* F = -------
|
||||
* mo + mw
|
||||
*
|
||||
*
|
||||
* dF dt mo 1
|
||||
* --- = -- · ------- ----- (dflux + mo*gflux)
|
||||
* dmw pv mw + mo mw + mo
|
||||
*
|
||||
* dt fo
|
||||
* = -- -- (dflux + mo*gflux)
|
||||
* pv mt
|
||||
*
|
||||
*
|
||||
* dF dt fw
|
||||
* --- = -1 * -- -- (dflux - mw*gflux)
|
||||
* dmo pv mt
|
||||
*
|
||||
*
|
||||
* and the chain rule
|
||||
*
|
||||
* dF dF dmw dF dmo
|
||||
* -- = --- --- + --- ---.
|
||||
* dS dmw dS dmo dS
|
||||
*
|
||||
*
|
||||
* If the gravity term is zero, we can use the simpler rule
|
||||
*
|
||||
* dF mo * dmw - mw * dmo
|
||||
* -- = -------------------
|
||||
* dS [ mw + mo]²
|
||||
*
|
||||
*
|
||||
*
|
||||
*/
|
||||
double
|
||||
spu_implicit(struct UnstructuredGrid *g, double *s0, double *s, double h, double x0, int ntab, double *tab,
|
||||
double *dflux, double *gflux, double *src, double dt,
|
||||
void (*linear_solver)(int, int*, int*, double *, double *, double *))
|
||||
{
|
||||
int nc = g->number_of_cells;
|
||||
int nhf = g->cell_facepos[nc]; /* Too much */
|
||||
double infnorm;
|
||||
double *b;
|
||||
double *x;
|
||||
double *mob, *dmob;
|
||||
int i;
|
||||
int it;
|
||||
|
||||
/* Allocate space for linear system etc.*/
|
||||
sparse_t *S = malloc(sizeof *S);
|
||||
if (S) {
|
||||
S->ia = malloc((nc+1) * sizeof *S->ia);
|
||||
S->ja = malloc( nhf * sizeof *S->ja);
|
||||
S->sa = malloc( nhf * sizeof *S->sa);
|
||||
S->n = nc;
|
||||
S->m = nc;
|
||||
}
|
||||
else {
|
||||
assert(0);
|
||||
}
|
||||
|
||||
b = malloc(nc * sizeof *b);
|
||||
x = malloc(nc * sizeof *x);
|
||||
|
||||
mob = malloc(g->number_of_cells *2* sizeof *mob);
|
||||
dmob = malloc(g->number_of_cells *2* sizeof *dmob);
|
||||
|
||||
infnorm = 1.0;
|
||||
it = 0;
|
||||
while (infnorm > 1e-9 && it++ < 20) {
|
||||
compute_mobilities(g->number_of_cells, s, mob, dmob, ntab, h, x0, tab);
|
||||
spu_implicit_assemble(g, s0, s, mob, dmob, dflux, gflux, src, dt, S, b);
|
||||
|
||||
/* Compute inf-norm of residual */
|
||||
infnorm = 0.0;
|
||||
for(i=0; i<nc; ++i) {
|
||||
infnorm = (infnorm > fabs(b[i]) ? infnorm : fabs(b[i]));
|
||||
}
|
||||
fprintf(stderr, " Max norm of residual: %e\n", infnorm);
|
||||
|
||||
/* Solve system */
|
||||
(*linear_solver)(S->m, S->ia, S->ja, S->sa, b, x);
|
||||
|
||||
/* Return x. */
|
||||
for(i=0; i<nc; ++i) {
|
||||
s[i] = s[i] - x[i];
|
||||
s[i] = s[i]>1.0 ? 1.0 : (s[i] < 0.0 ? 0.0 : s[i]);
|
||||
}
|
||||
}
|
||||
|
||||
free(mob);
|
||||
free(dmob);
|
||||
|
||||
free(b);
|
||||
free(x);
|
||||
free(S->ia);
|
||||
free(S->ja);
|
||||
free(S->sa);
|
||||
free(S);
|
||||
|
||||
return infnorm;
|
||||
}
|
||||
|
||||
|
||||
static void
|
||||
phase_upwind_mobility(double darcyflux, double gravityflux, int i, int c,
|
||||
double *mob, double *dmob, double *m, double *dm, int *cix)
|
||||
{
|
||||
/* ====================================================== */
|
||||
/* If darcy flux and gravity flux have same sign... */
|
||||
if ( (darcyflux>0.0 && gravityflux>0.0) ||
|
||||
(darcyflux<0.0 && gravityflux<0.0) ) {
|
||||
|
||||
/* Positive water phase flux */
|
||||
if (darcyflux>0) {
|
||||
m [0] = mob [2*i+0];
|
||||
dm[0] = dmob[2*i+0];
|
||||
cix[0] = i;
|
||||
}
|
||||
|
||||
/* Negative water phase flux */
|
||||
else {
|
||||
m [0] = mob [2*c+0];
|
||||
dm[0] = dmob[2*c+0];
|
||||
cix[0] = c;
|
||||
}
|
||||
|
||||
/* Positive oil phase flux */
|
||||
if (darcyflux - m[0]*gravityflux>0) {
|
||||
m [1] = mob[2*i+1];
|
||||
dm[1] = dmob[2*i+1];
|
||||
cix[1] = i;
|
||||
}
|
||||
|
||||
/* Negative oil phase flux */
|
||||
else {
|
||||
m [1] = mob[2*c+1];
|
||||
dm[1] = dmob[2*c+1];
|
||||
cix[1] = c;
|
||||
}
|
||||
}
|
||||
|
||||
/* ====================================================== */
|
||||
/* If Darcy flux and gravity flux have opposite sign... */
|
||||
else {
|
||||
|
||||
/* Positive oil phase flux */
|
||||
if (darcyflux>0) {
|
||||
m [1] = mob[2*i+1];
|
||||
dm[1] = dmob[2*i+1];
|
||||
cix[1] = i;
|
||||
}
|
||||
|
||||
/* Negative oil phase flux */
|
||||
else {
|
||||
m [1] = mob[2*c+1];
|
||||
dm[1] = dmob[2*c+1];
|
||||
cix[1] = c;
|
||||
}
|
||||
|
||||
/* Positive water phase flux */
|
||||
if (darcyflux+m[1]*gravityflux>0) {
|
||||
m [0] = mob[2*i+0];
|
||||
dm[0] = dmob[2*i+0];
|
||||
cix[0] = i;
|
||||
}
|
||||
|
||||
/* Negative water phase flux */
|
||||
else {
|
||||
m [0] = mob[2*c+0];
|
||||
dm[0] = dmob[2*c+0];
|
||||
cix[0] = c;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
spu_implicit_assemble(struct UnstructuredGrid *g, double *s0, double *s, double *mob, double *dmob,
|
||||
double *dflux, double *gflux, double *src, double dt, sparse_t *S,
|
||||
double *b)
|
||||
{
|
||||
int i, k, f, c1, c2, c;
|
||||
int nc = g->number_of_cells;
|
||||
|
||||
double m[6] = { 0, 0, 0, 0, 0, 0 };
|
||||
double dm[2] = { 0, 0 };
|
||||
int cix[2] = { 0, 0 };
|
||||
double m1, m2, dm1, dm2, mt2;
|
||||
|
||||
int *pja;
|
||||
double *psa;
|
||||
double *d;
|
||||
double sgn;
|
||||
double df, gf;
|
||||
|
||||
pja = S->ja;
|
||||
psa = S->sa;
|
||||
|
||||
/* Assemble system */
|
||||
S->ia[0] = 0;
|
||||
for (i=0; i<nc; ++i) {
|
||||
|
||||
/* Store position of diagonal element */
|
||||
d = psa++;
|
||||
*pja++ = i;
|
||||
*d = 0.0;
|
||||
|
||||
/* Accumulation term */
|
||||
b[i] = (s[i] - s0[i]);
|
||||
*d += 1.0;
|
||||
|
||||
/* Flux terms follows*/
|
||||
for (k=g->cell_facepos[i]; k<g->cell_facepos[i+1]; ++k) {
|
||||
f = g->cell_faces[k];
|
||||
c1 = g->face_cells[2*f+0];
|
||||
c2 = g->face_cells[2*f+1];
|
||||
|
||||
/* Skip all boundary terms (for now). */
|
||||
if (c1 == -1 || c2 == -1) { continue; }
|
||||
|
||||
/* Set cell index of other cell, set correct sign of fluxes. */
|
||||
c = (i==c1 ? c2 : c1);
|
||||
sgn = (i==c1)*2.0 - 1.0;
|
||||
df = sgn * dt*dflux[f];
|
||||
gf = sgn * dt*gflux[f];
|
||||
|
||||
phase_upwind_mobility(df, gf, i, c, mob, dmob, m, dm, cix);
|
||||
|
||||
/* Ensure we do not divide by zero. */
|
||||
if (m[0] + m[1] > 0.0) {
|
||||
|
||||
b[i] += m[0]/(m[0]+m[1])*(df + m[1]*gf);
|
||||
|
||||
mt2 = (m[0]+m[1])*(m[0]+m[1]);
|
||||
*psa = 0.0;
|
||||
*pja = c;
|
||||
|
||||
/* dFw/dmw·dmw/dsw */
|
||||
if (cix[0] == c ) { *psa += m[1]/mt2*(df + m[1]*gf)*dm[0]; }
|
||||
else { *d += m[1]/mt2*(df + m[1]*gf)*dm[0]; }
|
||||
|
||||
/* dFw/dmo·dmo/dsw */
|
||||
if (cix[1] == c) { *psa += -m[0]/mt2*(df - m[0]*gf)*dm[1]; }
|
||||
else { *d += -m[0]/mt2*(df - m[0]*gf)*dm[1]; }
|
||||
|
||||
|
||||
if (cix[0] == c || cix[1] == c) {
|
||||
++psa;
|
||||
++pja;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* Injection */
|
||||
if (src[i] > 0.0) {
|
||||
/* Assume sat==1.0 in source, and f(1.0)=1.0; */
|
||||
m1 = 1.0;
|
||||
m2 = 0.0;
|
||||
b[i] -= dt*src[i] * m1/(m1+m2);
|
||||
}
|
||||
|
||||
/* Production */
|
||||
else {
|
||||
m1 = mob [2*i+0];
|
||||
m2 = mob [2*i+1];
|
||||
dm1 = dmob[2*i+0];
|
||||
dm2 = dmob[2*i+1];
|
||||
|
||||
*d -= dt*src[i] *(m2*dm1-m1*dm2)/(m1+m2)/(m1+m2);
|
||||
b[i] -= dt*src[i] *m1/(m1+m2);
|
||||
|
||||
}
|
||||
S->ia[i+1] = pja - S->ja;
|
||||
}
|
||||
}
|
||||
@ -1,30 +0,0 @@
|
||||
/*
|
||||
* Copyright 2010 (c) SINTEF ICT, Applied Mathematics.
|
||||
* Jostein R. Natvig <Jostein.R.Natvig at sintef.no>
|
||||
*/
|
||||
|
||||
#ifndef SPU_IMPLICIT_H_INCLUDED
|
||||
#define SPU_IMPLICIT_H_INCLUDED
|
||||
|
||||
typedef struct Sparse {
|
||||
int m;
|
||||
int n;
|
||||
int *ia;
|
||||
int *ja;
|
||||
double *sa;
|
||||
} sparse_t;
|
||||
|
||||
|
||||
|
||||
void
|
||||
spu_implicit_assemble(struct UnstructuredGrid *g, double *s0, double *s, double *mob, double *dmob,
|
||||
double *dflux, double *gflux, double *src, double dt, sparse_t *S,
|
||||
double *b);
|
||||
|
||||
double
|
||||
spu_implicit(struct UnstructuredGrid *g, double *s0, double *s, double h, double x0, int ntab, double *tab,
|
||||
double *dflux, double *gflux, double *src, double dt,
|
||||
void (*linear_solver)(int, int*, int*, double *, double *, double *));
|
||||
|
||||
#endif /* SPU_IMPLICIT_H_INCLUDED */
|
||||
|
||||
Loading…
Reference in New Issue
Block a user