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

Merge e3362f3b01 into 9a0ccfd853

Browse Source
This commit is contained in:
Kai Bao 2023-11-26 11:16:36 +00:00 committed by GitHub
commit a2d0f38812
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
1 changed files with 40 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -297,10 +297,8 @@ protected:
typename Vector::field_type Kmin = K[0]; typename Vector::field_type Kmin = K[0];
typename Vector::field_type Kmax = K[0]; typename Vector::field_type Kmax = K[0];
for (int compIdx=1; compIdx<numComponents; ++compIdx){ for (int compIdx=1; compIdx<numComponents; ++compIdx){
if (K[compIdx] < Kmin) Kmin = std::min(Kmin, K[compIdx]);
Kmin = K[compIdx]; Kmax = std::max(Kmin, K[compIdx]);
else if (K[compIdx] >= Kmax)
Kmax = K[compIdx];
} }
// Lower and upper bound for solution // Lower and upper bound for solution
@ -308,11 +306,16 @@ protected:
auto Lmax = Kmax / (Kmax - 1); auto Lmax = Kmax / (Kmax - 1);
// Check if Lmin < Lmax, and switch if not // Check if Lmin < Lmax, and switch if not
if (Lmin > Lmax) if (Lmin > Lmax) {
{ std::swap(Lmin, Lmax);
auto Ltmp = Lmin; }
Lmin = Lmax;
Lmax = Ltmp; Lmin = std::clamp(Lmin, 0., 1.);
Lmax = std::clamp(Lmax, 0., 1.);
if (Lmin == Lmax) {
Lmin = 0.;
Lmax = 1.0;
} }
// Initial guess // Initial guess
@ -335,8 +338,7 @@ protected:
L -= delta; L -= delta;
// Check if L is within the bounds, and if not, we apply bisection method // 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 // Print info
if (verbosity == 3 || verbosity == 4) { if (verbosity == 3 || verbosity == 4) {
std::cout << "L is not within the the range [Lmin, Lmax], solve using Bisection method!" << std::endl; std::cout << "L is not within the the range [Lmin, Lmax], solve using Bisection method!" << std::endl;
@ -353,7 +355,7 @@ protected:
std::cout << "Rachford-Rice (Bisection) converged to final solution L = " << L << std::endl; std::cout << "Rachford-Rice (Bisection) converged to final solution L = " << L << std::endl;
} }
return L; return L;
} }
// Print iteration info // Print iteration info
if (verbosity == 3 || verbosity == 4) { if (verbosity == 3 || verbosity == 4) {
@ -602,27 +604,34 @@ protected:
// Calculate composition using nonlinear solver // Calculate composition using nonlinear solver
// Newton // Newton
bool converged = false;
if (flash_2p_method == "newton") { if (flash_2p_method == "newton") {
if (verbosity >= 1) { if (verbosity >= 1) {
std::cout << "Calculate composition using Newton." << std::endl; 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") { } else if (flash_2p_method == "ssi") {
// Successive substitution // Successive substitution
if (verbosity >= 1) { if (verbosity >= 1) {
std::cout << "Calculate composition using Succcessive Substitution." << std::endl; 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") { } else if (flash_2p_method == "ssi+newton") {
successiveSubstitutionComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, true, verbosity); converged = successiveSubstitutionComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, true, verbosity);
newtonComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, verbosity); if (!converged) {
converged = newtonComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, verbosity);
}
} else { } 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> 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, FlashFluidState& fluid_state, const ComponentVector& z,
int verbosity) int verbosity)
{ {
@ -780,6 +789,7 @@ protected:
} }
L = Opm::getValue(l); L = Opm::getValue(l);
fluid_state.setLvalue(L); fluid_state.setLvalue(L);
return converged;
} }
// TODO: the interface will need to refactor for later usage // TODO: the interface will need to refactor for later usage
@ -1121,11 +1131,11 @@ protected:
// TODO: or use typename FlashFluidState::Scalar // TODO: or use typename FlashFluidState::Scalar
template <class FlashFluidState, class ComponentVector> 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) 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. // 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 // Store cout format before manipulation
std::ios_base::fmtflags f(std::cout.flags()); std::ios_base::fmtflags f(std::cout.flags());
@ -1182,7 +1192,7 @@ protected:
} }
// Check convergence // Check convergence
if (convFugRatio.two_norm() < 1e-6) { if (convFugRatio.two_norm() < 1.e-6) {
// Restore cout format // Restore cout format
std::cout.flags(f); std::cout.flags(f);
@ -1190,7 +1200,7 @@ protected:
if (verbosity >= 1) { if (verbosity >= 1) {
std::cout << "Solution converged to the following result :" << std::endl; std::cout << "Solution converged to the following result :" << std::endl;
std::cout << "x = ["; std::cout << "x = [";
for (int compIdx=0; compIdx<numComponents; ++compIdx){ for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
if (compIdx < numComponents - 1) if (compIdx < numComponents - 1)
std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx) << " "; std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx) << " ";
else else
@ -1198,7 +1208,7 @@ protected:
} }
std::cout << "]" << std::endl; std::cout << "]" << std::endl;
std::cout << "y = ["; std::cout << "y = [";
for (int compIdx=0; compIdx<numComponents; ++compIdx){ for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
if (compIdx < numComponents - 1) if (compIdx < numComponents - 1)
std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx) << " "; std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx) << " ";
else else
@ -1209,11 +1219,10 @@ protected:
std::cout << "L = " << L << std::endl; std::cout << "L = " << L << std::endl;
} }
// Restore cout format format // Restore cout format format
return; return true;
} }
// If convergence is not met, K is updated in a successive substitution manner // If convergence is not met, K is updated in a successive substitution manner
else{ else {
// Update K // Update K
for (int compIdx=0; compIdx<numComponents; ++compIdx){ for (int compIdx=0; compIdx<numComponents; ++compIdx){
K[compIdx] *= newFugRatio[compIdx]; K[compIdx] *= newFugRatio[compIdx];
@ -1222,12 +1231,13 @@ protected:
// Solve Rachford-Rice to get L from updated K // Solve Rachford-Rice to get L from updated K
L = solveRachfordRice_g_(K, z, 0); L = solveRachfordRice_g_(K, z, 0);
} }
} }
if (!newton_afterwards) { // did not get converged. check whether will do more newton later afterward
throw std::runtime_error( if (newton_afterwards && verbosity > 0) {
"Successive substitution composition update did not converge within maxIterations " + std::to_string(maxIterations)); std::cout << "Successive substitution composition update did not converge within maxIterations "
<< maxIterations << std::endl;
} }
return false;
} }
};//end PTFlash };//end PTFlash