From 26b7c31ef01eccef66cf7822a48a34b566116139 Mon Sep 17 00:00:00 2001
From: Kilian <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Wed, 25 Mar 2020 18:08:48 +0100
Subject: [PATCH] [nonequilibrium] Cleanup localresidual

* use computeFlux of parent class
---
 .../nonequilibrium/localresidual.hh           | 95 ++++---------------
 1 file changed, 17 insertions(+), 78 deletions(-)

diff --git a/dumux/porousmediumflow/nonequilibrium/localresidual.hh b/dumux/porousmediumflow/nonequilibrium/localresidual.hh
index e0f26d871c..707d5cdae0 100644
--- a/dumux/porousmediumflow/nonequilibrium/localresidual.hh
+++ b/dumux/porousmediumflow/nonequilibrium/localresidual.hh
@@ -41,99 +41,35 @@ using NonEquilibriumLocalResidual = NonEquilibriumLocalResidualImplementation<Ty
 
 /*!
  * \ingroup NonEquilibriumModel
- * \brief The mass conservation part of the nonequilibrium model for a model without chemical non-equilibrium
+ * \brief The local residual for a model without chemical non-equilibrium
+ *        but potentially with thermal non-equilibrium
  */
 template<class TypeTag>
 class NonEquilibriumLocalResidualImplementation<TypeTag, false>: public GetPropType<TypeTag, Properties::EquilibriumLocalResidual>
 {
-    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using ParentType = GetPropType<TypeTag, Properties::EquilibriumLocalResidual>;
     using Problem = GetPropType<TypeTag, Properties::Problem>;
     using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
     using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
-    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
-    using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>;
-    using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
     using GridView = GetPropType<TypeTag, Properties::GridView>;
     using Element = typename GridView::template Codim<0>::Entity;
     using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
     using EnergyLocalResidual = GetPropType<TypeTag, Properties::EnergyLocalResidual>;
     using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
-    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
 
-    static constexpr int numPhases = ModelTraits::numFluidPhases();
-    static constexpr int numComponents = ModelTraits::numFluidComponents();
-    static constexpr auto conti0EqIdx = ModelTraits::Indices::conti0EqIdx;
 public:
     using ParentType::ParentType;
 
     /*!
      * \brief Calculates the source term of the equation.
      *
-     * \param problem The object specifying the problem which ought to be simulated
+     * \param problem The source term
      * \param element An element which contains part of the control volume
      * \param fvGeometry The finite-volume geometry
      * \param elemVolVars The volume variables of the current element
-     * \param scvf The sub control volume face to compute the flux on
-     * \param elemFluxVarsCache The cache related to flux compuation
-     *
-     * \note This is the default implementation for all models as sources are computed
-     *       in the user interface of the problem
-     *
+     * \param scv The sub-control volume over which we integrate the source term
      */
-
-    NumEqVector computeFlux(const Problem& problem,
-                            const Element& element,
-                            const FVElementGeometry& fvGeometry,
-                            const ElementVolumeVariables& elemVolVars,
-                            const SubControlVolumeFace& scvf,
-                            const ElementFluxVariablesCache& elemFluxVarsCache) const
-    {
-        FluxVariables fluxVars;
-        fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
-        NumEqVector flux(0.0);
-
-        const auto moleDensity = [](const auto& volVars, const int phaseIdx)
-        { return volVars.molarDensity(phaseIdx); };
-
-        const auto moleFraction= [](const auto& volVars, const int phaseIdx, const int compIdx)
-        { return  volVars.moleFraction(phaseIdx, compIdx); };
-
-        // advective fluxes
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-        {
-            const auto diffusiveFluxes = fluxVars.molecularDiffusionFlux(phaseIdx);
-            for (int compIdx = 0; compIdx < numComponents; ++compIdx)
-            {
-                // get equation index
-                const auto eqIdx = conti0EqIdx + compIdx;
-
-                // the physical quantities for which we perform upwinding
-                const auto upwindTerm = [&moleDensity, &moleFraction, phaseIdx, compIdx] (const auto& volVars)
-                { return moleDensity(volVars, phaseIdx)*moleFraction(volVars, phaseIdx, compIdx)*volVars.mobility(phaseIdx); };
-
-                    flux[eqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
-
-                // diffusive fluxes (only for the component balances)
-                //check for the reference system and adapt units of the diffusive flux accordingly.
-                if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
-                    flux[eqIdx] += diffusiveFluxes[compIdx]/FluidSystem::molarMass(compIdx);
-                else
-                    flux[eqIdx] += diffusiveFluxes[compIdx];
-            }
-
-            //! Add advective and diffusive phase energy fluxes
-            EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx);
-        }
-        //! Add diffusive energy fluxes. For isothermal model the contribution is zero.
-        EnergyLocalResidual::heatConductionFlux(flux, fluxVars);
-
-        return flux;
-
-
-    }
-
     NumEqVector computeSource(const Problem& problem,
                               const Element& element,
                               const FVElementGeometry& fvGeometry,
@@ -147,21 +83,26 @@ public:
 
         // add contribution from possible point sources
         source += problem.scvPointSources(element, fvGeometry, elemVolVars, scv);
+
         // Call the (kinetic) Energy module, for the source term.
         // it has to be called from here, because the mass transfered has to be known.
-        // TODO: check for thermal non-eq ?
-        EnergyLocalResidual::computeSourceEnergy(source,
-                                                 element,
-                                                 fvGeometry,
-                                                 elemVolVars,
-                                                 scv);
+        if constexpr(ModelTraits::enableThermalNonEquilibrium())
+        {
+            EnergyLocalResidual::computeSourceEnergy(source,
+                                                     element,
+                                                     fvGeometry,
+                                                     elemVolVars,
+                                                     scv);
+        }
+
         return source;
     }
 
 };
 
 /*!
- * \brief The mass conservation part of the nonequilibrium model for a model assuming chemical non-equilibrium
+ * \brief The local residual for a model assuming chemical non-equilibrium
+ *        and potentially thermal non-equilibrium
  */
 template<class TypeTag>
 class NonEquilibriumLocalResidualImplementation<TypeTag, true>: public GetPropType<TypeTag, Properties::EquilibriumLocalResidual>
@@ -187,10 +128,8 @@ class NonEquilibriumLocalResidualImplementation<TypeTag, true>: public GetPropTy
 
     static constexpr int numPhases = ModelTraits::numFluidPhases();
     static constexpr int numComponents = ModelTraits::numFluidComponents();
-    static constexpr bool useMoles = ModelTraits::useMoles();
 
     static constexpr auto conti0EqIdx = Indices::conti0EqIdx;
-    static constexpr auto comp1Idx = FluidSystem::comp1Idx;
     static constexpr auto comp0Idx = FluidSystem::comp0Idx;
     static constexpr auto phase0Idx = FluidSystem::phase0Idx;
     static constexpr auto phase1Idx = FluidSystem::phase1Idx;
@@ -199,7 +138,7 @@ class NonEquilibriumLocalResidualImplementation<TypeTag, true>: public GetPropTy
                   "chemical non-equlibrium only makes sense for multiple phases");
     static_assert(numPhases == numComponents,
                   "currently chemical non-equilibrium is only available when numPhases equals numComponents");
-    static_assert(useMoles == true,
+    static_assert(ModelTraits::useMoles(),
                   "chemical nonequilibrium can only be calculated based on mole fractions not mass fractions");
 
 public:
-- 
GitLab