From 832686a36a2d996bf696be701a3e3e350e761edf Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Thu, 30 Jul 2015 09:48:31 +0000
Subject: [PATCH] [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
---
dumux/freeflow/stokes/stokesfluxvariables.hh | 45 ++++-----
dumux/freeflow/stokes/stokeslocalresidual.hh | 4 +-
.../stokesnc/stokesncfluxvariables.hh | 93 +++++++++----------
.../stokesnc/stokesnclocalresidual.hh | 51 +++++-----
.../stokesncni/stokesncnifluxvariables.hh | 32 +++----
.../stokesncni/stokesncnilocalresidual.hh | 56 +++++------
dumux/freeflow/zeroeq/zeroeqfluxvariables.hh | 8 +-
7 files changed, 127 insertions(+), 162 deletions(-)
diff --git a/dumux/freeflow/stokes/stokesfluxvariables.hh b/dumux/freeflow/stokes/stokesfluxvariables.hh
index 31648f9eb6..1d118d17f8 100644
--- a/dumux/freeflow/stokes/stokesfluxvariables.hh
+++ b/dumux/freeflow/stokes/stokesfluxvariables.hh
@@ -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
diff --git a/dumux/freeflow/stokes/stokeslocalresidual.hh b/dumux/freeflow/stokes/stokeslocalresidual.hh
index ff855415c4..6f1beddf5e 100644
--- a/dumux/freeflow/stokes/stokeslocalresidual.hh
+++ b/dumux/freeflow/stokes/stokeslocalresidual.hh
@@ -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
diff --git a/dumux/freeflow/stokesnc/stokesncfluxvariables.hh b/dumux/freeflow/stokesnc/stokesncfluxvariables.hh
index ac9ad61f7c..def0d9efe3 100644
--- a/dumux/freeflow/stokesnc/stokesncfluxvariables.hh
+++ b/dumux/freeflow/stokesnc/stokesncfluxvariables.hh
@@ -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];
diff --git a/dumux/freeflow/stokesnc/stokesnclocalresidual.hh b/dumux/freeflow/stokesnc/stokesnclocalresidual.hh
index 0843d8510e..5376704c63 100644
--- a/dumux/freeflow/stokesnc/stokesnclocalresidual.hh
+++ b/dumux/freeflow/stokesnc/stokesnclocalresidual.hh
@@ -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]);
+ }
}
-
- }
+ }
}
};
diff --git a/dumux/freeflow/stokesncni/stokesncnifluxvariables.hh b/dumux/freeflow/stokesncni/stokesncnifluxvariables.hh
index b58cdfd135..75ee28d252 100644
--- a/dumux/freeflow/stokesncni/stokesncnifluxvariables.hh
+++ b/dumux/freeflow/stokesncni/stokesncnifluxvariables.hh
@@ -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]);
}
diff --git a/dumux/freeflow/stokesncni/stokesncnilocalresidual.hh b/dumux/freeflow/stokesncni/stokesncnilocalresidual.hh
index a2913e289d..9a8d49d49e 100644
--- a/dumux/freeflow/stokesncni/stokesncnilocalresidual.hh
+++ b/dumux/freeflow/stokesncni/stokesncnilocalresidual.hh
@@ -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]);
}
};
diff --git a/dumux/freeflow/zeroeq/zeroeqfluxvariables.hh b/dumux/freeflow/zeroeq/zeroeqfluxvariables.hh
index 98029292e1..0e93c587db 100644
--- a/dumux/freeflow/zeroeq/zeroeqfluxvariables.hh
+++ b/dumux/freeflow/zeroeq/zeroeqfluxvariables.hh
@@ -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;
--
GitLab