From 469f0a1909c24b4b00b909bf5cfbe4aa7fd9b61b Mon Sep 17 00:00:00 2001
From: Thomas Fetzer <thomas.fetzer@iws.uni-stuttgart.de>
Date: Tue, 15 Dec 2015 15:51:42 +0100
Subject: [PATCH] [multidomain,localoperator] Further cleanup and restructure
 of the localoperators.

This time:
  - deprecated old coupling12() and coupling21() and replaced them by one coupling()
  - try to group boundary condition types
  - unified naming and steps how fluxes are calculated (nothing new was added)
  - further replacing _s, _n by 1, 2

boundaryLayerModel: spurious changes in (~4% for 1e-11 and 1% for 1e-9 in velocity)
---
 .../2cnistokes2p2cnilocaloperator.hh          | 191 ++++++++++-
 .../2cstokes2p2c/2cstokes2p2clocaloperator.hh | 314 ++++++++++++++++--
 ...stokes2p2cniboundarylayer-pm-reference.vtu |  80 ++---
 3 files changed, 492 insertions(+), 93 deletions(-)

diff --git a/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh b/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh
index 3f184ee5af..2e96f1089a 100644
--- a/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh
+++ b/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh
@@ -169,7 +169,7 @@ public:
     };
 
     // Stokes
-    enum { numComponents2 = Stokes2cniIndices::numComponents };
+    enum { numComponents1 = Stokes2cniIndices::numComponents };
     enum { // equation indices
         energyEqIdx1 = Stokes2cniIndices::energyEqIdx             //!< Index of the energy balance equation
     };
@@ -201,9 +201,167 @@ public:
     { }
 
 public:
+    //! \copydoc Dumux::TwoCStokesTwoPTwoCLocalOperator::evalCoupling()
+    template<typename LFSU1, typename LFSU2, typename RES1, typename RES2, typename CParams>
+    void evalCoupling(const LFSU1& lfsu1, const LFSU2& lfsu2,
+                      const int vertInElem1, const int vertInElem2,
+                      const SDElement1& sdElement1, const SDElement2& sdElement2,
+                      const BoundaryVariables1& boundaryVars1, const BoundaryVariables2& boundaryVars2,
+                      const CParams &cParams,
+                      RES1& couplingRes1, RES2& couplingRes2) const
+    {
+        // evaluate coupling of mass and momentum balances
+        ParentType::evalCoupling(lfsu1, lfsu2,
+                                 vertInElem1, vertInElem2,
+                                 sdElement1, sdElement2,
+                                 boundaryVars1, boundaryVars2,
+                                 cParams,
+                                 couplingRes1, couplingRes2);
+
+        const GlobalPosition& globalPos1 = cParams.fvGeometry1.subContVol[vertInElem1].global;
+        const GlobalPosition& globalPos2 = cParams.fvGeometry2.subContVol[vertInElem2].global;
+
+        // ENERGY Balance
+        // Neumann-like conditions
+        if (cParams.boundaryTypes1.isCouplingNeumann(energyEqIdx1))
+        {
+            if (this->globalProblem().sdProblem2().isCornerPoint(globalPos2))
+            {
+                // convective energy flux (enthalpy is mass based, mass flux also needed for useMoles)
+                Scalar convectiveFlux = 0.0;
+                for (int phaseIdx=0; phaseIdx<numPhases2; ++phaseIdx)
+                {
+                    convectiveFlux -= boundaryVars2.volumeFlux(phaseIdx)
+                                      * cParams.elemVolVarsCur2[vertInElem2].density(phaseIdx)
+                                      * cParams.elemVolVarsCur2[vertInElem2].enthalpy(phaseIdx);
+                }
+
+                // conductive energy flux
+                Scalar conductiveFlux = boundaryVars2.normalMatrixHeatFlux();
+
+                couplingRes1.accumulate(lfsu1.child(energyEqIdx1), vertInElem1,
+                                        -(convectiveFlux - conductiveFlux));
+            }
+            else
+            {
+                couplingRes1.accumulate(lfsu1.child(energyEqIdx1), vertInElem1,
+                                        this->globalProblem().localResidual2().residual(vertInElem2)[energyEqIdx2]);
+            }
+        }
+
+        // TODO: unify the behavior for cParams.boundaryTypes2.isCouplingNeumann()
+        //       with the different one in the isothermal LOP
+        if (cParams.boundaryTypes2.isCouplingNeumann(energyEqIdx2))
+        {
+            const GlobalPosition& bfNormal1 = boundaryVars1.face().normal;
+            // only enter here, if a boundary layer model is used for the computation of the diffusive fluxes
+            if (ParentType::blModel_)
+            {
+                // convective energy flux (enthalpy is mass based, mass flux also needed for useMoles)
+                Scalar convectiveFlux = boundaryVars1.normalVelocity()
+                                        * cParams.elemVolVarsCur1[vertInElem1].density()
+                                        * cParams.elemVolVarsCur1[vertInElem1].enthalpy();
+
+                // conductive energy flux
+                Scalar conductiveFlux = bfNormal1.two_norm()
+                                        * evalBoundaryLayerTemperatureGradient(cParams, vertInElem1)
+                                        * (boundaryVars1.thermalConductivity()
+                                           + boundaryVars1.thermalEddyConductivity());
+
+                // enthalpy transported by diffusive fluxes
+                Scalar sumDiffusiveFluxes = 0.0;
+                Scalar sumDiffusiveEnergyFlux = 0.0;
+                for (int compIdx=0; compIdx < numComponents1; compIdx++)
+                {
+                    if (compIdx != phaseCompIdx1)
+                    {
+                        Scalar diffusiveFlux = bfNormal1.two_norm()
+                                               * ParentType::template evalBoundaryLayerConcentrationGradient<CParams>(cParams, vertInElem1)
+                                               * (boundaryVars1.diffusionCoeff(compIdx)
+                                                  + boundaryVars1.eddyDiffusivity())
+                                               * boundaryVars1.molarDensity()
+                                               * ParentType::template evalMassTransferCoefficient<CParams>(cParams, vertInElem1, vertInElem2);
+                        sumDiffusiveFluxes += diffusiveFlux;
+                        sumDiffusiveEnergyFlux += diffusiveFlux
+                                                  * boundaryVars1.componentEnthalpy(compIdx)
+                                                  * FluidSystem::molarMass(compIdx); // Multiplied by molarMass [kg/mol] to convert from [mol/m^3 s] to [kg/m^3 s]
+                    }
+                }
+                sumDiffusiveEnergyFlux -= sumDiffusiveFluxes
+                                          * boundaryVars1.componentEnthalpy(phaseCompIdx1)
+                                          * FluidSystem::molarMass(phaseCompIdx1);
+
+                // TODO: use mass transfer coefficient here?
+                couplingRes2.accumulate(lfsu2.child(energyEqIdx2), vertInElem2,
+                                        -(convectiveFlux - sumDiffusiveEnergyFlux - conductiveFlux));
+            }
+            else if (this->globalProblem().sdProblem1().isCornerPoint(globalPos1))
+            {
+                // convective energy flux (enthalpy is mass based, mass flux also needed for useMoles)
+                Scalar convectiveFlux = boundaryVars1.normalVelocity()
+                                        * cParams.elemVolVarsCur1[vertInElem1].density()
+                                        * cParams.elemVolVarsCur1[vertInElem1].enthalpy();
+
+                // conductive energy flux
+                Scalar conductiveFlux = bfNormal1
+                                        * boundaryVars1.temperatureGrad()
+                                        * (boundaryVars1.thermalConductivity()
+                                           + boundaryVars1.thermalEddyConductivity());
+
+                // enthalpy transported by diffusive fluxes
+                Scalar sumDiffusiveFluxes = 0.0;
+                Scalar sumDiffusiveEnergyFlux = 0.0;
+                for (int compIdx=0; compIdx < numComponents1; compIdx++)
+                {
+                    if (compIdx != phaseCompIdx1)
+                    {
+                        Scalar diffusiveFlux = bfNormal1
+                                               * boundaryVars1.moleFractionGrad(compIdx)
+                                               * (boundaryVars1.diffusionCoeff(compIdx)
+                                                  + boundaryVars1.eddyDiffusivity())
+                                               * boundaryVars1.molarDensity();
+                        sumDiffusiveFluxes += diffusiveFlux;
+                        sumDiffusiveEnergyFlux += diffusiveFlux
+                                                  * boundaryVars1.componentEnthalpy(compIdx)
+                                                  * FluidSystem::molarMass(compIdx); // Multiplied by molarMass [kg/mol] to convert from [mol/m^3 s] to [kg/m^3 s]
+                    }
+                }
+                sumDiffusiveEnergyFlux -= sumDiffusiveFluxes
+                                          * boundaryVars1.componentEnthalpy(phaseCompIdx1)
+                                          * FluidSystem::molarMass(phaseCompIdx1);
+
+                couplingRes2.accumulate(lfsu2.child(energyEqIdx2), vertInElem2,
+                                        -(convectiveFlux - sumDiffusiveEnergyFlux - conductiveFlux));
+            }
+            else
+            {
+                couplingRes2.accumulate(lfsu2.child(energyEqIdx2), vertInElem2,
+                                        this->globalProblem().localResidual1().residual(vertInElem1)[energyEqIdx1]);
+            }
+        }
+
+
+        // ENERGY Balance
+        // Dirichlet-like conditions
+        if (cParams.boundaryTypes1.isCouplingDirichlet(energyEqIdx1))
+        {
+            // set residualStokes[energyIdx1] = T in stokesncnicouplinglocalresidual.hh
+            couplingRes1.accumulate(lfsu1.child(energyEqIdx1), vertInElem1,
+                                    -cParams.elemVolVarsCur2[vertInElem2].temperature());
+        }
+
+        if (cParams.boundaryTypes2.isCouplingDirichlet(energyEqIdx2))
+        {
+            // set residualDarcy[energyEqIdx2] = T in 2p2cnicouplinglocalresidual.hh
+            couplingRes2.accumulate(lfsu2.child(energyEqIdx2), vertInElem2,
+                                    -cParams.elemVolVarsCur1[vertInElem1].temperature());
+        }
+    }
+
     //! \copydoc Dumux::TwoCStokesTwoPTwoCLocalOperator::evalCoupling12()
     template<typename LFSU1, typename LFSU2, typename RES1, typename RES2, typename CParams>
-    void evalCoupling12(const LFSU1& lfsu_s, const LFSU2& lfsu_n,
+    DUNE_DEPRECATED_MSG("evalCoupling12() is deprecated. Use evalCoupling() instead.")
+    void evalCoupling12(const LFSU1& lfsu1, const LFSU2& lfsu2,
                         const int vertInElem1, const int vertInElem2,
                         const SDElement1& sdElement1, const SDElement2& sdElement2,
                         const BoundaryVariables1& boundaryVars1, const BoundaryVariables2& boundaryVars2,
@@ -217,7 +375,7 @@ public:
         GlobalProblem& globalProblem = this->globalProblem();
 
         // evaluate coupling of mass and momentum balances
-        ParentType::evalCoupling12(lfsu_s, lfsu_n,
+        ParentType::evalCoupling12(lfsu1, lfsu2,
                                    vertInElem1, vertInElem2,
                                    sdElement1, sdElement2,
                                    boundaryVars1, boundaryVars2,
@@ -235,7 +393,7 @@ public:
 
                 // enthalpy transported by diffusive fluxes
                 // multiply the diffusive flux with the mass transfer coefficient
-                static_assert(numComponents2 == 2,
+                static_assert(numComponents1 == 2,
                               "This coupling condition is only implemented for two components.");
                 Scalar diffusiveEnergyFlux = 0.0;
                 Scalar diffusiveFlux = bfNormal1.two_norm()
@@ -256,9 +414,7 @@ public:
                                               * (boundaryVars1.thermalConductivity()
                                                  + boundaryVars1.thermalEddyConductivity());
 
-                // TODO: unify this behavior with the one in the isothermal LOP
-                // TODO: use mass transfer coefficient here?
-                couplingRes2.accumulate(lfsu_n.child(energyEqIdx2), vertInElem2,
+                couplingRes2.accumulate(lfsu2.child(energyEqIdx2), vertInElem2,
                                         -(convectiveFlux - diffusiveEnergyFlux - conductiveFlux));
             }
             else if (globalProblem.sdProblem1().isCornerPoint(globalPos1))
@@ -272,7 +428,7 @@ public:
                     (boundaryVars1.thermalConductivity() + boundaryVars1.thermalEddyConductivity());
                 Scalar sumDiffusiveFluxes = 0.0;
                 Scalar sumDiffusiveEnergyFlux = 0.0;
-                for (int compIdx=0; compIdx < numComponents2; compIdx++)
+                for (int compIdx=0; compIdx < numComponents1; compIdx++)
                 {
                     if (compIdx != phaseCompIdx1)
                     {
@@ -288,27 +444,28 @@ public:
                 }
                 sumDiffusiveEnergyFlux -= sumDiffusiveFluxes * boundaryVars1.componentEnthalpy(phaseCompIdx1)
                                           * FluidSystem::molarMass(phaseCompIdx1);
-                couplingRes2.accumulate(lfsu_n.child(energyEqIdx2), vertInElem2,
+                couplingRes2.accumulate(lfsu2.child(energyEqIdx2), vertInElem2,
                                         -(convectiveFlux - sumDiffusiveEnergyFlux - conductiveFlux));
             }
             else
             {
-                // the energy flux from the stokes domain
-                couplingRes2.accumulate(lfsu_n.child(energyEqIdx2), vertInElem2,
+                // the energy flux from the Stokes domain
+                couplingRes2.accumulate(lfsu2.child(energyEqIdx2), vertInElem2,
                                         globalProblem.localResidual1().residual(vertInElem1)[energyEqIdx1]);
             }
         }
         if (cParams.boundaryTypes2.isCouplingDirichlet(energyEqIdx2))
         {
             // set residualDarcy[energyEqIdx2] = T in 2p2cnilocalresidual.hh
-            couplingRes2.accumulate(lfsu_n.child(energyEqIdx2), vertInElem2,
+            couplingRes2.accumulate(lfsu2.child(energyEqIdx2), vertInElem2,
                                     -cParams.elemVolVarsCur1[vertInElem1].temperature());
         }
     }
 
     //! \copydoc Dumux::TwoCStokesTwoPTwoCLocalOperator::evalCoupling21()
     template<typename LFSU1, typename LFSU2, typename RES1, typename RES2, typename CParams>
-    void evalCoupling21(const LFSU1& lfsu_s, const LFSU2& lfsu_n,
+    DUNE_DEPRECATED_MSG("evalCoupling21() is deprecated. Use evalCoupling() instead.")
+    void evalCoupling21(const LFSU1& lfsu1, const LFSU2& lfsu2,
                         const int vertInElem1, const int vertInElem2,
                         const SDElement1& sdElement1, const SDElement2& sdElement2,
                         const BoundaryVariables1& boundaryVars1, const BoundaryVariables2& boundaryVars2,
@@ -318,7 +475,7 @@ public:
         GlobalProblem& globalProblem = this->globalProblem();
 
         // evaluate coupling of mass and momentum balances
-        ParentType::evalCoupling21(lfsu_s, lfsu_n,
+        ParentType::evalCoupling21(lfsu1, lfsu2,
                                    vertInElem1, vertInElem2,
                                    sdElement1, sdElement2,
                                    boundaryVars1, boundaryVars2,
@@ -337,7 +494,7 @@ public:
         if (cParams.boundaryTypes1.isCouplingDirichlet(energyEqIdx1))
         {
             // set residualStokes[energyIdx1] = T in stokes2cnilocalresidual.hh
-            couplingRes1.accumulate(lfsu_s.child(energyEqIdx1), vertInElem1,
+            couplingRes1.accumulate(lfsu1.child(energyEqIdx1), vertInElem1,
                                     -cParams.elemVolVarsCur2[vertInElem2].temperature());
         }
         if (cParams.boundaryTypes1.isCouplingNeumann(energyEqIdx1))
@@ -351,12 +508,12 @@ public:
                     cParams.elemVolVarsCur2[vertInElem2].enthalpy(wPhaseIdx2);
                 const Scalar conductiveFlux = boundaryVars2.normalMatrixHeatFlux();
 
-                couplingRes1.accumulate(lfsu_s.child(energyEqIdx1), vertInElem1,
+                couplingRes1.accumulate(lfsu1.child(energyEqIdx1), vertInElem1,
                                         -(convectiveFlux - conductiveFlux));
             }
             else
             {
-                couplingRes1.accumulate(lfsu_s.child(energyEqIdx1), vertInElem1,
+                couplingRes1.accumulate(lfsu1.child(energyEqIdx1), vertInElem1,
                                         globalProblem.localResidual2().residual(vertInElem2)[energyEqIdx2]);
             }
         }
diff --git a/dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh b/dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh
index 01d28ad2ed..e84e67ca1e 100644
--- a/dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh
+++ b/dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh
@@ -173,6 +173,7 @@ public:
 
     // Stokes
     enum { numEq1 = GET_PROP_VALUE(Stokes2cTypeTag, NumEq) };
+    enum { numComponents1 = Stokes2cIndices::numComponents };
     enum { // equation indices
         momentumXIdx1 = Stokes2cIndices::momentumXIdx,         //!< Index of the x-component of the momentum balance
         momentumYIdx1 = Stokes2cIndices::momentumYIdx,         //!< Index of the y-component of the momentum balance
@@ -244,7 +245,6 @@ public:
      * \param lfsv2 local basis for the test space of the Darcy domain
      * \param couplingRes1 the coupling residual from the Stokes domain
      * \param couplingRes2 the coupling residual from the Darcy domain
-     *
      */
     template<typename IntersectionGeom, typename LFSU1, typename LFSU2,
              typename X, typename LFSV1, typename LFSV2,typename RES>
@@ -308,18 +308,12 @@ public:
                                                    cParams.elemVolVarsCur2,
                                                    /*onBoundary=*/true);
 
-            asImp_()->evalCoupling12(lfsu1, lfsu2, // local function spaces
-                                     vertInElem1, vertInElem2,
-                                     sdElement1, sdElement2,
-                                     boundaryVars1, boundaryVars2,
-                                     cParams,
-                                     couplingRes1, couplingRes2);
-            asImp_()->evalCoupling21(lfsu1, lfsu2, // local function spaces
-                                     vertInElem1, vertInElem2,
-                                     sdElement1, sdElement2,
-                                     boundaryVars1, boundaryVars2,
-                                     cParams,
-                                     couplingRes1, couplingRes2);
+            asImp_()->evalCoupling(lfsu1, lfsu2,
+                                   vertInElem1, vertInElem2,
+                                   sdElement1, sdElement2,
+                                   boundaryVars1, boundaryVars2,
+                                   cParams,
+                                   couplingRes1, couplingRes2);
         }
     }
 
@@ -334,7 +328,6 @@ public:
      * \param sdElement1 the element in the Stokes domain
      * \param sdElement2 the element in the Darcy domain
      * \param cParams a parameter container
-     *
      */
     template<typename LFSU1, typename LFSU2, typename X, typename CParams>
     void updateElemVolVars(const LFSU1& lfsu1, const LFSU2& lfsu2,
@@ -377,6 +370,256 @@ public:
 
     }
 
+    /*!
+     * \brief Evaluation of the coupling between the Stokes (1) and Darcy (2).
+     *
+     * Dirichlet-like and Neumann-like conditions for the respective domain are evaluated.
+     *
+     * \param lfsu1 local basis for the trial space of the Stokes domain
+     * \param lfsu2 local basis for the trial space of the Darcy domain
+     * \param vertInElem1 local vertex index in element1
+     * \param vertInElem2 local vertex index in element2
+     * \param sdElement1 the element in the Stokes domain
+     * \param sdElement2 the element in the Darcy domain
+     * \param boundaryVars1 the boundary variables at the interface of the Stokes domain
+     * \param boundaryVars2 the boundary variables at the interface of the Darcy domain
+     * \param cParams a parameter container
+     * \param couplingRes1 the coupling residual from the Stokes domain
+     * \param couplingRes2 the coupling residual from the Darcy domain
+     */
+    template<typename LFSU1, typename LFSU2, typename RES1, typename RES2, typename CParams>
+    void evalCoupling(const LFSU1& lfsu1, const LFSU2& lfsu2,
+                      const int vertInElem1, const int vertInElem2,
+                      const SDElement1& sdElement1, const SDElement2& sdElement2,
+                      const BoundaryVariables1& boundaryVars1, const BoundaryVariables2& boundaryVars2,
+                      const CParams &cParams,
+                      RES1& couplingRes1, RES2& couplingRes2) const
+    {
+        const GlobalPosition& globalPos1 = cParams.fvGeometry1.subContVol[vertInElem1].global;
+        const GlobalPosition& globalPos2 = cParams.fvGeometry2.subContVol[vertInElem2].global;
+
+        const GlobalPosition& bfNormal1 = boundaryVars1.face().normal;
+        const Scalar normalMassFlux1 = boundaryVars1.normalVelocity()
+                                       * cParams.elemVolVarsCur1[vertInElem1].density();
+
+        // MASS Balance
+        // Neumann-like conditions
+        if (cParams.boundaryTypes2.isCouplingNeumann(massBalanceIdx2))
+        {
+            static_assert(!GET_PROP_VALUE(TwoPTwoCTypeTag, UseMoles),
+                          "This coupling condition is only implemented for mass fraction formulation.");
+
+            if (globalProblem_.sdProblem1().isCornerPoint(globalPos1))
+            {
+                couplingRes2.accumulate(lfsu2.child(massBalanceIdx2), vertInElem2,
+                                        -normalMassFlux1);
+            }
+            else
+            {
+                couplingRes2.accumulate(lfsu2.child(massBalanceIdx2), vertInElem2,
+                                        globalProblem_.localResidual1().residual(vertInElem1)[massBalanceIdx1]);
+            }
+        }
+
+        // MASS Balance
+        // Dirichlet-like
+        if (cParams.boundaryTypes2.isCouplingDirichlet(massBalanceIdx2))
+        {
+            couplingRes2.accumulate(lfsu2.child(massBalanceIdx2), vertInElem2,
+                                    globalProblem_.localResidual1().residual(vertInElem1)[momentumYIdx1]
+                                    -cParams.elemVolVarsCur1[vertInElem1].pressure());
+        }
+
+
+        // MOMENTUM_X Balance
+        SpatialParams spatialParams = globalProblem_.sdProblem1().spatialParams();
+        Scalar beaversJosephCoeff = spatialParams.beaversJosephCoeffAtPos(globalPos1);
+        assert(beaversJosephCoeff > 0);
+        beaversJosephCoeff /= std::sqrt(spatialParams.intrinsicPermeability(sdElement1, cParams.fvGeometry1, vertInElem1));
+        // Neumann-like conditions
+
+        // TODO revise comment
+        // Implementation as Dirichlet-like condition
+        // tangential component: vx = sqrt K /alpha * (grad v n(unity))t
+        if (cParams.boundaryTypes1.isCouplingNeumann(momentumXIdx1))
+        {
+            // TODO: Is not implemented anymore because curPrimaryVars_ is protected
+            //       and it could be that this part of code has never been checked
+            std::cerr << "The boundary condition isCouplingNeumann(momentumXIdx1) on the Stokes side is not implemented anymore." << std::endl;
+            exit(1);
+
+            // GlobalPosition tangentialVelGrad(0);
+            // boundaryVars1.velocityGrad().umv(elementUnitNormal, tangentialVelGrad);
+            // tangentialVelGrad /= -beaversJosephCoeff; // was - before
+            // couplingRes1.accumulate(lfsu1.child(momentumXIdx1), vertInElem1,
+            // this->residual_[vertInElem1][momentumXIdx1] =
+            //        tangentialVelGrad[momentumXIdx1] - globalProblem_.localResidual1().curPriVars_(vertInElem1)[momentumXIdx1]);
+        }
+
+        // MOMENTUM_X Balance
+        // Dirichlet-like conditions
+
+        // TODO revise comment
+        // Implementation as Neumann-like condition: (v.n)n
+        if (cParams.boundaryTypes1.isCouplingDirichlet(momentumXIdx1))
+        {
+            const Scalar normalComp = boundaryVars1.velocity()*bfNormal1;
+            GlobalPosition normalV = bfNormal1;
+            normalV *= normalComp; // v*n*n
+
+            // v_tau = v - (v.n)n
+            const GlobalPosition tangentialV = boundaryVars1.velocity() - normalV;
+
+            for (int dimIdx=0; dimIdx < dim; ++dimIdx)
+            {
+                couplingRes1.accumulate(lfsu1.child(momentumXIdx1), vertInElem1,
+                                        beaversJosephCoeff
+                                        * boundaryVars1.face().area
+                                        * tangentialV[dimIdx]
+                                        * (boundaryVars1.dynamicViscosity()
+                                          + boundaryVars1.dynamicEddyViscosity()));
+            }
+        }
+
+
+        // MOMENTUM_Y Balance
+        // Neumann-like conditions
+
+        // TODO revise comment
+        // v.n as Dirichlet-like condition for the Stokes domain
+        // set residualStokes[momentumYIdx1] = v_y in stokesnccouplinglocalresidual.hh
+        if (cParams.boundaryTypes1.isCouplingNeumann(momentumYIdx1))
+        {
+            if (globalProblem_.sdProblem2().isCornerPoint(globalPos2))
+            {
+                Scalar sumNormalPhaseFluxes = 0.0;
+                for (int phaseIdx=0; phaseIdx<numPhases2; ++phaseIdx)
+                {
+                    sumNormalPhaseFluxes -= boundaryVars2.volumeFlux(phaseIdx)
+                                            * cParams.elemVolVarsCur2[vertInElem2].density(phaseIdx);
+                }
+                couplingRes1.accumulate(lfsu1.child(momentumYIdx1), vertInElem1,
+                                        -sumNormalPhaseFluxes
+                                        / cParams.elemVolVarsCur1[vertInElem1].density());
+            }
+            else
+            {
+                // TODO revise comment
+                // v.n as DIRICHLET condition for the Stokes domain (negative sign!)
+                couplingRes1.accumulate(lfsu1.child(momentumYIdx1), vertInElem1,
+                                        globalProblem_.localResidual2().residual(vertInElem2)[massBalanceIdx2]
+                                        / cParams.elemVolVarsCur1[vertInElem1].density());
+            }
+        }
+
+        // MOMENTUM_Y Balance
+        // Dirichlet-like conditions
+
+        // TODO revise comments
+        //p*n as NEUMANN condition for free flow (set, if B&J defined as NEUMANN condition)
+        // p*A*n as NEUMANN condition for free flow (set, if B&J defined as NEUMANN condition)
+        if (cParams.boundaryTypes1.isCouplingDirichlet(momentumYIdx1))
+        {
+            // pressure correction is done in stokeslocalresidual.hh
+            couplingRes1.accumulate(lfsu1.child(momentumYIdx1), vertInElem1,
+                                    cParams.elemVolVarsCur2[vertInElem2].pressure(nPhaseIdx2) *
+                                    boundaryVars2.face().area);
+        }
+
+
+        // COMPONENT Balance
+        // Neumann-like conditions
+        if (cParams.boundaryTypes1.isCouplingNeumann(transportEqIdx1))
+        {
+            std::cerr << "The boundary condition isCouplingNeumann(transportEqIdx1) is not implemented";
+            std::cerr << "for the Stokes side for multicomponent systems." << std::endl;
+            exit(1);
+        }
+        if (cParams.boundaryTypes2.isCouplingNeumann(contiWEqIdx2))
+        {
+            // only enter here, if a boundary layer model is used for the computation of the diffusive fluxes
+            if (blModel_)
+            {
+                Scalar diffusiveFlux = bfNormal1.two_norm()
+                                       * evalBoundaryLayerConcentrationGradient<CParams>(cParams, vertInElem1)
+                                       * (boundaryVars1.diffusionCoeff(transportCompIdx1)
+                                         + boundaryVars1.eddyDiffusivity())
+                                       * boundaryVars1.molarDensity()
+                                       * FluidSystem::molarMass(transportCompIdx1);
+
+                Scalar advectiveFlux = normalMassFlux1
+                                       * cParams.elemVolVarsCur1[vertInElem1].massFraction(transportCompIdx1);
+
+                const Scalar massTransferCoeff = evalMassTransferCoefficient<CParams>(cParams, vertInElem1, vertInElem2);
+                // TODO: unify this behavior with the one in the non-isothermal LOP
+                if (globalProblem_.sdProblem1().isCornerPoint(globalPos1) && massTransferModel_)
+                {
+                    Scalar diffusiveFluxAtCorner = bfNormal1
+                                                   * boundaryVars1.moleFractionGrad(transportCompIdx1)
+                                                   * (boundaryVars1.diffusionCoeff(transportCompIdx1)
+                                                      + boundaryVars1.eddyDiffusivity())
+                                                   * boundaryVars1.molarDensity()
+                                                   * FluidSystem::molarMass(transportCompIdx1);
+
+                    couplingRes2.accumulate(lfsu2.child(contiWEqIdx2), vertInElem2,
+                                            -massTransferCoeff*(advectiveFlux - diffusiveFlux) -
+                                            (1.-massTransferCoeff)*(advectiveFlux - diffusiveFluxAtCorner));
+                }
+                else
+                {
+                    couplingRes2.accumulate(lfsu2.child(contiWEqIdx2), vertInElem2,
+                                            -massTransferCoeff*(advectiveFlux - diffusiveFlux) +
+                                            (1.-massTransferCoeff)*globalProblem_.localResidual1().residual(vertInElem1)[transportEqIdx1]);
+                }
+            }
+            else if (globalProblem_.sdProblem1().isCornerPoint(globalPos1))
+            {
+                static_assert(!GET_PROP_VALUE(TwoPTwoCTypeTag, UseMoles),
+                              "This coupling condition is only implemented for mass fraction formulation.");
+
+                Scalar advectiveFlux = normalMassFlux1
+                                       * cParams.elemVolVarsCur1[vertInElem1].massFraction(transportCompIdx1);
+                Scalar diffusiveFlux = bfNormal1
+                                       * boundaryVars1.moleFractionGrad(transportCompIdx1)
+                                       * (boundaryVars1.diffusionCoeff(transportCompIdx1)
+                                          + boundaryVars1.eddyDiffusivity())
+                                       * boundaryVars1.molarDensity()
+                                       * FluidSystem::molarMass(transportCompIdx1);
+
+                couplingRes2.accumulate(lfsu2.child(contiWEqIdx2), vertInElem2,
+                                        -(advectiveFlux - diffusiveFlux));
+            }
+            else
+            {
+                static_assert(GET_PROP_VALUE(Stokes2cTypeTag, UseMoles) == GET_PROP_VALUE(TwoPTwoCTypeTag, UseMoles),
+                              "This coupling condition is not implemented for different formulations (mass/mole) in the subdomains.");
+
+                // the component mass flux from the stokes domain
+                couplingRes2.accumulate(lfsu2.child(contiWEqIdx2), vertInElem2,
+                                        globalProblem_.localResidual1().residual(vertInElem1)[transportEqIdx1]);
+            }
+        }
+
+        // COMPONENT Balance
+        // Dirichlet-like conditions (coupling residual is added to "real" residual)
+        // TODO (stimmt der Kommentar): here each node is passed twice, hence only half of the dirichlet condition has to be set
+        if (cParams.boundaryTypes1.isCouplingDirichlet(transportEqIdx1))
+        {
+            // set residualStokes[transportEqIdx1] = x in stokesnccouplinglocalresidual.hh
+            static_assert(!GET_PROP_VALUE(Stokes2cTypeTag, UseMoles),
+                          "This coupling condition is only implemented for mass fraction formulation.");
+
+            couplingRes1.accumulate(lfsu1.child(transportEqIdx1), vertInElem1,
+                                    -cParams.elemVolVarsCur2[vertInElem2].massFraction(nPhaseIdx2, wCompIdx2));
+        }
+        if (cParams.boundaryTypes2.isCouplingDirichlet(contiWEqIdx2))
+        {
+            std::cerr << "The boundary condition isCouplingDirichlet(contiWEqIdx2) is not implemented";
+            std::cerr << "for the Darcy side for multicomponent systems." << std::endl;
+            exit(1);
+        }
+    }
+
     /*!
      * \brief Evaluation of the coupling from Stokes (1) to Darcy (2).
      *
@@ -393,6 +636,7 @@ public:
      * \param couplingRes2 the coupling residual from the Darcy domain
      */
     template<typename LFSU1, typename LFSU2, typename RES1, typename RES2, typename CParams>
+    DUNE_DEPRECATED_MSG("evalCoupling21() is deprecated. Use evalCoupling() instead.")
     void evalCoupling12(const LFSU1& lfsu1, const LFSU2& lfsu2,
                         const int vertInElem1, const int vertInElem2,
                         const SDElement1& sdElement1, const SDElement2& sdElement2,
@@ -446,7 +690,6 @@ public:
                 Scalar advectiveFlux = normalMassFlux * cParams.elemVolVarsCur1[vertInElem1].massFraction(transportCompIdx1);
 
                 const Scalar massTransferCoeff = evalMassTransferCoefficient<CParams>(cParams, vertInElem1, vertInElem2);
-                // TODO: unify this behavior with the one in the non-isothermal LOP
                 if (globalProblem_.sdProblem1().isCornerPoint(globalPos1) && massTransferModel_)
                 {
                     const Scalar diffusiveFluxAtCorner =
@@ -488,16 +731,18 @@ public:
             else
             {
                 static_assert(GET_PROP_VALUE(Stokes2cTypeTag, UseMoles) == GET_PROP_VALUE(TwoPTwoCTypeTag, UseMoles),
-                              "This coupling condition is only implemented or different formulations (mass/mole) in the subdomains.");
+                              "This coupling condition is not implemented dor different formulations (mass/mole) in the subdomains.");
 
                 // the component mass flux from the stokes domain
                 couplingRes2.accumulate(lfsu2.child(contiWEqIdx2), vertInElem2,
                                         globalProblem_.localResidual1().residual(vertInElem1)[transportEqIdx1]);
             }
         }
-        // TODO make this a static assert
         if (cParams.boundaryTypes2.isCouplingDirichlet(contiWEqIdx2))
-            std::cerr << "Upwind PM -> FF does not work for the transport equation for a 2-phase system!" << std::endl;
+        {
+            std::cerr << "The boundary condition isCouplingDirichlet(contiWEqIdx2) is not implemented for the Darcy side." << std::endl;
+            exit(1);
+        }
     }
 
     /*!
@@ -516,6 +761,7 @@ public:
      * \param couplingRes2 the coupling residual from the Darcy domain
      */
     template<typename LFSU1, typename LFSU2, typename RES1, typename RES2, typename CParams>
+    DUNE_DEPRECATED_MSG("evalCoupling21() is deprecated. Use evalCoupling() instead.")
     void evalCoupling21(const LFSU1& lfsu1, const LFSU2& lfsu2,
                         const int vertInElem1, const int vertInElem2,
                         const SDElement1& sdElement1, const SDElement2& sdElement2,
@@ -531,7 +777,6 @@ public:
             normalFlux2[phaseIdx] = -boundaryVars2.volumeFlux(phaseIdx)*
                 cParams.elemVolVarsCur2[vertInElem2].density(phaseIdx);
 
-        // TODO revise comment
         //p*n as NEUMANN condition for free flow (set, if B&J defined as NEUMANN condition)
         if (cParams.boundaryTypes1.isCouplingDirichlet(momentumYIdx1))
         {
@@ -543,8 +788,7 @@ public:
         }
         if (cParams.boundaryTypes1.isCouplingNeumann(momentumYIdx1))
         {
-            // TODO revise comment and move upwards
-            // v.n as Dirichlet condition for the Stokes domain
+            // v.n as Dirichlet-like condition for the Stokes domain
             // set residualStokes[momentumYIdx1] = vy in stokeslocalresidual.hh
             if (globalProblem_.sdProblem2().isCornerPoint(globalPos2))
             {
@@ -554,7 +798,6 @@ public:
             }
             else
             {
-                // TODO revise comment and move upwards
                 // v.n as DIRICHLET condition for the Stokes domain (negative sign!)
                 couplingRes1.accumulate(lfsu1.child(momentumYIdx1), vertInElem1,
                                         globalProblem_.localResidual2().residual(vertInElem2)[massBalanceIdx2]
@@ -573,8 +816,7 @@ public:
         beaversJosephCoeff /= std::sqrt(Kxx);
         const GlobalPosition& elementUnitNormal = boundaryVars1.face().normal;
 
-        // TODO revise comment
-        // Implementation as Neumann condition: (v.n)n
+        // Implementation as Neumann-like condition: (v.n)n
         if (cParams.boundaryTypes1.isCouplingDirichlet(momentumXIdx1))
         {
             const Scalar normalComp = boundaryVars1.velocity()*elementUnitNormal;
@@ -595,16 +837,15 @@ public:
                                           + boundaryVars1.dynamicEddyViscosity()));
             }
         }
-        // TODO revise comment
-        // Implementation as Dirichlet condition
+        // Implementation as Dirichlet-like condition
         // tangential component: vx = sqrt K /alpha * (grad v n(unity))t
         if (cParams.boundaryTypes1.isCouplingNeumann(momentumXIdx1))
         {
-            GlobalPosition tangentialVelGrad(0);
-            boundaryVars1.velocityGrad().umv(elementUnitNormal, tangentialVelGrad);
-            tangentialVelGrad /= -beaversJosephCoeff; // was - before
-            // TODO: 3 lines below not implemtented because curPrimaryVars_ is protected
-            // it could be that this part of code has never been checked
+            std::cerr << "The boundary conditionisCouplingNeumann(momentumXIdx1) on the Stokes side is not implemented anymore." << std::endl;
+            exit(1);
+            // GlobalPosition tangentialVelGrad(0);
+            // boundaryVars1.velocityGrad().umv(elementUnitNormal, tangentialVelGrad);
+            // tangentialVelGrad /= -beaversJosephCoeff; // was - before
             // couplingRes1.accumulate(lfsu1.child(momentumXIdx1), vertInElem1,
             // this->residual_[vertInElem1][momentumXIdx1] =
             //        tangentialVelGrad[momentumXIdx1] - globalProblem_.localResidual1().curPriVars_(vertInElem1)[momentumXIdx1]);
@@ -623,9 +864,11 @@ public:
             couplingRes1.accumulate(lfsu1.child(transportEqIdx1), vertInElem1,
                                     -cParams.elemVolVarsCur2[vertInElem2].massFraction(nPhaseIdx2, wCompIdx2));
         }
-        // TODO make this a static assert
         if (cParams.boundaryTypes1.isCouplingNeumann(transportEqIdx1))
-            std::cerr << "Upwind PM -> FF does not work for the transport equation for a 2-phase system!" << std::endl;
+        {
+            std::cerr << "The boundary condition isCouplingNeumann(transportEqIdx1) is not implemented for the Stokes side!" << std::endl;
+            exit(1);
+        }
     }
 
     /*!
@@ -642,9 +885,6 @@ public:
     template<typename CParams>
     BoundaryLayerModel<TypeTag> evalBoundaryLayerModel(CParams cParams, const int scvIdx1) const
     {
-        static_assert(!GET_PROP_VALUE(Stokes2cTypeTag, UseMoles),
-                      "Boundary layer and mass transfer models are only implemented for mass fraction formulation.");
-
         const Scalar velocity = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefVelocity);
         // current position + additional virtual runup distance
         const Scalar distance = cParams.fvGeometry1.subContVol[scvIdx1].global[0]
@@ -675,6 +915,8 @@ public:
     template<typename CParams>
     Scalar evalBoundaryLayerConcentrationGradient(CParams cParams, const int scvIdx1) const
     {
+        static_assert(numComponents1 == 2,
+                      "This coupling condition is only implemented for two components.");
         Scalar massFractionOut = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, FreeFlow, RefMassfrac);
         Scalar M1 = FluidSystem::molarMass(transportCompIdx1);
         Scalar M2 = FluidSystem::molarMass(phaseCompIdx1);
diff --git a/test/references/2cnistokes2p2cniboundarylayer-pm-reference.vtu b/test/references/2cnistokes2p2cniboundarylayer-pm-reference.vtu
index 62b6fee57c..d10fe9587d 100644
--- a/test/references/2cnistokes2p2cniboundarylayer-pm-reference.vtu
+++ b/test/references/2cnistokes2p2cniboundarylayer-pm-reference.vtu
@@ -24,13 +24,13 @@
           99017 99017 99017 99017 98715.9 98715.9 98715.9 98715.9 98715.9 98715.9 98715.9
         </DataArray>
         <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
-          17.4518 17.4517 63.3252 63.324 17.4518 63.323 17.4512 63.3198 17.4513 63.3176 17.4511 63.3157
-          17.4512 63.3149 532.36 532.359 532.359 532.358 532.358 532.358 532.357 983.379 983.379 983.378
+          17.4521 17.4516 63.3319 63.3243 17.452 63.3231 17.4518 63.3203 17.4514 63.3171 17.4514 63.3147
+          17.4512 63.3145 532.36 532.359 532.359 532.358 532.358 532.358 532.357 983.379 983.379 983.378
           983.378 983.378 983.378 983.377 1284.13 1284.13 1284.13 1284.13 1284.13 1284.13 1284.13
         </DataArray>
         <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
           997.056 997.056 997.076 997.076 997.056 997.076 997.056 997.076 997.056 997.076 997.056 997.076
-          997.056 997.076 997.18 997.18 997.18 997.18 997.18 997.18 997.18 997.401 997.402 997.403
+          997.056 997.076 997.179 997.18 997.18 997.18 997.18 997.18 997.18 997.401 997.402 997.403
           997.403 997.404 997.405 997.405 997.634 997.636 997.638 997.64 997.641 997.642 997.643
         </DataArray>
         <DataArray type="Float32" Name="rhoN" NumberOfComponents="1" format="ascii">
@@ -39,14 +39,14 @@
           1.16082 1.16084 1.16085 1.16085 1.16537 1.1654 1.16544 1.16547 1.1655 1.16552 1.16554
         </DataArray>
         <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
-          1123.57 1123.57 1121.58 1121.58 1123.57 1121.57 1123.57 1121.57 1123.57 1121.57 1123.57 1121.57
-          1123.57 1121.57 1110.06 1110.05 1110.04 1110.02 1110 1109.98 1109.98 1060.7 1060.65 1060.57
-          1060.49 1060.42 1060.37 1060.34 614.273 614.163 614.052 613.953 613.866 613.792 613.738
+          1123.57 1123.57 1121.57 1121.58 1123.57 1121.57 1123.57 1121.57 1123.57 1121.57 1123.57 1121.57
+          1123.57 1121.57 1110.08 1110.05 1110.04 1110.02 1110 1109.98 1109.98 1060.7 1060.65 1060.57
+          1060.49 1060.42 1060.37 1060.34 614.273 614.164 614.052 613.953 613.866 613.792 613.738
         </DataArray>
         <DataArray type="Float32" Name="mobN" NumberOfComponents="1" format="ascii">
-          0.000310649 0.000310643 0.00469836 0.00469816 0.000310648 0.00469801 0.000310627 0.0046975 0.000310629 0.00469715 0.000310623 0.00469685
-          0.000310626 0.00469672 0.604801 0.6048 0.604799 0.604799 0.604798 0.604797 0.604795 20.9933 20.9931 20.9931
-          20.9932 20.9932 20.9931 20.993 1292.53 1292.55 1292.56 1292.58 1292.59 1292.61 1292.61
+          0.000310659 0.000310639 0.00469941 0.0046982 0.000310657 0.00469802 0.00031065 0.00469757 0.000310631 0.00469708 0.000310632 0.0046967
+          0.000310627 0.00469666 0.6048 0.6048 0.604799 0.604798 0.604798 0.604797 0.604795 20.9933 20.9931 20.9931
+          20.9932 20.9932 20.9931 20.993 1292.53 1292.55 1292.56 1292.58 1292.59 1292.6 1292.61
         </DataArray>
         <DataArray type="Float32" Name="X_liquid^H2O" NumberOfComponents="1" format="ascii">
           0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978
@@ -54,18 +54,18 @@
           0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978
         </DataArray>
         <DataArray type="Float32" Name="X_liquid^Air" NumberOfComponents="1" format="ascii">
-          2.18092e-05 2.18092e-05 2.16226e-05 2.16226e-05 2.18092e-05 2.16227e-05 2.18092e-05 2.16227e-05 2.18092e-05 2.16227e-05 2.18092e-05 2.16228e-05
-          2.18092e-05 2.16228e-05 2.17329e-05 2.1733e-05 2.17333e-05 2.17336e-05 2.17339e-05 2.17341e-05 2.17341e-05 2.20832e-05 2.2084e-05 2.20853e-05
+          2.18092e-05 2.18092e-05 2.16228e-05 2.16227e-05 2.18092e-05 2.16227e-05 2.18092e-05 2.16227e-05 2.18092e-05 2.16227e-05 2.18092e-05 2.16228e-05
+          2.18092e-05 2.16228e-05 2.17327e-05 2.1733e-05 2.17333e-05 2.17336e-05 2.17339e-05 2.17341e-05 2.17341e-05 2.20832e-05 2.2084e-05 2.20853e-05
           2.20866e-05 2.20877e-05 2.20885e-05 2.2089e-05 2.24795e-05 2.24826e-05 2.24858e-05 2.24886e-05 2.24911e-05 2.24932e-05 2.24948e-05
         </DataArray>
         <DataArray type="Float32" Name="X_gas^H2O" NumberOfComponents="1" format="ascii">
-          0.0197216 0.0197216 0.019822 0.019822 0.0197216 0.0198219 0.0197216 0.0198218 0.0197216 0.0198217 0.0197216 0.0198216
-          0.0197216 0.0198216 0.0193764 0.0193759 0.0193752 0.0193743 0.0193734 0.0193728 0.0193726 0.0183488 0.0183467 0.018343
+          0.0197216 0.0197216 0.0198215 0.0198219 0.0197216 0.0198219 0.0197216 0.0198218 0.0197216 0.0198217 0.0197216 0.0198216
+          0.0197216 0.0198216 0.0193771 0.0193761 0.0193752 0.0193743 0.0193734 0.0193728 0.0193726 0.0183488 0.0183467 0.018343
           0.0183395 0.0183364 0.0183339 0.0183327 0.0172823 0.0172743 0.0172661 0.0172588 0.0172525 0.0172471 0.017243
         </DataArray>
         <DataArray type="Float32" Name="X_gas^Air" NumberOfComponents="1" format="ascii">
-          0.980278 0.980278 0.980178 0.980178 0.980278 0.980178 0.980278 0.980178 0.980278 0.980178 0.980278 0.980178
-          0.980278 0.980178 0.980624 0.980624 0.980625 0.980626 0.980627 0.980627 0.980627 0.981651 0.981653 0.981657
+          0.980278 0.980278 0.980179 0.980178 0.980278 0.980178 0.980278 0.980178 0.980278 0.980178 0.980278 0.980178
+          0.980278 0.980178 0.980623 0.980624 0.980625 0.980626 0.980627 0.980627 0.980627 0.981651 0.981653 0.981657
           0.98166 0.981664 0.981666 0.981667 0.982718 0.982726 0.982734 0.982741 0.982747 0.982753 0.982757
         </DataArray>
         <DataArray type="Float32" Name="x_liquid^H2O" NumberOfComponents="1" format="ascii">
@@ -74,18 +74,18 @@
           0.999986 0.999986 0.999986 0.999986 0.999986 0.999986 0.999986 0.999986 0.999986 0.999986 0.999986
         </DataArray>
         <DataArray type="Float32" Name="x_liquid^Air" NumberOfComponents="1" format="ascii">
-          1.3567e-05 1.3567e-05 1.34509e-05 1.34509e-05 1.3567e-05 1.3451e-05 1.3567e-05 1.3451e-05 1.3567e-05 1.3451e-05 1.3567e-05 1.3451e-05
-          1.3567e-05 1.3451e-05 1.35195e-05 1.35196e-05 1.35198e-05 1.352e-05 1.35201e-05 1.35202e-05 1.35203e-05 1.37375e-05 1.37379e-05 1.37387e-05
+          1.3567e-05 1.3567e-05 1.3451e-05 1.34509e-05 1.3567e-05 1.3451e-05 1.3567e-05 1.3451e-05 1.3567e-05 1.3451e-05 1.3567e-05 1.3451e-05
+          1.3567e-05 1.3451e-05 1.35194e-05 1.35196e-05 1.35198e-05 1.352e-05 1.35201e-05 1.35202e-05 1.35203e-05 1.37375e-05 1.37379e-05 1.37387e-05
           1.37395e-05 1.37402e-05 1.37408e-05 1.3741e-05 1.3984e-05 1.39859e-05 1.39879e-05 1.39896e-05 1.39912e-05 1.39925e-05 1.39935e-05
         </DataArray>
         <DataArray type="Float32" Name="x_gas^H2O" NumberOfComponents="1" format="ascii">
-          0.0313277 0.0313277 0.0314853 0.0314854 0.0313277 0.0314851 0.0313277 0.031485 0.0313277 0.0314849 0.0313277 0.0314848
-          0.0313277 0.0314847 0.0307857 0.0307851 0.0307839 0.0307824 0.0307811 0.0307802 0.0307798 0.0291711 0.0291679 0.0291621
+          0.0313277 0.0313277 0.0314846 0.0314852 0.0313277 0.0314851 0.0313277 0.031485 0.0313277 0.0314849 0.0313277 0.0314848
+          0.0313277 0.0314847 0.0307868 0.0307853 0.0307839 0.0307824 0.0307811 0.0307802 0.0307798 0.0291712 0.0291679 0.0291621
           0.0291565 0.0291516 0.0291477 0.0291458 0.0274932 0.0274805 0.0274677 0.0274563 0.0274463 0.0274377 0.0274314
         </DataArray>
         <DataArray type="Float32" Name="x_gas^Air" NumberOfComponents="1" format="ascii">
           0.968672 0.968672 0.968515 0.968515 0.968672 0.968515 0.968672 0.968515 0.968672 0.968515 0.968672 0.968515
-          0.968672 0.968515 0.969214 0.969215 0.969216 0.969218 0.969219 0.96922 0.96922 0.970829 0.970832 0.970838
+          0.968672 0.968515 0.969213 0.969215 0.969216 0.969218 0.969219 0.96922 0.96922 0.970829 0.970832 0.970838
           0.970843 0.970848 0.970852 0.970854 0.972507 0.972519 0.972532 0.972544 0.972554 0.972562 0.972569
         </DataArray>
         <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
@@ -94,8 +94,8 @@
           0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41
         </DataArray>
         <DataArray type="Float32" Name="temperature" NumberOfComponents="1" format="ascii">
-          298.15 298.15 298.073 298.073 298.15 298.073 298.15 298.072 298.15 298.072 298.15 298.072
-          298.15 298.072 297.662 297.661 297.661 297.66 297.659 297.659 297.658 296.764 296.762 296.759
+          298.15 298.15 298.072 298.073 298.15 298.073 298.15 298.072 298.15 298.072 298.15 298.072
+          298.15 298.072 297.662 297.661 297.661 297.66 297.659 297.659 297.658 296.764 296.763 296.759
           296.756 296.753 296.751 296.75 295.785 295.777 295.769 295.763 295.757 295.751 295.748
         </DataArray>
         <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
@@ -104,26 +104,26 @@
           3 3 3 3 3 3 3 3 3 3 3
         </DataArray>
         <DataArray type="Float32" Name="velocityW" NumberOfComponents="3" format="ascii">
-          -9.58175e-11 -4.2546e-10 0 -1.81911e-10 -4.49609e-10 0 -1.63525e-10 -1.92698e-09 0 -2.8296e-10 -2.03582e-09 0
-          -3.08578e-10 -5.63133e-10 0 -4.53465e-10 -2.28574e-09 0 -3.4457e-10 -6.7563e-10 0 -4.89316e-10 -2.53587e-09 0
-          -2.92234e-10 -7.83235e-10 0 -4.08145e-10 -2.75011e-09 0 -1.66588e-10 -8.58339e-10 0 -2.34515e-10 -2.91054e-09 0
-          -8.86966e-11 -8.85282e-10 0 -1.26835e-10 -2.97501e-09 0 -2.69547e-10 -2.39849e-09 0 -3.74168e-10 -2.80675e-09 0
-          -4.79685e-10 -3.34372e-09 0 -4.48735e-10 -3.78399e-09 0 -3.68919e-10 -4.13262e-09 0 -2.33728e-10 -4.43756e-09 0
-          -1.46506e-10 -4.61757e-09 0 2.98484e-11 1.2079e-08 0 2.02847e-10 1.12857e-08 0 4.53404e-10 1.06743e-08 0
-          5.20395e-10 1.0256e-08 0 4.43935e-10 9.92635e-09 0 2.51573e-10 9.61285e-09 0 1.25101e-10 9.332e-09 0
-          7.56158e-10 2.55265e-08 0 9.30427e-10 2.45628e-08 0 1.13512e-09 2.40277e-08 0 1.10738e-09 2.36838e-08 0
-          9.4487e-10 2.3401e-08 0 6.7723e-10 2.31381e-08 0 5.13928e-10 2.28344e-08 0
+          -6.28852e-11 -5.07559e-10 0 -1.54489e-10 -4.35918e-10 0 -1.932e-10 -1.94939e-09 0 -3.07104e-10 -2.03292e-09 0
+          -2.94124e-10 -5.44565e-10 0 -4.647e-10 -2.28172e-09 0 -3.39561e-10 -6.70764e-10 0 -4.91237e-10 -2.53483e-09 0
+          -2.90013e-10 -7.81291e-10 0 -4.07792e-10 -2.74885e-09 0 -1.65645e-10 -8.57174e-10 0 -2.33954e-10 -2.90881e-09 0
+          -8.82313e-11 -8.8423e-10 0 -1.26406e-10 -2.97319e-09 0 -2.87601e-10 -2.33343e-09 0 -3.8672e-10 -2.8227e-09 0
+          -4.81771e-10 -3.36184e-09 0 -4.45386e-10 -3.78779e-09 0 -3.65677e-10 -4.13146e-09 0 -2.31917e-10 -4.4343e-09 0
+          -1.4555e-10 -4.61386e-09 0 7.90956e-11 1.21311e-08 0 2.45312e-10 1.127e-08 0 4.79148e-10 1.06577e-08 0
+          5.32826e-10 1.02534e-08 0 4.5088e-10 9.92886e-09 0 2.54743e-10 9.61718e-09 0 1.26609e-10 9.33675e-09 0
+          7.846e-10 2.55378e-08 0 9.5871e-10 2.45554e-08 0 1.15555e-09 2.40202e-08 0 1.11695e-09 2.36836e-08 0
+          9.49721e-10 2.34042e-08 0 6.79347e-10 2.31425e-08 0 5.14868e-10 2.28391e-08 0
         </DataArray>
         <DataArray type="Float32" Name="velocityN" NumberOfComponents="3" format="ascii">
-          8.2456e-16 7.67906e-10 0 5.07216e-16 7.67892e-10 0 3.02122e-14 2.24876e-09 0 2.55256e-14 2.24867e-09 0
-          1.27307e-15 7.67906e-10 0 4.75043e-14 2.24861e-09 0 1.63519e-15 7.67854e-10 0 6.1894e-14 2.24836e-09 0
-          1.00968e-15 7.67862e-10 0 4.62091e-14 2.24821e-09 0 6.87663e-16 7.67848e-10 0 3.1429e-14 2.24807e-09 0
-          2.70063e-16 7.67856e-10 0 2.00584e-14 2.24802e-09 0 2.40055e-12 1.74641e-09 0 2.12686e-12 1.74973e-09 0
-          2.24766e-12 1.74615e-09 0 2.46124e-12 1.74471e-09 0 2.21766e-12 1.74218e-09 0 2.25344e-12 1.74032e-09 0
-          2.35189e-12 1.74301e-09 0 8.32294e-11 -1.49414e-08 0 6.54632e-11 -1.48446e-08 0 4.61632e-11 -1.48703e-08 0
-          4.49127e-11 -1.48792e-08 0 4.58852e-11 -1.48882e-08 0 6.65314e-11 -1.49122e-08 0 8.64886e-11 -1.48107e-08 0
-          5.4352e-09 -2.9646e-08 0 4.13295e-09 -2.94592e-08 0 2.80761e-09 -2.95036e-08 0 2.77917e-09 -2.95189e-08 0
-          2.78594e-09 -2.95322e-08 0 4.19428e-09 -2.95767e-08 0 5.5905e-09 -2.93793e-08 0
+          4.55767e-15 7.67927e-10 0 2.17238e-15 7.67881e-10 0 1.7326e-13 2.24923e-09 0 9.91557e-14 2.24869e-09 0
+          6.94714e-16 7.67928e-10 0 4.43153e-14 2.24862e-09 0 1.8945e-15 7.67911e-10 0 6.7715e-14 2.24842e-09 0
+          1.6398e-15 7.67868e-10 0 6.31986e-14 2.24818e-09 0 6.88169e-16 7.67872e-10 0 3.06159e-14 2.24801e-09 0
+          2.834e-16 7.67859e-10 0 6.68517e-15 2.24799e-09 0 4.76421e-12 1.74015e-09 0 3.45597e-12 1.75153e-09 0
+          2.31464e-12 1.74674e-09 0 2.57776e-12 1.74437e-09 0 2.53776e-12 1.74235e-09 0 2.25593e-12 1.74072e-09 0
+          2.11033e-12 1.74267e-09 0 7.91871e-11 -1.49529e-08 0 6.39269e-11 -1.48456e-08 0 4.67667e-11 -1.48691e-08 0
+          4.4907e-11 -1.4879e-08 0 4.56953e-11 -1.48867e-08 0 6.66551e-11 -1.49101e-08 0 8.6867e-11 -1.48095e-08 0
+          5.43433e-09 -2.96555e-08 0 4.13169e-09 -2.94647e-08 0 2.80698e-09 -2.95023e-08 0 2.77916e-09 -2.95178e-08 0
+          2.78568e-09 -2.95297e-08 0 4.19425e-09 -2.95735e-08 0 5.59055e-09 -2.93762e-08 0
         </DataArray>
       </PointData>
       <CellData Scalars="process rank">
-- 
GitLab