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

New updateState

Browse Source 
-- a new updateState is implemented based on dune vectors
-- the old is kept for comparision in this PR
-- the updateState is not identical.
Tested on spe1, spe9 and norne and it improves the convergence compares
to the old one.
This commit is contained in:
Tor Harald Sandve
2016-09-07 10:13:24 +02:00
parent 952ccf8338
commit 83ff3271af
4 changed files with 467 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

186
opm/autodiff/BlackoilModelEbos.hpp
View File

@@ -171,6 +171,7 @@ namespace Opm {
, terminal_output_ (terminal_output)
, current_relaxation_(1.0)
, isBeginReportStep_(false)
, dx_old_(AutoDiffGrid::numCells(grid_))
{
const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
const std::vector<double> pv(geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size());
@@ -218,7 +219,7 @@ namespace Opm {
// the mass balance for each active phase, the well flux and the well equations.
residual_norms_history_.clear();
current_relaxation_ = 1.0;
dx_old_ = V::Zero(sizeNonLinear());
dx_old_ = 0.0; //V::Zero(sizeNonLinear());
}
IterationReport iter_report = assemble(timer, iteration, reservoir_state, well_state);
std::vector<double> residual_norms;
@@ -231,7 +232,11 @@ namespace Opm {
//residual_.singlePrecision = (unit::convert::to(dt, unit::day) < 20.) ;
// Compute the nonlinear update.
V dx = solveJacobianSystem(result);
const int nc = AutoDiffGrid::numCells(grid_);
const int nw = wellModel().wells().number_of_wells;
BVector x(nc);
BVector xw(nw);
V dx = solveJacobianSystem(result, x, xw);
// Stabilize the nonlinear update.
bool isOscillate = false;
@@ -246,16 +251,38 @@ namespace Opm {
OpmLog::info(msg);
}
}
nonlinear_solver.stabilizeNonlinearUpdate(dx, dx_old_, current_relaxation_);
nonlinear_solver.stabilizeNonlinearUpdate(x, dx_old_, current_relaxation_);
// Apply the update, applying model-dependent
// limitations and chopping of the update.
updateState(dx, reservoir_state, well_state);
//auto well_state2 = well_state;
//auto reservoir_state2 = reservoir_state;
//updateState(dx, reservoir_state2, well_state);
updateState(x,reservoir_state);
wellModel().updateWellState(xw, well_state);
// for (int c = 0; c < nc; ++c) {
// printIf(c, reservoir_state.pressure()[c], reservoir_state2.pressure()[c], 1e5, "pressure");
// printIf(c, reservoir_state.rv()[c], reservoir_state2.rv()[c] , 1e-3, "rv");
// printIf(c, reservoir_state.gasoilratio()[c], reservoir_state2.gasoilratio()[c], 1, "rs");
// printIf(c, reservoir_state.saturation()[3*c+0], reservoir_state2.saturation()[3*c+0], 1e-3, "sw");
// printIf(c, reservoir_state.saturation()[3*c+1], reservoir_state2.saturation()[3*c+1], 1e-3, "so");
// printIf(c, reservoir_state.saturation()[3*c+2], reservoir_state2.saturation()[3*c+2], 1e-3, "sg");
// printIf(c, reservoir_state.hydroCarbonState()[c], reservoir_state2.hydroCarbonState()[c], 1e-6,"state");
// }
}
const bool failed = false; // Not needed in this model.
const int linear_iters = must_solve ? result.iterations : 0;
return IterationReport{ failed, converged, linear_iters, iter_report.well_iterations };
}
void printIf(int c, double x, double y, double eps, std::string type) {
if (std::abs(x-y) > eps) {
std::cout << type << " " <<c << ": "<<x << " " << y << std::endl;
}
}
/// Called once after each time step.
@@ -369,7 +396,7 @@ namespace Opm {
/// Solve the Jacobian system Jx = r where J is the Jacobian and
/// r is the residual.
V solveJacobianSystem(Dune::InverseOperatorResult& result) const
V solveJacobianSystem(Dune::InverseOperatorResult& result, BVector& x, BVector& xw) const
{
typedef double Scalar;
@@ -401,12 +428,13 @@ namespace Opm {
const int nc = AutoDiffGrid::numCells(grid_);
// Solve system.
BVector x(ebosJac.M());
//BVector x(ebosJac.M());
x = 0.0;
linsolve.apply(x, ebosResid, result);
const int nw = wellModel().wells().number_of_wells;
BVector xw(nw);
//BVector xw(nw);
xw = 0.0;
wellModel().recoverVariable(x, xw);
// convert to V;
@@ -432,6 +460,148 @@ namespace Opm {
return dx;
}
/// Apply an update to the primary variables, chopped if appropriate.
/// \param[in] dx updates to apply to primary variables
/// \param[in, out] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
void updateState(const BVector& dx,
ReservoirState& reservoir_state)
{
using namespace Opm::AutoDiffGrid;
const int np = fluid_.numPhases();
const int nc = numCells(grid_);
for (int cell_idx = 0; cell_idx < nc; ++cell_idx) {
double dp = dx[cell_idx][flowPhaseToEbosCompIdx(0)];
reservoir_state.pressure()[cell_idx] -= dp;
// Saturation updates.
const double dsw = active_[Water] ? dx[cell_idx][flowPhaseToEbosCompIdx(1)] : 0.0;
const int xvar_ind = active_[Water] ? 2 : 1;
const double dxvar = active_[Gas] ? dx[cell_idx][flowPhaseToEbosCompIdx(xvar_ind)] : 0.0;
double dso = 0.0;
double dsg = 0.0;
double drs = 0.0;
double drv = 0.0;
double maxVal = 0.0;
// water phase
maxVal = std::max(std::abs(dsw),maxVal);
dso -= dsw;
// gas phase
switch (reservoir_state.hydroCarbonState()[cell_idx]) {
case HydroCarbonState::GasAndOil:
dsg = dxvar;
break;
case HydroCarbonState::OilOnly:
drs = dxvar;
break;
case HydroCarbonState::GasOnly:
dsg -= dsw;
drv = dxvar;
break;
default:
OPM_THROW(std::logic_error, "Unknown primary variable enum value in cell " << cell_idx << ": " << reservoir_state.hydroCarbonState()[cell_idx]);
}
dso -= dsg;
// Appleyard chop process.
maxVal = std::max(std::abs(dsg),maxVal);
double step = dsMax()/maxVal;
step = std::min(step, 1.0);
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
if (active_[Water]) {
double& sw = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Water ]];
sw -= step * dsw;
}
if (active_[Gas]) {
double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]];
sg -= step * dsg;
}
double& so = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Oil ]];
so -= step * dso;
// const double drmaxrel = drMaxRel();
// Update rs and rv
if (has_disgas_) {
double& rs = reservoir_state.gasoilratio()[cell_idx];
rs -= drs;
}
if (has_vapoil_) {
double& rv = reservoir_state.rv()[cell_idx];
rv -= drv;
}
// Sg is used as primal variable for water only cells.
const double epsilon = 1e-4; //std::sqrt(std::numeric_limits<double>::epsilon());
// phase translation sg <-> rs
const HydroCarbonState hydroCarbonState = reservoir_state.hydroCarbonState()[cell_idx];
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& rs = reservoir_state.gasoilratio()[cell_idx];
rs = rsSat*(1-epsilon);
} else if (so <= 0.0 && has_vapoil_) {
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasOnly; // sg --> rv
so = 0;
sg = 1.0 - sw - so;
double& rv = reservoir_state.rv()[cell_idx];
double rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(0, 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 (rs > ( rsSat * (1+epsilon) ) ) {
reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
sg = epsilon;
so -= epsilon;
rs = rsSat;
}
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 (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;
}
default:
OPM_THROW(std::logic_error, "Unknown primary variable enum value in cell " << cell_idx << ": " << hydroCarbonState);
}
}
//wellModel().updateWellState(dwells, well_state);
// Update phase conditions used for property calculations.
updatePhaseCondFromPrimalVariable(reservoir_state);
}
/// Apply an update to the primary variables, chopped if appropriate.
/// \param[in] dx updates to apply to primary variables
/// \param[in, out] reservoir_state reservoir state variables
@@ -865,7 +1035,7 @@ namespace Opm {
std::vector<std::vector<double>> residual_norms_history_;
double current_relaxation_;
V dx_old_;
BVector dx_old_;

9
opm/autodiff/NonlinearSolver.hpp
View File

@@ -24,6 +24,9 @@
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/simulator/SimulatorTimerInterface.hpp>
#include <opm/autodiff/DuneMatrix.hpp>
#include <dune/common/fmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <memory>
namespace Opm {
@@ -40,6 +43,10 @@ namespace Opm {
typedef ADB::V V;
typedef ADB::M M;
typedef double Scalar;
typedef Dune::FieldVector<Scalar, 3 > VectorBlockType;
typedef Dune::BlockVector<VectorBlockType> BVector;
// Available relaxation scheme types.
enum RelaxType { DAMPEN, SOR };
@@ -138,6 +145,8 @@ namespace Opm {
/// Apply a stabilization to dx, depending on dxOld and relaxation parameters.
void stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega) const;
void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const;
/// The greatest relaxation factor (i.e. smallest factor) allowed.
double relaxMax() const { return param_.relax_max_; }

32
opm/autodiff/NonlinearSolver_impl.hpp
View File

@@ -278,6 +278,38 @@ namespace Opm
return;
}
template <class PhysicalModel>
void
NonlinearSolver<PhysicalModel>::stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const
{
// The dxOld is updated with dx.
// If omega is equal to 1., no relaxtion will be appiled.
BVector tempDxOld = dxOld;
dxOld = dx;
switch (relaxType()) {
case DAMPEN:
if (omega == 1.) {
return;
}
dx *= omega;
return;
case SOR:
if (omega == 1.) {
return;
}
dx *= omega;
tempDxOld *= (1.-omega);
dx += tempDxOld;
return;
default:
OPM_THROW(std::runtime_error, "Can only handle DAMPEN and SOR relaxation type.");
}
return;
}
} // namespace Opm

250
opm/autodiff/StandardWellsDense.hpp
View File

@@ -677,8 +677,16 @@ namespace Opm {
dx_new_eigen(idx) = dx_new[w][flowPhaseToEbosCompIdx(p)];
}
}
assert(dx.size() == total_residual_v.size());
updateWellState(dx_new_eigen.array(), well_state);
//auto well_state2 = well_state;
//updateWellState(dx_new_eigen.array(), well_state);
updateWellState(dx_new, well_state);
// for (int w = 0; w < nw; ++w) {
// printIf(w, well_state.wellSolutions()[w], well_state2.wellSolutions()[w], 1e-3 ,"well solution 1");
// printIf(w, well_state.wellSolutions()[nw + w], well_state2.wellSolutions()[nw + w], 1e-3 ,"well solution 2");
// printIf(w, well_state.wellSolutions()[2*nw + w], well_state2.wellSolutions()[2*nw + w], 1e-3 ,"well solution 3");
// }
updateWellControls(well_state);
setWellVariables(well_state);
}
@@ -693,6 +701,13 @@ namespace Opm {
return IterationReport{failed, converged, linear_iters, it};
}
void printIf(int c, double x, double y, double eps, std::string type) {
if (std::abs(x-y) > eps) {
std::cout << type << " " << c << ": "<<x << " " << y << std::endl;
}
}
std::vector<double> residual() {
const int np = numPhases();
@@ -880,6 +895,237 @@ namespace Opm {
}
}
template <class WellState>
void updateWellState(const BVector& dwells,
WellState& well_state)
{
if( localWellsActive() )
{
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
double dFLimit = 0.2;
double dBHPLimit = 2.0;
const Vector xvar_well_old = Eigen::Map<const Vector>(&well_state.wellSolutions()[0], nw*np);
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells().ctrls[w];
// The current control in the well state overrides
// the current control set in the Wells struct, which
// is instead treated as a default.
const int current = well_state.currentControls()[w];
const double target_rate = well_controls_iget_target(wc, current);
const double* distr = well_controls_iget_distr(wc, current);
std::vector<double> F(np,0.0);
const int sign2 = dwells[w][flowPhaseToEbosCompIdx(1)] > 0 ? 1: -1;
const double dx2_limited = sign2 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(1)]),dFLimit);
well_state.wellSolutions()[nw + w] = xvar_well_old[nw + w] - dx2_limited;
const int sign3 = dwells[w][flowPhaseToEbosCompIdx(2)] > 0 ? 1: -1;
const double dx3_limited = sign3 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(2)]),dFLimit);
well_state.wellSolutions()[2*nw + w] = xvar_well_old[2*nw + w] - dx3_limited;
F[Water] = well_state.wellSolutions()[nw + w];
F[Gas] = well_state.wellSolutions()[2*nw + w];
F[Oil] = 1.0 - F[Water] - F[Gas];
// const double dFw = dxvar_well[nw + w];
// const double dFg = dxvar_well[nw*2 + w];
// const double dFo = - dFw - dFg;
// //std::cout << w << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl;
// double step = dFLimit / std::max(std::abs(dFw),std::max(std::abs(dFg),std::abs(dFo))); //)) / dFLimit;
// step = std::min(step, 1.0);
// //std::cout << step << std::endl;
// F[Water] = xvar_well_old[nw + w] - step*dFw;
// F[Gas] = xvar_well_old[2*nw + w] - step*dFg;
// F[Oil] = (1.0 - xvar_well_old[2*nw + w] - xvar_well_old[nw + w]) - step * dFo;
if (F[Water] < 0.0) {
F[Gas] /= (1.0 - F[Water]);
F[Oil] /= (1.0 - F[Water]);
F[Water] = 0.0;
}
if (F[Gas] < 0.0) {
F[Water] /= (1.0 - F[Gas]);
F[Oil] /= (1.0 - F[Gas]);
F[Gas] = 0.0;
}
if (F[Oil] < 0.0) {
F[Water] /= (1.0 - F[Oil]);
F[Gas] /= (1.0 - F[Oil]);
F[Oil] = 0.0;
}
well_state.wellSolutions()[nw + w] = F[Water];
well_state.wellSolutions()[2*nw + w] = F[Gas];
//std::cout << wells().name[w] << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl;
std::vector<double> g = {1,1,0.01};
if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) {
for (int p = 0; p < np; ++p) {
F[p] /= distr[p];
}
} else {
for (int p = 0; p < np; ++p) {
F[p] /= g[p];
}
}
//std::cout << w << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl;
// const double dFw = dxvar_well[nw + w];
// const double dFg = dxvar_well[nw*2 + w];
// const double dFo = - dFw - dFg;
//std::cout << w << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl;
// double step = dFLimit / std::max(std::abs(dFw),std::max(std::abs(dFg),std::abs(dFo))); //)) / dFLimit;
// step = std::min(step, 1.0);
// std::cout << step << std::endl;
// F[Water] = xvar_well_old[nw + w] - step*dFw;
// F[Gas] = xvar_well_old[2*nw + w] - step*dFg;
// F[Oil] = (1.0 - xvar_well_old[2*nw + w] - xvar_well_old[nw + w]) - step * dFo;
// double sumF = F[Water]+F[Gas]+F[Oil];
// F[Water] /= sumF;
// F[Gas] /= sumF;
// F[Oil] /= sumF;
// well_state.wellSolutions()[nw + w] = F[Water];
// well_state.wellSolutions()[2 * nw + w] = F[Gas];
switch (well_controls_iget_type(wc, current)) {
case BHP:
{
//const int sign1 = dxvar_well[w] > 0 ? 1: -1;
//const double dx1_limited = sign1 * std::min(std::abs(dxvar_well[w]),std::abs(xvar_well_old[w])*dTotalRateLimit);
well_state.wellSolutions()[w] = xvar_well_old[w] - dwells[w][flowPhaseToEbosCompIdx(0)];
switch (wells().type[w]) {
case INJECTOR:
for (int p = 0; p < np; ++p) {
const double comp_frac = wells().comp_frac[np*w + p];
//if (comp_frac > 0) {
well_state.wellRates()[w*np + p] = comp_frac * well_state.wellSolutions()[w];
//}
}
break;
case PRODUCER:
for (int p = 0; p < np; ++p) {
well_state.wellRates()[w*np + p] = well_state.wellSolutions()[w] * F[p];
}
break;
}
}
break;
case SURFACE_RATE:
{
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;
//const int sign = (dxvar_well1[w] < 0) ? -1 : 1;
//well_state.bhp()[w] -= sign * std::min( std::abs(dxvar_well1[w]), std::abs(well_state.bhp()[w])*dpmaxrel) ;
well_state.bhp()[w] = well_state.wellSolutions()[w];
if (wells().type[w]==PRODUCER) {
double F_target = 0.0;
for (int p = 0; p < np; ++p) {
F_target += wells().comp_frac[np*w + p] * F[p];
}
for (int p = 0; p < np; ++p) {
//std::cout << F[p] << std::endl;
well_state.wellRates()[np*w + p] = F[p] * target_rate /F_target;
}
} else {
for (int p = 0; p < np; ++p) {
//std::cout << wells().comp_frac[np*w + p] << " " <<distr[p] << std::endl;
well_state.wellRates()[w*np + p] = wells().comp_frac[np*w + p] * target_rate;
}
}
}
break;
case RESERVOIR_RATE: {
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;
//const int sign = (dxvar_well1[w] < 0) ? -1 : 1;
//well_state.bhp()[w] -= sign * std::min( std::abs(dxvar_well1[w]), std::abs(well_state.bhp()[w])*dpmaxrel) ;
well_state.bhp()[w] = well_state.wellSolutions()[w];
for (int p = 0; p < np; ++p) {
well_state.wellRates()[np*w + p] = F[p] * target_rate;
}
}
break;
}
}
const Opm::PhaseUsage& pu = fluid_->phaseUsage();
//Loop over all wells
#pragma omp parallel for schedule(static)
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells().ctrls[w];
const int nwc = well_controls_get_num(wc);
//Loop over all controls until we find a THP control
//that specifies what we need...
//Will only update THP for wells with THP control
for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) {
if (well_controls_iget_type(wc, ctrl_index) == THP) {
double aqua = 0.0;
double liquid = 0.0;
double vapour = 0.0;
if ((*active_)[ Water ]) {
aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ];
}
if ((*active_)[ Oil ]) {
liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ];
}
if ((*active_)[ Gas ]) {
vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ];
}
double alq = well_controls_iget_alq(wc, ctrl_index);
int table_id = well_controls_iget_vfp(wc, ctrl_index);
const WellType& well_type = wells().type[w];
if (well_type == INJECTOR) {
double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(),
wellPerforationDensities(), gravity_);
well_state.thp()[w] = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, well_state.bhp()[w] + dp);
}
else if (well_type == PRODUCER) {
double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(),
wellPerforationDensities(), gravity_);
well_state.thp()[w] = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, well_state.bhp()[w] + dp, alq);
}
else {
OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well");
}
//Assume only one THP control specified for each well
break;
}
}
}
}
}
template <class WellState>