mirror of
https://github.com/OPM/opm-simulators.git
synced 2024-11-29 04:23:48 -06:00
Merge branch 'limited_variable_update' of github.com:GitPaean/opm-autodiff into oscillation_treatment_withlimitedupdate
Conflicts: opm/autodiff/FullyImplicitBlackoilSolver.hpp
This commit is contained in:
Kai Bao
commit
9557db4c15
@ -105,7 +105,7 @@ try
|
||||
Opm::LinearSolverFactory linsolver(param);
|
||||
Opm::NewtonIterationBlackoilSimple fis_solver(linsolver);
|
||||
|
||||
Opm::FullyImplicitBlackoilSolver<UnstructuredGrid> solver(*g, props, geo, 0, *wells, fis_solver);
|
||||
Opm::FullyImplicitBlackoilSolver<UnstructuredGrid> solver(param, *g, props, geo, 0, *wells, fis_solver);
|
||||
|
||||
Opm::BlackoilState state;
|
||||
initStateBasic(*g, props0, param, 0.0, state);
|
||||
|
||||
@ -31,6 +31,7 @@ struct Wells;
|
||||
|
||||
namespace Opm {
|
||||
|
||||
namespace parameter { class ParameterGroup; }
|
||||
class DerivedGeology;
|
||||
class RockCompressibility;
|
||||
class NewtonIterationBlackoilInterface;
|
||||
@ -56,13 +57,15 @@ namespace Opm {
|
||||
/// Construct a solver. It will retain references to the
|
||||
/// arguments of this functions, and they are expected to
|
||||
/// remain in scope for the lifetime of the solver.
|
||||
/// \param[in] param parameters
|
||||
/// \param[in] grid grid data structure
|
||||
/// \param[in] fluid fluid properties
|
||||
/// \param[in] geo rock properties
|
||||
/// \param[in] rock_comp_props if non-null, rock compressibility properties
|
||||
/// \param[in] wells well structure
|
||||
/// \param[in] linsolver linear solver
|
||||
FullyImplicitBlackoilSolver(const Grid& grid ,
|
||||
FullyImplicitBlackoilSolver(const parameter::ParameterGroup& param,
|
||||
const Grid& grid ,
|
||||
const BlackoilPropsAdInterface& fluid,
|
||||
const DerivedGeology& geo ,
|
||||
const RockCompressibility* rock_comp_props,
|
||||
@ -141,6 +144,9 @@ namespace Opm {
|
||||
HelperOps ops_;
|
||||
const WellOps wops_;
|
||||
const M grav_;
|
||||
double dp_max_rel_;
|
||||
double ds_max_;
|
||||
double drs_max_rel_;
|
||||
|
||||
std::vector<ReservoirResidualQuant> rq_;
|
||||
std::vector<PhasePresence> phaseCondition_;
|
||||
@ -278,6 +284,9 @@ namespace Opm {
|
||||
|
||||
void stablizeNewton(V &dx, V &dxOld, const bool &oscillate, const bool &stagnate, const double omega,
|
||||
const RelaxType relax_type) const;
|
||||
const double dpMaxRel () const { return dp_max_rel_; }
|
||||
const double dsMax () const { return ds_max_; }
|
||||
const double drsMaxRel () const { return drs_max_rel_; }
|
||||
|
||||
};
|
||||
} // namespace Opm
|
||||
|
||||
@ -35,6 +35,7 @@
|
||||
#include <opm/core/utility/Exceptions.hpp>
|
||||
#include <opm/core/utility/Units.hpp>
|
||||
#include <opm/core/well_controls.h>
|
||||
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
||||
|
||||
#include <cassert>
|
||||
#include <cmath>
|
||||
@ -212,7 +213,8 @@ namespace {
|
||||
|
||||
template<class T>
|
||||
FullyImplicitBlackoilSolver<T>::
|
||||
FullyImplicitBlackoilSolver(const Grid& grid ,
|
||||
FullyImplicitBlackoilSolver(const parameter::ParameterGroup& param,
|
||||
const Grid& grid ,
|
||||
const BlackoilPropsAdInterface& fluid,
|
||||
const DerivedGeology& geo ,
|
||||
const RockCompressibility* rock_comp_props,
|
||||
@ -230,15 +232,25 @@ namespace {
|
||||
, ops_ (grid)
|
||||
, wops_ (wells)
|
||||
, grav_ (gravityOperator(grid_, ops_, geo_))
|
||||
, dp_max_rel_ ( 1.0e9 )
|
||||
, ds_max_ ( 0.2 )
|
||||
, drs_max_rel_ ( 1.0e9 )
|
||||
, rq_ (fluid.numPhases())
|
||||
, phaseCondition_(AutoDiffGrid::numCells(grid))
|
||||
, residual_ ( { std::vector<ADB>(fluid.numPhases(), ADB::null()),
|
||||
ADB::null(),
|
||||
ADB::null() } )
|
||||
{
|
||||
if (param.has("dp_max_rel")) {
|
||||
dp_max_rel_ = param.get<double>(std::string("dp_max_rel"));
|
||||
}
|
||||
if (param.has("ds_max")) {
|
||||
ds_max_ = param.get<double>("ds_max");
|
||||
}
|
||||
if (param.has("drs_max_rel")) {
|
||||
drs_max_rel_ = param.get<double>("drs_max_rel");
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
template<class T>
|
||||
@ -1251,7 +1263,7 @@ namespace {
|
||||
assert(varstart == dx.size());
|
||||
|
||||
// Pressure update.
|
||||
const double dpmaxrel = 0.8;
|
||||
const double dpmaxrel = dpMaxRel();
|
||||
const V p_old = Eigen::Map<const V>(&state.pressure()[0], nc, 1);
|
||||
const V absdpmax = dpmaxrel*p_old.abs();
|
||||
const V dp_limited = sign(dp) * dp.abs().min(absdpmax);
|
||||
@ -1262,40 +1274,53 @@ namespace {
|
||||
// Saturation updates.
|
||||
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
|
||||
const DataBlock s_old = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
|
||||
const double dsmax = 0.3;
|
||||
const double dsmax = dsMax();
|
||||
V so = one;
|
||||
V sw;
|
||||
V sg;
|
||||
|
||||
{
|
||||
V maxVal = zero;
|
||||
V dso = zero;
|
||||
if (active_[Water]){
|
||||
maxVal = dsw.abs().max(maxVal);
|
||||
dso = dso - dsw;
|
||||
}
|
||||
|
||||
V dsg;
|
||||
if (active_[Gas]){
|
||||
dsg = isSg * dxvar - isRv * dsw;
|
||||
maxVal = dsg.abs().max(maxVal);
|
||||
dso = dso - dsg;
|
||||
}
|
||||
|
||||
maxVal = dso.abs().max(maxVal);
|
||||
|
||||
V step = dsmax/maxVal;
|
||||
step = step.min(1.);
|
||||
|
||||
if (active_[Water]) {
|
||||
const int pos = pu.phase_pos[ Water ];
|
||||
const V sw_old = s_old.col(pos);
|
||||
const V dsw_limited = sign(dsw) * dsw.abs().min(dsmax);
|
||||
sw = (sw_old - dsw_limited).unaryExpr(Chop01());
|
||||
sw = sw_old - step * dsw;
|
||||
so -= sw;
|
||||
for (int c = 0; c < nc; ++c) {
|
||||
state.saturation()[c*np + pos] = sw[c];
|
||||
}
|
||||
}
|
||||
|
||||
V sg;
|
||||
if (active_[Gas]) {
|
||||
const int pos = pu.phase_pos[ Gas ];
|
||||
const V sg_old = s_old.col(pos);
|
||||
const V dsg = isSg * dxvar - isRv * dsw;
|
||||
const V dsg_limited = sign(dsg) * dsg.abs().min(dsmax);
|
||||
sg = sg_old - dsg_limited;
|
||||
sg = sg_old - step * dsg;
|
||||
so -= sg;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
const double drsmax = 1e9;
|
||||
const double drsmaxrel = drsMaxRel();
|
||||
const double drvmax = 1e9;//% same as in Mrst
|
||||
V rs;
|
||||
if (disgas) {
|
||||
const V rs_old = Eigen::Map<const V>(&state.gasoilratio()[0], nc);
|
||||
const V drs = isRs * dxvar;
|
||||
const V drs_limited = sign(drs) * drs.abs().min(drsmax);
|
||||
const V drs_limited = sign(drs) * drs.abs().min(rs_old.abs()*drsmaxrel);
|
||||
rs = rs_old - drs_limited;
|
||||
}
|
||||
V rv;
|
||||
@ -1431,7 +1456,8 @@ namespace {
|
||||
|
||||
// Bhp update.
|
||||
const V bhp_old = Eigen::Map<const V>(&well_state.bhp()[0], nw, 1);
|
||||
const V bhp = bhp_old - dbhp;
|
||||
const V dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel);
|
||||
const V bhp = bhp_old - dbhp_limited;
|
||||
std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
|
||||
|
||||
}
|
||||
|
||||
@ -200,7 +200,7 @@ namespace Opm
|
||||
wells_(wells_manager.c_wells()),
|
||||
gravity_(gravity),
|
||||
geo_(grid_, props_, gravity_),
|
||||
solver_(grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver)
|
||||
solver_(param, grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver)
|
||||
/* param.getDefault("nl_pressure_residual_tolerance", 0.0),
|
||||
param.getDefault("nl_pressure_change_tolerance", 1.0),
|
||||
param.getDefault("nl_pressure_maxiter", 10),
|
||||
|
||||
Loading…
Reference in New Issue
Block a user