diff --git a/dumux/flux/box/fickslaw.hh b/dumux/flux/box/fickslaw.hh
index acd625fa0bdd1be73228ff70eb2a4a3431b508f6..8bc2cfa8aa8a0ff092d3ecf741a312e0c15d3e60 100644
--- a/dumux/flux/box/fickslaw.hh
+++ b/dumux/flux/box/fickslaw.hh
@@ -76,7 +76,10 @@ class FicksLawImplementation<TypeTag, DiscretizationMethod::box, referenceSystem
public:
template<int numPhases, int numComponents>
- using DiffusionCoefficientsContainer = FickianDiffusionCoefficients<Scalar, numPhases, numComponents>;
+ using DiffusionCoefficientsContainer = FickianDiffusionCoefficients<Scalar,
+ numPhases,
+ numComponents,
+ FluidSystem::isTracerFluidSystem()>;
//return the reference system
static constexpr ReferenceSystemFormulation referenceSystemFormulation()
diff --git a/dumux/flux/ccmpfa/fickslaw.hh b/dumux/flux/ccmpfa/fickslaw.hh
index 1b56a8a75d73fe86f228a958880be004c0e520f0..e4c69a4a7ee6df973235f9da319942e87d0aef0f 100644
--- a/dumux/flux/ccmpfa/fickslaw.hh
+++ b/dumux/flux/ccmpfa/fickslaw.hh
@@ -170,7 +170,10 @@ public:
using Cache = MpfaFicksLawCache;
//export the diffusion container
template<int numPhases, int numComponents>
- using DiffusionCoefficientsContainer = FickianDiffusionCoefficients<Scalar, numPhases, numComponents>;
+ using DiffusionCoefficientsContainer = FickianDiffusionCoefficients<Scalar,
+ numPhases,
+ numComponents,
+ FluidSystem::isTracerFluidSystem()>;
//! Compute the diffusive flux across an scvf
static ComponentFluxVector flux(const Problem& problem,
diff --git a/dumux/flux/cctpfa/fickslaw.hh b/dumux/flux/cctpfa/fickslaw.hh
index 0f94756dca8472a32d9f7039bc8036bfcab99092..b41c2ed11b345920d8094daf96ada305574ead96 100644
--- a/dumux/flux/cctpfa/fickslaw.hh
+++ b/dumux/flux/cctpfa/fickslaw.hh
@@ -128,7 +128,10 @@ public:
using Cache = TpfaFicksLawCache;
template<int numPhases, int numComponents>
- using DiffusionCoefficientsContainer = FickianDiffusionCoefficients<Scalar, numPhases, numComponents>;
+ using DiffusionCoefficientsContainer = FickianDiffusionCoefficients<Scalar,
+ numPhases,
+ numComponents,
+ FluidSystem::isTracerFluidSystem()>;
//! return diffusive fluxes for all components in a phase
static ComponentFluxVector flux(const Problem& problem,
diff --git a/dumux/flux/fickiandiffusioncoefficients.hh b/dumux/flux/fickiandiffusioncoefficients.hh
index f95429ce0ffa9c6f80ad5acf3f8a6b9201cacb15..cba9c93bda904c1fb98c9f7657ff01389a567d04 100644
--- a/dumux/flux/fickiandiffusioncoefficients.hh
+++ b/dumux/flux/fickiandiffusioncoefficients.hh
@@ -29,9 +29,27 @@
namespace Dumux {
-// General Case (mpnc)
+/*!
+ * \ingroup Flux
+ * \brief Container storing the diffusion coefficients required by Fick's law.
+ * Uses the minimal possible container size and provides unified access.
+ * \tparam Scalar The type used for scalar values
+ * \tparam numPhases Number of phases in the fluid composition
+ * \tparam numComponents Number of components in the fluid composition
+ * \tparam allTracerComponents If false, this means that the main component of
+ * a phase is part of the components. In this case,
+ * the storage container is optimized with respect to
+ * memory consumption as diffusion coefficients of the
+ * main component of a phase in itself are not stored.
+ * If true, all diffusion coefficients of all components
+ * are stored
+ */
+template <class Scalar, int numPhases, int numComponents, bool allTracerComponents>
+class FickianDiffusionCoefficients;
+
+//! General case (mpnc), for compositions containing the phases' main components
template <class Scalar, int numPhases, int numComponents>
-class FickianDiffusionCoefficients
+class FickianDiffusionCoefficients<Scalar, numPhases, numComponents, false>
{
public:
template<class DiffCoeffFunc>
@@ -63,9 +81,9 @@ private:
}
};
-// Specialization for 1pnc
+//! Specialization for 1pnc & compositions containing the phases' main components
template <class Scalar, int numComponents>
-class FickianDiffusionCoefficients<Scalar, 1, numComponents>
+class FickianDiffusionCoefficients<Scalar, 1, numComponents, false>
{
public:
template<class DiffCoeffFunc>
@@ -94,9 +112,9 @@ private:
{ return compJIdx-1; }
};
-// Specialization for 2p2c
+//! Specialization for 2p2c & compositions containing the phases' main components
template <class Scalar>
-class FickianDiffusionCoefficients<Scalar, 2, 2>
+class FickianDiffusionCoefficients<Scalar, 2, 2, false>
{
public:
template<class DiffCoeffFunc>
@@ -118,9 +136,9 @@ private:
{ return phaseIdx; }
};
-// Specialization for 3p2c
+//! Specialization for 3p2c & compositions containing the phases' main components
template <class Scalar>
-class FickianDiffusionCoefficients<Scalar, 3, 2>
+class FickianDiffusionCoefficients<Scalar, 3, 2, false>
{
public:
template<class DiffCoeffFunc>
@@ -143,6 +161,36 @@ private:
{ return phaseIdx; }
};
+//! Specialization for mpnc & compositions that only contain tracers
+template <class Scalar, int numPhases, int numComponents>
+class FickianDiffusionCoefficients<Scalar, numPhases, numComponents, true>
+{
+public:
+ template<class DiffCoeffFunc>
+ void update(DiffCoeffFunc& computeDiffCoeff)
+ {
+ for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx)
+ diffCoeff_[getIndex_(phaseIdx, compIdx)]
+ = computeDiffCoeff(phaseIdx, phaseIdx, compIdx);
+ }
+
+ const Scalar& operator()(int phaseIdx, int compIIdx, int compJIdx) const
+ {
+ assert(phaseIdx < numPhases);
+ assert(compIIdx != compJIdx);
+ assert(compIIdx < numComponents);
+ assert(compJIdx < numComponents);
+ assert(compIIdx < compJIdx);
+ return diffCoeff_[getIndex_(phaseIdx, compJIdx)];
+ }
+
+private:
+ std::array<Scalar, numPhases*numComponents> diffCoeff_;
+ int getIndex_(int phaseIdx, int compJIdx) const
+ { return phaseIdx*numComponents + compJIdx; }
+};
+
} // end namespace Dumux
#endif
diff --git a/dumux/flux/maxwellstefandiffusioncoefficients.hh b/dumux/flux/maxwellstefandiffusioncoefficients.hh
index 8e7ffe4b090756321f6fe275d1c2f43031d1ea6a..6d9a8d08ad9187178647b8bf8f0bdd1d3228d790 100644
--- a/dumux/flux/maxwellstefandiffusioncoefficients.hh
+++ b/dumux/flux/maxwellstefandiffusioncoefficients.hh
@@ -29,7 +29,15 @@
namespace Dumux {
-// General case
+/*!
+ * \ingroup Flux
+ * \brief Container storing the diffusion coefficients required by the Maxwell-
+ * Stefan diffusion law. Uses the minimal possible container size and
+ * provides unified access.
+ * \tparam Scalar The type used for scalar values
+ * \tparam numPhases Number of phases in the fluid composition
+ * \tparam numComponents Number of components in the fluid composition
+ */
template <class Scalar, int numPhases, int numComponents>
class MaxwellStefanDiffusionCoefficients
{
@@ -41,7 +49,8 @@ public:
for (unsigned int compIIdx = 0; compIIdx < numComponents; ++compIIdx)
for (unsigned int compJIdx = 0; compJIdx < numComponents; ++compJIdx)
if(compIIdx != compJIdx && compIIdx < compJIdx)
- diffCoeff_[getIndex_(phaseIdx, compIIdx, compJIdx)] = computeDiffCoeff(phaseIdx, compIIdx, compJIdx);
+ diffCoeff_[getIndex_(phaseIdx, compIIdx, compJIdx)]
+ = computeDiffCoeff(phaseIdx, compIIdx, compJIdx);
}
const Scalar& operator()(int phaseIdx, int compIIdx, int compJIdx) const
@@ -54,10 +63,10 @@ private:
std::array<Scalar, (numPhases * ((numComponents * numComponents - numComponents)/2))> diffCoeff_;
const int getIndex_(int phaseIdx, int compIIdx, int compJIdx) const
{
- return phaseIdx * ((numComponents * numComponents - numComponents) / 2)
- + compIIdx * numComponents
- - ((compIIdx * compIIdx + compIIdx) / 2)
- + compJIdx - (compIIdx +1);
+ return phaseIdx * ((numComponents * numComponents - numComponents) / 2)
+ + compIIdx * numComponents
+ - ((compIIdx * compIIdx + compIIdx) / 2)
+ + compJIdx - (compIIdx +1);
}
};
diff --git a/dumux/flux/staggered/freeflow/fickslaw.hh b/dumux/flux/staggered/freeflow/fickslaw.hh
index 6c0c426f6b9532200dada8b1d18941b871ab9ac2..e204c20075162ba10dc3adfae8340040aef27d70 100644
--- a/dumux/flux/staggered/freeflow/fickslaw.hh
+++ b/dumux/flux/staggered/freeflow/fickslaw.hh
@@ -79,7 +79,10 @@ public:
using Cache = FluxVariablesCaching::EmptyDiffusionCache;
template<int numPhases, int numComponents>
- using DiffusionCoefficientsContainer = FickianDiffusionCoefficients<Scalar, numPhases, numComponents>;
+ using DiffusionCoefficientsContainer = FickianDiffusionCoefficients<Scalar,
+ numPhases,
+ numComponents,
+ FluidSystem::isTracerFluidSystem()>;
template<class Problem, class ElementVolumeVariables>
static NumEqVector flux(const Problem& problem,
diff --git a/dumux/porousmediumflow/tracer/volumevariables.hh b/dumux/porousmediumflow/tracer/volumevariables.hh
index c90aeae1c876b4a7edfad99107d20e2c57d0ec33..fbd1ec6fbaf504ca95e5a18955008abcfe710736 100644
--- a/dumux/porousmediumflow/tracer/volumevariables.hh
+++ b/dumux/porousmediumflow/tracer/volumevariables.hh
@@ -59,7 +59,7 @@ class TracerVolumeVariables
static constexpr bool useMoles = Traits::ModelTraits::useMoles();
using EffDiffModel = typename Traits::EffectiveDiffusivityModel;
static constexpr int numFluidComps = ParentType::numFluidComponents();
- using DiffusionCoefficients = typename Traits::DiffusionType::template DiffusionCoefficientsContainer<1, numFluidComps+1>;
+ using DiffusionCoefficients = typename Traits::DiffusionType::template DiffusionCoefficientsContainer<1, numFluidComps>;
public:
//! Export fluid system type
@@ -98,13 +98,10 @@ public:
}
// Update the binary diffusion and effective diffusion coefficients.
- // The container starts counting from compJIdx = 1 as it assumes the
- // phase (main component) to have compIdx = 0. Tracer fluid systems
- // start counting from zero with the actual tracer components, so we
- // have to subtract 1 here before calling the fluid system
auto getDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
{
- return FluidSystem::binaryDiffusionCoefficient( compJIdx-1,
+ return FluidSystem::binaryDiffusionCoefficient( compIIdx,
+ compJIdx,
problem,
element,
scv);
@@ -112,7 +109,7 @@ public:
auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
{
- return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx-1, compJIdx-1);
+ return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
};
diffCoeff_.update(getDiffusionCoefficient);
@@ -207,27 +204,19 @@ public:
*/
[[deprecated("Signature deprecated. Use diffusionCoefficient(phaseIdx, compIIdx, compJIdx)!")]]
Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
- { return diffCoeff_(0, 0, compIdx+1); }
+ { return diffCoeff_(0, 0, compIdx); }
/*!
* \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
- * \note The container starts counting from compJIdx = 1 as it assumes the
- * phase (main component) to have compIdx = 0. Tracer fluid systems
- * start counting from zero with the actual tracer components, so we
- * have to increase compJIdx by 1 before calling the container.
*/
Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
- { return diffCoeff_(phaseIdx, compIIdx+1, compJIdx+1); }
+ { return diffCoeff_(phaseIdx, compIIdx, compJIdx); }
/*!
* \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$.
- * \note The container starts counting from compJIdx = 1 as it assumes the
- * phase (main component) to have compIdx = 0. Tracer fluid systems
- * start counting from zero with the actual tracer components, so we
- * have to increase compJIdx by 1 before calling the container.
*/
Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
- { return effectiveDiffCoeff_(phaseIdx, compIIdx+1, compJIdx+1); }
+ { return effectiveDiffCoeff_(phaseIdx, compIIdx, compJIdx); }
// /*!
// * \brief Returns the dispersivity of the fluid's streamlines.
diff --git a/test/porousmediumflow/tracer/2ptracer/problem_tracer.hh b/test/porousmediumflow/tracer/2ptracer/problem_tracer.hh
index eade9d78214aa32bb7be47ea628b2af2d5be2602..9bfa216b235964c158b1fdc14cc3a15f231b4890 100644
--- a/test/porousmediumflow/tracer/2ptracer/problem_tracer.hh
+++ b/test/porousmediumflow/tracer/2ptracer/problem_tracer.hh
@@ -105,7 +105,8 @@ public:
//! Binary diffusion coefficient
//! (might depend on spatial parameters like pressure / temperature)
- static Scalar binaryDiffusionCoefficient(unsigned int compIdx,
+ static Scalar binaryDiffusionCoefficient(unsigned int compIIdx,
+ unsigned int compJIdx,
const Problem& problem,
const Element& element,
const SubControlVolume& scv)
diff --git a/test/porousmediumflow/tracer/constvel/problem.hh b/test/porousmediumflow/tracer/constvel/problem.hh
index becabaa3da649c100e7efa1e00c17b6ee6e1ea2c..63e5bf5de0e1dd5e3d9c0d1735cee7b59e16a504 100644
--- a/test/porousmediumflow/tracer/constvel/problem.hh
+++ b/test/porousmediumflow/tracer/constvel/problem.hh
@@ -126,14 +126,15 @@ public:
//! Binary diffusion coefficient
//! (might depend on spatial parameters like pressure / temperature)
- static Scalar binaryDiffusionCoefficient(unsigned int compIdx,
+ static Scalar binaryDiffusionCoefficient(unsigned int compIIdx,
+ unsigned int compJIdx,
const Problem& problem,
const Element& element,
const SubControlVolume& scv)
{
static const Scalar D = getParam<Scalar>("Problem.D");
static const Scalar D2 = getParam<Scalar>("Problem.D2");
- if (compIdx == 0)
+ if (compJIdx == 0)
return D;
else
return D2;