updating based on review comments.

examples/problems/co2ptflashproblem.hh

@@ -147,7 +147,7 @@ struct EnableGravity<TypeTag, TTag::CO2PTBaseProblem> { static constexpr bool va
 template <class TypeTag>
 struct Initialpressure<TypeTag, TTag::CO2PTBaseProblem> {
     using type = GetPropType<TypeTag, Scalar>;
-    static constexpr type value = 1e5;
+    static constexpr type value = 75.e5;
 };
 
 template <class TypeTag>
@@ -278,8 +278,6 @@ struct CellsY<TypeTag, TTag::CO2PTBaseProblem> { static constexpr int value = 1;
 template<class TypeTag>
 struct CellsZ<TypeTag, TTag::CO2PTBaseProblem> { static constexpr int value = 1; };
 
-
-// compositional, with diffusion
 template <class TypeTag>
 struct EnableEnergy<TypeTag, TTag::CO2PTBaseProblem> {
     static constexpr bool value = false;
@@ -357,7 +355,6 @@ public:
         K_ = this->toDimMatrix_(9.869232667160131e-14);
 
         porosity_ = 0.1;
-    
     }
 
     template <class Context>
@@ -558,9 +555,8 @@ private:
         ComponentVector sat;
         sat[0] = 1.0;
         sat[1] = 1.0 - sat[0];
-        const Scalar temp = 423.25;
 
-        Scalar p0 = 75e5; // CONVERGENCE FAILURE WITH 75
+        Scalar p0 = EWOMS_GET_PARAM(TypeTag, Scalar, Initialpressure);
 
         //\Note, for an AD variable, if we multiply it with 2, the derivative will also be scalced with 2,
         //\Note, so we should not do it.
@@ -587,7 +583,7 @@ private:
         fs.setSaturation(FluidSystem::oilPhaseIdx, sat[0]);
         fs.setSaturation(FluidSystem::gasPhaseIdx, sat[1]);
 
-        fs.setTemperature(temp);
+        fs.setTemperature(temperature_);
 
         // ParameterCache paramCache;
         {

opm/models/ptflash/flashintensivequantities.hh

@@ -135,7 +135,6 @@ public:
 
         // Get initial K and L from storage initially (if enabled)
         const auto *hint = elemCtx.thermodynamicHint(dofIdx, timeIdx);
-        const auto *hint2 = elemCtx.thermodynamicHint(dofIdx, 1);
         if (hint) {
              for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                  const Evaluation& Ktmp = hint->fluidState().K(compIdx);
@@ -144,7 +143,9 @@ public:
              const Evaluation& Ltmp = hint->fluidState().L();
              fluidState_.setLvalue(Ltmp);
         }
-        else if (hint2) {
+        else if (timeIdx == 0 && elemCtx.thermodynamicHint(dofIdx, 1)) {
+             // checking the storage cache
+             const auto& hint2 = elemCtx.thermodynamicHint(dofIdx, 1);
              for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                  const Evaluation& Ktmp = hint2->fluidState().K(compIdx);
                  fluidState_.setKvalue(compIdx, Ktmp);
@@ -205,6 +206,7 @@ public:
 
 
         // Update saturation
+        // \Note: the current implementation assume oil-gas system.
         Evaluation L = fluidState_.L();
         Evaluation So = Opm::max((L * Z_L / ( L * Z_L + (1 - L) * Z_V)), 0.0);
         Evaluation Sg = Opm::max(1 - So, 0.0);

opm/models/ptflash/flashmodel.hh

@@ -284,40 +284,6 @@ public:
         return oss.str();
     }
 
-    /*!
-     * \copydoc FvBaseDiscretization::primaryVarWeight
-     */
-    Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
-    {
-        Scalar tmp = EnergyModule::primaryVarWeight(*this, globalDofIdx, pvIdx);
-        if (tmp > 0)
-            return tmp;
-
-        unsigned compIdx = pvIdx - Indices::conti0EqIdx;
-
-        // make all kg equal. also, divide the weight of all total
-        // compositions by 100 to make the relative errors more
-        // comparable to the ones of the other models (at 10% porosity
-        // the medium is fully saturated with water at atmospheric
-        // conditions if 100 kg/m^3 are present!)
-        return FluidSystem::molarMass(compIdx) / 100.0;
-    }
-
-    /*!
-     * \copydoc FvBaseDiscretization::eqWeight
-     */
-    Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
-    {
-        Scalar tmp = EnergyModule::eqWeight(*this, globalDofIdx, eqIdx);
-        if (tmp > 0)
-            return tmp;
-
-        unsigned compIdx = eqIdx - Indices::conti0EqIdx;
-
-        // make all kg equal
-        return FluidSystem::molarMass(compIdx);
-    }
-
     void registerOutputModules_()
     {
         ParentType::registerOutputModules_();

opm/models/ptflash/flashnewtonmethod.hh

@@ -91,9 +91,12 @@ protected:
         // Pressure updates
         ////
         // limit pressure reference change to 20% of the total value per iteration
+        constexpr Scalar max_percent_change = 0.2;
+        constexpr Scalar upper_bound = 1. + max_percent_change;
+        constexpr Scalar lower_bound = 1. - max_percent_change;
         clampValue_(nextValue[pressure0Idx],
-                    currentValue[pressure0Idx]*0.8,
+                    currentValue[pressure0Idx] * lower_bound,
-                    currentValue[pressure0Idx]*1.2);
+                    currentValue[pressure0Idx] * upper_bound);
 
         ////
         // z updates