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

Merge branch 'rock_comp_tpfa'

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen committed 2012-08-22 11:13:16 +02:00
2 parents ae66fe74e3 856c33af74
commit 51d4e04d36
10 files changed
+178 -58
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

examples/sim_wateroil.cpp
+15 -29
View File
@@ -141,7 +141,6 @@ main(int argc, char** argv)
std::cout << "--------------- Reading parameters ---------------" << std::endl;
// Reading various control parameters.
const bool use_reorder = param.getDefault("use_reorder", true);
const bool output = param.getDefault("output", true);
std::string output_dir;
int output_interval = 1;
@@ -229,28 +228,8 @@ main(int argc, char** argv)
}
}
bool use_segregation_split = false;
bool use_column_solver = false;
bool use_gauss_seidel_gravity = false;
if (use_gravity && use_reorder) {
if (use_gravity) {
use_segregation_split = param.getDefault("use_segregation_split", use_segregation_split);
if (use_segregation_split) {
use_column_solver = param.getDefault("use_column_solver", use_column_solver);
if (use_column_solver) {
use_gauss_seidel_gravity = param.getDefault("use_gauss_seidel_gravity", use_gauss_seidel_gravity);
}
}
}
// Check that rock compressibility is not used with solvers that do not handle it.
// int nl_pressure_maxiter = 0;
// double nl_pressure_tolerance = 0.0;
if (rock_comp->isActive()) {
THROW("No rock compressibility in comp. pressure solver yet.");
if (!use_reorder) {
THROW("Cannot run implicit (non-reordering) transport solver with rock compressibility yet.");
}
// nl_pressure_maxiter = param.getDefault("nl_pressure_maxiter", 10);
// nl_pressure_tolerance = param.getDefault("nl_pressure_tolerance", 1.0); // in Pascal
}
// Source-related variables init.
@@ -262,11 +241,11 @@ main(int argc, char** argv)
// Extra rock init.
std::vector<double> porevol;
if (rock_comp->isActive()) {
THROW("CompressibleTpfa solver does not handle this.");
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
} else {
computePorevolume(*grid->c_grid(), props->porosity(), porevol);
}
std::vector<double> initial_porevol = porevol;
double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
@@ -293,20 +272,20 @@ main(int argc, char** argv)
const double nl_press_change_tol = param.getDefault("nl_press_change_tol", 10.0);
const int nl_press_maxiter = param.getDefault("nl_press_maxiter", 20);
const double *grav = use_gravity ? &gravity[0] : 0;
Opm::CompressibleTpfa psolver(*grid->c_grid(), *props, linsolver,
Opm::CompressibleTpfa psolver(*grid->c_grid(), *props, rock_comp.get(), linsolver,
nl_press_res_tol, nl_press_change_tol, nl_press_maxiter,
grav, wells->c_wells());
// Reordering solver.
const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9);
const int nl_maxiter = param.getDefault("nl_maxiter", 30);
Opm::TransportModelCompressibleTwophase reorder_model(*grid->c_grid(), *props, nl_tolerance, nl_maxiter);
if (use_gauss_seidel_gravity) {
if (use_segregation_split) {
reorder_model.initGravity();
}
// Column-based gravity segregation solver.
std::vector<std::vector<int> > columns;
if (use_column_solver) {
if (use_segregation_split) {
Opm::extractColumn(*grid->c_grid(), columns);
}
@@ -409,6 +388,11 @@ main(int argc, char** argv)
Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0,
wells->c_wells(), well_state.perfRates(), reorder_src);
// Compute new porevolumes after pressure solve, if necessary.
if (rock_comp->isActive()) {
initial_porevol = porevol;
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
}
// Solve transport.
transport_timer.start();
double stepsize = simtimer.currentStepLength();
@@ -419,11 +403,13 @@ main(int argc, char** argv)
for (int tr_substep = 0; tr_substep < num_transport_substeps; ++tr_substep) {
// Note that for now we do not handle rock compressibility,
// although the transport solver should be able to.
reorder_model.solve(&state.faceflux()[0], &state.pressure()[0], &state.surfacevol()[0],
&porevol[0], &porevol[0], &reorder_src[0], stepsize, state.saturation());
reorder_model.solve(&state.faceflux()[0], &state.pressure()[0],
&porevol[0], &initial_porevol[0], &reorder_src[0], stepsize,
state.saturation(), state.surfacevol());
// Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced);
if (use_segregation_split) {
reorder_model.solveGravity(columns, &state.pressure()[0], &porevol[0], stepsize, grav, state.saturation());
reorder_model.solveGravity(columns, &state.pressure()[0], &initial_porevol[0],
stepsize, grav, state.saturation());
}
}
transport_timer.stop();
opm/core/fluid/RockFromDeck.cpp
+1 -1
View File
@@ -87,7 +87,7 @@ namespace Opm
double perm_threshold)
{
const int dim = 3;
const int num_global_cells = numGlobalCells(parser);
const int num_global_cells = grid.cartdims[0]*grid.cartdims[1]*grid.cartdims[2];
const int nc = grid.number_of_cells;
ASSERT (num_global_cells > 0);
opm/core/pressure/CompressibleTpfa.cpp
+34 -5
View File
@@ -30,6 +30,7 @@
#include <opm/core/newwells.h>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/fluid/RockCompressibility.hpp>
#include <algorithm>
#include <cmath>
@@ -58,6 +59,7 @@ namespace Opm
/// to change.
CompressibleTpfa::CompressibleTpfa(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
const LinearSolverInterface& linsolver,
const double residual_tol,
const double change_tol,
@@ -66,6 +68,7 @@ namespace Opm
const struct Wells* wells)
: grid_(grid),
props_(props),
rock_comp_props_(rock_comp_props),
linsolver_(linsolver),
residual_tol_(residual_tol),
change_tol_(change_tol),
@@ -74,7 +77,6 @@ namespace Opm
wells_(wells),
htrans_(grid.cell_facepos[ grid.number_of_cells ]),
trans_ (grid.number_of_faces),
porevol_(grid.number_of_cells),
allcells_(grid.number_of_cells)
{
if (wells_ && (wells_->number_of_phases != props.numPhases())) {
@@ -86,7 +88,12 @@ namespace Opm
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
tpfa_htrans_compute(gg, props.permeability(), &htrans_[0]);
tpfa_trans_compute(gg, &htrans_[0], &trans_[0]);
computePorevolume(grid_, props.porosity(), porevol_);
// If we have rock compressibility, pore volumes are updated
// in the compute*() methods, otherwise they are constant and
// hence may be computed here.
if (rock_comp_props_ == NULL || !rock_comp_props_->isActive()) {
computePorevolume(grid_, props.porosity(), porevol_);
}
for (int c = 0; c < grid.number_of_cells; ++c) {
allcells_[c] = c;
}
@@ -230,6 +237,9 @@ namespace Opm
const WellState& /*well_state*/)
{
computeWellPotentials(state);
if (rock_comp_props_ && rock_comp_props_->isActive()) {
computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), initial_porevol_);
}
}
@@ -252,6 +262,8 @@ namespace Opm
// std::vector<double> face_gravcap_;
// std::vector<double> wellperf_A_;
// std::vector<double> wellperf_phasemob_;
// std::vector<double> porevol_; // Only modified if rock_comp_props_ is non-null.
// std::vector<double> rock_comp_; // Empty unless rock_comp_props_ is non-null.
computeCellDynamicData(dt, state, well_state);
computeFaceDynamicData(dt, state, well_state);
computeWellDynamicData(dt, state, well_state);
@@ -273,6 +285,8 @@ namespace Opm
// std::vector<double> cell_viscosity_;
// std::vector<double> cell_phasemob_;
// std::vector<double> cell_voldisc_;
// std::vector<double> porevol_; // Only modified if rock_comp_props_ is non-null.
// std::vector<double> rock_comp_; // Empty unless rock_comp_props_ is non-null.
const int nc = grid_.number_of_cells;
const int np = props_.numPhases();
const double* cell_p = &state.pressure()[0];
@@ -296,6 +310,14 @@ namespace Opm
// TODO: Check this!
cell_voldisc_.clear();
cell_voldisc_.resize(nc, 0.0);
if (rock_comp_props_ && rock_comp_props_->isActive()) {
computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol_);
rock_comp_.resize(nc);
for (int cell = 0; cell < nc; ++cell) {
rock_comp_[cell] = rock_comp_props_->rockComp(state.pressure()[cell]);
}
}
}
@@ -465,9 +487,16 @@ namespace Opm
cq.Af = &face_A_[0];
cq.phasemobf = &face_phasemob_[0];
cq.voldiscr = &cell_voldisc_[0];
cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0],
&face_gravcap_[0], cell_press, well_bhp,
&porevol_[0], h_);
if (rock_comp_props_ == NULL || !rock_comp_props_->isActive()) {
cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0],
&face_gravcap_[0], cell_press, well_bhp,
&porevol_[0], h_);
} else {
cfs_tpfa_res_comprock_assemble(gg, dt, &forces, z, &cq, &trans_[0],
&face_gravcap_[0], cell_press, well_bhp,
&porevol_[0], &initial_porevol_[0],
&rock_comp_[0], h_);
}
}
opm/core/pressure/CompressibleTpfa.hpp
+22 -16
View File
@@ -33,6 +33,7 @@ namespace Opm
class BlackoilState;
class BlackoilPropertiesInterface;
class RockCompressibility;
class LinearSolverInterface;
class WellState;
@@ -44,23 +45,25 @@ namespace Opm
{
public:
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] props Rock and fluid properties.
/// \param[in] linsolver Linear solver to use.
/// \param[in] residual_tol Solution accepted if inf-norm of residual is smaller.
/// \param[in] change_tol Solution accepted if inf-norm of change in pressure is smaller.
/// \param[in] maxiter Maximum acceptable number of iterations.
/// \param[in] gravity Gravity vector. If non-null, the array should
/// have D elements.
/// \param[in] wells The wells argument. Will be used in solution,
/// is ignored if NULL.
/// Note: this class observes the well object, and
/// makes the assumption that the well topology
/// and completions does not change during the
/// run. However, controls (only) are allowed
/// to change.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] props Rock and fluid properties.
/// \param[in] rock_comp_props Rock compressibility properties. May be null.
/// \param[in] linsolver Linear solver to use.
/// \param[in] residual_tol Solution accepted if inf-norm of residual is smaller.
/// \param[in] change_tol Solution accepted if inf-norm of change in pressure is smaller.
/// \param[in] maxiter Maximum acceptable number of iterations.
/// \param[in] gravity Gravity vector. If non-null, the array should
/// have D elements.
/// \param[in] wells The wells argument. Will be used in solution,
/// is ignored if NULL.
/// Note: this class observes the well object, and
/// makes the assumption that the well topology
/// and completions does not change during the
/// run. However, controls (only) are allowed
/// to change.
CompressibleTpfa(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
const LinearSolverInterface& linsolver,
const double residual_tol,
const double change_tol,
@@ -107,6 +110,7 @@ namespace Opm
// ------ Data that will remain unmodified after construction. ------
const UnstructuredGrid& grid_;
const BlackoilPropertiesInterface& props_;
const RockCompressibility* rock_comp_props_;
const LinearSolverInterface& linsolver_;
const double residual_tol_;
const double change_tol_;
@@ -115,7 +119,6 @@ namespace Opm
const Wells* wells_; // May be NULL, outside may modify controls (only) between calls to solve().
std::vector<double> htrans_;
std::vector<double> trans_ ;
std::vector<double> porevol_;
std::vector<int> allcells_;
// ------ Internal data for the cfs_tpfa_res solver. ------
@@ -123,6 +126,7 @@ namespace Opm
// ------ Data that will be modified for every solve. ------
std::vector<double> wellperf_gpot_;
std::vector<double> initial_porevol_;
// ------ Data that will be modified for every solver iteration. ------
std::vector<double> cell_A_;
@@ -135,6 +139,8 @@ namespace Opm
std::vector<double> face_gravcap_;
std::vector<double> wellperf_A_;
std::vector<double> wellperf_phasemob_;
std::vector<double> porevol_; // Only modified if rock_comp_props_ is non-null.
std::vector<double> rock_comp_; // Empty unless rock_comp_props_ is non-null.
// The update to be applied to the pressures (cell and bhp).
std::vector<double> pressure_increment_;
opm/core/pressure/tpfa/cfs_tpfa_residual.c
+73
View File
@@ -1213,6 +1213,79 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G ,
}
/* ---------------------------------------------------------------------- */
void
cfs_tpfa_res_comprock_assemble(
struct UnstructuredGrid *G ,
double dt ,
struct cfs_tpfa_res_forces *forces ,
const double *zc ,
struct compr_quantities_gen *cq ,
const double *trans ,
const double *gravcap_f,
const double *cpress ,
const double *wpress ,
const double *porevol ,
const double *porevol0 ,
const double *rock_comp,
struct cfs_tpfa_res_data *h )
/* ---------------------------------------------------------------------- */
{
/* We want to add this term to the usual residual:
*
* (porevol(pressure)-porevol(initial_pressure))/dt.
*
* Its derivative (for the diagonal term of the Jacobian) is:
*
* porevol(pressure)*rock_comp(pressure)/dt
*/
int c, w, well_is_neumann, rock_is_incomp;
size_t j;
double dpv;
const struct Wells* W;
/* Assemble usual system (without rock compressibility). */
cfs_tpfa_res_assemble(G, dt, forces, zc, cq, trans, gravcap_f,
cpress, wpress, porevol0, h);
/* Check if we have only Neumann wells. */
well_is_neumann = 1;
W = forces->wells->W;
for (w = 0; well_is_neumann && w < W->number_of_wells; w++) {
if ((W->ctrls[w]->current >= 0) && /* OPEN? */
(W->ctrls[w]->type[ W->ctrls[w]->current ] == BHP)) {
well_is_neumann = 0;
}
}
/* If we made a singularity-removing adjustment in the
regular assembly, we undo it here. */
if (well_is_neumann && h->pimpl->is_incomp) {
h->J->sa[0] /= 2;
}
/* Add new terms to residual and Jacobian. */
rock_is_incomp = 1;
for (c = 0; c < G->number_of_cells; c++) {
j = csrmatrix_elm_index(c, c, h->J);
dpv = (porevol[c] - porevol0[c]);
if (dpv != 0.0 || rock_comp[c] != 0.0) {
rock_is_incomp = 0;
}
h->J->sa[j] += porevol[c] * rock_comp[c];
h->F[c] += dpv;
}
/* Re-do the singularity-removing adjustment if necessary */
if (rock_is_incomp && well_is_neumann && h->pimpl->is_incomp) {
h->J->sa[0] *= 2;
}
}
/* ---------------------------------------------------------------------- */
void
cfs_tpfa_res_flux(struct UnstructuredGrid *G ,
opm/core/pressure/tpfa/cfs_tpfa_residual.h
+16
View File
@@ -72,6 +72,22 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G,
const double *porevol,
struct cfs_tpfa_res_data *h);
void
cfs_tpfa_res_comprock_assemble(
struct UnstructuredGrid *G,
double dt,
struct cfs_tpfa_res_forces *forces,
const double *zc,
struct compr_quantities_gen *cq,
const double *trans,
const double *gravcap_f,
const double *cpress,
const double *wpress,
const double *porevol,
const double *porevol0,
const double *rock_comp,
struct cfs_tpfa_res_data *h);
void
cfs_tpfa_res_flux(struct UnstructuredGrid *G ,
struct cfs_tpfa_res_forces *forces ,
opm/core/pressure/tpfa/ifs_tpfa.c
+1 -1
View File
@@ -771,7 +771,7 @@ ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid *G ,
assemble_incompressible(G, F, trans, gpress, h, &system_singular, &ok);
/* We want to solve a Newton step for the residual
* (porevol(pressure)-porevol(initial_pressure))/dt + residual_for_imcompressible
* (porevol(pressure)-porevol(initial_pressure))/dt + residual_for_incompressible
*
*/
opm/core/simulator/SimulatorIncompTwophase.cpp
+1 -1
View File
@@ -150,7 +150,7 @@ namespace Opm
catch (...) {
THROW("Creating directories failed: " << fpath);
}
vtkfilename << "/" << std::setw(3) << std::setfill('0') << step << ".vtu";
vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
std::ofstream vtkfile(vtkfilename.str().c_str());
if (!vtkfile) {
THROW("Failed to open " << vtkfilename.str());
opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp
+13 -3
View File
@@ -75,15 +75,15 @@ namespace Opm
void TransportModelCompressibleTwophase::solve(const double* darcyflux,
const double* pressure,
const double* surfacevol0,
const double* porevolume0,
const double* porevolume,
const double* source,
const double dt,
std::vector<double>& saturation)
std::vector<double>& saturation,
std::vector<double>& surfacevol)
{
darcyflux_ = darcyflux;
surfacevol0_ = surfacevol0;
surfacevol0_ = &surfacevol[0];
porevolume0_ = porevolume0;
porevolume_ = porevolume;
source_ = source;
@@ -107,6 +107,15 @@ namespace Opm
&ia_downw_[0], &ja_downw_[0]);
reorderAndTransport(grid_, darcyflux);
toBothSat(saturation_, saturation);
// Compute surface volume as a postprocessing step from saturation and A_
surfacevol = saturation;
const int np = props_.numPhases();
for (int cell = 0; cell < grid_.number_of_cells; ++cell) {
for (int phase = 0; phase < np; ++phase) {
surfacevol[np*cell + phase] *= A_[np*np*cell + np*phase + phase];
}
}
}
// Residual function r(s) for a single-cell implicit Euler transport
@@ -381,6 +390,7 @@ namespace Opm
std::vector<double> htrans(grid_.cell_facepos[grid_.number_of_cells]);
const int nf = grid_.number_of_faces;
trans_.resize(nf);
gravflux_.resize(nf);
tpfa_htrans_compute(const_cast<UnstructuredGrid*>(&grid_), props_.permeability(), &htrans[0]);
tpfa_trans_compute(const_cast<UnstructuredGrid*>(&grid_), &htrans[0], &trans_[0]);
}
opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp
+2 -2
View File
@@ -54,12 +54,12 @@ namespace Opm
/// \param[in, out] saturation Phase saturations.
void solve(const double* darcyflux,
const double* pressure,
const double* surfacevol0,
const double* porevolume0,
const double* porevolume,
const double* source,
const double dt,
std::vector<double>& saturation);
std::vector<double>& saturation,
std::vector<double>& surfacevol);
/// Initialise quantities needed by gravity solver.
void initGravity();