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

@@ -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_<PrimaryFlashFluidState, PrimaryComponentVector, primary_num_pv, num_equations>
             (primary_fluid_state, primary_z, pri_jac, pri_res);
 
-        // not totally sure the following matrix operations are correct
+        //corresponds to julias J_s
-        pri_jac.invert();
+        std::cout << "sec_jac:" << std::endl;
-        sec_jac.template leftmultiply(pri_jac);
+        std::cout << "[" << sec_jac[0][0] << "  " << sec_jac[0][1] << "  " << sec_jac[0][2] << "  " << sec_jac[0][3] << "]  " << std::endl;
-        // TODO: then beginning from that point
+        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] << "  " <<xx[1][1] << "  " << xx[1][2] << "  " << xx[1][3] << "]  " << std::endl;
+        std::cout << "x3 = [" << x[2] << "  " << xx[2][0] << "  " << xx[2][1] << "  " << xx[2][2] << "  " << xx[2][3] << "]  " << std::endl;
+        std::cout << "y1 = [" << y[0] << "  " << xx[3][0] << "  " << xx[3][1] << "  " << xx[3][2] << "  " << xx[3][3] << "]  " << std::endl;
+        std::cout << "y2 = [" << y[1] << "  " << xx[4][0] << "  " << xx[4][1] << "  " << xx[4][2] << "  " << xx[4][3] << "]  " << std::endl;
+        std::cout << "y3 = [" << y[2] << "  " << xx[5][0] << "  " << xx[5][1] << "  " << xx[5][2] << "  " << xx[5][3] << "]  " << std::endl;
+        std::cout << "L = [" << L << "  " << xx[6][0] << "  " << xx[6][1] << "  " << xx[6][2] << "  " << xx[6][3] << "]  " << std::endl; 
+
+        // rewrite like Olav xx --> 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<Scalar, num_equations, secondary_num_pv-1>;
+        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] << "  " <<xxx[1][1] << "  " << xxx[1][2] <<  "]  " << std::endl;
+        std::cout << "x3 = [" << x[2] << "  " << xxx[2][0] << "  " << xxx[2][1] << "  " << xxx[2][2] <<  "]  " << std::endl;
+        std::cout << "y1 = [" << y[0] << "  " << xxx[3][0] << "  " << xxx[3][1] << "  " << xxx[3][2] <<  "]  " << std::endl;
+        std::cout << "y2 = [" << y[1] << "  " << xxx[4][0] << "  " << xxx[4][1] << "  " << xxx[4][2] <<  "]  " << std::endl;
+        std::cout << "y3 = [" << y[2] << "  " << xxx[5][0] << "  " << xxx[5][1] << "  " << xxx[5][2] <<  "]  " << std::endl;
+        std::cout << "L = [" << L << "  " << xxx[6][0] << "  " << xxx[6][1] << "  " << xxx[6][2] << "  " << xx[6][3] << "]  " << std::endl;
+
+        // TODO: then beginning from that point
+
+    }//end updateDerivative
 
     /* template <class Vector, class Matrix, class Eval, class ComponentVector>
     static void evalJacobian(const ComponentVector& globalComposition,

opm/material/eos/PengRobinsonMixture.hpp

@@ -94,21 +94,21 @@ public:
                                               unsigned phaseIdx,
                                               unsigned compIdx)
     {
-#warning also a HACK, should investigate why
+// #warning also a HACK, should investigate why
          auto fs = fs2;
-        double sumx = 0.0;
+        //  double sumx = 0.0;
-        for (int i = 0; i < FluidState::numComponents; ++i)
+        //  for (int i = 0; i < FluidState::numComponents; ++i)
-            sumx += Opm::scalarValue(fs.moleFraction(phaseIdx, i));
+        //      sumx += Opm::scalarValue(fs.moleFraction(phaseIdx, i));
-        if (sumx < 0.95)  {
+        //  if (sumx < 0.95)  {
-            double alpha = 0.95/sumx;
+        //      double alpha = 0.95/sumx;
-            std::cerr << "normalize: " << sumx
+        //      std::cerr << "normalize: " << sumx
-                      << " alpha: " << alpha << "\n";
+        //                << " alpha: " << alpha << "\n";
-            for (int i = 0; i < FluidState::numComponents; ++i)
+        //      for (int i = 0; i < FluidState::numComponents; ++i)
-                fs.setMoleFraction(phaseIdx, i, alpha*fs.moleFraction(phaseIdx, i));
+        //          fs.setMoleFraction(phaseIdx, i, alpha*fs.moleFraction(phaseIdx, i));
-        }
+        //  }
 
          auto params = params2;
-        params.updatePhase(fs, phaseIdx);
+        // 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);
+        //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))
+                 * sqrt(params.aPure(phaseIdx, compJIdx));
-                * (1.0 - StaticParameters::interactionCoefficient(compIdx, 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))) /
+            (Z + Bstar*m1) /
-            (2*Z + Bstar*(u - std::sqrt(u*u - 4*w)));
+            (Z + Bstar*m2);
         LhsEval expo =  Astar/(Bstar*std::sqrt(u*u - 4*w))*(bi_b - deltai);
 
-        LhsEval fugCoeff =
+        //julia:
-            exp(bi_b*(Z - 1))/max(1e-9, Z - Bstar) *
+        //α = -log(Z - B) + (B_i[i]/B)*(Z - 1)
-            pow(base, expo);
+        //β = 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:

opm/material/eos/PengRobinsonParamsMixture.hpp

@@ -66,6 +66,15 @@ class PengRobinsonParamsMixture
     typedef MathToolbox<Scalar> 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
+           // else
-                f_omega = 0.37464 + omega*(1.54226 - omega*0.26992);
+           //     f_omega = 0.37464 + omega*(1.54226 - omega*0.26992);
 
             Valgrind::CheckDefined(f_omega);
 
@@ -135,6 +144,7 @@ public:
     template <class FluidState>
     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;
+        FlashEval newA = 0;
-        Scalar newB = 0;
+        FlashEval newB = 0;
         for (unsigned compIIdx = 0; compIIdx < numComponents; ++compIIdx) {
-            const Scalar moleFracI = fs.moleFraction(phaseIdx, compIIdx);
+            const FlashEval moleFracI = fs.moleFraction(phaseIdx, compIIdx);
-            Scalar xi = max(0.0, min(1.0, moleFracI));
+            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 );
+                const FlashEval moleFracJ = fs.moleFraction(phaseIdx, compJIdx );
-                Scalar xj = max(0.0, min(1.0, moleFracJ));
+                FlashEval xj = max(0.0, min(1.0, moleFracJ));
                 Valgrind::CheckDefined(xj);
 
                 // mixing rule from Reid, page 82

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]
      *

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(
+            LhsEval phi = PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx); 
-                    PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
+            //Scalar phi = Opm::getValue(
+            //        PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
             return phi;
         }
 

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(
+                LhsEval phi = PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
-                PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
+            //Scalar phi = Opm::getValue(
+            //    PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
             return phi;