From 27e9c60ef64500d9b7fa8b9f67efda03cfdb07c0 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Mon, 5 Dec 2016 16:49:22 +0100
Subject: [PATCH] [staggeredGrid][localJacobian] Add derivative of cc dofs
w.r.t face dofs
*Same way as derivative of cc dofs w.r.t cc dofs
---
dumux/implicit/staggered/localjacobian.hh | 59 ++++++++++++++++++++++-
1 file changed, 58 insertions(+), 1 deletion(-)
diff --git a/dumux/implicit/staggered/localjacobian.hh b/dumux/implicit/staggered/localjacobian.hh
index 3302a894c8..7a2cbd3b0e 100644
--- a/dumux/implicit/staggered/localjacobian.hh
+++ b/dumux/implicit/staggered/localjacobian.hh
@@ -193,8 +193,12 @@ private:
// compute the derivatives of the cell center dofs with respect to cell center dofs
dCCdCC_(element, fvGeometry, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, residual, isGhost);
+ // compute the derivatives of the cell center dofs with respect to face dofs
+ dCCdFace_(element, fvGeometry, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, residual, isGhost);
+
printmatrix(std::cout, matrix[cellCenterIdx][cellCenterIdx], "A11 neu", "");
+ printmatrix(std::cout, matrix[cellCenterIdx][faceIdx], "A12 neu", "");
}
/*!
@@ -215,7 +219,6 @@ private:
// set the actual dof index
auto globalI = this->problem_().elementMapper().index(element);
-
// build derivatives with for cell center dofs w.r.t. cell center dofs
const auto& cellCenterToCellCenterStencil = this->model_().stencils(element).cellCenterToCellCenterStencil();
@@ -255,6 +258,60 @@ private:
}
}
+ /*!
+ * \brief Computes the derivatives of the cell center dofs with respect to face dofs
+ */
+ void dCCdFace_(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& prevElemVolVars,
+ const ElementVolumeVariables& curElemVolVars,
+ const GlobalFaceVars& prevGlobalFaceVars,
+ GlobalFaceVars& curGlobalFaceVars,
+ ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes& elemBcTypes,
+ JacobianMatrix& matrix,
+ SolutionVector& residual,
+ const bool isGhost)
+ {
+ // set the actual dof index
+ auto globalI = this->problem_().elementMapper().index(element);
+
+ // build derivatives with for cell center dofs w.r.t. face dofs
+ const auto& cellCenterToFaceStencil = this->model_().stencils(element).cellCenterToFaceStencil();
+
+ for(const auto& globalJ : cellCenterToFaceStencil)
+ {
+ // get the faceVars of the face with respect to which we are going to build the derivative
+
+ auto origFaceVars = curGlobalFaceVars.faceVars(globalJ);
+ auto& curFaceVars = curGlobalFaceVars.faceVars(globalJ);
+
+ FacePrimaryVariables priVars(this->model_().curSol()[faceIdx][globalJ]);
+
+ for(int pvIdx = 0; pvIdx < priVars.size(); ++pvIdx)
+ {
+ const Scalar eps = 1e-4; // TODO: do properly
+ priVars += eps;
+
+ curFaceVars.update(priVars);
+
+ this->localResidual().eval(element, fvGeometry,
+ prevElemVolVars, curElemVolVars,
+ prevGlobalFaceVars, curGlobalFaceVars,
+ elemBcTypes, elemFluxVarsCache);
+
+ auto partialDeriv = (this->localResidual().ccResidual() - residual[cellCenterIdx][globalI]);
+ partialDeriv /= eps;
+
+ // update the global jacobian matrix with the current partial derivatives
+ this->updateGlobalJacobian_(matrix[cellCenterIdx][faceIdx], globalI, globalJ, pvIdx, partialDeriv);
+
+ // restore the original faceVars
+ curFaceVars = origFaceVars;
+ }
+ }
+ }
+
/*!
* \brief Updates the current global Jacobian matrix with the
* partial derivatives of all equations in regard to the
--
GitLab