From 45fc9a6cf07fe29da751bd8405592270691c551c Mon Sep 17 00:00:00 2001
From: DennisGlaeser <dennis.glaeser@iws.uni-stuttgart.de>
Date: Fri, 11 Mar 2016 12:02:45 +0100
Subject: [PATCH] [2pdfm][volvars] default initialization fracture variables
The fracture variables and matrix saturations have to be initialized
by default in case the return functions are used for vertices that are
not connected to a fracture
---
.../2pdfm/implicit/volumevariables.hh | 30 ++++++++++++++-----
1 file changed, 23 insertions(+), 7 deletions(-)
diff --git a/dumux/porousmediumflow/2pdfm/implicit/volumevariables.hh b/dumux/porousmediumflow/2pdfm/implicit/volumevariables.hh
index 6f35b5e47e..aba771f9d0 100644
--- a/dumux/porousmediumflow/2pdfm/implicit/volumevariables.hh
+++ b/dumux/porousmediumflow/2pdfm/implicit/volumevariables.hh
@@ -94,13 +94,25 @@ public:
scvIdx,
isOldSol);
+ // initially we set the matrix saturations to the actual solution
+ satNMatrix_ = priVars[saturationIdx];
+ satWMatrix_ = 1.0 - satNMatrix_;
+
+ // fracture variables are set to zero for the case that no fracture is present in the scv
+ satNFracture_ = 0.0;
+ satWFracture_ = 0.0;
+ porosityFracture_ = 0.0;
+ permeabilityFracture_ = 0.0;
+ mobilityFracture_[wPhaseIdx] = 0;
+ mobilityFracture_[nPhaseIdx] = 0;
+
// calculate the matrix mobilities
this->completeFluidState(priVars, problem, element, fvGeometry, scvIdx, fluidState_);
const MaterialLawParams &materialParams = problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
mobilityMatrix_[wPhaseIdx] = MaterialLaw::krw(materialParams, fluidState_.saturation(wPhaseIdx)) / fluidState_.viscosity(wPhaseIdx);
mobilityMatrix_[nPhaseIdx] = MaterialLaw::krn(materialParams, fluidState_.saturation(wPhaseIdx)) / fluidState_.viscosity(nPhaseIdx);
- // update the fracture if we are on a fracture vertex
+ // update the fracture if fracture is present in the scv
if (problem.spatialParams().isVertexFracture(element, scvIdx))
asImp_().updateFracture(priVars, problem, element, fvGeometry, scvIdx, isOldSol);
}
@@ -122,13 +134,14 @@ public:
int scvIdx,
bool isOldSol)
{
+ // the scv is on a fracture, the actual solution corresponds to the solution in the fracture
+ // the matrix saturation will be modified below if interface condition is used
+ satNFracture_ = priVars[saturationIdx];
+ satWFracture_ = 1.0 - satNFracture_;
+
const MaterialLawParams &materialParamsMatrix =
problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
- // initially we set the matrix saturations equal to the fracture saturations
- satNMatrix_ = priVars[saturationIdx];
- satWMatrix_ = 1.0 - satNMatrix_;
-
// the fluid state of the fracture
this->completeFluidState(priVars, problem, element, fvGeometry, scvIdx, fluidStateFracture_);
@@ -243,9 +256,9 @@ public:
Scalar saturationFracture(int phaseIdx) const
{
if (phaseIdx == wPhaseIdx)
- return 1.0 - this->priVar(saturationIdx);
+ return satWFracture_;
else
- return this->priVar(saturationIdx);
+ return satNFracture_;
}
Scalar saturationMatrix(int phaseIdx) const
{
@@ -298,6 +311,9 @@ protected:
Scalar satWMatrix_;
Scalar satNMatrix_;
+ Scalar satWFracture_;
+ Scalar satNFracture_;
+
private:
Implementation &asImp_()
{ return *static_cast<Implementation*>(this); }
--
GitLab