Commit 832686a3 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

[stokes] remove unnecessary loops over coordinate index

Several explicit loops over the coordinate index can be replaced by
using the corresponding FieldVector functionality.

This implements FS#271.

Reviewed by fetzer.


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@15194 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 14bc7a84
......@@ -91,28 +91,27 @@ protected:
velocity_ = Scalar(0);
pressureGrad_ = Scalar(0);
velocityGrad_ = Scalar(0);
// velocityDiv_ = Scalar(0);
for (int idx = 0;
idx < fvGeometry_.numScv;
idx++) // loop over adjacent vertices
for (int scvIdx = 0;
scvIdx < fvGeometry_.numScv;
scvIdx++) // loop over adjacent vertices
{
// phase density and viscosity at IP
density_ += elemVolVars[idx].density() *
face().shapeValue[idx];
dynamicViscosity_ += elemVolVars[idx].dynamicViscosity() *
face().shapeValue[idx];
pressure_ += elemVolVars[idx].pressure() *
face().shapeValue[idx];
density_ += elemVolVars[scvIdx].density() *
face().shapeValue[scvIdx];
dynamicViscosity_ += elemVolVars[scvIdx].dynamicViscosity() *
face().shapeValue[scvIdx];
pressure_ += elemVolVars[scvIdx].pressure() *
face().shapeValue[scvIdx];
// velocity at the IP (fluxes)
DimVector velocityTimesShapeValue = elemVolVars[idx].velocity();
velocityTimesShapeValue *= face().shapeValue[idx];
DimVector velocityTimesShapeValue = elemVolVars[scvIdx].velocity();
velocityTimesShapeValue *= face().shapeValue[scvIdx];
velocity_ += velocityTimesShapeValue;
// the pressure gradient
tmp = face().grad[idx];
tmp *= elemVolVars[idx].pressure();
tmp = face().grad[scvIdx];
tmp *= elemVolVars[scvIdx].pressure();
pressureGrad_ += tmp;
// take gravity into account
tmp = problem.gravity();
......@@ -121,13 +120,11 @@ protected:
pressureGrad_ -= tmp;
// the velocity gradients and divergence
for (int dimIdx = 0; dimIdx<dim; ++dimIdx)
for (int dimIdx = 0; dimIdx < dim; ++dimIdx)
{
tmp = face().grad[idx];
tmp *= elemVolVars[idx].velocity()[dimIdx];
tmp = face().grad[scvIdx];
tmp *= elemVolVars[scvIdx].velocity()[dimIdx];
velocityGrad_[dimIdx] += tmp;
// velocityDiv_ += face().grad[idx][dimIdx]*elemVolVars[idx].velocity()[dimIdx];
}
}
......@@ -139,7 +136,6 @@ protected:
Valgrind::CheckDefined(velocity_);
Valgrind::CheckDefined(pressureGrad_);
Valgrind::CheckDefined(velocityGrad_);
// Valgrind::CheckDefined(velocityDiv_);
}
void determineUpwindDirection_(const ElementVolumeVariables &elemVolVars)
......@@ -244,14 +240,6 @@ public:
const Scalar kinematicEddyViscosity() const
{ return 0; }
// /*!
// * \brief Return the divergence of the normal velocity at the
// * integration point.
// */
// Scalar velocityDiv() const
// { return velocityDiv_; }
/*!
* \brief Return the local index of the upstream sub-control volume.
*/
......@@ -280,7 +268,6 @@ protected:
Scalar dynamicViscosity_;
Scalar pressure_;
Scalar normalvelocity_;
// Scalar velocityDiv_;
DimVector velocity_;
// gradients at the IPs
......
......@@ -324,7 +324,7 @@ protected:
// add the component of the pressure gradient to the respective part
// of the momentum equation and take the gravity term into account
// signs are inverted, since q is subtracted
for (int dimIdx=0; dimIdx<dim; ++dimIdx)
for (int dimIdx = 0; dimIdx < dim; ++dimIdx)
{
source[momentumXIdx + dimIdx] -= pressureGradAtSCVCenter[dimIdx];
source[momentumXIdx + dimIdx] += volVars.density()*this->problem_().gravity()[dimIdx];
......@@ -413,7 +413,7 @@ protected:
for (unsigned int i=0; i < this->residual_.size(); i++)
Valgrind::CheckDefined(this->residual_[i]);
for (int dimIdx=0; dimIdx < dim; ++dimIdx)
for (int dimIdx = 0; dimIdx < dim; ++dimIdx)
momentumResidual[dimIdx] = this->residual_[scvIdx][momentumXIdx+dimIdx];
//Sign is right!!!: boundary flux: -mu grad v n
......
......@@ -57,33 +57,33 @@ class StokesncFluxVariables : public StokesFluxVariables<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
//dimensions
enum { dim = GridView::dimension };
enum { dim = GridView::dimension };
//phase indices
enum { phaseIdx = Indices::phaseIdx };
enum { phaseIdx = Indices::phaseIdx };
//component indices
enum { phaseCompIdx = Indices::phaseCompIdx,
transportCompIdx = Indices::transportCompIdx };
enum { phaseCompIdx = Indices::phaseCompIdx,
transportCompIdx = Indices::transportCompIdx };
//number of components
enum { numComponents = Indices::numComponents };
enum { numComponents = Indices::numComponents };
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Entity Element;
typedef Dune::FieldVector<Scalar, dim> DimVector;
typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
public:
//Constrcutor calls ParentType function
StokesncFluxVariables(const Problem &problem,
//Constructor calls ParentType function
StokesncFluxVariables(const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry,
const int fIdx,
const ElementVolumeVariables &elemVolVars,
const bool onBoundary = false)
: ParentType(problem, element, fvGeometry, fIdx, elemVolVars, onBoundary)
: ParentType(problem, element, fvGeometry, fIdx, elemVolVars, onBoundary)
{
calculateValues_(problem, element, elemVolVars);
}
/*!
/*!
* \brief Return the molar density \f$ \mathrm{[mol/m^3]} \f$ at the integration point.
*/
const Scalar molarDensity() const
......@@ -113,7 +113,7 @@ public:
const DimVector &moleFractionGrad(int compIdx) const
{ return moleFractionGrad_[compIdx]; }
/*!
/*!
* \brief Return the eddy diffusivity (if implemented).
*/
const Scalar eddyDiffusivity() const
......@@ -125,51 +125,48 @@ protected:
const ElementVolumeVariables &elemVolVars)
{
// loop over all components
for (int compIdx=0; compIdx<numComponents; compIdx++){
// loop over all components
for (int compIdx=0; compIdx<numComponents; compIdx++){
massFraction_[compIdx] = Scalar(0.0);
moleFraction_[compIdx] = Scalar(0.0);
diffusionCoeff_[compIdx] = Scalar(0.0);
moleFractionGrad_[compIdx] = Scalar(0.0);
if (phaseCompIdx!=compIdx) //no transport equation parameters needed for the mass balance
{
if (phaseCompIdx!=compIdx) //no transport equation parameters needed for the mass balance
{
molarDensity_ = Scalar(0.0);
// calculate gradients and secondary variables at IPs
for (int scvIdx = 0;
scvIdx < this->fvGeometry_.numScv;
scvIdx++) // loop over vertices of the element
{
molarDensity_ += elemVolVars[scvIdx].molarDensity()*
this->face().shapeValue[scvIdx];
massFraction_[compIdx] += elemVolVars[scvIdx].massFraction(compIdx) *
this->face().shapeValue[scvIdx];
moleFraction_[compIdx] += elemVolVars[scvIdx].moleFraction(compIdx) *
this->face().shapeValue[scvIdx];
diffusionCoeff_[compIdx] += elemVolVars[scvIdx].diffusionCoeff(compIdx) *
this->face().shapeValue[scvIdx];
// the gradient of the mole fraction at the IP
for (int dimIdx=0; dimIdx<dim; ++dimIdx)
{
moleFractionGrad_[compIdx][dimIdx] +=
this->face().grad[scvIdx][dimIdx] *
elemVolVars[scvIdx].moleFraction(compIdx);
}
}
Valgrind::CheckDefined(molarDensity_);
Valgrind::CheckDefined(massFraction_[compIdx]);
Valgrind::CheckDefined(moleFraction_[compIdx]);
Valgrind::CheckDefined(diffusionCoeff_[compIdx]);
Valgrind::CheckDefined(moleFractionGrad_[compIdx]);
}
}
// calculate gradients and secondary variables at IPs
for (int scvIdx = 0;
scvIdx < this->fvGeometry_.numScv;
scvIdx++) // loop over vertices of the element
{
molarDensity_ += elemVolVars[scvIdx].molarDensity()*
this->face().shapeValue[scvIdx];
massFraction_[compIdx] += elemVolVars[scvIdx].massFraction(compIdx) *
this->face().shapeValue[scvIdx];
moleFraction_[compIdx] += elemVolVars[scvIdx].moleFraction(compIdx) *
this->face().shapeValue[scvIdx];
diffusionCoeff_[compIdx] += elemVolVars[scvIdx].diffusionCoeff(compIdx) *
this->face().shapeValue[scvIdx];
// the gradient of the mole fraction at the IP
DimVector grad = this->face().grad[scvIdx];
grad *= elemVolVars[scvIdx].moleFraction(compIdx);
moleFractionGrad_[compIdx] += grad;
}
Valgrind::CheckDefined(molarDensity_);
Valgrind::CheckDefined(massFraction_[compIdx]);
Valgrind::CheckDefined(moleFraction_[compIdx]);
Valgrind::CheckDefined(diffusionCoeff_[compIdx]);
Valgrind::CheckDefined(moleFractionGrad_[compIdx]);
}
}
}
Scalar molarDensity_;
Scalar molarDensity_;
Scalar massFraction_[numComponents];
Scalar moleFraction_[numComponents];
Scalar diffusionCoeff_[numComponents];
......
......@@ -201,36 +201,31 @@ public:
void computeDiffusiveFlux(PrimaryVariables &flux,
const FluxVariables &fluxVars) const
{
// diffusive component flux
for (int dimIdx = 0; dimIdx < dim; ++dimIdx)
// diffusive component flux
if(!useMoles)
{
if(!useMoles)
{
flux[transportEqIdx] -= fluxVars.moleFractionGrad(transportCompIdx)[dimIdx]
* fluxVars.face().normal[dimIdx]
*(fluxVars.diffusionCoeff(transportCompIdx) + fluxVars.eddyDiffusivity())
* fluxVars.molarDensity()
* FluidSystem::molarMass(transportCompIdx);// Multiplied by molarMass [kg/mol] to convert form [mol/m^3 s] to [kg/m^3 s]
Valgrind::CheckDefined(flux[transportEqIdx]);
}
else
{
//loop over secondary components
for (int compIdx=0; compIdx<numComponents; compIdx++)
{
if (conti0EqIdx+compIdx != massBalanceIdx)
{
flux[conti0EqIdx+compIdx] -= fluxVars.moleFractionGrad(compIdx)[dimIdx]
* fluxVars.face().normal[dimIdx]
*(fluxVars.diffusionCoeff(compIdx) + fluxVars.eddyDiffusivity())
* fluxVars.molarDensity();
Valgrind::CheckDefined(flux[conti0EqIdx+compIdx]);
}
}
flux[transportEqIdx] -= fluxVars.moleFractionGrad(transportCompIdx)
* fluxVars.face().normal
* (fluxVars.diffusionCoeff(transportCompIdx) + fluxVars.eddyDiffusivity())
* fluxVars.molarDensity()
* FluidSystem::molarMass(transportCompIdx);// Multiplied by molarMass [kg/mol] to convert form [mol/m^3 s] to [kg/m^3 s]
Valgrind::CheckDefined(flux[transportEqIdx]);
}
else
{
//loop over secondary components
for (int compIdx=0; compIdx<numComponents; compIdx++)
{
if (conti0EqIdx+compIdx != massBalanceIdx)
{
flux[conti0EqIdx+compIdx] -= fluxVars.moleFractionGrad(compIdx)
* fluxVars.face().normal
* (fluxVars.diffusionCoeff(compIdx) + fluxVars.eddyDiffusivity())
* fluxVars.molarDensity();
Valgrind::CheckDefined(flux[conti0EqIdx+compIdx]);
}
}
}
}
}
};
......
......@@ -126,25 +126,23 @@ protected:
temperatureGrad_ = Scalar(0);
// calculate gradients and secondary variables at IPs
DimVector tmp(0.0);
for (int idx = 0;
idx < this->fvGeometry_.numScv;
idx++) // loop over vertices of the element
for (int scvIdx = 0;
scvIdx < this->fvGeometry_.numScv;
scvIdx++) // loop over vertices of the element
{
temperature_ += elemVolVars[idx].temperature() *
this->face().shapeValue[idx];
temperature_ += elemVolVars[scvIdx].temperature() *
this->face().shapeValue[scvIdx];
thermalConductivity_ += elemVolVars[idx].thermalConductivity() *
this->face().shapeValue[idx];
thermalConductivity_ += elemVolVars[scvIdx].thermalConductivity() *
this->face().shapeValue[scvIdx];
heatCapacity_ += elemVolVars[idx].heatCapacity() *
this->face().shapeValue[idx];
heatCapacity_ += elemVolVars[scvIdx].heatCapacity() *
this->face().shapeValue[scvIdx];
// the gradient of the temperature at the IP
for (int dimIdx=0; dimIdx<dim; ++dimIdx)
temperatureGrad_[dimIdx] +=
this->face().grad[idx][dimIdx]*
elemVolVars[idx].temperature();
DimVector grad = this->face().grad[scvIdx];
grad *= elemVolVars[scvIdx].temperature();
temperatureGrad_ += grad;
}
Valgrind::CheckDefined(temperature_);
Valgrind::CheckDefined(thermalConductivity_);
......@@ -154,10 +152,10 @@ protected:
for (unsigned int i = 0; i < numComponents; ++i)
{
componentEnthalpy_[i] = Scalar(0.0);
for (int idx = 0; idx < this->fvGeometry_.numScv; idx++) // loop over vertices of the element
for (int scvIdx = 0; scvIdx < this->fvGeometry_.numScv; scvIdx++) // loop over vertices of the element
{
componentEnthalpy_[i] += elemVolVars[idx].componentEnthalpy(i)
* this->face().shapeValue[idx];
componentEnthalpy_[i] += elemVolVars[scvIdx].componentEnthalpy(i)
* this->face().shapeValue[scvIdx];
}
Valgrind::CheckDefined(componentEnthalpy_[i]);
}
......
......@@ -133,34 +133,32 @@ public:
computeConductiveFlux(flux, fluxVars);
// diffusive component energy flux
for (int dimIdx = 0; dimIdx < dim; ++dimIdx)
Scalar sumDiffusiveFluxes = 0;
for (int compIdx=0; compIdx<numComponents; compIdx++)
{
Scalar sumDiffusiveFluxes = 0;
for (int compIdx=0; compIdx<numComponents; compIdx++)
if (compIdx != phaseCompIdx)
{
if (compIdx != phaseCompIdx)
{
Valgrind::CheckDefined(fluxVars.moleFractionGrad(compIdx)[dimIdx]);
Valgrind::CheckDefined(fluxVars.face().normal[dimIdx]);
Valgrind::CheckDefined(fluxVars.diffusionCoeff(compIdx));
Valgrind::CheckDefined(fluxVars.eddyDiffusivity());
Valgrind::CheckDefined(fluxVars.molarDensity());
Valgrind::CheckDefined(FluidSystem::molarMass(compIdx));
Valgrind::CheckDefined(fluxVars.componentEnthalpy(compIdx));
Scalar diffusiveFlux = fluxVars.moleFractionGrad(compIdx)[dimIdx]
* fluxVars.face().normal[dimIdx]
*(fluxVars.diffusionCoeff(compIdx) + fluxVars.eddyDiffusivity())
* fluxVars.molarDensity();
sumDiffusiveFluxes += diffusiveFlux;
flux[energyEqIdx] -= diffusiveFlux * fluxVars.componentEnthalpy(compIdx)
* FluidSystem::molarMass(compIdx); // Multiplied by molarMass [kg/mol] to convert from [mol/m^3 s] to [kg/m^3 s];
}
Valgrind::CheckDefined(fluxVars.moleFractionGrad(compIdx));
Valgrind::CheckDefined(fluxVars.face().normal);
Valgrind::CheckDefined(fluxVars.diffusionCoeff(compIdx));
Valgrind::CheckDefined(fluxVars.eddyDiffusivity());
Valgrind::CheckDefined(fluxVars.molarDensity());
Valgrind::CheckDefined(FluidSystem::molarMass(compIdx));
Valgrind::CheckDefined(fluxVars.componentEnthalpy(compIdx));
Scalar diffusiveFlux = fluxVars.moleFractionGrad(compIdx)
* fluxVars.face().normal
* (fluxVars.diffusionCoeff(compIdx) + fluxVars.eddyDiffusivity())
* fluxVars.molarDensity();
sumDiffusiveFluxes += diffusiveFlux;
flux[energyEqIdx] -= diffusiveFlux * fluxVars.componentEnthalpy(compIdx)
* FluidSystem::molarMass(compIdx); // Multiplied by molarMass [kg/mol] to convert from [mol/m^3 s] to [kg/m^3 s];
}
// the diffusive flux of the phase component is the negative of the sum of the component fluxes
flux[energyEqIdx] += sumDiffusiveFluxes * fluxVars.componentEnthalpy(phaseCompIdx)
* FluidSystem::molarMass(phaseCompIdx); // Multiplied by molarMass [kg/mol] to convert from [mol/m^3 s] to [kg/m^3 s];
}
// the diffusive flux of the phase component is the negative of the sum of the component fluxes
flux[energyEqIdx] += sumDiffusiveFluxes * fluxVars.componentEnthalpy(phaseCompIdx)
* FluidSystem::molarMass(phaseCompIdx); // Multiplied by molarMass [kg/mol] to convert from [mol/m^3 s] to [kg/m^3 s];
Valgrind::CheckDefined(flux[energyEqIdx]);
}
......@@ -175,13 +173,9 @@ public:
const FluxVariables &fluxVars) const
{
// diffusive heat flux
for (int dimIdx = 0; dimIdx < dim; ++dimIdx)
{
flux[energyEqIdx] -=
fluxVars.temperatureGrad()[dimIdx] *
fluxVars.face().normal[dimIdx] *
(fluxVars.thermalConductivity() + fluxVars.thermalEddyConductivity());
}
flux[energyEqIdx] -=
fluxVars.temperatureGrad() * fluxVars.face().normal
* (fluxVars.thermalConductivity() + fluxVars.thermalEddyConductivity());
Valgrind::CheckDefined(flux[energyEqIdx]);
}
};
......
......@@ -170,13 +170,7 @@ protected:
fz_ = distanceToWallRough_ * omega * (1.0 - std::exp(-yPlusRough_ / aPlus ));
Scalar fMax = problem.model().wall[wallIdx_].fMax[posIdx_];
Scalar yMax = std::abs(problem.model().wall[wallIdx_].yMax[posIdx_]);
Scalar temp[2] = {0.0, 0.0};
for (int dimIdx = 0; dimIdx < dim; ++dimIdx)
{
temp[0] += maxVelocity_[dimIdx] * maxVelocity_[dimIdx];
temp[1] += minVelocity_[dimIdx] * minVelocity_[dimIdx];
}
Scalar uDiff = std::sqrt(temp[0]) - std::sqrt(temp[1]);
Scalar uDiff = maxVelocity_.two_norm() - minVelocity_.two_norm();
Scalar f1 = yMax * fMax;
Scalar f2 = cWK * yMax * uDiff * uDiff / fMax;
......
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