diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index 2536d97a7..3f2bec31f 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -122,7 +122,6 @@ namespace Opm using Eval = typename StdWellEval::Eval; using EvalWell = typename StdWellEval::EvalWell; using BVectorWell = typename StdWellEval::BVectorWell; - using DiagMatrixBlockWellType = typename StdWellEval::DiagMatrixBlockWellType; StandardWell(const Well& well, const ParallelWellInfo& pw_info, diff --git a/opm/simulators/wells/StandardWellEval.cpp b/opm/simulators/wells/StandardWellEval.cpp index cdfbd9648..4ee80ccb8 100644 --- a/opm/simulators/wells/StandardWellEval.cpp +++ b/opm/simulators/wells/StandardWellEval.cpp @@ -69,53 +69,6 @@ extendEval(const Eval& in) const return out; } -template -double -StandardWellEval:: -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::numActivePhases() > 1) { - if constexpr (has_wfrac_variable) { - const double relaxation_factor_w = StandardWellGeneric:: - relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]); - relaxation_factor = std::min(relaxation_factor, relaxation_factor_w); - } - - if constexpr (has_gfrac_variable) { - const double relaxation_factor_g = StandardWellGeneric:: - relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]); - relaxation_factor = std::min(relaxation_factor, relaxation_factor_g); - } - - - if constexpr (has_wfrac_variable && has_gfrac_variable) { - // 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; - // We only relax if fraction is above 1. - // The newton solver should handle the rest - const double epsilon = 0.001; - if (possible_updated_sum > 1.0 + epsilon) { - // since the orignal sum <= 1.0 the epsilon asserts that - // the relaxed_update is non trivial. - 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 void StandardWellEval:: @@ -166,7 +119,7 @@ updatePrimaryVariablesNewton(const BVectorWell& 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 [[maybe_unused]] const double relaxation_factor_fractions = - (baseif_.isProducer()) ? relaxationFactorFractionsProducer(old_primary_variables, dwells) + (baseif_.isProducer()) ? this->primary_variables_.relaxationFactorFractionsProducer(old_primary_variables, dwells) : 1.0; // update the second and third well variable (The flux fractions) diff --git a/opm/simulators/wells/StandardWellEval.hpp b/opm/simulators/wells/StandardWellEval.hpp index a78bd94c4..66940bd0b 100644 --- a/opm/simulators/wells/StandardWellEval.hpp +++ b/opm/simulators/wells/StandardWellEval.hpp @@ -62,7 +62,7 @@ protected: public: using EvalWell = typename PrimaryVariables::EvalWell; using Eval = DenseAd::Evaluation; - using BVectorWell = typename StandardWellGeneric::BVectorWell; + using BVectorWell = typename StandardWellEquations::BVectorWell; //! \brief Returns a const reference to equation system. const StandardWellEquations& linSys() const @@ -85,11 +85,6 @@ protected: EvalWell extendEval(const Eval& in) const; - // 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); - // computing the accumulation term for later use in well mass equations void computeAccumWell(); diff --git a/opm/simulators/wells/StandardWellGeneric.hpp b/opm/simulators/wells/StandardWellGeneric.hpp index 53247e717..8d7133106 100644 --- a/opm/simulators/wells/StandardWellGeneric.hpp +++ b/opm/simulators/wells/StandardWellGeneric.hpp @@ -48,25 +48,7 @@ class WellState; template class StandardWellGeneric { -protected: - // sparsity pattern for the matrices - //[A C^T [x = [ res - // B D ] x_well] res_well] - - // the vector type for the res_well and x_well - using VectorBlockWellType = Dune::DynamicVector; - using BVectorWell = Dune::BlockVector; - - // the matrix type for the diagonal matrix D - using DiagMatrixBlockWellType = Dune::DynamicMatrix; - using DiagMatWell = Dune::BCRSMatrix; - - // the matrix type for the non-diagonal matrix B and C^T - using OffDiagMatrixBlockWellType = Dune::DynamicMatrix; - using OffDiagMatWell = Dune::BCRSMatrix; - - StandardWellGeneric(const WellInterfaceGeneric& baseif); - +public: // calculate a relaxation factor to avoid overshoot of total rates static double relaxationFactorRate(const std::vector& primary_variables, const double newton_update); @@ -75,6 +57,9 @@ protected: static double relaxationFactorFraction(const double old_value, const double dx); +protected: + StandardWellGeneric(const WellInterfaceGeneric& baseif); + void computeConnectionPressureDelta(); // Base interface reference diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.cpp b/opm/simulators/wells/StandardWellPrimaryVariables.cpp index afa4324b6..8bab865e9 100644 --- a/opm/simulators/wells/StandardWellPrimaryVariables.cpp +++ b/opm/simulators/wells/StandardWellPrimaryVariables.cpp @@ -34,9 +34,12 @@ #include +#include #include #include +#include + namespace Opm { template @@ -524,6 +527,52 @@ processFractions() } } +template +double StandardWellPrimaryVariables:: +relaxationFactorFractionsProducer(const std::vector& primary_variables, + const BVectorWell& dwells) const +{ + // TODO: not considering solvent yet + // 0.95 is a experimental value, which remains to be optimized + double relaxation_factor = 1.0; + + if (FluidSystem::numActivePhases() > 1) { + if constexpr (has_wfrac_variable) { + const double relaxation_factor_w = StandardWellGeneric:: + relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]); + relaxation_factor = std::min(relaxation_factor, relaxation_factor_w); + } + + if constexpr (has_gfrac_variable) { + const double relaxation_factor_g = StandardWellGeneric:: + relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]); + relaxation_factor = std::min(relaxation_factor, relaxation_factor_g); + } + + + if constexpr (has_wfrac_variable && has_gfrac_variable) { + // 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; + // We only relax if fraction is above 1. + // The newton solver should handle the rest + const double epsilon = 0.001; + if (possible_updated_sum > 1.0 + epsilon) { + // since the orignal sum <= 1.0 the epsilon asserts that + // the relaxed_update is non trivial. + 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; +} + #define INSTANCE(...) \ template class StandardWellPrimaryVariables,__VA_ARGS__,double>; diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.hpp b/opm/simulators/wells/StandardWellPrimaryVariables.hpp index b1e3da68f..9b0a9d22e 100644 --- a/opm/simulators/wells/StandardWellPrimaryVariables.hpp +++ b/opm/simulators/wells/StandardWellPrimaryVariables.hpp @@ -132,6 +132,11 @@ public: //! \brief Handle non-reasonable fractions due to numerical overshoot. void processFractions(); + //! \brief Calculate a relaxation factor for producers. + //! \details To avoid overshoot of the fractions which might result in negative rates. + double relaxationFactorFractionsProducer(const std::vector& primary_variables, + const BVectorWell& dwells) const; + private: //! \brief Returns volume fraction for a component. EvalWell volumeFraction(const unsigned compIdx) const;