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 #144 from totto82/primalvariable

Browse Source 
Refactor primary variable switching
This commit is contained in:
Alf Birger Rustad 2014-06-13 16:18:51 +02:00
commit 2ffaeae466
8 changed files with 124 additions and 80 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
examples/sim_fibo_ad.cpp
View File

@ -218,7 +218,10 @@ try
rock_comp->isActive() ? rock_comp.get() : 0, rock_comp->isActive() ? rock_comp.get() : 0,
wells, wells,
*fis_solver, *fis_solver,
grav); grav,
deck->hasKeyword("DISGAS"),
deck->hasKeyword("VAPOIL") );
SimulatorReport episodeReport = simulator.run(simtimer, state, well_state); SimulatorReport episodeReport = simulator.run(simtimer, state, well_state);
++simtimer; ++simtimer;

2
examples/test_implicit_ad.cpp
View File

@ -104,7 +104,7 @@ try
Opm::NewtonIterationBlackoilSimple fis_solver(param); Opm::NewtonIterationBlackoilSimple fis_solver(param);
Opm::FullyImplicitBlackoilSolver<UnstructuredGrid> solver(param, *g, props, geo, 0, *wells, fis_solver); Opm::FullyImplicitBlackoilSolver<UnstructuredGrid> solver(param, *g, props, geo, 0, *wells, fis_solver, /*hasDisgas*/ true, /*hasVapoil=*/false);
Opm::BlackoilState state; Opm::BlackoilState state;
initStateBasic(*g, props0, param, 0.0, state); initStateBasic(*g, props0, param, 0.0, state);

2
opm/autodiff/BlackoilPropsAdFromDeck.cpp
View File

@ -94,6 +94,8 @@ namespace Opm
phase_usage_ = phaseUsageFromDeck(deck); phase_usage_ = phaseUsageFromDeck(deck);
// Surface densities. Accounting for different orders in eclipse and our code. // Surface densities. Accounting for different orders in eclipse and our code.
Opm::DeckKeywordConstPtr densityKeyword = deck->getKeyword("DENSITY"); Opm::DeckKeywordConstPtr densityKeyword = deck->getKeyword("DENSITY");
int numRegions = densityKeyword->size(); int numRegions = densityKeyword->size();

2
opm/autodiff/BlackoilPropsAdFromDeck.hpp
View File

@ -360,6 +360,8 @@ namespace Opm
std::unique_ptr<SaturationPropsInterface> satprops_; std::unique_ptr<SaturationPropsInterface> satprops_;
PhaseUsage phase_usage_; PhaseUsage phase_usage_;
bool has_vapoil_;
bool has_disgas_;
// The PVT region which is to be used for each cell // The PVT region which is to be used for each cell
std::vector<int> cellPvtRegionIdx_; std::vector<int> cellPvtRegionIdx_;

16
opm/autodiff/FullyImplicitBlackoilSolver.hpp
View File

@ -70,7 +70,9 @@ namespace Opm {
const DerivedGeology& geo , const DerivedGeology& geo ,
const RockCompressibility* rock_comp_props, const RockCompressibility* rock_comp_props,
const Wells& wells, const Wells& wells,
const NewtonIterationBlackoilInterface& linsolver); const NewtonIterationBlackoilInterface& linsolver,
const bool has_disgas,
const bool has_vapoil );
/// Take a single forward step, modifiying /// Take a single forward step, modifiying
/// state.pressure() /// state.pressure()
@ -127,7 +129,7 @@ namespace Opm {
// the Newton relaxation type // the Newton relaxation type
enum RelaxType { DAMPEN, SOR }; enum RelaxType { DAMPEN, SOR };
enum PrimalVariables { Sg = 0, RS = 1, RV = 2 };
// Member data // Member data
const Grid& grid_; const Grid& grid_;
@ -144,6 +146,8 @@ namespace Opm {
HelperOps ops_; HelperOps ops_;
const WellOps wops_; const WellOps wops_;
const M grav_; const M grav_;
const bool has_disgas_;
const bool has_vapoil_;
double dp_max_rel_; double dp_max_rel_;
double ds_max_; double ds_max_;
double drs_max_rel_; double drs_max_rel_;
@ -159,6 +163,8 @@ namespace Opm {
LinearisedBlackoilResidual residual_; LinearisedBlackoilResidual residual_;
std::vector<int> primalVariable_;
// Private methods. // Private methods.
SolutionState SolutionState
constantState(const BlackoilState& x, constantState(const BlackoilState& x,
@ -279,6 +285,12 @@ namespace Opm {
void void
classifyCondition(const BlackoilState& state); classifyCondition(const BlackoilState& state);
/// update the primal variable for Sg, Rv or Rs. The Gas phase must
/// be active to call this method.
void
updatePrimalVariableFromState(const BlackoilState& state);
/// Compute convergence based on total mass balance (tol_mb) and maximum /// Compute convergence based on total mass balance (tol_mb) and maximum
/// residual mass balance (tol_cnv). /// residual mass balance (tol_cnv).
bool getConvergence(const double dt); bool getConvergence(const double dt);

143
opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
View File

@ -219,7 +219,9 @@ namespace {
const DerivedGeology& geo , const DerivedGeology& geo ,
const RockCompressibility* rock_comp_props, const RockCompressibility* rock_comp_props,
const Wells& wells, const Wells& wells,
const NewtonIterationBlackoilInterface& linsolver) const NewtonIterationBlackoilInterface& linsolver,
const bool has_disgas,
const bool has_vapoil)
: grid_ (grid) : grid_ (grid)
, fluid_ (fluid) , fluid_ (fluid)
, geo_ (geo) , geo_ (geo)
@ -232,6 +234,8 @@ namespace {
, ops_ (grid) , ops_ (grid)
, wops_ (wells) , wops_ (wells)
, grav_ (gravityOperator(grid_, ops_, geo_)) , grav_ (gravityOperator(grid_, ops_, geo_))
, has_disgas_(has_disgas)
, has_vapoil_(has_vapoil)
, dp_max_rel_ (1.0e9) , dp_max_rel_ (1.0e9)
, ds_max_ (0.2) , ds_max_ (0.2)
, drs_max_rel_ (1.0e9) , drs_max_rel_ (1.0e9)
@ -272,7 +276,8 @@ namespace {
{ {
const V pvdt = geo_.poreVolume() / dt; const V pvdt = geo_.poreVolume() / dt;
classifyCondition(x); if (active_[Gas]) { updatePrimalVariableFromState(x); }
{ {
const SolutionState state = constantState(x, xw); const SolutionState state = constantState(x, xw);
computeAccum(state, 0); computeAccum(state, 0);
@ -463,26 +468,13 @@ namespace {
V isRs = V::Zero(nc,1); V isRs = V::Zero(nc,1);
V isRv = V::Zero(nc,1); V isRv = V::Zero(nc,1);
V isSg = V::Zero(nc,1); V isSg = V::Zero(nc,1);
bool disgas = false;
bool vapoil = false;
if (active_[ Gas ]){ if (active_[ Gas ]){
// this is a temporary hack to find if vapoil or disgas
// is a active component. Should be given directly from
// DISGAS and VAPOIL keywords in the deck.
for (int c = 0; c < nc ; c++ ) { for (int c = 0; c < nc ; c++ ) {
if(x.rv()[c] > 0) if ( primalVariable_[c] == PrimalVariables::RS ) {
vapoil = true;
if(x.gasoilratio ()[c] > 0)
disgas = true;
}
for (int c = 0; c < nc ; c++ ) {
const PhasePresence cond = phaseCondition()[c];
if ( (!cond.hasFreeGas()) && disgas ) {
isRs[c] = 1; isRs[c] = 1;
} }
else if ( (!cond.hasFreeOil()) && vapoil ) { else if ( primalVariable_[c] == PrimalVariables::RV ) {
isRv[c] = 1; isRv[c] = 1;
} }
else { else {
@ -546,12 +538,12 @@ namespace {
state.saturation[ pu.phase_pos[ Gas ] ] = sg; state.saturation[ pu.phase_pos[ Gas ] ] = sg;
so = so - sg; so = so - sg;
if (disgas) { if (has_disgas_) {
state.rs = (1-isRs) * rsSat + isRs*xvar; state.rs = (1-isRs) * rsSat + isRs*xvar;
} else { } else {
state.rs = rsSat; state.rs = rsSat;
} }
if (vapoil) { if (has_vapoil_) {
state.rv = (1-isRv) * rvSat + isRv*xvar; state.rv = (1-isRv) * rvSat + isRv*xvar;
} else { } else {
state.rv = rvSat; state.rv = rvSat;
@ -1232,33 +1224,19 @@ namespace {
V isRs = V::Zero(nc,1); V isRs = V::Zero(nc,1);
V isRv = V::Zero(nc,1); V isRv = V::Zero(nc,1);
V isSg = V::Zero(nc,1); V isSg = V::Zero(nc,1);
if (active_[Gas]) {
bool disgas = false;
bool vapoil = false;
// this is a temporary hack to find if vapoil or disgas
// is a active component. Should be given directly from
// DISGAS and VAPOIL keywords in the deck.
for (int c = 0; c < nc ; c++ ) { for (int c = 0; c < nc ; c++ ) {
if(state.rv()[c]>0) if ( primalVariable_[c] == PrimalVariables::RS ) {
vapoil = true;
if(state.gasoilratio()[c]>0)
disgas = true;
}
const std::vector<PhasePresence> conditions = phaseCondition();
for (int c = 0; c < nc; c++ ) {
const PhasePresence cond = conditions[c];
if ( (!cond.hasFreeGas()) && disgas ) {
isRs[c] = 1; isRs[c] = 1;
} }
else if ( (!cond.hasFreeOil()) && vapoil ) { else if ( primalVariable_[c] == PrimalVariables::RV ) {
isRv[c] = 1; isRv[c] = 1;
} }
else { else {
isSg[c] = 1; isSg[c] = 1;
} }
} }
}
// Extract parts of dx corresponding to each part. // Extract parts of dx corresponding to each part.
const V dp = subset(dx, Span(nc)); const V dp = subset(dx, Span(nc));
@ -1330,14 +1308,14 @@ namespace {
const double drsmaxrel = drsMaxRel(); const double drsmaxrel = drsMaxRel();
const double drvmax = 1e9;//% same as in Mrst const double drvmax = 1e9;//% same as in Mrst
V rs; V rs;
if (disgas) { if (has_disgas_) {
const V rs_old = Eigen::Map<const V>(&state.gasoilratio()[0], nc); const V rs_old = Eigen::Map<const V>(&state.gasoilratio()[0], nc);
const V drs = isRs * dxvar; const V drs = isRs * dxvar;
const V drs_limited = sign(drs) * drs.abs().min(rs_old.abs()*drsmaxrel); const V drs_limited = sign(drs) * drs.abs().min(rs_old.abs()*drsmaxrel);
rs = rs_old - drs_limited; rs = rs_old - drs_limited;
} }
V rv; V rv;
if (vapoil) { if (has_vapoil_) {
const V rv_old = Eigen::Map<const V>(&state.rv()[0], nc); const V rv_old = Eigen::Map<const V>(&state.rv()[0], nc);
const V drv = isRv * dxvar; const V drv = isRv * dxvar;
const V drv_limited = sign(drv) * drv.abs().min(drvmax); const V drv_limited = sign(drv) * drv.abs().min(drvmax);
@ -1356,23 +1334,23 @@ namespace {
// reset the phase conditions // reset the phase conditions
std::vector<PhasePresence> cond(nc); std::vector<PhasePresence> cond(nc);
if (disgas) { std::fill(primalVariable_.begin(), primalVariable_.end(), PrimalVariables::Sg);
// The obvioious case
auto ix0 = (sg > 0 && isRs == 0); if (has_disgas_) {
// The obvious case
auto hasGas = (sg > 0 && isRs == 0);
// keep oil saturated if previous sg is sufficient large: // keep oil saturated if previous sg is sufficient large:
const int pos = pu.phase_pos[ Gas ]; const int pos = pu.phase_pos[ Gas ];
auto ix1 = (sg < 0 && s_old.col(pos) > epsilon); auto hadGas = (sg < 0 && s_old.col(pos) > epsilon);
// Set oil saturated if previous rs is sufficiently large // Set oil saturated if previous rs is sufficiently large
const V rs_old = Eigen::Map<const V>(&state.gasoilratio()[0], nc); const V rs_old = Eigen::Map<const V>(&state.gasoilratio()[0], nc);
auto ix2 = ( (rs > rsSat * (1+epsilon) && isRs == 1 ) && (rs_old > rsSat0 * (1-epsilon)) ); auto gasVaporized = ( (rs > rsSat * (1+epsilon) && isRs == 1 ) && (rs_old > rsSat0 * (1-epsilon)) );
auto gasPresent = watOnly || ix0 || ix1 || ix2; auto useSg = watOnly || hasGas || hadGas || gasVaporized;
for (int c = 0; c < nc; ++c) { for (int c = 0; c < nc; ++c) {
if (gasPresent[c]) { if (useSg[c]) { rs[c] = rsSat[c];}
rs[c] = rsSat[c]; else { primalVariable_[c] = PrimalVariables::RS; }
cond[c].setFreeGas();
}
} }
} }
@ -1381,26 +1359,24 @@ namespace {
const V rvSat0 = fluidRvSat(p_old, cells_); const V rvSat0 = fluidRvSat(p_old, cells_);
const V rvSat = fluidRvSat(p, cells_); const V rvSat = fluidRvSat(p, cells_);
if (vapoil) { if (has_vapoil_) {
// The obvious case // The obvious case
auto ix0 = (so > 0 && isRv == 0); auto hasOil = (so > 0 && isRv == 0);
// keep oil saturated if previous sg is sufficient large: // keep oil saturated if previous so is sufficient large:
const int pos = pu.phase_pos[ Oil ]; const int pos = pu.phase_pos[ Oil ];
auto ix1 = (so < 0 && s_old.col(pos) > epsilon ); auto hadOil = (so < 0 && s_old.col(pos) > epsilon );
// Set oil saturated if previous rs is sufficiently large // Set oil saturated if previous rv is sufficiently large
const V rv_old = Eigen::Map<const V>(&state.rv()[0], nc); const V rv_old = Eigen::Map<const V>(&state.rv()[0], nc);
auto ix2 = ( (rv > rvSat * (1+epsilon) && isRv == 1) && (rv_old > rvSat0 * (1-epsilon)) ); auto oilCondensed = ( (rv > rvSat * (1+epsilon) && isRv == 1) && (rv_old > rvSat0 * (1-epsilon)) );
auto oilPresent = watOnly || ix0 || ix1 || ix2; auto useSg = watOnly || hasOil || hadOil || oilCondensed;
for (int c = 0; c < nc; ++c) { for (int c = 0; c < nc; ++c) {
if (oilPresent[c]) { if (useSg[c]) { rv[c] = rvSat[c]; }
rv[c] = rvSat[c]; else {primalVariable_[c] = PrimalVariables::RV; }
cond[c].setFreeOil();
}
} }
} }
std::copy(&cond[0], &cond[0] + nc, phaseCondition_.begin());
auto ixg = sg < 0; auto ixg = sg < 0;
for (int c = 0; c < nc; ++c) { for (int c = 0; c < nc; ++c) {
@ -1449,10 +1425,10 @@ namespace {
} }
// Rs and Rv updates // Rs and Rv updates
if (disgas) if (has_disgas_)
std::copy(&rs[0], &rs[0] + nc, state.gasoilratio().begin()); std::copy(&rs[0], &rs[0] + nc, state.gasoilratio().begin());
if (vapoil) if (has_vapoil_)
std::copy(&rv[0], &rv[0] + nc, state.rv().begin()); std::copy(&rv[0], &rv[0] + nc, state.rv().begin());
@ -2104,4 +2080,45 @@ namespace {
} }
template<class T>
void
FullyImplicitBlackoilSolver<T>::updatePrimalVariableFromState(const BlackoilState& state)
{
using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_);
const int np = state.numPhases();
const PhaseUsage& pu = fluid_.phaseUsage();
const DataBlock s = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
assert (active_[ Gas ]);
// reset the primary variables if RV and RS is not set Sg is used as primary variable.
primalVariable_.resize(nc);
std::fill(primalVariable_.begin(), primalVariable_.end(), PrimalVariables::Sg);
if (has_disgas_) {
// Oil/Gas or Water/Oil/Gas system
const V sg = s.col(pu.phase_pos[ Gas ]);
const V so = s.col(pu.phase_pos[ Oil ]);
for (V::Index c = 0, e = sg.size(); c != e; ++c) {
if ( sg[c] <= 0 && so[c] > 0 ) {primalVariable_[c] = PrimalVariables::RS; }
}
}
if (has_vapoil_) {
// Oil/Gas or Water/Oil/Gas system
const V sg = s.col(pu.phase_pos[ Gas ]);
const V so = s.col(pu.phase_pos[ Oil ]);
for (V::Index c = 0, e = so.size(); c != e; ++c) {
if (so[c] <= 0 && sg[c] > 0) {primalVariable_[c] = PrimalVariables::RV; }
}
}
}
} // namespace Opm } // namespace Opm

4
opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
View File

@ -74,7 +74,9 @@ namespace Opm
const RockCompressibility* rock_comp_props, const RockCompressibility* rock_comp_props,
WellsManager& wells_manager, WellsManager& wells_manager,
NewtonIterationBlackoilInterface& linsolver, NewtonIterationBlackoilInterface& linsolver,
const double* gravity); const double* gravity,
const bool disgas,
const bool vapoil );
/// Run the simulation. /// Run the simulation.
/// This will run succesive timesteps until timer.done() is true. It will /// This will run succesive timesteps until timer.done() is true. It will

16
opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp
View File

@ -66,7 +66,9 @@ namespace Opm
const RockCompressibility* rock_comp_props, const RockCompressibility* rock_comp_props,
WellsManager& wells_manager, WellsManager& wells_manager,
NewtonIterationBlackoilInterface& linsolver, NewtonIterationBlackoilInterface& linsolver,
const double* gravity); const double* gravity,
bool has_disgas,
bool has_vapoil );
SimulatorReport run(SimulatorTimer& timer, SimulatorReport run(SimulatorTimer& timer,
BlackoilState& state, BlackoilState& state,
@ -107,10 +109,12 @@ namespace Opm
const RockCompressibility* rock_comp_props, const RockCompressibility* rock_comp_props,
WellsManager& wells_manager, WellsManager& wells_manager,
NewtonIterationBlackoilInterface& linsolver, NewtonIterationBlackoilInterface& linsolver,
const double* gravity) const double* gravity,
const bool has_disgas,
const bool has_vapoil )
{ {
pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, linsolver, gravity)); pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, linsolver, gravity, has_disgas, has_vapoil));
} }
@ -192,7 +196,9 @@ namespace Opm
const RockCompressibility* rock_comp_props, const RockCompressibility* rock_comp_props,
WellsManager& wells_manager, WellsManager& wells_manager,
NewtonIterationBlackoilInterface& linsolver, NewtonIterationBlackoilInterface& linsolver,
const double* gravity) const double* gravity,
const bool has_disgas,
const bool has_vapoil)
: grid_(grid), : grid_(grid),
props_(props), props_(props),
rock_comp_props_(rock_comp_props), rock_comp_props_(rock_comp_props),
@ -200,7 +206,7 @@ namespace Opm
wells_(wells_manager.c_wells()), wells_(wells_manager.c_wells()),
gravity_(gravity), gravity_(gravity),
geo_(grid_, props_, gravity_), geo_(grid_, props_, gravity_),
solver_(param, grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver) solver_(param, grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver, has_disgas, has_vapoil)
/* param.getDefault("nl_pressure_residual_tolerance", 0.0), /* param.getDefault("nl_pressure_residual_tolerance", 0.0),
param.getDefault("nl_pressure_change_tolerance", 1.0), param.getDefault("nl_pressure_change_tolerance", 1.0),
param.getDefault("nl_pressure_maxiter", 10), param.getDefault("nl_pressure_maxiter", 10),