stokesncvolumevariables: introduced convenience functions for mass and mole fractions

stokesnclocalresidual: storage term is not set to zero here, this would lead to problems for a derived model model (e.g. nonisothermal stokesnc);
convenience functions are used for mass and mole fractions; small cleanup.

stokesncfluxvariables: use of convenience functions

Reviewed by Markus

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@12138 2fb0f335-1f38-0410-981e-8018bf24f1b0

dumux/freeflow/stokesnc/stokesncfluxvariables.hh

@@ -87,7 +87,7 @@ public:
      * \brief Return the molar density \f$ \mathrm{[mol/m^3]} \f$ at the integration point.
      */
     const Scalar molarDensity() const
-    { return this->molarDensity_; }
+    { return molarDensity_; }
 	
 	/*!
      * \brief Return the mass fraction of a transported component at the integration point.
@@ -121,9 +121,9 @@ protected:
 
 		// loop over all components
 		for (int compIdx=0; compIdx<numComponents; compIdx++){
-			if (phaseCompIdx!=compIdx) //no transport equationen parameters needed for the mass balance
+			if (phaseCompIdx!=compIdx) //no transport equation parameters needed for the mass balance
 			{
-				this->molarDensity_ = Scalar(0.0);  				
+				molarDensity_ = Scalar(0.0);
 				massFraction_[compIdx] = Scalar(0.0);
 				diffusionCoeff_[compIdx] = Scalar(0.0);
 				moleFractionGrad_[compIdx] = Scalar(0.0);
@@ -134,10 +134,9 @@ protected:
 					 scvIdx++) // loop over vertices of the element
 				{
             
-					this->molarDensity_ += elemVolVars[scvIdx].molarDensity()*
+					molarDensity_ += elemVolVars[scvIdx].molarDensity()*
 						this->face().shapeValue[scvIdx];
-					
-					massFraction_[compIdx] += elemVolVars[scvIdx].fluidState().massFraction(phaseIdx, compIdx) *
+					massFraction_[compIdx] += elemVolVars[scvIdx].massFraction(compIdx) *
 						this->face().shapeValue[scvIdx];
 					diffusionCoeff_[compIdx] += elemVolVars[scvIdx].diffusionCoeff(compIdx) *
 						this->face().shapeValue[scvIdx];
@@ -147,10 +146,11 @@ protected:
 					{
 						moleFractionGrad_[compIdx] +=
 							this->face().grad[scvIdx][dimIdx] *
-							elemVolVars[scvIdx].fluidState().moleFraction(phaseIdx, compIdx);
+							elemVolVars[scvIdx].moleFraction(compIdx);
 					}
 				}
 							
+				Valgrind::CheckDefined(molarDensity_);
 				Valgrind::CheckDefined(massFraction_[compIdx]);
 				Valgrind::CheckDefined(diffusionCoeff_[compIdx]);
 				Valgrind::CheckDefined(moleFractionGrad_[compIdx]);

dumux/freeflow/stokesnc/stokesnclocalresidual.hh

@@ -100,57 +100,44 @@ public:
      */
     void computeStorage(PrimaryVariables &storage, const int scvIdx, const bool usePrevSol) const
     {
-       
-		// if flag usePrevSol is set, the solution from the previous
+        // compute the storage term for the transport equation
+        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 &elemVolVars = usePrevSol ? this->prevVolVars_()
-		: this->curVolVars_();
+        const ElementVolumeVariables &elemVolVars = usePrevSol ?
+        			this->prevVolVars_() : this->curVolVars_();
         const VolumeVariables &volVars = elemVolVars[scvIdx];
 		
-        storage = 0.0;
         if (!useMoles)
 		{
 			/* works for a maximum of two components, for more components 
 			mole fractions must be used (property useMoles => true) */
 			
-			// mass and momentum balance
-			ParentType::computeStorage(storage, scvIdx, usePrevSol);
-			
 			//storage of transported component
 			storage[transportEqIdx] = volVars.density() 
-				* volVars.fluidState().massFraction(phaseIdx, transportCompIdx);
+				* volVars.massFraction(transportCompIdx);
 			
 			Valgrind::CheckDefined(volVars.density());
-			Valgrind::CheckDefined(volVars.fluidState().massFraction(phaseIdx, transportCompIdx));
+			Valgrind::CheckDefined(volVars.massFraction(transportCompIdx));
 			
 		}
         else
 		{
-        	/*//TODO call parent type function
-            // momentum balance
-            for (int momentumIdx = momentumXIdx; momentumIdx <= lastMomentumIdx; ++momentumIdx)
-                storage[momentumIdx] = volVars.molarDensity()
-					* volVars.velocity()[momentumIdx-momentumXIdx];*/
-			
-            // mass and momentum balance
-			ParentType::computeStorage(storage, scvIdx, usePrevSol);
-
-            
 			// mass balance and transport equations
 			for (int compIdx=0; compIdx<numComponents; compIdx++)
 			{
 				if (conti0EqIdx+compIdx != massBalanceIdx)
-					//storage[massBalanceIdx] = volVars.molarDensity();
 				//else // transport equations
 				{
 					storage[conti0EqIdx+compIdx] = volVars.molarDensity()
-						* volVars.fluidState().moleFraction(phaseIdx, compIdx);
+						* volVars.moleFraction(compIdx);
 					
 					Valgrind::CheckDefined(volVars.molarDensity());
-					Valgrind::CheckDefined(volVars.fluidState().moleFraction(phaseIdx, compIdx));
+					Valgrind::CheckDefined(volVars.moleFraction(compIdx));
 				}
 			}
 		}
@@ -169,53 +156,46 @@ public:
     void computeAdvectiveFlux(PrimaryVariables &flux,
                               const FluxVariables &fluxVars) const
     {
+      	// call ParentType function
+		ParentType::computeAdvectiveFlux(flux,fluxVars);
+            
         // data attached to upstream and the downstream vertices
 		const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx());
 		const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx());
-		Scalar tmp = 0.0;
+
+		Scalar tmp = fluxVars.normalVelocity();  
 		
 		if(!useMoles)
 		{
-        	// call ParentType function
-			ParentType::computeAdvectiveFlux(flux,fluxVars);
-			
 			// for transport equations
-			tmp = fluxVars.normalVelocity();  
-			
 			if (this->massUpwindWeight_ > 0.0)
 				tmp *=  this->massUpwindWeight_          // upwind data
 					* up.density()
-                    * up.fluidState().massFraction(phaseIdx, transportCompIdx);
+                    * up.massFraction(transportCompIdx);
 			
             if (this->massUpwindWeight_ < 1.0)
 				tmp += (1.0 - this->massUpwindWeight_)      // rest
 					* dn.density()
-                    * dn.fluidState().massFraction(phaseIdx, transportCompIdx);
+                    * dn.massFraction(transportCompIdx);
 			
 			flux[transportEqIdx] += tmp;
 			Valgrind::CheckDefined(flux[transportEqIdx]);
-			
 		}
         else
 		{
-        	// call ParentType function
-			ParentType::computeAdvectiveFlux(flux,fluxVars);
-            
 			//transport equations
 			for (int compIdx=0; compIdx<numComponents; compIdx++)
 			{
 				if (conti0EqIdx+compIdx != massBalanceIdx) //mass balance is calculated above
 				{
-					tmp = fluxVars.normalVelocity();
-				
 					if (this->massUpwindWeight_ > 0.0)
 						tmp *=  this->massUpwindWeight_          // upwind data
                             * up.molarDensity()
-                            * up.fluidState().moleFraction(phaseIdx, compIdx);
+                            * up.moleFraction(compIdx);
 					if (this->massUpwindWeight_ < 1.0)
 						tmp += (1.0 - this->massUpwindWeight_)      // rest
                             * dn.molarDensity()
-                            * dn.fluidState().moleFraction(phaseIdx, compIdx);
+                            * dn.moleFraction(compIdx);
 				
 					flux[conti0EqIdx+compIdx] += tmp;
 					Valgrind::CheckDefined(flux[conti0EqIdx+compIdx]);
@@ -236,15 +216,15 @@ public:
                               const FluxVariables &fluxVars) const
     {
 		// diffusive component flux
-        for (int dimIdx = 0; dimIdx < dim; ++dimIdx){
-
+        for (int dimIdx = 0; dimIdx < dim; ++dimIdx)
+        {
 			if(!useMoles)
 			{
                 flux[transportEqIdx] -= fluxVars.moleFractionGrad(transportCompIdx)[dimIdx]
 		                          * fluxVars.face().normal[dimIdx]
 		                          *(fluxVars.diffusionCoeff(transportCompIdx) + fluxVars.eddyDiffusivity())
 		                          * fluxVars.molarDensity()
-		                          * FluidSystem::molarMass(transportCompIdx);// Multipled by molarMass [kg/mol] to convert form [mol/m^3 s] to [kg/m^3 s]
+		                          * 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

dumux/freeflow/stokesnc/stokesncpropertydefaults.hh

@@ -49,7 +49,7 @@ namespace Properties
 // Properties
 //////////////////////////////////////////////////////////////////
 
-//!< Define that mass fractions are used in the balance equations
+//!< Define that mole fractions are used in the balance equations
 SET_BOOL_PROP(BoxStokesnc, UseMoles, true);
 		
 SET_PROP(BoxStokesnc, NumEq) //!< set the number of equations

dumux/freeflow/stokesnc/stokesncvolumevariables.hh

@@ -105,17 +105,20 @@ public:
         typename FluidSystem::ParameterCache paramCache;
         paramCache.updateAll(this->fluidState());
 		
-		for (int compIdx=0; compIdx<numComponents; compIdx++){
-			if (phaseCompIdx!=compIdx){
-				
+        for (int compIdx=0; compIdx<numComponents; compIdx++)
+        {
+			if (phaseCompIdx!=compIdx)
+			{
 				diffCoeff_[compIdx] = FluidSystem::binaryDiffusionCoefficient(this->fluidState(),
                                                              paramCache,
                                                              phaseIdx,
                                                              compIdx,
                                                              phaseCompIdx);
-
-				Valgrind::CheckDefined(diffCoeff_[compIdx]);
 			}
+			else
+				diffCoeff_[compIdx] = 0.0;
+
+			Valgrind::CheckDefined(diffCoeff_[compIdx]);
 		}
     };
 
@@ -163,19 +166,34 @@ public:
 			
 			//molefraction for the main component (no primary variable)
 			moleFracPhase = 1 - sumMoleFrac;
-			fluidState.setMoleFraction(phaseIdx, phaseIdx, moleFracPhase);
+			fluidState.setMoleFraction(phaseIdx, phaseCompIdx, moleFracPhase);
 		}
     }
 
+    /*!
+     * \brief Returns the mass fraction of a given component in the
+     * 		  given fluid phase within the control volume.
+     *
+     * \param compIdx The component index
+     */
+    Scalar massFraction(const int compIdx) const
+    { return this->fluidState_.massFraction(phaseIdx, compIdx); }
+
+    /*!
+     * \brief Returns the mass fraction of a given component in the
+     * 		  given fluid phase within the control volume.
+     *
+     * \param compIdx The component index
+     */
+    Scalar moleFraction(const int compIdx) const
+    { return this->fluidState_.moleFraction(phaseIdx, compIdx); }
+
     /*!
      * \brief Returns the molar density \f$\mathrm{[mol/m^3]}\f$ of the fluid within the
      *        sub-control volume.
      */
     Scalar molarDensity() const
-    { 
-		return this->fluidState_.density(phaseIdx) / this->fluidState_.averageMolarMass(phaseIdx); 
-		
-	}
+    { return this->fluidState_.density(phaseIdx) / this->fluidState_.averageMolarMass(phaseIdx); }
 
     /*!
      * \brief Returns the binary (mass) diffusion coefficient