187 const auto& gradCalc = elemCtx.gradientCalculator();
190 const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(faceIdx);
191 const auto& faceNormal = scvf.normal();
193 unsigned i = scvf.interiorIndex();
194 unsigned j = scvf.exteriorIndex();
195 interiorDofIdx_ =
static_cast<short>(i);
196 exteriorDofIdx_ =
static_cast<short>(j);
197 unsigned focusDofIdx = elemCtx.focusDofIndex();
200 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
201 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
202 Valgrind::SetUndefined(potentialGrad_[phaseIdx]);
207 gradCalc.calculateGradient(potentialGrad_[phaseIdx],
211 Valgrind::CheckDefined(potentialGrad_[phaseIdx]);
215 if (Parameters::Get<Parameters::EnableGravity>()) {
218 const auto& gIn = elemCtx.problem().gravity(elemCtx, i, timeIdx);
219 const auto& gEx = elemCtx.problem().gravity(elemCtx, j, timeIdx);
221 const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
222 const auto& intQuantsEx = elemCtx.intensiveQuantities(j, timeIdx);
224 const auto& posIn = elemCtx.pos(i, timeIdx);
225 const auto& posEx = elemCtx.pos(j, timeIdx);
226 const auto& posFace = scvf.integrationPos();
229 DimVector distVecIn(posIn);
230 DimVector distVecEx(posEx);
231 DimVector distVecTotal(posEx);
233 distVecIn -= posFace;
234 distVecEx -= posFace;
235 distVecTotal -= posIn;
236 Scalar absDistTotalSquared = distVecTotal.two_norm2();
237 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
238 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
244 if (std::is_same<Scalar, Evaluation>::value ||
245 interiorDofIdx_ ==
static_cast<int>(focusDofIdx))
247 const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
248 pStatIn = - rhoIn*(gIn*distVecIn);
251 Scalar rhoIn = Toolbox::value(intQuantsIn.fluidState().density(phaseIdx));
252 pStatIn = - rhoIn*(gIn*distVecIn);
259 if (std::is_same<Scalar, Evaluation>::value ||
260 exteriorDofIdx_ ==
static_cast<int>(focusDofIdx))
262 const Evaluation& rhoEx = intQuantsEx.fluidState().density(phaseIdx);
263 pStatEx = - rhoEx*(gEx*distVecEx);
266 Scalar rhoEx = Toolbox::value(intQuantsEx.fluidState().density(phaseIdx));
267 pStatEx = - rhoEx*(gEx*distVecEx);
274 Dune::FieldVector<Evaluation, dimWorld> f(distVecTotal);
275 f *= (pStatEx - pStatIn)/absDistTotalSquared;
278 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
279 potentialGrad_[phaseIdx][dimIdx] += f[dimIdx];
281 for (
unsigned dimIdx = 0; dimIdx < potentialGrad_[phaseIdx].size(); ++dimIdx) {
282 if (!isfinite(potentialGrad_[phaseIdx][dimIdx])) {
283 throw NumericalProblem(
"Non-finite potential gradient for phase '"
284 + std::string(FluidSystem::phaseName(phaseIdx))+
"'");
290 Valgrind::SetUndefined(K_);
291 elemCtx.problem().intersectionIntrinsicPermeability(K_, elemCtx, faceIdx, timeIdx);
292 Valgrind::CheckDefined(K_);
294 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
295 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
296 Valgrind::SetUndefined(potentialGrad_[phaseIdx]);
301 Evaluation tmp = 0.0;
302 for (
unsigned dimIdx = 0; dimIdx < faceNormal.size(); ++dimIdx)
303 tmp += potentialGrad_[phaseIdx][dimIdx]*faceNormal[dimIdx];
306 upstreamDofIdx_[phaseIdx] = exteriorDofIdx_;
307 downstreamDofIdx_[phaseIdx] = interiorDofIdx_;
310 upstreamDofIdx_[phaseIdx] = interiorDofIdx_;
311 downstreamDofIdx_[phaseIdx] = exteriorDofIdx_;
316 const auto& up = elemCtx.intensiveQuantities(upstreamDofIdx_[phaseIdx], timeIdx);
317 if (upstreamDofIdx_[phaseIdx] ==
static_cast<int>(focusDofIdx))
318 mobility_[phaseIdx] = up.mobility(phaseIdx);
320 mobility_[phaseIdx] = Toolbox::value(up.mobility(phaseIdx));
332 unsigned boundaryFaceIdx,
334 const FluidState& fluidState)
336 const auto& gradCalc = elemCtx.gradientCalculator();
340 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
341 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
342 Valgrind::SetUndefined(potentialGrad_[phaseIdx]);
347 gradCalc.calculateBoundaryGradient(potentialGrad_[phaseIdx],
351 Valgrind::CheckDefined(potentialGrad_[phaseIdx]);
354 const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(boundaryFaceIdx);
355 auto i = scvf.interiorIndex();
356 interiorDofIdx_ =
static_cast<short>(i);
357 exteriorDofIdx_ = -1;
358 int focusDofIdx = elemCtx.focusDofIndex();
361 const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
362 K_ = intQuantsIn.intrinsicPermeability();
365 if (Parameters::Get<Parameters::EnableGravity>()) {
368 const auto& gIn = elemCtx.problem().gravity(elemCtx, i, timeIdx);
369 const auto& posIn = elemCtx.pos(i, timeIdx);
370 const auto& posFace = scvf.integrationPos();
373 DimVector distVecIn(posIn);
374 distVecIn -= posFace;
375 Scalar absDistSquared = distVecIn.two_norm2();
376 Scalar gTimesDist = gIn*distVecIn;
378 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
379 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
383 Evaluation rhoIn = intQuantsIn.fluidState().density(phaseIdx);
384 Evaluation pStatIn = - gTimesDist*rhoIn;
386 Valgrind::CheckDefined(pStatIn);
392 EvalDimVector f(distVecIn);
393 f *= pStatIn / absDistSquared;
396 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
397 potentialGrad_[phaseIdx][dimIdx] += f[dimIdx];
399 Valgrind::CheckDefined(potentialGrad_[phaseIdx]);
400 for (
unsigned dimIdx = 0; dimIdx < potentialGrad_[phaseIdx].size(); ++dimIdx) {
401 if (!isfinite(potentialGrad_[phaseIdx][dimIdx])) {
402 throw NumericalProblem(
"Non-finite potential gradient for phase '"
403 + std::string(FluidSystem::phaseName(phaseIdx))+
"'");
410 const auto& faceNormal = scvf.normal();
412 const auto& matParams = elemCtx.problem().materialLawParams(elemCtx, i, timeIdx);
414 Scalar kr[numPhases];
415 MaterialLaw::relativePermeabilities(kr, matParams, fluidState);
416 Valgrind::CheckDefined(kr);
418 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
419 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
422 Evaluation tmp = 0.0;
423 for (
unsigned dimIdx = 0; dimIdx < faceNormal.size(); ++dimIdx)
424 tmp += potentialGrad_[phaseIdx][dimIdx]*faceNormal[dimIdx];
427 upstreamDofIdx_[phaseIdx] = exteriorDofIdx_;
428 downstreamDofIdx_[phaseIdx] = interiorDofIdx_;
431 upstreamDofIdx_[phaseIdx] = interiorDofIdx_;
432 downstreamDofIdx_[phaseIdx] = exteriorDofIdx_;
436 if (upstreamDofIdx_[phaseIdx] < 0) {
437 if (interiorDofIdx_ == focusDofIdx)
438 mobility_[phaseIdx] =
439 kr[phaseIdx] / fluidState.viscosity(phaseIdx);
441 mobility_[phaseIdx] =
442 Toolbox::value(kr[phaseIdx])
443 / Toolbox::value(fluidState.viscosity(phaseIdx));
445 else if (upstreamDofIdx_[phaseIdx] != focusDofIdx)
446 mobility_[phaseIdx] = Toolbox::value(intQuantsIn.mobility(phaseIdx));
448 mobility_[phaseIdx] = intQuantsIn.mobility(phaseIdx);
449 Valgrind::CheckDefined(mobility_[phaseIdx]);
461 const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
462 const DimVector& normal = scvf.normal();
463 Valgrind::CheckDefined(normal);
465 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
466 filterVelocity_[phaseIdx] = 0.0;
467 volumeFlux_[phaseIdx] = 0.0;
468 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
471 asImp_().calculateFilterVelocity_(phaseIdx);
472 Valgrind::CheckDefined(filterVelocity_[phaseIdx]);
474 volumeFlux_[phaseIdx] = 0.0;
475 for (
unsigned i = 0; i < normal.size(); ++i)
476 volumeFlux_[phaseIdx] += filterVelocity_[phaseIdx][i] * normal[i];
487 unsigned boundaryFaceIdx,
490 const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(boundaryFaceIdx);
491 const DimVector& normal = scvf.normal();
492 Valgrind::CheckDefined(normal);
494 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
495 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
496 filterVelocity_[phaseIdx] = 0.0;
497 volumeFlux_[phaseIdx] = 0.0;
501 asImp_().calculateFilterVelocity_(phaseIdx);
502 Valgrind::CheckDefined(filterVelocity_[phaseIdx]);
503 volumeFlux_[phaseIdx] = 0.0;
504 for (
unsigned i = 0; i < normal.size(); ++i)
505 volumeFlux_[phaseIdx] += filterVelocity_[phaseIdx][i] * normal[i];