diff --git a/dumux/discretization/box/maxwellstefanslaw.hh b/dumux/discretization/box/maxwellstefanslaw.hh
index ae4a2327202352c9fd44e3fe15da9716c98c0e69..8f1b0f7b91e6f0e4bdb8b3b032f149e2e34e98b2 100644
--- a/dumux/discretization/box/maxwellstefanslaw.hh
+++ b/dumux/discretization/box/maxwellstefanslaw.hh
@@ -151,86 +151,65 @@ private:
const int phaseIdx,
const ComponentFluxVector moleFrac)
{
-
ReducedComponentMatrix reducedDiffusionMatrix(0.0);
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
const auto insideScvIdx = scvf.insideScvIdx();
const auto outsideScvIdx = scvf.outsideScvIdx();
+
//this is to not devide by 0 if the saturation in 0 and the effectiveDiffusivity becomes zero due to that
- if(insideVolVars.saturation(phaseIdx) == 0 || outsideVolVars.saturation(phaseIdx) == 0)
+ if(insideVolVars.saturation(phaseIdx) == 0 || outsideVolVars.saturation(phaseIdx) == 0)
return reducedDiffusionMatrix;
- for (int compIIdx = 0; compIIdx < numComponents-1; compIIdx++)
- {
- // effective diffusion tensors
- using EffDiffModel = typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel);
-
- const auto xi = moleFrac[compIIdx];
+ for (int compIIdx = 0; compIIdx < numComponents-1; compIIdx++)
+ {
+ // effective diffusion tensors
+ using EffDiffModel = typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel);
- //calculate diffusivity for i,numComponents
+ const auto xi = moleFrac[compIIdx];
- auto tinInside = getDiffusionCoefficient(phaseIdx, compIIdx, numComponents-1, problem, element, insideVolVars, fvGeometry.scv(insideScvIdx));
- tinInside = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), tinInside);
- auto tinOutside = getDiffusionCoefficient(phaseIdx, compIIdx, numComponents-1, problem, element, outsideVolVars, fvGeometry.scv(outsideScvIdx));
- tinOutside = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), tinOutside);
+ //calculate diffusivity for i,numComponents
+ auto tinInside = getDiffusionCoefficient(phaseIdx, compIIdx, numComponents-1, problem, element, insideVolVars, fvGeometry.scv(insideScvIdx));
+ tinInside = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), tinInside);
+ auto tinOutside = getDiffusionCoefficient(phaseIdx, compIIdx, numComponents-1, problem, element, outsideVolVars, fvGeometry.scv(outsideScvIdx));
+ tinOutside = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), tinOutside);
// scale by extrusion factor
- tinInside *= insideVolVars.extrusionFactor();
- tinOutside *= outsideVolVars.extrusionFactor();
+ tinInside *= insideVolVars.extrusionFactor();
+ tinOutside *= outsideVolVars.extrusionFactor();
+
+ // the resulting averaged diffusion tensor
+ const auto tin = problem.spatialParams().harmonicMean(tinInside, tinOutside, scvf.unitOuterNormal());
+
+ //begin the entrys of the diffusion matrix of the diagonal
+ reducedDiffusionMatrix[compIIdx][compIIdx] += xi/tin;
+
+ // now set the rest of the entries (off-diagonal and additional entries for diagonal)
+ for (int compJIdx = 0; compJIdx < numComponents-1; compJIdx++)
+ {
+ //we don't want to calculate e.g. water in water diffusion
+ if (compIIdx == compJIdx)
+ continue;
+
+ //calculate diffusivity for compIIdx, compJIdx
+ const auto xj = moleFrac[compJIdx];
+ auto tijInside = getDiffusionCoefficient(phaseIdx, compIIdx, compJIdx, problem, element, insideVolVars, fvGeometry.scv(insideScvIdx));
+ tijInside = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), tijInside);
+ auto tijOutside = getDiffusionCoefficient(phaseIdx, compIIdx, compJIdx, problem, element, outsideVolVars, fvGeometry.scv(outsideScvIdx));
+ tijOutside = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), outsideVolVars.saturation(phaseIdx), tijOutside);
+
+ // scale by extrusion factor
+ tijInside *= insideVolVars.extrusionFactor();
+ tijOutside *= outsideVolVars.extrusionFactor();
// the resulting averaged diffusion tensor
- const auto tin = problem.spatialParams().harmonicMean(tinInside, tinOutside, scvf.unitOuterNormal());
-
- //set the entrys of the diffusion matrix of the diagonal
- reducedDiffusionMatrix[compIIdx][compIIdx] += xi/tin;
-
- for (int compkIdx = 0; compkIdx < numComponents; compkIdx++)
- {
- if (compkIdx == compIIdx)
- continue;
-
- const auto xk = moleFrac[compkIdx];
-
- auto tikInside = getDiffusionCoefficient(phaseIdx, compIIdx, compkIdx, problem, element, insideVolVars, fvGeometry.scv(insideScvIdx));
- tikInside = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), tikInside);
- auto tikOutside = getDiffusionCoefficient(phaseIdx, compIIdx, compkIdx, problem, element, outsideVolVars, fvGeometry.scv(outsideScvIdx));
- tikOutside = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), tikOutside);
-
- // scale by extrusion factor
- tikInside *= insideVolVars.extrusionFactor();
- tikOutside *= outsideVolVars.extrusionFactor();
-
- // the resulting averaged diffusion tensor
- const auto tik = problem.spatialParams().harmonicMean(tikInside, tikOutside, scvf.unitOuterNormal());
-
- reducedDiffusionMatrix[compIIdx][compIIdx] += xk/tik;
- }
-
- // now set the rest of the entries (off-diagonal)
- for (int compJIdx = 0; compJIdx < numComponents-1; compJIdx++)
- {
- //we don't want to calculate e.g. water in water diffusion
- if (compIIdx == compJIdx)
- continue;
- //calculate diffusivity for compIIdx, compJIdx
- auto tijInside = getDiffusionCoefficient(phaseIdx, compIIdx, compJIdx, problem, element, insideVolVars, fvGeometry.scv(insideScvIdx));
- tijInside = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), tijInside);
- auto tijOutside = getDiffusionCoefficient(phaseIdx, compIIdx, compJIdx, problem, element, outsideVolVars, fvGeometry.scv(outsideScvIdx));
- tijOutside = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), outsideVolVars.saturation(phaseIdx), tijOutside);
-
- // scale by extrusion factor
- tijInside *= insideVolVars.extrusionFactor();
- tijOutside *= outsideVolVars.extrusionFactor();
-
- // the resulting averaged diffusion tensor
- const auto tij = problem.spatialParams().harmonicMean(tijInside, tijOutside, scvf.unitOuterNormal());
-
- reducedDiffusionMatrix[compIIdx][compJIdx] +=xi*(1/tin - 1/tij);
- }
- }
+ const auto tij = problem.spatialParams().harmonicMean(tijInside, tijOutside, scvf.unitOuterNormal());
+ reducedDiffusionMatrix[compIIdx][compIIdx] += xj/tij;
+ reducedDiffusionMatrix[compIIdx][compJIdx] +=xi*(1/tin - 1/tij);
+ }
+ }
return reducedDiffusionMatrix;
}
diff --git a/dumux/discretization/cellcentered/tpfa/maxwellstefanslaw.hh b/dumux/discretization/cellcentered/tpfa/maxwellstefanslaw.hh
index cdf5a579f399ca76a66e7cbad21b3f153dacc087..ad226ea7baeaa86908cab08f1f95fc21ec12179a 100644
--- a/dumux/discretization/cellcentered/tpfa/maxwellstefanslaw.hh
+++ b/dumux/discretization/cellcentered/tpfa/maxwellstefanslaw.hh
@@ -120,8 +120,8 @@ public:
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideScv = fvGeometry.scv(insideScvIdx);
const Scalar omegai = calculateOmega_(scvf,
- insideScv,
- insideVolVars.extrusionFactor());
+ insideScv,
+ insideVolVars.extrusionFactor());
//now we have to do the tpfa: J_i = J_j which leads to: tij(xi -xj) = -rho Bi^-1 omegai(x*-xi) with x* = (omegai Bi^-1 + omegaj Bj^-1)^-1 (xi omegai Bi^-1 + xj omegaj Bj^-1) with i inside and j outside
reducedDiffusionMatrixInside = setupMSMatrix_(problem, element, fvGeometry, insideVolVars, insideScv, phaseIdx);
@@ -132,8 +132,6 @@ public:
moleFracOutside -= moleFracInside;
reducedDiffusionMatrixInside.solve(reducedFlux, moleFracOutside);
reducedFlux *= omegai;
-
-
}
else //we need outside cells as well if we are not on the boundary
{
@@ -146,12 +144,12 @@ public:
if (dim == dimWorld)
// assume the normal vector from outside is anti parallel so we save flipping a vector
omegaj = -1.0*calculateOmega_(scvf,
- outsideScv,
- outsideVolVars.extrusionFactor());
+ outsideScv,
+ outsideVolVars.extrusionFactor());
else
omegaj = calculateOmega_(fvGeometry.flipScvf(scvf.index()),
- outsideScv,
- outsideVolVars.extrusionFactor());
+ outsideScv,
+ outsideVolVars.extrusionFactor());
reducedDiffusionMatrixInside.invert();
reducedDiffusionMatrixOutside.invert();
@@ -204,11 +202,11 @@ private:
}
static ReducedComponentMatrix setupMSMatrix_(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const VolumeVariables& volVars,
- const SubControlVolume& scv,
- const int phaseIdx)
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const VolumeVariables& volVars,
+ const SubControlVolume& scv,
+ const int phaseIdx)
{
ReducedComponentMatrix reducedDiffusionMatrix(0.0);
@@ -219,33 +217,24 @@ private:
for (int compIIdx = 0; compIIdx < numComponents-1; compIIdx++)
{
const auto xi = volVars.moleFraction(phaseIdx, compIIdx);
-
- //calculate diffusivity for i,numComponents
Scalar tin = getDiffusionCoefficient(phaseIdx, compIIdx, numComponents-1, problem, element, volVars, scv);
tin = EffDiffModel::effectiveDiffusivity(volVars.porosity(), volVars.saturation(phaseIdx), tin);
- //set the entrys of the diffusion matrix of the diagonal
- reducedDiffusionMatrix[compIIdx][compIIdx] += xi/tin;
- for (int compkIdx = 0; compkIdx < numComponents; compkIdx++)
- {
- if (compkIdx == compIIdx)
- continue;
- const auto xk = volVars.moleFraction(phaseIdx, compkIdx);
- Scalar tik = getDiffusionCoefficient(phaseIdx, compIIdx, compkIdx, problem, element, volVars, scv);
- tik = EffDiffModel::effectiveDiffusivity(volVars.porosity(), volVars.saturation(phaseIdx), tik);
- reducedDiffusionMatrix[compIIdx][compIIdx] += xk/tik;
- }
+ // set the entries of the diffusion matrix of the diagonal
+ reducedDiffusionMatrix[compIIdx][compIIdx] += xi/tin;
- // now set the rest of the entries (off-diagonal)
- for (int compJIdx = 0; compJIdx < numComponents-1; compJIdx++)
+ // now set the rest of the entries (off-diagonal and additional entries for diagonal)
+ for (int compJIdx = 0; compJIdx < numComponents; compJIdx++)
{
- //we don't want to calculate e.g. water in water diffusion
- if (compIIdx == compJIdx)
+ // we don't want to calculate e.g. water in water diffusion
+ if (compJIdx == compIIdx)
continue;
- //calculate diffusivity for compIIdx, compJIdx
+
+ const auto xj = volVars.moleFraction(phaseIdx, compJIdx);
Scalar tij = getDiffusionCoefficient(phaseIdx, compIIdx, compJIdx, problem, element, volVars, scv);
tij = EffDiffModel::effectiveDiffusivity(volVars.porosity(), volVars.saturation(phaseIdx), tij);
- reducedDiffusionMatrix[compIIdx][compJIdx] +=xi*(1/tin - 1/tij);
+ reducedDiffusionMatrix[compIIdx][compIIdx] += xj/tij;
+ reducedDiffusionMatrix[compIIdx][compJIdx] += xi*(1/tin - 1/tij);
}
}
return reducedDiffusionMatrix;