From 83d4194876ba0f2ed1159cc06ab5ca65168b03df Mon Sep 17 00:00:00 2001 From: Trine Mykkeltvedt Date: Fri, 20 May 2022 13:13:11 +0200 Subject: [PATCH] fixed the derivatives in PengRobinsonMixture with a rewrite of a term to get pri_jac correct, also a getValue to get sec_jac equal, did rewriting and printing in ugly way jus t compare with Julia-code values --- opm/material/constraintsolvers/ChiFlash.hpp | 66 ++++++++++++- opm/material/eos/PengRobinsonMixture.hpp | 96 +++++++++++++------ .../eos/PengRobinsonParamsMixture.hpp | 30 ++++-- .../chifluid/ChiParameterCache.hpp | 21 ++++ .../juliathreecomponentfluidsystem.hh | 5 +- .../chifluid/twophasefluidsystem.hh | 5 +- 6 files changed, 177 insertions(+), 46 deletions(-) diff --git a/opm/material/constraintsolvers/ChiFlash.hpp b/opm/material/constraintsolvers/ChiFlash.hpp index f1a3dc497..9d0c97a38 100644 --- a/opm/material/constraintsolvers/ChiFlash.hpp +++ b/opm/material/constraintsolvers/ChiFlash.hpp @@ -955,6 +955,9 @@ protected: auto local_res = (fluid_state.fugacity(oilPhaseIdx, compIdx) - fluid_state.fugacity(gasPhaseIdx, compIdx)); res[compIdx + numComponents] = Opm::getValue(local_res); + //std::cout << "fugacity oil = " << local_res.derivative(i) << " gas = " << fluid_state.fugacity(gasPhaseIdx, compIdx) << " comp " << compIdx << std::endl; trine + + for (unsigned i = 0; i < num_primary; ++i) { jac[compIdx + numComponents][i] = local_res.derivative(i); } @@ -1084,11 +1087,66 @@ protected: assembleNewton_ (primary_fluid_state, primary_z, pri_jac, pri_res); - // not totally sure the following matrix operations are correct - pri_jac.invert(); - sec_jac.template leftmultiply(pri_jac); + //corresponds to julias J_s + std::cout << "sec_jac:" << std::endl; + std::cout << "[" << sec_jac[0][0] << " " << sec_jac[0][1] << " " << sec_jac[0][2] << " " << sec_jac[0][3] << "] " << std::endl; + std::cout << "[" << sec_jac[1][0] << " " << sec_jac[1][1] << " " << sec_jac[1][2] << " " << sec_jac[1][3] << "] " << std::endl; + std::cout << "[" << sec_jac[2][0] << " " << sec_jac[2][1] << " " << sec_jac[2][2] << " " << sec_jac[2][3] << "] " << std::endl; + std::cout << "[" << sec_jac[3][0] << " " << sec_jac[3][1] << " " << sec_jac[3][2] << " " << sec_jac[3][3] << "] " << std::endl; + std::cout << "[" << sec_jac[4][0] << " " << sec_jac[4][1] << " " << sec_jac[4][2] << " " << sec_jac[4][3] << "] " << std::endl; + std::cout << "[" << sec_jac[5][0] << " " << sec_jac[5][1] << " " << sec_jac[5][2] << " " << sec_jac[5][3] << "] " << std::endl; + std::cout << "[" << sec_jac[6][0] << " " << sec_jac[6][1] << " " << sec_jac[6][2] << " " << sec_jac[6][3] << "] " << std::endl; + + //corresponds to julias J_p (we miss d/dt, and have d/dL instead of d/dV) + std::cout << "pri_jac:" << std::endl; + std::cout << "[" << pri_jac[0][0] << " " << pri_jac[0][1] << " " << pri_jac[0][2] << " " << pri_jac[0][3] << " " << pri_jac[0][4]<< " " << pri_jac[0][5] << " " << pri_jac[0][6]<< "] " << std::endl; + std::cout << "[" << pri_jac[1][0] << " " << pri_jac[1][1] << " " << pri_jac[1][2] << " " << pri_jac[1][3] << " " << pri_jac[1][4]<< " " << pri_jac[1][5] << " " << pri_jac[1][6]<< "] " << std::endl; + std::cout << "[" << pri_jac[2][0] << " " << pri_jac[2][1] << " " << pri_jac[2][2] << " " << pri_jac[2][3] << " " << pri_jac[2][4]<< " " << pri_jac[2][5] << " " << pri_jac[2][6]<< "] " << std::endl; + std::cout << "[" << pri_jac[3][0] << " " << pri_jac[3][1] << " " << pri_jac[3][2] << " " << pri_jac[3][3] << " " << pri_jac[3][4]<< " " << pri_jac[3][5] << " " << pri_jac[3][6]<< "] " << std::endl; + std::cout << "[" << pri_jac[4][0] << " " << pri_jac[4][1] << " " << pri_jac[4][2] << " " << pri_jac[4][3] << " " << pri_jac[4][4]<< " " << pri_jac[4][5] << " " << pri_jac[4][6]<< "] " << std::endl; + std::cout << "[" << pri_jac[5][0] << " " << pri_jac[5][1] << " " << pri_jac[5][2] << " " << pri_jac[5][3] << " " << pri_jac[5][4]<< " " << pri_jac[5][5] << " " << pri_jac[5][6]<< "] " << std::endl; + std::cout << "[" << pri_jac[6][0] << " " << pri_jac[6][1] << " " << pri_jac[6][2] << " " << pri_jac[6][3] << " " << pri_jac[6][4]<< " " << pri_jac[6][5] << " " << pri_jac[6][6]<< "] " << std::endl; + + + SecondaryNewtonMatrix xx; + pri_jac.solve(xx,sec_jac); + std::cout << " corresponding to julia-code value and updated J_s " << std::endl; + std::cout << "x1 = [" << x[0] << " " << xx[0][0] << " " << xx[0][1] << " " << xx[0][2] << " " << xx[0][3] << "] " << std::endl; + std::cout << "x2 = [" << x[1] << " " << xx[1][0] << " " < xxx (do this properly to clean up =) + // z3 = 1 -z2 -z1; + // dx1/dp dx1/dt (dx1/dz1-dx1/dz3) (dx1/dz2 - dx1/dz3) + using TertiaryMatrix = Dune::FieldMatrix; + TertiaryMatrix xxx; + xxx[0][0] = xx[0][0]; + xxx[1][0] = xx[1][0]; + xxx[2][0] = xx[2][0]; + xxx[3][0] = xx[0][0]; + xxx[4][0] = xx[1][0]; + xxx[5][0] = xx[2][0]; + xxx[6][0] = xx[3][0]; + for (unsigned i = 0; i < primary_num_pv; ++i) { // 7 rekker + xxx[i][1] = Opm::getValue(xx[i][1])-Opm::getValue(xx[i][3]); + xxx[i][2] = Opm::getValue(xx[i][2])-Opm::getValue(xx[i][3]); + } + std::cout << " corresponding to julia-code value and derivatives listed in test_setup NB CHANGE SIGN, and we dont have d/dT " << std::endl; + std::cout << "x1 = [" << x[0] << " " << xxx[0][0] << " " << xxx[0][1] << " " << xxx[0][2] << "] " << std::endl; + std::cout << "x2 = [" << x[1] << " " << xxx[1][0] << " " < static void evalJacobian(const ComponentVector& globalComposition, diff --git a/opm/material/eos/PengRobinsonMixture.hpp b/opm/material/eos/PengRobinsonMixture.hpp index d93453999..850df6c56 100644 --- a/opm/material/eos/PengRobinsonMixture.hpp +++ b/opm/material/eos/PengRobinsonMixture.hpp @@ -94,21 +94,21 @@ public: unsigned phaseIdx, unsigned compIdx) { -#warning also a HACK, should investigate why - auto fs = fs2; - double sumx = 0.0; - for (int i = 0; i < FluidState::numComponents; ++i) - sumx += Opm::scalarValue(fs.moleFraction(phaseIdx, i)); - if (sumx < 0.95) { - double alpha = 0.95/sumx; - std::cerr << "normalize: " << sumx - << " alpha: " << alpha << "\n"; - for (int i = 0; i < FluidState::numComponents; ++i) - fs.setMoleFraction(phaseIdx, i, alpha*fs.moleFraction(phaseIdx, i)); - } +// #warning also a HACK, should investigate why + auto fs = fs2; + // double sumx = 0.0; + // for (int i = 0; i < FluidState::numComponents; ++i) + // sumx += Opm::scalarValue(fs.moleFraction(phaseIdx, i)); + // if (sumx < 0.95) { + // double alpha = 0.95/sumx; + // std::cerr << "normalize: " << sumx + // << " alpha: " << alpha << "\n"; + // for (int i = 0; i < FluidState::numComponents; ++i) + // fs.setMoleFraction(phaseIdx, i, alpha*fs.moleFraction(phaseIdx, i)); + // } - auto params = params2; - params.updatePhase(fs, phaseIdx); + auto params = params2; + // params.updatePhase(fs, phaseIdx); // note that we normalize the component mole fractions, so // that their sum is 100%. This increases numerical stability // considerably if the fluid state is not physical. @@ -126,29 +126,69 @@ public: LhsEval Astar = params.a(phaseIdx)*p/(RT*RT); LhsEval Bstar = params.b(phaseIdx)*p/(RT); + //MAYBE REMOVE THIS COMPLETELY ?? // calculate delta_i (see: Reid, p. 145) LhsEval sumMoleFractions = 0.0; for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) sumMoleFractions += fs.moleFraction(phaseIdx, compJIdx); LhsEval deltai = 2*sqrt(params.aPure(phaseIdx, compIdx))/params.a(phaseIdx); - LhsEval tmp = 0; - for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) { - tmp += - fs.moleFraction(phaseIdx, compJIdx) - / sumMoleFractions - * sqrt(params.aPure(phaseIdx, compJIdx)) - * (1.0 - StaticParameters::interactionCoefficient(compIdx, compJIdx)); - }; - deltai *= tmp; + //std::cout << "params.aPure [" << params.aPure(phaseIdx, compIdx) << "] " << std::endl; + LhsEval tmp = 0; + for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) { + tmp += + fs.moleFraction(phaseIdx, compJIdx) + / sumMoleFractions + * sqrt(params.aPure(phaseIdx, compJIdx)); + //* (1.0 - StaticParameters::interactionCoefficient(compIdx, compJIdx)); + }; + deltai *= tmp; + + // Calculate A_s LIKE IN JULIA CODE + + LhsEval A_s = 0.0; + for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) + { + //param + A_s += params.aCache(phaseIdx, compIdx, compJIdx) * fs.moleFraction(phaseIdx, compJIdx) * p / (RT * RT); + } + //LhsEval A_s = 0.5*Astar*deltai; + + + LhsEval julia_alpha; + LhsEval julia_betta; + LhsEval julia_gamma; + LhsEval ln_phi; + LhsEval fugCoeff; + + Scalar m1; + Scalar m2; + + m1 = 0.5*(u + std::sqrt(u*u - 4*w)); + m2 = 0.5*(u - std::sqrt(u*u - 4*w)); + LhsEval base = - (2*Z + Bstar*(u + std::sqrt(u*u - 4*w))) / - (2*Z + Bstar*(u - std::sqrt(u*u - 4*w))); + (Z + Bstar*m1) / + (Z + Bstar*m2); LhsEval expo = Astar/(Bstar*std::sqrt(u*u - 4*w))*(bi_b - deltai); - LhsEval fugCoeff = - exp(bi_b*(Z - 1))/max(1e-9, Z - Bstar) * - pow(base, expo); + //julia: + //α = -log(Z - B) + (B_i[i]/B)*(Z - 1) + //β = log((Z + m2*B)/(Z + m1*B))*A/(Δm*B) + //γ = (2/A)*A_s - B_i[i]/B + //return α + β*γ + + julia_alpha = -log(Z-Bstar)+bi_b*(Z - 1); + julia_betta = log((Z+m2*Bstar)/(Z+m1*Bstar))*Astar/((m1-m2)*Bstar); + //julia_gamma = -bi_b+deltai;//expo;//(2/Astar)*A_s -bi_b/Bstar; + julia_gamma = (2 / Astar )* A_s - bi_b; + + ln_phi = julia_alpha + (julia_betta*julia_gamma); + + + + fugCoeff = exp(ln_phi); + //exp(bi_b*(Z - 1))/max(1e-9, Z - Bstar) * pow(base, expo); //////// // limit the fugacity coefficient to a reasonable range: diff --git a/opm/material/eos/PengRobinsonParamsMixture.hpp b/opm/material/eos/PengRobinsonParamsMixture.hpp index b58876e28..1f02a12a7 100644 --- a/opm/material/eos/PengRobinsonParamsMixture.hpp +++ b/opm/material/eos/PengRobinsonParamsMixture.hpp @@ -66,6 +66,15 @@ class PengRobinsonParamsMixture typedef MathToolbox Toolbox; public: + + /*! + * \brief TODO + */ + Scalar getaCache(unsigned compIIdx, unsigned compJIdx ) + { + return aCache_[compIIdx][compJIdx]; + } + /*! * \brief Update Peng-Robinson parameters for the pure components. */ @@ -99,12 +108,12 @@ public: Scalar f_omega; - if (useSpe5Relations) { + // if (useSpe5Relations) { if (omega < 0.49) f_omega = 0.37464 + omega*(1.54226 + omega*(-0.26992)); else f_omega = 0.379642 + omega*(1.48503 + omega*(-0.164423 + omega*0.016666)); - } - else - f_omega = 0.37464 + omega*(1.54226 - omega*0.26992); + // } + // else + // f_omega = 0.37464 + omega*(1.54226 - omega*0.26992); Valgrind::CheckDefined(f_omega); @@ -135,6 +144,7 @@ public: template void updateMix(const FluidState& fs) { + using FlashEval = typename FluidState::Scalar; Scalar sumx = 0.0; for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) sumx += fs.moleFraction(phaseIdx, compIdx); @@ -144,16 +154,16 @@ public: // // See: R. Reid, et al.: The Properties of Gases and Liquids, // 4th edition, McGraw-Hill, 1987, p. 82 - Scalar newA = 0; - Scalar newB = 0; + FlashEval newA = 0; + FlashEval newB = 0; for (unsigned compIIdx = 0; compIIdx < numComponents; ++compIIdx) { - const Scalar moleFracI = fs.moleFraction(phaseIdx, compIIdx); - Scalar xi = max(0.0, min(1.0, moleFracI)); + const FlashEval moleFracI = fs.moleFraction(phaseIdx, compIIdx); + FlashEval xi = max(0.0, min(1.0, moleFracI)); Valgrind::CheckDefined(xi); for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) { - const Scalar moleFracJ = fs.moleFraction(phaseIdx, compJIdx ); - Scalar xj = max(0.0, min(1.0, moleFracJ)); + const FlashEval moleFracJ = fs.moleFraction(phaseIdx, compJIdx ); + FlashEval xj = max(0.0, min(1.0, moleFracJ)); Valgrind::CheckDefined(xj); // mixing rule from Reid, page 82 diff --git a/opm/material/fluidsystems/chifluid/ChiParameterCache.hpp b/opm/material/fluidsystems/chifluid/ChiParameterCache.hpp index 6a299597d..b16b712f5 100644 --- a/opm/material/fluidsystems/chifluid/ChiParameterCache.hpp +++ b/opm/material/fluidsystems/chifluid/ChiParameterCache.hpp @@ -177,6 +177,27 @@ public: }; } + + + /*! + * \brief TODO + * + * \param phaseIdx The fluid phase of interest + * \param compIdx The component phase of interest + * \param compJIdx Additional component index + */ + Scalar aCache(unsigned phaseIdx, unsigned compIdx, unsigned compJIdx) + { + switch (phaseIdx) + { + case oilPhaseIdx: return oilPhaseParams_.getaCache(compIdx,compJIdx); + case gasPhaseIdx: return gasPhaseParams_.getaCache(compIdx,compJIdx); + default: + throw std::logic_error("The aCache() parameter is only defined for " + "oil phase"); + }; + } + /*! * \brief Returns the molar volume of a phase [m^3/mol] * diff --git a/opm/material/fluidsystems/chifluid/juliathreecomponentfluidsystem.hh b/opm/material/fluidsystems/chifluid/juliathreecomponentfluidsystem.hh index fd9a04474..ce5ce8299 100644 --- a/opm/material/fluidsystems/chifluid/juliathreecomponentfluidsystem.hh +++ b/opm/material/fluidsystems/chifluid/juliathreecomponentfluidsystem.hh @@ -160,8 +160,9 @@ namespace Opm { // TODO: here the derivatives for the phi are dropped. Should we keep the derivatives against the pressure // and temperature? - Scalar phi = Opm::getValue( - PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx)); + LhsEval phi = PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx); + //Scalar phi = Opm::getValue( + // PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx)); return phi; } diff --git a/opm/material/fluidsystems/chifluid/twophasefluidsystem.hh b/opm/material/fluidsystems/chifluid/twophasefluidsystem.hh index aafb176e1..670b5b674 100644 --- a/opm/material/fluidsystems/chifluid/twophasefluidsystem.hh +++ b/opm/material/fluidsystems/chifluid/twophasefluidsystem.hh @@ -272,8 +272,9 @@ public: assert(0 <= phaseIdx && phaseIdx < numPhases); assert(0 <= compIdx && compIdx < numComponents); - Scalar phi = Opm::getValue( - PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx)); + LhsEval phi = PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx); + //Scalar phi = Opm::getValue( + // PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx)); return phi;