OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

improving the ssi+newton PTFlash method

Browse Source
This commit is contained in:
Kai Bao 2023-11-01 22:08:51 +01:00
parent aacf4a4d81
commit e3362f3b01
1 changed files with 29 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

53
opm/material/constraintsolvers/PTFlash.hpp
View File

@ -307,9 +307,7 @@ protected:
// Check if Lmin < Lmax, and switch if not
if (Lmin > Lmax) {
auto Ltmp = Lmin;
Lmin = Lmax;
Lmax = Ltmp;
std::swap(Lmin, Lmax);
}
Lmin = std::clamp(Lmin, 0., 1.);
@ -340,8 +338,7 @@ protected:
L -= delta;
// Check if L is within the bounds, and if not, we apply bisection method
if (L < Lmin || L > Lmax)
{
if (L < Lmin || L > Lmax) {
// Print info
if (verbosity == 3 || verbosity == 4) {
std::cout << "L is not within the the range [Lmin, Lmax], solve using Bisection method!" << std::endl;
@ -358,7 +355,7 @@ protected:
std::cout << "Rachford-Rice (Bisection) converged to final solution L = " << L << std::endl;
}
return L;
}
}
// Print iteration info
if (verbosity == 3 || verbosity == 4) {
@ -607,27 +604,34 @@ protected:
// Calculate composition using nonlinear solver
// Newton
bool converged = false;
if (flash_2p_method == "newton") {
if (verbosity >= 1) {
std::cout << "Calculate composition using Newton." << std::endl;
}
newtonComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, verbosity);
converged = newtonComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, verbosity);
} else if (flash_2p_method == "ssi") {
// Successive substitution
if (verbosity >= 1) {
std::cout << "Calculate composition using Succcessive Substitution." << std::endl;
}
successiveSubstitutionComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, false, verbosity);
converged = successiveSubstitutionComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, false, verbosity);
} else if (flash_2p_method == "ssi+newton") {
successiveSubstitutionComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, true, verbosity);
newtonComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, verbosity);
converged = successiveSubstitutionComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, true, verbosity);
if (!converged) {
converged = newtonComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, verbosity);
}
} else {
throw std::runtime_error("unknown two phase flash method " + flash_2p_method + " is specified");
throw std::logic_error("unknown two phase flash method " + flash_2p_method + " is specified");
}
if (!converged) {
throw std::runtime_error("flash calculation did not get converged with " + flash_2p_method);
}
}
template <class FlashFluidState, class ComponentVector>
static void newtonComposition_(ComponentVector& K, typename FlashFluidState::Scalar& L,
static bool newtonComposition_(ComponentVector& K, typename FlashFluidState::Scalar& L,
FlashFluidState& fluid_state, const ComponentVector& z,
int verbosity)
{
@ -785,6 +789,7 @@ protected:
}
L = Opm::getValue(l);
fluid_state.setLvalue(L);
return converged;
}
// TODO: the interface will need to refactor for later usage
@ -1126,11 +1131,11 @@ protected:
// TODO: or use typename FlashFluidState::Scalar
template <class FlashFluidState, class ComponentVector>
static void successiveSubstitutionComposition_(ComponentVector& K, typename ComponentVector::field_type& L, FlashFluidState& fluid_state, const ComponentVector& z,
static bool successiveSubstitutionComposition_(ComponentVector& K, typename ComponentVector::field_type& L, FlashFluidState& fluid_state, const ComponentVector& z,
const bool newton_afterwards, const int verbosity)
{
// Determine max. iterations based on if it will be used as a standalone flash or as a pre-process to Newton (or other) method.
const int maxIterations = newton_afterwards ? 3 : 100;
const int maxIterations = newton_afterwards ? 5 : 100;
// Store cout format before manipulation
std::ios_base::fmtflags f(std::cout.flags());
@ -1187,7 +1192,7 @@ protected:
}
// Check convergence
if (convFugRatio.two_norm() < 1e-6) {
if (convFugRatio.two_norm() < 1.e-6) {
// Restore cout format
std::cout.flags(f);
@ -1195,7 +1200,7 @@ protected:
if (verbosity >= 1) {
std::cout << "Solution converged to the following result :" << std::endl;
std::cout << "x = [";
for (int compIdx=0; compIdx<numComponents; ++compIdx){
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
if (compIdx < numComponents - 1)
std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx) << " ";
else
@ -1203,7 +1208,7 @@ protected:
}
std::cout << "]" << std::endl;
std::cout << "y = [";
for (int compIdx=0; compIdx<numComponents; ++compIdx){
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
if (compIdx < numComponents - 1)
std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx) << " ";
else
@ -1214,11 +1219,10 @@ protected:
std::cout << "L = " << L << std::endl;
}
// Restore cout format format
return;
return true;
}
// If convergence is not met, K is updated in a successive substitution manner
else{
else {
// Update K
for (int compIdx=0; compIdx<numComponents; ++compIdx){
K[compIdx] *= newFugRatio[compIdx];
@ -1227,12 +1231,13 @@ protected:
// Solve Rachford-Rice to get L from updated K
L = solveRachfordRice_g_(K, z, 0);
}
}
if (!newton_afterwards) {
throw std::runtime_error(
"Successive substitution composition update did not converge within maxIterations " + std::to_string(maxIterations));
// did not get converged. check whether will do more newton later afterward
if (newton_afterwards && verbosity > 0) {
std::cout << "Successive substitution composition update did not converge within maxIterations "
<< maxIterations << std::endl;
}
return false;
}
};//end PTFlash