diff --git a/dumux/freeflow/stokes/stokesfluxvariables.hh b/dumux/freeflow/stokes/stokesfluxvariables.hh index cd0937d634ad8abbaf0b14889b0e3a44e6e344eb..a688140812dcd01f64c32867c70cc2378211d77d 100644 --- a/dumux/freeflow/stokes/stokesfluxvariables.hh +++ b/dumux/freeflow/stokes/stokesfluxvariables.hh @@ -61,10 +61,8 @@ class StokesFluxVariables enum { dim = GridView::dimension }; typedef typename GridView::template Codim<0>::Entity Element; - typedef Dune::FieldVector<Scalar, dim> FieldVector; - typedef Dune::FieldVector<Scalar, dim> VelocityVector; - typedef Dune::FieldVector<Scalar, dim> ScalarGradient; - typedef Dune::FieldMatrix<Scalar, dim, dim> VectorGradient; + typedef Dune::FieldVector<Scalar, dim> DimVector; + typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix; typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; typedef typename FVElementGeometry::SubControlVolumeFace SCVFace; @@ -72,11 +70,11 @@ class StokesFluxVariables public: StokesFluxVariables(const Problem &problem, const Element &element, - const FVElementGeometry &elemGeom, - int faceIdx, + const FVElementGeometry &fvGeometry, + const int faceIdx, const ElementVolumeVariables &elemVolVars, - bool onBoundary = false) - : fvGeom_(elemGeom), onBoundary_(onBoundary), faceIdx_(faceIdx) + const bool onBoundary = false) + : fvGeometry_(fvGeometry), onBoundary_(onBoundary), faceIdx_(faceIdx) { calculateValues_(problem, element, elemVolVars); determineUpwindDirection_(elemVolVars); @@ -88,67 +86,67 @@ protected: const ElementVolumeVariables &elemVolVars) { // calculate gradients and secondary variables at IPs - FieldVector tmp(0.0); - - densityAtIP_ = Scalar(0); - molarDensityAtIP_ = Scalar(0); - viscosityAtIP_ = Scalar(0); - pressureAtIP_ = Scalar(0); - normalVelocityAtIP_ = Scalar(0); - velocityAtIP_ = Scalar(0); - pressureGradAtIP_ = Scalar(0); - velocityGradAtIP_ = Scalar(0); -// velocityDivAtIP_ = Scalar(0); + DimVector tmp(0.0); + + density_ = Scalar(0); + molarDensity_ = Scalar(0); + viscosity_ = Scalar(0); + pressure_ = Scalar(0); + normalvelocity_ = Scalar(0); + velocity_ = Scalar(0); + pressureGrad_ = Scalar(0); + velocityGrad_ = Scalar(0); +// velocityDiv_ = Scalar(0); for (int idx = 0; - idx < fvGeom_.numVertices; + idx < fvGeometry_.numVertices; idx++) // loop over adjacent vertices { // phase density and viscosity at IP - densityAtIP_ += elemVolVars[idx].density() * + density_ += elemVolVars[idx].density() * face().shapeValue[idx]; - molarDensityAtIP_ += elemVolVars[idx].molarDensity()* + molarDensity_ += elemVolVars[idx].molarDensity()* face().shapeValue[idx]; - viscosityAtIP_ += elemVolVars[idx].viscosity() * + viscosity_ += elemVolVars[idx].viscosity() * face().shapeValue[idx]; - pressureAtIP_ += elemVolVars[idx].pressure() * + pressure_ += elemVolVars[idx].pressure() * face().shapeValue[idx]; // velocity at the IP (fluxes) - VelocityVector velocityTimesShapeValue = elemVolVars[idx].velocity(); + DimVector velocityTimesShapeValue = elemVolVars[idx].velocity(); velocityTimesShapeValue *= face().shapeValue[idx]; - velocityAtIP_ += velocityTimesShapeValue; + velocity_ += velocityTimesShapeValue; // the pressure gradient tmp = face().grad[idx]; tmp *= elemVolVars[idx].pressure(); - pressureGradAtIP_ += tmp; + pressureGrad_ += tmp; // take gravity into account tmp = problem.gravity(); - tmp *= densityAtIP_; + tmp *= density_; // pressure gradient including influence of gravity - pressureGradAtIP_ -= tmp; + pressureGrad_ -= tmp; // the velocity gradients and divergence for (int dimIdx = 0; dimIdx<dim; ++dimIdx) { tmp = face().grad[idx]; tmp *= elemVolVars[idx].velocity()[dimIdx]; - velocityGradAtIP_[dimIdx] += tmp; + velocityGrad_[dimIdx] += tmp; -// velocityDivAtIP_ += face().grad[idx][dimIdx]*elemVolVars[idx].velocity()[dimIdx]; +// velocityDiv_ += face().grad[idx][dimIdx]*elemVolVars[idx].velocity()[dimIdx]; } } - normalVelocityAtIP_ = velocityAtIP_ * face().normal; + normalvelocity_ = velocity_ * face().normal; - Valgrind::CheckDefined(densityAtIP_); - Valgrind::CheckDefined(viscosityAtIP_); - Valgrind::CheckDefined(normalVelocityAtIP_); - Valgrind::CheckDefined(velocityAtIP_); - Valgrind::CheckDefined(pressureGradAtIP_); - Valgrind::CheckDefined(velocityGradAtIP_); -// Valgrind::CheckDefined(velocityDivAtIP_); + Valgrind::CheckDefined(density_); + Valgrind::CheckDefined(viscosity_); + Valgrind::CheckDefined(normalvelocity_); + Valgrind::CheckDefined(velocity_); + Valgrind::CheckDefined(pressureGrad_); + Valgrind::CheckDefined(velocityGrad_); +// Valgrind::CheckDefined(velocityDiv_); }; void determineUpwindDirection_(const ElementVolumeVariables &elemVolVars) @@ -158,7 +156,7 @@ protected: upstreamIdx_ = face().i; downstreamIdx_ = face().j; - if (normalVelocityAtIP() < 0) + if (normalVelocity() < 0) std::swap(upstreamIdx_, downstreamIdx_); }; @@ -170,9 +168,9 @@ public: const SCVFace &face() const { if (onBoundary_) - return fvGeom_.boundaryFace[faceIdx_]; + return fvGeometry_.boundaryFace[faceIdx_]; else - return fvGeom_.subContVolFace[faceIdx_]; + return fvGeometry_.subContVolFace[faceIdx_]; } /*! @@ -181,69 +179,130 @@ public: */ const Scalar averageSCVVolume() const { - return 0.5*(fvGeom_.subContVol[upstreamIdx_].volume + - fvGeom_.subContVol[downstreamIdx_].volume); + return 0.5*(fvGeometry_.subContVol[upstreamIdx_].volume + + fvGeometry_.subContVol[downstreamIdx_].volume); } /*! * \brief Return the pressure \f$\mathrm{[Pa]}\f$ at the integration * point. */ + Scalar pressure() const + { return pressure_; } + + /*! + * \brief Return the pressure \f$\mathrm{[Pa]}\f$ at the integration + * point. + */ + DUMUX_DEPRECATED_MSG("use pressure() instead") Scalar pressureAtIP() const - { return pressureAtIP_; } + { return pressure(); } + + /*! + * \brief Return the mass density \f$ \mathrm{[kg/m^3]} \f$ at the integration + * point. + */ + Scalar density() const + { return density_; } /*! * \brief Return the mass density \f$ \mathrm{[kg/m^3]} \f$ at the integration * point. */ + DUMUX_DEPRECATED_MSG("use density() instead") Scalar densityAtIP() const - { return densityAtIP_; } + { return density(); } /*! * \brief Return the molar density \f$ \mathrm{[mol/m^3]} \f$ at the integration point. */ + const Scalar molarDensity() const + { return molarDensity_; } + + /*! + * \brief Return the molar density \f$ \mathrm{[mol/m^3]} \f$ at the integration point. + */ + DUMUX_DEPRECATED_MSG("use molarDensity() instead") const Scalar molarDensityAtIP() const - { return molarDensityAtIP_; } + { return molarDensity(); } + + /*! + * \brief Return the viscosity \f$ \mathrm{[m^2/s]} \f$ at the integration + * point. + */ + Scalar viscosity() const + { return viscosity_; } /*! * \brief Return the viscosity \f$ \mathrm{[m^2/s]} \f$ at the integration * point. */ + DUMUX_DEPRECATED_MSG("use viscosity() instead") Scalar viscosityAtIP() const - { return viscosityAtIP_; } + { return viscosity(); } /*! * \brief Return the velocity \f$ \mathrm{[m/s]} \f$ at the integration * point multiplied by the normal and the area. */ + Scalar normalVelocity() const + { return normalvelocity_; } + + /*! + * \brief Return the velocity \f$ \mathrm{[m/s]} \f$ at the integration + * point multiplied by the normal and the area. + */ + DUMUX_DEPRECATED_MSG("use normalVelocity() instead") Scalar normalVelocityAtIP() const - { return normalVelocityAtIP_; } + { return normalVelocity(); } + + /*! + * \brief Return the pressure gradient at the integration point. + */ + const DimVector &pressureGrad() const + { return pressureGrad_; } /*! * \brief Return the pressure gradient at the integration point. */ - const ScalarGradient &pressureGradAtIP() const - { return pressureGradAtIP_; } + DUMUX_DEPRECATED_MSG("use pressureGrad() instead") + const DimVector &pressureGradAtIP() const + { return pressureGrad(); } /*! * \brief Return the velocity vector at the integration point. */ - const VelocityVector &velocityAtIP() const - { return velocityAtIP_; } + const DimVector &velocity() const + { return velocity_; } + + /*! + * \brief Return the velocity vector at the integration point. + */ + DUMUX_DEPRECATED_MSG("use velocity() instead") + const DimVector &velocityAtIP() const + { return velocity(); } + + /*! + * \brief Return the velocity gradient at the integration + * point of a face. + */ + const DimMatrix &velocityGrad() const + { return velocityGrad_; } /*! * \brief Return the velocity gradient at the integration * point of a face. */ - const VectorGradient &velocityGradAtIP() const - { return velocityGradAtIP_; } + DUMUX_DEPRECATED_MSG("use velocityGrad() instead") + const DimMatrix &velocityGradAtIP() const + { return velocityGrad(); } // /*! // * \brief Return the divergence of the normal velocity at the // * integration point. // */ -// Scalar velocityDivAtIP() const -// { return velocityDivAtIP_; } +// Scalar velocityDiv() const +// { return velocityDiv_; } /*! * \brief Return the local index of the upstream sub-control volume. @@ -265,21 +324,21 @@ public: { return onBoundary_; } protected: - const FVElementGeometry &fvGeom_; + const FVElementGeometry &fvGeometry_; const bool onBoundary_; // values at the integration point - Scalar densityAtIP_; - Scalar molarDensityAtIP_; - Scalar viscosityAtIP_; - Scalar pressureAtIP_; - Scalar normalVelocityAtIP_; -// Scalar velocityDivAtIP_; - VelocityVector velocityAtIP_; + Scalar density_; + Scalar molarDensity_; + Scalar viscosity_; + Scalar pressure_; + Scalar normalvelocity_; +// Scalar velocityDiv_; + DimVector velocity_; // gradients at the IPs - ScalarGradient pressureGradAtIP_; - VectorGradient velocityGradAtIP_; + DimVector pressureGrad_; + DimMatrix velocityGrad_; // local index of the upwind vertex int upstreamIdx_; diff --git a/dumux/freeflow/stokes/stokeslocaljacobian.hh b/dumux/freeflow/stokes/stokeslocaljacobian.hh index f934c83756f472a02da9e45dcd5c261daf48dba8..7db0bfd2b80d02a2fff6012937fbe15869867831 100644 --- a/dumux/freeflow/stokes/stokeslocaljacobian.hh +++ b/dumux/freeflow/stokes/stokeslocaljacobian.hh @@ -49,8 +49,8 @@ class StokesLocalJacobian : public BoxLocalJacobian<TypeTag> public: //! \copydoc BoxLocalJacobian::numericEpsilon() - Scalar numericEpsilon(int scvIdx, - int pvIdx) const + Scalar numericEpsilon(const int scvIdx, + const int pvIdx) const { Scalar pv = this->curVolVars_[scvIdx].primaryVars()[pvIdx]; if (pvIdx < GridView::dimension){ diff --git a/dumux/freeflow/stokes/stokeslocalresidual.hh b/dumux/freeflow/stokes/stokeslocalresidual.hh index 5daef4876d3b12c3675e7434066732d5aeba0aa1..01a0cfed342bbf7777cfa34a7e6b4e5f30ff78f6 100644 --- a/dumux/freeflow/stokes/stokeslocalresidual.hh +++ b/dumux/freeflow/stokes/stokeslocalresidual.hh @@ -80,8 +80,7 @@ protected: typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElements; typedef Dune::GenericReferenceElement<Scalar, dim> ReferenceElement; - typedef Dune::FieldVector<Scalar, dim> ScalarGradient; - typedef Dune::FieldVector<Scalar, dim> FieldVector; + typedef Dune::FieldVector<Scalar, dim> DimVector; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables; @@ -115,11 +114,11 @@ protected: * The result should be averaged over the volume (e.g. phase mass * inside a sub control volume divided by the volume) * - * \param result The mass of the component within the sub-control volume + * \param storage The mass of the component within the sub-control volume * \param scvIdx The SCV (sub-control-volume) index * \param usePrevSol Evaluate function with solution of current or previous time step */ - void computeStorage(PrimaryVariables &result, int scvIdx, bool usePrevSol) const + void computeStorage(PrimaryVariables &storage, const int scvIdx, const bool usePrevSol) const { // if flag usePrevSol is set, the solution from the previous // time step is used, otherwise the current solution is @@ -128,17 +127,17 @@ protected: // using the implicit Euler method. const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_() : this->curVolVars_(); - const VolumeVariables &vertexData = elemVolVars[scvIdx]; + const VolumeVariables &volVars = elemVolVars[scvIdx]; - result = 0.0; + storage = 0.0; // mass balance - result[massBalanceIdx] = vertexData.density(); + storage[massBalanceIdx] = volVars.density(); // momentum balance for (int momentumIdx = momentumXIdx; momentumIdx <= lastMomentumIdx; ++momentumIdx) - result[momentumIdx] = vertexData.density() - * vertexData.velocity()[momentumIdx-momentumXIdx]; + storage[momentumIdx] = volVars.density() + * volVars.velocity()[momentumIdx-momentumXIdx]; } /*! @@ -153,7 +152,7 @@ protected: * the created fluxVars object contains boundary variables evaluated at the IP of the * boundary face */ - void computeFlux(PrimaryVariables &flux, int faceIdx, bool onBoundary=false) const + void computeFlux(PrimaryVariables &flux, const int faceIdx, const bool onBoundary=false) const { const FluxVariables fluxVars(this->problem_(), this->element_(), @@ -190,7 +189,7 @@ protected: const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx()); // mass balance with upwinded density - FieldVector massBalanceResidual = fluxVars.velocityAtIP(); + DimVector massBalanceResidual = fluxVars.velocity(); if (massUpwindWeight_ == 1.0) // fully upwind massBalanceResidual *= up.density(); else @@ -201,7 +200,7 @@ protected: { // stabilization of the mass balance // with 0.5*alpha*(V_i + V_j)*grad P - FieldVector stabilizationTerm = fluxVars.pressureGradAtIP(); + DimVector stabilizationTerm = fluxVars.pressureGrad(); stabilizationTerm *= stabilizationAlpha_* fluxVars.averageSCVVolume(); massBalanceResidual += stabilizationTerm; @@ -216,26 +215,26 @@ protected: // compute symmetrized gradient for the momentum flux: // mu (grad v + (grad v)^t) - Dune::FieldMatrix<Scalar, dim, dim> symmVelGrad = fluxVars.velocityGradAtIP(); + Dune::FieldMatrix<Scalar, dim, dim> symmVelGrad = fluxVars.velocityGrad(); for (int i=0; i<dim; ++i) for (int j=0; j<dim; ++j) - symmVelGrad[i][j] += fluxVars.velocityGradAtIP()[j][i]; + symmVelGrad[i][j] += fluxVars.velocityGrad()[j][i]; - FieldVector velGradComp(0.); + DimVector velGradComp(0.); for (int velIdx = 0; velIdx < dim; ++velIdx) { velGradComp = symmVelGrad[velIdx]; // TODO: dilatation term has to be accounted for in outflow, coupling, neumann - // velGradComp[velIdx] += 2./3*fluxVars.velocityDivAtIP; + // velGradComp[velIdx] += 2./3*fluxVars.velocityDiv; - velGradComp *= fluxVars.viscosityAtIP(); + velGradComp *= fluxVars.viscosity(); flux[momentumXIdx + velIdx] -= velGradComp*fluxVars.face().normal; // gravity is accounted for in computeSource; alternatively: - // Scalar gravityTerm = fluxVars.densityAtIP * + // Scalar gravityTerm = fluxVars.density * // this->problem_().gravity()[dim-1] * // fluxVars.face().ipGlobal[dim-1]* // fluxVars.face().normal[velIdx]; @@ -251,7 +250,7 @@ protected: { for (int dimIndex = 0; dimIndex < dim; ++dimIndex) { - flux[momentumXIdx + dimIndex] += up.density() * up.velocity()[dimIndex] * fluxVars.normalVelocityAtIP(); + flux[momentumXIdx + dimIndex] += up.density() * up.velocity()[dimIndex] * fluxVars.normalVelocity(); } } } @@ -276,35 +275,35 @@ protected: * The pressure gradient at the center of a SCV is computed * and the gravity term evaluated. * - * \param q The source/sink in the sub control volume for each component - * \param localVertexIdx The local index of the sub-control volume + * \param source The source/sink in the sub control volume for each component + * \param scvIndex The local index of the sub-control volume */ - void computeSource(PrimaryVariables &q, int localVertexIdx) + void computeSource(PrimaryVariables &source, const int scvIdx) { const ElementVolumeVariables &elemVolVars = this->curVolVars_(); - const VolumeVariables &vertexData = elemVolVars[localVertexIdx]; + const VolumeVariables &volVars = elemVolVars[scvIdx]; // retrieve the source term intrinsic to the problem - this->problem_().source(q, + this->problem_().source(source, this->element_(), this->fvGeometry_(), - localVertexIdx); + scvIdx); // ATTENTION: The source term of the mass balance has to be chosen as // div (q_momentum) in the problem file const Scalar alphaH2 = stabilizationAlpha_* - this->fvGeometry_().subContVol[localVertexIdx].volume; - q[massBalanceIdx] *= alphaH2; // stabilization of the source term + this->fvGeometry_().subContVol[scvIdx].volume; + source[massBalanceIdx] *= alphaH2; // stabilization of the source term // pressure gradient at the center of the SCV, // the pressure is discretized as volume term, // while -mu grad v is calculated in computeFlux - ScalarGradient pressureGradAtSCVCenter(0.0); - ScalarGradient grad(0.0); + DimVector pressureGradAtSCVCenter(0.0); + DimVector grad(0.0); for (int vertexIdx = 0; vertexIdx < this->fvGeometry_().numVertices; vertexIdx++) { - grad = this->fvGeometry_().subContVol[localVertexIdx].gradCenter[vertexIdx]; + grad = this->fvGeometry_().subContVol[scvIdx].gradCenter[vertexIdx]; Valgrind::CheckDefined(grad); grad *= elemVolVars[vertexIdx].pressure(); @@ -316,8 +315,8 @@ protected: // signs are inverted, since q is subtracted for (int dimIdx=0; dimIdx<dim; ++dimIdx) { - q[momentumXIdx + dimIdx] -= pressureGradAtSCVCenter[dimIdx]; - q[momentumXIdx + dimIdx] += vertexData.density()*this->problem_().gravity()[dimIdx]; + source[momentumXIdx + dimIdx] -= pressureGradAtSCVCenter[dimIdx]; + source[momentumXIdx + dimIdx] += volVars.density()*this->problem_().gravity()[dimIdx]; } } @@ -338,8 +337,8 @@ protected: continue; // important at corners of the grid - FieldVector momentumResidual(0.0); - FieldVector averagedNormal(0.0); + DimVector momentumResidual(0.0); + DimVector averagedNormal(0.0); int numberOfOuterFaces = 0; // evaluate boundary conditions for the intersections of // the current element @@ -379,14 +378,14 @@ protected: // the fluxes at the boundary are added in the second step if (momentumBalanceDirichlet_(bcTypes)) { - FieldVector muGradVelNormal(0.); - const FieldVector &boundaryFaceNormal = + DimVector muGradVelNormal(0.); + const DimVector &boundaryFaceNormal = boundaryVars.face().normal; - boundaryVars.velocityGradAtIP().umv(boundaryFaceNormal, muGradVelNormal); - muGradVelNormal *= boundaryVars.viscosityAtIP(); + boundaryVars.velocityGrad().umv(boundaryFaceNormal, muGradVelNormal); + muGradVelNormal *= boundaryVars.viscosity(); - for (int i=0; i < this->residual_.size(); i++) + for (unsigned int i=0; i < this->residual_.size(); i++) Valgrind::CheckDefined(this->residual_[i]); for (int dimIdx=0; dimIdx < dim; ++dimIdx) momentumResidual[dimIdx] = this->residual_[vertexIdx][momentumXIdx+dimIdx]; @@ -450,7 +449,7 @@ protected: // mathematically Neumann BC: p n - mu grad v n = q // boundary terms: -mu grad v n // implement q * A (from evalBoundarySegment) - p n(unity) A - FieldVector pressureCorrection(boundaryVars.face().normal); + DimVector pressureCorrection(boundaryVars.face().normal); pressureCorrection *= this->curVolVars_(scvIdx).pressure(); for (int momentumIdx = momentumXIdx; momentumIdx <= lastMomentumIdx; momentumIdx++) if(bcTypes.isNeumann(momentumIdx)) @@ -460,15 +459,15 @@ protected: // in case of neumann conditions for the momentum equation; // calculate mu grad v t t // center in the face of the reference element - FieldVector tangent; + DimVector tangent; if (stabilizationBeta_ != 0) { - const FieldVector& elementUnitNormal = isIt->centerUnitOuterNormal(); + const DimVector& elementUnitNormal = isIt->centerUnitOuterNormal(); tangent[0] = elementUnitNormal[1]; //TODO: 3D tangent[1] = -elementUnitNormal[0]; - FieldVector tangentialVelGrad; - boundaryVars.velocityGradAtIP().mv(tangent, tangentialVelGrad); - tangentialVelGrad *= boundaryVars.viscosityAtIP(); + DimVector tangentialVelGrad; + boundaryVars.velocityGrad().mv(tangent, tangentialVelGrad); + tangentialVelGrad *= boundaryVars.viscosity(); this->residual_[scvIdx][massBalanceIdx] -= stabilizationBeta_*0.5* this->curVolVars_(scvIdx).pressure(); @@ -520,15 +519,15 @@ protected: { // calculate mu grad v t t for beta-stabilization // center in the face of the reference element - FieldVector tangent; - const FieldVector& elementUnitNormal = isIt->centerUnitOuterNormal(); + DimVector tangent; + const DimVector& elementUnitNormal = isIt->centerUnitOuterNormal(); tangent[0] = elementUnitNormal[1]; tangent[1] = -elementUnitNormal[0]; - FieldVector tangentialVelGrad; - boundaryVars.velocityGradAtIP().mv(tangent, tangentialVelGrad); + DimVector tangentialVelGrad; + boundaryVars.velocityGrad().mv(tangent, tangentialVelGrad); this->residual_[scvIdx][massBalanceIdx] -= 0.5*stabilizationBeta_ - * boundaryVars.viscosityAtIP() + * boundaryVars.viscosity() * (tangentialVelGrad*tangent); } } @@ -558,7 +557,7 @@ protected: const Scalar alphaH2 = stabilizationAlpha_* fluxVars.averageSCVVolume(); - Scalar stabilizationTerm = fluxVars.pressureGradAtIP() * + Scalar stabilizationTerm = fluxVars.pressureGrad() * this->fvGeometry_().subContVolFace[faceIdx].normal; stabilizationTerm *= alphaH2; @@ -610,8 +609,8 @@ protected: * \brief Replace the local residual of the mass balance equation by * the sum of the residuals of the momentum balance equation. */ - void replaceMassbalanceResidual_(const FieldVector& momentumResidual, - FieldVector& averagedNormal, + void replaceMassbalanceResidual_(const DimVector& momentumResidual, + DimVector& averagedNormal, const int vertexIdx) { assert(averagedNormal.two_norm() != 0.0); diff --git a/dumux/freeflow/stokes/stokesmodel.hh b/dumux/freeflow/stokes/stokesmodel.hh index a831b420746d08bdfa9a0947d81ce02fc7e655f7..a85505926e3d5a37da4589814d19d537ca3dcdbe 100644 --- a/dumux/freeflow/stokes/stokesmodel.hh +++ b/dumux/freeflow/stokes/stokesmodel.hh @@ -105,7 +105,7 @@ public: PrimaryVariables tmpFlux(0.0); int sign; - FVElementGeometry fvElemGeom; + FVElementGeometry fvGeometry; ElementVolumeVariables elemVolVars; // Loop over elements @@ -116,14 +116,14 @@ public: if (elemIt->partitionType() != Dune::InteriorEntity) continue; - fvElemGeom.update(this->gridView_(), *elemIt); - elemVolVars.update(this->problem_(), *elemIt, fvElemGeom); + fvGeometry.update(this->gridView_(), *elemIt); + elemVolVars.update(this->problem_(), *elemIt, fvGeometry); this->localResidual().evalFluxes(*elemIt, elemVolVars); bool hasLeft = false; bool hasRight = false; - for (int i = 0; i < fvElemGeom.numVertices; i++) { - const GlobalPosition &global = fvElemGeom.subContVol[i].global; + for (int i = 0; i < fvGeometry.numVertices; i++) { + const GlobalPosition &global = fvGeometry.subContVol[i].global; if (globalI[axis] < coordVal) hasLeft = true; else if (globalI[axis] >= coordVal) @@ -132,10 +132,10 @@ public: if (!hasLeft || !hasRight) continue; - for (int i = 0; i < fvElemGeom.numVertices; i++) { + for (int i = 0; i < fvGeometry.numVertices; i++) { const int &globalIdx = this->vertexMapper().map(*elemIt, i, dim); - const GlobalPosition &global = fvElemGeom.subContVol[i].global; + const GlobalPosition &global = fvGeometry.subContVol[i].global; if (globalI[axis] < coordVal) flux += this->localResidual().residual(i); } @@ -163,7 +163,7 @@ public: unsigned numElements = this->gridView_().size(0); ScalarField &rank = *writer.allocateManagedBuffer(numElements); - FVElementGeometry fvElemGeom; + FVElementGeometry fvGeometry; VolumeVariables volVars; ElementBoundaryTypes elemBcTypes; @@ -174,7 +174,7 @@ public: int idx = this->elementMapper().map(*elemIt); rank[idx] = this->gridView_().comm().rank(); - fvElemGeom.update(this->gridView_(), *elemIt); + fvGeometry.update(this->gridView_(), *elemIt); elemBcTypes.update(this->problem_(), *elemIt); int numLocalVerts = elemIt->template count<dim>(); @@ -184,7 +184,7 @@ public: volVars.update(sol[globalIdx], this->problem_(), *elemIt, - fvElemGeom, + fvGeometry, i, false); diff --git a/dumux/freeflow/stokes/stokesproblem.hh b/dumux/freeflow/stokes/stokesproblem.hh index 0f86e23a1c889310db4301f05f033bdd4c2afe8b..312f2bd488fc4ded2a185de6792767436c38cdf6 100644 --- a/dumux/freeflow/stokes/stokesproblem.hh +++ b/dumux/freeflow/stokes/stokesproblem.hh @@ -100,10 +100,10 @@ public: * \return (Scalar) permeability */ Scalar permeability(const Element &element, - const FVElementGeometry &fvElemGeom, + const FVElementGeometry &fvGeometry, const Intersection &is, - int scvIdx, - int boundaryFaceIdx) const + const int scvIdx, + const int boundaryFaceIdx) const { DUNE_THROW(Dune::NotImplemented, "permeability()"); } // \} diff --git a/dumux/freeflow/stokes/stokesvolumevariables.hh b/dumux/freeflow/stokes/stokesvolumevariables.hh index b0ddca47cc6611e242010807094ad3a55cc20c89..d6f6d000b2554f6721172945ce36d9b85a7344f0 100644 --- a/dumux/freeflow/stokes/stokesvolumevariables.hh +++ b/dumux/freeflow/stokes/stokesvolumevariables.hh @@ -70,7 +70,7 @@ class StokesVolumeVariables : public BoxVolumeVariables<TypeTag> typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState; - typedef Dune::FieldVector<Scalar, dim> VelocityVector; + typedef Dune::FieldVector<Scalar, dim> DimVector; public: /*! @@ -79,18 +79,18 @@ public: void update(const PrimaryVariables &priVars, const Problem &problem, const Element &element, - const FVElementGeometry &elemGeom, - int scvIdx, - bool isOldSol) + const FVElementGeometry &fvGeometry, + const int scvIdx, + const bool isOldSol) { ParentType::update(priVars, problem, element, - elemGeom, + fvGeometry, scvIdx, isOldSol); - completeFluidState(priVars, problem, element, elemGeom, scvIdx, fluidState_, isOldSol); + completeFluidState(priVars, problem, element, fvGeometry, scvIdx, fluidState_, isOldSol); for (int dimIdx=momentumXIdx; dimIdx<=lastMomentumIdx; ++dimIdx) velocity_[dimIdx] = priVars[dimIdx]; @@ -100,18 +100,18 @@ public: * \copydoc BoxModel::completeFluidState() * \param isOldSol Specifies whether this is the previous solution or the current one */ - static void completeFluidState(const PrimaryVariables& primaryVariables, + static void completeFluidState(const PrimaryVariables& priVars, const Problem& problem, const Element& element, - const FVElementGeometry& elementGeometry, - int scvIdx, + const FVElementGeometry& fvGeometry, + const int scvIdx, FluidState& fluidState, - bool isOldSol = false) + const bool isOldSol = false) { - Scalar temperature = Implementation::temperature_(primaryVariables, problem, - element, elementGeometry, scvIdx); + Scalar temperature = Implementation::temperature_(priVars, problem, + element, fvGeometry, scvIdx); fluidState.setTemperature(temperature); - fluidState.setPressure(phaseIdx, primaryVariables[pressureIdx]); + fluidState.setPressure(phaseIdx, priVars[pressureIdx]); // create NullParameterCache and do dummy update typename FluidSystem::ParameterCache paramCache; @@ -178,7 +178,7 @@ public: /*! * \brief Returns the velocity vector in the sub-control volume. */ - const VelocityVector &velocity() const + const DimVector &velocity() const { return velocity_; } protected: @@ -190,16 +190,16 @@ protected: return 0; } - static Scalar temperature_(const PrimaryVariables &primaryVars, + static Scalar temperature_(const PrimaryVariables &priVars, const Problem& problem, const Element &element, - const FVElementGeometry &elemGeom, - int scvIdx) + const FVElementGeometry &fvGeometry, + const int scvIdx) { - return problem.boxTemperature(element, elemGeom, scvIdx); + return problem.boxTemperature(element, fvGeometry, scvIdx); } - VelocityVector velocity_; + DimVector velocity_; FluidState fluidState_; private: diff --git a/dumux/freeflow/stokes2c/stokes2cfluxvariables.hh b/dumux/freeflow/stokes2c/stokes2cfluxvariables.hh index a457d1a35129b4d12ccb30e995743ab150b12313..dcb254fb0708f4bd9abe80e29bb584d755f9071c 100644 --- a/dumux/freeflow/stokes2c/stokes2cfluxvariables.hh +++ b/dumux/freeflow/stokes2c/stokes2cfluxvariables.hh @@ -66,21 +66,20 @@ class Stokes2cFluxVariables : public StokesFluxVariables<TypeTag> typedef typename GET_PROP_TYPE(TypeTag, Stokes2cIndices) Indices; enum { dim = GridView::dimension }; - enum { lCompIdx = Indices::lCompIdx }; + enum { comp1Idx = Indices::comp1Idx }; enum { phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIndex) }; typedef typename GridView::template Codim<0>::Entity Element; - typedef Dune::FieldVector<Scalar, dim> ScalarGradient; - + typedef Dune::FieldVector<Scalar, dim> DimVector; public: Stokes2cFluxVariables(const Problem &problem, const Element &element, - const FVElementGeometry &elemGeom, - int faceIdx, + const FVElementGeometry &fvGeometry, + const int faceIdx, const ElementVolumeVariables &elemVolVars, - bool onBoundary = false) - : ParentType(problem, element, elemGeom, faceIdx, elemVolVars, onBoundary) + const bool onBoundary = false) + : ParentType(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary) { calculateValues_(problem, element, elemVolVars); } @@ -88,20 +87,41 @@ public: /*! * \brief Return the mass fraction at the integration point. */ + Scalar massFraction() const + { return massFraction_; } + + /*! + * \brief Return the mass fraction at the integration point. + */ + DUMUX_DEPRECATED_MSG("use massFraction() instead") Scalar massFractionAtIP() const - { return massFractionAtIP_; } + { return massFraction(); } /*! * \brief Return the molar diffusion coefficient at the integration point. */ + Scalar diffusionCoeff() const + { return diffusionCoeff_; } + + /*! + * \brief Return the molar diffusion coefficient at the integration point. + */ + DUMUX_DEPRECATED_MSG("use diffusionCoeff() instead") Scalar diffusionCoeffAtIP() const - { return diffusionCoeffAtIP_; } + { return diffusionCoeff(); } + + /*! + * \brief Return the gradient of the mole fraction at the integration point. + */ + const DimVector &moleFractionGrad() const + { return moleFractionGrad_; } /*! * \brief Return the gradient of the mole fraction at the integration point. */ - const ScalarGradient &moleFractionGradAtIP() const - { return moleFractionGradAtIP_; } + DUMUX_DEPRECATED_MSG("use moleFractionGrad() instead") + const DimVector &moleFractionGradAtIP() const + { return moleFractionGrad(); } protected: @@ -109,37 +129,37 @@ protected: const Element &element, const ElementVolumeVariables &elemVolVars) { - massFractionAtIP_ = Scalar(0); - diffusionCoeffAtIP_ = Scalar(0); - moleFractionGradAtIP_ = Scalar(0); + massFraction_ = Scalar(0); + diffusionCoeff_ = Scalar(0); + moleFractionGrad_ = Scalar(0); // calculate gradients and secondary variables at IPs for (int idx = 0; - idx < this->fvGeom_.numVertices; + idx < this->fvGeometry_.numVertices; idx++) // loop over vertices of the element { - massFractionAtIP_ += elemVolVars[idx].fluidState().massFraction(phaseIdx, lCompIdx) * + massFraction_ += elemVolVars[idx].fluidState().massFraction(phaseIdx, comp1Idx) * this->face().shapeValue[idx]; - diffusionCoeffAtIP_ += elemVolVars[idx].diffusionCoeff() * + diffusionCoeff_ += elemVolVars[idx].diffusionCoeff() * this->face().shapeValue[idx]; // the gradient of the mass fraction at the IP for (int dimIdx=0; dimIdx<dim; ++dimIdx) { - moleFractionGradAtIP_ += + moleFractionGrad_ += this->face().grad[idx][dimIdx] * - elemVolVars[idx].fluidState().moleFraction(phaseIdx, lCompIdx); + elemVolVars[idx].fluidState().moleFraction(phaseIdx, comp1Idx); } }; - Valgrind::CheckDefined(massFractionAtIP_); - Valgrind::CheckDefined(diffusionCoeffAtIP_); - Valgrind::CheckDefined(moleFractionGradAtIP_); + Valgrind::CheckDefined(massFraction_); + Valgrind::CheckDefined(diffusionCoeff_); + Valgrind::CheckDefined(moleFractionGrad_); } - Scalar massFractionAtIP_; - Scalar diffusionCoeffAtIP_; - ScalarGradient moleFractionGradAtIP_; + Scalar massFraction_; + Scalar diffusionCoeff_; + DimVector moleFractionGrad_; }; } // end namespace diff --git a/dumux/freeflow/stokes2c/stokes2cindices.hh b/dumux/freeflow/stokes2c/stokes2cindices.hh index 1fd631b909fc3edc680c5b51a234526aaa4deede..61e144b764b7b8d9646a1f760333cb1962d499d0 100644 --- a/dumux/freeflow/stokes2c/stokes2cindices.hh +++ b/dumux/freeflow/stokes2c/stokes2cindices.hh @@ -49,19 +49,19 @@ struct Stokes2cCommonIndices : public StokesCommonIndices<TypeTag> public: // Phase indices - static const int lPhaseIdx = FluidSystem::lPhaseIdx; //!< Index of the liquid phase - static const int gPhaseIdx = FluidSystem::gPhaseIdx; //!< Index of the gas phase + static const int wPhaseIdx = FluidSystem::lPhaseIdx; //!< Index of the wetting phase + static const int nPhaseIdx = FluidSystem::gPhaseIdx; //!< Index of the non-wetting phase // Component indices - static const int lCompIdx = 0; //!< Index of the liquid's primary component - static const int gCompIdx = 1; //!< Index of the gas' primary component + static const int comp1Idx = 0; //!< Index of the wetting's primary component + static const int comp0Idx = 1; //!< Index of the non-wetting's primary component // equation and primary variable indices static const int dim = StokesCommonIndices<TypeTag>::dim; static const int transportIdx = PVOffset + dim+1; //! The index for the transport equation. static const int massOrMoleFracIndex = transportIdx; //! The index for the mass or mole fraction in primary variable vectors. - static const int phaseIdx = gPhaseIdx; //!< Index of the gas phase (required to use the same fluid system in coupled models) + static const int phaseIdx = nPhaseIdx; //!< Index of the non-wetting phase (required to use the same fluid system in coupled models) }; } // end namespace diff --git a/dumux/freeflow/stokes2c/stokes2clocalresidual.hh b/dumux/freeflow/stokes2c/stokes2clocalresidual.hh index f5c9d3d80d4e5a9bdce73bcdb1a4fe12713ba8f3..38828b5a0fa82ffdc880c6210c3dfd1417f6c2d8 100644 --- a/dumux/freeflow/stokes2c/stokes2clocalresidual.hh +++ b/dumux/freeflow/stokes2c/stokes2clocalresidual.hh @@ -54,7 +54,7 @@ class Stokes2cLocalResidual : public StokesLocalResidual<TypeTag> enum { dim = GridView::dimension }; enum { transportIdx = Indices::transportIdx }; //!< Index of the transport equation - enum { lCompIdx = Indices::lCompIdx }; //!< Index of the liquid component + enum { comp1Idx = Indices::comp1Idx }; //!< Index of the transported component enum { phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIndex)}; //!< Index of the considered phase (only of interest when using two-phase fluidsystems) typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; @@ -71,30 +71,30 @@ public: * The result should be averaged over the volume (e.g. phase mass * inside a sub control volume divided by the volume) * - * \param result The mass of the component within the sub-control volume + * \param storage The mass of the component within the sub-control volume * \param scvIdx The SCV (sub-control-volume) index * \param usePrevSol Evaluate function with solution of current or previous time step */ - void computeStorage(PrimaryVariables &result, int scvIdx, bool usePrevSol) const + void computeStorage(PrimaryVariables &storage, const int scvIdx, const bool usePrevSol) const { // compute the storage term for the transport equation - ParentType::computeStorage(result, scvIdx, usePrevSol); + ParentType::computeStorage(storage, scvIdx, usePrevSol); // if flag usePrevSol is set, the solution from the previous // time step is used, otherwise the current solution is // used. The secondary variables are used accordingly. This // is required to compute the derivative of the storage term // using the implicit euler method. - const ElementVolumeVariables &elemDat = usePrevSol ? this->prevVolVars_() : this->curVolVars_(); - const VolumeVariables &vertexDat = elemDat[scvIdx]; + const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_() : this->curVolVars_(); + const VolumeVariables &volVars = elemVolVars[scvIdx]; // compute the storage of the component - result[transportIdx] = - vertexDat.density() * - vertexDat.fluidState().massFraction(phaseIdx, lCompIdx); + storage[transportIdx] = + volVars.density() * + volVars.fluidState().massFraction(phaseIdx, comp1Idx); - Valgrind::CheckDefined(vertexDat.density()); - Valgrind::CheckDefined(vertexDat.fluidState().massFraction(phaseIdx, lCompIdx)); + Valgrind::CheckDefined(volVars.density()); + Valgrind::CheckDefined(volVars.fluidState().massFraction(phaseIdx, comp1Idx)); } /*! @@ -118,14 +118,14 @@ public: const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx()); const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx()); - Scalar tmp = fluxVars.normalVelocityAtIP(); + Scalar tmp = fluxVars.normalVelocity(); if (this->massUpwindWeight_ > 0.0) tmp *= this->massUpwindWeight_ * // upwind data - up.density() * up.fluidState().massFraction(phaseIdx, lCompIdx); + up.density() * up.fluidState().massFraction(phaseIdx, comp1Idx); if (this->massUpwindWeight_ < 1.0) tmp += (1.0 - this->massUpwindWeight_) * // rest - dn.density() * dn.fluidState().massFraction(phaseIdx, lCompIdx); + dn.density() * dn.fluidState().massFraction(phaseIdx, comp1Idx); flux[transportIdx] += tmp; Valgrind::CheckDefined(flux[transportIdx]); @@ -147,11 +147,11 @@ public: // diffusive component flux for (int dimIdx = 0; dimIdx < dim; ++dimIdx) flux[transportIdx] -= - fluxVars.moleFractionGradAtIP()[dimIdx] * + fluxVars.moleFractionGrad()[dimIdx] * fluxVars.face().normal[dimIdx] * - fluxVars.diffusionCoeffAtIP() * - fluxVars.molarDensityAtIP() * - FluidSystem::molarMass(lCompIdx); + fluxVars.diffusionCoeff() * + fluxVars.molarDensity() * + FluidSystem::molarMass(comp1Idx); Valgrind::CheckDefined(flux[transportIdx]); } diff --git a/dumux/freeflow/stokes2c/stokes2cmodel.hh b/dumux/freeflow/stokes2c/stokes2cmodel.hh index 3caeceec3efbcf68cbecc136a5dce33c4b8fe584..0bcaa270ee891213950e58189e7b192bcf90d79d 100644 --- a/dumux/freeflow/stokes2c/stokes2cmodel.hh +++ b/dumux/freeflow/stokes2c/stokes2cmodel.hh @@ -75,7 +75,7 @@ class Stokes2cModel : public StokesModel<TypeTag> typedef typename GET_PROP_TYPE(TypeTag, Stokes2cIndices) Indices; enum { dim = GridView::dimension }; - enum { lCompIdx = Indices::lCompIdx }; + enum { comp1Idx = Indices::comp1Idx }; enum { phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIndex) }; typedef typename GridView::template Codim<0>::Iterator ElementIterator; @@ -137,7 +137,7 @@ public: pN[globalIdx] = volVars.pressure()*scale_; delP[globalIdx] = volVars.pressure()*scale_ - 1e5; - Xw[globalIdx] = volVars.fluidState().massFraction(phaseIdx, lCompIdx); + Xw[globalIdx] = volVars.fluidState().massFraction(phaseIdx, comp1Idx); rho[globalIdx] = volVars.density()*scale_*scale_*scale_; mu[globalIdx] = volVars.viscosity()*scale_; velocity[globalIdx] = volVars.velocity(); @@ -147,7 +147,7 @@ public: writer.attachVertexData(pN, "P"); writer.attachVertexData(delP, "delP"); std::ostringstream outputNameX; - outputNameX << "X^" << FluidSystem::componentName(lCompIdx); + outputNameX << "X^" << FluidSystem::componentName(comp1Idx); writer.attachVertexData(Xw, outputNameX.str()); writer.attachVertexData(rho, "rho"); writer.attachVertexData(mu, "mu"); diff --git a/dumux/freeflow/stokes2c/stokes2cpropertydefaults.hh b/dumux/freeflow/stokes2c/stokes2cpropertydefaults.hh index 820e9e2c84e95055fa9556489fbeca4d4ca17b92..2775725bd402a2f73a4c59f9d06e151e007d92ac 100644 --- a/dumux/freeflow/stokes2c/stokes2cpropertydefaults.hh +++ b/dumux/freeflow/stokes2c/stokes2cpropertydefaults.hh @@ -96,7 +96,7 @@ SET_PROP(BoxStokes2c, PhaseIndex) { typedef typename GET_PROP_TYPE(TypeTag, Stokes2cIndices) Indices; public: - static constexpr int value = Indices::gPhaseIdx; + static constexpr int value = Indices::nPhaseIdx; }; } diff --git a/dumux/freeflow/stokes2c/stokes2cvolumevariables.hh b/dumux/freeflow/stokes2c/stokes2cvolumevariables.hh index c4bdbf2663658d296c1c99dabec1f43674338162..25d8853ee270424a4a5f6e6f1e953de5a2e3dbf7 100644 --- a/dumux/freeflow/stokes2c/stokes2cvolumevariables.hh +++ b/dumux/freeflow/stokes2c/stokes2cvolumevariables.hh @@ -57,8 +57,8 @@ class Stokes2cVolumeVariables : public StokesVolumeVariables<TypeTag> typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; enum { - lCompIdx = Indices::lCompIdx, - gCompIdx = Indices::gCompIdx + comp1Idx = Indices::comp1Idx, + comp0Idx = Indices::comp0Idx }; enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) }; enum { phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIndex) }; @@ -71,18 +71,18 @@ public: void update(const PrimaryVariables &priVars, const Problem &problem, const Element &element, - const FVElementGeometry &elemGeom, - int scvIdx, - bool isOldSol) + const FVElementGeometry &fvGeometry, + const int scvIdx, + const bool isOldSol) { // set the mole fractions first - completeFluidState(priVars, problem, element, elemGeom, scvIdx, this->fluidState(), isOldSol); + completeFluidState(priVars, problem, element, fvGeometry, scvIdx, this->fluidState(), isOldSol); // update vertex data for the mass and momentum balance ParentType::update(priVars, problem, element, - elemGeom, + fvGeometry, scvIdx, isOldSol); @@ -95,8 +95,8 @@ public: diffCoeff_ = FluidSystem::binaryDiffusionCoefficient(this->fluidState(), paramCache, phaseIdx, - lCompIdx, - gCompIdx); + comp1Idx, + comp0Idx); Valgrind::CheckDefined(diffCoeff_); }; @@ -105,27 +105,27 @@ public: * \copydoc BoxModel::completeFluidState() * \param isOldSol Specifies whether this is the previous solution or the current one */ - static void completeFluidState(const PrimaryVariables& primaryVariables, + static void completeFluidState(const PrimaryVariables& priVars, const Problem& problem, const Element& element, - const FVElementGeometry& elementGeometry, - int scvIdx, + const FVElementGeometry& fvGeometry, + const int scvIdx, FluidState& fluidState, - bool isOldSol = false) + const bool isOldSol = false) { Scalar massFraction[numComponents]; - massFraction[lCompIdx] = primaryVariables[transportIdx]; - massFraction[gCompIdx] = 1 - massFraction[lCompIdx]; + massFraction[comp1Idx] = priVars[transportIdx]; + massFraction[comp0Idx] = 1 - massFraction[comp1Idx]; // calculate average molar mass of the gas phase - Scalar M1 = FluidSystem::molarMass(lCompIdx); - Scalar M2 = FluidSystem::molarMass(gCompIdx); - Scalar X2 = massFraction[gCompIdx]; + Scalar M1 = FluidSystem::molarMass(comp1Idx); + Scalar M2 = FluidSystem::molarMass(comp0Idx); + Scalar X2 = massFraction[comp0Idx]; Scalar avgMolarMass = M1*M2/(M2 + X2*(M1 - M2)); // convert mass to mole fractions and set the fluid state - fluidState.setMoleFraction(phaseIdx, lCompIdx, massFraction[lCompIdx]*avgMolarMass/M1); - fluidState.setMoleFraction(phaseIdx, gCompIdx, massFraction[gCompIdx]*avgMolarMass/M2); + fluidState.setMoleFraction(phaseIdx, comp1Idx, massFraction[comp1Idx]*avgMolarMass/M1); + fluidState.setMoleFraction(phaseIdx, comp0Idx, massFraction[comp0Idx]*avgMolarMass/M2); } /*! diff --git a/dumux/freeflow/stokes2cni/stokes2cnifluxvariables.hh b/dumux/freeflow/stokes2cni/stokes2cnifluxvariables.hh index 40a519af1f96cf02405b5025593e9dce6040816f..aab4770b727c35018f0145c9c9ddfc72f867b95a 100644 --- a/dumux/freeflow/stokes2cni/stokes2cnifluxvariables.hh +++ b/dumux/freeflow/stokes2cni/stokes2cnifluxvariables.hh @@ -61,18 +61,18 @@ class Stokes2cniFluxVariables : public Stokes2cFluxVariables<TypeTag> enum { dim = GridView::dimension }; typedef typename GridView::template Codim<0>::Entity Element; - typedef Dune::FieldVector<Scalar, dim> ScalarGradient; + typedef Dune::FieldVector<Scalar, dim> DimVector; typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; public: Stokes2cniFluxVariables(const Problem &problem, const Element &element, - const FVElementGeometry &elemGeom, - int faceIdx, + const FVElementGeometry &fvGeometry, + const int faceIdx, const ElementVolumeVariables &elemVolVars, - bool onBoundary = false) - : ParentType(problem, element, elemGeom, faceIdx, elemVolVars, onBoundary) + const bool onBoundary = false) + : ParentType(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary) { calculateValues_(problem, element, elemVolVars); } @@ -80,44 +80,58 @@ public: /*! * \brief Returns the heat conductivity at the integration point. */ + Scalar heatConductivity() const + { return heatConductivity_; } + + /*! + * \brief Returns the heat conductivity at the integration point. + */ + DUMUX_DEPRECATED_MSG("use heatConductivity() instead") Scalar heatConductivityAtIP() const - { return heatConductivityAtIP_; } + { return heatConductivity(); } + + /*! + * \brief Returns the temperature gradient at the integration point. + */ + const DimVector &temperatureGrad() const + { return temperatureGrad_; } /*! * \brief Returns the temperature gradient at the integration point. */ - const ScalarGradient &temperatureGradAtIP() const - { return temperatureGradAtIP_; } + DUMUX_DEPRECATED_MSG("use temperatureGrad() instead") + const DimVector &temperatureGradAtIP() const + { return temperatureGrad(); } protected: void calculateValues_(const Problem &problem, const Element &element, const ElementVolumeVariables &elemVolVars) { - heatConductivityAtIP_ = Scalar(0); - temperatureGradAtIP_ = Scalar(0); + heatConductivity_ = Scalar(0); + temperatureGrad_ = Scalar(0); // calculate gradients and secondary variables at IPs - ScalarGradient tmp(0.0); + DimVector tmp(0.0); for (int idx = 0; - idx < this->fvGeom_.numVertices; + idx < this->fvGeometry_.numVertices; idx++) // loop over vertices of the element { - heatConductivityAtIP_ += elemVolVars[idx].heatConductivity() * + heatConductivity_ += elemVolVars[idx].heatConductivity() * this->face().shapeValue[idx]; // the gradient of the temperature at the IP for (int dimIdx=0; dimIdx<dim; ++dimIdx) - temperatureGradAtIP_ += + temperatureGrad_ += this->face().grad[idx][dimIdx]* elemVolVars[idx].temperature(); } - Valgrind::CheckDefined(heatConductivityAtIP_); - Valgrind::CheckDefined(temperatureGradAtIP_); + Valgrind::CheckDefined(heatConductivity_); + Valgrind::CheckDefined(temperatureGrad_); } - Scalar heatConductivityAtIP_; - ScalarGradient temperatureGradAtIP_; + Scalar heatConductivity_; + DimVector temperatureGrad_; }; } // end namespace diff --git a/dumux/freeflow/stokes2cni/stokes2cnilocalresidual.hh b/dumux/freeflow/stokes2cni/stokes2cnilocalresidual.hh index 9e668461cd5e691f1c1944c1d7db8aac48c17ae8..7aea7a619d07a11b8826ba838872485006dc56ec 100644 --- a/dumux/freeflow/stokes2cni/stokes2cnilocalresidual.hh +++ b/dumux/freeflow/stokes2cni/stokes2cnilocalresidual.hh @@ -43,10 +43,6 @@ namespace Dumux * \brief Element-wise calculation of the Jacobian matrix for problems * using the non-isothermal compositional Stokes box model. This class is derived * from the stokes2c box local residual and adds the energy balance equation. - * - * \param result The mass of the component within the sub-control volume - * \param scvIdx The SCV (sub-control-volume) index - * \param usePrevSol Evaluate function with solution of current or previous time step */ template<class TypeTag> class Stokes2cniLocalResidual : public Stokes2cLocalResidual<TypeTag> @@ -74,23 +70,23 @@ public: * The result should be averaged over the volume (e.g. phase mass * inside a sub control volume divided by the volume) */ - void computeStorage(PrimaryVariables &result, int scvIdx, bool usePrevSol) const + void computeStorage(PrimaryVariables &storage, const int scvIdx, const bool usePrevSol) const { // compute the storage term for the transport equation - ParentType::computeStorage(result, scvIdx, usePrevSol); + ParentType::computeStorage(storage, scvIdx, usePrevSol); // if flag usePrevSol is set, the solution from the previous // time step is used, otherwise the current solution is // used. The secondary variables are used accordingly. This // is required to compute the derivative of the storage term // using the implicit euler method. - const ElementVolumeVariables &elemDat = usePrevSol ? this->prevVolVars_() : this->curVolVars_(); - const VolumeVariables &vertexDat = elemDat[scvIdx]; + const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_() : this->curVolVars_(); + const VolumeVariables &volVars = elemVolVars[scvIdx]; // compute the storage of energy - result[energyIdx] = - vertexDat.density() * - vertexDat.internalEnergy(); + storage[energyIdx] = + volVars.density() * + volVars.internalEnergy(); } /*! @@ -112,7 +108,7 @@ public: const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx()); const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx()); - Scalar tmp = fluxVars.normalVelocityAtIP(); + Scalar tmp = fluxVars.normalVelocity(); tmp *= this->massUpwindWeight_ * // upwind data up.density() * up.enthalpy() + @@ -139,9 +135,9 @@ public: // diffusive heat flux for (int dimIdx = 0; dimIdx < dim; ++dimIdx) flux[energyIdx] -= - fluxVars.temperatureGradAtIP()[dimIdx] * + fluxVars.temperatureGrad()[dimIdx] * fluxVars.face().normal[dimIdx] * - fluxVars.heatConductivityAtIP(); + fluxVars.heatConductivity(); } }; diff --git a/dumux/freeflow/stokes2cni/stokes2cnimodel.hh b/dumux/freeflow/stokes2cni/stokes2cnimodel.hh index ebe36ca07103e1e7f65e971017ce956d12790976..894a07773dc86049dc45446b52471d5360795898 100644 --- a/dumux/freeflow/stokes2cni/stokes2cnimodel.hh +++ b/dumux/freeflow/stokes2cni/stokes2cnimodel.hh @@ -85,7 +85,7 @@ class Stokes2cniModel : public Stokes2cModel<TypeTag> typedef typename GET_PROP_TYPE(TypeTag, Stokes2cniIndices) Indices; enum { dim = GridView::dimension }; - enum { lCompIdx = Indices::lCompIdx }; + enum { comp1Idx = Indices::comp1Idx }; enum { phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIndex) }; typedef typename GridView::template Codim<0>::Iterator ElementIterator; @@ -148,7 +148,7 @@ public: pN [globalIdx] = volVars.pressure(); delP[globalIdx] = volVars.pressure() - 1e5; - Xw [globalIdx] = volVars.fluidState().massFraction(phaseIdx, lCompIdx); + Xw [globalIdx] = volVars.fluidState().massFraction(phaseIdx, comp1Idx); T [globalIdx] = volVars.temperature(); rho [globalIdx] = volVars.density(); mu [globalIdx] = volVars.viscosity(); @@ -161,7 +161,7 @@ public: writer.attachVertexData(delP, "delP"); // writer.attachVertexData(D, "Dwg"); std::ostringstream outputNameX; - outputNameX << "X^" << FluidSystem::componentName(lCompIdx); + outputNameX << "X^" << FluidSystem::componentName(comp1Idx); writer.attachVertexData(Xw, outputNameX.str()); writer.attachVertexData(T, "temperature"); writer.attachVertexData(rho, "rhoG"); diff --git a/dumux/freeflow/stokes2cni/stokes2cnivolumevariables.hh b/dumux/freeflow/stokes2cni/stokes2cnivolumevariables.hh index 02c854a656f10596ad5fd13a250eb2408b77405b..89044c9a1178f5d85a844accd231acf6f0f405cf 100644 --- a/dumux/freeflow/stokes2cni/stokes2cnivolumevariables.hh +++ b/dumux/freeflow/stokes2cni/stokes2cnivolumevariables.hh @@ -69,7 +69,7 @@ public: const Problem &problem, const Element &element, const FVElementGeometry &elemGeom, - int scvIdx, + const int scvIdx, bool isOldSol) { // the internal energies and the enthalpies @@ -116,7 +116,7 @@ protected: const Problem& problem, const Element &element, const FVElementGeometry &elemGeom, - int scvIdx) + const int scvIdx) { return priVars[energyIdx]; } @@ -124,7 +124,7 @@ protected: template<class ParameterCache> static Scalar enthalpy_(const FluidState& fluidState, const ParameterCache& paramCache, - int phaseIdx) + const int phaseIdx) { return FluidSystem::enthalpy(fluidState, paramCache, phaseIdx); }