diff --git a/examples/problems/fingerproblem.hh b/examples/problems/fingerproblem.hh index 0eaff67fd..2514060d0 100644 --- a/examples/problems/fingerproblem.hh +++ b/examples/problems/fingerproblem.hh @@ -177,6 +177,7 @@ class FingerProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) enum { // number of phases + numPhases = FluidSystem::numPhases, // phase indices wettingPhaseIdx = FluidSystem::wettingPhaseIdx, @@ -520,6 +521,14 @@ private: Scalar pn = 1e5; fs.setPressure(nonWettingPhaseIdx, pn); fs.setPressure(wettingPhaseIdx, pn); + + typename FluidSystem::template ParameterCache paramCache; + paramCache.updateAll(fs); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + fs.setDensity(phaseIdx, FluidSystem::density(fs, paramCache, phaseIdx)); + fs.setViscosity(phaseIdx, FluidSystem::viscosity(fs, paramCache, phaseIdx)); + } + } DimMatrix K_; diff --git a/examples/problems/fractureproblem.hh b/examples/problems/fractureproblem.hh index c9194de70..e37297be9 100644 --- a/examples/problems/fractureproblem.hh +++ b/examples/problems/fractureproblem.hh @@ -470,6 +470,15 @@ public: fluidState.setPressure(wettingPhaseIdx, 1e5); fluidState.setPressure(nonWettingPhaseIdx, fluidState.pressure(wettingPhaseIdx)); + typename FluidSystem::template ParameterCache paramCache; + paramCache.updateAll(fluidState); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + fluidState.setDensity(phaseIdx, + FluidSystem::density(fluidState, paramCache, phaseIdx)); + fluidState.setViscosity(phaseIdx, + FluidSystem::viscosity(fluidState, paramCache, phaseIdx)); + } + // set a free flow (i.e. Dirichlet) boundary values.setFreeFlow(context, spaceIdx, timeIdx, fluidState); } diff --git a/examples/problems/groundwaterproblem.hh b/examples/problems/groundwaterproblem.hh index 00708b385..e0a6a2357 100644 --- a/examples/problems/groundwaterproblem.hh +++ b/examples/problems/groundwaterproblem.hh @@ -128,10 +128,13 @@ class GroundWaterProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; // copy some indices for convenience typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; enum { + numPhases = FluidSystem::numPhases, + // Grid and world dimension dim = GridView::dimension, dimWorld = GridView::dimensionworld, @@ -141,7 +144,6 @@ class GroundWaterProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) }; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; - typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector; typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector; @@ -320,6 +322,13 @@ public: fs.setPressure(/*phaseIdx=*/0, pressure); fs.setTemperature(T); + typename FluidSystem::template ParameterCache paramCache; + paramCache.updateAll(fs); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + fs.setDensity(phaseIdx, FluidSystem::density(fs, paramCache, phaseIdx)); + fs.setViscosity(phaseIdx, FluidSystem::viscosity(fs, paramCache, phaseIdx)); + } + // impose an freeflow boundary condition values.setFreeFlow(context, spaceIdx, timeIdx, fs); } diff --git a/examples/problems/infiltrationproblem.hh b/examples/problems/infiltrationproblem.hh index 3cdaeb896..655878c8a 100644 --- a/examples/problems/infiltrationproblem.hh +++ b/examples/problems/infiltrationproblem.hh @@ -449,7 +449,7 @@ private: typedef Opm::ComputeFromReferencePhase CFRP; typename FluidSystem::template ParameterCache paramCache; CFRP::solve(fs, paramCache, gasPhaseIdx, - /*setViscosity=*/false, + /*setViscosity=*/true, /*setEnthalpy=*/false); fs.setMoleFraction(waterPhaseIdx, H2OIdx, diff --git a/examples/problems/lensproblem.hh b/examples/problems/lensproblem.hh index 928201b18..21d975541 100644 --- a/examples/problems/lensproblem.hh +++ b/examples/problems/lensproblem.hh @@ -426,9 +426,9 @@ public: const GlobalPosition& pos = context.pos(spaceIdx, timeIdx); if (onLeftBoundary_(pos) || onRightBoundary_(pos)) { - // free flow boundary - Scalar densityW = WettingPhase::density(temperature_, - /*pressure=*/Scalar(1e5)); + // free flow boundary. we assume incompressible fluids + Scalar densityW = WettingPhase::density(temperature_, /*pressure=*/Scalar(1e5)); + Scalar densityN = NonwettingPhase::density(temperature_, /*pressure=*/Scalar(1e5)); Scalar T = temperature(context, spaceIdx, timeIdx); Scalar pw, Sw; @@ -465,6 +465,12 @@ public: fs.setPressure(wettingPhaseIdx, pw); fs.setPressure(nonWettingPhaseIdx, pw + pC[nonWettingPhaseIdx] - pC[wettingPhaseIdx]); + fs.setDensity(wettingPhaseIdx, densityW); + fs.setDensity(nonWettingPhaseIdx, densityN); + + fs.setViscosity(wettingPhaseIdx, WettingPhase::viscosity(temperature_, fs.pressure(wettingPhaseIdx))); + fs.setViscosity(nonWettingPhaseIdx, NonwettingPhase::viscosity(temperature_, fs.pressure(nonWettingPhaseIdx))); + // impose an freeflow boundary condition values.setFreeFlow(context, spaceIdx, timeIdx, fs); } diff --git a/examples/problems/obstacleproblem.hh b/examples/problems/obstacleproblem.hh index 0fc665963..7c0847bce 100644 --- a/examples/problems/obstacleproblem.hh +++ b/examples/problems/obstacleproblem.hh @@ -528,12 +528,11 @@ private: // make the fluid state consistent with local thermodynamic // equilibrium - typedef Opm::ComputeFromReferencePhase - ComputeFromReferencePhase; + typedef Opm::ComputeFromReferencePhase ComputeFromReferencePhase; typename FluidSystem::template ParameterCache paramCache; ComputeFromReferencePhase::solve(fs, paramCache, refPhaseIdx, - /*setViscosity=*/false, + /*setViscosity=*/true, /*setEnthalpy=*/false); } diff --git a/examples/problems/outflowproblem.hh b/examples/problems/outflowproblem.hh index 6449a1b3d..5d4dda9c1 100644 --- a/examples/problems/outflowproblem.hh +++ b/examples/problems/outflowproblem.hh @@ -122,6 +122,8 @@ class OutflowProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) dim = GridView::dimension, dimWorld = GridView::dimensionworld, + numPhases = FluidSystem::numPhases, + // component indices H2OIdx = FluidSystem::H2OIdx, N2Idx = FluidSystem::N2Idx @@ -267,6 +269,13 @@ public: fs.setMoleFraction(/*phaseIdx=*/0, N2Idx, xlN2); fs.setMoleFraction(/*phaseIdx=*/0, H2OIdx, 1 - xlN2); + typename FluidSystem::template ParameterCache paramCache; + paramCache.updateAll(fs); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + fs.setDensity(phaseIdx, FluidSystem::density(fs, paramCache, phaseIdx)); + fs.setViscosity(phaseIdx, FluidSystem::viscosity(fs, paramCache, phaseIdx)); + } + // impose an freeflow boundary condition values.setFreeFlow(context, spaceIdx, timeIdx, fs); } @@ -343,6 +352,13 @@ private: fs.setMoleFraction(/*phaseIdx=*/0, H2OIdx, 1.0); fs.setMoleFraction(/*phaseIdx=*/0, N2Idx, 0); fs.setTemperature(T); + + typename FluidSystem::template ParameterCache paramCache; + paramCache.updateAll(fs); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + fs.setDensity(phaseIdx, FluidSystem::density(fs, paramCache, phaseIdx)); + fs.setViscosity(phaseIdx, FluidSystem::viscosity(fs, paramCache, phaseIdx)); + } } const Scalar eps_; diff --git a/examples/problems/powerinjectionproblem.hh b/examples/problems/powerinjectionproblem.hh index bf59a6b2c..ed91e577b 100644 --- a/examples/problems/powerinjectionproblem.hh +++ b/examples/problems/powerinjectionproblem.hh @@ -170,6 +170,7 @@ class PowerInjectionProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) enum { // number of phases + numPhases = FluidSystem::numPhases, // phase indices wettingPhaseIdx = FluidSystem::wettingPhaseIdx, @@ -409,6 +410,15 @@ private: Scalar p = 1e5; initialFluidState_.setPressure(wettingPhaseIdx, p); initialFluidState_.setPressure(nonWettingPhaseIdx, p); + + typename FluidSystem::template ParameterCache paramCache; + paramCache.updateAll(initialFluidState_); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + initialFluidState_.setDensity(phaseIdx, + FluidSystem::density(initialFluidState_, paramCache, phaseIdx)); + initialFluidState_.setViscosity(phaseIdx, + FluidSystem::viscosity(initialFluidState_, paramCache, phaseIdx)); + } } DimMatrix K_; diff --git a/examples/problems/reservoirproblem.hh b/examples/problems/reservoirproblem.hh index 7d0f748d4..7ad661dd7 100644 --- a/examples/problems/reservoirproblem.hh +++ b/examples/problems/reservoirproblem.hh @@ -663,7 +663,7 @@ private: CFRP::solve(injFs, paramCache, /*refPhaseIdx=*/waterPhaseIdx, - /*setViscosities=*/false, + /*setViscosities=*/true, /*setEnthalpies=*/false); // set up the fluid state used for the producer @@ -681,7 +681,7 @@ private: CFRP::solve(prodFs, paramCache, /*refPhaseIdx=*/oilPhaseIdx, - /*setViscosities=*/false, + /*setViscosities=*/true, /*setEnthalpies=*/false); } diff --git a/examples/problems/richardslensproblem.hh b/examples/problems/richardslensproblem.hh index 014d7b66b..2152cabdd 100644 --- a/examples/problems/richardslensproblem.hh +++ b/examples/problems/richardslensproblem.hh @@ -368,6 +368,14 @@ public: fs.setPressure(wettingPhaseIdx, pnRef_ + pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]); fs.setPressure(nonWettingPhaseIdx, pnRef_); + typename FluidSystem::template ParameterCache paramCache; + paramCache.updateAll(fs); + fs.setDensity(wettingPhaseIdx, FluidSystem::density(fs, paramCache, wettingPhaseIdx)); + //fs.setDensity(nonWettingPhaseIdx, FluidSystem::density(fs, paramCache, nonWettingPhaseIdx)); + + fs.setViscosity(wettingPhaseIdx, FluidSystem::viscosity(fs, paramCache, wettingPhaseIdx)); + //fs.setViscosity(nonWettingPhaseIdx, FluidSystem::viscosity(fs, paramCache, nonWettingPhaseIdx)); + values.setFreeFlow(context, spaceIdx, timeIdx, fs); } else if (onInlet_(pos)) { diff --git a/examples/problems/waterairproblem.hh b/examples/problems/waterairproblem.hh index e0697df07..e152a4be3 100644 --- a/examples/problems/waterairproblem.hh +++ b/examples/problems/waterairproblem.hh @@ -527,7 +527,7 @@ private: typename FluidSystem::template ParameterCache paramCache; typedef Opm::ComputeFromReferencePhase CFRP; - CFRP::solve(fs, paramCache, liquidPhaseIdx, /*setViscosity=*/false, /*setEnthalpy=*/true); + CFRP::solve(fs, paramCache, liquidPhaseIdx, /*setViscosity=*/true, /*setEnthalpy=*/true); } void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro) diff --git a/examples/tutorial1problem.hh b/examples/tutorial1problem.hh index 997e93ce7..e76e1e67c 100644 --- a/examples/tutorial1problem.hh +++ b/examples/tutorial1problem.hh @@ -243,6 +243,13 @@ public: fs.setPressure(wettingPhaseIdx, 200e3); fs.setPressure(nonWettingPhaseIdx, 200e3 + pC[nonWettingPhaseIdx] - pC[nonWettingPhaseIdx]); + typename FluidSystem::template ParameterCache paramCache; + paramCache.updateAll(fs); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + fs.setDensity(phaseIdx, FluidSystem::density(fs, paramCache, phaseIdx)); + fs.setViscosity(phaseIdx, FluidSystem::viscosity(fs, paramCache, phaseIdx)); + } + values.setFreeFlow(context, spaceIdx, timeIdx, fs); } else if (pos[0] > this->boundingBoxMax()[0] - eps_) {