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

WIP in getting the derivatives correct after flash calculation

Browse Source
This commit is contained in:
Kai Bao 2022-01-18 11:29:56 +01:00 committed by Trine Mykkeltvedt
parent 3368256ecc
commit a216d61377
2 changed files with 18 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -143,11 +143,6 @@ public:
fluid_state_scalar.setPressure(FluidSystem::gasPhaseIdx, fluid_state_scalar.setPressure(FluidSystem::gasPhaseIdx,
Opm::getValue(fluid_state.pressure(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))); fluid_state_scalar.setTemperature(Opm::getValue(fluid_state.temperature(0)));
// Rachford Rice equation to get initial L for composition solver // Rachford Rice equation to get initial L for composition solver
@ -169,6 +164,9 @@ public:
// Update the composition if cell is two-phase // Update the composition if cell is two-phase
if ( !isStable) { if ( !isStable) {
flash_2ph(z_scalar, twoPhaseMethod, K_scalar, L_scalar, fluid_state_scalar, verbosity); flash_2ph(z_scalar, twoPhaseMethod, K_scalar, L_scalar, fluid_state_scalar, verbosity);
// TODO: update the fluid_state, and also get the derivatives correct.
} 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_(fluid_state, z, verbosity); L = li_single_phase_label_(fluid_state, z, verbosity);
@ -719,7 +717,7 @@ protected:
if (verbosity >= 1) { if (verbosity >= 1) {
std::cout << "Calculate composition using Newton." << std::endl; std::cout << "Calculate composition using Newton." << std::endl;
} }
newtonCompositionUpdate_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, verbosity); 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) {
@ -728,20 +726,23 @@ protected:
successiveSubstitutionComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, false, verbosity); 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); successiveSubstitutionComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, true, verbosity);
newtonCompositionUpdate_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, verbosity); 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::runtime_error("unknown two phase flash method " + flash_2p_method + " is specified");
} }
} }
template <class FlashFluidState, class ComponentVector> template <class FlashFluidState, class ComponentVector>
static void newtonCompositionUpdate_(ComponentVector& K, typename FlashFluidState::Scalar& L, static void newtonComposition_(ComponentVector& K, typename FlashFluidState::Scalar& L,
FlashFluidState& fluidState, const ComponentVector& globalComposition, FlashFluidState& fluidState, const ComponentVector& globalComposition,
int verbosity) int verbosity)
{ {
// Note: due to the need for inverse flash update for derivatives, the following two can be different
// Looking for a good way to organize them
constexpr size_t num_equations = numMisciblePhases * numMiscibleComponents + 1; constexpr size_t num_equations = numMisciblePhases * numMiscibleComponents + 1;
constexpr size_t num_primary_variables = numMisciblePhases * numMiscibleComponents + 1;
using NewtonVector = Dune::FieldVector<Scalar, num_equations>; using NewtonVector = Dune::FieldVector<Scalar, num_equations>;
using NewtonMatrix = Dune::FieldMatrix<Scalar, num_equations, num_equations>; using NewtonMatrix = Dune::FieldMatrix<Scalar, num_equations, num_primary_variables>;
constexpr Scalar tolerance = 1.e-8; constexpr Scalar tolerance = 1.e-8;
NewtonVector soln(0.); NewtonVector soln(0.);
@ -777,9 +778,9 @@ protected:
} }
// AD type // AD type
using Eval = DenseAd::Evaluation<Scalar, num_equations>; using Eval = DenseAd::Evaluation<Scalar, num_primary_variables>;
// TODO: we might need to use numMiscibleComponents here // TODO: we might need to use numMiscibleComponents here
std::vector<Eval> x(num_equations), y(num_equations); std::vector<Eval> x(numComponents), y(numComponents);
Eval l; Eval l;
// TODO: I might not need to set soln anything here. // TODO: I might not need to set soln anything here.
@ -818,6 +819,7 @@ protected:
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
paramCache.updatePhase(flash_fluid_state, phaseIdx); paramCache.updatePhase(flash_fluid_state, phaseIdx);
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
// TODO: will phi here carry the correct derivatives?
Eval phi = FluidSystem::fugacityCoefficient(flash_fluid_state, paramCache, phaseIdx, compIdx); Eval phi = FluidSystem::fugacityCoefficient(flash_fluid_state, paramCache, phaseIdx, compIdx);
flash_fluid_state.setFugacityCoefficient(phaseIdx, compIdx, phi); flash_fluid_state.setFugacityCoefficient(phaseIdx, compIdx, phi);
} }
@ -825,7 +827,7 @@ protected:
bool converged = false; bool converged = false;
unsigned iter = 0; unsigned iter = 0;
constexpr unsigned max_iter = 1000; constexpr unsigned max_iter = 1000;
while (!converged && iter < max_iter) { while (iter < max_iter) {
// assembling the Jacobian and residuals // assembling the Jacobian and residuals
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
@ -833,7 +835,7 @@ protected:
// z - L*x - (1-L) * y // z - L*x - (1-L) * y
auto local_res = -Opm::getValue(globalComposition[compIdx]) + l * x[compIdx] + (1 - l) * y[compIdx]; auto local_res = -Opm::getValue(globalComposition[compIdx]) + l * x[compIdx] + (1 - l) * y[compIdx];
res[compIdx] = Opm::getValue(local_res); res[compIdx] = Opm::getValue(local_res);
for (unsigned i = 0; i < num_equations; ++i) { for (unsigned i = 0; i < num_primary_variables; ++i) {
jac[compIdx][i] = local_res.derivative(i); jac[compIdx][i] = local_res.derivative(i);
} }
} }

2
tests/test_chiflash.cpp
View File

@ -78,7 +78,7 @@ void testChiFlash()
fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp1Idx, comp[1]); fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp1Idx, comp[1]);
fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp2Idx, comp[2]); fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp2Idx, comp[2]);
// TODO: why we need this one? But without this, it caused problem for the ParamCache // It is used here only for calculate the z
fs.setSaturation(FluidSystem::oilPhaseIdx, sat[0]); fs.setSaturation(FluidSystem::oilPhaseIdx, sat[0]);
fs.setSaturation(FluidSystem::gasPhaseIdx, sat[1]); fs.setSaturation(FluidSystem::gasPhaseIdx, sat[1]);