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 pull request #847 from totto82/frankenstein_fixModel2

Browse Source 
Frankenstein fixes for model2
This commit is contained in:
Andreas Lauser
2016-09-30 11:16:03 +02:00
committed by GitHub
parent 7d7e619333 8e2657ce4d
commit f9bf09e832
5 changed files with 62 additions and 25 deletions
Download Patch File Download Diff File Expand all files Collapse all files

51
opm/autodiff/BlackoilModelEbos.hpp
View File

@@ -456,7 +456,10 @@ namespace Opm {
for (int cell_idx = 0; cell_idx < nc; ++cell_idx) {
double dp = dx[cell_idx][flowPhaseToEbosCompIdx(0)];
reservoir_state.pressure()[cell_idx] -= dp;
//reservoir_state.pressure()[cell_idx] -= dp;
double& p = reservoir_state.pressure()[cell_idx];
p -= dp;
p = std::max(p, 1e5);
// Saturation updates.
const double dsw = active_[Water] ? dx[cell_idx][flowPhaseToEbosCompIdx(1)] : 0.0;
@@ -512,29 +515,38 @@ namespace Opm {
if (has_disgas_) {
double& rs = reservoir_state.gasoilratio()[cell_idx];
rs -= drs;
rs = std::max(rs, 0.0);
}
if (has_vapoil_) {
double& rv = reservoir_state.rv()[cell_idx];
rv -= drv;
rv = std::max(rv, 0.0);
}
// Sg is used as primal variable for water only cells.
const double epsilon = 1e-4; //std::sqrt(std::numeric_limits<double>::epsilon());
double& sw = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Water ]];
double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]];
double& rs = reservoir_state.gasoilratio()[cell_idx];
double& rv = reservoir_state.rv()[cell_idx];
// phase translation sg <-> rs
const HydroCarbonState hydroCarbonState = reservoir_state.hydroCarbonState()[cell_idx];
const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
const auto& fs = intQuants.fluidState();
switch (hydroCarbonState) {
case HydroCarbonState::GasAndOil: {
double& sw = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Water ]];
double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]];
if (sw > (1.0 - epsilon)) // water only i.e. do nothing
break;
if (sg <= 0.0 && has_disgas_) {
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::OilOnly; // sg --> rs
sg = 0;
so = 1.0 - sw - sg;
double rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(0, reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
double rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
double& rs = reservoir_state.gasoilratio()[cell_idx];
rs = rsSat*(1-epsilon);
} else if (so <= 0.0 && has_vapoil_) {
@@ -543,16 +555,25 @@ namespace Opm {
sg = 1.0 - sw - so;
double& rv = reservoir_state.rv()[cell_idx];
// use gas pressure?
double rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(0, reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
double rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
rv = rvSat*(1-epsilon);
}
break;
}
case HydroCarbonState::OilOnly: {
double& rs = reservoir_state.gasoilratio()[cell_idx];
double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]];
// TODO:: not hardcode pvtRegion = 0
double rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(0, reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
if (sw > (1.0 - epsilon)) {
// water only change to Sg
rs = 0;
rv = 0;
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
//std::cout << "watonly rv -> sg" << cell_idx << std::endl;
break;
}
double rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
if (rs > ( rsSat * (1+epsilon) ) ) {
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
sg = epsilon;
@@ -562,13 +583,19 @@ namespace Opm {
break;
}
case HydroCarbonState::GasOnly: {
double& rv = reservoir_state.rv()[cell_idx];
double rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(0, reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
if (sw > (1.0 - epsilon)) {
// water only change to Sg
rs = 0;
rv = 0;
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
//std::cout << "watonly rv -> sg" << cell_idx << std::endl;
break;
}
double rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
if (rv > rvSat * (1+epsilon) ) {
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
so = epsilon;
rv = rvSat;
double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]];
sg -= epsilon;
}
break;

2
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@@ -239,7 +239,7 @@ public:
// write the inital state at the report stage
if (timer.initialStep()) {
output_writer_.writeTimeStep( timer, state, well_state, solver->model() );
//output_writer_.writeTimeStep( timer, state, well_state, solver->model() );
}
if( terminal_output_ )

2
opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp
View File

@@ -167,7 +167,7 @@ namespace Opm
substep,
timer.simulationTimeElapsed(),
simToSolution( state, phaseUsage_ ),
wellState.report(),
wellState.report(phaseUsage_),
simProps);
}
}

28
opm/autodiff/StandardWellsDense.hpp
View File

@@ -226,6 +226,7 @@ namespace Opm {
const double volume = 0.002831684659200; // 0.1 cu ft;
for (int w = 0; w < nw; ++w) {
double total_flux = 0.0;
for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) {
@@ -243,6 +244,7 @@ namespace Opm {
// subtract sum of phase fluxes in the well equations.
resWell_[w][flowPhaseToEbosCompIdx(p1)] -= cq_s[p1].value;
total_flux += cq_s[p1].value;
// assemble the jacobians
for (int p2 = 0; p2 < np; ++p2) {
@@ -260,15 +262,23 @@ namespace Opm {
// Store the perforation pressure for later usage.
well_state.perfPress()[perf] = well_state.bhp()[w] + wellPerforationPressureDiffs()[perf];
}
// add vol * dF/dt + Q to the well equations;
for (int p1 = 0; p1 < np; ++p1) {
EvalWell resWell_loc = (wellVolumeFraction(w, p1) - F0_[w + nw*p1]) * volume / dt;
resWell_loc += getQs(w, p1);
for (int p2 = 0; p2 < np; ++p2) {
invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] += resWell_loc.derivatives[p2+3];
if (std::abs(total_flux) == 0) { // add a trivial equation for the case when there is no flow
for (int p1 = 0; p1 < np; ++p1) {
invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p1)] += 1;
resWell_[w][flowPhaseToEbosCompIdx(p1)] += 0;
}
}
else {
// add vol * dF/dt + Q to the well equations;
for (int p1 = 0; p1 < np; ++p1) {
EvalWell resWell_loc = (wellVolumeFraction(w, p1) - F0_[w + nw*p1]) * volume / dt;
resWell_loc += getQs(w, p1);
for (int p2 = 0; p2 < np; ++p2) {
invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] += resWell_loc.derivatives[p2+3];
}
resWell_[w][flowPhaseToEbosCompIdx(p1)] += resWell_loc.value;
}
resWell_[w][flowPhaseToEbosCompIdx(p1)] += resWell_loc.value;
}
}
@@ -997,7 +1007,7 @@ namespace Opm {
{
const int sign1 = dwells[w][flowPhaseToEbosCompIdx(0)] > 0 ? 1: -1;
const double dx1_limited = sign1 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(0)]),std::abs(xvar_well_old[w])*dBHPLimit);
well_state.wellSolutions()[w] = xvar_well_old[w] - dx1_limited;
well_state.wellSolutions()[w] = std::max(xvar_well_old[w] - dx1_limited,1e5);
well_state.bhp()[w] = well_state.wellSolutions()[w];
if (well_controls_iget_type(wc, current) == SURFACE_RATE) {

4
opm/autodiff/WellStateFullyImplicitBlackoil.hpp
View File

@@ -276,8 +276,8 @@ namespace Opm
std::vector<double>& wellSolutions() { return well_solutions_; }
const std::vector<double>& wellSolutions() const { return well_solutions_; }
data::Wells report() const override {
data::Wells res = WellState::report();
data::Wells report(const PhaseUsage& pu) const override {
data::Wells res = WellState::report(pu);
const int nw = this->numWells();
// If there are now wells numPhases throws a floating point