updateNewton: avoid copying

This commit is contained in:
Arne Morten Kvarving
2022-11-09 11:10:59 +01:00
parent e5e6e21122
commit bc3680f822

View File

@@ -244,12 +244,13 @@ updateNewton(const BVectorWell& dwells,
[[maybe_unused]] const double dFLimit,
const double dBHPLimit)
{
const std::vector<double> old_primary_variables = value_;
const double relaxation_factor_rate = relaxationFactorRate(value_[WQTotal],
dwells[0][WQTotal]);
// 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 =
well_.isProducer() ? relaxationFactorFractionsProducer(old_primary_variables, dwells)
well_.isProducer() ? relaxationFactorFractionsProducer(value_, dwells)
: 1.0;
// update the second and third well variable (The flux fractions)
@@ -257,35 +258,32 @@ updateNewton(const BVectorWell& dwells,
if constexpr (has_wfrac_variable) {
const int sign2 = dwells[0][WFrac] > 0 ? 1: -1;
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;
value_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
value_[WFrac] = value_[WFrac] - dx2_limited;
}
if constexpr (has_gfrac_variable) {
const int sign3 = dwells[0][GFrac] > 0 ? 1: -1;
const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac] * relaxation_factor_fractions), dFLimit);
value_[GFrac] = old_primary_variables[GFrac] - dx3_limited;
value_[GFrac] = value_[GFrac] - dx3_limited;
}
if constexpr (Indices::enableSolvent) {
const int sign4 = dwells[0][SFrac] > 0 ? 1: -1;
const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]) * relaxation_factor_fractions, dFLimit);
value_[SFrac] = old_primary_variables[SFrac] - dx4_limited;
value_[SFrac] = value_[SFrac] - dx4_limited;
}
this->processFractions();
// updating the total rates Q_t
const double relaxation_factor_rate = relaxationFactorRate(old_primary_variables[WQTotal],
dwells[0][WQTotal]);
value_[WQTotal] = old_primary_variables[WQTotal] - dwells[0][WQTotal] * relaxation_factor_rate;
value_[WQTotal] = value_[WQTotal] - dwells[0][WQTotal] * relaxation_factor_rate;
// updating the bottom hole pressure
const int sign1 = dwells[0][Bhp] > 0 ? 1: -1;
const double dx1_limited = sign1 * std::min(std::abs(dwells[0][Bhp]),
std::abs(old_primary_variables[Bhp]) * dBHPLimit);
std::abs(value_[Bhp]) * dBHPLimit);
// 1e5 to make sure bhp will not be below 1bar
value_[Bhp] = std::max(old_primary_variables[Bhp] - dx1_limited, 1e5);
value_[Bhp] = std::max(value_[Bhp] - dx1_limited, 1e5);
}
template<class FluidSystem, class Indices, class Scalar>