From 1f9b5d5652f597303e8903780ce1da98400dec7b Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Mon, 29 Jun 2020 14:19:41 +0200
Subject: [PATCH] [cclocalassembler] Implement internal Dirichlet constraints
---
dumux/assembly/cclocalassembler.hh | 124 ++++++++++++++++++++++++++---
1 file changed, 114 insertions(+), 10 deletions(-)
diff --git a/dumux/assembly/cclocalassembler.hh b/dumux/assembly/cclocalassembler.hh
index 004bb7e8f0..37e51257d3 100644
--- a/dumux/assembly/cclocalassembler.hh
+++ b/dumux/assembly/cclocalassembler.hh
@@ -63,8 +63,6 @@ class CCLocalAssemblerBase : public FVLocalAssemblerBase<TypeTag, Assembler, Imp
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
- static_assert(!Assembler::Problem::enableInternalDirichletConstraints(), "Internal Dirichlet constraints are currently not implemented for cc-methods!");
-
public:
using ParentType::ParentType;
@@ -141,6 +139,7 @@ class CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/tru
using FVElementGeometry = typename GridGeometry::LocalView;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
+ using Problem = typename GridVariables::GridVolumeVariables::Problem;
enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };
enum { dim = GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimension };
@@ -275,15 +274,54 @@ public:
// add the current partial derivatives to the global jacobian matrix
// no special treatment is needed if globalJ is a ghost because then derivatives have been assembled to 0 above
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ if constexpr (Problem::enableInternalDirichletConstraints())
{
- // the diagonal entries
- A[globalI][globalI][eqIdx][pvIdx] += partialDerivs[0][eqIdx];
+ // check if own element has internal Dirichlet constraint
+ const auto internalDirichletConstraintsOwnElement = this->problem().hasInternalDirichletConstraint(this->element(), scv);
+ const auto dirichletValues = this->problem().internalDirichlet(this->element(), scv);
+
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (internalDirichletConstraintsOwnElement[eqIdx])
+ {
+ origResiduals[0][eqIdx] = this->curElemVolVars()[scv].priVars()[eqIdx] - dirichletValues[eqIdx];
+ A[globalI][globalI][eqIdx][pvIdx] = (eqIdx == pvIdx) ? 1.0 : 0.0;
+ }
+ else
+ A[globalI][globalI][eqIdx][pvIdx] += partialDerivs[0][eqIdx];
+ }
// off-diagonal entries
j = 1;
for (const auto& dataJ : connectivityMap[globalI])
- A[dataJ.globalJ][globalI][eqIdx][pvIdx] += partialDerivs[j++][eqIdx];
+ {
+ const auto& neighborElement = neighborElements[j-1];
+ const auto& neighborScv = fvGeometry.scv(dataJ.globalJ);
+ const auto internalDirichletConstraintsNeighbor = this->problem().hasInternalDirichletConstraint(neighborElement, neighborScv);
+
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (internalDirichletConstraintsNeighbor[eqIdx])
+ A[dataJ.globalJ][globalI][eqIdx][pvIdx] = 0.0;
+ else
+ A[dataJ.globalJ][globalI][eqIdx][pvIdx] += partialDerivs[j][eqIdx];
+ }
+
+ ++j;
+ }
+ }
+ else // no internal Dirichlet constraints specified
+ {
+ for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ {
+ // the diagonal entries
+ A[globalI][globalI][eqIdx][pvIdx] += partialDerivs[0][eqIdx];
+
+ // off-diagonal entries
+ j = 1;
+ for (const auto& dataJ : connectivityMap[globalI])
+ A[dataJ.globalJ][globalI][eqIdx][pvIdx] += partialDerivs[j++][eqIdx];
+ }
}
// restore the original state of the scv's volume variables
@@ -325,6 +363,7 @@ class CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/fal
using Element = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView::template Codim<0>::Entity;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
+ using Problem = typename GridVariables::GridVolumeVariables::Problem;
enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };
@@ -343,7 +382,7 @@ public:
DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
// assemble the undeflected residual
- const auto residual = this->evalLocalResidual()[0];
+ auto residual = this->evalLocalResidual()[0];
const auto storageResidual = this->evalLocalStorageResidual()[0];
//////////////////////////////////////////////////////////////////////////////////////////////////
@@ -405,8 +444,28 @@ public:
// add the current partial derivatives to the global jacobian matrix
// only diagonal entries for explicit jacobians
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
- A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ if constexpr (Problem::enableInternalDirichletConstraints())
+ {
+ // check if own element has internal Dirichlet constraint
+ const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv);
+ const auto dirichletValues = this->problem().internalDirichlet(this->element(), scv);
+
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (internalDirichletConstraints[eqIdx])
+ {
+ residual[eqIdx] = this->curElemVolVars()[scv].priVars()[eqIdx] - dirichletValues[eqIdx];
+ A[globalI][globalI][eqIdx][pvIdx] = (eqIdx == pvIdx) ? 1.0 : 0.0;
+ }
+ else
+ A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ }
+ }
+ else
+ {
+ for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ }
// restore the original state of the scv's volume variables
curVolVars = origVolVars;
@@ -435,6 +494,7 @@ class CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/tr
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ using Problem = typename GridVariables::GridVolumeVariables::Problem;
enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };
@@ -461,7 +521,7 @@ public:
}
// assemble the undeflected residual
- const auto residual = this->evalLocalResidual()[0];
+ auto residual = this->evalLocalResidual()[0];
// get some aliases for convenience
const auto& problem = this->problem();
@@ -507,6 +567,30 @@ public:
}
}
+ if constexpr (Problem::enableInternalDirichletConstraints())
+ {
+ // check if own element has internal Dirichlet constraint
+ const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv);
+ const auto dirichletValues = this->problem().internalDirichlet(this->element(), scv);
+
+ for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
+ {
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (internalDirichletConstraints[eqIdx])
+ {
+ residual[eqIdx] = this->curElemVolVars()[scv].priVars()[eqIdx] - dirichletValues[eqIdx];
+ A[globalI][globalI][eqIdx][pvIdx] = (eqIdx == pvIdx) ? 1.0 : 0.0;
+
+ // inner faces
+ for (const auto& scvf : scvfs(fvGeometry))
+ if (!scvf.boundary())
+ A[globalI][fvGeometry.scv(scvf.outsideScvIdx()).dofIndex()][eqIdx][pvIdx] = 0.0;
+ }
+ }
+ }
+ }
+
// return element residual
return residual;
}
@@ -527,6 +611,7 @@ class CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/fa
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ using Problem = typename GridVariables::GridVolumeVariables::Problem;
enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };
@@ -563,6 +648,25 @@ public:
// add hand-code derivative of storage term
this->localResidual().addStorageDerivatives(A[globalI][globalI], this->problem(), this->element(), this->fvGeometry(), volVars, scv);
+ if constexpr (Problem::enableInternalDirichletConstraints())
+ {
+ // check if own element has internal Dirichlet constraint
+ const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv);
+ const auto dirichletValues = this->problem().internalDirichlet(this->element(), scv);
+
+ for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
+ {
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (internalDirichletConstraints[eqIdx])
+ {
+ residual[eqIdx] = this->curElemVolVars()[scv].priVars()[eqIdx] - dirichletValues[eqIdx];
+ A[globalI][globalI][eqIdx][pvIdx] = (eqIdx == pvIdx) ? 1.0 : 0.0;
+ }
+ }
+ }
+ }
+
// return the original residual
return residual;
}
--
GitLab