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 };
86 static constexpr auto liquidPhaseIdx = FluidSystem::phase0Idx;
87 static constexpr auto gasPhaseIdx = FluidSystem::phase1Idx;
88 static constexpr auto liquidCompIdx = FluidSystem::comp0Idx;
90 static constexpr bool enableWaterDiffusionInAir
91 = getPropValue<TypeTag, Properties::EnableWaterDiffusionInAir>();
97 double operator[] (
const Index i)
const
98 {
static_assert(
AlwaysFalse<Index>::value,
"Solution-dependent material parameters not supported with analytical differentiation");
return 0.0; }
102 using ParentType::ParentType;
116 const SubControlVolume& scv,
117 const VolumeVariables& volVars)
const
120 NumEqVector storage(0.0);
121 storage[conti0EqIdx] = volVars.porosity()
122 * volVars.density(liquidPhaseIdx)
123 * volVars.saturation(liquidPhaseIdx);
126 if constexpr (enableWaterDiffusionInAir)
127 storage[conti0EqIdx] += volVars.porosity()
128 * volVars.molarDensity(gasPhaseIdx)
129 * volVars.moleFraction(gasPhaseIdx, liquidCompIdx)
130 * FluidSystem::molarMass(liquidCompIdx)
131 * volVars.saturation(gasPhaseIdx);
134 EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, liquidPhaseIdx);
135 EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, gasPhaseIdx);
136 EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
153 const Element& element,
154 const FVElementGeometry& fvGeometry,
155 const ElementVolumeVariables& elemVolVars,
156 const SubControlVolumeFace& scvf,
157 const ElementFluxVariablesCache& elemFluxVarsCache)
const
159 FluxVariables fluxVars;
160 fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
162 NumEqVector flux(0.0);
164 auto upwindTerm = [](
const auto& volVars)
165 {
return volVars.density(liquidPhaseIdx)*volVars.mobility(liquidPhaseIdx); };
167 flux[conti0EqIdx] = fluxVars.advectiveFlux(liquidPhaseIdx, upwindTerm);
170 if constexpr (enableWaterDiffusionInAir)
174 flux[conti0EqIdx] += fluxVars.molecularDiffusionFlux(gasPhaseIdx)[liquidCompIdx];
176 flux[conti0EqIdx] += fluxVars.molecularDiffusionFlux(gasPhaseIdx)[liquidCompIdx]*FluidSystem::molarMass(liquidCompIdx);
180 EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, liquidPhaseIdx);
184 EnergyLocalResidual::heatConductionFlux(flux, fluxVars);
199 template<
class PartialDerivativeMatrix>
201 const Problem& problem,
202 const Element& element,
203 const FVElementGeometry& fvGeometry,
204 const VolumeVariables& curVolVars,
205 const SubControlVolume& scv)
const
207 static_assert(!enableWaterDiffusionInAir,
208 "richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
209 static_assert(!FluidSystem::isCompressible(0),
210 "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
212 const auto poreVolume =
Extrusion::volume(scv)*curVolVars.porosity()*curVolVars.extrusionFactor();
213 static const auto rho = curVolVars.density(0);
217 const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, InvalidElemSol{});
218 partialDerivatives[conti0EqIdx][0] += -rho*poreVolume/this->timeLoop().timeStepSize()*fluidMatrixInteraction.dsw_dpc(curVolVars.capillaryPressure());
233 template<
class PartialDerivativeMatrix>
235 const Problem& problem,
236 const Element& element,
237 const FVElementGeometry& fvGeometry,
238 const VolumeVariables& curVolVars,
239 const SubControlVolume& scv)
const
256 template<
class PartialDerivativeMatrices,
class T = TypeTag>
259 const Problem& problem,
260 const Element& element,
261 const FVElementGeometry& fvGeometry,
262 const ElementVolumeVariables& curElemVolVars,
263 const ElementFluxVariablesCache& elemFluxVarsCache,
264 const SubControlVolumeFace& scvf)
const
266 static_assert(!enableWaterDiffusionInAir,
267 "richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
268 static_assert(!FluidSystem::isCompressible(0),
269 "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
270 static_assert(FluidSystem::viscosityIsConstant(0),
271 "richards/localresidual.hh: Analytic Jacobian only supports fluids with constant viscosity!");
274 const auto insideScvIdx = scvf.insideScvIdx();
275 const auto outsideScvIdx = scvf.outsideScvIdx();
276 const auto outsideElement = fvGeometry.gridGeometry().element(outsideScvIdx);
277 const auto& insideScv = fvGeometry.scv(insideScvIdx);
278 const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
279 const auto& insideVolVars = curElemVolVars[insideScvIdx];
280 const auto& outsideVolVars = curElemVolVars[outsideScvIdx];
283 static const auto rho = insideVolVars.density(0);
284 static const auto mu = insideVolVars.viscosity(0);
285 static const auto rho_mu = rho/mu;
290 static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(),
"Flux.UpwindWeight");
291 const auto flux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache);
292 const auto insideWeight = std::signbit(flux) ? (1.0 - upwindWeight) : upwindWeight;
293 const auto outsideWeight = 1.0 - insideWeight;
294 const auto upwindTerm = rho*insideVolVars.mobility(0)*insideWeight + rho*outsideVolVars.mobility(0)*outsideWeight;
296 const auto insideFluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, insideScv, InvalidElemSol{});
297 const auto outsideFluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(outsideElement, outsideScv, InvalidElemSol{});
300 const auto insideSw = insideVolVars.saturation(0);
301 const auto outsideSw = outsideVolVars.saturation(0);
302 const auto insidePc = insideVolVars.capillaryPressure();
303 const auto outsidePc = outsideVolVars.capillaryPressure();
304 const auto dkrw_dsw_inside = insideFluidMatrixInteraction.dkrw_dsw(insideSw);
305 const auto dkrw_dsw_outside = outsideFluidMatrixInteraction.dkrw_dsw(outsideSw);
306 const auto dsw_dpw_inside = -insideFluidMatrixInteraction.dsw_dpc(insidePc);
307 const auto dsw_dpw_outside = -outsideFluidMatrixInteraction.dsw_dpc(outsidePc);
310 const auto tij = elemFluxVarsCache[scvf].advectionTij();
313 auto& dI_dI = derivativeMatrices[insideScvIdx];
314 auto& dI_dJ = derivativeMatrices[outsideScvIdx];
317 dI_dI[conti0EqIdx][0] += tij*upwindTerm + rho_mu*flux*insideWeight*dkrw_dsw_inside*dsw_dpw_inside;
318 dI_dJ[conti0EqIdx][0] += -tij*upwindTerm + rho_mu*flux*outsideWeight*dkrw_dsw_outside*dsw_dpw_outside;
332 template<
class JacobianMatrix,
class T = TypeTag>
335 const Problem& problem,
336 const Element& element,
337 const FVElementGeometry& fvGeometry,
338 const ElementVolumeVariables& curElemVolVars,
339 const ElementFluxVariablesCache& elemFluxVarsCache,
340 const SubControlVolumeFace& scvf)
const
342 static_assert(!enableWaterDiffusionInAir,
343 "richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
344 static_assert(!FluidSystem::isCompressible(0),
345 "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
346 static_assert(FluidSystem::viscosityIsConstant(0),
347 "richards/localresidual.hh: Analytic Jacobian only supports fluids with constant viscosity!");
350 const auto insideScvIdx = scvf.insideScvIdx();
351 const auto outsideScvIdx = scvf.outsideScvIdx();
352 const auto& insideScv = fvGeometry.scv(insideScvIdx);
353 const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
354 const auto& insideVolVars = curElemVolVars[insideScvIdx];
355 const auto& outsideVolVars = curElemVolVars[outsideScvIdx];
358 static const auto rho = insideVolVars.density(0);
359 static const auto mu = insideVolVars.viscosity(0);
360 static const auto rho_mu = rho/mu;
365 static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(),
"Flux.UpwindWeight");
366 const auto flux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache);
367 const auto insideWeight = std::signbit(flux) ? (1.0 - upwindWeight) : upwindWeight;
368 const auto outsideWeight = 1.0 - insideWeight;
369 const auto upwindTerm = rho*insideVolVars.mobility(0)*insideWeight + rho*outsideVolVars.mobility(0)*outsideWeight;
371 const auto insideFluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, insideScv, InvalidElemSol{});
372 const auto outsideFluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, outsideScv, InvalidElemSol{});
375 const auto insideSw = insideVolVars.saturation(0);
376 const auto outsideSw = outsideVolVars.saturation(0);
377 const auto insidePc = insideVolVars.capillaryPressure();
378 const auto outsidePc = outsideVolVars.capillaryPressure();
379 const auto dkrw_dsw_inside = insideFluidMatrixInteraction.dkrw_dsw(insideSw);
380 const auto dkrw_dsw_outside = outsideFluidMatrixInteraction.dkrw_dsw(outsideSw);
381 const auto dsw_dpw_inside = -insideFluidMatrixInteraction.dsw_dpc(insidePc);
382 const auto dsw_dpw_outside = -outsideFluidMatrixInteraction.dsw_dpc(outsidePc);
386 const auto ti = AdvectionType::calculateTransmissibilities(
387 problem, element, fvGeometry, curElemVolVars, scvf, elemFluxVarsCache[scvf]
391 auto& dI_dJ_inside = A[insideScv.dofIndex()];
392 auto& dI_dJ_outside = A[outsideScv.dofIndex()];
395 for (
const auto& scvJ : scvs(fvGeometry))
398 const auto& tj = ti[scvJ.indexInElement()];
399 const auto globalJ = scvJ.dofIndex();
402 dI_dJ_inside[globalJ][conti0EqIdx][0] += tj*upwindTerm;
403 dI_dJ_outside[globalJ][conti0EqIdx][0] += -tj*upwindTerm;
406 const auto upwindContributionInside = rho_mu*flux*insideWeight*dkrw_dsw_inside*dsw_dpw_inside;
407 const auto upwindContributionOutside = rho_mu*flux*outsideWeight*dkrw_dsw_outside*dsw_dpw_outside;
410 A[insideScv.dofIndex()][insideScv.dofIndex()][conti0EqIdx][0] += upwindContributionInside;
411 A[insideScv.dofIndex()][outsideScv.dofIndex()][conti0EqIdx][0] += upwindContributionOutside;
413 A[outsideScv.dofIndex()][insideScv.dofIndex()][conti0EqIdx][0] -= upwindContributionInside;
414 A[outsideScv.dofIndex()][outsideScv.dofIndex()][conti0EqIdx][0] -= upwindContributionOutside;
431 template<
class PartialDerivativeMatrices>
433 const Problem& problem,
434 const Element& element,
435 const FVElementGeometry& fvGeometry,
436 const ElementVolumeVariables& curElemVolVars,
437 const ElementFluxVariablesCache& elemFluxVarsCache,
438 const SubControlVolumeFace& scvf)
const
440 static_assert(!enableWaterDiffusionInAir,
441 "richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
442 static_assert(!FluidSystem::isCompressible(0),
443 "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
444 static_assert(FluidSystem::viscosityIsConstant(0),
445 "richards/localresidual.hh: Analytic Jacobian only supports fluids with constant viscosity!");
449 const auto insideScvIdx = scvf.insideScvIdx();
450 const auto& insideScv = fvGeometry.scv(insideScvIdx);
451 const auto& insideVolVars = curElemVolVars[insideScvIdx];
452 const auto& outsideVolVars = curElemVolVars[scvf.outsideScvIdx()];
453 const auto insideFluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, insideScv, InvalidElemSol{});
456 static const auto rho = insideVolVars.density(0);
457 static const auto mu = insideVolVars.viscosity(0);
458 static const auto rho_mu = rho/mu;
463 static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(),
"Flux.UpwindWeight");
464 const auto flux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache);
465 const auto insideWeight = std::signbit(flux) ? (1.0 - upwindWeight) : upwindWeight;
466 const auto outsideWeight = 1.0 - insideWeight;
467 const auto upwindTerm = rho*insideVolVars.mobility(0)*insideWeight + rho*outsideVolVars.mobility(0)*outsideWeight;
470 const auto insideSw = insideVolVars.saturation(0);
471 const auto insidePc = insideVolVars.capillaryPressure();
472 const auto dkrw_dsw_inside = insideFluidMatrixInteraction.dkrw_dsw(insideSw);
473 const auto dsw_dpw_inside = -insideFluidMatrixInteraction.dsw_dpc(insidePc);
476 const auto tij = elemFluxVarsCache[scvf].advectionTij();
479 derivativeMatrices[insideScvIdx][conti0EqIdx][0] += tij*upwindTerm + rho_mu*flux*insideWeight*dkrw_dsw_inside*dsw_dpw_inside;
493 template<
class PartialDerivativeMatrices>
495 const Problem& problem,
496 const Element& element,
497 const FVElementGeometry& fvGeometry,
498 const ElementVolumeVariables& curElemVolVars,
499 const ElementFluxVariablesCache& elemFluxVarsCache,
500 const SubControlVolumeFace& scvf)
const
503 PartialDerivativeMatrices&, Element, FVElementGeometry,
504 ElementVolumeVariables, ElementFluxVariablesCache, SubControlVolumeFace>()
506 problem.addRobinFluxDerivatives(derivativeMatrices, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);
510 Implementation *asImp_()
511 {
return static_cast<Implementation *
> (
this); }
513 const Implementation *asImp_()
const
514 {
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.
Helper classes to compute the integration elements.
The available discretization methods in Dumux.
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
Distance implementation details.
Definition: fclocalassembler.hh:42
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
constexpr CCTpfa cctpfa
Definition: method.hh:137
constexpr Box box
Definition: method.hh:139
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:35
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:234
std::enable_if_t< GetPropType< T, Properties::GridGeometry >::discMethod==DiscretizationMethods::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:258
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:432
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:115
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:200
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:152
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:494
std::enable_if_t< GetPropType< T, Properties::GridGeometry >::discMethod==DiscretizationMethods::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:334
Declares all properties used in Dumux.