Commit 20e7be33 authored by Timo Koch's avatar Timo Koch
Browse files

[md][embedded][kernel] Replace couplmanager with a recent kernel coupling scheme

Coupling methox from Koch et al 2020 JCP (https://doi.org/10.1016/j.jcp.2020.109370)
parent 88edf0eb
......@@ -18,9 +18,10 @@ Differences Between DuMu<sup>x</sup> 3.3 and DuMu<sup>x</sup> 3.2
- __Quadmath__: Dumux::Quad has been removed without deprecation. Use Dune::Float128 instead.
- Within the RANS group, two additional runtime parameters have been included 'IsFlatWallBounded' and 'WriteFlatWallBoundedFields'.
For both the K-Epsilon and Zero-eq RANS models the 'IsFlatWallBounded' runtime parameter should be set as True,
as wall topology is not supported for these models with our geometric contraints. If not set as true, the geometry
will be checked before the model is run. If either the runtime parameter or the geometry check indicate non-flat walls,
as wall topology is not supported for these models with our geometric contraints. If not set as true, the geometry
will be checked before the model is run. If either the runtime parameter or the geometry check indicate non-flat walls,
the model will terminate. To add FlatWallBounded specific output to the vtk output, WriteFlatWallBoundedFields can be set as True.
- __1d3d coupling__: The kernel coupling manager has been replaced with the one from Koch et al (2020) JCP https://doi.org/10.1016/j.jcp.2020.109370
### Deprecated properties/classes/functions/files, to be removed after 3.3:
......
[MixedDimension]
NumCircleSegments = 100
IntegrationOrder = 2
KernelWidth = 0.3
KernelWidthFactor = 1.0
[Vessel.Grid]
LowerLeft = 0 0 0
......
......@@ -256,9 +256,7 @@ public:
* the absolute rate mass generated or annihilated in kg/s. Positive values mean
* that mass is created, negative ones mean that it vanishes.
*/
template<class ElementVolumeVariables,
bool enable = (CouplingManager::couplingMode == Embedded1d3dCouplingMode::kernel),
std::enable_if_t<!enable, int> = 0>
template<class ElementVolumeVariables>
void pointSource(PointSource& source,
const Element &element,
const FVElementGeometry& fvGeometry,
......@@ -272,60 +270,13 @@ public:
// calculate the source
const Scalar radius = this->couplingManager().radius(source.id());
const Scalar beta = 2*M_PI/(2*M_PI + std::log(radius));
const Scalar sourceValue = beta*(pressure3D - pressure1D);//*bulkVolVars.density();
Scalar sourceValue = beta*(pressure3D - pressure1D);
if constexpr(CouplingManager::couplingMode == EmbeddedCouplingMode::kernel)
sourceValue *= this->couplingManager().fluxScalingFactor(source.id());
source = sourceValue*source.quadratureWeight()*source.integrationElement();
}
//! Specialization for kernel method
template<class ElementVolumeVariables,
bool enable = (CouplingManager::couplingMode == Embedded1d3dCouplingMode::kernel),
std::enable_if_t<enable, int> = 0>
void pointSource(PointSource& source,
const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolume &scv) const
{
static const Scalar kernelWidth = getParam<Scalar>("MixedDimension.KernelWidth");
static const Scalar kernelFactor = std::log(kernelWidth)-9.0/10.0;
// compute source at every integration point
const Scalar pressure1D = this->couplingManager().lowDimPriVars(source.id())[Indices::pressureIdx];
const Scalar pressure3D = this->couplingManager().bulkPriVars(source.id())[Indices::pressureIdx];
// calculate the source
const Scalar radius = this->couplingManager().radius(source.id());
const Scalar beta = 2*M_PI/(2*M_PI + std::log(radius));
const Scalar sourceValue = beta*(pressure3D / kernelFactor * std::log(radius) - pressure1D);//*bulkVolVars.density();
source = sourceValue*source.quadratureWeight()*source.integrationElement();
}
//! Evaluate coupling residual for the derivative bulk DOF with respect to low dim DOF
//! We only need to evaluate the part of the residual that will be influence by the low dim DOF
template<class MatrixBlock, class VolumeVariables>
void addSourceDerivatives(MatrixBlock& block,
const Element& element,
const FVElementGeometry& fvGeometry,
const VolumeVariables& curElemVolVars,
const SubControlVolume& scv) const
{
const auto eIdx = this->gridGeometry().elementMapper().index(element);
auto key = std::make_pair(eIdx, 0);
if (this->pointSourceMap().count(key))
{
// call the solDependent function. Herein the user might fill/add values to the point sources
// we make a copy of the local point sources here
auto pointSources = this->pointSourceMap().at(key);
// add the point source values to the local residual (negative sign is convention for source term)
for (const auto& source : pointSources)
block[0][0] -= this->couplingManager().pointSourceDerivative(source, Dune::index_constant<1>{}, Dune::index_constant<1>{});
}
}
/*!
* \brief Evaluates the initial value for a control volume.
*
......
......@@ -247,15 +247,18 @@ public:
const ElementVolumeVariables& elemVolVars,
const SubControlVolume &scv) const
{
// compute source at every integration point
const Scalar pressure3D = this->couplingManager().bulkPriVars(source.id())[Indices::pressureIdx];
const Scalar pressure1D = this->couplingManager().lowDimPriVars(source.id())[Indices::pressureIdx];
// calculate the source
const Scalar radius = this->couplingManager().radius(source.id());
const Scalar beta = 2*M_PI/(2*M_PI + std::log(radius));
const Scalar sourceValue = beta*(pressure1D - pressure3D);
source = sourceValue*source.quadratureWeight()*source.integrationElement();
if constexpr (CouplingManager::couplingMode != Embedded1d3dCouplingMode::kernel)
{
// compute source at every integration point
const Scalar pressure3D = this->couplingManager().bulkPriVars(source.id())[Indices::pressureIdx];
const Scalar pressure1D = this->couplingManager().lowDimPriVars(source.id())[Indices::pressureIdx];
// calculate the source
const Scalar radius = this->couplingManager().radius(source.id());
const Scalar beta = 2*M_PI/(2*M_PI + std::log(radius));
const Scalar sourceValue = beta*(pressure1D - pressure3D);
source = sourceValue*source.quadratureWeight()*source.integrationElement();
}
}
/*!
......@@ -276,77 +279,51 @@ public:
* that the conserved quantity is created, negative ones mean that it vanishes.
* E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$.
*/
template<class ElementVolumeVariables,
bool enable = (CouplingManager::couplingMode == Embedded1d3dCouplingMode::kernel),
std::enable_if_t<enable, int> = 0>
template<class ElementVolumeVariables>
NumEqVector source(const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolume &scv) const
{
static const Scalar kernelWidth = getParam<Scalar>("MixedDimension.KernelWidth");
static const Scalar kernelFactor = std::log(kernelWidth)-9.0/10.0;
NumEqVector source(0.0);
const auto eIdx = this->gridGeometry().elementMapper().index(element);
const auto& sourceIds = this->couplingManager().bulkSourceIds(eIdx);
const auto& sourceWeights = this->couplingManager().bulkSourceWeights(eIdx);
for (int i = 0; i < sourceIds.size(); ++i)
if constexpr(CouplingManager::couplingMode == EmbeddedCouplingMode::kernel)
{
const auto id = sourceIds[i];
const auto weight = sourceWeights[i];
const Scalar radius = this->couplingManager().radius(id);
const Scalar pressure3D = this->couplingManager().bulkPriVars(id)[Indices::pressureIdx];
const Scalar pressure1D = this->couplingManager().lowDimPriVars(id)[Indices::pressureIdx];
const auto eIdx = this->gridGeometry().elementMapper().index(element);
const auto& sourceIds = this->couplingManager().bulkSourceIds(eIdx, scv.indexInElement());
const auto& sourceWeights = this->couplingManager().bulkSourceWeights(eIdx, scv.indexInElement());
// calculate the source
const Scalar beta = 2*M_PI/(2*M_PI + std::log(radius));
const Scalar sourceValue = beta*(pressure1D - pressure3D / kernelFactor * std::log(radius));
for (int i = 0; i < sourceIds.size(); ++i)
{
const auto id = sourceIds[i];
const auto weight = sourceWeights[i];
const auto xi = this->couplingManager().fluxScalingFactor(id);
const Scalar radius = this->couplingManager().radius(id);
const Scalar pressure3D = this->couplingManager().bulkPriVars(id)[Indices::pressureIdx];
const Scalar pressure1D = this->couplingManager().lowDimPriVars(id)[Indices::pressureIdx];
// calculate the source
static const Scalar beta = 2*M_PI/(2*M_PI + std::log(radius));
const Scalar sourceValue = beta*(pressure1D - pressure3D)*xi;
source[Indices::conti0EqIdx] += sourceValue*weight;
}
source[Indices::conti0EqIdx] += sourceValue*weight;
const auto volume = scv.volume()*elemVolVars[scv].extrusionFactor();
source[Indices::conti0EqIdx] /= volume;
}
const auto volume = scv.volume()*elemVolVars[scv].extrusionFactor();
source[Indices::conti0EqIdx] /= volume;
return source;
}
//! Other methods
template<class ElementVolumeVariables,
bool enable = (CouplingManager::couplingMode == Embedded1d3dCouplingMode::kernel),
std::enable_if_t<!enable, int> = 0>
NumEqVector source(const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolume &scv) const
{ return NumEqVector(0.0); }
//! Evaluate coupling residual for the derivative bulk DOF with respect to low dim DOF
//! We only need to evaluate the part of the residual that will be influence by the low dim DOF
template<class MatrixBlock, class VolumeVariables>
void addSourceDerivatives(MatrixBlock& block,
const Element& element,
const FVElementGeometry& fvGeometry,
const VolumeVariables& curElemVolVars,
const SubControlVolume& scv) const
//! compute the flux scaling factor (xi) for a distance r for a vessel with radius R and kernel width rho
Scalar fluxScalingFactor(const Scalar r, const Scalar R, const Scalar rho) const
{
const auto eIdx = this->gridGeometry().elementMapper().index(element);
auto key = std::make_pair(eIdx, 0);
if (this->pointSourceMap().count(key))
{
// call the solDependent function. Herein the user might fill/add values to the point sources
// we make a copy of the local point sources here
auto pointSources = this->pointSourceMap().at(key);
// add the point source values to the local residual (negative sign is convention for source term)
for (const auto& source : pointSources)
block[0][0] -= this->couplingManager().pointSourceDerivative(source, Dune::index_constant<0>{}, Dune::index_constant<0>{});
}
using std::log;
static const Scalar beta = 2*M_PI/(2*M_PI + log(R));
static const Scalar kernelWidthFactorLog = log(rho/R);
return r < rho ? 1.0/(1.0 + R*beta/(2*M_PI*R)*(r*r/(2*rho*rho) + kernelWidthFactorLog - 0.5))
: 1.0/(1.0 + R*beta/(2*M_PI*R)*log(r/R));
}
/*!
......@@ -360,33 +337,32 @@ public:
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{ return PrimaryVariables(0.0); }
Scalar p1DExact(const GlobalPosition &globalPos) const
{ return 1.0 + globalPos[2]; }
//! The exact solution
Scalar exactSolution(const GlobalPosition &globalPos) const
{
Dune::FieldVector<double, 2> xy({globalPos[0], globalPos[1]});
if (CouplingManager::couplingMode == Embedded1d3dCouplingMode::surface)
const auto r = std::sqrt(globalPos[0]*globalPos[0] + globalPos[1]*globalPos[1]);
static const auto R = getParam<Scalar>("SpatialParams.Radius");
if (CouplingManager::couplingMode == Embedded1d3dCouplingMode::kernel)
{
static const auto R = getParam<Scalar>("SpatialParams.Radius");
if (xy.two_norm() > R)
return -1.0*(1+globalPos[2])/(2*M_PI)*std::log(xy.two_norm());
static const auto rho = getParam<Scalar>("MixedDimension.KernelWidthFactor")*R;
if (r > rho)
return -1.0*p1DExact(globalPos)/(2*M_PI)*std::log(r);
else
return -1.0*(1+globalPos[2])/(2*M_PI)*std::log(R);
return -1.0*p1DExact(globalPos)/(2*M_PI)*(r*r/(2.0*rho*rho) + std::log(rho) - 0.5);
}
else if (CouplingManager::couplingMode == Embedded1d3dCouplingMode::kernel)
else if (CouplingManager::couplingMode == Embedded1d3dCouplingMode::surface)
{
static const auto rho = getParam<Scalar>("MixedDimension.KernelWidth");
const auto& r = xy.two_norm();
if (xy.two_norm() >= rho)
return -1.0*(1+globalPos[2])/(2*M_PI)*std::log(xy.two_norm());
if (r > R)
return -1.0*p1DExact(globalPos)/(2*M_PI)*std::log(r);
else
return -1.0*(1+globalPos[2])/(2*M_PI)*(11.0/15.0*std::pow(r/rho, 5)
- 3.0/2.0*std::pow(r/rho, 4)
+ 5.0/3.0*std::pow(r/rho, 2)
- 9.0/10.0 + std::log(rho));
return -1.0*p1DExact(globalPos)/(2*M_PI)*std::log(R);
}
else
return -1.0*(1+globalPos[2])/(2*M_PI)*std::log(xy.two_norm());
return -1.0*p1DExact(globalPos)/(2*M_PI)*std::log(r);
}
//! Called after every time step
......
Markdown is supported
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