making the types for LBCJulia more proper

We will need to handle AD type with this function.
This commit is contained in:
Kai Bao 2021-12-13 09:59:21 +01:00 committed by Trine Mykkeltvedt
parent bb86f443c7
commit 30c0ce7fac

View File

@ -215,20 +215,20 @@ public:
constexpr Scalar rankine = 5. / 9.;
constexpr Scalar psia = 6.894757293168360e+03;
constexpr Scalar R = 8.3144598;
const Scalar T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const Scalar P = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
const Scalar Z = Opm::decay<LhsEval>(fluidState.compressFactor(phaseIdx));
const Scalar rho = P / (R * T * Z);
const LhsEval T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval P = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
const LhsEval Z = Opm::decay<LhsEval>(fluidState.compressFactor(phaseIdx));
const LhsEval rho = P / (R * T * Z);
Scalar P_pc = 0.;
Scalar T_pc = 0.;
Scalar Vc = 0.;
Scalar mwc = 0.;
Scalar a = 0.;
Scalar b = 0.;
LhsEval P_pc = 0.;
LhsEval T_pc = 0.;
LhsEval Vc = 0.;
LhsEval mwc = 0.;
LhsEval a = 0.;
LhsEval b = 0.;
for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
const Scalar mol_frac = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
const Scalar mol_weight = FluidSystem::molarMass(compIdx); // TODO: check values
const LhsEval mol_frac = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
const LhsEval mol_weight = FluidSystem::molarMass(compIdx); // TODO: check values
const Scalar p_c = FluidSystem::criticalPressure(compIdx);
const Scalar T_c = FluidSystem::criticalTemperature(compIdx);
const Scalar v_c = FluidSystem::criticalVolume(compIdx);
@ -237,14 +237,14 @@ public:
T_pc += mol_frac * T_c;
Vc += mol_frac * v_c;
const Scalar tr = T / T_c;
const LhsEval tr = T / T_c;
const Scalar Tc = T_c / rankine;
const Scalar Pc = p_c / psia;
const Scalar mwi = Opm::sqrt(mol_factor * mol_weight);
const Scalar e_i = 5.4402 * Opm::pow(Tc, 1./6.) / (mwi * Opm::pow(Pc, 2./3.) * 1.e-3);
const LhsEval mwi = Opm::sqrt(mol_factor * mol_weight);
const LhsEval e_i = 5.4402 * Opm::pow(Tc, 1./6.) / (mwi * Opm::pow(Pc, 2./3.) * 1.e-3);
Scalar mu_i;
LhsEval mu_i;
if (tr > 1.5) {
mu_i = 17.78e-5 * Opm::pow(4.58*tr - 1.67, 0.625) / e_i;
} else {
@ -253,12 +253,12 @@ public:
a += mol_frac * mu_i * mwi;
b += mol_frac * mwi;
}
const Scalar mu_atm = a / b;
const Scalar e_mix = 5.4402 * Opm::pow(T_pc/rankine, 1./6.) /
const LhsEval mu_atm = a / b;
const LhsEval e_mix = 5.4402 * Opm::pow(T_pc/rankine, 1./6.) /
(Opm::sqrt(mol_factor * mwc) * Opm::pow(P_pc/psia, 2./3.) * (1e-3));
const Scalar rhor = Vc * rho;
const LhsEval rhor = Vc * rho;
Scalar corr = 0.;
LhsEval corr = 0.;
const std::vector<Scalar> LBC{0.10230, 0.023364, 0.058533, -0.040758, 0.0093324};
const Scalar shift = -1.e-4;
for (unsigned i = 0; i < 5; ++i) {