some shallow cleaning up

opm/material/constraintsolvers/PTFlash.hpp

@@ -65,13 +65,13 @@ class PTFlash
     enum { numComponents = FluidSystem::numComponents };
     enum { oilPhaseIdx = FluidSystem::oilPhaseIdx};
     enum { gasPhaseIdx = FluidSystem::gasPhaseIdx};
-    enum { numMiscibleComponents = FluidSystem::numMiscibleComponents}; //octane, co2 // should be brine instead of brine here.
+    enum { numMiscibleComponents = FluidSystem::numMiscibleComponents};
     enum { numMisciblePhases = FluidSystem::numMisciblePhases}; //oil, gas
     enum {
         numEq =
         numMisciblePhases+
         numMisciblePhases*numMiscibleComponents
-    };  //pressure, saturation, composition
+    };
 
 public:
     /*!
@@ -82,9 +82,9 @@ public:
     static void solve(FluidState& fluid_state,
                       const Dune::FieldVector<typename FluidState::Scalar, numComponents>& z,
                       int spatialIdx,
-                      int verbosity,
                       std::string twoPhaseMethod,
-                      Scalar tolerance = -1.)
+                      Scalar tolerance = -1.,
+                      int verbosity = 0)
     {
 
         using InputEval = typename FluidState::Scalar;
@@ -163,20 +163,16 @@ public:
             // 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);
         }
+        fluid_state_scalar.setLvalue(L_scalar);
 
         // Print footer
         if (verbosity >= 1) {
             std::cout << "********" << std::endl;
         }
 
-        // all the solution should be processed in scalar form
+        // the flash solution process were performed in scalar form, after the flash calculation finishes,
-        // now we should update the derivatives
+        // we transform the derivatives to the final solution
-        // TODO: should be able the handle the derivatives directly, which will affect the final organization
+        // ensure that things in fluid_state_scalar is transformed to fluid_state
-        // of the code
-        // TODO: Does fluid_state_scalar contain z with derivatives?
-        fluid_state_scalar.setLvalue(L_scalar);
-        fluid_state.setLvalue(L_scalar);
-        // ensure that things in fluid_state_scalar is transformed to fluid_state 
         for (int compIdx=0; compIdx<numComponents; ++compIdx){
                 const auto x_i = fluid_state_scalar.moleFraction(oilPhaseIdx, compIdx);
                 fluid_state.setMoleFraction(oilPhaseIdx, compIdx, x_i);
@@ -188,6 +184,7 @@ public:
             fluid_state.setKvalue(compIdx, K_scalar[compIdx]);
             fluid_state_scalar.setKvalue(compIdx, K_scalar[compIdx]);
         }
+        fluid_state.setLvalue(L_scalar);
         updateDerivatives_(fluid_state_scalar, z, fluid_state, is_single_phase);
     }//end solve
 
@@ -279,7 +276,6 @@ protected:
         return L;
     }
 
-    // TODO: not totally sure whether ValueType can be obtained from Vector::field_type
     template <class Vector>
     static typename Vector::field_type rachfordRice_g_(const Vector& K, typename Vector::field_type L, const Vector& z)
     {
@@ -394,12 +390,13 @@ protected:
         typename Vector::field_type gLmin = rachfordRice_g_(K, Lmin, z);
       
         // Print new header
-        if (verbosity == 3 || verbosity == 4) {
+        if (verbosity >= 3) {
                 std::cout << std::setw(10) << "Iteration" << std::setw(16) << "g(Lmid)" << std::setw(16) << "L" << std::endl;
         }
-            
+
+        constexpr int max_it = 100;
         // Bisection loop
-        for (int iteration=1; iteration<100; ++iteration){
+        for (int iteration = 0; iteration < max_it; ++iteration){
             // New midpoint
             auto L = (Lmin + Lmax) / 2;
             auto gMid = rachfordRice_g_(K, L, z);
@@ -423,7 +420,7 @@ protected:
                 gLmin = gMid;
             }
         }
-        throw std::runtime_error(" Rachford-Rice with bisection failed!");
+        throw std::runtime_error(" Rachford-Rice with bisection failed with " + std::to_string(max_it) + " iterations!");
     }
 
     template <class FlashFluidState, class ComponentVector>
@@ -599,9 +596,6 @@ protected:
         }
     }
 
-    // TODO: refactoring the template parameter later
-    // For the function flash_2ph, we should carry the derivatives in, and
-    // return with the correct and updated derivatives
     template <class FluidState, class ComponentVector>
     static void flash_2ph(const ComponentVector& z_scalar,
                           const std::string& flash_2p_method,
@@ -653,7 +647,9 @@ protected:
 
         // Compute x and y from K, L and Z
         computeLiquidVapor_(fluid_state, L, K, z);
-        std::cout << " the current L is " << Opm::getValue(L) << std::endl;
+        if (verbosity >= 1) {
+            std::cout << " the current L is " << Opm::getValue(L) << std::endl;
+        }
 
         // Print initial condition
         if (verbosity >= 1) {
@@ -732,7 +728,9 @@ protected:
         while (iter < max_iter) {
             assembleNewton_<CompositionalFluidState<Eval, FluidSystem>, ComponentVector, num_primary_variables, num_equations>
                     (flash_fluid_state, z, jac, res);
-            std::cout << " newton residual is " << res.two_norm() << std::endl;
+            if (verbosity >= 1) {
+                std::cout << " newton residual is " << res.two_norm() << std::endl;
+            }
             converged = res.two_norm() < tolerance;
             if (converged) {
                 break;
@@ -763,15 +761,17 @@ protected:
                 }
             }
         }
-        for (unsigned i = 0; i < num_equations; ++i) {
+        if (verbosity >= 1) {
-            for (unsigned j = 0; j < num_primary_variables; ++j) {
+            for (unsigned i = 0; i < num_equations; ++i) {
-                std::cout << " " << jac[i][j] ;
+                for (unsigned j = 0; j < num_primary_variables; ++j) {
+                    std::cout << " " << jac[i][j];
+                }
+                std::cout << std::endl;
             }
             std::cout << std::endl;
         }
-        std::cout << std::endl;
         if (!converged) {
-            throw std::runtime_error(" Newton composition update did not converge within maxIterations");
+            throw std::runtime_error(" Newton composition update did not converge within maxIterations " + std::to_string(max_iter));
         }
 
         // fluid_state is scalar, we need to update all the values for fluid_state here

opm/material/eos/PengRobinsonMixture.hpp

@@ -86,14 +86,11 @@ public:
       * 4th edition, McGraw-Hill, 1987, pp. 42-44, 143-145
       */
     template <class FluidState, class Params, class LhsEval = typename FluidState::Scalar>
-    static LhsEval computeFugacityCoefficient(const FluidState& fs_arg,
+    static LhsEval computeFugacityCoefficient(const FluidState& fs,
-                                              const Params& params_arg,
+                                              const Params& params,
                                               unsigned phaseIdx,
                                               unsigned compIdx)
     {
-        auto fs = fs_arg;
-        auto params = params_arg;
-
         // 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.

opm/material/eos/PengRobinsonParamsMixture.hpp

@@ -70,7 +70,7 @@ public:
     /*!
      * \brief TODO
      */
-    Scalar getaCache(unsigned compIIdx, unsigned compJIdx )
+    Scalar getaCache(unsigned compIIdx, unsigned compJIdx ) const
     {
         return aCache_[compIIdx][compJIdx];
     }

opm/material/fluidsystems/PTFlashParameterCache.hpp

@@ -2,7 +2,8 @@
 // vi: set et ts=4 sw=4 sts=4:
 /*
   Copyright 2022 NORCE.
-  
+  Copyright 2022 SINTEF Digital, Mathematics and Cybernetics.
+
   This file is part of the Open Porous Media project (OPM).
 
   OPM is free software: you can redistribute it and/or modify
@@ -73,16 +74,8 @@ public:
                      unsigned phaseIdx,
                      int exceptQuantities = ParentType::None)
     {
-        // if (phaseIdx != oilPhaseIdx)
-        //     return;
-
         updateEosParams(fluidState, phaseIdx, exceptQuantities);
 
-        // if we don't need to recalculate the molar volume, we exit
-        // here
-        // if (VmUpToDate_[phaseIdx])
-        //     return;
-
         // update the phase's molar volume
         updateMolarVolume_(fluidState, phaseIdx);
     }
@@ -117,7 +110,7 @@ public:
         case gasPhaseIdx: return gasPhaseParams_.a();
         default:
             throw std::logic_error("The a() parameter is only defined for "
-                                   "oil phases");
+                                   "oil and gas phases");
         };
     }
 
@@ -134,7 +127,7 @@ public:
         case gasPhaseIdx: return gasPhaseParams_.b();
         default:
             throw std::logic_error("The b() parameter is only defined for "
-                                   "oil phase");
+                                   "oil and gas phase");
         };
     }
 
@@ -154,7 +147,7 @@ public:
         case gasPhaseIdx: return gasPhaseParams_.pureParams(compIdx).a();
         default:
             throw std::logic_error("The a() parameter is only defined for "
-                                   "oil phase");
+                                   "oil and gas phase");
         };
     }
 
@@ -173,7 +166,7 @@ public:
         case gasPhaseIdx: return gasPhaseParams_.pureParams(compIdx).b();
         default:
             throw std::logic_error("The b() parameter is only defined for "
-                                   "oil phase");
+                                   "oil and gas phase");
         };
     }
 
@@ -186,7 +179,7 @@ public:
      * \param compIdx The component phase of interest
      * \param compJIdx Additional component index
      */
-    Scalar aCache(unsigned phaseIdx, unsigned compIdx, unsigned compJIdx)
+    Scalar aCache(unsigned phaseIdx, unsigned compIdx, unsigned compJIdx) const
     {
         switch (phaseIdx)
         {
@@ -194,7 +187,7 @@ public:
         case gasPhaseIdx: return gasPhaseParams_.getaCache(compIdx,compJIdx);
         default:
             throw std::logic_error("The aCache() parameter is only defined for "
-                                   "oil phase");
+                                   "oil and gas phase");
         };
     }
 
@@ -235,9 +228,6 @@ public:
                          unsigned phaseIdx,
                          int exceptQuantities = ParentType::None)
     {
-        // if (phaseIdx != oilPhaseIdx)
-        //     return;
-
         if (!(exceptQuantities & ParentType::Temperature))
         {
             updatePure_(fluidState, phaseIdx);

opm/material/fluidsystems/Spe5ParameterCache.hpp

@@ -182,7 +182,7 @@ public:
      * \param compIdx The component phase of interest
      * \param compJIdx Additional component index
      */
-    Scalar aCache(unsigned phaseIdx, unsigned compIdx, unsigned compJIdx)
+    Scalar aCache(unsigned phaseIdx, unsigned compIdx, unsigned compJIdx) const
     {
         switch (phaseIdx)
         {

opm/material/fluidsystems/ThreeComponentFluidSystem.hh

@@ -181,7 +181,6 @@ namespace Opm {
 
             LhsEval dens;
             if (phaseIdx == oilPhaseIdx || phaseIdx == gasPhaseIdx) {
-                // paramCache.updatePhase(fluidState, phaseIdx);
                 dens = fluidState.averageMolarMass(phaseIdx) / paramCache.molarVolume(phaseIdx);
             }
             return dens;
@@ -210,11 +209,7 @@ namespace Opm {
             assert(phaseIdx < numPhases);
             assert(compIdx < numComponents);
 
-            // TODO: here the derivatives for the phi are dropped. Should we keep the derivatives against the pressure
+            LhsEval phi = PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
-            // and temperature?
-            LhsEval phi = PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx); 
-            //Scalar phi = Opm::getValue(
-            //        PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
             return phi;
         }
 

tests/test_co2brine_ptflash.cpp

@@ -127,7 +127,7 @@ bool testPTFlash(const std::string& flash_twophase_method)
 
     const int spatialIdx = 0;
     using Flash = Opm::PTFlash<double, FluidSystem>;
-    Flash::solve(fluid_state, z, spatialIdx, flash_verbosity, flash_twophase_method, flash_tolerance);
+    Flash::solve(fluid_state, z, spatialIdx, flash_twophase_method, flash_tolerance, flash_verbosity);
 
     return result_okay(fluid_state);
 }
@@ -229,4 +229,4 @@ int main(int argc, char **argv)
     }
 
     return 0;
-}
+}

tests/test_threecomponents_ptflash.cpp

@@ -117,7 +117,7 @@ bool testPTFlash(const std::string& flash_twophase_method)
     }
 
     const double flash_tolerance = 1.e-12; // just to test the setup in co2-compositional
-    const int flash_verbosity = 1;
+    const int flash_verbosity = 0;
 
     // TODO: should we set these?
     // Set initial K and L
@@ -130,7 +130,7 @@ bool testPTFlash(const std::string& flash_twophase_method)
 
     const int spatialIdx = 0;
     using Flash = Opm::PTFlash<double, FluidSystem>;
-    Flash::solve(fluid_state, z, spatialIdx, flash_verbosity, flash_twophase_method, flash_tolerance);
+    Flash::solve(fluid_state, z, spatialIdx, flash_twophase_method, flash_tolerance, flash_verbosity);
 
     return result_okay(fluid_state);
 }
@@ -242,7 +242,9 @@ int main(int argc, char **argv)
     }
 
     if (!test_passed) {
-        throw std::runtime_error(" PTFlash tests failed");
+        throw std::runtime_error(" test_threecomponents_ptflash tests FAILED");
+    } else {
+        std::cout << "test_threecomponents_ptflash testing is SUCCESSFUL";
     }
 
     return 0;