From 7e17a60c58c74a86cfa13f74ede07066b91c94c6 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 16 Oct 2018 13:17:13 +0200 Subject: [PATCH 1/3] relaxation to avoid overshoot during Newton update in StandardWell Solvent model is not handled yet. --- opm/autodiff/StandardWell.hpp | 12 +++ opm/autodiff/StandardWell_impl.hpp | 136 ++++++++++++++++++++++++++++- 2 files changed, 144 insertions(+), 4 deletions(-) diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index da83ef86a..921eeaa15 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 + // which might result in negative rates + static double determineRelaxationFactorProducer(const std::vector& primary_variables, + const BVectorWell& dwells); + + // calcualte a relaxation factor to avoid overshoot for injectors + static double determineRelaxationFactorInjector(const std::vector& primary_variables, + const BVectorWell& dwells); }; } diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index a94c56612..5ae33590a 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -870,29 +870,34 @@ namespace Opm const std::vector old_primary_variables = primary_variables_; + const double relaxation_factor = (well_type_ == PRODUCER) ? + determineRelaxationFactorProducer(old_primary_variables, dwells) + : determineRelaxationFactorInjector(old_primary_variables, dwells); + // 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), 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), 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, 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]; + primary_variables_[WQTotal] = old_primary_variables[WQTotal] - dwells[0][WQTotal] * relaxation_factor; // updating the bottom hole pressure { @@ -2215,4 +2220,127 @@ namespace Opm } } } + + + + + template + double + StandardWell:: + relaxationFactorFraction(const double old_value, + const double dx) { + // may adding a tolerance? + assert(old_value >= 0. && old_value <= 1.0); + + // to avoid to get negative value or value above 1 + double relaxation_factor = 1.; + + bool relaxation_avoid_negative = false; + bool relaxation_avoid_over_one = false; + { + 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) { + return relaxation_factor; // 1.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; + } else if (relaxation_avoid_over_one) { + relaxation_factor = std::abs((1. - old_value) / dx) * 0.95; + } + + assert(relaxation_factor >= 0. && relaxation_factor <= 1.); + return relaxation_factor; + } + + + + + + template + double + StandardWell:: + determineRelaxationFactorProducer(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]); + if (relaxation_factor_w < relaxation_factor) { + relaxation_factor = relaxation_factor_w; + } + } + + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const double relaxation_factor_g = relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]); + if (relaxation_factor_g < relaxation_factor) { + 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; + } + } + + + // relaxation factor for the total rates + { + const double original_total_rate = primary_variables[WQTotal]; + const double relaxed_update = dwells[0][WQTotal] * relaxation_factor; + const double possible_update_total_rate = primary_variables[WQTotal] - relaxed_update; + + if (original_total_rate * possible_update_total_rate < 0.) { // sign changed + const double further_relaxation_factor = std::abs(original_total_rate / relaxed_update) * 0.95; + relaxation_factor *= further_relaxation_factor; + } + } + + return relaxation_factor; + } + + + + + + template + double + StandardWell:: + determineRelaxationFactorInjector(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 (original_total_rate * possible_update_total_rate < 0.) { // sign changed + relaxation_factor = std::abs(original_total_rate / newton_update) * 0.8; + } + + return relaxation_factor; + } } From 179a505b834fc2069c2128cf969057635a3013ed Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 6 Nov 2018 21:15:37 +0100 Subject: [PATCH 2/3] addressing the review comments #1597 --- opm/autodiff/StandardWell.hpp | 2 +- opm/autodiff/StandardWell_impl.hpp | 38 ++++++++++-------------------- 2 files changed, 13 insertions(+), 27 deletions(-) diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index 921eeaa15..4f4369050 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -363,7 +363,7 @@ namespace Opm static double determineRelaxationFactorProducer(const std::vector& primary_variables, 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& primary_variables, const BVectorWell& dwells); }; diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index 5ae33590a..db768ae6e 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -2228,35 +2228,26 @@ namespace Opm double StandardWell:: relaxationFactorFraction(const double old_value, - const double dx) { - // may adding a tolerance? + const double dx) + { assert(old_value >= 0. && old_value <= 1.0); - // to avoid to get negative value or value above 1 double relaxation_factor = 1.; - bool relaxation_avoid_negative = false; - bool relaxation_avoid_over_one = false; - { - const double possible_updated_value = old_value - dx; - relaxation_avoid_negative = possible_updated_value < 0.; - relaxation_avoid_over_one = possible_updated_value > 1.0; + // updated values without relaxation factor + const double possible_updated_value = old_value - dx; - if (!relaxation_avoid_negative && !relaxation_avoid_over_one) { - return relaxation_factor; // 1.0 - } - } - - // we will not be this step if dx == 0. - assert(dx != 0.0); - - if (relaxation_avoid_negative) { + // 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 (relaxation_avoid_over_one) { + } 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; } @@ -2276,16 +2267,12 @@ namespace Opm if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { const double relaxation_factor_w = relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]); - if (relaxation_factor_w < relaxation_factor) { - relaxation_factor = relaxation_factor_w; - } + 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]); - if (relaxation_factor_g < relaxation_factor) { - relaxation_factor = relaxation_factor_g; - } + relaxation_factor = std::min(relaxation_factor, relaxation_factor_g); } if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { @@ -2303,7 +2290,6 @@ namespace Opm } } - // relaxation factor for the total rates { const double original_total_rate = primary_variables[WQTotal]; From 5e736fc4c78119f9956072bbb0663912068010de Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 8 Nov 2018 13:37:17 +0100 Subject: [PATCH 3/3] decoupe the rate and fraction relaxation the relaxation for these two have different goals and problems, it is wise to decuple them. --- opm/autodiff/StandardWell.hpp | 10 +++---- opm/autodiff/StandardWell_impl.hpp | 42 ++++++++++++++---------------- 2 files changed, 25 insertions(+), 27 deletions(-) diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index 4f4369050..234dd6a4e 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -358,14 +358,14 @@ namespace Opm static double relaxationFactorFraction(const double old_value, const double dx); - // calculate a relaxation factor to avoid overshoot + // calculate a relaxation factor to avoid overshoot of the fractions for producers // which might result in negative rates - static double determineRelaxationFactorProducer(const std::vector& primary_variables, + static double relaxationFactorFractionsProducer(const std::vector& primary_variables, const BVectorWell& dwells); - // calculate a relaxation factor to avoid overshoot for injectors - static double determineRelaxationFactorInjector(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 db768ae6e..aa25485e8 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -870,34 +870,37 @@ namespace Opm const std::vector old_primary_variables = primary_variables_; - const double relaxation_factor = (well_type_ == PRODUCER) ? - determineRelaxationFactorProducer(old_primary_variables, dwells) - : determineRelaxationFactorInjector(old_primary_variables, dwells); + // 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] * relaxation_factor), 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] * relaxation_factor), 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]) * relaxation_factor, 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] * relaxation_factor; + // 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 { @@ -2224,6 +2227,7 @@ namespace Opm + template double StandardWell:: @@ -2258,7 +2262,7 @@ namespace Opm template double StandardWell:: - determineRelaxationFactorProducer(const std::vector& primary_variables, + relaxationFactorFractionsProducer(const std::vector& primary_variables, const BVectorWell& dwells) { // TODO: not considering solvent yet @@ -2290,17 +2294,7 @@ namespace Opm } } - // relaxation factor for the total rates - { - const double original_total_rate = primary_variables[WQTotal]; - const double relaxed_update = dwells[0][WQTotal] * relaxation_factor; - const double possible_update_total_rate = primary_variables[WQTotal] - relaxed_update; - - if (original_total_rate * possible_update_total_rate < 0.) { // sign changed - const double further_relaxation_factor = std::abs(original_total_rate / relaxed_update) * 0.95; - relaxation_factor *= further_relaxation_factor; - } - } + assert(relaxation_factor >= 0.0 && relaxation_factor <= 1.0); return relaxation_factor; } @@ -2312,8 +2306,8 @@ namespace Opm template double StandardWell:: - determineRelaxationFactorInjector(const std::vector& primary_variables, - const BVectorWell& dwells) + relaxationFactorRate(const std::vector& primary_variables, + const BVectorWell& dwells) { double relaxation_factor = 1.0; @@ -2323,10 +2317,14 @@ namespace Opm 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; } }