Skip to content
GitLab
Commit 85ecc5ec authored by Kilian Weishaupt's avatar Kilian Weishaupt Committed by Ned Coltman
Browse files

[stokesdarcy][couplingdata] Fix index bug in FicksLaw

parent 4bb084ed
Hide whitespace changes
Inline Side-by-side
dumux/flux/box/fourierslawnonequilibrium.hh
View file @ 85ecc5ec
......@@ -59,7 +59,9 @@ class FouriersLawNonEquilibriumImplementation<TypeTag, DiscretizationMethod::box
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using Element = typename GridView::template Codim<0>::Entity;
static constexpr auto numEnergyEqSolid = getPropValue<TypeTag, Properties::NumEnergyEqSolid>();
static constexpr auto numEnergyEqFluid = getPropValue<TypeTag, Properties::NumEnergyEqFluid>();
static constexpr auto numEnergyEq = numEnergyEqSolid + numEnergyEqFluid;
static constexpr auto sPhaseIdx = ModelTraits::numFluidPhases();
public:
......@@ -76,29 +78,19 @@ public:
const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
const auto& insideVolVars = elemVolVars[insideScv];
const auto& outsideVolVars = elemVolVars[outsideScv];
Scalar insideLambda = 0.0;
Scalar outsideLambda = 0.0;
// effective diffusion tensors
if (phaseIdx != sPhaseIdx)
{
//when number of energyEq for the fluid are smaller than numPhases that means that we need an effecitve law
if (numEnergyEqFluid < ModelTraits::numFluidPhases())
{
insideLambda += insideVolVars.effectiveThermalConductivity();
outsideLambda += outsideVolVars.effectiveThermalConductivity();
}
const auto computeLambda = [&](const auto& v){
if constexpr (numEnergyEq == 1)
return v.effectiveThermalConductivity();
else if constexpr (numEnergyEqFluid == 1)
return (phaseIdx != sPhaseIdx)
? v.effectiveFluidThermalConductivity()
: v.effectiveSolidThermalConductivity();
else
{
insideLambda += insideVolVars.fluidThermalConductivity(phaseIdx)*insideVolVars.saturation(phaseIdx)*insideVolVars.porosity();
outsideLambda += outsideVolVars.fluidThermalConductivity(phaseIdx)*outsideVolVars.saturation(phaseIdx)*outsideVolVars.porosity();
}
}
//solid phase
else
{
insideLambda += insideVolVars.solidThermalConductivity()*(1.0-insideVolVars.porosity());
outsideLambda += outsideVolVars.solidThermalConductivity()*(1.0-outsideVolVars.porosity());
}
return v.effectivePhaseThermalConductivity(phaseIdx);
};
auto insideLambda = computeLambda(insideVolVars);
auto outsideLambda = computeLambda(outsideVolVars);
// scale by extrusion factor
insideLambda *= insideVolVars.extrusionFactor();
......
dumux/flux/cctpfa/fourierslawnonequilibrium.hh
View file @ 85ecc5ec
......@@ -59,7 +59,9 @@ class FouriersLawNonEquilibriumImplementation<TypeTag, DiscretizationMethod::cct
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
static constexpr auto numEnergyEqSolid = getPropValue<TypeTag, Properties::NumEnergyEqSolid>();
static constexpr auto numEnergyEqFluid = getPropValue<TypeTag, Properties::NumEnergyEqFluid>();
static constexpr auto numEnergyEq = numEnergyEqSolid + numEnergyEqFluid;
static constexpr auto sPhaseIdx = ModelTraits::numFluidPhases();
public:
......@@ -103,48 +105,32 @@ public:
const SubControlVolumeFace& scvf,
const int phaseIdx)
{
Scalar tij;
const auto insideScvIdx = scvf.insideScvIdx();
const auto& insideScv = fvGeometry.scv(insideScvIdx);
const auto& insideVolVars = elemVolVars[insideScvIdx];
Scalar insideLambda = 0.0;
Scalar outsideLambda = 0.0;
// effective diffusion tensors
if (phaseIdx != sPhaseIdx)
{
//when number of energyEq for the fluid are smaller than numPhases that means that we need an effecitve law
if (numEnergyEqFluid < ModelTraits::numFluidPhases())
insideLambda += insideVolVars.effectiveThermalConductivity();
else // numEnergyEqFluid >1
insideLambda += insideVolVars.fluidThermalConductivity(phaseIdx)*insideVolVars.saturation(phaseIdx)*insideVolVars.porosity();
}
else // solid phase
insideLambda += insideVolVars.solidThermalConductivity()*(1.0-insideVolVars.porosity());
const auto computeLambda = [&](const auto& v){
if constexpr (numEnergyEq == 1)
return v.effectiveThermalConductivity();
else if constexpr (numEnergyEqFluid == 1)
return (phaseIdx != sPhaseIdx)
? v.effectiveFluidThermalConductivity()
: v.effectiveSolidThermalConductivity();
else
return v.effectivePhaseThermalConductivity(phaseIdx);
};
const auto insideLambda = computeLambda(insideVolVars);
const Scalar ti = computeTpfaTransmissibility(scvf, insideScv, insideLambda, insideVolVars.extrusionFactor());
// for the boundary (dirichlet) or at branching points we only need ti
if (scvf.boundary() || scvf.numOutsideScvs() > 1)
tij = scvf.area()*ti;
return scvf.area()*ti;
else // otherwise we compute a tpfa harmonic mean
{
const auto outsideScvIdx = scvf.outsideScvIdx();
const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
const auto& outsideVolVars = elemVolVars[outsideScvIdx];
// effective diffusion tensors
if (phaseIdx != sPhaseIdx)
{
//when number of energyEq for the fluid are smaller than numPhases that means that we need an effecitve law
if (numEnergyEqFluid < ModelTraits::numFluidPhases())
outsideLambda += outsideVolVars.effectiveThermalConductivity();
else
outsideLambda += outsideVolVars.fluidThermalConductivity(phaseIdx)*outsideVolVars.saturation(phaseIdx)*outsideVolVars.porosity();
}
else //solid phase
outsideLambda +=outsideVolVars.solidThermalConductivity()*(1.0-outsideVolVars.porosity());
const auto outsideLambda = computeLambda(outsideVolVars);
Scalar tj;
if (dim == dimWorld)
......@@ -155,11 +141,10 @@ public:
// check for division by zero!
if (ti*tj <= 0.0)
tij = 0;
return 0.0;
else
tij = scvf.area()*(ti * tj)/(ti + tj);
return scvf.area()*(ti * tj)/(ti + tj);
}
return tij;
}
};
......
dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
View file @ 85ecc5ec
......@@ -1063,6 +1063,8 @@ protected:
for (int compIdx = 1; compIdx < numComponents; ++compIdx)
{
const int domainIMainCompIdx = couplingPhaseIdx(domainI);
const int domainJMainCompIdx = couplingPhaseIdx(domainJ);
const int domainICompIdx = couplingCompIdx(domainI, compIdx);
const int domainJCompIdx = couplingCompIdx(domainJ, compIdx);
......@@ -1076,8 +1078,8 @@ protected:
domainJ,
insideDistance,
outsideDistance,
Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(volVarsI, couplingPhaseIdx(domainI), domainICompIdx, domainICompIdx),
Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(volVarsJ, couplingPhaseIdx(domainJ), domainICompIdx, domainJCompIdx),
Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(volVarsI, couplingPhaseIdx(domainI), domainIMainCompIdx, domainICompIdx),
Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(volVarsJ, couplingPhaseIdx(domainJ), domainJMainCompIdx, domainJCompIdx),
diffCoeffAvgType);
diffusiveFlux[domainICompIdx] += -avgDensity * tij * deltaMassOrMoleFrac;
}
......
dumux/porousmediumflow/nonisothermal/volumevariables.hh
View file @ 85ecc5ec
......@@ -149,8 +149,11 @@ class EnergyVolumeVariablesImplementation<Traits, Impl, true>
using ParentType = PorousMediumFlowVolumeVariables<Traits>;
using EffCondModel = typename Traits::EffectiveThermalConductivityModel;
static const int temperatureIdx = Idx::temperatureIdx;
static const int numEnergyEq = Traits::ModelTraits::numEnergyEq();
static constexpr int temperatureIdx = Idx::temperatureIdx;
static constexpr int numEnergyEq = Traits::ModelTraits::numEnergyEq();
static constexpr bool fullThermalEquilibrium = (numEnergyEq == 1);
static constexpr bool fluidThermalEquilibrium = (numEnergyEq == 2);
public:
// export the fluidstate
......@@ -169,7 +172,7 @@ public:
FluidState& fluidState,
SolidState& solidState)
{
if (numEnergyEq == 1)
if constexpr (fullThermalEquilibrium)
{
// retrieve temperature from solution vector, all phases have the same temperature
const Scalar T = elemSol[scv.localDofIndex()][temperatureIdx];
......@@ -182,8 +185,8 @@ public:
else
{
// if numEnergyEq == 2 this means we have 1 temp for fluid phase, one for solid
if (numEnergyEq == 2)
// this means we have 1 temp for fluid phase, one for solid
if constexpr (fluidThermalEquilibrium)
{
const Scalar T = elemSol[scv.localDofIndex()][temperatureIdx];
for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
......@@ -226,7 +229,28 @@ public:
// updates the effective thermal conductivity
void updateEffectiveThermalConductivity()
{
lambdaEff_ = EffCondModel::effectiveThermalConductivity(asImp_());
if constexpr (fullThermalEquilibrium)
{
// Full Thermal Equilibirum: One effectiveThermalConductivity value for all phases (solid & fluid).
lambdaEff_[0] = EffCondModel::effectiveThermalConductivity(asImp_());
}
else if constexpr (fluidThermalEquilibrium)
{
// Fluid Thermal Equilibrium (Partial Nonequilibrium): One effectiveThermalConductivity for the fluids, one for the solids.
Scalar fluidLambda = 0.0;
for (int phaseIdx = 0; phaseIdx < FluidSystem::numPhases; phaseIdx++)
fluidLambda += fluidThermalConductivity(phaseIdx) * asImp_().saturation(phaseIdx) * asImp_().porosity();
lambdaEff_[0] = fluidLambda;
lambdaEff_[numEnergyEq-1] = solidThermalConductivity() * (1.0 - asImp_().porosity());
}
else
{
// Full Thermal Nonequilibrium: One effectiveThermal Conductivity per phase (solid & fluid).
for (int phaseIdx = 0; phaseIdx < FluidSystem::numPhases; phaseIdx++)
lambdaEff_[phaseIdx] = fluidThermalConductivity(phaseIdx) * asImp_().saturation(phaseIdx) * asImp_().porosity();
lambdaEff_[numEnergyEq-1] = solidThermalConductivity() * (1.0 - asImp_().porosity());
}
}
/*!
......@@ -278,24 +302,56 @@ public:
{ return asImp_().solidState().density(); }
/*!
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$ of the solid phase in the sub-control volume.
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
* of the solid phase in the sub-control volume.
*/
Scalar solidThermalConductivity() const
{ return asImp_().solidState().thermalConductivity(); }
/*!
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$ of a fluid phase in
* the sub-control volume.
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
* of a fluid phase in the sub-control volume.
*/
Scalar fluidThermalConductivity(const int phaseIdx) const
{ return FluidSystem::thermalConductivity(asImp_().fluidState(), phaseIdx); }
/*!
* \brief Returns the effective thermal conductivity \f$\mathrm{[W/(m*K)]}\f$ in
* the sub-control volume.
* the sub-control volume. Specific to equilibirum models (case fullThermalEquilibrium).
*/
template< bool enable = fullThermalEquilibrium,
std::enable_if_t<enable, int> = 0>
Scalar effectiveThermalConductivity() const
{ return lambdaEff_; }
{ return lambdaEff_[0]; }
/*!
* \brief Returns the effective thermal conductivity \f$\mathrm{[W/(m*K)]}\f$ of the fluids in
* the sub-control volume. Specific to partially nonequilibrium models (case fluidThermalEquilibrium).
*/
template< bool enable = fluidThermalEquilibrium,
std::enable_if_t<enable, int> = 0>
Scalar effectiveFluidThermalConductivity() const
{ return lambdaEff_[0]; }
/*!
* \brief Returns the effective thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
* of the solid phase in the sub-control volume.
* Specific to partially nonequilibrium models (case fluidThermalEquilibrium)
*/
template< bool enable = fluidThermalEquilibrium,
std::enable_if_t<enable, int> = 0>
Scalar effectiveSolidThermalConductivity() const
{ return lambdaEff_[numEnergyEq-1]; }
/*!
* \brief Returns the effective thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
* per fluid phase in the sub-control volume.
* Specific to nonequilibrium models (case full non-equilibrium)
*/
template< bool enable = (!fullThermalEquilibrium && !fluidThermalEquilibrium),
std::enable_if_t<enable, int> = 0>
Scalar effectivePhaseThermalConductivity(const int phaseIdx) const
{ return lambdaEff_[phaseIdx]; }
//! The phase enthalpy is zero for isothermal models
//! This is needed for completing the fluid state
......@@ -473,7 +529,7 @@ private:
return problem.spatialParams().solidThermalConductivity(element, scv, elemSol, solidState);
}
Scalar lambdaEff_;
std::array<Scalar, numEnergyEq> lambdaEff_;
// \}
};
......
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