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 #1597 from GitPaean/avoiding_overshoot_fractions

Browse Source 
relaxation to avoid overshoot of fractions of StandardWell
This commit is contained in:
Atgeirr Flø Rasmussen 2018-11-12 14:20:52 +01:00 committed by GitHub
commit d562df9d8a
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 129 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
opm/autodiff/StandardWell.hpp
View File

@ -354,6 +354,18 @@ namespace Opm
// handle the non reasonable fractions due to numerical overshoot
void processFractions() const;
// relaxation factor considering only one fraction value
static double relaxationFactorFraction(const double old_value,
const double dx);
// calculate a relaxation factor to avoid overshoot of the fractions for producers
// which might result in negative rates
static double relaxationFactorFractionsProducer(const std::vector<double>& primary_variables,
const BVectorWell& dwells);
// calculate a relaxation factor to avoid overshoot of total rates
static double relaxationFactorRate(const std::vector<double>& primary_variables,
const BVectorWell& dwells);
};
}

122
opm/autodiff/StandardWell_impl.hpp
View File

@ -870,29 +870,37 @@ namespace Opm
const std::vector<double> old_primary_variables = primary_variables_;
// for injectors, very typical one of the fractions will be one, and it is easy to get zero value
// fractions. not sure what is the best way to handle it yet, so we just use 1.0 here
const double relaxation_factor_fractions = (well_type_ == PRODUCER) ?
relaxationFactorFractionsProducer(old_primary_variables, dwells)
: 1.0;
// update the second and third well variable (The flux fractions)
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int sign2 = dwells[0][WFrac] > 0 ? 1: -1;
const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac]),dFLimit);
const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac] * relaxation_factor_fractions), dFLimit);
// primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int sign3 = dwells[0][GFrac] > 0 ? 1: -1;
const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac]),dFLimit);
const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac] * relaxation_factor_fractions), dFLimit);
primary_variables_[GFrac] = old_primary_variables[GFrac] - dx3_limited;
}
if (has_solvent) {
const int sign4 = dwells[0][SFrac] > 0 ? 1: -1;
const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]),dFLimit);
const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]) * relaxation_factor_fractions, dFLimit);
primary_variables_[SFrac] = old_primary_variables[SFrac] - dx4_limited;
}
processFractions();
// updating the total rates G_t
primary_variables_[WQTotal] = old_primary_variables[WQTotal] - dwells[0][WQTotal];
// updating the total rates Q_t
const double relaxation_factor_rate = relaxationFactorRate(old_primary_variables, dwells);
primary_variables_[WQTotal] = old_primary_variables[WQTotal] - dwells[0][WQTotal] * relaxation_factor_rate;
// updating the bottom hole pressure
{
@ -2215,4 +2223,108 @@ namespace Opm
}
}
}
template<typename TypeTag>
double
StandardWell<TypeTag>::
relaxationFactorFraction(const double old_value,
const double dx)
{
assert(old_value >= 0. && old_value <= 1.0);
double relaxation_factor = 1.;
// updated values without relaxation factor
const double possible_updated_value = old_value - dx;
// 0.95 is an experimental value remains to be optimized
if (possible_updated_value < 0.0) {
relaxation_factor = std::abs(old_value / dx) * 0.95;
} else if (possible_updated_value > 1.0) {
relaxation_factor = std::abs((1. - old_value) / dx) * 0.95;
}
// if possible_updated_value is between 0. and 1.0, then relaxation_factor
// remains to be one
assert(relaxation_factor >= 0. && relaxation_factor <= 1.);
return relaxation_factor;
}
template<typename TypeTag>
double
StandardWell<TypeTag>::
relaxationFactorFractionsProducer(const std::vector<double>& primary_variables,
const BVectorWell& dwells)
{
// TODO: not considering solvent yet
// 0.95 is a experimental value, which remains to be optimized
double relaxation_factor = 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const double relaxation_factor_w = relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]);
relaxation_factor = std::min(relaxation_factor, relaxation_factor_w);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const double relaxation_factor_g = relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]);
relaxation_factor = std::min(relaxation_factor, relaxation_factor_g);
}
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
// We need to make sure the even with the relaxation_factor, the sum of F_w and F_g is below one, so there will
// not be negative oil fraction later
const double original_sum = primary_variables[WFrac] + primary_variables[GFrac];
const double relaxed_update = (dwells[0][WFrac] + dwells[0][GFrac]) * relaxation_factor;
const double possible_updated_sum = original_sum - relaxed_update;
if (possible_updated_sum > 1.0) {
assert(relaxed_update != 0.);
const double further_relaxation_factor = std::abs((1. - original_sum) / relaxed_update) * 0.95;
relaxation_factor *= further_relaxation_factor;
}
}
assert(relaxation_factor >= 0.0 && relaxation_factor <= 1.0);
return relaxation_factor;
}
template<typename TypeTag>
double
StandardWell<TypeTag>::
relaxationFactorRate(const std::vector<double>& primary_variables,
const BVectorWell& dwells)
{
double relaxation_factor = 1.0;
// For injector, we only check the total rates to avoid sign change of rates
const double original_total_rate = primary_variables[WQTotal];
const double newton_update = dwells[0][WQTotal];
const double possible_update_total_rate = primary_variables[WQTotal] - newton_update;
// 0.8 here is a experimental value, which remains to be optimized
// if the original rate is zero or possible_update_total_rate is zero, relaxation_factor will
// always be 1.0, more thoughts might be needed.
if (original_total_rate * possible_update_total_rate < 0.) { // sign changed
relaxation_factor = std::abs(original_total_rate / newton_update) * 0.8;
}
assert(relaxation_factor >= 0.0 && relaxation_factor <= 1.0);
return relaxation_factor;
}
}