Merge branch 'feature/improve-staggered-localJacobian' into 'next'

Feature/improve staggered local jacobian

See merge request !388

dumux/freeflow/staggered/localresidual.hh

@@ -114,6 +114,7 @@ class StaggeredNavierStokesResidualImpl<TypeTag, false, false> : public Dumux::S
     using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
 
     static constexpr bool navierStokes = GET_PROP_VALUE(TypeTag, EnableInertiaTerms);
+    static constexpr bool normalizePressure = GET_PROP_VALUE(TypeTag, NormalizePressure);
 
 public:
     // copying the local residual class is not a good idea
@@ -321,12 +322,14 @@ private:
         const auto insideScvIdx = scvf.insideScvIdx();
         const auto& insideVolVars = elemVolVars[insideScvIdx];
 
-        Scalar result = insideVolVars.pressure() * scvf.area() * -1.0 * sign(scvf.outerNormalScalar());
+        const Scalar deltaP = normalizePressure ? this->problem().initialAtPos(scvf.center())[cellCenterIdx][pressureIdx] : 0.0;
+
+        Scalar result = (insideVolVars.pressure() - deltaP) * scvf.area() * -1.0 * sign(scvf.outerNormalScalar());
 
         // treat outflow BCs
         if(scvf.boundary())
         {
-            const Scalar pressure = this->problem().dirichletAtPos(scvf.center())[cellCenterIdx][pressureIdx];
+            const Scalar pressure = this->problem().dirichletAtPos(scvf.center())[cellCenterIdx][pressureIdx] - deltaP;
             result += pressure * scvf.area() * sign(scvf.outerNormalScalar());
         }
         return result;

dumux/freeflow/staggered/properties.hh

@@ -65,6 +65,7 @@ NEW_PROP_TAG(BoundaryValues); //!< Type to set values on the boundary
 NEW_PROP_TAG(EnableComponentTransport); //!< Returns whether to consider component transport or not
 NEW_PROP_TAG(EnableEnergyTransport); //!<  Returns whether to consider energy transport or not
 NEW_PROP_TAG(FaceVariables); //!<  Returns whether to consider energy transport or not
+NEW_PROP_TAG(NormalizePressure); //!<  Returns whether to normalize the pressure term in the momentum balance or not
 // \}
 }
 

dumux/freeflow/staggered/propertydefaults.hh

@@ -166,6 +166,9 @@ SET_BOOL_PROP(NavierStokes, EnableEnergyTransport, false);
 
 SET_BOOL_PROP(NavierStokes, EnableComponentTransport, false);
 
+//! Normalize the pressure term in the momentum balance or not
+SET_BOOL_PROP(NavierStokes, NormalizePressure, true);
+
 SET_PROP(NavierStokes, BoundaryValues)
 {
 private:

dumux/implicit/staggered/localjacobian.hh

@@ -279,7 +279,7 @@ private:
                 PrimaryVariables priVars(CellCenterPrimaryVariables(this->model_().curSol()[cellCenterIdx][globalJ]),
                                          FacePrimaryVariables(0.0));
 
-                const Scalar eps = numericEpsilon(priVars[pvIdx]);
+                const Scalar eps = numericEpsilon(priVars[pvIdx], cellCenterIdx, cellCenterIdx);
                 priVars[pvIdx] += eps;
                 ElementSolutionVector elemSol{std::move(priVars)};
                 curVolVars.update(elemSol, this->problem_(), elementJ, scvJ);
@@ -327,7 +327,7 @@ private:
             for(auto pvIdx : PriVarIndices(faceIdx))
             {
                 PrimaryVariables priVars(CellCenterPrimaryVariables(0.0), FacePrimaryVariables(this->model_().curSol()[faceIdx][globalJ]));
-                const Scalar eps = numericEpsilon(priVars[pvIdx]);
+                const Scalar eps = numericEpsilon(priVars[pvIdx], cellCenterIdx, faceIdx);
                 priVars[pvIdx] += eps;
                 curFaceVars.update(priVars[faceIdx]);
 
@@ -381,7 +381,7 @@ private:
                     PrimaryVariables priVars(CellCenterPrimaryVariables(this->model_().curSol()[cellCenterIdx][globalJ]),
                                              FacePrimaryVariables(0.0));
 
-                    const Scalar eps = numericEpsilon(priVars[pvIdx]);
+                    const Scalar eps = numericEpsilon(priVars[pvIdx], faceIdx, cellCenterIdx);
                     priVars[pvIdx] += eps;
                     ElementSolutionVector elemSol{std::move(priVars)};
                     curVolVars.update(elemSol, this->problem_(), elementJ, scvJ);
@@ -433,7 +433,7 @@ private:
                 {
                     PrimaryVariables priVars(CellCenterPrimaryVariables(0.0), FacePrimaryVariables(this->model_().curSol()[faceIdx][globalJ]));
 
-                    const Scalar eps = numericEpsilon(priVars[pvIdx]);
+                    const Scalar eps = numericEpsilon(priVars[pvIdx], faceIdx, faceIdx);
                     priVars[pvIdx] += eps;
                     curFaceVars.update(priVars[faceIdx]);
 
@@ -497,7 +497,14 @@ private:
     const JacobianAssembler &jacAsm_() const
     { return model_().jacobianAssembler(); }
 
-    Scalar numericEpsilon(const Scalar priVar) const
+     /*!
+     * \brief Return the epsilon used to calculate the numeric derivative of a localResidual w.r.t a certain priVar
+     *
+     * \param priVar The the value of primary varible w.r.t which to derivative of the localResidual is calculated
+     * \param idx1 Indicates whether the the derivative is build for a cellCenter or face localResidual
+     * \param idx2 Indicates whether the the derivative is build w.r.t a priVar living on a cellCenter or face
+     */
+    Scalar numericEpsilon(const Scalar priVar, const int idx1, const int idx2) const
     {
         // define the base epsilon as the geometric mean of 1 and the
         // resolution of the scalar type. E.g. for standard 64 bit
@@ -507,12 +514,14 @@ private:
         static const Scalar baseEps
             = Dumux::geometricMean<Scalar>(std::numeric_limits<Scalar>::epsilon(), 1.0);
         */
-        static const Scalar baseEps = 1e-10;
+
+        static const Scalar baseEps = baseEps_[idx1][idx2];
         assert(std::numeric_limits<Scalar>::epsilon()*1e4 < baseEps);
         // the epsilon value used for the numeric differentiation is
         // now scaled by the absolute value of the primary variable...
         return baseEps*(std::abs(priVar) + 1.0);
     }
+
     /*!
      * \brief Updates the current global Jacobian matrix with the
      *        partial derivatives of all equations in regard to the
@@ -561,6 +570,8 @@ private:
     LocalResidual localResidual_;
 
     AssemblyMap assemblyMap_;
+
+    static constexpr auto baseEps_{GET_PROP_VALUE(TypeTag, BaseEpsilon)};
 };
 
 }

dumux/implicit/staggered/properties.hh

@@ -52,6 +52,8 @@ NEW_PROP_TAG(FaceSolutionVector); //!< Vector containing all face primary variab
 
 NEW_PROP_TAG(EnableInteriorBoundaries); //!< For compatibility
 
+NEW_PROP_TAG(BaseEpsilon); //!< Set one or different base epsilons for the calculations of the localJacobian's derivatives
+
 NEW_PROP_TAG(NumEqCellCenter); //!< Number of equations per cell center dof
 NEW_PROP_TAG(NumEqFace); //!< Number of equations per face dof
 NEW_PROP_TAG(DofTypeIndices); //!< Indices to choose between cell center and face dofs

dumux/implicit/staggered/propertydefaults.hh

@@ -248,8 +248,24 @@ SET_BOOL_PROP(StaggeredModel, VtkWriteFaceData, true);
 //! For compatibility
 SET_BOOL_PROP(StaggeredModel, EnableInteriorBoundaries, false);
 
+//! Set staggered vtk output module
 SET_TYPE_PROP(StaggeredModel, VtkOutputModule, StaggeredVtkOutputModule<TypeTag>);
 
+//! Set one or different base epsilons for the calculations of the localJacobian's derivatives
+SET_PROP(StaggeredModel, BaseEpsilon)
+{
+private:
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    static constexpr Scalar dCCdCC = 1e-8;
+    static constexpr Scalar dCCdFace = 1e-8;
+    static constexpr Scalar dFacedCC = 1e-8;
+    static constexpr Scalar dFacedFace = 1e-8;
+
+public:
+    static constexpr std::array<std::array<Scalar, 2>, 2> value{{{dCCdCC,   dCCdFace},
+                                                                 {dFacedCC, dFacedFace}}};
+};
+
 } // namespace Properties
 
 } // namespace Dumux