[porousmediumflow][velocity] Improve calculation in case of Neumann BCs

• check if we have non-compositional model
• estimate the volume flow by neumannFlow/density

dumux/porousmediumflow/velocity.hh

@@ -26,6 +26,7 @@
 #define DUMUX_POROUSMEDIUMFLOW_VELOCITY_HH
 
 #include <vector>
+#include <type_traits>
 
 #include <dune/common/fvector.hh>
 #include <dune/common/float_cmp.hh>
@@ -38,6 +39,27 @@
 
 namespace Dumux {
 
+#ifndef DOXYGEN
+namespace Detail {
+// helper structs and functions detecting if the model is compositional
+
+template <typename T, typename ...Ts>
+using MoleFractionDetector = decltype(std::declval<T>().moleFraction(std::declval<Ts>()...));
+
+template<class T, typename ...Args>
+static constexpr bool hasMoleFraction()
+{ return Dune::Std::is_detected<MoleFractionDetector, T, Args...>::value; }
+
+template <typename T, typename ...Ts>
+using MassFractionDetector = decltype(std::declval<T>().massFraction(std::declval<Ts>()...));
+
+template<class T, typename ...Args>
+static constexpr bool hasMassFraction()
+{ return Dune::Std::is_detected<MassFractionDetector, T, Args...>::value; }
+
+} // end namespace Detail
+#endif // DOXYGEN
+
 /*!
  * \ingroup PorousmediumflowModels
  * \brief Velocity computation for implicit (porous media) models.
@@ -70,6 +92,11 @@ class PorousMediumFlowVelocity
 
     using Problem = typename GridVolumeVariables::Problem;
 
+    static constexpr bool modelIsCompositional = Detail::hasMoleFraction<typename GridVolumeVariables::VolumeVariables, int, int>() ||
+                                                 Detail::hasMassFraction<typename GridVolumeVariables::VolumeVariables, int, int>();
+
+    static_assert(VolumeVariables::numFluidPhases() >= 1, "Velocity output only makes sense for models with fluid phases.");
+
 public:
     static constexpr int numFluidPhases = VolumeVariables::numFluidPhases();
     using VelocityVector = std::vector<Dune::FieldVector<Scalar, dimWorld>>;
@@ -303,13 +330,26 @@ public:
                             if (Dune::FloatCmp::eq<NumEqVector, Dune::FloatCmp::CmpStyle::absolute>(neumannFlux, NumEqVector(0.0), 1e-30))
                                 scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0;
 
-                            // otherwise, we try some reconstruction (TODO: Can this be improved?)
+                            // otherwise, we try some reconstruction
                             // for cubes
                             else if (dim == 1 || geomType.isCube())
                             {
                                 const auto fIdx = scvfIndexInInside[localScvfIdx];
-                                const auto fIdxOpposite = fIdx%2 ? fIdx-1 : fIdx+1;
-                                scvfFluxes[fIdx] = -scvfFluxes[fIdxOpposite];
+
+                                if constexpr (!modelIsCompositional)
+                                {
+                                    // We assume that the density at the face equals the one at the cell center and reconstruct a volume flux from the Neumann mass flux.
+                                    const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+                                    const auto eqIdx = VolumeVariables::Indices::conti0EqIdx + phaseIdx;
+                                    scvfFluxes[fIdx] += neumannFlux[eqIdx] / insideVolVars.density(phaseIdx) * scvf.area() * insideVolVars.extrusionFactor();
+                                }
+                                else
+                                {
+                                    // For compositional models, we generally can't reconstruct the volume flow from the Neumann flux (which is a species flux rather
+                                    // than a phase flux here). Instead, we use the velocity of the opposing face.
+                                    const auto fIdxOpposite = fIdx%2 ? fIdx-1 : fIdx+1;
+                                    scvfFluxes[fIdx] = -scvfFluxes[fIdxOpposite];
+                                }
                             }
 
                             // for simplices