Merge pull request #5862 from GitPaean/simplying_flash_usage

adapting to the interface change for PTFlash::solve()

examples/problems/co2ptflashproblem.hh

@@ -457,9 +457,6 @@ private:
             comp[1] = Evaluation::createVariable(0.01 - 1e-3, 2);
             comp[2] = 1. - comp[0] - comp[1];
         }
-        ComponentVector sat;
-        sat[0] = 1.0;
-        sat[1] = 1.0 - sat[0];
 
         Scalar p0 = Parameters::Get<Parameters::Initialpressure<Scalar>>();
 
@@ -476,46 +473,10 @@ private:
         fs.setPressure(FluidSystem::oilPhaseIdx, p_init);
         fs.setPressure(FluidSystem::gasPhaseIdx, p_init);
 
-        for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
-            fs.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, comp[compIdx]);
-            fs.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, comp[compIdx]);
-        }
-
-        // It is used here only for calculate the z
-        fs.setSaturation(FluidSystem::oilPhaseIdx, sat[0]);
-        fs.setSaturation(FluidSystem::gasPhaseIdx, sat[1]);
-
         fs.setTemperature(temperature_);
 
-        // ParameterCache paramCache;
-        {
-            typename FluidSystem::template ParameterCache<Evaluation> paramCache;
-            paramCache.updatePhase(fs, FluidSystem::oilPhaseIdx);
-            paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx);
-            fs.setDensity(FluidSystem::oilPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::oilPhaseIdx));
-            fs.setDensity(FluidSystem::gasPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx));
-            fs.setViscosity(FluidSystem::oilPhaseIdx, FluidSystem::viscosity(fs, paramCache, FluidSystem::oilPhaseIdx));
-            fs.setViscosity(FluidSystem::gasPhaseIdx, FluidSystem::viscosity(fs, paramCache, FluidSystem::gasPhaseIdx));
-        }
-
-        // determine the component total fractions
-        // TODO: duplicated code here, while should be refactored out when we swithing
-        // to starting from total mole fractions
-        Dune::FieldVector<Scalar, numComponents> z(0.0);
-        Scalar sumMoles = 0.0;
-        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-            const auto saturation = getValue(fs.saturation(phaseIdx));
-            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
-                Scalar tmp = getValue(fs.molarity(phaseIdx, compIdx)) * saturation;
-                tmp = max(tmp, 1.e-8);
-                z[compIdx] += tmp;
-                sumMoles += tmp;
-            }
-        }
-        z /= sumMoles;
-
         for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
-            fs.setMoleFraction(compIdx, z[compIdx]);
+            fs.setMoleFraction(compIdx, comp[compIdx]);
         }
 
         // Set initial K and L

opm/models/ptflash/flashintensivequantities.hh

@@ -178,7 +178,7 @@ public:
             const int spatialIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
             std::cout << " updating the intensive quantities for Cell " << spatialIdx << std::endl;
         }
-        FlashSolver::solve(fluidState_, z, flashTwoPhaseMethod, flashTolerance, flashVerbosity);
+        FlashSolver::solve(fluidState_, flashTwoPhaseMethod, flashTolerance, flashVerbosity);
 
         if (flashVerbosity >= 5) {
             // printing of flash result after solve

opm/models/ptflash/flashmodel.hh

@@ -231,7 +231,7 @@ public:
             ("Method for solving vapor-liquid composition. Available options include: "
              "ssi, newton, ssi+newton");
 
-        Parameters::SetDefault<Parameters::FlashTolerance<Scalar>>(1e-12);
+        Parameters::SetDefault<Parameters::FlashTolerance<Scalar>>(1.e-8);
         Parameters::SetDefault<Parameters::EnableIntensiveQuantityCache>(true);
 
         // since thermodynamic hints are basically free if the cache for intensive quantities is

opm/simulators/flow/FlowProblemComp.hpp

@@ -324,56 +324,54 @@ public:
         const unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
         const auto& initial_fs = initialFluidStates_[globalDofIdx];
         Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
-        // TODO: the current approach is assuming we begin with XMF and YMF.
-        // TODO: maybe we should make it begin with ZMF
-        using ComponentVector = Dune::FieldVector<Scalar, numComponents>;
         for (unsigned p = 0; p < numPhases; ++p) { // TODO: assuming the phaseidx continuous
-            ComponentVector vals;
-            auto& last_eval = vals[numComponents - 1];
-            last_eval = 1.;
-            for (unsigned c = 0; c < numComponents - 1; ++c) {
-                const auto val = initial_fs.moleFraction(p, c);
-                vals[c] = val;
-                last_eval -= val;
-            }
-            for (unsigned c = 0; c < numComponents; ++c) {
-                fs.setMoleFraction(p, c, vals[c]);
-            }
-
             // pressure
-            const auto p_val = initial_fs.pressure(p);
+            fs.setPressure(p, initial_fs.pressure(p));
-            fs.setPressure(p, p_val);
 
-            const auto sat_val = initial_fs.saturation(p);
+            // saturation
-            fs.setSaturation(p, sat_val);
+            fs.setSaturation(p, initial_fs.saturation(p));
 
-            const auto temp_val = initial_fs.temperature(p);
+            // temperature
-            fs.setTemperature(temp_val);
+            fs.setTemperature(initial_fs.temperature(p));
         }
 
-        {
-            typename FluidSystem::template ParameterCache<Scalar> paramCache;
-            paramCache.updatePhase(fs, FluidSystem::oilPhaseIdx);
-            paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx);
-            fs.setDensity(FluidSystem::oilPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::oilPhaseIdx));
-            fs.setDensity(FluidSystem::gasPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx));
-            fs.setViscosity(FluidSystem::oilPhaseIdx, FluidSystem::viscosity(fs, paramCache, FluidSystem::oilPhaseIdx));
-            fs.setViscosity(FluidSystem::gasPhaseIdx, FluidSystem::viscosity(fs, paramCache, FluidSystem::gasPhaseIdx));
-        }
 
-        // determine the component fractions
+        if (!zmf_initialization_) {
-        Dune::FieldVector<Scalar, numComponents> z(0.0);
+            for (unsigned p = 0; p < numPhases; ++p) {
-        Scalar sumMoles = 0.0;
+                for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
-        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+                    fs.setMoleFraction(p, compIdx, initial_fs.moleFraction(p, compIdx));
-            const auto saturation = getValue(fs.saturation(phaseIdx));
+                }
+            }
+
+            {
+                typename FluidSystem::template ParameterCache<Scalar> paramCache;
+                paramCache.updatePhase(fs, FluidSystem::oilPhaseIdx);
+                paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx);
+                fs.setDensity(FluidSystem::oilPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::oilPhaseIdx));
+                fs.setDensity(FluidSystem::gasPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx));
+            }
+            // determine the component fractions
+            Dune::FieldVector<Scalar, numComponents> z(0.0);
+            Scalar sumMoles = 0.0;
+            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+                const auto saturation = fs.saturation(phaseIdx);
+                for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                    Scalar tmp = fs.molarity(phaseIdx, compIdx) * saturation;
+                    tmp = max(tmp, 1e-8);
+                    z[compIdx] += tmp;
+                    sumMoles += tmp;
+                }
+            }
+            z /= sumMoles;
             for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
-                Scalar tmp = getValue(fs.molarity(phaseIdx, compIdx)) * saturation;
+                fs.setMoleFraction(compIdx, z[compIdx]);
-                tmp = max(tmp, 1e-8);
+            }
-                z[compIdx] += tmp;
+        } else {
-                sumMoles += tmp;
+            // TODO: should we normalize the input?
+            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                fs.setMoleFraction(compIdx, initial_fs.moleFraction(compIdx));
             }
         }
-        z /= sumMoles;
 
         // Set initial K and L
         for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
@@ -381,10 +379,6 @@ public:
             fs.setKvalue(compIdx, Ktmp);
         }
 
-        for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
-            fs.setMoleFraction(compIdx, z[compIdx]);
-        }
-
         const Scalar& Ltmp = -1.0;
         fs.setLvalue(Ltmp);
 
@@ -560,6 +554,7 @@ protected:
             }
 
             if (has_zmf) {
+                zmf_initialization_ = true;
                 const auto& zmfData = fp.get_double("ZMF");
                 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                     const std::size_t data_idx = compIdx * numDof + dofIdx;
@@ -595,6 +590,8 @@ private:
 
     std::vector<InitialFluidState> initialFluidStates_;
 
+    bool zmf_initialization_ {false};
+
     bool enableEclOutput_{false};
     std::unique_ptr<EclWriterType> eclWriter_;
 };