26#ifndef DUMUX_RICHARDS_LOCAL_RESIDUAL_HH
27#define DUMUX_RICHARDS_LOCAL_RESIDUAL_HH
45template<
class TypeTag>
59 using FVElementGeometry =
typename GridGeometry::LocalView;
61 using SubControlVolume =
typename GridGeometry::SubControlVolume;
62 using SubControlVolumeFace =
typename GridGeometry::SubControlVolumeFace;
63 using GridView =
typename GridGeometry::GridView;
64 using Element =
typename GridView::template Codim<0>::Entity;
69 enum { conti0EqIdx = Indices::conti0EqIdx };
73 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
74 gasPhaseIdx = FluidSystem::gasPhaseIdx,
75 liquidCompIdx = FluidSystem::liquidCompIdx
78 static constexpr bool enableWaterDiffusionInAir
79 = getPropValue<TypeTag, Properties::EnableWaterDiffusionInAir>();
85 double operator[] (
const Index i)
const
86 {
static_assert(
AlwaysFalse<Index>::value,
"Solution-dependent material parameters not supported with analytical differentiation");
return 0.0; }
90 using ParentType::ParentType;
104 const SubControlVolume& scv,
105 const VolumeVariables& volVars)
const
108 NumEqVector storage(0.0);
109 storage[conti0EqIdx] = volVars.porosity()
110 * volVars.density(liquidPhaseIdx)
111 * volVars.saturation(liquidPhaseIdx);
114 if (enableWaterDiffusionInAir)
115 storage[conti0EqIdx] += volVars.porosity()
116 * volVars.molarDensity(gasPhaseIdx)
117 * volVars.moleFraction(gasPhaseIdx, liquidCompIdx)
118 * FluidSystem::molarMass(liquidCompIdx)
119 * volVars.saturation(gasPhaseIdx);
122 EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, liquidPhaseIdx);
123 EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, gasPhaseIdx);
124 EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
141 const Element& element,
142 const FVElementGeometry& fvGeometry,
143 const ElementVolumeVariables& elemVolVars,
144 const SubControlVolumeFace& scvf,
145 const ElementFluxVariablesCache& elemFluxVarsCache)
const
147 FluxVariables fluxVars;
148 fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
150 NumEqVector flux(0.0);
152 auto upwindTerm = [](
const auto& volVars)
153 {
return volVars.density(liquidPhaseIdx)*volVars.mobility(liquidPhaseIdx); };
155 flux[conti0EqIdx] = fluxVars.advectiveFlux(liquidPhaseIdx, upwindTerm);
158 if (enableWaterDiffusionInAir)
162 flux[conti0EqIdx] += fluxVars.molecularDiffusionFlux(gasPhaseIdx)[liquidCompIdx];
164 flux[conti0EqIdx] += fluxVars.molecularDiffusionFlux(gasPhaseIdx)[liquidCompIdx]*FluidSystem::molarMass(liquidCompIdx);
168 EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, liquidPhaseIdx);
172 EnergyLocalResidual::heatConductionFlux(flux, fluxVars);
187 template<
class PartialDerivativeMatrix>
189 const Problem& problem,
190 const Element& element,
191 const FVElementGeometry& fvGeometry,
192 const VolumeVariables& curVolVars,
193 const SubControlVolume& scv)
const
195 static_assert(!enableWaterDiffusionInAir,
196 "richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
197 static_assert(!FluidSystem::isCompressible(0),
198 "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
200 const auto poreVolume = Extrusion::volume(scv)*curVolVars.porosity();
201 static const auto rho = curVolVars.density(0);
206 const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem.spatialParams(), element, scv, InvalidElemSol{});
210 partialDerivatives[conti0EqIdx][0] += -rho*poreVolume/this->timeLoop().timeStepSize()*fluidMatrixInteraction.dsw_dpc(curVolVars.capillaryPressure());
225 template<
class PartialDerivativeMatrix>
227 const Problem& problem,
228 const Element& element,
229 const FVElementGeometry& fvGeometry,
230 const VolumeVariables& curVolVars,
231 const SubControlVolume& scv)
const
248 template<
class PartialDerivativeMatrices,
class T = TypeTag>
251 const Problem& problem,
252 const Element& element,
253 const FVElementGeometry& fvGeometry,
254 const ElementVolumeVariables& curElemVolVars,
255 const ElementFluxVariablesCache& elemFluxVarsCache,
256 const SubControlVolumeFace& scvf)
const
258 static_assert(!enableWaterDiffusionInAir,
259 "richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
260 static_assert(!FluidSystem::isCompressible(0),
261 "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
262 static_assert(!FluidSystem::viscosityIsConstant(0),
263 "richards/localresidual.hh: Analytic Jacobian only supports fluids with constant viscosity!");
266 const auto insideScvIdx = scvf.insideScvIdx();
267 const auto outsideScvIdx = scvf.outsideScvIdx();
268 const auto outsideElement = fvGeometry.gridGeometry().element(outsideScvIdx);
269 const auto& insideScv = fvGeometry.scv(insideScvIdx);
270 const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
271 const auto& insideVolVars = curElemVolVars[insideScvIdx];
272 const auto& outsideVolVars = curElemVolVars[outsideScvIdx];
275 static const auto rho = insideVolVars.density(0);
276 static const auto mu = insideVolVars.viscosity(0);
277 static const auto rho_mu = rho/mu;
282 static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(),
"Flux.UpwindWeight");
283 const auto flux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache);
284 const auto insideWeight = std::signbit(flux) ? (1.0 - upwindWeight) : upwindWeight;
285 const auto outsideWeight = 1.0 - insideWeight;
286 const auto upwindTerm = rho*insideVolVars.mobility(0)*insideWeight + rho*outsideVolVars.mobility(0)*outsideWeight;
291 const auto insideFluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem.spatialParams(), element, insideScv, InvalidElemSol{});
292 const auto outsideFluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem.spatialParams(), outsideElement, outsideScv, InvalidElemSol{});
295 const auto insideSw = insideVolVars.saturation(0);
296 const auto outsideSw = outsideVolVars.saturation(0);
297 const auto insidePc = insideVolVars.capillaryPressure();
298 const auto outsidePc = outsideVolVars.capillaryPressure();
299 const auto dkrw_dsw_inside = insideFluidMatrixInteraction.dkrw_dsw(insideSw);
300 const auto dkrw_dsw_outside = outsideFluidMatrixInteraction.dkrw_dsw(outsideSw);
301 const auto dsw_dpw_inside = -insideFluidMatrixInteraction.dsw_dpc(insidePc);
302 const auto dsw_dpw_outside = -outsideFluidMatrixInteraction.dsw_dpc(outsidePc);
305 const auto tij = elemFluxVarsCache[scvf].advectionTij();
308 auto& dI_dI = derivativeMatrices[insideScvIdx];
309 auto& dI_dJ = derivativeMatrices[outsideScvIdx];
312 dI_dI[conti0EqIdx][0] += tij*upwindTerm + rho_mu*flux*insideWeight*dkrw_dsw_inside*dsw_dpw_inside;
313 dI_dJ[conti0EqIdx][0] += -tij*upwindTerm + rho_mu*flux*outsideWeight*dkrw_dsw_outside*dsw_dpw_outside;
327 template<
class JacobianMatrix,
class T = TypeTag>
330 const Problem& problem,
331 const Element& element,
332 const FVElementGeometry& fvGeometry,
333 const ElementVolumeVariables& curElemVolVars,
334 const ElementFluxVariablesCache& elemFluxVarsCache,
335 const SubControlVolumeFace& scvf)
const
337 static_assert(!enableWaterDiffusionInAir,
338 "richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
339 static_assert(!FluidSystem::isCompressible(0),
340 "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
341 static_assert(!FluidSystem::viscosityIsConstant(0),
342 "richards/localresidual.hh: Analytic Jacobian only supports fluids with constant viscosity!");
345 const auto insideScvIdx = scvf.insideScvIdx();
346 const auto outsideScvIdx = scvf.outsideScvIdx();
347 const auto& insideScv = fvGeometry.scv(insideScvIdx);
348 const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
349 const auto& insideVolVars = curElemVolVars[insideScvIdx];
350 const auto& outsideVolVars = curElemVolVars[outsideScvIdx];
353 static const auto rho = insideVolVars.density(0);
354 static const auto mu = insideVolVars.viscosity(0);
355 static const auto rho_mu = rho/mu;
360 static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(),
"Flux.UpwindWeight");
361 const auto flux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache);
362 const auto insideWeight = std::signbit(flux) ? (1.0 - upwindWeight) : upwindWeight;
363 const auto outsideWeight = 1.0 - insideWeight;
364 const auto upwindTerm = rho*insideVolVars.mobility(0)*insideWeight + rho*outsideVolVars.mobility(0)*outsideWeight;
369 const auto insideFluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem.spatialParams(), element, insideScv, InvalidElemSol{});
370 const auto outsideFluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem.spatialParams(), element, outsideScv, InvalidElemSol{});
373 const auto insideSw = insideVolVars.saturation(0);
374 const auto outsideSw = outsideVolVars.saturation(0);
375 const auto insidePc = insideVolVars.capillaryPressure();
376 const auto outsidePc = outsideVolVars.capillaryPressure();
377 const auto dkrw_dsw_inside = insideFluidMatrixInteraction.dkrw_dsw(insideSw);
378 const auto dkrw_dsw_outside = outsideFluidMatrixInteraction.dkrw_dsw(outsideSw);
379 const auto dsw_dpw_inside = -insideFluidMatrixInteraction.dsw_dpc(insidePc);
380 const auto dsw_dpw_outside = -outsideFluidMatrixInteraction.dsw_dpc(outsidePc);
384 const auto ti = AdvectionType::calculateTransmissibilities(
385 problem, element, fvGeometry, curElemVolVars, scvf, elemFluxVarsCache[scvf]
389 auto& dI_dJ_inside = A[insideScv.dofIndex()];
390 auto& dI_dJ_outside = A[outsideScv.dofIndex()];
393 for (
const auto& scvJ : scvs(fvGeometry))
396 const auto& tj = ti[scvJ.indexInElement()];
397 const auto globalJ = scvJ.dofIndex();
400 dI_dJ_inside[globalJ][conti0EqIdx][0] += tj*upwindTerm;
401 dI_dJ_outside[globalJ][conti0EqIdx][0] += -tj*upwindTerm;
404 const auto upwindContributionInside = rho_mu*flux*insideWeight*dkrw_dsw_inside*dsw_dpw_inside;
405 const auto upwindContributionOutside = rho_mu*flux*outsideWeight*dkrw_dsw_outside*dsw_dpw_outside;
408 A[insideScv.dofIndex()][insideScv.dofIndex()][conti0EqIdx][0] += upwindContributionInside;
409 A[insideScv.dofIndex()][outsideScv.dofIndex()][conti0EqIdx][0] += upwindContributionOutside;
411 A[outsideScv.dofIndex()][insideScv.dofIndex()][conti0EqIdx][0] -= upwindContributionInside;
412 A[outsideScv.dofIndex()][outsideScv.dofIndex()][conti0EqIdx][0] -= upwindContributionOutside;
429 template<
class PartialDerivativeMatrices>
431 const Problem& problem,
432 const Element& element,
433 const FVElementGeometry& fvGeometry,
434 const ElementVolumeVariables& curElemVolVars,
435 const ElementFluxVariablesCache& elemFluxVarsCache,
436 const SubControlVolumeFace& scvf)
const
438 static_assert(!enableWaterDiffusionInAir,
439 "richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
440 static_assert(!FluidSystem::isCompressible(0),
441 "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
442 static_assert(!FluidSystem::viscosityIsConstant(0),
443 "richards/localresidual.hh: Analytic Jacobian only supports fluids with constant viscosity!");
447 const auto insideScvIdx = scvf.insideScvIdx();
448 const auto& insideScv = fvGeometry.scv(insideScvIdx);
449 const auto& insideVolVars = curElemVolVars[insideScvIdx];
450 const auto& outsideVolVars = curElemVolVars[scvf.outsideScvIdx()];
455 const auto insideFluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem.spatialParams(), element, insideScv, InvalidElemSol{});
458 static const auto rho = insideVolVars.density(0);
459 static const auto mu = insideVolVars.viscosity(0);
460 static const auto rho_mu = rho/mu;
465 static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(),
"Flux.UpwindWeight");
466 const auto flux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache);
467 const auto insideWeight = std::signbit(flux) ? (1.0 - upwindWeight) : upwindWeight;
468 const auto outsideWeight = 1.0 - insideWeight;
469 const auto upwindTerm = rho*insideVolVars.mobility(0)*insideWeight + rho*outsideVolVars.mobility(0)*outsideWeight;
472 const auto insideSw = insideVolVars.saturation(0);
473 const auto insidePc = insideVolVars.capillaryPressure();
474 const auto dkrw_dsw_inside = insideFluidMatrixInteraction.dkrw_dsw(insideSw);
475 const auto dsw_dpw_inside = -insideFluidMatrixInteraction.dsw_dpc(insidePc);
478 const auto tij = elemFluxVarsCache[scvf].advectionTij();
481 derivativeMatrices[insideScvIdx][conti0EqIdx][0] += tij*upwindTerm + rho_mu*flux*insideWeight*dkrw_dsw_inside*dsw_dpw_inside;
495 template<
class PartialDerivativeMatrices>
497 const Problem& problem,
498 const Element& element,
499 const FVElementGeometry& fvGeometry,
500 const ElementVolumeVariables& curElemVolVars,
501 const ElementFluxVariablesCache& elemFluxVarsCache,
502 const SubControlVolumeFace& scvf)
const
506 Implementation *asImp_()
507 {
return static_cast<Implementation *
> (
this); }
509 const Implementation *asImp_()
const
510 {
return static_cast<const Implementation *
> (
this); }
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
The available discretization methods in Dumux.
Helper classes to compute the integration elements.
The reference frameworks and formulations available for splitting total fluxes into a advective and d...
typename Extrusion< T >::type Extrusion_t
Convenience alias for obtaining the extrusion type.
Definition: extrusion.hh:177
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
Template which always yields a false value.
Definition: typetraits.hh:36
Element-wise calculation of the Jacobian matrix for problems using the Richards fully implicit models...
Definition: porousmediumflow/richards/localresidual.hh:47
void addSourceDerivatives(PartialDerivativeMatrix &partialDerivatives, const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const VolumeVariables &curVolVars, const SubControlVolume &scv) const
Adds source derivatives for wetting and nonwetting phase.
Definition: porousmediumflow/richards/localresidual.hh:226
void addCCDirichletFluxDerivatives(PartialDerivativeMatrices &derivativeMatrices, const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &curElemVolVars, const ElementFluxVariablesCache &elemFluxVarsCache, const SubControlVolumeFace &scvf) const
Adds cell-centered Dirichlet flux derivatives for wetting and nonwetting phase.
Definition: porousmediumflow/richards/localresidual.hh:430
NumEqVector computeStorage(const Problem &problem, const SubControlVolume &scv, const VolumeVariables &volVars) const
Evaluates the rate of change of all conservation quantites (e.g. phase mass) within a sub-control vol...
Definition: porousmediumflow/richards/localresidual.hh:103
void addStorageDerivatives(PartialDerivativeMatrix &partialDerivatives, const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const VolumeVariables &curVolVars, const SubControlVolume &scv) const
Adds the storage derivative.
Definition: porousmediumflow/richards/localresidual.hh:188
NumEqVector computeFlux(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolumeFace &scvf, const ElementFluxVariablesCache &elemFluxVarsCache) const
Evaluates the mass flux over a face of a sub control volume.
Definition: porousmediumflow/richards/localresidual.hh:140
std::enable_if_t< GetPropType< T, Properties::GridGeometry >::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
Adds flux derivatives for wetting and nonwetting phase for cell-centered FVM using TPFA.
Definition: porousmediumflow/richards/localresidual.hh:250
void addRobinFluxDerivatives(PartialDerivativeMatrices &derivativeMatrices, const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &curElemVolVars, const ElementFluxVariablesCache &elemFluxVarsCache, const SubControlVolumeFace &scvf) const
Adds Robin flux derivatives for wetting and nonwetting phase.
Definition: porousmediumflow/richards/localresidual.hh:496
std::enable_if_t< GetPropType< T, Properties::GridGeometry >::discMethod==DiscretizationMethod::box, void > addFluxDerivatives(JacobianMatrix &A, const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &curElemVolVars, const ElementFluxVariablesCache &elemFluxVarsCache, const SubControlVolumeFace &scvf) const
Adds flux derivatives for box method.
Definition: porousmediumflow/richards/localresidual.hh:329
Declares all properties used in Dumux.