Commit 879eeac2 authored by Dennis Gläser's avatar Dennis Gläser Committed by Ned Coltman
Browse files

[tracer][diffusion] make tracer model support both fickian and ms diffusion

parent 877db73e
......@@ -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()
......
......@@ -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,
......
......@@ -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,
......
......@@ -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
......@@ -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);
}
};
......
......@@ -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,
......
......@@ -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.
......
......@@ -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)
......
......@@ -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;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment