Update surface volume in transport solver.

This commit is contained in:
Xavier Raynaud
2012-08-22 10:11:32 +02:00
parent a8eefa4a4c
commit 6fd14bd76a
3 changed files with 17 additions and 7 deletions

View File

@@ -403,8 +403,9 @@ 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], &initial_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], &initial_porevol[0],

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

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();