Skip to content
GitLab
Commit aa24c42c authored by Martin Schneider's avatar Martin Schneider
Browse files

[2phybrid] First working implementation

parent f6f8f99a
Hide whitespace changes
Inline Side-by-side
dumux/porousmediumflow/2pff/implicit/darcyslaw.hh
View file @ aa24c42c
......@@ -98,7 +98,7 @@ public:
const auto& insideVolVars = elemVolVars[insideScv];
const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
if (eqIdx == transportEqIdx)
if (eqIdx == transportEqIdx || eqIdx == totalvelocityEqIdx)
{
//////////////////////////////////////////////////////////////
// we return three separate fluxes: the viscous, gravity and capillary flux
......@@ -108,8 +108,6 @@ public:
using VelocityVector = Dune::FieldVector<Scalar, GridView::dimensionworld>;
// static const VelocityVector v_t = getParamFromGroup<VelocityVector>(problem.paramGroup(), "Problem.TotalVelocity");
if (eqIdx == transportEqIdx)
{
// viscous Flu
//TODO: calculate the velocity (no more constant)
......@@ -176,13 +174,13 @@ public:
const auto pcOutside = outsideVolVars.capillaryPressure();
Scalar WeightedCapillaryDp = WeightedCapillaryMobility*(pcOutside - pcInside);
static VelocityVector v_t = fluxVarsCache.advectionTij()*(mobT_WA*(pnInside - pnOutside)+ WeightedGravityDp + WeightedCapillaryDp);
fluxes[viscousFluxIdx] = fluxVarsCache.advectionTij()*(mobT_WA*(pnInside - pnOutside)+ WeightedGravityDp + WeightedCapillaryDp);
//////////////////////////////////////////////////////////////
// The viscous flux before upwinding is the total velocity
//////////////////////////////////////////////////////////////
fluxes[viscousFluxIdx] = (v_t*scvf.unitOuterNormal())*scvf.area(); // TODO extrusion factor
//fluxes[viscousFluxIdx] = (v_t*scvf.unitOuterNormal())*scvf.area(); // TODO extrusion factor
//////////////////////////////////////////////////////////////
// Compute the capillary flux before upwinding
......@@ -245,6 +243,6 @@ public:
}
};
}; // end namespace Dumux
} // end namespace Dumux
#endif
dumux/porousmediumflow/2pff/implicit/localresidual.hh
View file @ aa24c42c
......@@ -108,7 +108,7 @@ public:
// this is just a dummy we do upwinding internally in the upwind scheme class
auto upwindTermTransport = [](const auto& volVars, const int phaseIdx = 0){ return 0.0; };
flux[transportEqIdx] = fluxVars.advectiveFlux(transportEqIdx, upwindTermTransport);
flux[totalvelocityEqIdx] = fluxVars.advetiveFlux(totalvelocityEqIdx,upwindTermTransport) +fluxVars.advectiveFlux(transportEqIdx, upwindTermTransport);
flux[totalvelocityEqIdx] = fluxVars.advectiveFlux(totalvelocityEqIdx,upwindTermTransport);
//! Add advective phase energy fluxes. For isothermal model the contribution is zero.
// EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, wPhaseIdx);
......
dumux/porousmediumflow/2pff/implicit/upwindscheme.hh
View file @ aa24c42c
......@@ -48,6 +48,12 @@ class TwoPFractionalFlowUpwindScheme
capillaryFluxIdx = 2
};
enum
{
transportEqIdx = 0,
totalvelocityEqIdx = 1
};
public:
// applies a simple upwind scheme to the precalculated advective flux
template<class FluxVariables, class UpwindTermFunction, class Flux>
......@@ -60,7 +66,6 @@ public:
using MaterialLaw = typename std::decay_t<decltype(fluxVars.problem().spatialParams())>::MaterialLaw;
static const bool useHybridUpwinding = getParamFromGroup<bool>(fluxVars.problem().paramGroup(), "Problem.UseHybridUpwinding");
static const bool enableGravity = getParamFromGroup<bool>(fluxVars.problem().paramGroup(), "Problem.EnableGravity");
static constexpr int transportEqIdx = 0;
if (useHybridUpwinding)
{
......@@ -71,9 +76,9 @@ public:
{
// viscous Flux
const auto mobW_V = std::signbit(flux[viscousFluxIdx]) ? outsideVolVars.mobility(wPhaseIdx)
: insideVolVars.mobility(wPhaseIdx);
: insideVolVars.mobility(wPhaseIdx);
const auto mobN_V = std::signbit(flux[viscousFluxIdx]) ? outsideVolVars.mobility(nPhaseIdx)
: insideVolVars.mobility(nPhaseIdx);
: insideVolVars.mobility(nPhaseIdx);
const auto mobT_V = std::max(mobW_V + mobN_V, 1.0e-30);
const Scalar viscousFlux = mobW_V/mobT_V*flux[viscousFluxIdx];
......@@ -131,25 +136,24 @@ public:
const auto& materialLaws = fluxVars.problem().spatialParams().materialLawParamsAtPos(scv.center());
// compute D_max
unsigned int numIntervals = 10;
for(int k=0; k <= numIntervals; k++)
{
Scalar Sw = S_min + (S_max-S_min)*k/numIntervals;
Scalar mobW = MaterialLaw::krw(materialLaws, Sw)/insideVolVars.viscosity(wPhaseIdx);
Scalar mobN = MaterialLaw::krn(materialLaws, Sw)/insideVolVars.viscosity(nPhaseIdx);
Scalar dPc_dSw = MaterialLaw::dpc_dsw(materialLaws, Sw);
D_max = std::max(D_max, -(mobW*mobN)/(mobW + mobN)*dPc_dSw);
}
// compute D_max
unsigned int numIntervals = 10;
for(int k=0; k <= numIntervals; k++)
{
Scalar Sw = S_min + (S_max-S_min)*k/numIntervals;
Scalar mobW = MaterialLaw::krw(materialLaws, Sw)/insideVolVars.viscosity(wPhaseIdx);
Scalar mobN = MaterialLaw::krn(materialLaws, Sw)/insideVolVars.viscosity(nPhaseIdx);
Scalar dPc_dSw = MaterialLaw::dpc_dsw(materialLaws, Sw);
D_max = std::max(D_max, -(mobW*mobN)/(mobW + mobN)*dPc_dSw);
}
return viscousFlux
+ D_max*flux[capillaryFluxIdx]
+ gravFlux;
+ D_max*flux[capillaryFluxIdx]
+ gravFlux;
}
else if (eqIdx == totalvelocityEqIdx)
{
return flux[viscousFluxIdx];
}
else
{
......@@ -175,122 +179,119 @@ public:
Scalar wPhasepotentialdifference,nPhasepotentialdifference =0;
if (enableGravity)
{
const auto& scvf = fluxVars.scvFace();
const auto& insideScv = fluxVars.fvGeometry.scv(scvf.insideScvIdx());
// do averaging for the density over all neighboring elements
const auto rhow = scvf.boundary() ? outsideVolVars.density(wPhaseIdx)
: (insideVolVars.density(wPhaseIdx) + outsideVolVars.density(wPhaseIdx))*0.5;
const auto rhon = scvf.boundary() ? outsideVolVars.density(nPhaseIdx)
: (insideVolVars.density(nPhaseIdx) + outsideVolVars.density(nPhaseIdx))*0.5;
// Obtain inside and outside pressures
const auto pwInside = insideVolVars.pressure(wPhaseIdx);
const auto pwOutside = outsideVolVars.pressure(wPhaseIdx);
const auto pnInside = insideVolVars.pressure(nPhaseIdx);
const auto pnOutside = outsideVolVars.pressure(nPhaseIdx);
const auto& g = gravityAtPos(scvf.ipGlobal());
const auto xInside = insideScv.center();
const auto xOutside = scvf.boundary() ? scvf.ipGlobal()
: fluxVars.fvGeometry.scv(scvf.outsideScvIdx()).center();
wPhasepotentialdifference = (pwInside - pwOutside) + rhow*g*(xInside-xOutside);
nPhasepotentialdifference = (pnInside - pnOutside) + rhon*g*(xInside-xOutside);
}
else
{
// Obtain inside and outside pressures
const auto pwInside = insideVolVars.pressure(wPhaseIdx);
const auto pwOutside = outsideVolVars.pressure(wPhaseIdx);
const auto pnInside = insideVolVars.pressure(nPhaseIdx);
const auto pnOutside = outsideVolVars.pressure(nPhaseIdx);
wPhasepotentialdifference = (pwInside - pwOutside);
nPhasepotentialdifference = (pnInside - pnOutside);
Scalar wPhasepotentialdifference,nPhasepotentialdifference = 0;
if (enableGravity)
{
const auto& scvf = fluxVars.scvFace();
const auto& insideScv = fluxVars.fvGeometry().scv(scvf.insideScvIdx());
// do averaging for the density over all neighboring elements
const auto rhow = scvf.boundary() ? outsideVolVars.density(wPhaseIdx)
: (insideVolVars.density(wPhaseIdx) + outsideVolVars.density(wPhaseIdx))*0.5;
const auto rhon = scvf.boundary() ? outsideVolVars.density(nPhaseIdx)
: (insideVolVars.density(nPhaseIdx) + outsideVolVars.density(nPhaseIdx))*0.5;
}
// Obtain inside and outside pressures
const auto pwInside = insideVolVars.pressure(wPhaseIdx);
const auto pwOutside = outsideVolVars.pressure(wPhaseIdx);
const auto pnInside = insideVolVars.pressure(nPhaseIdx);
const auto pnOutside = outsideVolVars.pressure(nPhaseIdx);
const auto& g = fluxVars.problem().gravityAtPos(scvf.ipGlobal());
const auto xInside = insideScv.center();
const auto xOutside = scvf.boundary() ? scvf.ipGlobal()
: fluxVars.fvGeometry().scv(scvf.outsideScvIdx()).center();
wPhasepotentialdifference = (pwInside - pwOutside) + rhow*(g*(xInside-xOutside));
nPhasepotentialdifference = (pnInside - pnOutside) + rhon*(g*(xInside-xOutside));
//upwind according to the sign of Phasepotentialdifference
const auto mobW = std::signbit(wPhasepotentialdifference) ? outsideVolVars.mobility(wPhaseIdx)
: insideVolVars.mobility(wPhaseIdx);
const auto mobN = std::signbit(nPhasepotentialdifference) ? outsideVolVars.mobility(nPhaseIdx)
: insideVolVars.mobility(nPhaseIdx);
const auto mobT= std::max(mobW + mobN, 1.0e-30);
}
else
{
// Obtain inside and outside pressures
const auto pwInside = insideVolVars.pressure(wPhaseIdx);
const auto pwOutside = outsideVolVars.pressure(wPhaseIdx);
const auto pnInside = insideVolVars.pressure(nPhaseIdx);
const auto pnOutside = outsideVolVars.pressure(nPhaseIdx);
wPhasepotentialdifference = (pwInside - pwOutside);
nPhasepotentialdifference = (pnInside - pnOutside);
}
//upwind according to the sign of Phasepotentialdifference
const auto mobW = std::signbit(wPhasepotentialdifference) ? outsideVolVars.mobility(wPhaseIdx)
: insideVolVars.mobility(wPhaseIdx);
const auto mobN = std::signbit(nPhasepotentialdifference) ? outsideVolVars.mobility(nPhaseIdx)
: insideVolVars.mobility(nPhaseIdx);
const auto mobT= std::max(mobW + mobN, 1.0e-30);
if (eqIdx == transportEqIdx)
{
// viscous Flux
// const auto mobW_V = std::signbit(flux[viscousFluxIdx]) ? outsideVolVars.mobility(wPhaseIdx)
// : insideVolVars.mobility(wPhaseIdx);
// const auto mobN_V = std::signbit(flux[viscousFluxIdx]) ? outsideVolVars.mobility(nPhaseIdx)
// : insideVolVars.mobility(nPhaseIdx);
// const auto mobT_V = std::max(mobW_V + mobN_V, 1.0e-30);
// const auto mobW_V = std::signbit(flux[viscousFluxIdx]) ? outsideVolVars.mobility(wPhaseIdx)
// : insideVolVars.mobility(wPhaseIdx);
// const auto mobN_V = std::signbit(flux[viscousFluxIdx]) ? outsideVolVars.mobility(nPhaseIdx)
// : insideVolVars.mobility(nPhaseIdx);
// const auto mobT_V = std::max(mobW_V + mobN_V, 1.0e-30);
// const Scalar viscousFlux = mobW_V/mobT_V*flux[viscousFluxIdx];
// const Scalar viscousFlux = mobW_V/mobT_V*flux[viscousFluxIdx];
const Scalar viscousFlux = mobW/mobT*flux[viscousFluxIdx];
// Calculate gravity flux
Scalar gravFlux = 0.0;
// if (enableGravity)
// {
// const auto& scvf = fluxVars.scvFace();
// // do averaging for the density over all neighboring elements
// const auto rho = [&](int phaseIdx)
// {
// // boundaries
// if (scvf.boundary())
// return insideVolVars.density(phaseIdx);
//
// // inner faces with two neighboring elements
// else
// return (insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx))*0.5;
// };
// if(rho(nPhaseIdx) - rho(wPhaseIdx) < 0)
// {
// const auto mobW_G = std::signbit(flux[gravityFluxIdx]) ? outsideVolVars.mobility(wPhaseIdx)
// : insideVolVars.mobility(wPhaseIdx);
// const auto mobN_G = std::signbit(flux[gravityFluxIdx]) ? insideVolVars.mobility(nPhaseIdx)
// : outsideVolVars.mobility(nPhaseIdx);
// const auto mobT_G = std::max(mobW_G + mobN_G, 1.0e-30);
// gravFlux = mobW_G*mobN_G/mobT_G * flux[gravityFluxIdx];
// }
// else
// {
// const auto mobW_G = std::signbit(flux[gravityFluxIdx]) ? insideVolVars.mobility(wPhaseIdx)
// : outsideVolVars.mobility(wPhaseIdx);
// const auto mobN_G = std::signbit(flux[gravityFluxIdx]) ? outsideVolVars.mobility(nPhaseIdx)
// : insideVolVars.mobility(nPhaseIdx);
// const auto mobT_G = std::max(mobW_G + mobN_G, 1.0e-30);
// gravFlux = mobW_G*mobN_G/mobT_G * flux[gravityFluxIdx];
// if (enableGravity)
// {
// const auto& scvf = fluxVars.scvFace();
// // do averaging for the density over all neighboring elements
// const auto rho = [&](int phaseIdx)
// {
// // boundaries
// if (scvf.boundary())
// return insideVolVars.density(phaseIdx);
//
// // inner faces with two neighboring elements
// else
// return (insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx))*0.5;
// };
// if(rho(nPhaseIdx) - rho(wPhaseIdx) < 0)
// {
// const auto mobW_G = std::signbit(flux[gravityFluxIdx]) ? outsideVolVars.mobility(wPhaseIdx)
// : insideVolVars.mobility(wPhaseIdx);
// const auto mobN_G = std::signbit(flux[gravityFluxIdx]) ? insideVolVars.mobility(nPhaseIdx)
// : outsideVolVars.mobility(nPhaseIdx);
// const auto mobT_G = std::max(mobW_G + mobN_G, 1.0e-30);
// gravFlux = mobW_G*mobN_G/mobT_G * flux[gravityFluxIdx];
// }
// else
// {
// const auto mobW_G = std::signbit(flux[gravityFluxIdx]) ? insideVolVars.mobility(wPhaseIdx)
// : outsideVolVars.mobility(wPhaseIdx);
// const auto mobN_G = std::signbit(flux[gravityFluxIdx]) ? outsideVolVars.mobility(nPhaseIdx)
// : insideVolVars.mobility(nPhaseIdx);
// const auto mobT_G = std::max(mobW_G + mobN_G, 1.0e-30);
// gravFlux = mobW_G*mobN_G/mobT_G * flux[gravityFluxIdx];
gravFlux = mobW*mobN/mobT * flux[gravityFluxIdx];
// }
// }
// }
// }
// capillary flux
// const auto mobW_C = std::signbit(flux[capillaryFluxIdx]) ? outsideVolVars.mobility(wPhaseIdx)
// : insideVolVars.mobility(wPhaseIdx);
// const auto mobN_C = std::signbit(flux[capillaryFluxIdx]) ? insideVolVars.mobility(nPhaseIdx)
// : outsideVolVars.mobility(nPhaseIdx);
// const auto mobT_C = std::max(mobW_C + mobN_C, 1.0e-30);
// const auto mobW_C = std::signbit(flux[capillaryFluxIdx]) ? outsideVolVars.mobility(wPhaseIdx)
// : insideVolVars.mobility(wPhaseIdx);
// const auto mobN_C = std::signbit(flux[capillaryFluxIdx]) ? insideVolVars.mobility(nPhaseIdx)
// : outsideVolVars.mobility(nPhaseIdx);
// const auto mobT_C = std::max(mobW_C + mobN_C, 1.0e-30);
// const Scalar capillaryFlux = mobW_C*mobN_C/mobT_C*flux[capillaryFluxIdx];
// const Scalar capillaryFlux = mobW_C*mobN_C/mobT_C*flux[capillaryFluxIdx];
const Scalar capillaryFlux = mobW*mobN/mobT*flux[capillaryFluxIdx];
return viscousFlux
+ capillaryFlux
+ gravFlux;
+ capillaryFlux
+ gravFlux;
// for PPU the flux[capillaryFluxIdx] should contain tij*(pc_i - pc_j)
......@@ -304,7 +305,7 @@ Scalar wPhasepotentialdifference,nPhasepotentialdifference =0;
DUNE_THROW(Dune::InvalidStateException, "Unknown equation index!");
}
} // end phase potential upwinding
};
}
};
} // end namespace Dumux
......
test/porousmediumflow/2pff/implicit/constvelproblem.hh
View file @ aa24c42c
......@@ -98,8 +98,11 @@ class ConstVelProblem : public PorousMediumFlowProblem<TypeTag>
enum {
// equation indices
pressureIdx = Indices::pressureIdx,
saturationIdx = Indices::saturationIdx,
transportEqIdx = 0,
saturationIdx = 0,
totalvelocityEqIdx = 1,
// phase indices
wPhaseIdx = FluidSystem::phase0Idx,
......@@ -158,13 +161,14 @@ public:
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes values;
// if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos) || onLowerBoundary_(globalPos)) {
// values.setAllNeumann();
// }
// else {
// values.setAllDirichlet();
// }
values.setAllNeumann();
if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos) || onLowerBoundary_(globalPos)) {
values.setAllNeumann();
}
else
{
values.setAllDirichlet();
}
// values.setAllNeumann();
return values;
}
......@@ -179,7 +183,8 @@ public:
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values(0.0);
values[saturationIdx] = 0.9;
values[saturationIdx] = 0.0;
values[pressureIdx] = 1.0e5;
return values;
}
......@@ -190,28 +195,35 @@ public:
const SubControlVolumeFace& scvf) const
{
NeumannFluxes values(0.0);
// const auto globalPos = scvf.ipGlobal();
//
// static const GlobalPosition v_t = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, GlobalPosition, Problem, TotalVelocity);
//
// if(onLowerBoundary_(globalPos))
// {
// //we assume constant saturation at the lower boundary
// const auto& volVars = elemVolvars[scvf.insideScvIdx()];
//
// Scalar mobW = volVars.mobility(wPhaseIdx);
// Scalar mobN = volVars.mobility(nPhaseIdx);
// Scalar mobT = mobW + mobN;
// auto K = this->spatialParams().permeabilityAtPos(globalPos);
//
// Scalar densityW = volVars.density(wPhaseIdx);
// Scalar densityN = volVars.density(nPhaseIdx);
//
// //ToDo use fluxVars for calculation
// Scalar fracW = mobW/mobT;
// Scalar velW = fracW*(v_t*scvf.unitOuterNormal()) + fracW*mobN*K*(densityN-densityW)*this->gravity()[dimWorld-1];
// values[transportEqIdx] = velW;
// }
const auto globalPos = scvf.ipGlobal();
static const GlobalPosition v_t = getParam<GlobalPosition>("Problem.TotalVelocity");
if(onLowerBoundary_(globalPos))
{
//we assume constant saturation at the lower boundary
const auto& volVars = elemVolvars[scvf.insideScvIdx()];
Scalar mobW = volVars.mobility(wPhaseIdx);
Scalar mobN = volVars.mobility(nPhaseIdx);
Scalar mobT = mobW + mobN;
auto K = this->spatialParams().permeabilityAtPos(globalPos);
Scalar densityW = volVars.density(wPhaseIdx);
Scalar densityN = volVars.density(nPhaseIdx);
//ToDo use fluxVars for calculation
Scalar fracW = mobW/mobT;
Scalar fracN = mobN/mobT;
Scalar fluxW = fracW*(v_t*scvf.unitOuterNormal()) + fracW*mobN*K*(densityN-densityW)*this->gravity()[dimWorld-1];
Scalar fluxN = fracN*(v_t*scvf.unitOuterNormal()) - fracN*mobW*K*(densityN-densityW)*this->gravity()[dimWorld-1];
fluxW *= scvf.area();
fluxN *= scvf.area();
values[transportEqIdx] = fluxW;
values[totalvelocityEqIdx] = fluxW + fluxN;
}
return values;
}
// \}
......@@ -232,6 +244,7 @@ public:
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values(0.0);
values[pressureIdx] = 1.0e5;
if(globalPos[1] > 0.65*this->fvGridGeometry().bBoxMax()[1] + eps_)
values[saturationIdx] = initSaturations_[0];
else
......@@ -244,22 +257,22 @@ private:
bool onLeftBoundary_(const GlobalPosition &globalPos) const
{
return globalPos[0] < this->bBoxMin()[0] + eps_;
return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_;
}
bool onRightBoundary_(const GlobalPosition &globalPos) const
{
return globalPos[0] > this->bBoxMax()[0] - eps_;
return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_;
}
bool onLowerBoundary_(const GlobalPosition &globalPos) const
{
return globalPos[1] < this->bBoxMin()[1] + eps_;
return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_;
}
bool onUpperBoundary_(const GlobalPosition &globalPos) const
{
return globalPos[1] > this->bBoxMax()[1] - eps_;
return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_;
}
Scalar temperature_;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment