225 enum { dimWorld = GridView::dimensionworld };
228 using Toolbox = MathToolbox<Evaluation>;
230 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
231 using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
232 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
233 using DimEvalMatrix = Dune::FieldMatrix<Evaluation, dimWorld, dimWorld>;
240 {
return ergunCoefficient_; }
249 {
return mobilityPassabilityRatio_[phaseIdx]; }
252 void calculateGradients_(
const ElementContext& elemCtx,
258 auto focusDofIdx = elemCtx.focusDofIndex();
259 unsigned i =
static_cast<unsigned>(this->interiorDofIdx_);
260 unsigned j =
static_cast<unsigned>(this->exteriorDofIdx_);
261 const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
262 const auto& intQuantsEx = elemCtx.intensiveQuantities(j, timeIdx);
267 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
268 sqrtK_[dimIdx] = std::sqrt(this->K_[dimIdx][dimIdx]);
271 if (focusDofIdx == i) {
273 (intQuantsIn.ergunCoefficient() +
274 getValue(intQuantsEx.ergunCoefficient()))/2;
276 else if (focusDofIdx == j)
278 (getValue(intQuantsIn.ergunCoefficient()) +
279 intQuantsEx.ergunCoefficient())/2;
282 (getValue(intQuantsIn.ergunCoefficient()) +
283 getValue(intQuantsEx.ergunCoefficient()))/2;
286 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
287 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
290 unsigned upIdx =
static_cast<unsigned>(this->upstreamIndex_(phaseIdx));
291 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
293 if (focusDofIdx == upIdx) {
295 up.fluidState().density(phaseIdx);
296 mobilityPassabilityRatio_[phaseIdx] =
297 up.mobilityPassabilityRatio(phaseIdx);
301 getValue(up.fluidState().density(phaseIdx));
302 mobilityPassabilityRatio_[phaseIdx] =
303 getValue(up.mobilityPassabilityRatio(phaseIdx));
308 template <
class Flu
idState>
309 void calculateBoundaryGradients_(
const ElementContext& elemCtx,
310 unsigned boundaryFaceIdx,
312 const FluidState& fluidState)
319 auto focusDofIdx = elemCtx.focusDofIndex();
320 unsigned i =
static_cast<unsigned>(this->interiorDofIdx_);
321 const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
325 if (focusDofIdx == i)
326 ergunCoefficient_ = intQuantsIn.ergunCoefficient();
328 ergunCoefficient_ = getValue(intQuantsIn.ergunCoefficient());
333 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
334 sqrtK_[dimIdx] = std::sqrt(this->K_[dimIdx][dimIdx]);
336 for (
unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
337 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
340 if (focusDofIdx == i) {
341 density_[phaseIdx] = intQuantsIn.fluidState().density(phaseIdx);
342 mobilityPassabilityRatio_[phaseIdx] = intQuantsIn.mobilityPassabilityRatio(phaseIdx);
346 getValue(intQuantsIn.fluidState().density(phaseIdx));
347 mobilityPassabilityRatio_[phaseIdx] =
348 getValue(intQuantsIn.mobilityPassabilityRatio(phaseIdx));
361 auto focusDofIdx = elemCtx.focusDofIndex();
362 auto i = asImp_().interiorIndex();
363 auto j = asImp_().exteriorIndex();
364 const auto& intQuantsI = elemCtx.intensiveQuantities(i, timeIdx);
365 const auto& intQuantsJ = elemCtx.intensiveQuantities(j, timeIdx);
367 const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
368 const auto& normal = scvf.normal();
369 Valgrind::CheckDefined(normal);
373 if (focusDofIdx == i)
375 (intQuantsI.ergunCoefficient() +
376 getValue(intQuantsJ.ergunCoefficient())) / 2;
377 else if (focusDofIdx == j)
379 (getValue(intQuantsI.ergunCoefficient()) +
380 intQuantsJ.ergunCoefficient()) / 2;
383 (getValue(intQuantsI.ergunCoefficient()) +
384 getValue(intQuantsJ.ergunCoefficient())) / 2;
389 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; phaseIdx++) {
390 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
391 this->filterVelocity_[phaseIdx] = 0.0;
392 this->volumeFlux_[phaseIdx] = 0.0;
396 calculateForchheimerFlux_(phaseIdx);
398 this->volumeFlux_[phaseIdx] = 0.0;
399 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++ dimIdx)
400 this->volumeFlux_[phaseIdx] +=
401 this->filterVelocity_[phaseIdx][dimIdx]*normal[dimIdx];
412 const auto& boundaryFace = elemCtx.stencil(timeIdx).boundaryFace(bfIdx);
413 const auto& normal = boundaryFace.normal();
418 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; phaseIdx++) {
419 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
420 this->filterVelocity_[phaseIdx] = 0.0;
421 this->volumeFlux_[phaseIdx] = 0.0;
425 calculateForchheimerFlux_(phaseIdx);
427 this->volumeFlux_[phaseIdx] = 0.0;
428 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
429 this->volumeFlux_[phaseIdx] +=
430 this->filterVelocity_[phaseIdx][dimIdx]*normal[dimIdx];
434 void calculateForchheimerFlux_(
unsigned phaseIdx)
437 DimEvalVector& velocity = this->filterVelocity_[phaseIdx];
441 DimEvalVector deltaV(1e5);
444 DimEvalVector residual;
446 DimEvalMatrix gradResid;
449 unsigned newtonIter = 0;
450 while (deltaV.one_norm() > 1e-11) {
451 if (newtonIter >= 50)
452 throw NumericalProblem(
"Could not determine Forchheimer velocity within "
453 + std::to_string(newtonIter)+
" iterations");
457 gradForchheimerResid_(residual, gradResid, phaseIdx);
460 gradResid.solve(deltaV, residual);
465 void forchheimerResid_(DimEvalVector& residual,
unsigned phaseIdx)
const
467 const DimEvalVector& velocity = this->filterVelocity_[phaseIdx];
470 const auto& mobility = this->mobility_[phaseIdx];
471 const auto& density = density_[phaseIdx];
475 const auto& pGrad = this->potentialGrad_[phaseIdx];
483 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++ dimIdx)
484 residual[dimIdx] += mobility*pGrad[dimIdx]*this->K_[dimIdx][dimIdx];
497 Evaluation absVel = 0.0;
498 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
499 absVel += velocity[dimIdx]*velocity[dimIdx];
505 absVel = Toolbox::sqrt(absVel);
507 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
508 residual[dimIdx] += sqrtK_[dimIdx]*alpha*velocity[dimIdx];
509 Valgrind::CheckDefined(residual);
512 void gradForchheimerResid_(DimEvalVector& residual,
513 DimEvalMatrix& gradResid,
517 DimEvalVector& velocity = this->filterVelocity_[phaseIdx];
518 forchheimerResid_(residual, phaseIdx);
522 for (
unsigned i = 0; i < dimWorld; ++i) {
523 Scalar coordEps = std::max(eps, Toolbox::scalarValue(velocity[i]) * (1 + eps));
524 velocity[i] += coordEps;
525 forchheimerResid_(tmp, phaseIdx);
529 velocity[i] -= coordEps;
542 for (
unsigned i = 0; i < dimWorld; i++) {
543 for (
unsigned j = 0; j < dimWorld; j++) {
547 if (std::abs(K[i][j]) > 1e-25)
555 Implementation& asImp_()
556 {
return *
static_cast<Implementation *
>(
this); }
558 const Implementation& asImp_()
const
559 {
return *
static_cast<const Implementation *
>(
this); }
566 Evaluation ergunCoefficient_;
569 Evaluation mobilityPassabilityRatio_[numPhases];
572 Evaluation density_[numPhases];