diff --git a/dumux/common/math.hh b/dumux/common/math.hh
index 8222a2ca9d92d9dfd9c333aae517891c0673931c..cc5eb1d58ac87311357f16dc56cd35e37b8b9eb6 100644
--- a/dumux/common/math.hh
+++ b/dumux/common/math.hh
@@ -553,6 +553,26 @@ Dune::FieldMatrix<Scalar, n, m> getTransposed(const Dune::FieldMatrix<Scalar, m,
return T;
}
+/*!
+ * \brief Transpose a DynamicMatrix
+ *
+ * \param M The matrix to be transposed
+ */
+template <class Scalar>
+Dune::DynamicMatrix<Scalar> getTransposed(const Dune::DynamicMatrix<Scalar>& M)
+{
+ std::size_t rows_T = M.M();
+ std::size_t cols_T = M.N();
+
+ Dune::DynamicMatrix<Scalar> M_T(rows_T, cols_T, 0.0);
+
+ for (std::size_t i = 0; i < rows_T; ++i)
+ for (std::size_t j = 0; j < cols_T; ++j)
+ M_T[i][j] = M[j][i];
+
+ return M_T;
+}
+
/*!
* \brief Multiply two dynamic matrices
*
@@ -577,6 +597,27 @@ Dune::DynamicMatrix<Scalar> multiplyMatrices(const Dune::DynamicMatrix<Scalar> &
return result;
}
+/*!
+ * \brief Trace of dynamic matrix
+ *
+ * \param M The dynamic matrix
+ */
+template <class Scalar>
+Scalar trace(const Dune::DynamicMatrix<Scalar>& M)
+{
+ std::size_t rows_T = M.M();
+ std::size_t cols_T = M.N();
+
+ DUNE_ASSERT_BOUNDS(rows_T == cols_T);
+
+ Scalar trace = 0.0;
+
+ for (std::size_t i = 0; i < rows_T; ++i)
+ trace += M[i][i];
+
+ return trace;
+}
+
} // end namespace Dumux
#endif
diff --git a/dumux/freeflow/staggered/indices.hh b/dumux/freeflow/staggered/indices.hh
index 9ff3adb4c10c970915da31a37fd9fc41c3a350d4..704590dac4c5db4a160bf7a1ea16c0bd31838cb7 100644
--- a/dumux/freeflow/staggered/indices.hh
+++ b/dumux/freeflow/staggered/indices.hh
@@ -50,7 +50,7 @@ struct NavierStokesCommonIndices
static constexpr int momentumXBalanceIdx = momentumBalanceIdx; //!< Index of the momentum balance equation
static constexpr int momentumYBalanceIdx = momentumBalanceIdx + 1; //!< Index of the momentum balance equation
static constexpr int momentumZBalanceIdx = momentumBalanceIdx + 2; //!< Index of the momentum balance equation
- static constexpr int velocityIdx = PVOffset + momentumBalanceIdx; //!< Index of the velocity in a solution vector (NOTE: This PV lives in a different vector than the pressure)
+ static constexpr int velocityIdx = momentumBalanceIdx; //!< Index of the velocity in a solution vector (NOTE: This PV lives in a different vector than the pressure)
static constexpr int velocityXIdx = velocityIdx; //!< Index of the velocity in a solution vector (NOTE: This PV lives in a different vector than the pressure)
static constexpr int velocityYIdx = velocityIdx + 1; //!< Index of the velocity in a solution vector (NOTE: This PV lives in a different vector than the pressure)
static constexpr int velocityZIdx = velocityIdx + 2; //!< Index of the velocity in a solution vector (NOTE: This PV lives in a different vector than the pressure)
diff --git a/dumux/freeflow/staggered/localresidual.hh b/dumux/freeflow/staggered/localresidual.hh
index c4006e0119ad701dfbbec6c98a554eff726c9b61..b5a440999c8b8edf478c3524fe469beb5b04c30a 100644
--- a/dumux/freeflow/staggered/localresidual.hh
+++ b/dumux/freeflow/staggered/localresidual.hh
@@ -127,10 +127,10 @@ public:
const ElementVolumeVariables& elemVolVars,
const GlobalFaceVars& globalFaceVars,
const SubControlVolumeFace &scvf,
- const FluxVariablesCache& fluxVarsCache)
+ const ElementFluxVariablesCache& elemFluxVarsCache)
{
FluxVariables fluxVars;
- return fluxVars.computeFluxForCellCenter(element, fvGeometry, elemVolVars, globalFaceVars, scvf, fluxVarsCache);
+ return fluxVars.computeFluxForCellCenter(element, fvGeometry, elemVolVars, globalFaceVars, scvf, elemFluxVarsCache[scvf]);
}
CellCenterPrimaryVariables computeSourceForCellCenter(const Element &element,
@@ -206,7 +206,8 @@ public:
const SubControlVolumeFace& scvf,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
- const GlobalFaceVars& globalFaceVars)
+ const GlobalFaceVars& globalFaceVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache)
{
FacePrimaryVariables flux(0.0);
FluxVariables fluxVars;
@@ -250,7 +251,7 @@ protected:
if (scvf.boundary())
{
// For the mass-balance residual, do the same as if the face was not on a boundary.This might need to be changed sometime...
- this->ccResidual_ += computeFluxForCellCenter(element, fvGeometry, elemVolVars, faceVars, scvf, elemFluxVarsCache[scvf]);
+ this->ccResidual_ += computeFluxForCellCenter(element, fvGeometry, elemVolVars, faceVars, scvf, elemFluxVarsCache);
// handle the actual boundary conditions:
const auto bcTypes = this->problem().boundaryTypes(element, scvf);
@@ -294,7 +295,7 @@ protected:
if(bcTypes.isOutflow(momentumBalanceIdx))
{
if(bcTypes.isDirichlet(massBalanceIdx))
- this->faceResiduals_[scvf.localFaceIdx()] += computeFluxForFace(element, scvf, fvGeometry, elemVolVars, faceVars);
+ this->faceResiduals_[scvf.localFaceIdx()] += computeFluxForFace(element, scvf, fvGeometry, elemVolVars, faceVars, elemFluxVarsCache);
else
DUNE_THROW(Dune::InvalidStateException, "Face at " << scvf.center() << " has an outflow BC for the momentum balance but no Dirichlet BC for the pressure!");
}
diff --git a/dumux/implicit/staggered/localjacobian.hh b/dumux/implicit/staggered/localjacobian.hh
index 28e12e50a926e85c189571866ed3874cdf107797..d54bcc9787cd4145528c2a832dc7851d6c75b960 100644
--- a/dumux/implicit/staggered/localjacobian.hh
+++ b/dumux/implicit/staggered/localjacobian.hh
@@ -193,7 +193,17 @@ public:
this->model_().updatePVWeights(fvGeometry);
// calculate derivatives of all dofs in stencil with respect to the dofs in the element
- evalPartialDerivatives_(element, fvGeometry, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, residual[cellCenterIdx][ccGlobalI_], faceResidualCache);
+ evalPartialDerivatives_(element,
+ fvGeometry,
+ prevElemVolVars,
+ curElemVolVars,
+ prevGlobalFaceVars,
+ curGlobalFaceVars,
+ elemFluxVarsCache,
+ elemBcTypes,
+ matrix,
+ residual[cellCenterIdx][ccGlobalI_],
+ faceResidualCache);
}
/*!
diff --git a/dumux/implicit/staggered/localresidual.hh b/dumux/implicit/staggered/localresidual.hh
index a2d11b36c628f627be5fc3d58eb9164ae5864aa1..deb06209fb725b43c1d105ad717c9795bbe44acc 100644
--- a/dumux/implicit/staggered/localresidual.hh
+++ b/dumux/implicit/staggered/localresidual.hh
@@ -57,7 +57,6 @@ class StaggeredLocalResidual
using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
- using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
@@ -297,7 +296,7 @@ protected:
for (auto&& scvf : scvfs(fvGeometry))
{
if(!scvf.boundary())
- ccResidual_ += asImp_().computeFluxForCellCenter(element, fvGeometry, elemVolVars, faceVars, scvf, elemFluxVarsCache[scvf]);
+ ccResidual_ += asImp_().computeFluxForCellCenter(element, fvGeometry, elemVolVars, faceVars, scvf, elemFluxVarsCache);
}
}
@@ -313,7 +312,7 @@ protected:
const ElementFluxVariablesCache& elemFluxVarsCache)
{
if(!scvf.boundary())
- faceResiduals_[scvf.localFaceIdx()] += asImp_().computeFluxForFace(element, scvf, fvGeometry, elemVolVars, globalFaceVars);
+ faceResiduals_[scvf.localFaceIdx()] += asImp_().computeFluxForFace(element, scvf, fvGeometry, elemVolVars, globalFaceVars, elemFluxVarsCache);
}
/*!
diff --git a/dumux/implicit/staggered/propertydefaults.hh b/dumux/implicit/staggered/propertydefaults.hh
index 3917d92c883135b3ed37806e885dd6c21d9f310b..5e5f215a59dba5b9dec31177d9d94d049268fe5e 100644
--- a/dumux/implicit/staggered/propertydefaults.hh
+++ b/dumux/implicit/staggered/propertydefaults.hh
@@ -49,6 +49,8 @@
#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+
#include "assembler.hh"
#include "localresidual.hh"
#include "localjacobian.hh"
@@ -246,6 +248,7 @@ SET_BOOL_PROP(StaggeredModel, VtkWriteFaceData, true);
//! For compatibility
SET_BOOL_PROP(StaggeredModel, EnableInteriorBoundaries, false);
+SET_TYPE_PROP(StaggeredModel, VtkOutputModule, StaggeredVtkOutputModule<TypeTag>);
} // namespace Properties