diff --git a/dumux/freeflow/staggered/localresidual.hh b/dumux/freeflow/staggered/localresidual.hh
index b5a440999c8b8edf478c3524fe469beb5b04c30a..e6210b546e063233ce84cb7175d56375439d099f 100644
--- a/dumux/freeflow/staggered/localresidual.hh
+++ b/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;
diff --git a/dumux/freeflow/staggered/properties.hh b/dumux/freeflow/staggered/properties.hh
index b516c753178967ca07cbd4b09104f149a437fca4..3d9a591cd30a3765b81e9f4adc09d1847673a5cc 100644
--- a/dumux/freeflow/staggered/properties.hh
+++ b/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
// \}
}
diff --git a/dumux/freeflow/staggered/propertydefaults.hh b/dumux/freeflow/staggered/propertydefaults.hh
index dfa075570fe407e087b63adae964b4b0d7664f9f..9852a5e6b2db9161a2db8fff1d15dc7bf2a3134b 100644
--- a/dumux/freeflow/staggered/propertydefaults.hh
+++ b/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:
diff --git a/dumux/implicit/staggered/localjacobian.hh b/dumux/implicit/staggered/localjacobian.hh
index d54bcc9787cd4145528c2a832dc7851d6c75b960..d2da798d54d4adc7ae07e2eac54e7651a12d8760 100644
--- a/dumux/implicit/staggered/localjacobian.hh
+++ b/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)};
};
}
diff --git a/dumux/implicit/staggered/properties.hh b/dumux/implicit/staggered/properties.hh
index 1c49343c05c80d9e369e093cdbd69a268b8b4c9e..30f49ad0d214bf5f40b240c2bc2c85deb1303640 100644
--- a/dumux/implicit/staggered/properties.hh
+++ b/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
diff --git a/dumux/implicit/staggered/propertydefaults.hh b/dumux/implicit/staggered/propertydefaults.hh
index 5e5f215a59dba5b9dec31177d9d94d049268fe5e..d9c8d6679a4cff44774389662a886d8d74dea559 100644
--- a/dumux/implicit/staggered/propertydefaults.hh
+++ b/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