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

added Kais startup cleaning, no functionality change

Browse Source
This commit is contained in:
Trine Mykkeltvedt 2022-06-20 09:54:15 +02:00
parent 5b2845bf9a
commit 82c6bfe610
2 changed files with 64 additions and 77 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -58,7 +58,6 @@ namespace Opm {
template <class Scalar, class FluidSystem>
class ChiFlash
{
//using Problem = GetPropType<TypeTag, Properties::Problem>;
enum { numPhases = FluidSystem::numPhases };
enum { numComponents = FluidSystem::numComponents };
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx};
@ -69,9 +68,7 @@ class ChiFlash
numEq =
numMisciblePhases+
numMisciblePhases*numMiscibleComponents
};//pressure, saturation, composition
}; //pressure, saturation, composition
public:
/*!
@ -84,7 +81,7 @@ public:
int spatialIdx,
int verbosity,
std::string twoPhaseMethod,
Scalar tolerance)
Scalar tolerance = -1.)
{
using InputEval = typename FluidState::Scalar;
@ -93,11 +90,11 @@ public:
#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7)
Dune::FMatrixPrecision<InputEval>::set_singular_limit(1e-35);
#endif
if (tolerance <= 0) {
tolerance = std::min<Scalar>(1e-3, 1e8 * std::numeric_limits<Scalar>::epsilon());
}
if (tolerance <= 0)
tolerance = std::min<Scalar>(1e-3, 1e8*std::numeric_limits<Scalar>::epsilon());
//K and L from previous timestep (wilson and -1 initially)
// K and L from previous timestep (wilson and -1 initially)
ComponentVector K;
for(int compIdx = 0; compIdx < numComponents; ++compIdx) {
K[compIdx] = fluid_state.K(compIdx);
@ -146,7 +143,6 @@ public:
std::cout << "Perform stability test (L <= 0 or L == 1)!" << std::endl;
}
phaseStabilityTest_(isStable, K_scalar, fluid_state_scalar, z_scalar, verbosity);
}
if (verbosity >= 1) {
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;
@ -164,11 +160,9 @@ public:
} else {
// Cell is one-phase. Use Li's phase labeling method to see if it's liquid or vapor
L_scalar = li_single_phase_label_(fluid_state_scalar, z_scalar, verbosity);
single = true;
single = true;
}
// Print footer
if (verbosity >= 1) {
std::cout << "********" << std::endl;
@ -193,9 +187,7 @@ public:
fluid_state.setKvalue(compIdx, K_scalar[compIdx]);
fluid_state_scalar.setKvalue(compIdx, K_scalar[compIdx]);
}
updateDerivatives_(fluid_state_scalar, z, fluid_state, single);
}//end solve
/*!
@ -433,7 +425,6 @@ protected:
throw std::runtime_error(" Rachford-Rice with bisection failed!");
}
template <class FlashFluidState, class ComponentVector>
static void phaseStabilityTest_(bool& isStable, ComponentVector& K, FlashFluidState& fluidState, const ComponentVector& globalComposition, int verbosity)
{
@ -475,6 +466,7 @@ protected:
}
}
}
template <class FlashFluidState, class ComponentVector>
static void checkStability_(const FlashFluidState& fluidState, bool& isTrivial, ComponentVector& K, ComponentVector& xy_loc,
typename FlashFluidState::Scalar& S_loc, const ComponentVector& globalComposition, bool isGas, int verbosity)
@ -537,6 +529,7 @@ protected:
fluidState_fake.setFugacityCoefficient(phaseIdx, compIdx, phiFake);
fluidState_global.setFugacityCoefficient(phaseIdx2, compIdx, phiGlobal);
}
ComponentVector R;
for (int compIdx=0; compIdx<numComponents; ++compIdx){
@ -771,7 +764,13 @@ protected:
}
}
}
for (unsigned i = 0; i < num_equations; ++i) {
for (unsigned j = 0; j < num_primary_variables; ++j) {
std::cout << " " << jac[i][j] ;
}
std::cout << std::endl;
}
std::cout << std::endl;
if (!converged) {
throw std::runtime_error(" Newton composition update did not converge within maxIterations");
}
@ -788,7 +787,8 @@ protected:
K[idx] = K_i;
}
L = Opm::getValue(l);
fluidState.setLvalue(L); }
fluidState.setLvalue(L);
}
template <class DefectVector>
static void updateCurrentSol_(DefectVector& x, DefectVector& d)
@ -1021,7 +1021,6 @@ protected:
for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
primary_z[comp_idx] = Opm::getValue(z[comp_idx]);
}
for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
const auto x_ii = PrimaryEval(fluid_state_scalar.moleFraction(oilPhaseIdx, comp_idx), comp_idx);
primary_fluid_state.setMoleFraction(oilPhaseIdx, comp_idx, x_ii);
@ -1106,24 +1105,20 @@ protected:
deri[idx] += pz * zi.derivative(idx);
}
}
for (unsigned idx = 0; idx < num_deri; ++idx) {
x[compIdx].setDerivative(idx, deri[idx]);
}
}
// handling y
for (unsigned idx = 0; idx < num_deri; ++idx) {
deri[idx] = - xx[compIdx + numComponents][0]* p_v.derivative(idx);
}
}
for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
const double pz = -xx[compIdx + numComponents][cIdx + 1];
const auto& zi = z[cIdx];
for (unsigned idx = 0; idx < num_deri; ++idx) {
deri[idx] += pz * zi.derivative(idx);
}
}
}
for (unsigned idx = 0; idx < num_deri; ++idx) {
y[compIdx].setDerivative(idx, deri[idx]);
}
@ -1133,7 +1128,6 @@ protected:
for (unsigned idx = 0; idx < num_deri; ++idx) {
deriL[idx] = - xx[2*numComponents][0] * p_v.derivative(idx);
}
for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
const double pz = -xx[2*numComponents][cIdx + 1];
const auto& zi = z[cIdx];
@ -1169,6 +1163,7 @@ protected:
fluidState.setMoleFraction(gasPhaseIdx, compIdx, x[compIdx + numMiscibleComponents]);
}
// Compute fugacities
using ValueType = typename FluidState::Scalar;
using ParamCache = typename FluidSystem::template ParameterCache<typename FluidState::Scalar>;
@ -1198,7 +1193,7 @@ protected:
// sum(x) - sum(y) = 0
b[numMiscibleComponents*numMisciblePhases] += -x[compIdx] + x[compIdx + numMiscibleComponents];
}
}//end evalDefect
}//end valDefect
template <class FluidState, class DefectVector, class DefectMatrix, class ComponentVector>
static void evalJacobian_(DefectMatrix& A,
@ -1349,4 +1344,4 @@ protected:
} // namespace Opm
#endif
#endif

88
tests/test_chiflash.cpp
View File

@ -30,21 +30,17 @@
#include <opm/material/constraintsolvers/ChiFlash.hpp>
#include <opm/material/fluidsystems/chifluid/juliathreecomponentfluidsystem.hh>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
#include <opm/material/fluidstates/CompositionalFluidState.hpp>
#include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
#include <dune/common/parallel/mpihelper.hh>
// the following include should be removed later
// #include <opm/material/fluidsystems/chifluid/chiwoms.h>
void testChiFlash()
{
using Scalar = double;
// TODO: the name Julia should not be there, remaining to be changed
using FluidSystem = Opm::JuliaThreeComponentFluidSystem<Scalar>;
constexpr auto numComponents = FluidSystem::numComponents;
@ -52,89 +48,85 @@ void testChiFlash()
typedef Dune::FieldVector<Evaluation, numComponents> ComponentVector;
typedef Opm::CompositionalFluidState<Evaluation, FluidSystem> FluidState;
// input
// It is a three component system
// Initial: the primary variables are, pressure, molar fractions of the first and second component
Evaluation p_init = Evaluation::createVariable(10e5, 0); // 10 bar
ComponentVector comp;
comp[0] = Evaluation::createVariable(0.5, 1);
comp[1] = Evaluation::createVariable(0.3, 2);
comp[2] = 1. - comp[0] - comp[1];
ComponentVector sat;
sat[0] = 1.0; sat[1] = 1.0-sat[0];
// TODO: should we put the derivative against the temperature here?
Scalar temp = 300.0;
// From co2-compositional branch, it uses
// typedef typename FluidSystem::template ParameterCache<Scalar> ParameterCache;
FluidState fs;
// TODO: no capillary pressure for now
fs.setPressure(FluidSystem::oilPhaseIdx, p_init);
fs.setPressure(FluidSystem::gasPhaseIdx, p_init);
// TODO: not sure whether the saturation matter here.
ComponentVector sat;
// We assume that currently everything is in the oil phase
sat[0] = 1.0; sat[1] = 1.0-sat[0];
Scalar temp = 300.0;
fs.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp0Idx, comp[0]);
fs.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp1Idx, comp[1]);
fs.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp2Idx, comp[2]);
// FluidState will be the input for the flash calculation
FluidState fluid_state;
fluid_state.setPressure(FluidSystem::oilPhaseIdx, p_init);
fluid_state.setPressure(FluidSystem::gasPhaseIdx, p_init);
fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp0Idx, comp[0]);
fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp1Idx, comp[1]);
fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp2Idx, comp[2]);
fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp0Idx, comp[0]);
fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp1Idx, comp[1]);
fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp2Idx, comp[2]);
fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp0Idx, comp[0]);
fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp1Idx, comp[1]);
fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp2Idx, comp[2]);
// It is used here only for calculate the z
fs.setSaturation(FluidSystem::oilPhaseIdx, sat[0]);
fs.setSaturation(FluidSystem::gasPhaseIdx, sat[1]);
fluid_state.setSaturation(FluidSystem::oilPhaseIdx, sat[0]);
fluid_state.setSaturation(FluidSystem::gasPhaseIdx, sat[1]);
fs.setTemperature(temp);
fluid_state.setTemperature(temp);
// ParameterCache paramCache;
{
typename FluidSystem::template ParameterCache<Evaluation> paramCache;
paramCache.updatePhase(fs, FluidSystem::oilPhaseIdx);
paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx);
fs.setDensity(FluidSystem::oilPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::oilPhaseIdx));
fs.setDensity(FluidSystem::gasPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx));
paramCache.updatePhase(fluid_state, FluidSystem::oilPhaseIdx);
paramCache.updatePhase(fluid_state, FluidSystem::gasPhaseIdx);
fluid_state.setDensity(FluidSystem::oilPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::oilPhaseIdx));
fluid_state.setDensity(FluidSystem::gasPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx));
}
ComponentVector zInit(0.); // TODO; zInit needs to be normalized.
ComponentVector z(0.); // TODO; z needs to be normalized.
{
Scalar sumMoles = 0.0;
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
Scalar tmp = Opm::getValue(fs.molarity(phaseIdx, compIdx) * fs.saturation(phaseIdx));
zInit[compIdx] += Opm::max(tmp, 1e-8);
Scalar tmp = Opm::getValue(fluid_state.molarity(phaseIdx, compIdx) * fluid_state.saturation(phaseIdx));
z[compIdx] += Opm::max(tmp, 1e-8);
sumMoles += tmp;
}
}
zInit /= sumMoles;
// initialize the derivatives
// TODO: the derivative eventually should be from the reservoir flow equations
z /= sumMoles;
// p And z is the primary variables
Evaluation z_last = 1.;
for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
zInit[compIdx] = Evaluation::createVariable(Opm::getValue(zInit[compIdx]), compIdx + 1);
z_last -= zInit[compIdx];
z[compIdx] = Evaluation::createVariable(Opm::getValue(z[compIdx]), compIdx + 1);
z_last -= z[compIdx];
}
zInit[numComponents - 1] = z_last;
z[numComponents - 1] = z_last;
}
// TODO: only, p, z need the derivatives.
const double flash_tolerance = 1.e-12; // just to test the setup in co2-compositional
const int flash_verbosity = 1;
//const std::string flash_twophase_method = "ssi"; // "ssi"
const std::string flash_twophase_method = "newton";
// const std::string flash_twophase_method = "ssi+newton";
// TODO: should we set these?
// Set initial K and L
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
const Evaluation Ktmp = fs.wilsonK_(compIdx);
fs.setKvalue(compIdx, Ktmp);
const Evaluation Ktmp = fluid_state.wilsonK_(compIdx);
fluid_state.setKvalue(compIdx, Ktmp);
}
const Evaluation Ltmp = 1.;
fs.setLvalue(Ltmp);
fluid_state.setLvalue(Ltmp);
const int spatialIdx = 0;
using Flash = Opm::ChiFlash<double, FluidSystem>;
// TODO: here the zInit does not have the proper derivatives
Flash::solve(fs, zInit, spatialIdx, flash_verbosity, flash_twophase_method, flash_tolerance);
Flash::solve(fluid_state, z, spatialIdx, flash_verbosity, flash_twophase_method, flash_tolerance);
}
@ -145,4 +137,4 @@ int main(int argc, char **argv)
testChiFlash();
return 0;
}
}