26#ifndef DUMUX_RICHARDS_LOCAL_RESIDUAL_HH
27#define DUMUX_RICHARDS_LOCAL_RESIDUAL_HH
40template <
typename T,
typename ...Ts>
41using RobinDerivDetector =
decltype(std::declval<T>().addRobinFluxDerivatives(std::declval<Ts>()...));
43template<
class T,
typename ...Args>
56template<
class TypeTag>
70 using FVElementGeometry =
typename GridGeometry::LocalView;
72 using SubControlVolume =
typename GridGeometry::SubControlVolume;
73 using SubControlVolumeFace =
typename GridGeometry::SubControlVolumeFace;
74 using GridView =
typename GridGeometry::GridView;
75 using Element =
typename GridView::template Codim<0>::Entity;
80 enum { conti0EqIdx = Indices::conti0EqIdx };
84 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
85 gasPhaseIdx = FluidSystem::gasPhaseIdx,
86 liquidCompIdx = FluidSystem::liquidCompIdx
89 static constexpr bool enableWaterDiffusionInAir
90 = getPropValue<TypeTag, Properties::EnableWaterDiffusionInAir>();
96 double operator[] (
const Index i)
const
97 {
static_assert(
AlwaysFalse<Index>::value,
"Solution-dependent material parameters not supported with analytical differentiation");
return 0.0; }
101 using ParentType::ParentType;
115 const SubControlVolume& scv,
116 const VolumeVariables& volVars)
const
119 NumEqVector storage(0.0);
120 storage[conti0EqIdx] = volVars.porosity()
121 * volVars.density(liquidPhaseIdx)
122 * volVars.saturation(liquidPhaseIdx);
125 if constexpr (enableWaterDiffusionInAir)
126 storage[conti0EqIdx] += volVars.porosity()
127 * volVars.molarDensity(gasPhaseIdx)
128 * volVars.moleFraction(gasPhaseIdx, liquidCompIdx)
129 * FluidSystem::molarMass(liquidCompIdx)
130 * volVars.saturation(gasPhaseIdx);
133 EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, liquidPhaseIdx);
134 EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, gasPhaseIdx);
135 EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
152 const Element& element,
153 const FVElementGeometry& fvGeometry,
154 const ElementVolumeVariables& elemVolVars,
155 const SubControlVolumeFace& scvf,
156 const ElementFluxVariablesCache& elemFluxVarsCache)
const
158 FluxVariables fluxVars;
159 fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
161 NumEqVector flux(0.0);
163 auto upwindTerm = [](
const auto& volVars)
164 {
return volVars.density(liquidPhaseIdx)*volVars.mobility(liquidPhaseIdx); };
166 flux[conti0EqIdx] = fluxVars.advectiveFlux(liquidPhaseIdx, upwindTerm);
169 if constexpr (enableWaterDiffusionInAir)
173 flux[conti0EqIdx] += fluxVars.molecularDiffusionFlux(gasPhaseIdx)[liquidCompIdx];
175 flux[conti0EqIdx] += fluxVars.molecularDiffusionFlux(gasPhaseIdx)[liquidCompIdx]*FluidSystem::molarMass(liquidCompIdx);
179 EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, liquidPhaseIdx);
183 EnergyLocalResidual::heatConductionFlux(flux, fluxVars);
198 template<
class PartialDerivativeMatrix>
200 const Problem& problem,
201 const Element& element,
202 const FVElementGeometry& fvGeometry,
203 const VolumeVariables& curVolVars,
204 const SubControlVolume& scv)
const
206 static_assert(!enableWaterDiffusionInAir,
207 "richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
208 static_assert(!FluidSystem::isCompressible(0),
209 "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
211 const auto poreVolume =
Extrusion::volume(scv)*curVolVars.porosity()*curVolVars.extrusionFactor();
212 static const auto rho = curVolVars.density(0);
216 const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, InvalidElemSol{});
217 partialDerivatives[conti0EqIdx][0] += -rho*poreVolume/this->timeLoop().timeStepSize()*fluidMatrixInteraction.dsw_dpc(curVolVars.capillaryPressure());
232 template<
class PartialDerivativeMatrix>
234 const Problem& problem,
235 const Element& element,
236 const FVElementGeometry& fvGeometry,
237 const VolumeVariables& curVolVars,
238 const SubControlVolume& scv)
const
255 template<
class PartialDerivativeMatrices,
class T = TypeTag>
258 const Problem& problem,
259 const Element& element,
260 const FVElementGeometry& fvGeometry,
261 const ElementVolumeVariables& curElemVolVars,
262 const ElementFluxVariablesCache& elemFluxVarsCache,
263 const SubControlVolumeFace& scvf)
const
265 static_assert(!enableWaterDiffusionInAir,
266 "richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
267 static_assert(!FluidSystem::isCompressible(0),
268 "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
269 static_assert(!FluidSystem::viscosityIsConstant(0),
270 "richards/localresidual.hh: Analytic Jacobian only supports fluids with constant viscosity!");
273 const auto insideScvIdx = scvf.insideScvIdx();
274 const auto outsideScvIdx = scvf.outsideScvIdx();
275 const auto outsideElement = fvGeometry.gridGeometry().element(outsideScvIdx);
276 const auto& insideScv = fvGeometry.scv(insideScvIdx);
277 const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
278 const auto& insideVolVars = curElemVolVars[insideScvIdx];
279 const auto& outsideVolVars = curElemVolVars[outsideScvIdx];
282 static const auto rho = insideVolVars.density(0);
283 static const auto mu = insideVolVars.viscosity(0);
284 static const auto rho_mu = rho/mu;
289 static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(),
"Flux.UpwindWeight");
290 const auto flux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache);
291 const auto insideWeight = std::signbit(flux) ? (1.0 - upwindWeight) : upwindWeight;
292 const auto outsideWeight = 1.0 - insideWeight;
293 const auto upwindTerm = rho*insideVolVars.mobility(0)*insideWeight + rho*outsideVolVars.mobility(0)*outsideWeight;
295 const auto insideFluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, insideScv, InvalidElemSol{});
296 const auto outsideFluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(outsideElement, outsideScv, InvalidElemSol{});
299 const auto insideSw = insideVolVars.saturation(0);
300 const auto outsideSw = outsideVolVars.saturation(0);
301 const auto insidePc = insideVolVars.capillaryPressure();
302 const auto outsidePc = outsideVolVars.capillaryPressure();
303 const auto dkrw_dsw_inside = insideFluidMatrixInteraction.dkrw_dsw(insideSw);
304 const auto dkrw_dsw_outside = outsideFluidMatrixInteraction.dkrw_dsw(outsideSw);
305 const auto dsw_dpw_inside = -insideFluidMatrixInteraction.dsw_dpc(insidePc);
306 const auto dsw_dpw_outside = -outsideFluidMatrixInteraction.dsw_dpc(outsidePc);
309 const auto tij = elemFluxVarsCache[scvf].advectionTij();
312 auto& dI_dI = derivativeMatrices[insideScvIdx];
313 auto& dI_dJ = derivativeMatrices[outsideScvIdx];
316 dI_dI[conti0EqIdx][0] += tij*upwindTerm + rho_mu*flux*insideWeight*dkrw_dsw_inside*dsw_dpw_inside;
317 dI_dJ[conti0EqIdx][0] += -tij*upwindTerm + rho_mu*flux*outsideWeight*dkrw_dsw_outside*dsw_dpw_outside;
331 template<
class JacobianMatrix,
class T = TypeTag>
334 const Problem& problem,
335 const Element& element,
336 const FVElementGeometry& fvGeometry,
337 const ElementVolumeVariables& curElemVolVars,
338 const ElementFluxVariablesCache& elemFluxVarsCache,
339 const SubControlVolumeFace& scvf)
const
341 static_assert(!enableWaterDiffusionInAir,
342 "richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
343 static_assert(!FluidSystem::isCompressible(0),
344 "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
345 static_assert(!FluidSystem::viscosityIsConstant(0),
346 "richards/localresidual.hh: Analytic Jacobian only supports fluids with constant viscosity!");
349 const auto insideScvIdx = scvf.insideScvIdx();
350 const auto outsideScvIdx = scvf.outsideScvIdx();
351 const auto& insideScv = fvGeometry.scv(insideScvIdx);
352 const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
353 const auto& insideVolVars = curElemVolVars[insideScvIdx];
354 const auto& outsideVolVars = curElemVolVars[outsideScvIdx];
357 static const auto rho = insideVolVars.density(0);
358 static const auto mu = insideVolVars.viscosity(0);
359 static const auto rho_mu = rho/mu;
364 static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(),
"Flux.UpwindWeight");
365 const auto flux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache);
366 const auto insideWeight = std::signbit(flux) ? (1.0 - upwindWeight) : upwindWeight;
367 const auto outsideWeight = 1.0 - insideWeight;
368 const auto upwindTerm = rho*insideVolVars.mobility(0)*insideWeight + rho*outsideVolVars.mobility(0)*outsideWeight;
370 const auto insideFluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, insideScv, InvalidElemSol{});
371 const auto outsideFluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, outsideScv, InvalidElemSol{});
374 const auto insideSw = insideVolVars.saturation(0);
375 const auto outsideSw = outsideVolVars.saturation(0);
376 const auto insidePc = insideVolVars.capillaryPressure();
377 const auto outsidePc = outsideVolVars.capillaryPressure();
378 const auto dkrw_dsw_inside = insideFluidMatrixInteraction.dkrw_dsw(insideSw);
379 const auto dkrw_dsw_outside = outsideFluidMatrixInteraction.dkrw_dsw(outsideSw);
380 const auto dsw_dpw_inside = -insideFluidMatrixInteraction.dsw_dpc(insidePc);
381 const auto dsw_dpw_outside = -outsideFluidMatrixInteraction.dsw_dpc(outsidePc);
385 const auto ti = AdvectionType::calculateTransmissibilities(
386 problem, element, fvGeometry, curElemVolVars, scvf, elemFluxVarsCache[scvf]
390 auto& dI_dJ_inside = A[insideScv.dofIndex()];
391 auto& dI_dJ_outside = A[outsideScv.dofIndex()];
394 for (
const auto& scvJ : scvs(fvGeometry))
397 const auto& tj = ti[scvJ.indexInElement()];
398 const auto globalJ = scvJ.dofIndex();
401 dI_dJ_inside[globalJ][conti0EqIdx][0] += tj*upwindTerm;
402 dI_dJ_outside[globalJ][conti0EqIdx][0] += -tj*upwindTerm;
405 const auto upwindContributionInside = rho_mu*flux*insideWeight*dkrw_dsw_inside*dsw_dpw_inside;
406 const auto upwindContributionOutside = rho_mu*flux*outsideWeight*dkrw_dsw_outside*dsw_dpw_outside;
409 A[insideScv.dofIndex()][insideScv.dofIndex()][conti0EqIdx][0] += upwindContributionInside;
410 A[insideScv.dofIndex()][outsideScv.dofIndex()][conti0EqIdx][0] += upwindContributionOutside;
412 A[outsideScv.dofIndex()][insideScv.dofIndex()][conti0EqIdx][0] -= upwindContributionInside;
413 A[outsideScv.dofIndex()][outsideScv.dofIndex()][conti0EqIdx][0] -= upwindContributionOutside;
430 template<
class PartialDerivativeMatrices>
432 const Problem& problem,
433 const Element& element,
434 const FVElementGeometry& fvGeometry,
435 const ElementVolumeVariables& curElemVolVars,
436 const ElementFluxVariablesCache& elemFluxVarsCache,
437 const SubControlVolumeFace& scvf)
const
439 static_assert(!enableWaterDiffusionInAir,
440 "richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
441 static_assert(!FluidSystem::isCompressible(0),
442 "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
443 static_assert(!FluidSystem::viscosityIsConstant(0),
444 "richards/localresidual.hh: Analytic Jacobian only supports fluids with constant viscosity!");
448 const auto insideScvIdx = scvf.insideScvIdx();
449 const auto& insideScv = fvGeometry.scv(insideScvIdx);
450 const auto& insideVolVars = curElemVolVars[insideScvIdx];
451 const auto& outsideVolVars = curElemVolVars[scvf.outsideScvIdx()];
452 const auto insideFluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, insideScv, InvalidElemSol{});
455 static const auto rho = insideVolVars.density(0);
456 static const auto mu = insideVolVars.viscosity(0);
457 static const auto rho_mu = rho/mu;
462 static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(),
"Flux.UpwindWeight");
463 const auto flux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache);
464 const auto insideWeight = std::signbit(flux) ? (1.0 - upwindWeight) : upwindWeight;
465 const auto outsideWeight = 1.0 - insideWeight;
466 const auto upwindTerm = rho*insideVolVars.mobility(0)*insideWeight + rho*outsideVolVars.mobility(0)*outsideWeight;
469 const auto insideSw = insideVolVars.saturation(0);
470 const auto insidePc = insideVolVars.capillaryPressure();
471 const auto dkrw_dsw_inside = insideFluidMatrixInteraction.dkrw_dsw(insideSw);
472 const auto dsw_dpw_inside = -insideFluidMatrixInteraction.dsw_dpc(insidePc);
475 const auto tij = elemFluxVarsCache[scvf].advectionTij();
478 derivativeMatrices[insideScvIdx][conti0EqIdx][0] += tij*upwindTerm + rho_mu*flux*insideWeight*dkrw_dsw_inside*dsw_dpw_inside;
492 template<
class PartialDerivativeMatrices>
494 const Problem& problem,
495 const Element& element,
496 const FVElementGeometry& fvGeometry,
497 const ElementVolumeVariables& curElemVolVars,
498 const ElementFluxVariablesCache& elemFluxVarsCache,
499 const SubControlVolumeFace& scvf)
const
502 PartialDerivativeMatrices&, Element, FVElementGeometry,
503 ElementVolumeVariables, ElementFluxVariablesCache, SubControlVolumeFace>()
505 problem.addRobinFluxDerivatives(derivativeMatrices, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);
509 Implementation *asImp_()
510 {
return static_cast<Implementation *
> (
this); }
512 const Implementation *asImp_()
const
513 {
return static_cast<const Implementation *
> (
this); }
A helper to deduce a vector with the same size as numbers of equations.
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 NumEqVectorTraits< PrimaryVariables >::type NumEqVector
A vector with the same size as numbers of equations This is the default implementation and has to be ...
Definition: numeqvector.hh:46
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
Definition: propertysystem.hh:150
Definition: gridcapabilities.hh:60
decltype(std::declval< T >().addRobinFluxDerivatives(std::declval< Ts >()...)) RobinDerivDetector
Definition: porousmediumflow/richards/localresidual.hh:41
static constexpr bool hasAddRobinFluxDerivatives()
Definition: porousmediumflow/richards/localresidual.hh:44
Scalar volume(Shape shape, Scalar inscribedRadius)
Returns the volume of a given geometry based on the inscribed radius.
Definition: poreproperties.hh:73
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:58
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:233
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:431
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:114
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:199
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:151
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:257
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:493
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:333
Declares all properties used in Dumux.