OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'upstream/master'

Browse Source
This commit is contained in:
Markus Blatt 2012-11-22 12:48:02 +01:00
commit 87a5ba847e
20 changed files with 2129 additions and 149 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

177
examples/compute_tof_from_files.cpp Normal file
View File

@ -0,0 +1,177 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/grid.h>
#include <opm/core/GridManager.hpp>
#include <opm/core/newwells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/initState.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/transport/reorder/TransportModelTracerTof.hpp>
#include <opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/filesystem.hpp>
#include <algorithm>
#include <iostream>
#include <vector>
#include <numeric>
#include <iterator>
namespace
{
void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
{
if (param.anyUnused()) {
std::cout << "-------------------- Warning: unused parameters: --------------------\n";
param.displayUsage();
std::cout << "-------------------------------------------------------------------------" << std::endl;
}
}
} // anon namespace
// ----------------- Main program -----------------
int
main(int argc, char** argv)
{
using namespace Opm;
parameter::ParameterGroup param(argc, argv, false);
// Read grid.
GridManager grid_manager(param.get<std::string>("grid_filename"));
const UnstructuredGrid& grid = *grid_manager.c_grid();
// Read porosity, compute pore volume.
std::vector<double> porevol;
{
std::ifstream poro_stream(param.get<std::string>("poro_filename").c_str());
std::istream_iterator<double> beg(poro_stream);
std::istream_iterator<double> end;
porevol.assign(beg, end); // Now contains poro.
if (int(porevol.size()) != grid.number_of_cells) {
THROW("Size of porosity field differs from number of cells.");
}
for (int i = 0; i < grid.number_of_cells; ++i) {
porevol[i] *= grid.cell_volumes[i];
}
}
// Read flux.
std::vector<double> flux;
{
std::ifstream flux_stream(param.get<std::string>("flux_filename").c_str());
std::istream_iterator<double> beg(flux_stream);
std::istream_iterator<double> end;
flux.assign(beg, end);
if (int(flux.size()) != grid.number_of_faces) {
THROW("Size of flux field differs from number of faces.");
}
}
// Read source terms.
std::vector<double> src;
{
std::ifstream src_stream(param.get<std::string>("src_filename").c_str());
std::istream_iterator<double> beg(src_stream);
std::istream_iterator<double> end;
src.assign(beg, end);
if (int(src.size()) != grid.number_of_cells) {
THROW("Size of source term field differs from number of cells.");
}
}
// Choice of tof solver.
bool use_dg = param.getDefault("use_dg", false);
int dg_degree = -1;
bool use_cvi = false;
bool use_multidim_upwind = false;
if (use_dg) {
dg_degree = param.getDefault("dg_degree", 0);
use_cvi = param.getDefault("use_cvi", false);
} else {
use_multidim_upwind = param.getDefault("use_multidim_upwind", false);
}
// Write parameters used for later reference.
bool output = param.getDefault("output", true);
std::ofstream epoch_os;
std::string output_dir;
if (output) {
output_dir =
param.getDefault("output_dir", std::string("output"));
boost::filesystem::path fpath(output_dir);
try {
create_directories(fpath);
}
catch (...) {
THROW("Creating directories failed: " << fpath);
}
param.writeParam(output_dir + "/simulation.param");
}
// Issue a warning if any parameters were unused.
warnIfUnusedParams(param);
// Solve time-of-flight.
Opm::time::StopWatch transport_timer;
transport_timer.start();
std::vector<double> tof;
if (use_dg) {
Opm::TransportModelTracerTofDiscGal tofsolver(grid, use_cvi);
tofsolver.solveTof(&flux[0], &porevol[0], &src[0], dg_degree, tof);
} else {
Opm::TransportModelTracerTof tofsolver(grid, use_multidim_upwind);
tofsolver.solveTof(&flux[0], &porevol[0], &src[0], tof);
}
transport_timer.stop();
double tt = transport_timer.secsSinceStart();
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
// Output.
if (output) {
std::string tof_filename = output_dir + "/tof.txt";
std::ofstream tof_stream(tof_filename.c_str());
tof_stream.precision(16);
std::copy(tof.begin(), tof.end(), std::ostream_iterator<double>(tof_stream, "\n"));
}
}

2
opm/core/fluid/IncompPropertiesInterface.hpp
View File

@ -62,7 +62,7 @@ namespace Opm
/// \return Array of P viscosity values.
virtual const double* viscosity() const = 0;
/// Densities of fluid phases at surface conditions.
/// Densities of fluid phases at reservoir conditions.
/// \return Array of P density values.
virtual const double* density() const = 0;

225
opm/core/newwells.h
View File

@ -33,97 +33,177 @@
extern "C" {
#endif
/** Well type indicates desired/expected well behaviour. */
enum WellType { INJECTOR, PRODUCER };
/** Type of well control equation or inequality constraint.
* BHP -> Well constrained by bottom-hole pressure target.
* RESERVOIR_RATE -> Well constrained by reservoir volume flow rates.
* SURFACE_RATE -> Well constrained by surface volume flow rates.
/**
* Well type indicates desired/expected well behaviour.
*/
enum WellControlType { BHP, RESERVOIR_RATE, SURFACE_RATE };
enum WellType {
INJECTOR, /**< Well is an injector */
PRODUCER /**< Well is a producer */
};
/** Controls for a single well.
* Each control specifies a well rate or bottom-hole pressure. Only
* one control can be active at a time, indicated by current. The
* meaning of each control's target value depends on the control type:
* BHP -> target pressure in Pascal.
* RESERVOIR_RATE -> target reservoir volume rate in cubic(meter)/second
* SURFACE_RATE -> target surface volume rate in cubic(meter)/second
* The sign convention for RATE targets is as follows:
* (+) Fluid flowing into reservoir, i.e. injecting.
* (-) Fluid flowing out of reservoir, i.e. producing.
* For *_RATE controls, the distribution of phases used for the control
* is also needed. For example, a total rate control should have 1.0
* for each phase, whereas a control on oil rate should have 1.0 for
* the oil phase and 0.0 for the rest. For BHP controls, this is unused.
* The active control acts as an equality constraint, whereas the
* non-active controls should be interpreted as inequality
* constraints (upper or lower bounds). For instance, a PRODUCER's
* BHP constraint defines a minimum acceptable bottom-hole pressure
* value for the well.
/**
* Type of well control equation or inequality constraint.
*/
enum WellControlType {
BHP, /**< Well constrained by BHP target */
RESERVOIR_RATE, /**< Well constrained by reservoir volume flow rate */
SURFACE_RATE /**< Well constrained by surface volume flow rate */
};
/**
* Controls for a single well.
* Each control specifies a well rate or bottom-hole pressure. Only
* one control can be active at a time, indicated by current. The
* meaning of each control's target value depends on the control type:
*
* - BHP -> target pressure in Pascal.
* - RESERVOIR_RATE -> target reservoir volume rate in cubic(meter)/second
* - SURFACE_RATE -> target surface volume rate in cubic(meter)/second
*
* The sign convention for RATE targets is as follows:
*
* - (+) Fluid flowing into reservoir, i.e. injecting.
* - (-) Fluid flowing out of reservoir, i.e. producing.
*
* For *_RATE controls, the distribution of phases used for the control
* is also needed. For example, a total rate control should have 1.0
* for each phase, whereas a control on oil rate should have 1.0 for
* the oil phase and 0.0 for the rest. For BHP controls, this is unused.
* The active control acts as an equality constraint, whereas the
* non-active controls should be interpreted as inequality
* constraints (upper or lower bounds). For instance, a PRODUCER's
* BHP constraint defines a minimum acceptable bottom-hole pressure
* value for the well.
*/
struct WellControls
{
int num; /** Number of controls. */
enum WellControlType *type; /** Array of control types.*/
double *target; /** Array of control targets */
double *distr; /** Array of rate control distributions,
Wells::number_of_phases numbers per control */
int current; /** Index of current active control. */
/**
* Number of controls.
*/
int num;
void *data; /** Internal management structure. */
/**
* Array of control types.
*/
enum WellControlType *type;
/**
* Array of control targets
*/
double *target;
/**
* Array of rate control distributions,
* <CODE>Wells::number_of_phases</CODE> numbers for each control
*/
double *distr;
/**
* Index of current active control.
*/
int current;
/**
* Internal management structure.
*/
void *data;
};
/** Data structure aggregating static information about all wells in a scenario. */
/**
* Data structure aggregating static information about all wells in a scenario.
*/
struct Wells
{
int number_of_wells; /** Number of wells. */
int number_of_phases; /** Number of phases. */
int number_of_wells; /**< Number of wells. */
int number_of_phases; /**< Number of phases. */
enum WellType *type; /** Array of well types. */
double *depth_ref; /** Array of well bhp reference depths. */
double *comp_frac; /** Component fractions for each well, size is (number_of_wells*number_of_phases).
* This is intended to be used for injection wells. For production wells
* the component fractions will vary and cannot be specified a priori.
*/
int *well_connpos; /** Array of indices into well_cells (and WI).
* For a well w, well_connpos[w] and well_connpos[w+1] yield
* start and one-beyond-end indices into the well_cells array
* for accessing w's perforation cell indices.
*/
int *well_cells; /** Array of perforation cell indices.
* Size is number of perforations (== well_connpos[number_of_wells]).
*/
double *WI; /** Well productivity index, same size and structure as well_cells. */
struct WellControls **ctrls; /** Well controls, one set of controls for each well. */
/**
* Array of well types.
*/
enum WellType *type;
char **name; /** Well names. One string for each well. */
/**
* Array of well types.
*/
double *depth_ref;
void *data; /** Internal management structure. */
/**
* Component fractions for each well. Array of size
* <CODE>number_of_wells * number_of_phases</CODE>.
* This is intended to be used for injection wells. For production wells
* the component fractions will vary and cannot be specified a priori.
*/
double *comp_frac;
/**
* Array of indices into well_cells (and WI). For a well @c w,
* <CODE>well_connpos[w]</CODE> and <CODE>well_connpos[w+1]</CODE> are start
* and one-beyond-end indices into the @c well_cells array for accessing
* @c w's perforation cell indices.
*/
int *well_connpos;
/**
* Array of perforation cell indices.
* Size is number of perforations (== well_connpos[number_of_wells]).
*/
int *well_cells;
/**
* Well productivity index, same size and structure as well_cells.
*/
double *WI;
/**
* Well controls, one set of controls for each well.
*/
struct WellControls **ctrls;
/**
* Well names. One string for each well.
*/
char **name;
/**
* Internal management structure.
*/
void *data;
};
/** Data structure aggregating dynamic information about all wells in a scenario.
* All arrays in this structure contain data for each perforation,
* ordered the same as Wells::well_cells and Wells:WI. The array
* sizes are, respectively,
/**
* Data structure aggregating dynamic information about all wells in a scenario.
* All arrays in this structure contain data for each perforation, ordered the
* same as Wells::well_cells and Wells:WI. The array sizes are, respectively,
*
* gpot n*NP
* wdp NP
* A n²*NP (matrix in column-major (i.e., Fortran) order).
* phasemob n*NP
*
* in which "n" denotes the number of active fluid phases (and
* constituent components) and "NP" is the total number of
* perforations, <CODE>well_connpos[ number_of_wells ]</CODE>.
* in which "n" denotes the number of active fluid phases (and constituent
* components) and "NP" is the total number of perforations,
* <CODE>well_connpos[ number_of_wells ]</CODE>.
*/
struct CompletionData
{
double *gpot; /** Gravity potentials. */
double *A; /** Volumes to surface-components matrix, A = RB^{-1}. */
double *phasemob; /** Phase mobilities. */
/**
* Gravity potentials.
*/
double *wdp;
/**
* Volumes to surface-components matrix, A = RB^{-1}.
*/
double *A;
/**
* Phase mobilities for all perforations, stored consecutively with the
* phase index cycling the most rapidly.
*/
double *phasemob;
};
/**
@ -235,6 +315,19 @@ void
destroy_wells(struct Wells *W);
/**
* Create a deep-copy (i.e., clone) of an existing Wells object, including its
* controls.
*
* @param[in] W Existing Wells object.
* @return Complete clone of the input object. Dispose of resources using
* function destroy_wells() when no longer needed. Returns @c NULL in case of
* allocation failure.
*/
struct Wells *
clone_wells(const struct Wells *W);
#ifdef __cplusplus
}
#endif

32
opm/core/pressure/CompressibleTpfa.cpp
View File

@ -216,9 +216,9 @@ namespace Opm
const int nperf = wells_->well_connpos[nw];
const int dim = grid_.dimensions;
const double grav = gravity_ ? gravity_[dim - 1] : 0.0;
wellperf_gpot_.clear();
wellperf_gpot_.resize(np*nperf, 0.0);
if (grav == 0.0) {
wellperf_wdp_.clear();
wellperf_wdp_.resize(nperf, 0.0);
if (not (std::abs(grav) > 0.0)) {
return;
}
@ -228,9 +228,9 @@ namespace Opm
// Main loop, iterate over all perforations,
// using the following formula (by phase):
// gpot(perf) = g*(perf_z - well_ref_z)*rho(perf)
// where the phase densities rho(perf) are taken to be
// the densities in the perforation cell.
// wdp(perf) = g*(perf_z - well_ref_z)*rho(perf)
// where the total density rho(perf) is taken to be
// sum_p (rho_p*saturation_p) in the perforation cell.
for (int w = 0; w < nw; ++w) {
const double ref_depth = wells_->depth_ref[w];
for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w + 1]; ++j) {
@ -239,7 +239,8 @@ namespace Opm
props_.matrix(1, &state.pressure()[cell], &state.surfacevol()[np*cell], &cell, &A[0], 0);
props_.density(1, &A[0], &rho[0]);
for (int phase = 0; phase < np; ++phase) {
wellperf_gpot_[np*j + phase] = rho[phase]*grav*(cell_depth - ref_depth);
const double s_phase = state.saturation()[np*cell + phase];
wellperf_wdp_[j] += s_phase*rho[phase]*grav*(cell_depth - ref_depth);
}
}
}
@ -478,10 +479,7 @@ namespace Opm
std::copy(cM, cM + np, wpM);
} else {
const double bhp = well_state.bhp()[w];
double perf_p = bhp;
for (int phase = 0; phase < np; ++phase) {
perf_p += wellperf_gpot_[np*j + phase]*comp_frac[phase];
}
double perf_p = bhp + wellperf_wdp_[j];
// Hack warning: comp_frac is used as a component
// surface-volume variable in calls to matrix() and
// viscosity(), but as a saturation in the call to
@ -512,7 +510,7 @@ namespace Opm
const double* z = &state.surfacevol()[0];
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
CompletionData completion_data;
completion_data.gpot = ! wellperf_gpot_.empty() ? &wellperf_gpot_[0] : 0;
completion_data.wdp = ! wellperf_wdp_.empty() ? &wellperf_wdp_[0] : 0;
completion_data.A = ! wellperf_A_.empty() ? &wellperf_A_[0] : 0;
completion_data.phasemob = ! wellperf_phasemob_.empty() ? &wellperf_phasemob_[0] : 0;
cfs_tpfa_res_wells wells_tmp;
@ -599,7 +597,7 @@ namespace Opm
{
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
CompletionData completion_data;
completion_data.gpot = ! wellperf_gpot_.empty() ? const_cast<double*>(&wellperf_gpot_[0]) : 0;
completion_data.wdp = ! wellperf_wdp_.empty() ? const_cast<double*>(&wellperf_wdp_[0]) : 0;
completion_data.A = ! wellperf_A_.empty() ? const_cast<double*>(&wellperf_A_[0]) : 0;
completion_data.phasemob = ! wellperf_phasemob_.empty() ? const_cast<double*>(&wellperf_phasemob_[0]) : 0;
cfs_tpfa_res_wells wells_tmp;
@ -636,16 +634,10 @@ namespace Opm
// Compute well perforation pressures (not done by the C code).
if (wells_ != 0) {
const int nw = wells_->number_of_wells;
const int np = props_.numPhases();
for (int w = 0; w < nw; ++w) {
const double* comp_frac = &wells_->comp_frac[np*w];
for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) {
const double bhp = well_state.bhp()[w];
double perf_p = bhp;
for (int phase = 0; phase < np; ++phase) {
perf_p += wellperf_gpot_[np*j + phase]*comp_frac[phase];
}
well_state.perfPress()[j] = perf_p;
well_state.perfPress()[j] = bhp + wellperf_wdp_[j];
}
}
}

2
opm/core/pressure/CompressibleTpfa.hpp
View File

@ -134,7 +134,7 @@ namespace Opm
struct cfs_tpfa_res_data* h_;
// ------ Data that will be modified for every solve. ------
std::vector<double> wellperf_gpot_;
std::vector<double> wellperf_wdp_;
std::vector<double> initial_porevol_;
// ------ Data that will be modified for every solver iteration. ------

32
opm/core/pressure/tpfa/cfs_tpfa_residual.c
View File

@ -420,15 +420,15 @@ compute_well_compflux_and_deriv(struct cfs_tpfa_res_wells *wells ,
const double *wpress,
struct cfs_tpfa_res_impl *pimpl )
{
int c, i, w, np2;
double pw, dp;
double *WI, *gpot, *Ap, *pmobp;
double *pflux, *dpflux;
int c, i, w, np2;
double pw, dp;
const double *WI, *wdp, *Ap, *pmobp;
double *pflux, *dpflux, gpot[3] = { 0.0 };
struct Wells *W;
assert (wells->W != NULL);
assert (wells->W->number_of_phases <= 3);
W = wells->W;
@ -436,7 +436,7 @@ compute_well_compflux_and_deriv(struct cfs_tpfa_res_wells *wells ,
assert (W->data != NULL);
WI = W->WI;
gpot = wells->data->gpot;
wdp = wells->data->wdp;
Ap = wells->data->A;
pmobp = wells->data->phasemob;
@ -449,11 +449,11 @@ compute_well_compflux_and_deriv(struct cfs_tpfa_res_wells *wells ,
pw = wpress[w];
for (; i < W->well_connpos[w + 1]; i++,
gpot += np, Ap += np2, pmobp += np,
Ap += np2, pmobp += np,
pflux += np, dpflux += 2 * np) {
c = W->well_cells[i];
dp = pw - cpress[c];
dp = pw + wdp[i]- cpress[c];
compute_darcyflux_and_deriv(np, WI[i], dp, pmobp, gpot,
pimpl->flux_work,
@ -909,8 +909,8 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
{
int w, i, c, np, np2, nc;
int is_neumann, is_open;
double pw, dp;
double *WI, *gpot, *pmobp;
double pw, dp, gpot[3] = { 0.0 };
const double *WI, *wdp, *pmobp;
const double *Ac, *dAc;
struct Wells *W;
@ -923,7 +923,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
W = wells->W;
WI = W->WI;
gpot = wells->data->gpot;
wdp = wells->data->wdp;
pmobp = wells->data->phasemob;
is_neumann = 1;
@ -932,14 +932,13 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
pw = wpress[ w ];
is_open = W->ctrls[w]->current >= 0;
for (; i < W->well_connpos[w + 1];
i++, gpot += np, pmobp += np) {
for (; i < W->well_connpos[w + 1]; i++, pmobp += np) {
c = W->well_cells[ i ];
Ac = cq->Ac + (c * np2);
dAc = cq->dAc + (c * np2);
dp = pw - cpress[ c ];
dp = pw + wdp[i] - cpress[ c ];
init_completion_contrib(i, np, Ac, dAc, h->pimpl);
@ -1091,7 +1090,7 @@ compute_wflux(int np ,
{
int w, c, i, p;
double pw, dp, t;
const double *pmob;
const double *pmob, *wdp;
struct Wells *W;
struct CompletionData *cdata;
@ -1102,13 +1101,14 @@ compute_wflux(int np ,
W = wells->W;
cdata = wells->data;
wdp = cdata->wdp;
for (w = i = 0; w < W->number_of_wells; w++) {
pw = wpress[w];
for (; i < W->well_connpos[w + 1]; i++) {
c = W->well_cells[ i ];
dp = pw - cpress[c];
dp = pw + wdp[ i ] - cpress[c];
if (dp > 0) { pmob = cdata->phasemob + (i * np); } /* w->c */
else { pmob = pmobc + (c * np); } /* c->w */

238
opm/core/pressure/tpfa/cfs_tpfa_residual.h
View File

@ -25,6 +25,37 @@
#include <opm/core/pressure/tpfa/compr_source.h>
/**
* \file
* Public interface to assembler for (compressible) discrete pressure system
* based on two-point flux approximation method. The assembler implements a
* residual formulation that for a single cell \f$i\f$ reads
* \f[
* \mathit{pv}_i\cdot (1 - \sum_\alpha A_i^{-1}(p^{n+1},z^n)z^n) +
* \Delta t\sum_\alpha A_i^{-1}(p^{n+1},z^n) \Big(\sum_j A_{ij}v_{ij}^{n+1} -
* q_i^{n+1}\Big) = 0
* \f]
* in which \f$\mathit{pv}_i\f$ is the (constant or pressure-dependent)
* pore-volume of cell \f$i\f$. Moreover, \f$\Delta t\f$ is the time step size,
* \f$n\f$ denotes the time level, \f$A\f$ is the pressure and mass-dependent
* fluid matrix that converts phase volumes at reservoir conditions into
* component volumes at surface conditions and \f$v_{ij}\f$ is the vector of
* outward (with respect to cell \f$i\f$) phase fluxes across the \f$ij\f$
* cell-interface.
*
* This module's usage model is intended to be
* -# Construct assembler
* -# for (each time step)
* -# while (pressure not converged)
* -# Assemble appropriate (Jacobian) system of linear equations
* -# Solve Jacobian system to derive pressure increments
* -# Include increments into current state
* -# Check convergence
* -# Derive fluxes and, optionally, interface pressures
* -# Solve transport by some means
* -# Destroy assembler
*/
#ifdef __cplusplus
extern "C" {
#endif
@ -33,36 +64,121 @@ struct cfs_tpfa_res_impl;
struct CSRMatrix;
struct compr_quantities_gen;
/**
* Type encapsulating well topology and completion data (e.g., phase mobilities
* per connection (perforation)).
*/
struct cfs_tpfa_res_wells {
/**
* All wells pertaining to a particular linear system assembly.
* Must include current controls/targets and complete well topology.
*/
struct Wells *W ;
/**
* Completion data describing the fluid state at the current time level.
*/
struct CompletionData *data;
};
/**
* Type encapsulating all driving forces affecting the discrete pressure system.
*/
struct cfs_tpfa_res_forces {
struct cfs_tpfa_res_wells *wells;
struct compr_src *src ;
struct cfs_tpfa_res_wells *wells; /**< Wells */
struct compr_src *src ; /**< Explicit source terms */
};
/**
* Result structure that presents the fully assembled system of linear
* equations, linearised around the current pressure point.
*/
struct cfs_tpfa_res_data {
struct CSRMatrix *J;
double *F;
struct CSRMatrix *J; /**< Jacobian matrix */
double *F; /**< Residual vector (right-hand side) */
struct cfs_tpfa_res_impl *pimpl;
struct cfs_tpfa_res_impl *pimpl; /**< Internal management structure */
};
/**
* Construct assembler for system of linear equations.
*
* @param[in] G Grid
* @param[in] wells Well description. @c NULL in case of no wells.
* For backwards compatibility, the constructor also
* interprets <CODE>(wells != NULL) &&
* (wells->W == NULL)</CODE> as "no wells present", but new
* code should use <CODE>wells == NULL</CODE> to signify
* "no wells".
* @param[in] nphases Number of active fluid phases in this simulation run.
* Needed to correctly size various internal work arrays.
* @return Fully formed assembler structure suitable for forming systems of
* linear equations using, e.g., function cfs_tpfa_res_assemble(). @c NULL in
* case of allocation failure. Must be destroyed using function
* cfs_tpfa_res_destroy().
*/
struct cfs_tpfa_res_data *
cfs_tpfa_res_construct(struct UnstructuredGrid *G ,
struct cfs_tpfa_res_wells *wells ,
int nphases);
/**
* Destroy assembler for system of linear equations.
*
* Disposes of all resources acquired in a previous call to construction
* function cfs_tpfa_res_construct(). Note that the statement
* <CODE>
* cfs_tpfa_res_destroy(NULL)
* </CODE>
* is supported and benign (i.e., behaves like <CODE>free(NULL)</CODE>).
*
* @param[in,out] h On input - assembler obtained through a previous call to
* construction function cfs_tpfa_res_construct(). On output - invalid pointer.
*/
void
cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h);
/* Return value is 1 if the assembled matrix was adjusted to remove a
singularity. This happens if all fluids are incompressible and
there are no pressure conditions on wells or boundaries.
Otherwise return 0. */
/**
* Assemble system of linear equations by linearising the residual around the
* current pressure point. Assume incompressible rock (i.e., that the
* pore-volume is independent of pressure).
*
* The fully assembled system is presented in <CODE>h->J</CODE> and
* <CODE>h->F</CODE> and must be solved separately using external software.
*
* @param[in] G Grid.
* @param[in] dt Time step size \f$\Delta t\f$.
* @param[in] forces Driving forces.
* @param[in] zc Component volumes, per pore-volume, at surface
* conditions for all components in all cells stored
* consecutively per cell. Array of size
* <CODE>G->number_of_cells * cq->nphases</CODE>.
* @param[in] cq Compressible quantities describing the current fluid
* state. Fields @c Ac, @c dAc, @c Af, and
* @c phasemobf must be valid.
* @param[in] trans Background transmissibilities as defined by function
* tpfa_trans_compute().
* @param[in] gravcap_f Discrete gravity and capillary forces.
* @param[in] cpress Cell pressures. One scalar value per grid cell.
* @param[in] wpress Well (bottom-hole) pressures. One scalar value per
* well. @c NULL in case of no wells.
* @param[in] porevol Pore-volumes. One (positive) scalar value for each
* grid cell.
* @param[in,out] h On input-a valid (non-@c NULL) assembler obtained
* from a previous call to constructor function
* cfs_tpfa_res_construct(). On output-valid assembler
* that includes the new system of linear equations in
* its @c J and @c F fields.
*
* @return 1 if the assembled matrix was adjusted to remove a singularity. This
* happens if all fluids are incompressible and there are no pressure conditions
* on wells or boundaries. Otherwise return 0.
*/
int
cfs_tpfa_res_assemble(struct UnstructuredGrid *G,
double dt,
@ -76,11 +192,48 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G,
const double *porevol,
struct cfs_tpfa_res_data *h);
/* Return value is 1 if the assembled matrix was adjusted to remove a
singularity. This happens if all fluids are incompressible, the
rock is incompressible, and there are no pressure conditions on
wells or boundaries.
Otherwise return 0. */
/**
* Assemble system of linear equations by linearising the residual around the
* current pressure point. Assume compressible rock (i.e., that the pore-volume
* depends on pressure).
*
* The fully assembled system is presented in <CODE>h->J</CODE> and
* <CODE>h->F</CODE> and must be solved separately using external software.
*
* @param[in] G Grid.
* @param[in] dt Time step size \f$\Delta t\f$.
* @param[in] forces Driving forces.
* @param[in] zc Component volumes, per pore-volume, at surface
* conditions for all components in all cells stored
* consecutively per cell. Array of size
* <CODE>G->number_of_cells * cq->nphases</CODE>.
* @param[in] cq Compressible quantities describing the current fluid
* state. Fields @c Ac, @c dAc, @c Af, and
* @c phasemobf must be valid.
* @param[in] trans Background transmissibilities as defined by function
* tpfa_trans_compute().
* @param[in] gravcap_f Discrete gravity and capillary forces.
* @param[in] cpress Cell pressures. One scalar value per grid cell.
* @param[in] wpress Well (bottom-hole) pressures. One scalar value per
* well. @c NULL in case of no wells.
* @param[in] porevol Pore-volumes. One (positive) scalar value for each
* grid cell.
* @param[in] porevol0 Pore-volumes at start of time step (i.e., at time
* level \f$n\f$). One (positive) scalar value for
* each grid cell.
* @param[in] rock_comp Rock compressibility. One non-negative scalar for
* each grid cell.
* @param[in,out] h On input-a valid (non-@c NULL) assembler obtained
* from a previous call to constructor function
* cfs_tpfa_res_construct(). On output-valid assembler
* that includes the new system of linear equations in
* its @c J and @c F fields.
*
* @return 1 if the assembled matrix was adjusted to remove a singularity. This
* happens if all fluids are incompressible, the rock is incompressible, and
* there are no pressure conditions on wells or boundaries. Otherwise return 0.
*/
int
cfs_tpfa_res_comprock_assemble(
struct UnstructuredGrid *G,
@ -97,6 +250,38 @@ cfs_tpfa_res_comprock_assemble(
const double *rock_comp,
struct cfs_tpfa_res_data *h);
/**
* Derive interface (total) Darcy fluxes from (converged) pressure solution.
*
* @param[in] G Grid
* @param[in] forces Driving forces. Must correspond to those used when
* forming the system of linear equations, e.g., in the
* call to function cfs_tpfa_res_assemble().
* @param[in] np Number of fluid phases (and components).
* @param[in] trans Background transmissibilities as defined by function
* tpfa_trans_compute(). Must correspond to equally
* named parameter of the assembly functions.
* @param[in] pmobc Phase mobilities stored consecutively per cell with
* phase index cycling the most rapidly. Array of size
* <CODE>G->number_of_cells * np</CODE>.
* @param[in] pmobf Phase mobilities stored consecutively per interface
* with phase index cycling the most rapidly. Array of
* size <CODE>G->number_of_faces * np</CODE>.
* @param[in] gravcap_f Discrete gravity and capillary forces.
* @param[in] cpress Cell pressure values. One (positive) scalar for each
* grid cell.
* @param[in] wpress Well (bottom-hole) pressure values. One (positive)
* scalar value for each well. @c NULL in case of no
* wells.
* @param[out] fflux Total Darcy interface fluxes. One scalar value for
* each interface (inter-cell connection). Array of size
* <CODE>G->number_of_faces</CODE>.
* @param[out] wflux Total Darcy well connection fluxes. One scalar value
* for each well connection (perforation). Array of size
* <CODE>forces->wells->W->well_connpos
* [forces->wells->W->number_of_wells]</CODE>.
*/
void
cfs_tpfa_res_flux(struct UnstructuredGrid *G ,
struct cfs_tpfa_res_forces *forces ,
@ -110,6 +295,31 @@ cfs_tpfa_res_flux(struct UnstructuredGrid *G ,
double *fflux ,
double *wflux );
/**
* Derive interface pressures from (converged) pressure solution.
*
* @param[in] G Grid
* @param[in] np Number of fluid phases (and components).
* @param[in] htrans Background one-sided ("half") transmissibilities as
* defined by function tpfa_htrans_compute().
* @param[in] pmobf Phase mobilities stored consecutively per interface
* with phase index cycling the most rapidly. Array of
* size <CODE>G->number_of_faces * np</CODE>.
* @param[in] gravcap_f Discrete gravity and capillary forces.
* @param[in] h System assembler. Must correspond to the assembler
* state used to form the final system of linear equations
* from which the converged pressure solution was derived.
* @param[in] cpress Cell pressure values. One (positive) scalar for each
* grid cell.
* @param[in] fflux Total Darcy interface fluxes. One scalar value for
* each interface (inter-cell connection). Array of size
* <CODE>G->number_of_faces</CODE>. Typically computed
* using function cfs_tpfa_res_flux().
* @param[out] fpress Interface pressure values. One (positive) scalar for
* each interface. Array of size
* <CODE>G->number_of_faces</CODE>.
*/
void
cfs_tpfa_res_fpress(struct UnstructuredGrid *G,
int np,

98
opm/core/pressure/tpfa/compr_quant_general.h
View File

@ -22,23 +22,109 @@
#include <stddef.h>
/**
* \file
* Module defining derived fluid quantities needed to discretise compressible
* and miscible pressure (flow) problems.
*/
#ifdef __cplusplus
extern "C" {
#endif
/**
* Aggregate structure that represents an atomic view of the current fluid
* state. These quantities are used directly in the cfs_tpfa_residual module to
* discretise a particular, linearised flow problem.
*/
struct compr_quantities_gen {
int nphases; /* Number of phases/components */
/**
* Number of fluid phases. The pressure solvers also assume that the number
* of fluid components (at surface conditions) equals the number of fluid
* phases.
*/
int nphases;
double *Ac; /* RB^{-1} per cell */
double *dAc; /* d/dp (RB^{-1}) per cell */
double *Af; /* RB^{-1} per face */
double *phasemobf; /* Phase mobility per face */
double *voldiscr; /* Volume discrepancy per cell */
/**
* Pressure and mass-dependent fluid matrix that converts phase volumes at
* reservoir conditions into component volumes at surface conditions. Obeys
* the defining formula
* \f[
* A = RB^{-1}
* \f]
* in which \f$R\f$ denotes the miscibility ratios (i.e., the dissolved
* gas-oil ratio, \f$R_s\f$ and the evaporated oil-gas ratio, \f$R_v\f$)
* collected as a matrix and \f$B\f$ is the diagonal matrix of
* formation-volume expansion factors. The function is sampled in each grid
* cell. Array of size <CODE>nphases * nphases * nc</CODE>.
*/
double *Ac;
/**
* Derivative of \f$A\f$ with respect to pressure,
* \f[
* \frac{\partial A}{\partial p} = \frac{\partial R}{\partial p}B^{-1} +
* R\frac{\partial B^{-1}}{\partial p} = (\frac{\partial R}{\partial p} -
* A\frac{\partial B}{\partial p})B^{-1}
* \f]
* sampled in each grid cell. Array of size
* <CODE>nphases * nphases * nc</CODE>.
*/
double *dAc;
/**
* Fluid matrix sampled at each interface. Possibly as a result of an
* upwind calculation. Array of size <CODE>nphases * nphases * nf</CODE>.
*/
double *Af;
/**
* Phase mobility per interface. Possibly defined through an upwind
* calculation. Array of size <CODE>nphases * nf</CODE>.
*/
double *phasemobf;
/**
* Deceptively named "volume-discrepancy" term. Array of size @c nc.
* Unused by the cfs_tpfa_residual module.
*/
double *voldiscr;
};
/**
* Create management structure capable of storing derived fluid quantities
* pertaining to a reservoir model of a particular size.
*
* The resources acquired using function compr_quantities_gen_allocate() must
* be released using the destructor function compr_quantities_gen_deallocate().
*
* @param[in] nc Number of grid cells
* @param[in] nf Number of grid interfaces
* @param[in] np Number of fluid phases (and components)
*
* @return Fully formed management structure. Returns @c NULL in case of
* allocation failure.
*/
struct compr_quantities_gen *
compr_quantities_gen_allocate(size_t nc, size_t nf, int np);
/**
* Release resources acquired in a previous call to constructor function
* compr_quantities_gen_allocate().
*
* Note that
* <CODE>
* compr_quantities_gen_deallocate(NULL)
* </CODE>
* is supported and benign (i.e., this statement behaves as
* <CODE>free(NULL)</CODE>).
*
* @param[in,out] cq On input - compressible quantity management structure
* obtained through a previous call to construction function
* compr_quantities_gen_allocate(). On output - invalid pointer.
*/
void
compr_quantities_gen_deallocate(struct compr_quantities_gen *cq);

101
opm/core/pressure/tpfa/compr_source.h
View File

@ -36,28 +36,112 @@
#ifndef OPM_COMPR_SOURCE_H_HEADER
#define OPM_COMPR_SOURCE_H_HEADER
/**
* \file
* Data structures and support routines needed to represent explicit,
* compressible source terms.
*/
#ifdef __cplusplus
extern "C" {
#endif
/**
* Collection of explicit, compressible source terms.
*/
struct compr_src {
int nsrc;
int cpty;
/**
* Number of source terms.
*/
int nsrc;
int nphases;
/**
* Source term capacity. Client code should treat this member as read-only.
* The field is used in internal memory management.
*/
int cpty;
int *cell;
/**
* Number of fluid phases.
*/
int nphases;
/**
* Cells influenced by explicit source terms. Array of size @c cpty, the
* @c nsrc first elements (only) of which are valid.
*/
int *cell;
/**
* Total Darcy rate of inflow (measured at reservoir conditions) of each
* individual source term. Sign convention: Positive rate into reservoir
* (i.e., injection) and negative rate out of reservoir (production).
* Array of size @c cpty, the @c nsrc first elements (only) of which are
* valid.
*/
double *flux;
/**
* Injection composition for all explicit source terms. Not referenced for
* production sources (i.e., those terms for which <CODE>->flux[]</CODE> is
* negative). Array of size <CODE>nphases * cpty</CODE>, the
* <CODE>nphases * nsrc</CODE> of which (only) are valid.
*/
double *saturation;
};
/**
* Create a management structure that is, initially, capable of storing a
* specified number of sources defined by a particular number a fluid phases.
*
* @param[in] np Number of fluid phases. Must be positive to actually
* allocate any sources.
* @param[in] nsrc Initial management capacity. If positive, attempt to
* allocate that number of source terms. Otherwise, the
* initial capacity is treated as (and set to) zero.
* @return Fully formed management structure if <CODE>np > 0</CODE> and
* allocation success. @c NULL otherwise. The resources must be released using
* destructor function compr_src_deallocate().
*/
struct compr_src *
compr_src_allocate(int np, int nsrc);
/**
* Release memory resources acquired in a previous call to constructor function
* compr_src_allocate() and, possibly, source term insertion function
* append_compr_source_term().
*
* @param[in,out] src On input - source term management structure obtained
* through a previous call to construction function compr_src_allocate().
* On output - invalid pointer.
*/
void
compr_src_deallocate(struct compr_src *src);
/**
* Insert a new explicit source term into an existing collection.
*
* @param[in] c Cell influenced by this source term.
* @param[in] np Number of fluid phases. Used for consistency checking
* only.
* @param[in] v Source term total reservoir volume Darcy flux. Positive
* if the source term is to be interpreted as an injection
* source and negative otherwise.
* @param[in] sat Injection composition for this source term. Array of size
* @c np. Copied to internal storage, but the actual numbers
* are not inspected unless <CODE>v > 0.0</CODE>.
* @param[in,out] src On input - source term management structure obtained
* through a previous call to construction function compr_src_allocate() and,
* possibly, another call to function append_compr_source_term(). On output -
* source term collection that includes the new source term if successful and
* unchanged if not.
*
* @return One (1, true) if successful (i.e., if the new source term could be
* included in the existing collection) and zero (0, false) if not.
*/
int
append_compr_source_term(int c ,
int np ,
@ -65,6 +149,15 @@ append_compr_source_term(int c ,
const double *sat,
struct compr_src *src);
/**
* Empty source term collection while maintaining existing capacity.
*
* @param[in,out] src On input - an existing source term collection with a
* given number of sources and source capacity. On output - an empty source
* term collection (i.e., <CODE>src->nsrc == 0</CODE>) with an unchanged
* capacity.
*/
void
clear_compr_source_term(struct compr_src *src);

5
opm/core/transport/reorder/TransportModelInterface.cpp
View File

@ -20,6 +20,7 @@
#include <opm/core/transport/reorder/TransportModelInterface.hpp>
#include <opm/core/transport/reorder/reordersequence.h>
#include <opm/core/grid.h>
#include <opm/core/utility/StopWatch.hpp>
#include <vector>
#include <cassert>
@ -31,7 +32,11 @@ void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& g
std::vector<int> sequence(grid.number_of_cells);
std::vector<int> components(grid.number_of_cells + 1);
int ncomponents;
time::StopWatch clock;
clock.start();
compute_sequence(&grid, darcyflux, &sequence[0], &components[0], &ncomponents);
clock.stop();
std::cout << "Topological sort took: " << clock.secsSinceStart() << " seconds." << std::endl;
// Invoke appropriate solve method for each interdependent component.
for (int comp = 0; comp < ncomponents; ++comp) {

258
opm/core/transport/reorder/TransportModelTracerTof.cpp Normal file
View File

@ -0,0 +1,258 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/transport/reorder/TransportModelTracerTof.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <algorithm>
#include <numeric>
#include <cmath>
namespace Opm
{
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
TransportModelTracerTof::TransportModelTracerTof(const UnstructuredGrid& grid,
const bool use_multidim_upwind)
: grid_(grid), use_multidim_upwind_(use_multidim_upwind)
{
}
/// Solve for time-of-flight.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[out] tof Array of time-of-flight values.
void TransportModelTracerTof::solveTof(const double* darcyflux,
const double* porevolume,
const double* source,
std::vector<double>& tof)
{
darcyflux_ = darcyflux;
porevolume_ = porevolume;
source_ = source;
#ifndef NDEBUG
// Sanity check for sources.
const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0);
if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) {
THROW("Sources do not sum to zero: " << cum_src);
}
#endif
tof.resize(grid_.number_of_cells);
std::fill(tof.begin(), tof.end(), 0.0);
tof_ = &tof[0];
if (use_multidim_upwind_) {
face_tof_.resize(grid_.number_of_faces);
std::fill(face_tof_.begin(), face_tof_.end(), 0.0);
}
reorderAndTransport(grid_, darcyflux);
}
void TransportModelTracerTof::solveSingleCell(const int cell)
{
if (use_multidim_upwind_) {
solveSingleCellMultidimUpwind(cell);
return;
}
// Compute flux terms.
// Sources have zero tof, and therefore do not contribute
// to upwind_term. Sinks on the other hand, must be added
// to the downwind_flux (note sign change resulting from
// different sign conventions: pos. source is injection,
// pos. flux is outflow).
double upwind_term = 0.0;
double downwind_flux = std::max(-source_[cell], 0.0);
for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) {
int f = grid_.cell_faces[i];
double flux;
int other;
// Compute cell flux
if (cell == grid_.face_cells[2*f]) {
flux = darcyflux_[f];
other = grid_.face_cells[2*f+1];
} else {
flux =-darcyflux_[f];
other = grid_.face_cells[2*f];
}
// Add flux to upwind_term or downwind_flux, if interior.
if (other != -1) {
if (flux < 0.0) {
upwind_term += flux*tof_[other];
} else {
downwind_flux += flux;
}
}
}
// Compute tof.
tof_[cell] = (porevolume_[cell] - upwind_term)/downwind_flux;
}
void TransportModelTracerTof::solveSingleCellMultidimUpwind(const int cell)
{
// Compute flux terms.
// Sources have zero tof, and therefore do not contribute
// to upwind_term. Sinks on the other hand, must be added
// to the downwind terms (note sign change resulting from
// different sign conventions: pos. source is injection,
// pos. flux is outflow).
double upwind_term = 0.0;
double downwind_term_cell_factor = std::max(-source_[cell], 0.0);
double downwind_term_face = 0.0;
for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) {
int f = grid_.cell_faces[i];
double flux;
int other;
// Compute cell flux
if (cell == grid_.face_cells[2*f]) {
flux = darcyflux_[f];
other = grid_.face_cells[2*f+1];
} else {
flux =-darcyflux_[f];
other = grid_.face_cells[2*f];
}
// Add flux to upwind_term or downwind_flux, if interior.
if (other != -1) {
if (flux < 0.0) {
upwind_term += flux*face_tof_[f];
} else {
double fterm, cterm_factor;
multidimUpwindTerms(f, cell, fterm, cterm_factor);
downwind_term_face += fterm*flux;
downwind_term_cell_factor += cterm_factor*flux;
}
}
}
// Compute tof for cell.
tof_[cell] = (porevolume_[cell] - upwind_term - downwind_term_face)/downwind_term_cell_factor; // }
// Compute tof for downwind faces.
for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) {
int f = grid_.cell_faces[i];
const double outflux_f = (grid_.face_cells[2*f] == cell) ? darcyflux_[f] : -darcyflux_[f];
if (outflux_f > 0.0) {
double fterm, cterm_factor;
multidimUpwindTerms(f, cell, fterm, cterm_factor);
face_tof_[f] = fterm + cterm_factor*tof_[cell];
}
}
}
void TransportModelTracerTof::solveMultiCell(const int num_cells, const int* cells)
{
std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl;
for (int i = 0; i < num_cells; ++i) {
solveSingleCell(cells[i]);
}
}
// Assumes that face_tof_[f] is known for all upstream faces f of upwind_cell.
// Assumes that darcyflux_[face] is != 0.0.
// This function returns factors to compute the tof for 'face':
// tof(face) = face_term + cell_term_factor*tof(upwind_cell).
// It is not computed here, since these factors are needed to
// compute the tof(upwind_cell) itself.
void TransportModelTracerTof::multidimUpwindTerms(const int face,
const int upwind_cell,
double& face_term,
double& cell_term_factor) const
{
// Implements multidim upwind according to
// "Multidimensional upstream weighting for multiphase transport on general grids"
// by Keilegavlen, Kozdon, Mallison.
// However, that article does not give a 3d extension other than noting that using
// multidimensional upwinding in the XY-plane and not in the Z-direction may be
// a good idea. We have here attempted some generalization, by looking at all
// face-neighbours across edges as upwind candidates, and giving them all uniform weight.
// This will over-weight the immediate upstream cell value in an extruded 2d grid with
// one layer (top and bottom no-flow faces will enter the computation) compared to the
// original 2d case. Improvements are welcome.
// Identify the adjacent faces of the upwind cell.
const int* face_nodes_beg = grid_.face_nodes + grid_.face_nodepos[face];
const int* face_nodes_end = grid_.face_nodes + grid_.face_nodepos[face + 1];
ASSERT(face_nodes_end - face_nodes_beg == 2 || grid_.dimensions != 2);
adj_faces_.clear();
for (int hf = grid_.cell_facepos[upwind_cell]; hf < grid_.cell_facepos[upwind_cell + 1]; ++hf) {
const int f = grid_.cell_faces[hf];
if (f != face) {
const int* f_nodes_beg = grid_.face_nodes + grid_.face_nodepos[f];
const int* f_nodes_end = grid_.face_nodes + grid_.face_nodepos[f + 1];
// Find out how many vertices they have in common.
// Using simple linear searches since sets are small.
int num_common = 0;
for (const int* f_iter = f_nodes_beg; f_iter < f_nodes_end; ++f_iter) {
num_common += std::count(face_nodes_beg, face_nodes_end, *f_iter);
}
if (num_common == grid_.dimensions - 1) {
// Neighbours over an edge (3d) or vertex (2d).
adj_faces_.push_back(f);
} else {
ASSERT(num_common == 0);
}
}
}
// Indentify adjacent faces with inflows, compute omega_star, omega,
// add up contributions.
const int num_adj = adj_faces_.size();
ASSERT(num_adj == face_nodes_end - face_nodes_beg);
const double flux_face = std::fabs(darcyflux_[face]);
face_term = 0.0;
cell_term_factor = 0.0;
for (int ii = 0; ii < num_adj; ++ii) {
const int f = adj_faces_[ii];
const double influx_f = (grid_.face_cells[2*f] == upwind_cell) ? -darcyflux_[f] : darcyflux_[f];
const double omega_star = influx_f/flux_face;
// SPU
// const double omega = 0.0;
// TMU
// const double omega = omega_star > 0.0 ? std::min(omega_star, 1.0) : 0.0;
// SMU
const double omega = omega_star > 0.0 ? omega_star/(1.0 + omega_star) : 0.0;
face_term += omega*face_tof_[f];
cell_term_factor += (1.0 - omega);
}
face_term /= double(num_adj);
cell_term_factor /= double(num_adj);
}
} // namespace Opm

83
opm/core/transport/reorder/TransportModelTracerTof.hpp Normal file
View File

@ -0,0 +1,83 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED
#define OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED
#include <opm/core/transport/reorder/TransportModelInterface.hpp>
#include <vector>
#include <map>
#include <ostream>
struct UnstructuredGrid;
namespace Opm
{
class IncompPropertiesInterface;
/// Implements a first-order finite volume solver for
/// (single-phase) time-of-flight using reordering.
/// The equation solved is:
/// v \cdot \grad\tau = \phi
/// where v is the fluid velocity, \tau is time-of-flight and
/// \phi is the porosity. This is a boundary value problem, where
/// \tau is specified to be zero on all inflow boundaries.
class TransportModelTracerTof : public TransportModelInterface
{
public:
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] use_multidim_upwind If true, use multidimensional tof upwinding.
TransportModelTracerTof(const UnstructuredGrid& grid,
const bool use_multidim_upwind = false);
/// Solve for time-of-flight.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[out] tof Array of time-of-flight values.
void solveTof(const double* darcyflux,
const double* porevolume,
const double* source,
std::vector<double>& tof);
private:
virtual void solveSingleCell(const int cell);
void solveSingleCellMultidimUpwind(const int cell);
virtual void solveMultiCell(const int num_cells, const int* cells);
void multidimUpwindTerms(const int face, const int upwind_cell,
double& face_term, double& cell_term_factor) const;
private:
const UnstructuredGrid& grid_;
const double* darcyflux_; // one flux per grid face
const double* porevolume_; // one volume per cell
const double* source_; // one volumetric source term per cell
double* tof_;
bool use_multidim_upwind_;
std::vector<double> face_tof_; // For multidim upwind face tofs.
mutable std::vector<int> adj_faces_; // For multidim upwind logic.
};
} // namespace Opm
#endif // OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED

657
opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp Normal file
View File

@ -0,0 +1,657 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/VelocityInterpolation.hpp>
#include <opm/core/linalg/blas_lapack.h>
#include <algorithm>
#include <cmath>
#include <numeric>
namespace Opm
{
// --------------- Helpers for TransportModelTracerTofDiscGal ---------------
/// A class providing discontinuous Galerkin basis functions.
struct DGBasis
{
static int numBasisFunc(const int dimensions,
const int degree)
{
switch (dimensions) {
case 1:
return degree + 1;
case 2:
return (degree + 2)*(degree + 1)/2;
case 3:
return (degree + 3)*(degree + 2)*(degree + 1)/6;
default:
THROW("Dimensions must be 1, 2 or 3.");
}
}
/// Evaluate all nonzero basis functions at x,
/// writing to f_x. The array f_x must have
/// size numBasisFunc(grid.dimensions, degree).
///
/// The basis functions are the following
/// Degree 0: 1.
/// Degree 1: x - xc, y - yc, z - zc etc.
/// Further degrees await development.
static void eval(const UnstructuredGrid& grid,
const int cell,
const int degree,
const double* x,
double* f_x)
{
const int dim = grid.dimensions;
const double* cc = grid.cell_centroids + dim*cell;
// Note intentional fallthrough in this switch statement!
switch (degree) {
case 1:
for (int ix = 0; ix < dim; ++ix) {
f_x[1 + ix] = x[ix] - cc[ix];
}
case 0:
f_x[0] = 1;
break;
default:
THROW("Maximum degree is 1 for now.");
}
}
/// Evaluate gradients of all nonzero basis functions at x,
/// writing to grad_f_x. The array grad_f_x must have size
/// numBasisFunc(grid.dimensions, degree) * grid.dimensions.
/// The <grid.dimensions> components of the first basis function
/// gradient come before the components of the second etc.
static void evalGrad(const UnstructuredGrid& grid,
const int /*cell*/,
const int degree,
const double* /*x*/,
double* grad_f_x)
{
const int dim = grid.dimensions;
const int num_basis = numBasisFunc(dim, degree);
std::fill(grad_f_x, grad_f_x + num_basis*dim, 0.0);
if (degree > 1) {
THROW("Maximum degree is 1 for now.");
} else if (degree == 1) {
for (int ix = 0; ix < dim; ++ix) {
grad_f_x[dim*(ix + 1) + ix] = 1.0;
}
}
}
};
static void cross(const double* a, const double* b, double* res)
{
res[0] = a[1]*b[2] - a[2]*b[1];
res[1] = a[2]*b[0] - a[0]*b[2];
res[2] = a[0]*b[1] - a[1]*b[0];
}
static double triangleArea3d(const double* p0,
const double* p1,
const double* p2)
{
double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] };
double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] };
double cr[3];
cross(a, b, cr);
return 0.5*std::sqrt(cr[0]*cr[0] + cr[1]*cr[1] + cr[2]*cr[2]);
}
/// Calculates the determinant of a 3 x 3 matrix, represented as
/// three three-dimensional arrays.
static double determinantOf(const double* a0,
const double* a1,
const double* a2)
{
return
a0[0] * (a1[1] * a2[2] - a2[1] * a1[2]) -
a0[1] * (a1[0] * a2[2] - a2[0] * a1[2]) +
a0[2] * (a1[0] * a2[1] - a2[0] * a1[1]);
}
/// Computes the volume of a tetrahedron consisting of 4 vertices
/// with 3-dimensional coordinates
static double tetVolume(const double* p0,
const double* p1,
const double* p2,
const double* p3)
{
double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] };
double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] };
double c[3] = { p3[0] - p0[0], p3[1] - p0[1], p3[2] - p0[2] };
return std::fabs(determinantOf(a, b, c) / 6.0);
}
/// A class providing numerical quadrature for cells.
/// In general: \int_{cell} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i).
/// Note that this class does multiply weights by cell volume,
/// so weights always sum to cell volume.
/// Degree 1 method:
/// Midpoint (centroid) method.
/// n = 1, w_0 = cell volume, x_0 = cell centroid
/// Degree 2 method:
/// Based on subdivision of each cell face into triangles
/// with the face centroid as a common vertex, and then
/// subdividing the cell into tetrahedra with the cell
/// centroid as a common vertex. Then apply the tetrahedron
/// rule with the following 4 nodes (uniform weights):
/// a = 0.138196601125010515179541316563436
/// x_i has all barycentric coordinates = a, except for
/// the i'th coordinate which is = 1 - 3a.
/// This rule is from http://nines.cs.kuleuven.be/ecf,
/// it is the second degree 2 4-point rule for tets,
/// referenced to Stroud(1971).
/// The tetrahedra are numbered T_{i,j}, and are given by the
/// cell centroid C, the face centroid FC_i, and two nodes
/// of face i: FN_{i,j}, FN_{i,j+1}.
class CellQuadrature
{
public:
CellQuadrature(const UnstructuredGrid& grid,
const int cell,
const int degree)
: grid_(grid), cell_(cell), degree_(degree)
{
if (degree > 2) {
THROW("CellQuadrature exact for polynomial degrees > 1 not implemented.");
}
if (degree == 2) {
// Prepare subdivision.
}
}
int numQuadPts() const
{
if (degree_ < 2) {
return 1;
}
// Degree 2 case.
int sumnodes = 0;
for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
const int face = grid_.cell_faces[hf];
sumnodes += grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
}
return 4*sumnodes;
}
void quadPtCoord(const int index, double* coord) const
{
const int dim = grid_.dimensions;
const double* cc = grid_.cell_centroids + dim*cell_;
if (degree_ < 2) {
std::copy(cc, cc + dim, coord);
return;
}
// Degree 2 case.
int tetindex = index / 4;
const int subindex = index % 4;
const double* nc = grid_.node_coordinates;
for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
const int face = grid_.cell_faces[hf];
const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
if (nfn <= tetindex) {
// Our tet is not associated with this face.
tetindex -= nfn;
continue;
}
const double* fc = grid_.face_centroids + dim*face;
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face];
const int node0 = fnodes[tetindex];
const int node1 = fnodes[(tetindex + 1) % nfn];
const double* n0c = nc + dim*node0;
const double* n1c = nc + dim*node1;
const double a = 0.138196601125010515179541316563436;
// Barycentric coordinates of our point in the tet.
double baryc[4] = { a, a, a, a };
baryc[subindex] = 1.0 - 3.0*a;
for (int dd = 0; dd < dim; ++dd) {
coord[dd] = baryc[0]*cc[dd] + baryc[1]*fc[dd] + baryc[2]*n0c[dd] + baryc[3]*n1c[dd];
}
return;
}
THROW("Should never reach this point.");
}
double quadPtWeight(const int index) const
{
if (degree_ < 2) {
return grid_.cell_volumes[cell_];
}
// Degree 2 case.
const int dim = grid_.dimensions;
const double* cc = grid_.cell_centroids + dim*cell_;
int tetindex = index / 4;
const double* nc = grid_.node_coordinates;
for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
const int face = grid_.cell_faces[hf];
const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
if (nfn <= tetindex) {
// Our tet is not associated with this face.
tetindex -= nfn;
continue;
}
const double* fc = grid_.face_centroids + dim*face;
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face];
const int node0 = fnodes[tetindex];
const int node1 = fnodes[(tetindex + 1) % nfn];
const double* n0c = nc + dim*node0;
const double* n1c = nc + dim*node1;
return 0.25*tetVolume(cc, fc, n0c, n1c);
}
THROW("Should never reach this point.");
}
private:
const UnstructuredGrid& grid_;
const int cell_;
const int degree_;
};
/// A class providing numerical quadrature for faces.
/// In general: \int_{face} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i).
/// Note that this class does multiply weights by face area,
/// so weights always sum to face area.
/// Degree 1 method:
/// Midpoint (centroid) method.
/// n = 1, w_0 = face area, x_0 = face centroid
/// Degree 2 method:
/// Based on subdivision of the face into triangles,
/// with the centroid as a common vertex, and the triangle
/// edge midpoint rule.
/// Triangle i consists of the centroid C, nodes N_i and N_{i+1}.
/// Its area is A_i.
/// n = 2 * nn (nn = num nodes in face)
/// For i = 0..(nn-1):
/// w_i = 1/3 A_i.
/// w_{nn+i} = 1/3 A_{i-1} + 1/3 A_i
/// x_i = (N_i + N_{i+1})/2
/// x_{nn+i} = (C + N_i)/2
/// All N and A indices are interpreted cyclic, modulus nn.
class FaceQuadrature
{
public:
FaceQuadrature(const UnstructuredGrid& grid,
const int face,
const int degree)
: grid_(grid), face_(face), degree_(degree)
{
if (grid_.dimensions != 3) {
THROW("FaceQuadrature only implemented for 3D case.");
}
if (degree_ > 2) {
THROW("FaceQuadrature exact for polynomial degrees > 2 not implemented.");
}
}
int numQuadPts() const
{
if (degree_ < 2) {
return 1;
}
// Degree 2 case.
return 2 * (grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]);
}
void quadPtCoord(const int index, double* coord) const
{
const int dim = grid_.dimensions;
const double* fc = grid_.face_centroids + dim*face_;
if (degree_ < 2) {
std::copy(fc, fc + dim, coord);
return;
}
// Degree 2 case.
const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_];
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_];
const double* nc = grid_.node_coordinates;
if (index < nn) {
// Boundary edge midpoint.
const int node0 = fnodes[index];
const int node1 = fnodes[(index + 1)%nn];
for (int dd = 0; dd < dim; ++dd) {
coord[dd] = 0.5*(nc[dim*node0 + dd] + nc[dim*node1 + dd]);
}
} else {
// Interiour edge midpoint.
// Recall that index is now in [nn, 2*nn).
const int node = fnodes[index - nn];
for (int dd = 0; dd < dim; ++dd) {
coord[dd] = 0.5*(nc[dim*node + dd] + fc[dd]);
}
}
}
double quadPtWeight(const int index) const
{
if (degree_ < 2) {
return grid_.face_areas[face_];
}
// Degree 2 case.
const int dim = grid_.dimensions;
const double* fc = grid_.face_centroids + dim*face_;
const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_];
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_];
const double* nc = grid_.node_coordinates;
if (index < nn) {
// Boundary edge midpoint.
const int node0 = fnodes[index];
const int node1 = fnodes[(index + 1)%nn];
const double area = triangleArea3d(nc + dim*node1, nc + dim*node0, fc);
return area / 3.0;
} else {
// Interiour edge midpoint.
// Recall that index is now in [nn, 2*nn).
const int node0 = fnodes[(index - 1) % nn];
const int node1 = fnodes[index - nn];
const int node2 = fnodes[(index + 1) % nn];
const double area0 = triangleArea3d(nc + dim*node1, nc + dim*node0, fc);
const double area1 = triangleArea3d(nc + dim*node2, nc + dim*node1, fc);
return (area0 + area1) / 3.0;
}
}
private:
const UnstructuredGrid& grid_;
const int face_;
const int degree_;
};
// --------------- Methods of TransportModelTracerTofDiscGal ---------------
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] use_cvi If true, use corner point velocity interpolation.
/// Otherwise, use the basic constant interpolation.
TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid,
const bool use_cvi)
: grid_(grid),
coord_(grid.dimensions),
velocity_(grid.dimensions)
{
if (use_cvi) {
velocity_interpolation_.reset(new VelocityInterpolationECVI(grid));
} else {
velocity_interpolation_.reset(new VelocityInterpolationConstant(grid));
}
}
/// Solve for time-of-flight.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[in] degree Polynomial degree of DG basis functions used.
/// \param[out] tof_coeff Array of time-of-flight solution coefficients.
/// The values are ordered by cell, meaning that
/// the K coefficients corresponding to the first
/// cell comes before the K coefficients corresponding
/// to the second cell etc.
/// K depends on degree and grid dimension.
void TransportModelTracerTofDiscGal::solveTof(const double* darcyflux,
const double* porevolume,
const double* source,
const int degree,
std::vector<double>& tof_coeff)
{
darcyflux_ = darcyflux;
porevolume_ = porevolume;
source_ = source;
#ifndef NDEBUG
// Sanity check for sources.
const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0);
if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) {
THROW("Sources do not sum to zero: " << cum_src);
}
#endif
degree_ = degree;
const int num_basis = DGBasis::numBasisFunc(grid_.dimensions, degree_);
tof_coeff.resize(num_basis*grid_.number_of_cells);
std::fill(tof_coeff.begin(), tof_coeff.end(), 0.0);
tof_coeff_ = &tof_coeff[0];
rhs_.resize(num_basis);
jac_.resize(num_basis*num_basis);
orig_jac_.resize(num_basis*num_basis);
basis_.resize(num_basis);
basis_nb_.resize(num_basis);
grad_basis_.resize(num_basis*grid_.dimensions);
velocity_interpolation_->setupFluxes(darcyflux);
reorderAndTransport(grid_, darcyflux);
}
void TransportModelTracerTofDiscGal::solveSingleCell(const int cell)
{
// Residual:
// For each cell K, basis function b_j (spanning V_h),
// writing the solution u_h|K = \sum_i c_i b_i
// Res = - \int_K \sum_i c_i b_i v(x) \cdot \grad b_j dx
// + \int_{\partial K} F(u_h, u_h^{ext}, v(x) \cdot n) b_j ds
// - \int_K \phi b_j
// This is linear in c_i, so we do not need any nonlinear iterations.
// We assemble the jacobian and the right-hand side. The residual is
// equal to Res = Jac*c - rhs, and we compute rhs directly.
const int dim = grid_.dimensions;
const int num_basis = DGBasis::numBasisFunc(dim, degree_);
std::fill(rhs_.begin(), rhs_.end(), 0.0);
std::fill(jac_.begin(), jac_.end(), 0.0);
// Compute cell residual contribution.
// Note: Assumes that \int_K b_j = 0 for all j > 0
rhs_[0] += porevolume_[cell];
// Compute upstream residual contribution.
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
const int face = grid_.cell_faces[hface];
double flux = 0.0;
int upstream_cell = -1;
if (cell == grid_.face_cells[2*face]) {
flux = darcyflux_[face];
upstream_cell = grid_.face_cells[2*face+1];
} else {
flux = -darcyflux_[face];
upstream_cell = grid_.face_cells[2*face];
}
if (upstream_cell < 0) {
// This is an outer boundary. Assumed tof = 0 on inflow, so no contribution.
continue;
}
if (flux >= 0.0) {
// This is an outflow boundary.
continue;
}
// Do quadrature over the face to compute
// \int_{\partial K} u_h^{ext} (v(x) \cdot n) b_j ds
// (where u_h^{ext} is the upstream unknown (tof)).
const double normal_velocity = flux / grid_.face_areas[face];
FaceQuadrature quad(grid_, face, degree_);
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
quad.quadPtCoord(quad_pt, &coord_[0]);
DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
DGBasis::eval(grid_, upstream_cell, degree_, &coord_[0], &basis_nb_[0]);
const double tof_upstream = std::inner_product(basis_nb_.begin(), basis_nb_.end(),
tof_coeff_ + num_basis*upstream_cell, 0.0);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
rhs_[j] -= w * tof_upstream * normal_velocity * basis_[j];
}
}
}
// Compute cell jacobian contribution. We use Fortran ordering
// for jac_, i.e. rows cycling fastest.
{
CellQuadrature quad(grid_, cell, 2*degree_ - 1);
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
// b_i (v \cdot \grad b_j)
quad.quadPtCoord(quad_pt, &coord_[0]);
DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
DGBasis::evalGrad(grid_, cell, degree_, &coord_[0], &grad_basis_[0]);
velocity_interpolation_->interpolate(cell, &coord_[0], &velocity_[0]);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
for (int i = 0; i < num_basis; ++i) {
for (int dd = 0; dd < dim; ++dd) {
jac_[j*num_basis + i] -= w * basis_[j] * grad_basis_[dim*i + dd] * velocity_[dd];
}
}
}
}
}
// Compute downstream jacobian contribution from faces.
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
const int face = grid_.cell_faces[hface];
double flux = 0.0;
if (cell == grid_.face_cells[2*face]) {
flux = darcyflux_[face];
} else {
flux = -darcyflux_[face];
}
if (flux <= 0.0) {
// This is an inflow boundary.
continue;
}
// Do quadrature over the face to compute
// \int_{\partial K} b_i (v(x) \cdot n) b_j ds
const double normal_velocity = flux / grid_.face_areas[face];
FaceQuadrature quad(grid_, face, 2*degree_);
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
// u^ext flux B (B = {b_j})
quad.quadPtCoord(quad_pt, &coord_[0]);
DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
for (int i = 0; i < num_basis; ++i) {
jac_[j*num_basis + i] += w * basis_[i] * normal_velocity * basis_[j];
}
}
}
}
// Compute downstream jacobian contribution from sink terms.
// Contribution from inflow sources would be
// similar to the contribution from upstream faces, but
// it is zero since we let all external inflow be associated
// with a zero tof.
if (source_[cell] < 0.0) {
// A sink.
const double flux = -source_[cell]; // Sign convention for flux: outflux > 0.
const double flux_density = flux / grid_.cell_volumes[cell];
// Do quadrature over the cell to compute
// \int_{K} b_i flux b_j dx
CellQuadrature quad(grid_, cell, 2*degree_);
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
quad.quadPtCoord(quad_pt, &coord_[0]);
DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
for (int i = 0; i < num_basis; ++i) {
jac_[j*num_basis + i] += w * basis_[i] * flux_density * basis_[j];
}
}
}
}
// Solve linear equation.
MAT_SIZE_T n = num_basis;
MAT_SIZE_T nrhs = 1;
MAT_SIZE_T lda = num_basis;
std::vector<MAT_SIZE_T> piv(num_basis);
MAT_SIZE_T ldb = num_basis;
MAT_SIZE_T info = 0;
orig_jac_ = jac_;
dgesv_(&n, &nrhs, &jac_[0], &lda, &piv[0], &rhs_[0], &ldb, &info);
if (info != 0) {
// Print the local matrix and rhs.
std::cerr << "Failed solving single-cell system Ax = b in cell " << cell
<< " with A = \n";
for (int row = 0; row < n; ++row) {
for (int col = 0; col < n; ++col) {
std::cerr << " " << orig_jac_[row + n*col];
}
std::cerr << '\n';
}
std::cerr << "and b = \n";
for (int row = 0; row < n; ++row) {
std::cerr << " " << rhs_[row] << '\n';
}
THROW("Lapack error: " << info << " encountered in cell " << cell);
}
// The solution ends up in rhs_, so we must copy it.
std::copy(rhs_.begin(), rhs_.end(), tof_coeff_ + num_basis*cell);
}
void TransportModelTracerTofDiscGal::solveMultiCell(const int num_cells, const int* cells)
{
std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl;
for (int i = 0; i < num_cells; ++i) {
solveSingleCell(cells[i]);
}
}
} // namespace Opm

105
opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp Normal file
View File

@ -0,0 +1,105 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED
#define OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED
#include <opm/core/transport/reorder/TransportModelInterface.hpp>
#include <boost/shared_ptr.hpp>
#include <vector>
#include <map>
#include <ostream>
struct UnstructuredGrid;
namespace Opm
{
class IncompPropertiesInterface;
class VelocityInterpolationInterface;
/// Implements a discontinuous Galerkin solver for
/// (single-phase) time-of-flight using reordering.
/// The equation solved is:
/// v \cdot \grad\tau = \phi
/// where v is the fluid velocity, \tau is time-of-flight and
/// \phi is the porosity. This is a boundary value problem, where
/// \tau is specified to be zero on all inflow boundaries.
/// The user may specify the polynomial degree of the basis function space
/// used, but only degrees 0 and 1 are supported so far.
class TransportModelTracerTofDiscGal : public TransportModelInterface
{
public:
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] use_cvi If true, use corner point velocity interpolation.
/// Otherwise, use the basic constant interpolation.
TransportModelTracerTofDiscGal(const UnstructuredGrid& grid,
const bool use_cvi);
/// Solve for time-of-flight.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[in] degree Polynomial degree of DG basis functions used.
/// \param[out] tof_coeff Array of time-of-flight solution coefficients.
/// The values are ordered by cell, meaning that
/// the K coefficients corresponding to the first
/// cell comes before the K coefficients corresponding
/// to the second cell etc.
/// K depends on degree and grid dimension.
void solveTof(const double* darcyflux,
const double* porevolume,
const double* source,
const int degree,
std::vector<double>& tof_coeff);
private:
virtual void solveSingleCell(const int cell);
virtual void solveMultiCell(const int num_cells, const int* cells);
private:
// Disable copying and assignment.
TransportModelTracerTofDiscGal(const TransportModelTracerTofDiscGal&);
TransportModelTracerTofDiscGal& operator=(const TransportModelTracerTofDiscGal&);
const UnstructuredGrid& grid_;
boost::shared_ptr<VelocityInterpolationInterface> velocity_interpolation_;
const double* darcyflux_; // one flux per grid face
const double* porevolume_; // one volume per cell
const double* source_; // one volumetric source term per cell
int degree_;
double* tof_coeff_;
std::vector<double> rhs_; // single-cell right-hand-side
std::vector<double> jac_; // single-cell jacobian
std::vector<double> orig_jac_; // single-cell jacobian (copy)
// Below: storage for quantities needed by solveSingleCell().
std::vector<double> coord_;
std::vector<double> basis_;
std::vector<double> basis_nb_;
std::vector<double> grad_basis_;
std::vector<double> velocity_;
};
} // namespace Opm
#endif // OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED

1
opm/core/utility/miscUtilities.cpp
View File

@ -553,6 +553,7 @@ namespace Opm
{
const int np = wells.number_of_phases;
const int nw = wells.number_of_wells;
ASSERT(int(flow_rates_per_well_cell.size()) == wells.well_connpos[nw]);
phase_flow_per_well.resize(nw * np);
for (int wix = 0; wix < nw; ++wix) {
for (int phase = 0; phase < np; ++phase) {

14
opm/core/utility/writeVtkData.hpp
View File

@ -26,25 +26,23 @@
#include <vector>
#include <tr1/array>
#include <iosfwd>
#include <opm/core/utility/DataMap.hpp>
struct UnstructuredGrid;
namespace Opm
{
/// Intended to map strings (giving the output field names) to data.
typedef std::map<std::string, const std::vector<double>*> DataMap;
/// Vtk output for cartesian grids.
void writeVtkData(const std::tr1::array<int, 3>& dims,
const std::tr1::array<double, 3>& cell_size,
const DataMap& data,
std::ostream& os);
const std::tr1::array<double, 3>& cell_size,
const DataMap& data,
std::ostream& os);
/// Vtk output for general grids.
void writeVtkData(const UnstructuredGrid& grid,
const DataMap& data,
std::ostream& os);
const DataMap& data,
std::ostream& os);
} // namespace Opm
#endif // OPM_WRITEVTKDATA_HEADER_INCLUDED

26
opm/core/wells/WellsGroup.cpp
View File

@ -572,6 +572,7 @@ namespace Opm
break;
}
const double total_produced = getTotalProductionFlow(well_surfacerates_phase, phase);
const double total_reinjected = - total_produced; // Production negative, injection positive
const double my_guide_rate = injectionGuideRate(true);
for (size_t i = 0; i < children_.size(); ++i) {
// Apply for all children.
@ -580,11 +581,11 @@ namespace Opm
const double children_guide_rate = children_[i]->injectionGuideRate(true);
#ifdef DIRTY_WELLCTRL_HACK
children_[i]->applyInjGroupControl(InjectionSpecification::RESV,
(children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_,
(children_guide_rate / my_guide_rate) * total_reinjected * injSpec().reinjection_fraction_target_,
false);
#else
children_[i]->applyInjGroupControl(InjectionSpecification::RATE,
(children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_,
(children_guide_rate / my_guide_rate) * total_reinjected * injSpec().reinjection_fraction_target_,
false);
#endif
}
@ -600,15 +601,15 @@ namespace Opm
if (phaseUsage().phase_used[BlackoilPhases::Vapour]) {
total_produced += getTotalProductionFlow(well_reservoirrates_phase, BlackoilPhases::Vapour);
}
const double my_guide_rate = injectionGuideRate(true);
const double total_reinjected = - total_produced; // Production negative, injection positive
const double my_guide_rate = injectionGuideRate(true);
for (size_t i = 0; i < children_.size(); ++i) {
// Apply for all children.
// Note, we do _not_ want to call the applyProdGroupControl in this object,
// as that would check if we're under group control, something we're not.
const double children_guide_rate = children_[i]->injectionGuideRate(true);
children_[i]->applyInjGroupControl(InjectionSpecification::RESV,
(children_guide_rate / my_guide_rate) * total_produced * injSpec().voidage_replacment_fraction_,
(children_guide_rate / my_guide_rate) * total_reinjected * injSpec().voidage_replacment_fraction_,
false);
}
@ -760,16 +761,17 @@ namespace Opm
wells_->ctrls[self_index_]->current = ~ wells_->ctrls[self_index_]->current;
}
}
std::pair<WellNode*, double> WellNode::getWorstOffending(const std::vector<double>& well_reservoirrates_phase,
const std::vector<double>& well_surfacerates_phase,
ProductionSpecification::ControlMode mode)
{
const int np = phaseUsage().num_phases;
const int index = self_index_*np;
return std::make_pair<WellNode*, double>(this, rateByMode(&well_reservoirrates_phase[index],
&well_surfacerates_phase[index],
mode));
return std::pair<WellNode*, double>(this,
rateByMode(&well_reservoirrates_phase[index],
&well_surfacerates_phase[index],
mode));
}
void WellNode::applyInjGroupControl(const InjectionSpecification::ControlMode control_mode,
@ -781,7 +783,7 @@ namespace Opm
&& (injSpec().control_mode_ != InjectionSpecification::GRUP && injSpec().control_mode_ != InjectionSpecification::NONE)) {
return;
}
if (!wells_->type[self_index_] == INJECTOR) {
if (wells_->type[self_index_] != INJECTOR) {
ASSERT(target == 0.0);
return;
}
@ -858,12 +860,12 @@ namespace Opm
std::cout << "Returning" << std::endl;
return;
}
if (!wells_->type[self_index_] == PRODUCER) {
if (wells_->type[self_index_] != PRODUCER) {
ASSERT(target == 0.0);
return;
}
// We're a producer, so we need to negate the input
double ntarget = target;
double ntarget = -target;
double distr[3] = { 0.0, 0.0, 0.0 };
const int* phase_pos = phaseUsage().phase_pos;

8
opm/core/wells/WellsManager.cpp
View File

@ -226,6 +226,14 @@ namespace Opm
/// Construct from existing wells object.
WellsManager::WellsManager(struct Wells* W)
: w_(clone_wells(W))
{
}
/// Construct wells from deck.
WellsManager::WellsManager(const Opm::EclipseGridParser& deck,
const UnstructuredGrid& grid,

9
opm/core/wells/WellsManager.hpp
View File

@ -41,7 +41,14 @@ namespace Opm
{
public:
/// Default constructor -- no wells.
WellsManager();
WellsManager();
/// Construct from existing wells object.
/// WellsManager is not properly initialised in the sense that the logic to
/// manage control switching does not exist.
///
/// @param[in] W Existing wells object.
WellsManager(struct Wells* W);
/// Construct from input deck and grid.
/// The permeability argument may be zero if the input contain

205
tests/test_wells.cpp Normal file
View File

@ -0,0 +1,205 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define NVERBOSE // Suppress own messages when throw()ing
#define BOOST_TEST_MODULE WellsModuleTest
#include <boost/test/unit_test.hpp>
#include <opm/core/newwells.h>
#include <iostream>
#include <vector>
#include <boost/shared_ptr.hpp>
BOOST_AUTO_TEST_CASE(Construction)
{
const int nphases = 2;
const int nwells = 2;
const int nperfs = 2;
boost::shared_ptr<Wells> W(create_wells(nphases, nwells, nperfs),
destroy_wells);
if (W) {
int cells[] = { 0, 9 };
double WI = 1.0;
const double ifrac[] = { 1.0, 0.0 };
const bool ok0 = add_well(INJECTOR, 0.0, 1, &ifrac[0], &cells[0],
&WI, "INJECTOR", W.get());
const double pfrac[] = { 0.0, 0.0 };
const bool ok1 = add_well(PRODUCER, 0.0, 1, &pfrac[0], &cells[1],
&WI, "PRODUCER", W.get());
if (ok0 && ok1) {
BOOST_CHECK_EQUAL(W->number_of_phases, nphases);
BOOST_CHECK_EQUAL(W->number_of_wells , nwells );
BOOST_CHECK_EQUAL(W->well_connpos[0], 0);
BOOST_CHECK_EQUAL(W->well_connpos[1], 1);
BOOST_CHECK_EQUAL(W->well_connpos[W->number_of_wells], nperfs);
BOOST_CHECK_EQUAL(W->well_cells[W->well_connpos[0]], cells[0]);
BOOST_CHECK_EQUAL(W->well_cells[W->well_connpos[1]], cells[1]);
BOOST_CHECK_EQUAL(W->WI[W->well_connpos[0]], WI);
BOOST_CHECK_EQUAL(W->WI[W->well_connpos[1]], WI);
using std::string;
BOOST_CHECK_EQUAL(string(W->name[0]), string("INJECTOR"));
BOOST_CHECK_EQUAL(string(W->name[1]), string("PRODUCER"));
}
}
}
BOOST_AUTO_TEST_CASE(Controls)
{
const int nphases = 2;
const int nwells = 1;
const int nperfs = 2;
boost::shared_ptr<Wells> W(create_wells(nphases, nwells, nperfs),
destroy_wells);
if (W) {
int cells[] = { 0 , 9 };
double WI [] = { 1.0, 1.0 };
const double ifrac[] = { 1.0, 0.0 };
const bool ok = add_well(INJECTOR, 0.0, nperfs, &ifrac[0], &cells[0],
&WI[0], "INJECTOR", W.get());
if (ok) {
const double distr[] = { 1.0, 0.0 };
const bool ok1 = append_well_controls(BHP, 1, &distr[0],
0, W.get());
const bool ok2 = append_well_controls(SURFACE_RATE, 1,
&distr[0], 0, W.get());
if (ok1 && ok2) {
WellControls* ctrls = W->ctrls[0];
BOOST_CHECK_EQUAL(ctrls->num , 2);
BOOST_CHECK_EQUAL(ctrls->current, -1);
set_current_control(0, 0, W.get());
BOOST_CHECK_EQUAL(ctrls->current, 0);
set_current_control(0, 1, W.get());
BOOST_CHECK_EQUAL(ctrls->current, 1);
BOOST_CHECK_EQUAL(ctrls->type[0], BHP);
BOOST_CHECK_EQUAL(ctrls->type[1], SURFACE_RATE);
BOOST_CHECK_EQUAL(ctrls->target[0], 1.0);
BOOST_CHECK_EQUAL(ctrls->target[1], 1.0);
}
}
}
}
BOOST_AUTO_TEST_CASE(Copy)
{
const int nphases = 2;
const int nwells = 2;
const int nperfs = 2;
boost::shared_ptr<Wells> W1(create_wells(nphases, nwells, nperfs),
destroy_wells);
boost::shared_ptr<Wells> W2;
if (W1) {
int cells[] = { 0, 9 };
const double WI = 1.0;
const double ifrac[] = { 1.0, 0.0 };
const bool ok0 = add_well(INJECTOR, 0.0, 1, &ifrac[0], &cells[0],
&WI, "INJECTOR", W1.get());
const double pfrac[] = { 0.0, 0.0 };
const bool ok1 = add_well(PRODUCER, 0.0, 1, &pfrac[0], &cells[1],
&WI, "PRODUCER", W1.get());
bool ok = ok0 && ok1;
for (int w = 0; ok && (w < W1->number_of_wells); ++w) {
const double distr[] = { 1.0, 0.0 };
const bool okc1 = append_well_controls(BHP, 1, &distr[0],
w, W1.get());
const bool okc2 = append_well_controls(SURFACE_RATE, 1,
&distr[0], w,
W1.get());
ok = okc1 && okc2;
}
if (ok) {
W2.reset(clone_wells(W1.get()), destroy_wells);
}
}
if (W2) {
BOOST_CHECK_EQUAL(W2->number_of_phases, W1->number_of_phases);
BOOST_CHECK_EQUAL(W2->number_of_wells , W1->number_of_wells );
BOOST_CHECK_EQUAL(W2->well_connpos[0] , W1->well_connpos[0] );
for (int w = 0; w < W1->number_of_wells; ++w) {
using std::string;
BOOST_CHECK_EQUAL(string(W2->name[w]), string(W1->name[w]));
BOOST_CHECK_EQUAL( W2->type[w] , W1->type[w] );
BOOST_CHECK_EQUAL(W2->well_connpos[w + 1],
W1->well_connpos[w + 1]);
for (int j = W1->well_connpos[w];
j < W1->well_connpos[w + 1]; ++j) {
BOOST_CHECK_EQUAL(W2->well_cells[j], W1->well_cells[j]);
BOOST_CHECK_EQUAL(W2->WI [j], W1->WI [j]);
}
BOOST_CHECK(W1->ctrls[w] != 0);
BOOST_CHECK(W2->ctrls[w] != 0);
WellControls* c1 = W1->ctrls[w];
WellControls* c2 = W2->ctrls[w];
BOOST_CHECK_EQUAL(c2->num , c1->num );
BOOST_CHECK_EQUAL(c2->current, c1->current);
for (int c = 0; c < c1->num; ++c) {
BOOST_CHECK_EQUAL(c2->type [c], c1->type [c]);
BOOST_CHECK_EQUAL(c2->target[c], c1->target[c]);
for (int p = 0; p < W1->number_of_phases; ++p) {
BOOST_CHECK_EQUAL(c2->distr[c*W1->number_of_phases + p],
c1->distr[c*W1->number_of_phases + p]);
}
}
}
}
}