Commit 88154e73 authored by Martin Schneider's avatar Martin Schneider
Browse files

[IMPES] Implementation of analytical derivatives for pressure eq

parent f5973da9
......@@ -59,6 +59,8 @@ class PressureLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
static constexpr int numPhases = ModelTraits::numPhases();
static constexpr int pressureEqIdx = ModelTraits::Indices::pressureEqIdx; //!< first index for the mass balance
static constexpr int pressureIdx = ModelTraits::Indices::pressureIdx; //!< first index for the mass balance
public:
using ParentType::ParentType;
......@@ -119,26 +121,150 @@ public:
return flux;
}
/*!
* \brief TODO docme!
*
* \param partialDerivatives TODO docme!
* \param problem TODO docme!
* \param element The element
* \param fvGeometry The finite volume geometry context
* \param curVolVars TODO docme!
* \param scv The sub control volume.
*/
template<class PartialDerivativeMatrix>
void addStorageDerivatives(PartialDerivativeMatrix& partialDerivatives,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const VolumeVariables& curVolVars,
const SubControlVolume& scv) const
const SubControlVolume& scv) const {}
template<class PartialDerivativeMatrix>
void addSourceDerivatives(PartialDerivativeMatrix& partialDerivatives,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const VolumeVariables& curVolVars,
const SubControlVolume& scv) const
{
// Todo what about a slightly compressible case?
//problem.addSourceDerivatives(partialDerivatives, element, fvGeometry, curVolVars, scv);
}
//! flux derivatives for the cell-centered tpfa scheme
template<class PartialDerivativeMatrices, class T = TypeTag>
std::enable_if_t<GET_PROP_TYPE(T, FVGridGeometry)::discMethod == DiscretizationMethod::cctpfa, void>
addFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& curElemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
static_assert(!FluidSystem::isCompressible(0),
"2p/sequential/pressure/localresidual.hh: Only incompressible fluids are allowed!");
static_assert(!FluidSystem::isCompressible(1),
"2p/sequential/pressure/localresidual.hh: Only incompressible fluids are allowed!");
static_assert(FluidSystem::viscosityIsConstant(0),
"2p/sequential/pressure/localresidual.hh: Only fluids with constant viscosities are allowed!");
static_assert(FluidSystem::viscosityIsConstant(1),
"2p/sequential/pressure/localresidual.hh: Only fluids with constant viscosities are allowed!");
static_assert(ModelTraits::numPhases() == 2,
"2p/sequential/pressure/localresidual.hh: Only two-phase models are allowed!");
//static_assert(ModelTraits::priVarFormulation() == TwoPFormulation::p0s1,
// "2p/sequential/pressure/incompressiblelocalresidual.hh: Analytic differentiation has to be checked for p1-s0 formulation!");
using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
// evaluate the current wetting phase Darcy flux and resulting upwind weights
const auto flux_w = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache);
const auto flux_n = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 1, elemFluxVarsCache);
// get references to the two participating vol vars & parameters
const auto insideScvIdx = scvf.insideScvIdx();
const auto outsideScvIdx = scvf.outsideScvIdx();
const auto outsideElement = fvGeometry.fvGridGeometry().element(outsideScvIdx);
const auto& insideScv = fvGeometry.scv(insideScvIdx);
const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
const auto& insideVolVars = curElemVolVars[insideScvIdx];
const auto& outsideVolVars = curElemVolVars[outsideScvIdx];
const auto& insideMaterialParams = problem.spatialParams().materialLawParams(element,
insideScv,
elementSolution<FVElementGeometry>(insideVolVars.priVars()));
const auto& outsideMaterialParams = problem.spatialParams().materialLawParams(outsideElement,
outsideScv,
elementSolution<FVElementGeometry>(outsideVolVars.priVars()));
// some quantities to be reused (rho & mu are constant and thus equal for all cells)
const auto mobW_inside = insideVolVars.mobility(0);
const auto mobN_inside = insideVolVars.mobility(1);
const auto mobW_outside = outsideVolVars.mobility(0);
const auto mobN_outside = outsideVolVars.mobility(1);
const auto tij = elemFluxVarsCache[scvf].advectionTij();
// precalculate values
const auto mobup_w = (flux_w > 0.) ? mobW_inside : mobW_outside;
const auto mobup_n = (flux_n > 0.) ? mobN_inside : mobN_outside;
const auto tij_up_t = tij*(mobup_w + mobup_n);
// partial derivative of the wetting phase flux w.r.t. p_w
// add partial derivatives to the respective given matrices
derivativeMatrices[scvf.insideScvIdx()][pressureEqIdx][pressureIdx] += tij_up_t;
derivativeMatrices[scvf.outsideScvIdx()][pressureEqIdx][pressureIdx] -= tij_up_t;
}
//! Dirichlet flux derivatives for the cell-centered tpfa scheme
template<class PartialDerivativeMatrices, class T = TypeTag>
std::enable_if_t<GET_PROP_TYPE(T, FVGridGeometry)::discMethod == DiscretizationMethod::cctpfa, void>
addCCDirichletFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& curElemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
// evaluate the current wetting phase Darcy flux and resulting upwind weights
const auto flux_w = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache);
const auto flux_n = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 1, elemFluxVarsCache);
// get references to the two participating vol vars & parameters
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideScv = fvGeometry.scv(insideScvIdx);
const auto& insideVolVars = curElemVolVars[insideScvIdx];
const auto& outsideVolVars = curElemVolVars[scvf.outsideScvIdx()];
const auto& insideMaterialParams = problem.spatialParams().materialLawParams(element,
insideScv,
elementSolution<FVElementGeometry>(insideVolVars.priVars()));
// some quantities to be reused (rho & mu are constant and thus equal for all cells)
const auto mobW_inside = insideVolVars.mobility(0);
const auto mobN_inside = insideVolVars.mobility(1);
const auto mobW_outside = outsideVolVars.mobility(0);
const auto mobN_outside = outsideVolVars.mobility(1);
const auto tij = elemFluxVarsCache[scvf].advectionTij();
// precalculate values
const auto mobup_w = (flux_w > 0.) ? mobW_inside : mobW_outside;
const auto mobup_n = (flux_n > 0.) ? mobN_inside : mobN_outside;
const auto tij_up_t = tij*(mobup_w + mobup_n);
// partial derivative of the wetting phase flux w.r.t. p_w
// add partial derivatives to the respective given matrices
derivativeMatrices[scvf.insideScvIdx()][pressureEqIdx][pressureIdx] += tij_up_t;
}
//! Robin-type flux derivatives
template<class PartialDerivativeMatrices>
void addRobinFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& curElemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
//! Robin-type boundary conditions are problem-specific.
//! We can't put a general implementation here - users defining Robin-type BCs
//! while using analytical Jacobian assembly must overload this function!
}
};
} // end namespace Dumux
......
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