Merge pull request #276 from andlaus/boundary_condition_refactoring

refactor the boundary condition handling slightly
This commit is contained in:
Andreas Lauser 2018-01-22 13:02:54 +01:00 committed by GitHub
commit 922036da89
12 changed files with 84 additions and 11 deletions

View File

@ -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<Scalar> 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_;

View File

@ -470,6 +470,15 @@ public:
fluidState.setPressure(wettingPhaseIdx, 1e5);
fluidState.setPressure(nonWettingPhaseIdx, fluidState.pressure(wettingPhaseIdx));
typename FluidSystem::template ParameterCache<Scalar> 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);
}

View File

@ -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<Scalar> 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);
}

View File

@ -449,7 +449,7 @@ private:
typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem> CFRP;
typename FluidSystem::template ParameterCache<Scalar> paramCache;
CFRP::solve(fs, paramCache, gasPhaseIdx,
/*setViscosity=*/false,
/*setViscosity=*/true,
/*setEnthalpy=*/false);
fs.setMoleFraction(waterPhaseIdx, H2OIdx,

View File

@ -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);
}

View File

@ -528,12 +528,11 @@ private:
// make the fluid state consistent with local thermodynamic
// equilibrium
typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem>
ComputeFromReferencePhase;
typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem> ComputeFromReferencePhase;
typename FluidSystem::template ParameterCache<Scalar> paramCache;
ComputeFromReferencePhase::solve(fs, paramCache, refPhaseIdx,
/*setViscosity=*/false,
/*setViscosity=*/true,
/*setEnthalpy=*/false);
}

View File

@ -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<Scalar> 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<Scalar> 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_;

View File

@ -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<Scalar> 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_;

View File

@ -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);
}

View File

@ -368,6 +368,14 @@ public:
fs.setPressure(wettingPhaseIdx, pnRef_ + pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]);
fs.setPressure(nonWettingPhaseIdx, pnRef_);
typename FluidSystem::template ParameterCache<Scalar> 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)) {

View File

@ -527,7 +527,7 @@ private:
typename FluidSystem::template ParameterCache<Scalar> paramCache;
typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem> 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)

View File

@ -243,6 +243,13 @@ public:
fs.setPressure(wettingPhaseIdx, 200e3);
fs.setPressure(nonWettingPhaseIdx, 200e3 + pC[nonWettingPhaseIdx] - pC[nonWettingPhaseIdx]);
typename FluidSystem::template ParameterCache<Scalar> 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_) {