26#ifndef DUMUX_RICHARDS_LOCAL_RESIDUAL_HH
27#define DUMUX_RICHARDS_LOCAL_RESIDUAL_HH
42template<
class TypeTag>
56 using SubControlVolume =
typename FVElementGeometry::SubControlVolume;
57 using SubControlVolumeFace =
typename FVElementGeometry::SubControlVolumeFace;
59 using Element =
typename GridView::template Codim<0>::Entity;
64 enum { conti0EqIdx = Indices::conti0EqIdx };
68 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
69 gasPhaseIdx = FluidSystem::gasPhaseIdx,
70 liquidCompIdx = FluidSystem::liquidCompIdx
73 static constexpr bool enableWaterDiffusionInAir
74 = getPropValue<TypeTag, Properties::EnableWaterDiffusionInAir>();
80 double operator[] (
const Index i)
const
81 {
static_assert(
AlwaysFalse<Index>::value,
"Solution-dependent material parameters not supported with analytical differentiation");
return 0.0; }
85 using ParentType::ParentType;
99 const SubControlVolume& scv,
100 const VolumeVariables& volVars)
const
103 NumEqVector storage(0.0);
104 storage[conti0EqIdx] = volVars.porosity()
105 * volVars.density(liquidPhaseIdx)
106 * volVars.saturation(liquidPhaseIdx);
109 if (enableWaterDiffusionInAir)
110 storage[conti0EqIdx] += volVars.porosity()
111 * volVars.molarDensity(gasPhaseIdx)
112 * volVars.moleFraction(gasPhaseIdx, liquidCompIdx)
113 * FluidSystem::molarMass(liquidCompIdx)
114 * volVars.saturation(gasPhaseIdx);
117 EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, liquidPhaseIdx);
118 EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, gasPhaseIdx);
119 EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
136 const Element& element,
137 const FVElementGeometry& fvGeometry,
138 const ElementVolumeVariables& elemVolVars,
139 const SubControlVolumeFace& scvf,
140 const ElementFluxVariablesCache& elemFluxVarsCache)
const
142 FluxVariables fluxVars;
143 fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
145 NumEqVector flux(0.0);
147 auto upwindTerm = [](
const auto& volVars)
148 {
return volVars.density(liquidPhaseIdx)*volVars.mobility(liquidPhaseIdx); };
150 flux[conti0EqIdx] = fluxVars.advectiveFlux(liquidPhaseIdx, upwindTerm);
153 if (enableWaterDiffusionInAir)
157 flux[conti0EqIdx] += fluxVars.molecularDiffusionFlux(gasPhaseIdx)[liquidCompIdx];
159 flux[conti0EqIdx] += fluxVars.molecularDiffusionFlux(gasPhaseIdx)[liquidCompIdx]*FluidSystem::molarMass(liquidCompIdx);
163 EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, liquidPhaseIdx);
167 EnergyLocalResidual::heatConductionFlux(flux, fluxVars);
182 template<
class PartialDerivativeMatrix>
184 const Problem& problem,
185 const Element& element,
186 const FVElementGeometry& fvGeometry,
187 const VolumeVariables& curVolVars,
188 const SubControlVolume& scv)
const
190 static_assert(!enableWaterDiffusionInAir,
191 "richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
192 static_assert(!FluidSystem::isCompressible(0),
193 "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
195 const auto poreVolume = scv.volume()*curVolVars.porosity();
196 static const auto rho = curVolVars.density(0);
199 const auto& matParams = problem.spatialParams().materialLawParams(element, scv, InvalidElemSol{});
203 using MaterialLaw =
typename Problem::SpatialParams::MaterialLaw;
204 partialDerivatives[conti0EqIdx][0] += -rho*poreVolume/this->timeLoop().timeStepSize()*MaterialLaw::dsw_dpc(matParams, curVolVars.capillaryPressure());
219 template<
class PartialDerivativeMatrix>
221 const Problem& problem,
222 const Element& element,
223 const FVElementGeometry& fvGeometry,
224 const VolumeVariables& curVolVars,
225 const SubControlVolume& scv)
const
242 template<
class PartialDerivativeMatrices,
class T = TypeTag>
245 const Problem& problem,
246 const Element& element,
247 const FVElementGeometry& fvGeometry,
248 const ElementVolumeVariables& curElemVolVars,
249 const ElementFluxVariablesCache& elemFluxVarsCache,
250 const SubControlVolumeFace& scvf)
const
252 static_assert(!enableWaterDiffusionInAir,
253 "richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
254 static_assert(!FluidSystem::isCompressible(0),
255 "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
256 static_assert(!FluidSystem::viscosityIsConstant(0),
257 "richards/localresidual.hh: Analytic Jacobian only supports fluids with constant viscosity!");
260 const auto insideScvIdx = scvf.insideScvIdx();
261 const auto outsideScvIdx = scvf.outsideScvIdx();
262 const auto outsideElement = fvGeometry.gridGeometry().element(outsideScvIdx);
263 const auto& insideScv = fvGeometry.scv(insideScvIdx);
264 const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
265 const auto& insideVolVars = curElemVolVars[insideScvIdx];
266 const auto& outsideVolVars = curElemVolVars[outsideScvIdx];
267 const auto& insideMaterialParams = problem.spatialParams().materialLawParams(element, insideScv, InvalidElemSol{});
268 const auto& outsideMaterialParams = problem.spatialParams().materialLawParams(outsideElement, outsideScv, InvalidElemSol{});
271 static const auto rho = insideVolVars.density(0);
272 static const auto mu = insideVolVars.viscosity(0);
273 static const auto rho_mu = rho/mu;
278 static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(),
"Flux.UpwindWeight");
279 const auto flux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache);
280 const auto insideWeight = std::signbit(flux) ? (1.0 - upwindWeight) : upwindWeight;
281 const auto outsideWeight = 1.0 - insideWeight;
282 const auto upwindTerm = rho*insideVolVars.mobility(0)*insideWeight + rho*outsideVolVars.mobility(0)*outsideWeight;
285 const auto insideSw = insideVolVars.saturation(0);
286 const auto outsideSw = outsideVolVars.saturation(0);
287 const auto insidePc = insideVolVars.capillaryPressure();
288 const auto outsidePc = outsideVolVars.capillaryPressure();
289 using MaterialLaw =
typename Problem::SpatialParams::MaterialLaw;
290 const auto dkrw_dsw_inside = MaterialLaw::dkrw_dsw(insideMaterialParams, insideSw);
291 const auto dkrw_dsw_outside = MaterialLaw::dkrw_dsw(outsideMaterialParams, outsideSw);
292 const auto dsw_dpw_inside = -MaterialLaw::dsw_dpc(insideMaterialParams, insidePc);
293 const auto dsw_dpw_outside = -MaterialLaw::dsw_dpc(outsideMaterialParams, outsidePc);
296 const auto tij = elemFluxVarsCache[scvf].advectionTij();
299 auto& dI_dI = derivativeMatrices[insideScvIdx];
300 auto& dI_dJ = derivativeMatrices[outsideScvIdx];
303 dI_dI[conti0EqIdx][0] += tij*upwindTerm + rho_mu*flux*insideWeight*dkrw_dsw_inside*dsw_dpw_inside;
304 dI_dJ[conti0EqIdx][0] += -tij*upwindTerm + rho_mu*flux*outsideWeight*dkrw_dsw_outside*dsw_dpw_outside;
318 template<
class JacobianMatrix,
class T = TypeTag>
321 const Problem& problem,
322 const Element& element,
323 const FVElementGeometry& fvGeometry,
324 const ElementVolumeVariables& curElemVolVars,
325 const ElementFluxVariablesCache& elemFluxVarsCache,
326 const SubControlVolumeFace& scvf)
const
328 static_assert(!enableWaterDiffusionInAir,
329 "richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
330 static_assert(!FluidSystem::isCompressible(0),
331 "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
332 static_assert(!FluidSystem::viscosityIsConstant(0),
333 "richards/localresidual.hh: Analytic Jacobian only supports fluids with constant viscosity!");
336 const auto insideScvIdx = scvf.insideScvIdx();
337 const auto outsideScvIdx = scvf.outsideScvIdx();
338 const auto& insideScv = fvGeometry.scv(insideScvIdx);
339 const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
340 const auto& insideVolVars = curElemVolVars[insideScvIdx];
341 const auto& outsideVolVars = curElemVolVars[outsideScvIdx];
342 const auto& insideMaterialParams = problem.spatialParams().materialLawParams(element, insideScv, InvalidElemSol{});
343 const auto& outsideMaterialParams = problem.spatialParams().materialLawParams(element, outsideScv, InvalidElemSol{});
346 static const auto rho = insideVolVars.density(0);
347 static const auto mu = insideVolVars.viscosity(0);
348 static const auto rho_mu = rho/mu;
353 static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(),
"Flux.UpwindWeight");
354 const auto flux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache);
355 const auto insideWeight = std::signbit(flux) ? (1.0 - upwindWeight) : upwindWeight;
356 const auto outsideWeight = 1.0 - insideWeight;
357 const auto upwindTerm = rho*insideVolVars.mobility(0)*insideWeight + rho*outsideVolVars.mobility(0)*outsideWeight;
360 const auto insideSw = insideVolVars.saturation(0);
361 const auto outsideSw = outsideVolVars.saturation(0);
362 const auto insidePc = insideVolVars.capillaryPressure();
363 const auto outsidePc = outsideVolVars.capillaryPressure();
364 using MaterialLaw =
typename Problem::SpatialParams::MaterialLaw;
365 const auto dkrw_dsw_inside = MaterialLaw::dkrw_dsw(insideMaterialParams, insideSw);
366 const auto dkrw_dsw_outside = MaterialLaw::dkrw_dsw(outsideMaterialParams, outsideSw);
367 const auto dsw_dpw_inside = -MaterialLaw::dsw_dpc(insideMaterialParams, insidePc);
368 const auto dsw_dpw_outside = -MaterialLaw::dsw_dpc(outsideMaterialParams, outsidePc);
372 const auto ti = AdvectionType::calculateTransmissibilities(
373 problem, element, fvGeometry, curElemVolVars, scvf, elemFluxVarsCache[scvf]
377 auto& dI_dJ_inside = A[insideScv.dofIndex()];
378 auto& dI_dJ_outside = A[outsideScv.dofIndex()];
381 for (
const auto& scvJ : scvs(fvGeometry))
384 const auto& tj = ti[scvJ.indexInElement()];
385 const auto globalJ = scvJ.dofIndex();
388 dI_dJ_inside[globalJ][conti0EqIdx][0] += tj*upwindTerm;
389 dI_dJ_outside[globalJ][conti0EqIdx][0] += -tj*upwindTerm;
392 const auto upwindContributionInside = rho_mu*flux*insideWeight*dkrw_dsw_inside*dsw_dpw_inside;
393 const auto upwindContributionOutside = rho_mu*flux*outsideWeight*dkrw_dsw_outside*dsw_dpw_outside;
396 A[insideScv.dofIndex()][insideScv.dofIndex()][conti0EqIdx][0] += upwindContributionInside;
397 A[insideScv.dofIndex()][outsideScv.dofIndex()][conti0EqIdx][0] += upwindContributionOutside;
399 A[outsideScv.dofIndex()][insideScv.dofIndex()][conti0EqIdx][0] -= upwindContributionInside;
400 A[outsideScv.dofIndex()][outsideScv.dofIndex()][conti0EqIdx][0] -= upwindContributionOutside;
417 template<
class PartialDerivativeMatrices>
419 const Problem& problem,
420 const Element& element,
421 const FVElementGeometry& fvGeometry,
422 const ElementVolumeVariables& curElemVolVars,
423 const ElementFluxVariablesCache& elemFluxVarsCache,
424 const SubControlVolumeFace& scvf)
const
426 static_assert(!enableWaterDiffusionInAir,
427 "richards/localresidual.hh: Analytic Jacobian not implemented for the water diffusion in air version!");
428 static_assert(!FluidSystem::isCompressible(0),
429 "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!");
430 static_assert(!FluidSystem::viscosityIsConstant(0),
431 "richards/localresidual.hh: Analytic Jacobian only supports fluids with constant viscosity!");
435 const auto insideScvIdx = scvf.insideScvIdx();
436 const auto& insideScv = fvGeometry.scv(insideScvIdx);
437 const auto& insideVolVars = curElemVolVars[insideScvIdx];
438 const auto& outsideVolVars = curElemVolVars[scvf.outsideScvIdx()];
439 const auto& insideMaterialParams = problem.spatialParams().materialLawParams(element, insideScv, InvalidElemSol{});
442 static const auto rho = insideVolVars.density(0);
443 static const auto mu = insideVolVars.viscosity(0);
444 static const auto rho_mu = rho/mu;
449 static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(),
"Flux.UpwindWeight");
450 const auto flux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache);
451 const auto insideWeight = std::signbit(flux) ? (1.0 - upwindWeight) : upwindWeight;
452 const auto outsideWeight = 1.0 - insideWeight;
453 const auto upwindTerm = rho*insideVolVars.mobility(0)*insideWeight + rho*outsideVolVars.mobility(0)*outsideWeight;
456 const auto insideSw = insideVolVars.saturation(0);
457 const auto insidePc = insideVolVars.capillaryPressure();
458 using MaterialLaw =
typename Problem::SpatialParams::MaterialLaw;
459 const auto dkrw_dsw_inside = MaterialLaw::dkrw_dsw(insideMaterialParams, insideSw);
460 const auto dsw_dpw_inside = -MaterialLaw::dsw_dpc(insideMaterialParams, insidePc);
463 const auto tij = elemFluxVarsCache[scvf].advectionTij();
466 derivativeMatrices[insideScvIdx][conti0EqIdx][0] += tij*upwindTerm + rho_mu*flux*insideWeight*dkrw_dsw_inside*dsw_dpw_inside;
480 template<
class PartialDerivativeMatrices>
482 const Problem& problem,
483 const Element& element,
484 const FVElementGeometry& fvGeometry,
485 const ElementVolumeVariables& curElemVolVars,
486 const ElementFluxVariablesCache& elemFluxVarsCache,
487 const SubControlVolumeFace& scvf)
const
491 Implementation *asImp_()
492 {
return static_cast<Implementation *
> (
this); }
494 const Implementation *asImp_()
const
495 {
return static_cast<const Implementation *
> (
this); }
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
The available discretization methods in Dumux.
The reference frameworks and formulations available for splitting total fluxes into a advective and d...
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:37
Element-wise calculation of the Jacobian matrix for problems using the Richards fully implicit models...
Definition: porousmediumflow/richards/localresidual.hh:44
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 non-wetting phase.
Definition: porousmediumflow/richards/localresidual.hh:220
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 non-wetting phase.
Definition: porousmediumflow/richards/localresidual.hh:418
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:98
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:183
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:135
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 non-wetting phase for cell-centered FVM using TPFA.
Definition: porousmediumflow/richards/localresidual.hh:244
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 non-wetting phase.
Definition: porousmediumflow/richards/localresidual.hh:481
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:320
Declares all properties used in Dumux.