diff --git a/dumux/flux/box/darcyslaw.hh b/dumux/flux/box/darcyslaw.hh index a1326195e31839e3725ce92242ced0086b24ad44..5662119778b735f70ef23d23ab0c9740cafde970 100644 --- a/dumux/flux/box/darcyslaw.hh +++ b/dumux/flux/box/darcyslaw.hh @@ -33,6 +33,7 @@ #include #include #include +#include namespace Dumux { @@ -105,7 +106,7 @@ public: insideK *= insideVolVars.extrusionFactor(); outsideK *= outsideVolVars.extrusionFactor(); - const auto K = problem.spatialParams().harmonicMean(insideK, outsideK, scvf.unitOuterNormal()); + const auto K = faceTensorAverage(insideK, outsideK, scvf.unitOuterNormal()); static const bool enableGravity = getParamFromGroup(problem.paramGroup(), "Problem.EnableGravity"); const auto& shapeValues = fluxVarCache.shapeValues(); @@ -152,7 +153,7 @@ public: insideK *= insideVolVars.extrusionFactor(); outsideK *= outsideVolVars.extrusionFactor(); - const auto K = problem.spatialParams().harmonicMean(insideK, outsideK, scvf.unitOuterNormal()); + const auto K = faceTensorAverage(insideK, outsideK, scvf.unitOuterNormal()); std::vector ti(fvGeometry.numScv()); for (const auto& scv : scvs(fvGeometry)) diff --git a/dumux/flux/box/fickslaw.hh b/dumux/flux/box/fickslaw.hh index 0349e2549cab98b22f64b3412b803f72fc6b5f04..f802f6ec7f9afeaab4475c825aa504f8fa12ef02 100644 --- a/dumux/flux/box/fickslaw.hh +++ b/dumux/flux/box/fickslaw.hh @@ -35,6 +35,7 @@ #include #include +#include namespace Dumux { @@ -180,7 +181,7 @@ private: outsideD *= outsideVV.extrusionFactor(); // the resulting averaged diffusion tensor - return problem.spatialParams().harmonicMean(insideD, outsideD, scvf.unitOuterNormal()); + return faceTensorAverage(insideD, outsideD, scvf.unitOuterNormal()); } static std::pair diff --git a/dumux/flux/box/forchheimerslaw.hh b/dumux/flux/box/forchheimerslaw.hh index 0ce29b5f0e92c4d60659a4668d6c9e27cb699c0a..3ea0759f6d8d300053066622d1e728319597fbe4 100644 --- a/dumux/flux/box/forchheimerslaw.hh +++ b/dumux/flux/box/forchheimerslaw.hh @@ -35,6 +35,7 @@ #include #include #include +#include namespace Dumux { @@ -121,7 +122,7 @@ public: insideK *= insideVolVars.extrusionFactor(); outsideK *= outsideVolVars.extrusionFactor(); - const auto K = problem.spatialParams().harmonicMean(insideK, outsideK, scvf.unitOuterNormal()); + const auto K = faceTensorAverage(insideK, outsideK, scvf.unitOuterNormal()); static const bool enableGravity = getParamFromGroup(problem.paramGroup(), "Problem.EnableGravity"); const auto& shapeValues = fluxVarCache.shapeValues(); diff --git a/dumux/flux/box/fourierslaw.hh b/dumux/flux/box/fourierslaw.hh index 6a825e69464a0a4ec4a5fb42a94c06651202bab1..0a556891128ff337583fea50016ee8a022332c95 100644 --- a/dumux/flux/box/fourierslaw.hh +++ b/dumux/flux/box/fourierslaw.hh @@ -29,6 +29,7 @@ #include #include #include +#include namespace Dumux { @@ -84,7 +85,7 @@ public: outsideLambda *= outsideVolVars.extrusionFactor(); // the resulting averaged diffusion tensor - const auto lambda = problem.spatialParams().harmonicMean(insideLambda, outsideLambda, scvf.unitOuterNormal()); + const auto lambda = faceTensorAverage(insideLambda, outsideLambda, scvf.unitOuterNormal()); // evaluate gradTemp at integration point const auto& fluxVarsCache = elemFluxVarsCache[scvf]; diff --git a/dumux/flux/box/fourierslawnonequilibrium.hh b/dumux/flux/box/fourierslawnonequilibrium.hh index 22ea175d80443dd893be1e4722cb37f87ba2f2db..c6b33d9597346fea21faf33055ba024353a5929b 100644 --- a/dumux/flux/box/fourierslawnonequilibrium.hh +++ b/dumux/flux/box/fourierslawnonequilibrium.hh @@ -33,6 +33,8 @@ #include #include +#include + namespace Dumux { // forward declaration @@ -101,7 +103,7 @@ public: outsideLambda *= outsideVolVars.extrusionFactor(); // the resulting averaged diffusion tensor - const auto lambda = problem.spatialParams().harmonicMean(insideLambda, outsideLambda, scvf.unitOuterNormal()); + const auto lambda = faceTensorAverage(insideLambda, outsideLambda, scvf.unitOuterNormal()); // evaluate gradTemp at integration point const auto& fluxVarsCache = elemFluxVarsCache[scvf]; diff --git a/dumux/flux/box/maxwellstefanslaw.hh b/dumux/flux/box/maxwellstefanslaw.hh index fa2a845e9a7109e87200e7986ef7230c177cba63..758e80e18cc6710c8fd727b100c1d7ef84c0639a 100644 --- a/dumux/flux/box/maxwellstefanslaw.hh +++ b/dumux/flux/box/maxwellstefanslaw.hh @@ -36,6 +36,7 @@ #include #include #include +#include namespace Dumux { @@ -181,7 +182,7 @@ private: tinOutside *= outsideVolVars.extrusionFactor(); // the resulting averaged diffusion tensor - const auto tin = problem.spatialParams().harmonicMean(tinInside, tinOutside, scvf.unitOuterNormal()); + const auto tin = faceTensorAverage(tinInside, tinOutside, scvf.unitOuterNormal()); //begin the entrys of the diffusion matrix of the diagonal reducedDiffusionMatrix[compIIdx][compIIdx] += xi*avgMolarMass/(tin*Mn); @@ -207,7 +208,7 @@ private: tijOutside *= outsideVolVars.extrusionFactor(); // the resulting averaged diffusion tensor - const auto tij = problem.spatialParams().harmonicMean(tijInside, tijOutside, scvf.unitOuterNormal()); + const auto tij = faceTensorAverage(tijInside, tijOutside, scvf.unitOuterNormal()); reducedDiffusionMatrix[compIIdx][compIIdx] += xj*avgMolarMass/(tij*Mi); if (compJIdx < numComponents-1) diff --git a/dumux/flux/facetensoraverage.hh b/dumux/flux/facetensoraverage.hh new file mode 100644 index 0000000000000000000000000000000000000000..38a1f5b9f838714dcf2cc106bd0ba19b7659ec77 --- /dev/null +++ b/dumux/flux/facetensoraverage.hh @@ -0,0 +1,87 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \ingroup Flux + * \brief A free function to average a Tensor at an interface. + */ +#ifndef DUMUX_FACE_TENSOR_AVERAGE_HH +#define DUMUX_FACE_TENSOR_AVERAGE_HH + +#include +#include + +namespace Dumux { + +/*! + * \brief Average of a discontinuous scalar field at discontinuity interface + * (for compatibility reasons with the function below) + * \return the harmonic average of the scalars + * \param T1 first scalar parameter + * \param T2 second scalar parameter + * \param normal The unit normal vector of the interface + */ +template +Scalar faceTensorAverage(const Scalar T1, + const Scalar T2, + const Dune::FieldVector& normal) +{ return Dumux::harmonicMean(T1, T2); } + +/*! + * \brief Average of a discontinuous tensorial field at discontinuity interface + * \note We do a harmonic average of the part normal to the interface (alpha*I) and + * an arithmetic average of the tangential part (T - alpha*I). + * \return the averaged tensor + * \param T1 first tensor + * \param T2 second tensor + * \param normal The unit normal vector of the interface + */ +template +Dune::FieldMatrix faceTensorAverage(const Dune::FieldMatrix& T1, + const Dune::FieldMatrix& T2, + const Dune::FieldVector& normal) +{ + // determine nT*k*n + Dune::FieldVector tmp; + Dune::FieldVector tmp2; + T1.mv(normal, tmp); + T2.mv(normal, tmp2); + const Scalar alpha1 = tmp*normal; + const Scalar alpha2 = tmp2*normal; + + const Scalar alphaHarmonic = Dumux::harmonicMean(alpha1, alpha2); + const Scalar alphaAverage = 0.5*(alpha1 + alpha2); + + Dune::FieldMatrix T(0.0); + for (int i = 0; i < dim; ++i) + { + for (int j = 0; j < dim; ++j) + { + T[i][j] += 0.5*(T1[i][j] + T2[i][j]); + if (i == j) + T[i][j] += alphaHarmonic - alphaAverage; + } + } + + return T; +} + +} // end namespace Dumux + +#endif diff --git a/dumux/flux/forchheimervelocity.hh b/dumux/flux/forchheimervelocity.hh index bcf1d19cc8bd4d92198eb4a50a014b48419f69a5..7ea9ae297486465a6d0d10752419a3f6ac5a40e7 100644 --- a/dumux/flux/forchheimervelocity.hh +++ b/dumux/flux/forchheimervelocity.hh @@ -34,6 +34,8 @@ #include #include +#include + namespace Dumux { /*! @@ -207,7 +209,7 @@ class ForchheimerVelocity const auto& outsideVolVars = elemVolVars[outsideScvIdx]; const Scalar Kj = getPermeability_(problem, outsideVolVars, scvf.ipGlobal()); const Scalar sqrtKj = sqrt(Kj); - harmonicMeanSqrtK = problem.spatialParams().harmonicMean(sqrtKi, sqrtKj, scvf.unitOuterNormal()); + harmonicMeanSqrtK = faceTensorAverage(sqrtKi, sqrtKj, scvf.unitOuterNormal()); } else harmonicMeanSqrtK = sqrtKi; @@ -259,7 +261,7 @@ class ForchheimerVelocity for (int i = 0; i < dim; ++i) sqrtKj[i][i] = sqrt(Kj[i][i]); - harmonicMeanSqrtK = problem.spatialParams().harmonicMean(sqrtKi, sqrtKj, scvf.unitOuterNormal()); + harmonicMeanSqrtK = faceTensorAverage(sqrtKi, sqrtKj, scvf.unitOuterNormal()); } else harmonicMeanSqrtK = sqrtKi; diff --git a/dumux/material/fluidmatrixinteractions/dispersiontensors/scheidegger.hh b/dumux/material/fluidmatrixinteractions/dispersiontensors/scheidegger.hh index 20773b2ca84d696b2f14f67b2bdfe05a3337772a..d3e1ad30f8c122380ef8a1bf5ba0f553df7215f1 100644 --- a/dumux/material/fluidmatrixinteractions/dispersiontensors/scheidegger.hh +++ b/dumux/material/fluidmatrixinteractions/dispersiontensors/scheidegger.hh @@ -27,6 +27,7 @@ #include #include #include +#include namespace Dumux { @@ -86,9 +87,9 @@ public: // get inside and outside permeability tensors and calculate the harmonic mean const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()]; - const auto K = problem.spatialParams().harmonicMean(insideVolVars.permeability(), - outsideVolVars.permeability(), - scvf.unitOuterNormal()); + const auto K = faceTensorAverage(insideVolVars.permeability(), + outsideVolVars.permeability(), + scvf.unitOuterNormal()); // evaluate gradP - rho*g at integration point Dune::FieldVector gradP(0.0); diff --git a/dumux/porousmediumflow/fvspatialparams.hh b/dumux/porousmediumflow/fvspatialparams.hh index d685ef98151421a0ba4985ca58439b3956b47cd1..824e5b165053371d8820e5f0b54aa973fa0f9107 100644 --- a/dumux/porousmediumflow/fvspatialparams.hh +++ b/dumux/porousmediumflow/fvspatialparams.hh @@ -29,6 +29,7 @@ #include #include +#include // remove include after release 3.5 #include #include @@ -90,6 +91,7 @@ public: * \param T2 second scalar parameter * \param normal The unit normal vector of the interface */ + [[deprecated("Use Dumux::faceTensorAverage from dumux/flux/facetensoraverage.hh instead. This function will be removed after 3.5.")]] Scalar harmonicMean(const Scalar T1, const Scalar T2, const GlobalPosition& normal) const @@ -104,33 +106,12 @@ public: * \param T2 second tensor * \param normal The unit normal vector of the interface */ + [[deprecated("Use Dumux::faceTensorAverage from dumux/flux/facetensoraverage.hh instead. This function will be removed after 3.5.")]] DimWorldMatrix harmonicMean(const DimWorldMatrix& T1, const DimWorldMatrix& T2, const GlobalPosition& normal) const { - // determine nT*k*n - GlobalPosition tmp; - GlobalPosition tmp2; - T1.mv(normal, tmp); - T2.mv(normal, tmp2); - const Scalar alpha1 = tmp*normal; - const Scalar alpha2 = tmp2*normal; - - const Scalar alphaHarmonic = Dumux::harmonicMean(alpha1, alpha2); - const Scalar alphaAverage = 0.5*(alpha1 + alpha2); - - DimWorldMatrix T(0.0); - for (int i = 0; i < dimWorld; ++i) - { - for (int j = 0; j < dimWorld; ++j) - { - T[i][j] += 0.5*(T1[i][j] + T2[i][j]); - if (i == j) - T[i][j] += alphaHarmonic - alphaAverage; - } - } - - return T; + return faceTensorAverage(T1, T2, normal); } /*!