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

some refactoring for the Newton flash solver.

Browse Source
This commit is contained in:
Kai Bao 2021-12-22 11:30:04 +01:00 committed by Trine Mykkeltvedt
parent 6066a8dfa7
commit 3368256ecc
2 changed files with 98 additions and 67 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

162
opm/material/constraintsolvers/ChiFlash.hpp
View File

@ -86,8 +86,8 @@ public:
* *
*/ */
template <class FluidState> template <class FluidState>
static void solve(FluidState& fluidState, static void solve(FluidState& fluid_state,
const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalComposition, const Dune::FieldVector<typename FluidState::Scalar, numComponents>& z,
int spatialIdx, int spatialIdx,
int verbosity, int verbosity,
std::string twoPhaseMethod, std::string twoPhaseMethod,
@ -107,113 +107,146 @@ public:
//K and L from previous timestep (wilson and -1 initially) //K and L from previous timestep (wilson and -1 initially)
ComponentVector K; ComponentVector K;
for(int compIdx = 0; compIdx < numComponents; ++compIdx) { for(int compIdx = 0; compIdx < numComponents; ++compIdx) {
K[compIdx] = fluidState.K(compIdx); K[compIdx] = fluid_state.K(compIdx);
} }
InputEval L; InputEval L;
L = fluidState.L(0); L = fluid_state.L(0);
// Print header // Print header
if (verbosity >= 1) { if (verbosity >= 1) {
std::cout << "********" << std::endl; std::cout << "********" << std::endl;
std::cout << "Flash calculations on Cell " << spatialIdx << std::endl; std::cout << "Flash calculations on Cell " << spatialIdx << std::endl;
std::cout << "Inputs are K = [" << K << "], L = [" << L << "], z = [" << globalComposition << "], P = " << fluidState.pressure(0) << ", and T = " << fluidState.temperature(0) << std::endl; std::cout << "Inputs are K = [" << K << "], L = [" << L << "], z = [" << z << "], P = " << fluid_state.pressure(0) << ", and T = " << fluid_state.temperature(0) << std::endl;
} }
using ScalarVector = Dune::FieldVector<Scalar, numComponents>;
Scalar L_scalar = Opm::getValue(L);
ScalarVector z_scalar;
ScalarVector K_scalar;
for (unsigned i = 0; i < numComponents; ++i) {
z_scalar[i] = Opm::getValue(z[i]);
K_scalar[i] = Opm::getValue(K[i]);
}
using ScalarFluidState = CompositionalFluidState<Scalar, FluidSystem>;
ScalarFluidState fluid_state_scalar;
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
fluid_state_scalar.setMoleFraction(oilPhaseIdx, compIdx, Opm::getValue(fluid_state.moleFraction(oilPhaseIdx, compIdx)));
fluid_state_scalar.setMoleFraction(gasPhaseIdx, compIdx, Opm::getValue(fluid_state.moleFraction(gasPhaseIdx, compIdx)));
fluid_state_scalar.setKvalue(compIdx, Opm::getValue(fluid_state.K(compIdx)));
}
fluid_state_scalar.setLvalue(L_scalar);
// other values need to be Scalar, but I guess the fluidstate does not support it yet.
fluid_state_scalar.setPressure(FluidSystem::oilPhaseIdx,
Opm::getValue(fluid_state.pressure(FluidSystem::oilPhaseIdx)));
fluid_state_scalar.setPressure(FluidSystem::gasPhaseIdx,
Opm::getValue(fluid_state.pressure(FluidSystem::gasPhaseIdx)));
// TODO: not sure whether we need to set the saturations
fluid_state_scalar.setSaturation(FluidSystem::gasPhaseIdx,
Opm::getValue(fluid_state_scalar.saturation(FluidSystem::gasPhaseIdx)));
fluid_state_scalar.setSaturation(FluidSystem::oilPhaseIdx,
Opm::getValue(fluid_state_scalar.saturation(FluidSystem::oilPhaseIdx)));
fluid_state_scalar.setTemperature(Opm::getValue(fluid_state.temperature(0)));
// Rachford Rice equation to get initial L for composition solver
L_scalar = solveRachfordRice_g_(K_scalar, z_scalar, verbosity);
// Do a stability test to check if cell is single-phase (do for all cells the first time). // Do a stability test to check if cell is single-phase (do for all cells the first time).
bool isStable = false; bool isStable = false;
if ( L <= 0 || L == 1 ) { if ( L <= 0 || L == 1 ) {
if (verbosity >= 1) { if (verbosity >= 1) {
std::cout << "Perform stability test (L <= 0 or L == 1)!" << std::endl; std::cout << "Perform stability test (L <= 0 or L == 1)!" << std::endl;
} }
phaseStabilityTest_(isStable, K, fluidState, globalComposition, verbosity); phaseStabilityTest_(isStable, K_scalar, fluid_state_scalar, z_scalar, verbosity);
} }
if (verbosity >= 1) { if (verbosity >= 1) {
std::cout << "Inputs after stability test are K = [" << K << "], L = [" << L << "], z = [" << globalComposition << "], P = " << fluidState.pressure(0) << ", and T = " << fluidState.temperature(0) << std::endl; std::cout << "Inputs after stability test are K = [" << K_scalar << "], L = [" << L_scalar << "], z = [" << z_scalar << "], P = " << fluid_state.pressure(0) << ", and T = " << fluid_state.temperature(0) << std::endl;
} }
// Update the composition if cell is two-phase // Update the composition if cell is two-phase
if ( !isStable) { if ( !isStable) {
flash_2ph(globalComposition, twoPhaseMethod, K, L, fluidState, verbosity); flash_2ph(z_scalar, twoPhaseMethod, K_scalar, L_scalar, fluid_state_scalar, verbosity);
} else { } else {
// Cell is one-phase. Use Li's phase labeling method to see if it's liquid or vapor // Cell is one-phase. Use Li's phase labeling method to see if it's liquid or vapor
L = li_single_phase_label_(fluidState, globalComposition, verbosity); L = li_single_phase_label_(fluid_state, z, verbosity);
} }
newtonCompositionUpdate_(K, L, fluidState, globalComposition, verbosity);
// Print footer // Print footer
if (verbosity >= 1) { if (verbosity >= 1) {
std::cout << "********" << std::endl; std::cout << "********" << std::endl;
} }
// Update phases // Update phases
typename FluidSystem::template ParameterCache<InputEval> paramCache; /* typename FluidSystem::template ParameterCache<InputEval> paramCache;
paramCache.updatePhase(fluidState, oilPhaseIdx); paramCache.updatePhase(fluid_state, oilPhaseIdx);
paramCache.updatePhase(fluidState, gasPhaseIdx); paramCache.updatePhase(fluid_state, gasPhaseIdx); */
// Calculate compressibility factor /* // Calculate compressibility factor
const Scalar R = Opm::Constants<Scalar>::R; const Scalar R = Opm::Constants<Scalar>::R;
InputEval Z_L = (paramCache.molarVolume(oilPhaseIdx) * fluidState.pressure(oilPhaseIdx) )/ InputEval Z_L = (paramCache.molarVolume(oilPhaseIdx) * fluid_state.pressure(oilPhaseIdx) ) /
(R * fluidState.temperature(oilPhaseIdx)); (R * fluid_state.temperature(oilPhaseIdx));
InputEval Z_V = (paramCache.molarVolume(gasPhaseIdx) * fluidState.pressure(gasPhaseIdx) )/ InputEval Z_V = (paramCache.molarVolume(gasPhaseIdx) * fluid_state.pressure(gasPhaseIdx) ) /
(R * fluidState.temperature(gasPhaseIdx)); (R * fluid_state.temperature(gasPhaseIdx));
std::cout << " the type of InputEval here is " << Dune::className<InputEval>() << std::endl; std::cout << " the type of InputEval here is " << Dune::className<InputEval>() << std::endl;
// Update saturation // Update saturation
InputEval So = L*Z_L/(L*Z_L+(1-L)*Z_V); InputEval So = L*Z_L/(L*Z_L+(1-L)*Z_V);
InputEval Sg = 1-So; InputEval Sg = 1-So;
fluidState.setSaturation(oilPhaseIdx, So); fluid_state.setSaturation(oilPhaseIdx, So);
fluidState.setSaturation(gasPhaseIdx, Sg); fluid_state.setSaturation(gasPhaseIdx, Sg);
//Update L and K to the problem for the next flash //Update L and K to the problem for the next flash
for (int compIdx = 0; compIdx < numComponents; ++compIdx){ for (int compIdx = 0; compIdx < numComponents; ++compIdx){
fluidState.setKvalue(compIdx, K[compIdx]); fluid_state.setKvalue(compIdx, K[compIdx]);
} }
fluidState.setCompressFactor(oilPhaseIdx, Z_L); fluid_state.setCompressFactor(oilPhaseIdx, Z_L);
fluidState.setCompressFactor(gasPhaseIdx, Z_V); fluid_state.setCompressFactor(gasPhaseIdx, Z_V);
fluidState.setLvalue(L); fluid_state.setLvalue(L); */
// Print saturation // Print saturation
/* std::cout << " output the molefraction derivatives" << std::endl; /* std::cout << " output the molefraction derivatives" << std::endl;
std::cout << " for vapor comp 1 " << std::endl; std::cout << " for vapor comp 1 " << std::endl;
fluidState.moleFraction(gasPhaseIdx, 0).print(); fluid_state.moleFraction(gasPhaseIdx, 0).print();
std::cout << std::endl << " for vapor comp 2 " << std::endl; std::cout << std::endl << " for vapor comp 2 " << std::endl;
fluidState.moleFraction(gasPhaseIdx, 1).print(); fluid_state.moleFraction(gasPhaseIdx, 1).print();
std::cout << std::endl << " for vapor comp 3 " << std::endl; std::cout << std::endl << " for vapor comp 3 " << std::endl;
fluidState.moleFraction(gasPhaseIdx, 2).print(); fluid_state.moleFraction(gasPhaseIdx, 2).print();
std::cout << std::endl; std::cout << std::endl;
std::cout << " for oil comp 1 " << std::endl; std::cout << " for oil comp 1 " << std::endl;
fluidState.moleFraction(oilPhaseIdx, 0).print(); fluid_state.moleFraction(oilPhaseIdx, 0).print();
std::cout << std::endl << " for oil comp 2 " << std::endl; std::cout << std::endl << " for oil comp 2 " << std::endl;
fluidState.moleFraction(oilPhaseIdx, 1).print(); fluid_state.moleFraction(oilPhaseIdx, 1).print();
std::cout << std::endl << " for oil comp 3 " << std::endl; std::cout << std::endl << " for oil comp 3 " << std::endl;
fluidState.moleFraction(oilPhaseIdx, 2).print(); fluid_state.moleFraction(oilPhaseIdx, 2).print();
std::cout << std::endl; std::cout << std::endl;
std::cout << " for pressure " << std::endl; std::cout << " for pressure " << std::endl;
fluidState.pressure(0).print(); fluid_state.pressure(0).print();
std::cout<< std::endl; std::cout<< std::endl;
fluidState.pressure(1).print(); fluid_state.pressure(1).print();
std::cout<< std::endl; */ std::cout<< std::endl; */
// Update densities // Update densities
fluidState.setDensity(oilPhaseIdx, FluidSystem::density(fluidState, paramCache, oilPhaseIdx)); // fluid_state.setDensity(oilPhaseIdx, FluidSystem::density(fluid_state, paramCache, oilPhaseIdx));
fluidState.setDensity(gasPhaseIdx, FluidSystem::density(fluidState, paramCache, gasPhaseIdx)); // fluid_state.setDensity(gasPhaseIdx, FluidSystem::density(fluid_state, paramCache, gasPhaseIdx));
// check the residuals of the equations // check the residuals of the equations
using NewtonVector = Dune::FieldVector<InputEval, numMisciblePhases*numMiscibleComponents+1>; /* using NewtonVector = Dune::FieldVector<InputEval, numMisciblePhases*numMiscibleComponents+1>;
NewtonVector newtonX; NewtonVector newtonX;
NewtonVector newtonB; NewtonVector newtonB;
for (int compIdx=0; compIdx<numComponents; ++compIdx){ for (int compIdx=0; compIdx<numComponents; ++compIdx){
newtonX[compIdx] = Opm::getValue(fluidState.moleFraction(oilPhaseIdx, compIdx)); newtonX[compIdx] = Opm::getValue(fluid_state.moleFraction(oilPhaseIdx, compIdx));
newtonX[compIdx + numMiscibleComponents] = Opm::getValue(fluidState.moleFraction(gasPhaseIdx, compIdx)); newtonX[compIdx + numMiscibleComponents] = Opm::getValue(fluid_state.moleFraction(gasPhaseIdx, compIdx));
} }
newtonX[numMisciblePhases*numMiscibleComponents] = Opm::getValue(L); newtonX[numMisciblePhases*numMiscibleComponents] = Opm::getValue(L);
evalDefect_(newtonB, newtonX, fluidState, globalComposition); evalDefect_(newtonB, newtonX, fluid_state, z); */
std::cout << " the residuals of the equations is " << std::endl; /* std::cout << " the residuals of the equations is " << std::endl;
for (unsigned i = 0; i < newtonB.N(); ++i) { for (unsigned i = 0; i < newtonB.N(); ++i) {
std::cout << newtonB[i] << std::endl; std::cout << newtonB[i] << std::endl;
} }
@ -221,11 +254,11 @@ public:
if (verbosity >= 5) { if (verbosity >= 5) {
std::cout << " mole fractions for oil " << std::endl; std::cout << " mole fractions for oil " << std::endl;
for (int i = 0; i < numComponents; ++i) { for (int i = 0; i < numComponents; ++i) {
std::cout << " i " << i << " " << fluidState.moleFraction(oilPhaseIdx, i) << std::endl; std::cout << " i " << i << " " << fluid_state.moleFraction(oilPhaseIdx, i) << std::endl;
} }
std::cout << " mole fractions for gas " << std::endl; std::cout << " mole fractions for gas " << std::endl;
for (int i = 0; i < numComponents; ++i) { for (int i = 0; i < numComponents; ++i) {
std::cout << " i " << i << " " << fluidState.moleFraction(gasPhaseIdx, i) << std::endl; std::cout << " i " << i << " " << fluid_state.moleFraction(gasPhaseIdx, i) << std::endl;
} }
std::cout << " K " << std::endl; std::cout << " K " << std::endl;
for (int i = 0; i < numComponents; ++i) { for (int i = 0; i < numComponents; ++i) {
@ -235,25 +268,25 @@ public:
std::cout << "Z_V = " << Z_V <<std::endl; std::cout << "Z_V = " << Z_V <<std::endl;
std::cout << " fugacity for oil phase " << std::endl; std::cout << " fugacity for oil phase " << std::endl;
for (int i = 0; i < numComponents; ++i) { for (int i = 0; i < numComponents; ++i) {
auto phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, oilPhaseIdx, i); auto phi = FluidSystem::fugacityCoefficient(fluid_state, paramCache, oilPhaseIdx, i);
std::cout << " i " << i << " " << phi << std::endl; std::cout << " i " << i << " " << phi << std::endl;
} }
std::cout << " fugacity for gas phase " << std::endl; std::cout << " fugacity for gas phase " << std::endl;
for (int i = 0; i < numComponents; ++i) { for (int i = 0; i < numComponents; ++i) {
auto phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, gasPhaseIdx, i); auto phi = FluidSystem::fugacityCoefficient(fluid_state, paramCache, gasPhaseIdx, i);
std::cout << " i " << i << " " << phi << std::endl; std::cout << " i " << i << " " << phi << std::endl;
} }
std::cout << " density for oil phase " << std::endl; std::cout << " density for oil phase " << std::endl;
std::cout << FluidSystem::density(fluidState, paramCache, oilPhaseIdx) << std::endl; std::cout << FluidSystem::density(fluid_state, paramCache, oilPhaseIdx) << std::endl;
std::cout << " density for gas phase " << std::endl; std::cout << " density for gas phase " << std::endl;
std::cout << FluidSystem::density(fluidState, paramCache, gasPhaseIdx) << std::endl; std::cout << FluidSystem::density(fluid_state, paramCache, gasPhaseIdx) << std::endl;
std::cout << " viscosities for oil " << std::endl; std::cout << " viscosities for oil " << std::endl;
std::cout << FluidSystem::viscosity(fluidState, paramCache, oilPhaseIdx) << std::endl; std::cout << FluidSystem::viscosity(fluid_state, paramCache, oilPhaseIdx) << std::endl;
std::cout << " viscosities for gas " << std::endl; std::cout << " viscosities for gas " << std::endl;
std::cout << FluidSystem::viscosity(fluidState, paramCache, gasPhaseIdx) << std::endl; std::cout << FluidSystem::viscosity(fluid_state, paramCache, gasPhaseIdx) << std::endl;
std::cout << "So = " << So <<std::endl; std::cout << "So = " << So <<std::endl;
std::cout << "Sg = " << Sg <<std::endl; std::cout << "Sg = " << Sg <<std::endl;
} } */
}//end solve }//end solve
/*! /*!
@ -670,32 +703,32 @@ protected:
// For the function flash_2ph, we should carry the derivatives in, and // For the function flash_2ph, we should carry the derivatives in, and
// return with the correct and updated derivatives // return with the correct and updated derivatives
template <class FluidState, class ComponentVector> template <class FluidState, class ComponentVector>
static void flash_2ph(const ComponentVector& z, static void flash_2ph(const ComponentVector& z_scalar,
const std::string& flash_2p_method, const std::string& flash_2p_method,
ComponentVector& K, ComponentVector& K_scalar,
typename FluidState::Scalar& L, typename FluidState::Scalar& L_scalar,
FluidState& fluid_state, FluidState& fluid_state_scalar,
int verbosity = 0) { int verbosity = 0) {
if (verbosity >= 1) { if (verbosity >= 1) {
std::cout << "Cell is two-phase! Solve Rachford-Rice with initial K = [" << K << "]" << std::endl; std::cout << "Cell is two-phase! Solve Rachford-Rice with initial K = [" << K_scalar << "]" << std::endl;
} }
// Rachford Rice equation to get initial L for composition solver
L = solveRachfordRice_g_(K, z, verbosity);
// Calculate composition using nonlinear solver // Calculate composition using nonlinear solver
// Newton // Newton
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;
} }
newtonCompositionUpdate_(K, L, fluid_state, z, verbosity); newtonCompositionUpdate_(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, L, fluid_state, z, /*standAlone=*/true, verbosity); 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);
newtonCompositionUpdate_(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::runtime_error("unknown two phase flash method " + flash_2p_method + " is specified");
} }
@ -717,6 +750,7 @@ protected:
// Compute x and y from K, L and Z // Compute x and y from K, L and Z
computeLiquidVapor_(fluidState, L, K, globalComposition); computeLiquidVapor_(fluidState, L, K, globalComposition);
std::cout << " the current L is " << Opm::getValue(L) << std::endl;
// Print initial condition // Print initial condition
if (verbosity >= 1) { if (verbosity >= 1) {
@ -1034,15 +1068,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& fluidState, const ComponentVector& globalComposition, static void successiveSubstitutionComposition_(ComponentVector& K, typename ComponentVector::field_type& L, FlashFluidState& fluidState, const ComponentVector& globalComposition,
bool standAlone, int verbosity) const bool newton_afterwards, const int verbosity)
{ {
// std::cout << " Evaluation in successiveSubstitutionComposition_ is " << Dune::className(L) << std::endl; // std::cout << " Evaluation in successiveSubstitutionComposition_ is " << Dune::className(L) << std::endl;
// 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.
int maxIterations; const int maxIterations = newton_afterwards ? 3 : 10;
if (standAlone == true)
maxIterations = 100;
else
maxIterations = 10;
// 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());

3
tests/test_chiflash.cpp
View File

@ -110,7 +110,8 @@ void testChiFlash()
const double flash_tolerance = 1.e-12; // just to test the setup in co2-compositional const double flash_tolerance = 1.e-12; // just to test the setup in co2-compositional
const int flash_verbosity = 1; const int flash_verbosity = 1;
// const std::string flash_twophase_method = "newton"; // "ssi" // const std::string flash_twophase_method = "newton"; // "ssi"
const std::string flash_twophase_method = "ssi"; // const std::string flash_twophase_method = "ssi";
const std::string flash_twophase_method = "ssi+newton";
// TODO: should we set these? // TODO: should we set these?
// Set initial K and L // Set initial K and L