From 0406d6780fd09bcce5c27671b77d59e6ef33aa5d Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Mon, 22 Jan 2018 10:33:55 +0100 Subject: [PATCH] refactor the boundary condition handling slightly instead of passing a "minimal" fluid state that defines the thermodynamic conditions on the domain boundary and the models calculating everything they need based on this, it is now assumed that all quantities needed by the code that computes the boundary fluxes are defined. This simplifies the boundary flux computation code, it allows to get rid of the `paramCache` argument for these methods and to potentially speed things up because quantities do not get re-calculated unconditionally. on the flipside, this requires slightly more effort to define the conditions at the boundary on the problem level and it makes it less obvious which quantities are actually used. That said, one now has the freedom to shoot oneself into the foot more easily when specifying boundary conditions and also tools like valgrind or ASAN will normally complain about undefined quantities if this happens. --- examples/problems/fingerproblem.hh | 9 +++++++++ examples/problems/fractureproblem.hh | 9 +++++++++ examples/problems/groundwaterproblem.hh | 11 ++++++++++- examples/problems/infiltrationproblem.hh | 2 +- examples/problems/lensproblem.hh | 12 +++++++++--- examples/problems/obstacleproblem.hh | 5 ++--- examples/problems/outflowproblem.hh | 16 ++++++++++++++++ examples/problems/powerinjectionproblem.hh | 10 ++++++++++ examples/problems/reservoirproblem.hh | 4 ++-- examples/problems/richardslensproblem.hh | 8 ++++++++ examples/problems/waterairproblem.hh | 2 +- examples/tutorial1problem.hh | 7 +++++++ 12 files changed, 84 insertions(+), 11 deletions(-) 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_) {