[staggeredGrid][freeflow] Adapt nonisothermal localresidual and fluxvariables

* no test problem yet

dumux/freeflow/staggeredni/fluxvariables.hh

@@ -65,13 +65,14 @@ class FreeFlowFluxVariablesImpl<TypeTag, false, true> : public FreeFlowFluxVaria
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
     using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+    using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables); // TODO ?
     using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
     using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
     using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
     using IndexType = typename GridView::IndexSet::IndexType;
     using Stencil = std::vector<IndexType>;
 
-    using MolecularDiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
+//    using MolecularDiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
 
     static constexpr bool navierStokes = GET_PROP_VALUE(TypeTag, EnableInertiaTerms);
 //    static constexpr auto numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
@@ -108,11 +109,11 @@ public:
                                                         const SubControlVolumeFace &scvf,
                                                         const FluxVariablesCache& fluxVarsCache)
     {
-//        CellCenterPrimaryVariables flux(0.0);
-//
-//        flux += advectiveFluxForCellCenter_(problem, fvGeometry, elemVolVars, globalFaceVars, scvf);
-//        flux += MolecularDiffusionType::diffusiveFluxForCellCenter(problem, fvGeometry, elemVolVars, scvf);
-//        return flux;
+        CellCenterPrimaryVariables flux(0.0);
+
+        flux += advectiveFluxForCellCenter_(problem, fvGeometry, elemVolVars, globalFaceVars, scvf);
+        flux += diffusiveFluxForCellCenter_(problem, fvGeometry, elemVolVars, scvf);
+        return flux;
     }
 
 private:
@@ -123,53 +124,60 @@ private:
                                                           const GlobalFaceVars& globalFaceVars,
                                                           const SubControlVolumeFace &scvf)
     {
-//        CellCenterPrimaryVariables flux(0.0);
-//
-//        const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
-//        const auto& insideVolVars = elemVolVars[insideScv];
-//
-//        // if we are on an inflow/outflow boundary, use the volVars of the element itself
-//        // TODO: catch neumann and outflow in localResidual's evalBoundary_()
-//        bool isOutflow = false;
-//        if(scvf.boundary())
-//        {
-//            const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
-//                if(bcTypes.isOutflow(momentumBalanceIdx))
-//                    isOutflow = true;
-//        }
-//
-//        const auto& outsideVolVars = isOutflow ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
-//
-//        const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndex()).velocity();
-//
-//        const bool insideIsUpstream = sign(scvf.outerNormalScalar()) == sign(velocity);
-//        const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
-//        const auto& downstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
-//
-//        const Scalar upWindWeight = GET_PROP_VALUE(TypeTag, ImplicitUpwindWeight);
-//        const Scalar upstreamDensity = useMoles ? upstreamVolVars.molarDensity() : upstreamVolVars.density();
-//        const Scalar downstreamDensity = useMoles ? downstreamVolVars.molarDensity() : downstreamVolVars.density();
-//
-//        for (int compIdx = 0; compIdx < numComponents; ++compIdx)
-//        {
-//            // get equation index
-//            const auto eqIdx = conti0EqIdx + compIdx;
-//            if (eqIdx == replaceCompEqIdx)
-//                continue;
-//
-//            const Scalar upstreamFraction = useMoles ? upstreamVolVars.moleFraction(phaseIdx, compIdx) : upstreamVolVars.massFraction(phaseIdx, compIdx);
-//            const Scalar downstreamFraction = useMoles ? downstreamVolVars.moleFraction(phaseIdx, compIdx) : downstreamVolVars.massFraction(phaseIdx, compIdx);
-//
-//            flux[eqIdx] = (upWindWeight * upstreamDensity * upstreamFraction +
-//                          (1.0 - upWindWeight) * downstreamDensity * downstreamFraction)
-//                          * velocity;
-//        }
-//        // in case one balance is substituted by the total mass balance
-//        if (replaceCompEqIdx < numComponents)
-//            flux[replaceCompEqIdx] = (upWindWeight * upstreamDensity + (1.0 - upWindWeight) * downstreamDensity) * velocity;
-//
-//        flux *= scvf.area() * sign(scvf.outerNormalScalar());
-//        return flux;
+        CellCenterPrimaryVariables flux(0.0);
+
+        const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+        const auto& insideVolVars = elemVolVars[insideScv];
+
+        // if we are on an inflow/outflow boundary, use the volVars of the element itself
+        // TODO: catch neumann and outflow in localResidual's evalBoundary_()
+        bool isOutflow = false;
+        if(scvf.boundary())
+        {
+            const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
+                if(bcTypes.isOutflow(momentumBalanceIdx))
+                    isOutflow = true;
+        }
+
+        const auto& outsideVolVars = isOutflow ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
+
+        const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndex()).velocity();
+
+        const bool insideIsUpstream = sign(scvf.outerNormalScalar()) == sign(velocity);
+        const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
+        const auto& downstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
+
+        const Scalar upWindWeight = GET_PROP_VALUE(TypeTag, ImplicitUpwindWeight);
+        const Scalar upstreamDensity = useMoles ? upstreamVolVars.molarDensity() : upstreamVolVars.density();
+        const Scalar downstreamDensity = useMoles ? downstreamVolVars.molarDensity() : downstreamVolVars.density();
+        const Scalar upstreamEnthalpy = upstreamVolVars.enthalpy();
+        const Scalar downstreamEnthalpy = downstreamVolVars.enthalpy();
+
+//        flux[massBalanceIdx] = TODO??
+        flux[energyBalanceIdx] = (upWindWeight * upstreamDensity * upstreamEnthalpy
+                                 + (1.0 - upWindWeight) * downstreamDensity * downstreamEnthalpy)
+                                 * velocity;
+
+        flux *= scvf.area() * sign(scvf.outerNormalScalar());
+        return flux;
+    }
+
+    CellCenterPrimaryVariables diffusiveFluxForCellCenter_(const FluxVariables& fluxVars)
+    {
+        CellCenterPrimaryVariables flux(0.0);
+
+        // compute diffusive flux --> no diffusive flux (only 1 component)
+
+        // compute conductive flux
+        computeConductiveFlux_(flux, fluxVars);
+
+        return flux;
+    }
+
+    void computeConductiveFlux_(CellCenterPrimaryVariables& flux, FluxVariables& fluxVars)
+    {
+        flux[energyBalanceIdx] -= fluxVars.temperatureGrad() * fluxVars.face().normal
+                   * (fluxVars.thermalConductivity() + fluxVars.thermalEddyConductivity());
     }
 };
 

dumux/freeflow/staggeredni/localresidual.hh

@@ -78,10 +78,6 @@ class StaggeredNavierStokesResidualImpl<TypeTag, false, true> : public Staggered
     typename DofTypeIndices::FaceIdx faceIdx;
 
     enum { // TODO adapt
-         // grid and world dimension
-        dim = GridView::dimension,
-        dimWorld = GridView::dimensionworld,
-
         pressureIdx = Indices::pressureIdx,
         velocityIdx = Indices::velocityIdx,
 
@@ -110,105 +106,13 @@ public:
 //        const Scalar density = useMoles? volVars.molarDensity() : volVars.density();
 
         // compute storage of mass
-        storage[massBalanceIdx] = volVars.density(0);
+        storage[massBalanceIdx] = volVars.density(0); // TODO ParentType?
 
         // compute the storage of energy
         storage[energyBalanceIdx] = volVars.density(0) * volVars.internalEnergy(0);
 
         return storage;
     }
-
-    // TODO implement advectiveFlux, conductiveFlux
-
-    /*!
-     * \brief Evaluates the convective energy flux
-     *        over a face of a sub-control volume and writes the result in
-     *        the flux vector. This method is called by computeFlux in the base class.
-     *
-     * \param flux The vector for the fluxes over the SCV/boundary face
-     * \param fluxVars The flux variables at the current SCV/boundary face
-     */
-    void computeAdvectiveFlux(PrimaryVariables &flux,
-                              const FluxVariables &fluxVars) const
-    {
-        // call computation of the advective fluxes of the stokes model
-        // (momentum and mass fluxes)
-        ParentType::computeAdvectiveFlux(flux, fluxVars);
-
-        // vertex data of the upstream and the downstream vertices
-        const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx());
-        const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx());
-
-        Scalar tmp = fluxVars.normalVelocity();
-        tmp *= (this->massUpwindWeight_ * up.density() * up.enthalpy()
-                + (1.0 - this->massUpwindWeight_) * dn.density() * dn.enthalpy());
-
-        flux[energyEqIdx] += tmp;
-        Valgrind::CheckDefined(flux[energyEqIdx]);
-    }
-
-    /*!
-     * \brief Evaluates the diffusive component energy flux
-     *        over the face of a sub-control volume.
-     *
-     * \param flux The vector for the fluxes over the SCV face
-     * \param fluxVars The flux variables at the current SCV face
-     */
-    void computeDiffusiveFlux(PrimaryVariables &flux,
-                              const FluxVariables &fluxVars) const
-    {
-        // diffusive mass flux
-        ParentType::computeDiffusiveFlux(flux, fluxVars);
-
-        // conductive energy flux
-        computeConductiveFlux(flux, fluxVars);
-
-        // diffusive component energy flux
-        Scalar sumDiffusiveFluxes = 0;
-        for (int compIdx=0; compIdx<numComponents; compIdx++)
-        {
-            if (compIdx != phaseCompIdx)
-            {
-                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];
-
-        Valgrind::CheckDefined(flux[energyEqIdx]);
-    }
-
-    /*!
-     * \brief Evaluates the conductive energy flux
-     *        over the face of a sub-control volume.
-     *
-     * \param flux The vector for the fluxes over the SCV face
-     * \param fluxVars The flux variables at the current SCV face
-     */
-    void computeConductiveFlux(PrimaryVariables &flux,
-                               const FluxVariables &fluxVars) const
-    {
-        // diffusive heat flux
-        flux[energyEqIdx] -=
-            fluxVars.temperatureGrad() * fluxVars.face().normal
-            * (fluxVars.thermalConductivity() + fluxVars.thermalEddyConductivity());
-        Valgrind::CheckDefined(flux[energyEqIdx]);
-    }
 };
 
 } // end namespace