diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index da83ef86a..234dd6a4e 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -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& primary_variables, + const BVectorWell& dwells); + + // calculate a relaxation factor to avoid overshoot of total rates + static double relaxationFactorRate(const std::vector& primary_variables, + const BVectorWell& dwells); }; } diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index a94c56612..aa25485e8 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -870,29 +870,37 @@ namespace Opm const std::vector 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 + double + StandardWell:: + 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 + double + StandardWell:: + relaxationFactorFractionsProducer(const std::vector& 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 + double + StandardWell:: + relaxationFactorRate(const std::vector& 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; + } }