[fix] fix bugs introduces by solidsystem

dumux/porousmediumflow/nonisothermal/volumevariables.hh

@@ -135,13 +135,22 @@ public:
 
         else
         {
-            for(int phaseIdx=0; phaseIdx < numEnergyEq-1; ++phaseIdx)
+            //if numEnergyEq == 2 this means we have 1 temp for fluid phase, one for solid
+            if (numEnergyEq == 2)
             {
-                // retrieve temperatures from solution vector, phases might have different temperature
-                const Scalar T = elemSol[scv.localDofIndex()][temperatureIdx + phaseIdx];
-                fluidState.setTemperature(phaseIdx, T);
+                const Scalar T = elemSol[scv.localDofIndex()][temperatureIdx];
+                fluidState.setTemperature(T);
+            }
+            //this is for numEnergyEqFluid > 1
+            else
+            {
+                for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
+                {
+                    // retrieve temperatures from solution vector, phases might have different temperature
+                    const Scalar T = elemSol[scv.localDofIndex()][temperatureIdx + phaseIdx];
+                    fluidState.setTemperature(phaseIdx, T);
+                }
             }
-
             const Scalar solidTemperature = elemSol[scv.localDofIndex()][temperatureIdx+numEnergyEq-1];
             solidState.setTemperature(solidTemperature);
         }

test/porousmediumflow/mpnc/implicit/CMakeLists.txt

@@ -37,7 +37,7 @@ dune_add_test(COMPILE_ONLY # since it currently fails miserably with very differ
               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMD_ARGS  --script fuzzy
                         --files ${CMAKE_SOURCE_DIR}/test/references/combustion-reference.vtp
-                                ${CMAKE_CURRENT_BINARY_DIR}/test_boxmpncthermalnonequil-00084.vtp
+                                ${CMAKE_CURRENT_BINARY_DIR}/test_boxmpncthermalnonequil-00073.vtp
                         --command "${CMAKE_CURRENT_BINARY_DIR}/test_mpncthermalnonequil_box test_boxmpncthermalnonequil.input -Problem.Name test_boxmpncthermalnonequil")
 
 

test/porousmediumflow/mpnc/implicit/combustionproblem1c.hh

@@ -37,10 +37,10 @@
 
 #include <dumux/material/solidstates/compositionalsolidstate.hh>
 #include <dumux/material/solidsystems/compositionalsolidphase.hh>
+#include <dumux/material/components/constant.hh>
 
 #include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysimplefluidlumping.hh>
 #include <dumux/material/constraintsolvers/computefromreferencephase.hh>
-#include <dumux/material/components/constant.hh>
 
 #include "combustionspatialparams.hh"
 #include "combustionfluidsystem.hh"
@@ -69,8 +69,7 @@ SET_TYPE_PROP(CombustionOneComponentTypeTag,
               Problem,
               CombustionProblemOneComponent<TypeTag>);
 
-// Set the spatial parameters
-SET_TYPE_PROP(CombustionOneComponentTypeTag, SpatialParams, CombustionSpatialParams<TypeTag>);
+
 
 SET_TYPE_PROP(CombustionOneComponentTypeTag,
               FluidSystem,
@@ -116,6 +115,7 @@ SET_INT_PROP(CombustionOneComponentTypeTag, NumEnergyEqSolid, 1);
 // by default chemical non equilibrium is enabled in the nonequil model, switch that off here
 SET_BOOL_PROP(CombustionOneComponentTypeTag, EnableChemicalNonEquilibrium, false);
 //#################
+
 SET_PROP(CombustionOneComponentTypeTag, SolidSystem)
 {
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
@@ -124,6 +124,19 @@ SET_PROP(CombustionOneComponentTypeTag, SolidSystem)
     static constexpr int numInertComponents = 2;
     using type = SolidSystems::CompositionalSolidPhase<Scalar, ComponentOne, ComponentTwo, numInertComponents>;
 };
+
+SET_PROP(CombustionOneComponentTypeTag, SolidState)
+{
+private:
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using SolidSystem = typename GET_PROP_TYPE(TypeTag, SolidSystem);
+public:
+    using type = CompositionalSolidState<Scalar, SolidSystem>;
+};
+
+// Set the spatial parameters
+SET_TYPE_PROP(CombustionOneComponentTypeTag, SpatialParams, CombustionSpatialParams<TypeTag>);
+
 }
 /*!
  * \ingroup MPNCTests
@@ -327,7 +340,7 @@ public:
         fluidState.setPressure(nPhaseIdx, pn);
         fluidState.setPressure(wPhaseIdx, pw);
 
-        fluidState.setTemperature(TBoundary_ );
+        fluidState.setTemperature(TBoundary_);
         ParameterCache dummyCache;
         fluidState.setMoleFraction(wPhaseIdx, nCompIdx, 0.0);
         fluidState.setMoleFraction(wPhaseIdx, wCompIdx, 1.0);

test/porousmediumflow/mpnc/implicit/combustionspatialparams.hh

@@ -118,23 +118,45 @@ public:
      * \param elemSol The solution at the dofs connected to the element.
      * \return the porosity
      */
-    template<class SolidSystem>
-    Scalar inertVolumeFractionAtPos(const GlobalPosition& globalPos,
-                                    int compIdx) const
+    template<class ElementSolution>
+    Scalar porosity(const Element& element,
+                    const SubControlVolume& scv,
+                    const ElementSolution& elemSol) const
+    {
+        if (inOutFlow(scv.dofPosition()))
+            return porosityOutFlow_;
+        else
+            return porosity_;
+    }
+
+    /*!
+     * \brief Function for defining the porosity.
+     *        That is possibly solution dependent.
+     *
+     * \param element The current element
+     * \param scv The sub-control volume inside the element.
+     * \param elemSol The solution at the dofs connected to the element.
+     * \return the porosity
+     */
+    template<class SolidSystem, class ElementSolution>
+    Scalar inertVolumeFraction(const Element& element,
+                               const SubControlVolume& scv,
+                               const ElementSolution& elemSol,
+                               int compIdx) const
     {
         if (compIdx == SolidSystem::comp0Idx)
         {
-            if (inOutFlow(globalPos))
-                return 1-porosityOutFlow_;
+            if (inOutFlow(scv.dofPosition()))
+                return 1.0-porosityOutFlow_;
             else
                 return 0;
         }
         else
         {
-            if (inOutFlow(globalPos))
+            if (inOutFlow(scv.dofPosition()))
                 return 0;
             else
-                return 1-porosity_;
+                return 1.0-porosity_;
         }
     }
 

test/porousmediumflow/mpnc/implicit/test_boxmpncthermalnonequil.input

@@ -31,13 +31,8 @@ factorMassTransfer = 1 #
 [SpatialParams.Outflow]
 permeabilityOutFlow = 1e-6 # m^2
 porosityOutFlow = 0.35 #
-soilThermalConductivityOutFlow = 0.01#
 
 [SpatialParams.soil]
-density = 2600. # kg/m^3
-thermalConductivity = 30 # W / (m K) soil:3, metal:30
-heatCapacity = 466 # J / (kg K) steel:466 granite:800
-
 Swr =  0.0 #  5e-3
 Snr = 0 #
 
@@ -53,12 +48,14 @@ nRestart = 10000 # after so many timesteps a restart file should be written
 AddVelocity = 1 # enable velocity output
 
 [1.Component]
-SolidDensity = 2600
+SolidDensity = 2600.
 SolidThermalConductivity = 0.01
 SolidHeatCapacity = 466
 
 [2.Component]
-SolidDensity = 2600
+SolidDensity = 2600.
 SolidThermalConductivity = 30
 SolidHeatCapacity = 466
 
+[Newton]
+MaxRelativeShift = 1e-7