diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 569d4958..bd6ecda5 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -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 diff --git a/opm/core/pressure/cfsh.c b/opm/core/pressure/cfsh.c deleted file mode 100644 index 8ca7e9f7..00000000 --- a/opm/core/pressure/cfsh.c +++ /dev/null @@ -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 . -*/ - -#include "config.h" -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - - - - -/* ---------------------------------------------------------------------- */ -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 */ - } -} diff --git a/opm/core/pressure/fsh.c b/opm/core/pressure/fsh.c deleted file mode 100644 index 0dbd3fcd..00000000 --- a/opm/core/pressure/fsh.c +++ /dev/null @@ -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 . -*/ - -#include "config.h" -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - - - - -/* ---------------------------------------------------------------------- */ -/* 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; - } -} diff --git a/opm/core/pressure/fsh.h b/opm/core/pressure/fsh.h deleted file mode 100644 index 08cfebc5..00000000 --- a/opm/core/pressure/fsh.h +++ /dev/null @@ -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 . -*/ - -#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 -#include -#include - -#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 FSH object. - * - * @param[in,out] h FSH 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 - * h->A and h->b. - * - * @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 - * h->A and h->b. - * - * @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, h->x. - * 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 */ diff --git a/opm/core/pressure/fsh_common_impl.c b/opm/core/pressure/fsh_common_impl.c deleted file mode 100644 index 463d5455..00000000 --- a/opm/core/pressure/fsh_common_impl.c +++ /dev/null @@ -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 . -*/ - -#include "config.h" -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -#include -#include - -#include -#include - -#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 ]; - } -} diff --git a/opm/core/pressure/fsh_common_impl.h b/opm/core/pressure/fsh_common_impl.h deleted file mode 100644 index 3147f4f9..00000000 --- a/opm/core/pressure/fsh_common_impl.h +++ /dev/null @@ -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 . -*/ - -#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 */ diff --git a/opm/core/pressure/ifsh.c b/opm/core/pressure/ifsh.c deleted file mode 100644 index a174bda2..00000000 --- a/opm/core/pressure/ifsh.c +++ /dev/null @@ -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 . -*/ - -#include "config.h" -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - - -/* ---------------------------------------------------------------------- */ -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 */ - } -} diff --git a/opm/core/pressure/legacy_well.c b/opm/core/pressure/legacy_well.c deleted file mode 100644 index 94dae9d8..00000000 --- a/opm/core/pressure/legacy_well.c +++ /dev/null @@ -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 . -*/ - -#include "config.h" -#include - -#include - - -/* 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; - } - } -} diff --git a/opm/core/pressure/legacy_well.h b/opm/core/pressure/legacy_well.h index ccf93089..39f40201 100644 --- a/opm/core/pressure/legacy_well.h +++ b/opm/core/pressure/legacy_well.h @@ -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 - * nc + 1 if successful. - * @param[out] cwells Cell-to-well mapping. Points to array - * of size W->well_connpos[ - * W->number_of_wells] if successful. - * @return Positive number (size of *cwells) - * 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 W->well_cells). - * - * @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 diff --git a/opm/core/pressure/mimetic/hybsys.c b/opm/core/pressure/mimetic/hybsys.c deleted file mode 100644 index 694c5616..00000000 --- a/opm/core/pressure/mimetic/hybsys.c +++ /dev/null @@ -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 . -*/ - -#include "config.h" -#include -#include -#include -#include - -#include -#include - - -#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); -} diff --git a/opm/core/pressure/mimetic/hybsys.h b/opm/core/pressure/mimetic/hybsys.h deleted file mode 100644 index 12ace45e..00000000 --- a/opm/core/pressure/mimetic/hybsys.h +++ /dev/null @@ -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 . -*/ - -#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. hybsys_free(NULL) 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. hybsys_well_free(NULL) 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 sys->L and wsys->F1. - * - * @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 cwpos[nc]. 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 sys->F1) 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 pwconn[nc]. 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 - * pwconn[nc]. - * @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 - * pconn[nc]. - * @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 nw > 0, - * 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 pwconn[nc]. 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 - * pwconn[nw]. - * @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 */ diff --git a/opm/core/pressure/mimetic/hybsys_global.c b/opm/core/pressure/mimetic/hybsys_global.c deleted file mode 100644 index 33a352e0..00000000 --- a/opm/core/pressure/mimetic/hybsys_global.c +++ /dev/null @@ -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 . -*/ - -#include "config.h" -#include -#include -#include - -/* ---> Layering violation! */ -#include -/* <--- Layering violation! */ - -#include - - -#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]; - } -} diff --git a/opm/core/pressure/mimetic/hybsys_global.h b/opm/core/pressure/mimetic/hybsys_global.h deleted file mode 100644 index 1d98cd0d..00000000 --- a/opm/core/pressure/mimetic/hybsys_global.h +++ /dev/null @@ -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 . -*/ - -#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 -#include -#include - - -/** - * 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 (G->number_of_faces) - * plus the number of wells (W->number_of_wells). - */ -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 - * G->number_of_faces - * 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 */ diff --git a/opm/core/pressure/msmfem/coarse_conn.c b/opm/core/pressure/msmfem/coarse_conn.c deleted file mode 100644 index 065bece5..00000000 --- a/opm/core/pressure/msmfem/coarse_conn.c +++ /dev/null @@ -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 . -*/ - -#include "config.h" -#include -#include -#include -#include -#include - -#include -#include - -#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); -} diff --git a/opm/core/pressure/msmfem/coarse_conn.h b/opm/core/pressure/msmfem/coarse_conn.h deleted file mode 100644 index a35b9f65..00000000 --- a/opm/core/pressure/msmfem/coarse_conn.h +++ /dev/null @@ -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 . -*/ - -#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 */ diff --git a/opm/core/pressure/msmfem/coarse_sys.c b/opm/core/pressure/msmfem/coarse_sys.c deleted file mode 100644 index ec746e19..00000000 --- a/opm/core/pressure/msmfem/coarse_sys.c +++ /dev/null @@ -1,1865 +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 . -*/ - -#include "config.h" -#include -#include -#include -#include -#include -#include - - -#ifndef DEBUG_OUTPUT -#define DEBUG_OUTPUT 0 -#endif - -#if DEBUG_OUTPUT -#include -#endif - - -#include -#include - -#include -#include -#include - -#include -#include -#include - - -#if defined(MAX) -#undef MAX -#endif -#define MAX(a,b) (((a) > (b)) ? (a) : (b)) - - -/* ====================================================================== - * Data structures - * ====================================================================== */ - -struct coarse_sys_meta { - size_t max_ngconn; /* max(diff(face_pos)) */ - size_t sum_ngconn2; /* sum(diff(face_pos)^2) */ - - size_t max_blk_cells; /* max(accumarray(p,1)) */ - size_t max_blk_nhf; /* max(accumarray(p,diff(face_pos))) */ - size_t max_blk_nfsf; /* Maximum # block fine-scale faces */ - size_t max_blk_sum_nhf2; /* max_i \sum_{j\in\Omega_i} ncf(j)^2 */ - size_t max_cf_nf; /* Maximum # fs faces in a coarse face */ - size_t n_act_bf; /* Number of active coarse connections */ - - int *blk_nhf; /* Number of fs hfaces per block */ - int *blk_nfsf; /* Number of fs faces per block */ - - int *loc_fno; /* Local (fs) face numbering */ - int *ncf; /* diff(face_pos) */ - int *pconn2; /* cumsum([0; diff(face_pos).^2]) */ - - int *pb2c, *b2c; /* Block->cell mapping (CSR packed) */ - - int *bfno; /* Active basis function numbering */ - int *loc_dofno; /* Block-local DOF number of a CF */ - - /* -------------------------------------------------------------- */ - int *data; /* Linear storage. Don't touch. */ -}; - - -struct bf_asm_data { - struct hybsys *fsys; /* Fine-scale hybrid system contributions */ - - struct CSRMatrix *A; /* BF coefficient matrix */ - double *b; /* BF system RHS */ - double *x; /* BF system solution */ - - double *v; /* BF (hf) flux. */ - double *p; /* BF pressure. */ - double *flux; /* BF flux. Symmetrised. */ - - double *gpress; /* BF gravity contrib. (== 0) */ - - double *work; /* Back-substitution work array */ - - int *pdof; /* Indirection pointer to linearised DOF */ - int *dof; /* Linearised DOFs per BF */ - int *fcount; /* Flux symmetrisation face count. */ - - /* -------------------------------------------------------------- */ - int *idata; /* Linear (integer) storage */ - double *ddata; /* Linear (double) storage */ -}; - - -/* ====================================================================== - * Memory management - * ====================================================================== */ - - -/* ---------------------------------------------------------------------- */ -static void -coarse_sys_meta_destroy(struct coarse_sys_meta *m) -/* ---------------------------------------------------------------------- */ -{ - if (m != NULL) { - free(m->data); - } - - free(m); -} - - -/* ---------------------------------------------------------------------- */ -static struct coarse_sys_meta * -coarse_sys_meta_allocate(size_t nblocks, size_t nfaces_c, - size_t nc , size_t nfaces_f) -/* ---------------------------------------------------------------------- */ -{ - size_t i, alloc_sz; - struct coarse_sys_meta *new; - - new = malloc(1 * sizeof *new); - - if (new != NULL) { - alloc_sz = nblocks; /* blk_nhf */ - alloc_sz += nblocks; /* blk_nfsf */ - alloc_sz += nc; /* ncf */ - alloc_sz += nc + 1; /* pconn2 */ - alloc_sz += nfaces_f; /* loc_fno */ - alloc_sz += nblocks + 1; /* pb2c */ - alloc_sz += nc; /* b2c */ - alloc_sz += nfaces_c; /* bfno */ - alloc_sz += 2*nfaces_c; /* loc_dofno */ - - new->data = malloc(alloc_sz * sizeof *new->data); - - if (new->data == NULL) { - coarse_sys_meta_destroy(new); - - new = NULL; - } else { - for (i = 0; i < alloc_sz; i++) { - new->data[i] = 0; - } - - new->blk_nhf = new->data; - new->blk_nfsf = new->blk_nhf + nblocks; - new->ncf = new->blk_nfsf + nblocks; - new->pconn2 = new->ncf + nc; - new->loc_fno = new->pconn2 + nc + 1; - - new->pb2c = new->loc_fno + nfaces_f; - new->b2c = new->pb2c + nblocks + 1; - - new->bfno = new->b2c + nc; - new->loc_dofno = new->bfno + nfaces_c; - } - } - - return new; -} - - -/* ---------------------------------------------------------------------- */ -static void -bf_asm_data_deallocate(struct bf_asm_data *data) -/* ---------------------------------------------------------------------- */ -{ - if (data != NULL) { - free (data->ddata); - free (data->idata); - csrmatrix_delete(data->A); - hybsys_free (data->fsys); - } - - free(data); -} - - -/* ---------------------------------------------------------------------- */ -static struct bf_asm_data * -bf_asm_data_allocate(struct UnstructuredGrid *g, - struct coarse_sys_meta *m) -/* ---------------------------------------------------------------------- */ -{ - int tot_ngconn; - size_t max_nhf, max_cells, max_faces, nnz; - size_t alloc_sz; - struct bf_asm_data *new; - - new = malloc(1 * sizeof *new); - - if (new != NULL) { - max_nhf = 2 * m->max_blk_nhf; - max_cells = 2 * m->max_blk_cells; - max_faces = 2 * m->max_blk_nfsf; - nnz = 2 * m->max_blk_sum_nhf2; - - tot_ngconn = g->cell_facepos[ g->number_of_cells ]; - - new->fsys = hybsys_allocate_symm((int) m->max_ngconn, - g->number_of_cells, - tot_ngconn); - - new->A = csrmatrix_new_known_nnz(max_faces, nnz); - - alloc_sz = max_cells + 1; /* pdof */ - alloc_sz += max_nhf; /* dof */ - alloc_sz += max_faces; /* fcount */ - - new->idata = malloc(alloc_sz * sizeof *new->idata); - - alloc_sz = 2 * max_faces; /* b, x */ - alloc_sz += 1 * max_nhf; /* v */ - alloc_sz += 1 * max_cells; /* p */ - alloc_sz += 1 * max_faces; /* flux */ - alloc_sz += tot_ngconn; /* gpress */ - alloc_sz += m->max_ngconn; /* work */ - - new->ddata = malloc(alloc_sz * sizeof *new->ddata); - - if ((new->fsys == NULL) || (new->A == NULL) || - (new->idata == NULL) || (new->ddata == NULL)) { - bf_asm_data_deallocate(new); - new = NULL; - } else { - new->pdof = new->idata; - new->dof = new->pdof + max_cells + 1; - new->fcount = new->dof + max_nhf; - - new->b = new->ddata; - new->x = new->b + max_faces; - new->v = new->x + max_faces; - new->p = new->v + max_nhf; - - new->flux = new->p + max_cells; - - new->gpress = new->flux + max_faces; - new->work = new->gpress + tot_ngconn; - - hybsys_init((int) m->max_ngconn, new->fsys); - } - } - - return new; -} - - -/* ====================================================================== - * Utilities. - * ====================================================================== */ - - -/* max(diff(p(1:n))) */ -/* ---------------------------------------------------------------------- */ -static int -max_diff(int n, int *p) -/* ---------------------------------------------------------------------- */ -{ - int i, d, ret; - - assert ((n > 0) && (p != NULL)); - - ret = p[1] - p[0]; assert (ret >= 0); - for (i = 1; i < n; i++) { - d = p[i + 1] - p[i]; - ret = (d > ret) ? d : ret; - } - - return ret; -} - - -/* Enumerate active basis functions/coarse connetions according to - * block proximity. - * - * Returns number of active connections. */ -/* ---------------------------------------------------------------------- */ -static int -enumerate_active_bf(struct coarse_topology *ct, - struct coarse_sys_meta *m) -/* ---------------------------------------------------------------------- */ -{ - int act, cf, b_in, b_out, b1, b2, p; - - for (p = 0; p < ct->nfaces; p++) { - m->bfno[p] = -1; - } - - act = 0; - - for (b_in = p = 0; b_in < ct->nblocks; b_in++) { - for (; p < ct->blkfacepos[b_in + 1]; p++) { - cf = ct->blkfaces[p]; - - if (m->bfno[cf] < 0) { /* Previously undiscovered */ - b1 = ct->neighbours[2*cf + 0]; - b2 = ct->neighbours[2*cf + 1]; - - assert (b1 != b2); - - b_out = (b1 == b_in) ? b2 : b1; - - if (b_out >= 0) { - /* Restricted to internal cf's for now. */ - m->bfno[cf] = act++; - } - } - } - } - - return act; -} - - -/* ---------------------------------------------------------------------- */ -static void -compute_loc_dofno(struct coarse_topology *ct, - struct coarse_sys_meta *m) -/* ---------------------------------------------------------------------- */ -{ - int b, nb, p, cf, locno, col_off; - - for (p = 0; p < 2 * ct->nfaces; p++) { - m->loc_dofno[p] = -1; - } - - nb = ct->nblocks; - - for (b = p = 0; b < nb; b++) { - locno = 0; - - for (; p < ct->blkfacepos[b + 1]; p++) { - cf = ct->blkfaces[p]; - - if (m->bfno[cf] >= 0) { - col_off = 1 - (ct->neighbours[2*cf + 0] == b); - - assert (m->loc_dofno[2*cf + col_off] == -1); - - m->loc_dofno[2*cf + col_off] = locno++; - } - } - } -} - - -/* ---------------------------------------------------------------------- */ -static void -coarse_sys_meta_fill(int nc, const int *pgconn, - size_t nneigh, const int *neigh, - const int *p, - struct coarse_topology *ct, - struct coarse_sys_meta *m) -/* ---------------------------------------------------------------------- */ -{ - int c1, b1, c2, b2, i, n; - size_t f, blk_sum_nhf2; - - m->max_blk_nhf = m->max_blk_nfsf = 0; - - for (f = 0; f < nneigh; f++) { - c1 = neigh[2*f + 0]; b1 = (c1 >= 0) ? p[c1] : -1; - c2 = neigh[2*f + 1]; b2 = (c2 >= 0) ? p[c2] : -1; - - assert ((b1 >= 0) || (b2 >= 0)); - - if (b1 >= 0) { - m->blk_nhf[b1] += 1; - m->max_blk_nhf = MAX(m->max_blk_nhf, - (size_t) m->blk_nhf[b1]); - } - - if (b2 >= 0) { - m->blk_nhf[b2] += 1; - m->max_blk_nhf = MAX(m->max_blk_nhf, - (size_t) m->blk_nhf[b2]); - } - - if (b1 >= 0) { - m->blk_nfsf[b1] += 1; - m->max_blk_nfsf = MAX(m->max_blk_nfsf, - (size_t) m->blk_nfsf[b1]); - - if ((b2 >= 0) && (b2 != b1)) { - m->blk_nfsf[b2] += 1; - m->max_blk_nfsf = MAX(m->max_blk_nfsf, - (size_t) m->blk_nfsf[b2]); - } - } else { - m->blk_nfsf[b2] += 1; - m->max_blk_nfsf = MAX(m->max_blk_nfsf, - (size_t) m->blk_nfsf[b2]); - } - } - - for (f = 0; f < nneigh; f++) { - m->loc_fno[f] = -1; - } - - m->max_cf_nf = 0; - - for (f = 0; f < (size_t) ct->nfaces; f++) { - m->max_cf_nf = MAX(m->max_cf_nf, - (size_t) (ct->subfacepos[f + 1] - - ct->subfacepos[f])); - } - - m->max_ngconn = m->sum_ngconn2 = 0; - for (c1 = 0; c1 < nc; c1++) { - n = pgconn[c1 + 1] - pgconn[c1]; - - m->max_ngconn = MAX(m->max_ngconn, (size_t) n); - m->sum_ngconn2 += n * n; - - m->ncf[c1] = n; - m->pconn2[c1 + 1] = m->pconn2[c1] + (n * n); - } - - partition_invert(nc, p, m->pb2c, m->b2c); - - m->max_blk_cells = 0; - for (b1 = 0; b1 < ct->nblocks; b1++) { - m->max_blk_cells = MAX(m->max_blk_cells, - (size_t) (m->pb2c[b1 + 1] - - m->pb2c[b1 + 0])); - } - - m->max_blk_sum_nhf2 = 0; - for (b1 = 0; b1 < ct->nblocks; b1++) { - blk_sum_nhf2 = 0; - - for (i = m->pb2c[b1]; i < m->pb2c[b1 + 1]; i++) { - c1 = m->b2c[i]; - blk_sum_nhf2 += m->pconn2[c1 + 1] - m->pconn2[c1]; - } - - m->max_blk_sum_nhf2 = MAX(m->max_blk_sum_nhf2, blk_sum_nhf2); - } - - m->n_act_bf = enumerate_active_bf(ct, m); - - compute_loc_dofno(ct, m); -} - - -/* ---------------------------------------------------------------------- */ -static struct coarse_sys_meta * -coarse_sys_meta_construct(struct UnstructuredGrid *g, const int *p, - struct coarse_topology *ct) -/* ---------------------------------------------------------------------- */ -{ - struct coarse_sys_meta *m; - - m = coarse_sys_meta_allocate(ct->nblocks, ct->nfaces, - g->number_of_cells, g->number_of_faces); - - if (m != NULL) { - coarse_sys_meta_fill(g->number_of_cells, - g->cell_facepos, - g->number_of_faces, - g->face_cells, p, ct, m); - } - - return m; -} - -#define USE_MIM_IP_SIMPLE 0 -#define USE_MIM_IP_TPFA 1 -#define USE_MIM_IP_QFAMILY 0 - -#if USE_MIM_IP_SIMPLE -/* ---------------------------------------------------------------------- */ -static double * -compute_fs_ip(struct UnstructuredGrid *g, const double *perm, - const struct coarse_sys_meta *m) -/* ---------------------------------------------------------------------- */ -{ - double *Binv; - - Binv = malloc(m->sum_ngconn2 * sizeof *Binv); - - if (Binv != NULL) { - mim_ip_simple_all(g->number_of_cells, g->dimensions, m->max_ngconn, - g->cell_facepos, g->cell_faces, - g->face_cells, g->face_centroids, g->face_normals, - g->face_areas, g->cell_centroids, g->cell_volumes, - (double *) perm, Binv); /* const_cast<>() */ - } - - return Binv; -} -#elif USE_MIM_IP_TPFA -#include -/* ---------------------------------------------------------------------- */ -static double * -compute_fs_ip(struct UnstructuredGrid *g, const double *perm, - const struct coarse_sys_meta *m) -/* ---------------------------------------------------------------------- */ -{ - size_t c, nc, nconn, p1, p2, i, j; - double *Binv, *htrans; - - nc = g->number_of_cells; - - Binv = malloc(m->sum_ngconn2 * sizeof *Binv); - htrans = malloc(g->cell_facepos[ nc ] * sizeof *htrans); - - if ((Binv != NULL) && (htrans != NULL)) { - tpfa_htrans_compute(g, perm, htrans); - - for (c = p2 = 0; c < nc; c++) { - p1 = g->cell_facepos[c + 0] ; - nconn = g->cell_facepos[c + 1] - p1; - - for (i = 0; i < nconn; i++) { - Binv[p2 + i*(nconn + 1)] = htrans[p1 + i]; - - for (j = (i + 1) % nconn; j != i; j = (j + 1) % nconn) { - Binv[p2 + i*nconn + j] = 0.0; - } - } - - p2 += nconn * nconn; - } - } - - free(htrans); - - return Binv; -} -#endif - - -/* Create basis function weighting source term (unsigned) based on - * trace of permeability. - * - * Returns valid pointer (one scalar per cell in underlying fine-scale - * grid model) if succesful, and NULL otherwise. */ -/* ---------------------------------------------------------------------- */ -static double * -perm_weighting(size_t nc, size_t nd, - const double *perm, - const double *cvol) -/* ---------------------------------------------------------------------- */ -{ - size_t c, d, off; - double t; - double *w; - - w = malloc(nc * sizeof *w); - - if (w != NULL) { - for (c = off = 0; c < nc; c++, off += nd*nd) { - t = 0.0; - - for (d = 0; d < nd; d++) { - t += perm[off + d*(nd + 1)]; - } - - w[c] = t * cvol[c]; - } - } - - return w; -} - - -/* Use prescribed sources if applicable (replace synthetic source term). - * - * Returns one (1) if successful (needs to allocate internal data) and - * zero (0) otherwise. */ -/* ---------------------------------------------------------------------- */ -static int -enforce_explicit_source(size_t nc, size_t nb, const int *p, - const double *src, struct coarse_sys_meta *m, - double *w) -/* ---------------------------------------------------------------------- */ -{ - int i, ret; - size_t c, b; - int *has_src; - - ret = 0; - has_src = malloc(nb * sizeof *has_src); - - if (has_src != NULL) { - for (b = 0; b < nb; b++) { has_src[b] = 0; } - - for (c = 0; c < nc; c++) { - has_src[p[c]] += fabs(src[c]) > 0.0; - } - - for (b = 0; b < nb; b++) { - if (has_src[b]) { - for (i = m->pb2c[b]; i < m->pb2c[b + 1]; i++) { - w[m->b2c[i]] = 0.0; - } - } - } - - for (c = 0; c < nc; c++) { - if (fabs(src[c]) > 0.0) { - w[c] = src[c]; - } - } - - ret = 1; - } - - free(has_src); - - return ret; -} - - -/* Enforce \int_{\Omega_i} w(x) dx == 1 for all blocks, \Omega_i. - * - * Allocates internal work array. - * - * Returns one (1) if successful, and zero (0) otherwise. */ -/* ---------------------------------------------------------------------- */ -static int -normalize_weighting(size_t nc, size_t nb, const int *p, double *w) -/* ---------------------------------------------------------------------- */ -{ - int ret; - size_t c, b; - double *bw; - - ret = 0; - bw = malloc(nb * sizeof *bw); - - if (bw != NULL) { - for (b = 0; b < nb; b++) { bw[ b ] = 0.0; } - for (c = 0; c < nc; c++) { bw[p[c]] += w[c]; } - - for (c = 0; c < nc; c++) { - assert (fabs(bw[p[c]]) > 0.0); - - w[c] /= bw[p[c]]; - } - - ret = 1; - } - - free(bw); - - return ret; -} - - -/* Create basis function weighting term (unsigned), one scalar per - * grid cell in the underlying fine-scale grid. Integral one per - * block. Satisfies any prescribed external source terms. - * - * Returns valid ponter if successful and NULL if not. */ -/* ---------------------------------------------------------------------- */ -static double * -coarse_weight(struct UnstructuredGrid *g, size_t nb, - const int *p, - struct coarse_sys_meta *m, - const double *perm, const double *src) -/* ---------------------------------------------------------------------- */ -{ - int ok; - double *w; - - ok = 0; - w = perm_weighting(g->number_of_cells, - g->dimensions, perm, g->cell_volumes); - - if (w != NULL) { - ok = enforce_explicit_source(g->number_of_cells, - nb, p, src, m, w); - } - - if (ok) { - ok = normalize_weighting(g->number_of_cells, nb, p, w); - } - - if (!ok) { free(w); w = NULL; } - - return w; -} - - -/* ---------------------------------------------------------------------- */ -/* Determine which, and how many, degrees of freedom (i.e., active - * basis functions) are connected to every coarse block. Allocates - * and fills the CSR array pair (sys->blkdof_pos,sys->blkdof). - * - * Returns 1, and sets ->blkdof_pos,->blkdof to valid arrays if successful. - * Returns 0, and sets ->blkdof_pos,->blkdof to NULL if not. - * - * Uses the standard two-pass CSR push-back building strategy. */ -/* ---------------------------------------------------------------------- */ -static int -blkdof_fill(struct coarse_topology *ct, - struct coarse_sys_meta *m, - struct coarse_sys *sys) -/* ---------------------------------------------------------------------- */ -{ - int p, ret = 0, dof; - size_t b, nb; - - nb = ct->nblocks; - - assert(nb > 0); - - sys->blkdof_pos = malloc((nb + 1) * sizeof *sys->blkdof_pos); - - if (sys->blkdof_pos != NULL) { - for (b = 0; b <= nb; b++) { sys->blkdof_pos[b] = 0; } - - /* Count number of active BFs per block */ - for (b = 0, p = 0; b < nb; b++) { - for (; p < ct->blkfacepos[b + 1]; p++) { - dof = m->bfno[ ct->blkfaces[p] ]; - - sys->blkdof_pos[b + 1] += dof >= 0; - } - } - - /* Derive CSR start pointers */ - for (b = 1; b <= nb; b++) { - sys->blkdof_pos[0] += sys->blkdof_pos[b]; - sys->blkdof_pos[b] = sys->blkdof_pos[0] - sys->blkdof_pos[b]; - } - - sys->blkdof = malloc(sys->blkdof_pos[0] * sizeof *sys->blkdof); - - /* Fill each block's active BFs if we can allocate ->blkdof */ - if (sys->blkdof == NULL) { - free(sys->blkdof_pos); - sys->blkdof_pos = NULL; - } else { - sys->blkdof_pos[0] = 0; - - for (b = 0, p = 0; b < nb; b++) { - for (; p < ct->blkfacepos[b + 1]; p++) { - dof = m->bfno[ ct->blkfaces[p] ]; - - if (dof >= 0) { - sys->blkdof[ sys->blkdof_pos[b + 1] ++ ] = dof; - } - } - } - - ret = 1; - } - } - - return ret; -} - - -/* ---------------------------------------------------------------------- */ -/* Compute aggregate allocation sizes for ->basis, ->cell_ip, and ->Binv. - * - * sizeof(->basis) == \sum_i nbf(i)*nhf(i) - * sizeof(->cell_ip) == \sum_i npairs(i)*ncells(i) - * sizeof(->Binv) == \sum_i nbf(i)^2 - * - * Does not fail. */ -/* ---------------------------------------------------------------------- */ -static void -compute_alloc_sizes(size_t nb, - struct coarse_sys_meta *m, - struct coarse_sys *sys, - size_t *bf_asz, size_t *ip_asz, size_t *Binv_asz) -/* ---------------------------------------------------------------------- */ -{ - size_t nf, nc, b; - - *bf_asz = *ip_asz = *Binv_asz = 0; - - for (b = 0; b < nb; b++) { - nf = sys->blkdof_pos[b + 1] - sys->blkdof_pos[b]; /* # coarse dof */ - nc = m ->pb2c [b + 1] - m ->pb2c [b]; /* # block c */ - - *bf_asz += nf * m->blk_nhf[b]; - *ip_asz += nf * (nf + 1) / 2 * nc; - *Binv_asz += nf * nf; - } -} - - -/* ---------------------------------------------------------------------- */ -/* Allocate a coarse sys structure as well as suitably sized - * individual data arrays within this structure. - * - * Returns fully allocated structure, with ->blkdof_pos and ->blkdof - * fully constructed if successful, and NULL if not. */ -/* ---------------------------------------------------------------------- */ -static struct coarse_sys * -coarse_sys_allocate(struct coarse_topology *ct, - struct coarse_sys_meta *m) -/* ---------------------------------------------------------------------- */ -{ - int alloc_ok; - size_t i, nb; - size_t bf_asz, ip_asz, Binv_asz; /* Allocation sizes */ - - struct coarse_sys *new; - - new = malloc(1 * sizeof *new); - - if (new != NULL) { - nb = ct->nblocks; - - assert(nb > 0); - - alloc_ok = blkdof_fill(ct, m, new); - - if (alloc_ok) { - compute_alloc_sizes(nb, m, new, &bf_asz, &ip_asz, &Binv_asz); - - new->dof2conn = malloc(m->n_act_bf * sizeof *new->dof2conn); - - new->basis_pos = malloc((nb + 1) * sizeof *new->basis_pos ); - new->cell_ip_pos = malloc((nb + 1) * sizeof *new->cell_ip_pos); - - new->basis = malloc(bf_asz * sizeof *new->basis ); - new->cell_ip = malloc(ip_asz * sizeof *new->cell_ip); - new->Binv = malloc(Binv_asz * sizeof *new->Binv ); - - alloc_ok += new->dof2conn != NULL; - alloc_ok += new->basis_pos != NULL; - alloc_ok += new->cell_ip_pos != NULL; - alloc_ok += new->basis != NULL; - alloc_ok += new->cell_ip != NULL; - alloc_ok += new->Binv != NULL; - } - - if (alloc_ok < 7) { - coarse_sys_destroy(new); - new = NULL; - } else { - for (i = 0; i <= nb; i++) { - new->basis_pos [i] = 0; - new->cell_ip_pos[i] = 0; - } - } - } - - return new; -} - - -/* ---------------------------------------------------------------------- */ -/* Map degree of freedom back to global grid connection (coarse face). - * In other words, fills previously allocated ->dof2conn array. - * - * Does not fail. */ -/* ---------------------------------------------------------------------- */ -static void -map_dof_to_conn(struct coarse_topology *ct, - struct coarse_sys_meta *m , - struct coarse_sys *sys) -/* ---------------------------------------------------------------------- */ -{ - int f, dof; - - for (f = 0; f < ct->nfaces; f++) { - dof = m->bfno[f]; - - if (dof >= 0) { - sys->dof2conn[ dof ] = f; - } - } -} - - -/* ---------------------------------------------------------------------- */ -/* Fill previously allocated ->basis_pos and ->cell_ip_pos arrays. - * See compute_alloc_sizes(). - * - * Does not fail. */ -/* ---------------------------------------------------------------------- */ -static void -set_csys_block_pointers(struct coarse_topology *ct, - struct coarse_sys_meta *m , - struct coarse_sys *sys) -/* ---------------------------------------------------------------------- */ -{ - int b, nb, nc, ndof, npairs; - - nb = ct->nblocks; - - sys->basis_pos[0] = sys->cell_ip_pos[0] = 0; - - for (b = 0; b < nb; b++) { - ndof = sys->blkdof_pos[b + 1] - sys->blkdof_pos[b]; - nc = m ->pb2c [b + 1] - m ->pb2c [b]; - - npairs = ndof * (ndof + 1) / 2; - - sys->basis_pos[b + 1] = sys->basis_pos[b] + ndof*m->blk_nhf[b]; - - sys->cell_ip_pos[b + 1] = sys->cell_ip_pos[b] + npairs*nc; - } -} - - -/* Create local numbering of the fine-scale faces contained in a pair - * of blocks denoted by 'cf'. - * - * Precondition: m->loc_fno[0 .. g->number_of_faces-1] < 0 - * - * Returns the number of local fine-scale faces. */ -/* ---------------------------------------------------------------------- */ -static int -enumerate_local_dofs(size_t cf, - struct UnstructuredGrid *g , - struct coarse_topology *ct, - struct coarse_sys_meta *m) -/* ---------------------------------------------------------------------- */ -{ - int *b, *c, i, f, loc_no; - - loc_no = 0; - - for (b = ct->neighbours + 2*(cf + 0); - b != ct->neighbours + 2*(cf + 1); b++) { - - if (*b >= 0) { - for (c = m->b2c + m->pb2c[*b + 0]; - c != m->b2c + m->pb2c[*b + 1]; c++) { - - for (i = g->cell_facepos[*c + 0]; - i < g->cell_facepos[*c + 1]; i++) { - - f = g->cell_faces[i]; - - if (m->loc_fno[f] < 0) { - m->loc_fno[f] = loc_no++; - } - } - } - } - } - - assert (loc_no > 0); - - return loc_no; -} - - -/* Destroy local numbering of fine-scale faces contained in a pair of - * blocks denoted by 'cf'. - * - * This is the inverse of enumerate_local_dofs(), and is needed to - * prepare assembly of the discrete local system defining the next - * basis function. */ -/* ---------------------------------------------------------------------- */ -static void -unenumerate_local_dofs(size_t cf, - struct UnstructuredGrid *g , - struct coarse_topology *ct, - struct coarse_sys_meta *m) -/* ---------------------------------------------------------------------- */ -{ - int *b, *c, i; - - for (b = ct->neighbours + 2*(cf + 0); - b != ct->neighbours + 2*(cf + 1); b++) { - if (*b >= 0) { - for (c = m->b2c + m->pb2c[*b + 0]; - c != m->b2c + m->pb2c[*b + 1]; c++) { - - for (i = g->cell_facepos[*c + 0]; - i < g->cell_facepos[*c + 1]; i++) { - - m->loc_fno[ g->cell_faces[i] ] = -1; - } - } - } - } -} - - -/* ---------------------------------------------------------------------- */ -/* Define local (to a single BF) pdof/dof CSR table. - * - * Precondition: m->loc_fno valid for BF (i.e., called after - * enumerate_local_dofs()). - * - * Does not fail. */ -/* ---------------------------------------------------------------------- */ -static void -linearise_local_dof(size_t cf, - struct UnstructuredGrid *g , - struct coarse_topology *ct, - struct coarse_sys_meta *m , - struct bf_asm_data *bf_asm) -/* ---------------------------------------------------------------------- */ -{ - int *b, *c, i; - int *pdof, *dof; - - pdof = bf_asm->pdof; - dof = bf_asm->dof; - - pdof[0] = 0; - - for (b = ct->neighbours + 2*(cf + 0); - b != ct->neighbours + 2*(cf + 1); b++) { - if (*b >= 0) { - for (c = m->b2c + m->pb2c[*b + 0]; - c != m->b2c + m->pb2c[*b + 1]; c++) { - - for (i = g->cell_facepos[*c + 0]; - i < g->cell_facepos[*c + 1]; i++) { - *dof++ = m->loc_fno[ g->cell_faces[i] ]; - } - - *++pdof = dof - bf_asm->dof; - } - } - } -} - - -/* ---------------------------------------------------------------------- */ -/* Construct coefficient matrix sparsity structure for single BF. - * - * Precondition: bf_asm->pdof and bf_asm->dof valid (i.e., called - * after linearise_local_dof()). - * - * Does not fail. */ -/* ---------------------------------------------------------------------- */ -static void -define_csr_sparsity(size_t nc, size_t m, struct bf_asm_data *bf_asm) -/* ---------------------------------------------------------------------- */ -{ - int p, *dof; - size_t c, i; - - MAT_SIZE_T n, i1, i2, dof1, dof2; - - struct CSRMatrix *A = bf_asm->A; - - /* ------------------------------------------------------------------ - * Count connections (O(m)) - * ------------------------------------------------------------------ */ - for (i = 0; i < m; i++) { A->ia[ i + 1 ] = 1; } /* Count self */ - - for (c = 0, p = 0; c < nc; c++) { - n = bf_asm->pdof[c + 1] - bf_asm->pdof[c]; - - for (; p < bf_asm->pdof[c + 1]; p++) { - assert ((size_t) bf_asm->dof[p] < m); - - A->ia[ bf_asm->dof[p] + 1 ] += n - 1; /* Count other */ - } - } - - /* ------------------------------------------------------------------ - * Define start pointers (O(m)) - * ------------------------------------------------------------------ */ - A->ia[0] = 0; - for (i = 1; i <= m; i++) { - A->ia[0] += A->ia[i]; - A->ia[i] = A->ia[0] - A->ia[i]; - } - - /* ------------------------------------------------------------------ - * Set sparsity (O(nnz)) - * ------------------------------------------------------------------ */ - A->ia[0] = 0; - for (i = 0; i < m; i++) { - A->ja[ A->ia[i + 1] ++ ] = (MAT_SIZE_T) i; /* Set self */ - } - - for (c = 0; c < nc; c++) { - n = bf_asm->pdof[c + 1] - bf_asm->pdof[c]; - - dof = bf_asm->dof + bf_asm->pdof[c]; - - for (i1 = 0; i1 < n; i1++) { - dof1 = dof[i1]; - - for (i2 = (i1 + 1) % n; i2 != i1; i2 = (i2 + 1) % n) { - dof2 = dof[i2]; - - A->ja[ A->ia[ dof1 + 1 ] ++ ] = dof2; - } - } - } - - A->m = m; - A->nnz = A->ia[m]; - - /* Enforce sorted connection structure per row */ - csrmatrix_sortrows(A); -} - - -/* ---------------------------------------------------------------------- */ -/* Assemble system of linear equations corresponding to local - * discretisation of flow problem on domain connected to coarse face - * 'cf'. The domain has a total of 'nlocf' fine-scale interfaces, and - * the BF weighting function 'w' is pre-calculated using function - * coarse_weight(). - * - * Does not fail. */ -/* ---------------------------------------------------------------------- */ -static void -assemble_local_system(size_t cf , - size_t nlocf, - struct UnstructuredGrid *g , - const double *Binv , - double *w , - struct coarse_topology *ct , - struct coarse_sys_meta *m , - struct bf_asm_data *bf_asm) -/* ---------------------------------------------------------------------- */ -{ - int c, i, p1, p2, ndof; - int *b, *dof; - size_t nc; - - double sgn; - - linearise_local_dof(cf, g, ct, m, bf_asm); - - nc = 0; - for (b = ct->neighbours + 2*(cf + 0); - b != ct->neighbours + 2*(cf + 1); b++) { - if (*b >= 0) { - nc += m->pb2c[*b + 1] - m->pb2c[*b]; - } - } - - define_csr_sparsity(nc, nlocf, bf_asm); - - csrmatrix_zero( bf_asm->A); - vector_zero (nlocf, bf_asm->b); - - sgn = 1.0; - dof = bf_asm->dof; - for (b = ct->neighbours + 2*(cf + 0); - b != ct->neighbours + 2*(cf + 1); b++) { - - if (*b >= 0) { - for (i = m->pb2c[*b]; i < m->pb2c[*b + 1]; i++) { - c = m->b2c[i]; - p1 = g->cell_facepos[c]; - p2 = m->pconn2[c]; - ndof = g->cell_facepos[c + 1] - p1; - - w[c] *= sgn; /* Set w-sign according to source/sink */ - - hybsys_cellcontrib_symm(c, ndof, p1, p2, bf_asm->gpress, - w, Binv, bf_asm->fsys); - - hybsys_global_assemble_cell(ndof, dof, bf_asm->fsys->S, - bf_asm->fsys->r, bf_asm->A, - bf_asm->b); - - w[c] *= sgn; /* Restore original weighting sign */ - - dof += ndof; - } - - sgn = -sgn; - } - } - - /* Account for zero eigenvalue of pure Neumann problem */ - bf_asm->A->sa[0] *= 2; -} - - -/* ---------------------------------------------------------------------- */ -/* Scale the fine-scale (inverse) inner product 'Binv' by the - * corresponding cell's total mobility. This includes mobility - * effects in the resulting BFs. */ -/* ---------------------------------------------------------------------- */ -static void -Binv_scale_mobility(int nc, struct coarse_sys_meta *m, - const double *totmob, double *Binv) -/* ---------------------------------------------------------------------- */ -{ - int c, i; - - for (c = i = 0; c < nc; c++) { - for (; i < m->pconn2[c]; i++) { - Binv[i] *= totmob[c]; - } - } -} - - -/* ---------------------------------------------------------------------- */ -/* Project BF flux values for coarse face 'cf' onto continuous flux - * field. */ -/* ---------------------------------------------------------------------- */ -static void -symmetrise_flux(size_t cf, struct UnstructuredGrid *g, struct coarse_topology *ct, - struct coarse_sys_meta *m, struct bf_asm_data *bf_asm) -/* ---------------------------------------------------------------------- */ -{ - int i, j, ndof, p1l, p1g, *b, *c, *dof, *cnt; - size_t f; - double s, blk_sgn, *flux; - - dof = bf_asm->dof; - flux = bf_asm->flux; - cnt = bf_asm->fcount; - - vector_zero(bf_asm->A->m, flux); - for (f = 0; f < bf_asm->A->m; f++) { - cnt[f] = 0; - } - - /* Accumulate (and count) number of fine-scale flux contributions - * from this particular BF. */ - i = 0; - for (b = ct->neighbours + 2*(cf + 0); - b != ct->neighbours + 2*(cf + 1); b++) { - - if (*b >= 0) { - for (c = m->b2c + m->pb2c[*b + 0]; - c != m->b2c + m->pb2c[*b + 1]; c++, i++) { - - p1l = bf_asm->pdof [ i]; - p1g = g->cell_facepos[*c]; - - ndof = bf_asm->pdof[i + 1] - p1l; - - assert (g->cell_facepos[*c + 1] - p1g == ndof); - - for (j = 0; j < ndof; j++) { - f = g->cell_faces[p1g + j]; - s = 2.0*(g->face_cells[2*f + 0] == *c) - 1.0; - - flux[dof[p1l + j]] += s * bf_asm->v[p1l + j]; - cnt [dof[p1l + j]] += 1; - } - } - } - } - - /* Arithmetic average. */ - for (f = 0; f < bf_asm->A->m; f++) { - flux[f] /= cnt[f]; - } - - /* Store symmetrised flux values back to bf_asm->v, with - * additional block outflow flux sign. */ - i = 0; blk_sgn = 1.0; - for (b = ct->neighbours + 2*(cf + 0); - b != ct->neighbours + 2*(cf + 1); b++) { - - if (*b >= 0) { - for (c = m->b2c + m->pb2c[*b + 0]; - c != m->b2c + m->pb2c[*b + 1]; c++, i++) { - - p1l = bf_asm->pdof [ i]; - p1g = g->cell_facepos[*c]; - - ndof = bf_asm->pdof[i + 1] - p1l; - - for (j = 0; j < ndof; j++) { - f = g->cell_faces[p1g + j]; - s = 2.0*(g->face_cells[2*f + 0] == *c) - 1.0; - - s *= blk_sgn; - - bf_asm->v[p1l + j] = s * flux[dof[p1l + j]]; - } - } - - blk_sgn = -blk_sgn; - } - } -} - - -/* ---------------------------------------------------------------------- */ -/* Solve local system of linear equations to derive interface - * pressures (Lagrange multipliers), then perform back-substitution to - * derive cell pressures and interface fluxes. The fluxes are the BF - * values on the coarse face denoted by 'cf'. - * - * Does not fail. */ -/* ---------------------------------------------------------------------- */ -static void -solve_local_system(size_t cf , - struct UnstructuredGrid *g , - const double *Binv , - struct coarse_topology *ct , - struct coarse_sys_meta *m , - struct bf_asm_data *bf_asm, - LocalSolver linsolve) -/* ---------------------------------------------------------------------- */ -{ - int i, j, ndof, p1l, p1g, *b, *c; - double a1, a2; - - MAT_SIZE_T incx, incy, nrows, ncols, lda; - - /* Compute interface pressures (solve system of lin. eqns.) */ - (*linsolve)(bf_asm->A, bf_asm->b, bf_asm->x); - - /* Back substitution to derive cell pressure/interface fluxes */ - incx = incy = 1; - - a1 = 1.0; - a2 = 0.0; - - i = p1l = p1g = 0; - for (b = ct->neighbours + 2*(cf + 0); - b != ct->neighbours + 2*(cf + 1); b++) { - - if (*b >= 0) { - for (c = m->b2c + m->pb2c[*b + 0]; - c != m->b2c + m->pb2c[*b + 1]; c++, i++) { - p1l = bf_asm->pdof[i + 0]; - ndof = bf_asm->pdof[i + 1] - p1l; - - nrows = ncols = lda = ndof; - - for (j = 0; j < ndof; j++) { - bf_asm->work[j] = bf_asm->x[ bf_asm->dof[ p1l + j ] ]; - } - - p1g = g->cell_facepos[*c]; - - bf_asm->p[i] = bf_asm->fsys->q[*c]; - bf_asm->p[i] += ddot_(&nrows, &bf_asm->fsys->F2[p1g], - &incx, bf_asm->work, &incy); - bf_asm->p[i] /= bf_asm->fsys->L[*c]; - - for (j = 0; j < ndof; j++) { - bf_asm->work[j] = bf_asm->p[i] - bf_asm->work[j]; - } - - dgemv_("No Transpose", &nrows, &ncols, - &a1, &Binv[m->pconn2[*c]], &lda, bf_asm->work, &incx, - &a2, &bf_asm->v[p1l] , &incy); - } - } - } - - symmetrise_flux(cf, g, ct, m, bf_asm); -} - - -/* ---------------------------------------------------------------------- */ -/* Store BF values at appropriate offsets into sys->basis. - * - * Must be called after solve_local_system(). - * - * Does not fail. */ -/* ---------------------------------------------------------------------- */ -static void -store_basis_function(size_t cf , - struct coarse_topology *ct , - struct coarse_sys_meta *m , - struct bf_asm_data *bf_asm, - struct coarse_sys *sys) -/* ---------------------------------------------------------------------- */ -{ - int i, loc_dofno, *b, *loc_dof; - ptrdiff_t sstart, dstart, nhf; - - loc_dof = m->loc_dofno + 2*(cf + 0); - - i = 0; - sstart = 0; - for (b = ct->neighbours + 2*(cf + 0); - b != ct->neighbours + 2*(cf + 1); b++, i++) { - - if (*b >= 0) { - loc_dofno = loc_dof[i]; - nhf = m->blk_nhf[*b]; - dstart = sys->basis_pos[*b] + loc_dofno*nhf; - - memcpy(sys->basis + dstart, - bf_asm->v + sstart, nhf * sizeof *sys->basis); - - sstart += nhf; - } - } -} - - -/* ====================================================================== - * Public interfaces below. - * ====================================================================== */ - - -/* ---------------------------------------------------------------------- */ -/* Construct coarse system from fine-scale grid (g), partition vector - * (p), coarse topology (ct), fine-scale permeability tensor (perm), - * fine-scale source terms (src), and fine-scale (total) mobility - * field (totmob). - * - * Uses 'linsolve' to resolve local systems of linear equations. - * - * Returns fully constructed coarse system if successful (i.e., if all - * internal allocations succeed and all BFs can be constructed), and - * NULL if not. */ -/* ---------------------------------------------------------------------- */ -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) -/* ---------------------------------------------------------------------- */ -{ - int cf; - size_t nlocf; - double *Binv, *w; - struct coarse_sys_meta *m; - struct bf_asm_data *bf_asm; - struct coarse_sys *sys; - - sys = NULL; bf_asm = NULL; Binv = NULL; w = NULL; - - m = coarse_sys_meta_construct(g, p, ct); - - if (m != NULL) { - Binv = compute_fs_ip(g, perm, m); - w = coarse_weight(g, ct->nblocks, p, m, perm, src); - bf_asm = bf_asm_data_allocate(g, m); - sys = coarse_sys_allocate(ct, m); - } - - if ((Binv != NULL) && (w != NULL) && - (bf_asm != NULL) && (sys != NULL)) { - - /* Provide reverse BF->face mapping for fs flux reconstruction */ - map_dof_to_conn(ct, m, sys); - - /* Prepare storage tables */ - set_csys_block_pointers(ct, m, sys); - - /* Exclude effects of gravity */ - vector_zero(g->cell_facepos[ g->number_of_cells ], bf_asm->gpress); - - /* Include mobility effects (multiple phases) */ - Binv_scale_mobility(g->number_of_cells, m, totmob, Binv); - - /* Discretise flow equation on fine scale */ - hybsys_schur_comp_symm(g->number_of_cells, g->cell_facepos, - Binv, bf_asm->fsys); - - for (cf = 0; cf < ct->nfaces; cf++) { - if (m->bfno[cf] >= 0) { - nlocf = enumerate_local_dofs(cf, g, ct, m); - - assemble_local_system(cf, nlocf, g, Binv, w, - ct, m, bf_asm); - - solve_local_system(cf, g, Binv, ct, m, - bf_asm, linsolve); - - store_basis_function(cf, ct, m, bf_asm, sys); - - unenumerate_local_dofs(cf, g, ct, m); - } - } - - coarse_sys_compute_cell_ip(g->number_of_cells, - m->max_ngconn, - ct->nblocks, - g->cell_facepos, - Binv, - m->pb2c, m->b2c, - sys); - } - - bf_asm_data_deallocate(bf_asm); - - free(w); free(Binv); - coarse_sys_meta_destroy(m); - - return sys; -} - - -/* ---------------------------------------------------------------------- */ -/* Release dynamic memory resources for coarse system data structure. */ -/* ---------------------------------------------------------------------- */ -void -coarse_sys_destroy(struct coarse_sys *sys) -/* ---------------------------------------------------------------------- */ -{ - if (sys != NULL) { - free(sys->Binv); - free(sys->cell_ip); - free(sys->basis); - - free(sys->cell_ip_pos); - free(sys->basis_pos); - - free(sys->blkdof); - free(sys->blkdof_pos); - - free(sys->dof2conn); - } - - free(sys); -} - - -/* ---------------------------------------------------------------------- */ -/* Compute \Psi'_i * B * \Psi_j for all basis function pairs (i,j) for - * all cells. Inverts inv(B) (i.e., Binv) in each cell. Iterates - * over blocks (CSR representation b2c_pos, b2c). Result store in - * sys->cell_ip, a packed representation of the IP pairs (one col per - * cell per block). - * - * Allocates work arrays and may fail. Does currently not report failure.*/ -/* ---------------------------------------------------------------------- */ -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) -/* ---------------------------------------------------------------------- */ -{ - int i, i1, i2, b, c, n, bf, *pconn2; - int max_nbf, nbf, loc_nc; - - size_t p, nbf_pairs, bf_off, bf_sz; - - MAT_SIZE_T mm, nn, kk, nrhs, ld1, ld2, ld3, info; - - double a1, a2; - double *work, *BI, *Psi, *IP; - -#if DEBUG_OUTPUT - FILE *fp; -#endif - - - max_nbf = max_diff(nb, sys->blkdof_pos); - - pconn2 = malloc((nc + 1) * sizeof *pconn2); - work = malloc(((max_nconn * max_nconn) + /* BI */ - (max_nconn * max_nbf ) + /* Psi */ - (max_nbf * max_nbf )) /* IP */ - * sizeof *work); - - if ((pconn2 != NULL) && (work != NULL)) { - BI = work + 0 ; - Psi = BI + (max_nconn * max_nconn); - IP = Psi + (max_nconn * max_nbf ); - - pconn2[0] = 0; - - for (i = 1; i <= nc; i++) { - n = pconn[i] - pconn[i - 1]; - pconn2[i] = pconn2[i - 1] + (n * n); - } - -#if DEBUG_OUTPUT - fp = fopen("debug_out.m", "wt"); -#endif - - for (b = 0; b < nb; b++) { - loc_nc = b2c_pos[b + 1] - b2c_pos[b]; - bf_off = 0; - nbf = sys->blkdof_pos[b + 1] - sys->blkdof_pos[b]; - - assert ((sys->basis_pos[b + 1] - - sys->basis_pos[b + 0]) % nbf == 0); - - bf_sz = (sys->basis_pos[b + 1] - sys->basis_pos[b]) / nbf; - - nbf_pairs = nbf * (nbf + 1) / 2; - - for (i = 0; i < loc_nc; i++) { - c = b2c[b2c_pos[b] + i]; - n = pconn[c + 1] - pconn[c]; - - /* Linearise (collect) BF values for cell */ - p = sys->basis_pos[b] + bf_off; - for (bf = 0; bf < nbf; bf++, p += bf_sz) { - memcpy(Psi + bf*n, &sys->basis[p], n * sizeof *Psi); - } -#if DEBUG_OUTPUT - fprintf(fp, "Psi{%d,%d} = [\n", b+1, i+1); - for (i2 = 0; i2 < nbf; i2++) { - for (i1 = 0; i1 < n; i1++) { - fprintf(fp, "%22.15e ", Psi[i1 + i2*n]); - } - fprintf(fp, ";\n"); - } - fprintf(fp, "].';\n"); -#endif - - /* Extract cell's inv(B) values... */ - memcpy(BI, &Binv[pconn2[c]], n * n * sizeof *BI); -#if DEBUG_OUTPUT - fprintf(fp, "BI{%d,%d} = [\n", b+1, i+1); - for (i2 = 0; i2 < n; i2++) { - for (i1 = 0; i1 < n; i1++) { - fprintf(fp, "%22.15e ", BI[i1 + i2*n]); - } - fprintf(fp, ";\n"); - } - fprintf(fp, "].';\n"); -#endif - - /* ...and (Cholesky) factor it... */ - nn = n; - ld1 = n; - dpotrf_("Upper Triangular", &nn, BI, &ld1, &info); - - /* ...and solve BI*X = Psi (=> Psi (=X) <- B*Psi) */ - nrhs = nbf; - ld2 = n; - dpotrs_("Upper Triangular", &nn, &nrhs, BI, &ld1, - Psi, &ld2, &info); -#if DEBUG_OUTPUT - fprintf(fp, "BPsi{%d,%d} = [\n", b+1, i+1); - for (i2 = 0; i2 < nbf; i2++) { - for (i1 = 0; i1 < n; i1++) { - fprintf(fp, "%22.15e ", Psi[i1 + i2*n]); - } - fprintf(fp, ";\n"); - } - fprintf(fp, "].';\n"); -#endif - - /* Finally, compute IP = Psi'*X = Psi'*B*Psi... */ - mm = nn = nbf; - kk = n; - ld1 = bf_sz; ld2 = n; ld3 = nbf; - a1 = 1.0; a2 = 0.0; - dgemm_("Transpose", "No Transpose", &mm, &nn, &kk, - &a1, &sys->basis[sys->basis_pos[b] + bf_off], &ld1, - Psi, &ld2, &a2, IP, &ld3); -#if DEBUG_OUTPUT - fprintf(fp, "PsiTBPsi{%d,%d} = [\n", b+1, i+1); - for (i2 = 0; i2 < nbf; i2++) { - for (i1 = 0; i1 < nbf; i1++) { - fprintf(fp, "%22.15e ", IP[i1 + i2*nbf]); - } - fprintf(fp, ";\n"); - } - fprintf(fp, "].';\n"); -#endif - - /* ...and fill results into ip-contrib for this cell... */ - p = sys->cell_ip_pos[b] + i*nbf_pairs; - for (i2 = 0; i2 < nbf; i2++) { - for (i1 = 0; i1 <= i2; i1++, p++) { - sys->cell_ip[p] = IP[i1 + i2*nbf]; - } - } - - /* ...and prepare for next cell. */ - bf_off += n; - } - } -#if DEBUG_OUTPUT - fclose(fp); -#endif - } - - free(work); free(pconn2); -} - - -/* ---------------------------------------------------------------------- */ -/* Compute inv(B) on coarse scale from fine-scale contributions. - * Specifically, this function computes the inverse of - * B=sum(1/lambda_c * B_c, c\in Blk_j) for all blocks, 'j'. The - * fine-scale B contributions are computed in - * coarse_sys_compute_cell_ip() above. - * - * Does not fail. */ -/* ---------------------------------------------------------------------- */ -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) -/* ---------------------------------------------------------------------- */ -{ - int b, i, i1, i2, nbf_pairs, loc_nc, nbf; - - double a1, a2; - double *Lti, *B; - - size_t p1, p2; - - MAT_SIZE_T mm, nn, ld, incx, incy, info; - -#if DEBUG_OUTPUT - FILE *fp; -#endif - - Lti = work + 0; - B = work + max_bcells; - - incx = incy = 1; - p2 = 0; - for (b = 0; b < nb; b++) { - loc_nc = b2c_pos[b + 1] - b2c_pos[b]; - - for (i = 0; i < loc_nc; i++) { - Lti[i] = 1.0 / totmob[b2c[b2c_pos[b] + i]]; - } - - /* Form coarse inner-product matrix for block 'b' as (inverse) - * mobility weighted sum of cell contributions */ - nbf = sys->blkdof_pos[b + 1] - sys->blkdof_pos[b]; - nbf_pairs = nbf * (nbf + 1) / 2; - - mm = ld = nbf_pairs; - nn = loc_nc; - a1 = 1.0; - a2 = 0.0; - dgemv_("No Transpose", &mm, &nn, &a1, - &sys->cell_ip[sys->cell_ip_pos[b]], &ld, - Lti, &incx, &a2, B, &incy); - - /* Factor (packed) SPD inner-product matrix... */ - mm = nbf; - dpptrf_("Upper Triangular", &mm, B, &info); - if (info == 0) { - /* ...and invert it... */ - dpptri_("Upper Triangular", &mm, B, &info); - } else { -#if DEBUG_OUTPUT - fp = fopen("debug_out.m", "at"); - mm = ld = nbf_pairs; - nn = loc_nc; - a1 = 1.0; - a2 = 0.0; - dgemv_("No Transpose", &mm, &nn, &a1, - &sys->cell_ip[sys->cell_ip_pos[b]], &ld, - Lti, &incx, &a2, B, &incy); - fprintf(fp, "IP{%d} = [\n", b + 1); - for (i2 = 0; i2 < nn; i2++) { - for (i1 = 0; i1 < mm; i1++) { - fprintf(fp, "%18.12e ", - sys->cell_ip[sys->cell_ip_pos[b] + i1 + i2*mm]); - } - fprintf(fp, ";\n"); - } - fprintf(fp, "].';\n"); - - fprintf(fp, "B{%d} = [ %% info = %lu\n", b + 1, - (unsigned long) info); - for (i1 = 0; i1 < nbf_pairs; i1++) - fprintf(fp, "%22.15e ", B[i1]); - fprintf(fp, "].';\n"); - fclose(fp); -#endif - } - - /* ...and write result to permanent storage suitable for - * reduction functions hybsys_schur_comp*() */ - p1 = 0; - for (i2 = 0; i2 < nbf; i2++) { - for (i1 = 0; i1 <= i2; i1++, p1++) { - sys->Binv[p2 + i1 + i2*nbf] = B[p1]; /* col i2 */ - sys->Binv[p2 + i2 + i1*nbf] = B[p1]; /* row i2 */ - } - } - - p2 += nbf * nbf; - } -} - - -/* ---------------------------------------------------------------------- */ -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) -/* ---------------------------------------------------------------------- */ -{ - int b, c1, c2, f, i, n; - const int *c; - - MAT_SIZE_T nrows, ncols, lda, incx, incy; - double a1, a2, s; - - vector_zero(G->number_of_faces, flux); - - incx = incy = 1; - a1 = 1.0; - a2 = 0.0; - - for (b = 0; b < ct->nblocks; b++) { - /* Construct fs hc fluxes in block (\Psi_b * v_c) */ - ncols = sys->blkdof_pos[b + 1] - - sys->blkdof_pos[b + 0]; /* ndof in block */ - nrows = sys->basis_pos [b + 1] - - sys->basis_pos [b + 0]; /* NUMEL(Psi(:)) */ - nrows /= ncols; - - lda = nrows; - dgemv_("No Transpose", &nrows, &ncols, &a1, - sys->basis + sys->basis_pos[b], - &lda, v_c + sys->blkdof_pos[b], - &incx, &a2, work, &incy); - - /* Accumulate fs interface fluxes (internal interfaces visited - * twice). */ - n = 0; - for (c = b2c + b2c_pos[b + 0]; - c != b2c + b2c_pos[b + 1]; c++) { - - for (i = G->cell_facepos[*c + 0]; - i < G->cell_facepos[*c + 1]; i++, n++) { - - f = G->cell_faces[i]; - s = 2.0*(G->face_cells[2*f + 0] == *c) - 1.0; - - flux[ f ] += s * work[ n ]; - } - } - } - - /* Symmetrise fine-scale flux */ - for (f = 0; f < G->number_of_faces; f++) { - c1 = G->face_cells[2*f + 0]; - c2 = G->face_cells[2*f + 1]; - - flux[f] /= (c1 >= 0) + (c2 >= 0); - } -} diff --git a/opm/core/pressure/msmfem/coarse_sys.h b/opm/core/pressure/msmfem/coarse_sys.h deleted file mode 100644 index 146a7df4..00000000 --- a/opm/core/pressure/msmfem/coarse_sys.h +++ /dev/null @@ -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 . -*/ - -#ifndef OPM_COARSE_SYS_HEADER_INCLUDED -#define OPM_COARSE_SYS_HEADER_INCLUDED - -#include - -#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 */ diff --git a/opm/core/pressure/msmfem/hash_set.c b/opm/core/pressure/msmfem/hash_set.c deleted file mode 100644 index 720ec241..00000000 --- a/opm/core/pressure/msmfem/hash_set.c +++ /dev/null @@ -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 . -*/ - -#include "config.h" -#include -#include -#include -#include -#include - -#include - -/* ====================================================================== - * 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<> 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; -} diff --git a/opm/core/pressure/msmfem/hash_set.h b/opm/core/pressure/msmfem/hash_set.h deleted file mode 100644 index 1edfd1ce..00000000 --- a/opm/core/pressure/msmfem/hash_set.h +++ /dev/null @@ -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 . -*/ - -#ifndef OPM_HASH_SET_HEADER_INCLUDED -#define OPM_HASH_SET_HEADER_INCLUDED - -#ifdef __cplusplus -extern "C" { -#endif - -#include - - -/* ---------------------------------------------------------------------- */ -/* Poor-man's unordered set (ind. key insert/all key extract only). */ -/* ---------------------------------------------------------------------- */ -struct hash_set { - size_t m; /* Table/set capacity (1<

. -*/ - -#include "config.h" -#include -#include -#include - -#include -#include - -#include -#include -#include - -#include -#include -#include -#include - - -#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 */ -} diff --git a/opm/core/pressure/msmfem/ifsh_ms.h b/opm/core/pressure/msmfem/ifsh_ms.h deleted file mode 100644 index 681d9231..00000000 --- a/opm/core/pressure/msmfem/ifsh_ms.h +++ /dev/null @@ -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 . -*/ - -#ifndef OPM_IFSH_MS_HEADER_INCLUDED -#define OPM_IFSH_MS_HEADER_INCLUDED - -#ifdef __cplusplus -extern "C" { -#endif - -#include - -#include -#include - -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 */ diff --git a/opm/core/pressure/tpfa/TransTpfa.cpp b/opm/core/pressure/tpfa/TransTpfa.cpp deleted file mode 100644 index 13e9f2d4..00000000 --- a/opm/core/pressure/tpfa/TransTpfa.cpp +++ /dev/null @@ -1,10 +0,0 @@ -#include "TransTpfa.hpp" -#include - -/* 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*); diff --git a/opm/core/pressure/tpfa/compr_bc.c b/opm/core/pressure/tpfa/compr_bc.c deleted file mode 100644 index 32c5bc13..00000000 --- a/opm/core/pressure/tpfa/compr_bc.c +++ /dev/null @@ -1,183 +0,0 @@ -/*=========================================================================== -// -// File: compr_bc.c -// -// Created: 2011-10-24 16:07:17+0200 -// -// Authors: Ingeborg S. Ligaarden -// Jostein R. Natvig -// Halvor M. Nilsen -// Atgeirr F. Rasmussen -// Bård Skaflestad -// -//==========================================================================*/ - - -/* - 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 . -*/ - -#include "config.h" -#include -#include -#include - -#include - -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; -} diff --git a/opm/core/pressure/tpfa/compr_bc.h b/opm/core/pressure/tpfa/compr_bc.h deleted file mode 100644 index e50ee4a3..00000000 --- a/opm/core/pressure/tpfa/compr_bc.h +++ /dev/null @@ -1,83 +0,0 @@ -/*=========================================================================== -// -// File: compr_bc.h -// -// Created: 2011-10-24 15:58:16+0200 -// -// Authors: Ingeborg S. Ligaarden -// Jostein R. Natvig -// Halvor M. Nilsen -// Atgeirr F. Rasmussen -// Bård Skaflestad -// -//==========================================================================*/ - - -/* - 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 . -*/ - -#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 */ diff --git a/opm/core/pressure/tpfa/compr_quant_general.c b/opm/core/pressure/tpfa/compr_quant_general.c deleted file mode 100644 index 4042e445..00000000 --- a/opm/core/pressure/tpfa/compr_quant_general.c +++ /dev/null @@ -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 . -*/ - -#include "config.h" -#include - -#include -#include - - -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; -} diff --git a/opm/core/pressure/tpfa/compr_quant_general.h b/opm/core/pressure/tpfa/compr_quant_general.h index 7027ee39..b74d09b8 100644 --- a/opm/core/pressure/tpfa/compr_quant_general.h +++ b/opm/core/pressure/tpfa/compr_quant_general.h @@ -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 - * - * compr_quantities_gen_deallocate(NULL) - * - * is supported and benign (i.e., this statement behaves as - * free(NULL)). - * - * @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 diff --git a/opm/core/pressure/tpfa/compr_source.c b/opm/core/pressure/tpfa/compr_source.c deleted file mode 100644 index eef23bb7..00000000 --- a/opm/core/pressure/tpfa/compr_source.c +++ /dev/null @@ -1,170 +0,0 @@ -/*=========================================================================== -// -// File: compr_source.c -// -// Created: 2011-10-19 19:22:24+0200 -// -// Authors: Ingeborg S. Ligaarden -// Jostein R. Natvig -// Halvor M. Nilsen -// Atgeirr F. Rasmussen -// Bård Skaflestad -// -//==========================================================================*/ - - -/* - 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 . -*/ - -#include "config.h" -#include -#include -#include - -#include - -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; -} diff --git a/opm/core/pressure/tpfa/compr_source.h b/opm/core/pressure/tpfa/compr_source.h index 9beb0964..37b05df3 100644 --- a/opm/core/pressure/tpfa/compr_source.h +++ b/opm/core/pressure/tpfa/compr_source.h @@ -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 np > 0 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 v > 0.0. - * @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., src->nsrc == 0) with an unchanged - * capacity. - */ -void -clear_compr_source_term(struct compr_src *src); - - #ifdef __cplusplus } #endif diff --git a/opm/core/transport/minimal/spu_explicit.c b/opm/core/transport/minimal/spu_explicit.c deleted file mode 100644 index 38a62b29..00000000 --- a/opm/core/transport/minimal/spu_explicit.c +++ /dev/null @@ -1,93 +0,0 @@ -/* - * Copyright 2010 (c) SINTEF ICT, Applied Mathematics. - * Jostein R. Natvig - */ - -#include "config.h" -#include -#include - -#include -#include - - -/* 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 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; fface_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; - } - } -} diff --git a/opm/core/transport/minimal/spu_explicit.h b/opm/core/transport/minimal/spu_explicit.h deleted file mode 100644 index 3d169f70..00000000 --- a/opm/core/transport/minimal/spu_explicit.h +++ /dev/null @@ -1,19 +0,0 @@ -/* - * Copyright 2010 (c) SINTEF ICT, Applied Mathematics. - * Jostein R. Natvig - */ - -#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 */ - diff --git a/opm/core/transport/minimal/spu_implicit.c b/opm/core/transport/minimal/spu_implicit.c deleted file mode 100644 index 1f1143c3..00000000 --- a/opm/core/transport/minimal/spu_implicit.c +++ /dev/null @@ -1,372 +0,0 @@ -/* - * Copyright 2010 (c) SINTEF ICT, Applied Mathematics. - * Jostein R. Natvig - */ - - -#include "config.h" -#include -#include -#include -#include -#include - -#include -#include - - - - -/* 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; inumber_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 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; i1.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; icell_facepos[i]; kcell_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; - } -} diff --git a/opm/core/transport/minimal/spu_implicit.h b/opm/core/transport/minimal/spu_implicit.h deleted file mode 100644 index 2851f920..00000000 --- a/opm/core/transport/minimal/spu_implicit.h +++ /dev/null @@ -1,30 +0,0 @@ -/* - * Copyright 2010 (c) SINTEF ICT, Applied Mathematics. - * Jostein R. Natvig - */ - -#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 */ -