From 5e72f03e6ba73007d4f6f91097c30fdcdefeb4f1 Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Thu, 31 Mar 2022 10:30:36 +0200 Subject: [PATCH 1/2] [flux] extract tensor average at interface --- dumux/flux/facetensoraverage.hh | 87 +++++++++++++++++++++++ dumux/porousmediumflow/fvspatialparams.hh | 27 ++----- 2 files changed, 91 insertions(+), 23 deletions(-) create mode 100644 dumux/flux/facetensoraverage.hh diff --git a/dumux/flux/facetensoraverage.hh b/dumux/flux/facetensoraverage.hh new file mode 100644 index 0000000000..38a1f5b9f8 --- /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/porousmediumflow/fvspatialparams.hh b/dumux/porousmediumflow/fvspatialparams.hh index d685ef9815..824e5b1650 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); } /*! -- GitLab From 92fcdc6ed90f9b46ca0790f7465e65b5bdfa52b4 Mon Sep 17 00:00:00 2001 From: Mathis Kelm Date: Thu, 31 Mar 2022 11:04:56 +0200 Subject: [PATCH 2/2] [flux] adjust deprecated interface average calls --- dumux/flux/box/darcyslaw.hh | 5 +++-- dumux/flux/box/fickslaw.hh | 3 ++- dumux/flux/box/forchheimerslaw.hh | 3 ++- dumux/flux/box/fourierslaw.hh | 3 ++- dumux/flux/box/fourierslawnonequilibrium.hh | 4 +++- dumux/flux/box/maxwellstefanslaw.hh | 5 +++-- dumux/flux/forchheimervelocity.hh | 6 ++++-- .../dispersiontensors/scheidegger.hh | 7 ++++--- 8 files changed, 23 insertions(+), 13 deletions(-) diff --git a/dumux/flux/box/darcyslaw.hh b/dumux/flux/box/darcyslaw.hh index a1326195e3..5662119778 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 0349e2549c..f802f6ec7f 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 0ce29b5f0e..3ea0759f6d 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 6a825e6946..0a55689112 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 22ea175d80..c6b33d9597 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 fa2a845e9a..758e80e18c 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/forchheimervelocity.hh b/dumux/flux/forchheimervelocity.hh index bcf1d19cc8..7ea9ae2974 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 20773b2ca8..d3e1ad30f8 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); -- GitLab