addressing the review comments #1597

This commit is contained in:
Kai Bao
2018-11-06 21:15:37 +01:00
parent 7e17a60c58
commit 179a505b83
2 changed files with 13 additions and 27 deletions

View File

@@ -363,7 +363,7 @@ namespace Opm
static double determineRelaxationFactorProducer(const std::vector<double>& primary_variables, static double determineRelaxationFactorProducer(const std::vector<double>& primary_variables,
const BVectorWell& dwells); const BVectorWell& dwells);
// calcualte a relaxation factor to avoid overshoot for injectors // calculate a relaxation factor to avoid overshoot for injectors
static double determineRelaxationFactorInjector(const std::vector<double>& primary_variables, static double determineRelaxationFactorInjector(const std::vector<double>& primary_variables,
const BVectorWell& dwells); const BVectorWell& dwells);
}; };

View File

@@ -2228,35 +2228,26 @@ namespace Opm
double double
StandardWell<TypeTag>:: StandardWell<TypeTag>::
relaxationFactorFraction(const double old_value, relaxationFactorFraction(const double old_value,
const double dx) { const double dx)
// may adding a tolerance? {
assert(old_value >= 0. && old_value <= 1.0); assert(old_value >= 0. && old_value <= 1.0);
// to avoid to get negative value or value above 1
double relaxation_factor = 1.; double relaxation_factor = 1.;
bool relaxation_avoid_negative = false; // updated values without relaxation factor
bool relaxation_avoid_over_one = false; const double possible_updated_value = old_value - dx;
{
const double possible_updated_value = old_value - dx;
relaxation_avoid_negative = possible_updated_value < 0.;
relaxation_avoid_over_one = possible_updated_value > 1.0;
if (!relaxation_avoid_negative && !relaxation_avoid_over_one) { // 0.95 is an experimental value remains to be optimized
return relaxation_factor; // 1.0 if (possible_updated_value < 0.0) {
}
}
// we will not be this step if dx == 0.
assert(dx != 0.0);
if (relaxation_avoid_negative) {
relaxation_factor = std::abs(old_value / dx) * 0.95; relaxation_factor = std::abs(old_value / dx) * 0.95;
} else if (relaxation_avoid_over_one) { } else if (possible_updated_value > 1.0) {
relaxation_factor = std::abs((1. - old_value) / dx) * 0.95; 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.); assert(relaxation_factor >= 0. && relaxation_factor <= 1.);
return relaxation_factor; return relaxation_factor;
} }
@@ -2276,16 +2267,12 @@ namespace Opm
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const double relaxation_factor_w = relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]); const double relaxation_factor_w = relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]);
if (relaxation_factor_w < relaxation_factor) { relaxation_factor = std::min(relaxation_factor, relaxation_factor_w);
relaxation_factor = relaxation_factor_w;
}
} }
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const double relaxation_factor_g = relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]); const double relaxation_factor_g = relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]);
if (relaxation_factor_g < relaxation_factor) { relaxation_factor = std::min(relaxation_factor, relaxation_factor_g);
relaxation_factor = relaxation_factor_g;
}
} }
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
@@ -2303,7 +2290,6 @@ namespace Opm
} }
} }
// relaxation factor for the total rates // relaxation factor for the total rates
{ {
const double original_total_rate = primary_variables[WQTotal]; const double original_total_rate = primary_variables[WQTotal];