Commit 7969e43a authored by Ned Coltman's avatar Ned Coltman Committed by Melanie Lipp

[rans][problem] migrate dynamic variable calculations to private functions

parent 8aaabe31
......@@ -228,9 +228,6 @@ public:
*/
void updateDynamicWallProperties(const SolutionVector& curSol)
{
using std::abs;
using std::max;
using std::min;
std::cout << "Update dynamic wall properties." << std::endl;
if (!calledUpdateStaticWallProperties)
DUNE_THROW(Dune::InvalidStateException,
......@@ -240,6 +237,161 @@ public:
velocityMaximum_.assign(this->gridGeometry().elementMapper().size(), DimVector(1e-16));
velocityMinimum_.assign(this->gridGeometry().elementMapper().size(), DimVector(std::numeric_limits<Scalar>::max()));
calculateCCVelocities_(curSol);
calculateCCVelocityGradients_();
calculateMaxMinVelocities_();
calculateStressTensor_();
calculateVorticityTensor_();
storeKinematicViscosity_(curSol);
}
/*!
* \brief Returns whether a wall function should be used at a given face
*
* \param element The element.
* \param scvf The sub control volume face.
* \param eqIdx The equation index.
*/
bool useWallFunction(const Element& element,
const SubControlVolumeFace& scvf,
const int& eqIdx) const
{ return false; }
/*!
* \brief Returns an additional wall function momentum flux
*/
template<class ElementVolumeVariables, class ElementFaceVariables>
FacePrimaryVariables wallFunction(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const SubControlVolumeFace& scvf,
const SubControlVolumeFace& lateralBoundaryFace) const
{ return FacePrimaryVariables(0.0); }
/*!
* \brief Returns an additional wall function flux for cell-centered quantities
*/
template<class ElementVolumeVariables, class ElementFaceVariables>
CellCenterPrimaryVariables wallFunction(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const SubControlVolumeFace& scvf) const
{ return CellCenterPrimaryVariables(0.0); }
/*!
* \brief Returns whether a given sub control volume face is on a wall
*
* \param scvf The sub control volume face.
*/
bool isOnWall(const SubControlVolumeFace& scvf) const
{
return asImp_().isOnWallAtPos(scvf.center());
}
/*!
* \brief Returns whether a given point is on a wall
*
* \param globalPos The position in global coordinates.
*/
bool isOnWallAtPos(const GlobalPosition &globalPos) const
{
// Throw an exception if no walls are implemented
DUNE_THROW(Dune::InvalidStateException,
"The problem does not provide an isOnWall() method.");
}
bool hasChannelGeometry() const
{
static const bool hasChannelGeometry = getParamFromGroup<bool>(this->paramGroup(), "RANS.HasChannelGeometry");
return hasChannelGeometry;
}
/*!
* \brief Returns the Karman constant
*/
const Scalar karmanConstant() const
{ return 0.41; }
//! \brief Returns the \f$ \beta_{\omega} \f$ constant
const Scalar betaOmega() const
{
return 0.0708;
}
/*!
* \brief Return the turbulent Prandtl number \f$ [-] \f$ which is used to convert
* the eddy viscosity to an eddy thermal conductivity
*/
Scalar turbulentPrandtlNumber() const
{
static const Scalar turbulentPrandtlNumber
= getParamFromGroup<Scalar>(this->paramGroup(), "RANS.TurbulentPrandtlNumber", 1.0);
return turbulentPrandtlNumber;
}
/*!
* \brief Return the turbulent Schmidt number \f$ [-] \f$ which is used to convert
* the eddy viscosity to an eddy diffusivity
*/
Scalar turbulentSchmidtNumber() const
{
static const Scalar turbulentSchmidtNumber
= getParamFromGroup<Scalar>(this->paramGroup(), "RANS.TurbulentSchmidtNumber", 1.0);
return turbulentSchmidtNumber;
}
int wallNormalAxis(const int elementIdx) const
{ return wallNormalAxis_[elementIdx]; }
int flowNormalAxis(const int elementIdx) const
{ return flowNormalAxis_[elementIdx]; }
unsigned int wallElementIndex(const int elementIdx) const
{ return wallElementIdx_[elementIdx]; }
Scalar wallDistance(const int elementIdx) const
{ return wallDistance_[elementIdx]; }
GlobalPosition cellCenter(const int elementIdx) const
{ return cellCenter_[elementIdx]; }
unsigned int neighborIndex(const int elementIdx, const int dimIdx, const int sideIdx) const
{ return neighborIdx_[elementIdx][dimIdx][sideIdx];}
DimVector ccVelocityVector(const int elementIdx) const
{ return velocity_[elementIdx]; }
Scalar ccVelocity(const int elementIdx, const int dimIdx) const
{ return velocity_[elementIdx][dimIdx]; }
DimVector profileVelocityMaximum(const int elementIdx) const
{ return velocityMaximum_[elementIdx]; }
DimVector profileVelocityMinimum(const int elementIdx) const
{ return velocityMinimum_[elementIdx]; }
DimMatrix velocityGradientTensor(const int elementIdx) const
{ return velocityGradients_[elementIdx]; }
Scalar velocityGradient(const int elementIdx, const int dimIdx, const int velIdx) const
{ return velocityGradients_[elementIdx][dimIdx][velIdx]; }
Scalar stressTensorScalarProduct(const int elementIdx) const
{ return stressTensorScalarProduct_[elementIdx]; }
Scalar vorticityTensorScalarProduct(const int elementIdx) const
{ return vorticityTensorScalarProduct_[elementIdx]; }
Scalar kinematicViscosity(const int elementIdx) const
{ return kinematicViscosity_[elementIdx]; }
bool calledUpdateStaticWallProperties = false;
private:
void calculateCCVelocities_(const SolutionVector& curSol)
{
// calculate cell-center-averaged velocities
for (const auto& element : elements(this->gridGeometry().gridView()))
{
......@@ -258,14 +410,16 @@ public:
for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
velocity_[elementIdx][dimIdx] = velocityTemp[dimIdx] * 0.5; // faces are equidistant to cell center
}
}
void calculateCCVelocityGradients_()
{
using std::abs;
// calculate cell-center-averaged velocity gradients, maximum, and minimum values
for (const auto& element : elements(this->gridGeometry().gridView()))
{
const unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
const unsigned int wallElementIdx = wallElementIdx_[elementIdx];
Scalar maxVelocity = 0.0;
for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
{
for (unsigned int velIdx = 0; velIdx < dim; ++velIdx)
......@@ -280,18 +434,6 @@ public:
if (abs(cellCenter(neighborIndex1)[dimIdx] - cellCenter(neighborIndex0)[dimIdx]) < 1e-8)
velocityGradients_[elementIdx][velIdx][dimIdx] = 0.0;
}
if (abs(ccVelocity(elementIdx, dimIdx)) > abs(velocityMaximum_[wallElementIdx][dimIdx]))
velocityMaximum_[wallElementIdx][dimIdx] = ccVelocity(elementIdx, dimIdx);
if (abs(ccVelocity(elementIdx, dimIdx)) < abs(velocityMinimum_[wallElementIdx][dimIdx]))
velocityMinimum_[wallElementIdx][dimIdx] = ccVelocity(elementIdx, dimIdx);
if ((hasParam("RANS.FlowNormalAxis") != 1) && (maxVelocity) < abs(ccVelocity(elementIdx, dimIdx)))
{
maxVelocity = abs(ccVelocity(elementIdx, dimIdx));
flowNormalAxis_[elementIdx] = dimIdx;
}
}
auto fvGeometry = localView(this->gridGeometry());
......@@ -362,12 +504,67 @@ public:
}
}
}
}
// calculate or call all secondary variables
void calculateMaxMinVelocities_()
{
using std::abs;
if (hasChannelGeometry())
{
// re-initialize min and max values
velocityMaximum_.assign(this->gridGeometry().elementMapper().size(), DimVector(1e-16));
velocityMinimum_.assign(this->gridGeometry().elementMapper().size(), DimVector(std::numeric_limits<Scalar>::max()));
// For each profile perpendicular to the channel wall, find the max and minimum velocities
for (const auto& element : elements(this->gridGeometry().gridView()))
{
const unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
Scalar maxVelocity = 0.0;
const unsigned int wallElementIdx = wallElementIdx_[elementIdx];
for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
{
if (abs(ccVelocity(elementIdx, dimIdx)) > abs(velocityMaximum_[wallElementIdx][dimIdx]))
velocityMaximum_[wallElementIdx][dimIdx] = ccVelocity(elementIdx, dimIdx);
if (abs(ccVelocity(elementIdx, dimIdx)) < abs(velocityMinimum_[wallElementIdx][dimIdx]))
velocityMinimum_[wallElementIdx][dimIdx] = ccVelocity(elementIdx, dimIdx);
if ((hasParam("RANS.FlowNormalAxis") != 1) && (maxVelocity) < abs(ccVelocity(elementIdx, dimIdx)))
{
maxVelocity = abs(ccVelocity(elementIdx, dimIdx));
flowNormalAxis_[elementIdx] = dimIdx;
}
}
}
}
else
{
DimVector maxVelocity(0.0);
DimVector minVelocity(std::numeric_limits<Scalar>::max());
// Find the max and minimum velocities in the full domain
for (const auto& element : elements(this->gridGeometry().gridView()))
{
const unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
{
if (abs(ccVelocity(elementIdx, dimIdx)) > abs(maxVelocity[dimIdx]))
maxVelocity[dimIdx] = ccVelocity(elementIdx, dimIdx);
if (abs(ccVelocity(elementIdx, dimIdx)) < abs(minVelocity[dimIdx]))
minVelocity[dimIdx] = ccVelocity(elementIdx, dimIdx);
}
}
velocityMaximum_.assign(this->gridGeometry().elementMapper().size(), maxVelocity);
velocityMinimum_.assign(this->gridGeometry().elementMapper().size(), minVelocity);
}
}
void calculateStressTensor_()
{
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
Dune::FieldMatrix<Scalar, GridView::dimension, GridView::dimension> stressTensor(0.0);
for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
{
......@@ -385,7 +582,14 @@ public:
stressTensorScalarProduct_[elementIdx] += stressTensor[dimIdx][velIdx] * stressTensor[dimIdx][velIdx];
}
}
}
}
void calculateVorticityTensor_()
{
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
Dune::FieldMatrix<Scalar, GridView::dimension, GridView::dimension> vorticityTensor(0.0);
for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx)
{
......@@ -403,13 +607,20 @@ public:
vorticityTensorScalarProduct_[elementIdx] += vorticityTensor[dimIdx][velIdx] * vorticityTensor[dimIdx][velIdx];
}
}
}
}
void storeKinematicViscosity_(const SolutionVector& curSol)
{
// calculate or call all secondary variables
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
auto fvGeometry = localView(this->gridGeometry());
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
const int dofIdx = scv.dofIndex();
// construct a privars object from the cell center solution vector
const auto& cellCenterPriVars = curSol[GridGeometry::cellCenterIdx()][dofIdx];
PrimaryVariables priVars = makePriVarsFromCellCenterPriVars<PrimaryVariables>(cellCenterPriVars);
......@@ -422,151 +633,6 @@ public:
}
}
/*!
* \brief Returns whether a wall function should be used at a given face
*
* \param element The element.
* \param scvf The sub control volume face.
* \param eqIdx The equation index.
*/
bool useWallFunction(const Element& element,
const SubControlVolumeFace& scvf,
const int& eqIdx) const
{ return false; }
/*!
* \brief Returns an additional wall function momentum flux
*/
template<class ElementVolumeVariables, class ElementFaceVariables>
FacePrimaryVariables wallFunction(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const SubControlVolumeFace& scvf,
const SubControlVolumeFace& lateralBoundaryFace) const
{ return FacePrimaryVariables(0.0); }
/*!
* \brief Returns an additional wall function flux for cell-centered quantities
*/
template<class ElementVolumeVariables, class ElementFaceVariables>
CellCenterPrimaryVariables wallFunction(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const SubControlVolumeFace& scvf) const
{ return CellCenterPrimaryVariables(0.0); }
/*!
* \brief Returns whether a given sub control volume face is on a wall
*
* \param scvf The sub control volume face.
*/
bool isOnWall(const SubControlVolumeFace& scvf) const
{
return asImp_().isOnWallAtPos(scvf.center());
}
/*!
* \brief Returns whether a given point is on a wall
*
* \param globalPos The position in global coordinates.
*/
bool isOnWallAtPos(const GlobalPosition &globalPos) const
{
// Throw an exception if no walls are implemented
DUNE_THROW(Dune::InvalidStateException,
"The problem does not provide an isOnWall() method.");
}
bool hasChannelGeometry() const
{
static const bool hasChannelGeometry = getParamFromGroup<bool>(this->paramGroup(), "RANS.HasChannelGeometry");
return hasChannelGeometry;
}
/*!
* \brief Returns the Karman constant
*/
const Scalar karmanConstant() const
{ return 0.41; }
//! \brief Returns the \f$ \beta_{\omega} \f$ constant
const Scalar betaOmega() const
{
return 0.0708;
}
/*!
* \brief Return the turbulent Prandtl number \f$ [-] \f$ which is used to convert
* the eddy viscosity to an eddy thermal conductivity
*/
Scalar turbulentPrandtlNumber() const
{
static const Scalar turbulentPrandtlNumber
= getParamFromGroup<Scalar>(this->paramGroup(), "RANS.TurbulentPrandtlNumber", 1.0);
return turbulentPrandtlNumber;
}
/*!
* \brief Return the turbulent Schmidt number \f$ [-] \f$ which is used to convert
* the eddy viscosity to an eddy diffusivity
*/
Scalar turbulentSchmidtNumber() const
{
static const Scalar turbulentSchmidtNumber
= getParamFromGroup<Scalar>(this->paramGroup(), "RANS.TurbulentSchmidtNumber", 1.0);
return turbulentSchmidtNumber;
}
int wallNormalAxis(const int elementIdx) const
{ return wallNormalAxis_[elementIdx]; }
int flowNormalAxis(const int elementIdx) const
{ return flowNormalAxis_[elementIdx]; }
unsigned int wallElementIndex(const int elementIdx) const
{ return wallElementIdx_[elementIdx]; }
Scalar wallDistance(const int elementIdx) const
{ return wallDistance_[elementIdx]; }
GlobalPosition cellCenter(const int elementIdx) const
{ return cellCenter_[elementIdx]; }
unsigned int neighborIndex(const int elementIdx, const int dimIdx, const int sideIdx) const
{ return neighborIdx_[elementIdx][dimIdx][sideIdx];}
DimVector ccVelocityVector(const int elementIdx) const
{ return velocity_[elementIdx]; }
Scalar ccVelocity(const int elementIdx, const int dimIdx) const
{ return velocity_[elementIdx][dimIdx]; }
DimVector profileVelocityMaximum(const int elementIdx) const
{ return velocityMaximum_[elementIdx]; }
DimVector profileVelocityMinimum(const int elementIdx) const
{ return velocityMinimum_[elementIdx]; }
DimMatrix velocityGradientTensor(const int elementIdx) const
{ return velocityGradients_[elementIdx]; }
Scalar velocityGradient(const int elementIdx, const int dimIdx, const int velIdx) const
{ return velocityGradients_[elementIdx][dimIdx][velIdx]; }
Scalar stressTensorScalarProduct(const int elementIdx) const
{ return stressTensorScalarProduct_[elementIdx]; }
Scalar vorticityTensorScalarProduct(const int elementIdx) const
{ return vorticityTensorScalarProduct_[elementIdx]; }
Scalar kinematicViscosity(const int elementIdx) const
{ return kinematicViscosity_[elementIdx]; }
bool calledUpdateStaticWallProperties = false;
private:
const int fixedFlowNormalAxis_ = getParam<int>("RANS.FlowNormalAxis", 0);
const int fixedWallNormalAxis_ = getParam<int>("RANS.WallNormalAxis", 1);
......
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