Skip to content
Snippets Groups Projects
Commit f88374d6 authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'feature/add-coupled-assembly-for-box' into 'next' 

Feature/add coupled assembly for box

See merge request !433
parents b41e2628 00c478f1
No related branches found
No related tags found
2 merge requests!617[WIP] Next,!433Feature/add coupled assembly for box
Hide whitespace changes
Inline Side-by-side
dumux/mixeddimension/assembler.hh
+ 66
− 64
View file @ f88374d6
......@@ -182,7 +182,8 @@ protected:
residual_[lowDimIdx]);
}
void buildBulkMatrixBlocksCC_(BulkMatrixBlock& m, BulkCouplingMatrixBlock& cm)
template <class T = BulkProblemTypeTag, typename std::enable_if<((!GET_PROP_VALUE(T, ImplicitIsBox))), bool>::type = 0>
void buildBulkMatrixBlocks_(BulkMatrixBlock& m, BulkCouplingMatrixBlock& cm)
{
// get occupation pattern of the matrix
Dune::MatrixIndexSet bulkPattern, bulkCouplingPattern;
......@@ -213,7 +214,8 @@ protected:
bulkCouplingPattern.exportIdx(cm);
}
void buildLowDimMatrixBlocksCC_(LowDimMatrixBlock& m, LowDimCouplingMatrixBlock& cm)
template <class T = LowDimProblemTypeTag, typename std::enable_if<((!GET_PROP_VALUE(T, ImplicitIsBox))), bool>::type = 0>
void buildLowDimMatrixBlocks_(LowDimMatrixBlock& m, LowDimCouplingMatrixBlock& cm)
{
// get occupation pattern of the matrix
Dune::MatrixIndexSet lowDimPattern, lowDimCouplingPattern;
......@@ -244,63 +246,73 @@ protected:
lowDimCouplingPattern.exportIdx(cm);
}
// void buildBulkMatrixBlocksBox_(BulkMatrixBlock& m, BulkCouplingMatrixBlock& cm)
// {
// // get occupation pattern of the matrix
// Dune::MatrixIndexSet bulkPattern, bulkCouplingPattern;
// bulkPattern.resize(m.N(), m.M());
// bulkCouplingPattern.resize(cm.N(), cm.M());
// for (const auto& element : elements(problem_().bulkGridView()))
// {
// const auto& stencil = problem_().couplingManager().stencil(element);
// const auto& couplingStencil = problem_().couplingManager().couplingStencil(element);
// for (unsigned int scvIdx = 0; scvIdx < element.subEntities(bulkDim); ++scvIdx)
// {
// auto globalI = problem_().model().bulkVertexMapper().subIndex(element, scvIdx, bulkDim);
// for (auto&& globalJ : stencil)
// bulkPattern.add(globalI, globalJ);
// for (auto&& globalJ : couplingStencil)
// bulkCouplingPattern.add(globalI, globalJ);
// }
// }
template <class T = BulkProblemTypeTag, typename std::enable_if<((GET_PROP_VALUE(T, ImplicitIsBox))), bool>::type = 0>
void buildBulkMatrixBlocks_(BulkMatrixBlock& m, BulkCouplingMatrixBlock& cm)
{
// get occupation pattern of the matrix
Dune::MatrixIndexSet bulkPattern, bulkCouplingPattern;
bulkPattern.resize(m.N(), m.M());
bulkCouplingPattern.resize(cm.N(), cm.M());
// // export occupation patterns to matrices
// bulkPattern.exportIdx(m);
// bulkCouplingPattern.exportIdx(cm);
// }
for (const auto& element : elements(problem_().bulkGridView()))
{
const auto& couplingStencil = problem_().couplingManager().couplingStencil(element);
// void buildLowDimMatrixBlocksBox_(LowDimMatrixBlock& m, LowDimCouplingMatrixBlock& cm)
// {
// // get occupation pattern of the matrix
// Dune::MatrixIndexSet lowDimPattern, lowDimCouplingPattern;
// lowDimPattern.resize(m.N(), m.M());
// lowDimCouplingPattern.resize(cm.N(), cm.M());
for (unsigned int vIdx = 0; vIdx < element.subEntities(bulkDim); ++vIdx)
{
const auto globalI = problem_().model().bulkVertexMapper().subIndex(element, vIdx, bulkDim);
for (unsigned int vIdx2 = vIdx; vIdx2 < element.subEntities(bulkDim); ++vIdx2)
{
const auto globalJ = problem_().model().bulkVertexMapper().subIndex(element, vIdx2, bulkDim);
bulkPattern.add(globalI, globalJ);
bulkPattern.add(globalJ, globalI);
}
for (auto&& globalJ : couplingStencil)
bulkCouplingPattern.add(globalI, globalJ);
//TODO: additionalDofDependencies
}
}
// for (const auto& element : elements(problem_().lowDimGridView()))
// {
// const auto& stencil = problem_().couplingManager().stencil(element);
// const auto& couplingStencil = problem_().couplingManager().couplingStencil(element);
// export occupation patterns to matrices
bulkPattern.exportIdx(m);
bulkCouplingPattern.exportIdx(cm);
}
// for (unsigned int scvIdx = 0; scvIdx < element.subEntities(lowDimDim); ++scvIdx)
// {
// auto globalI = problem_().model().lowDimVertexMapper().subIndex(element, scvIdx, lowDimDim);
template <class T = LowDimProblemTypeTag, typename std::enable_if<((GET_PROP_VALUE(T, ImplicitIsBox))), bool>::type = 0>
void buildLowDimMatrixBlocks_(LowDimMatrixBlock& m, LowDimCouplingMatrixBlock& cm)
{
// get occupation pattern of the matrix
Dune::MatrixIndexSet lowDimPattern, lowDimCouplingPattern;
lowDimPattern.resize(m.N(), m.M());
lowDimCouplingPattern.resize(cm.N(), cm.M());
// for (auto&& globalJ : stencil)
// lowDimPattern.add(globalI, globalJ);
for (const auto& element : elements(problem_().lowDimGridView()))
{
const auto& couplingStencil = problem_().couplingManager().couplingStencil(element);
// for (auto&& globalJ : couplingStencil)
// lowDimCouplingPattern.add(globalI, globalJ);
// }
// }
for (unsigned int vIdx = 0; vIdx < element.subEntities(lowDimDim); ++vIdx)
{
const auto globalI = problem_().model().lowDimVertexMapper().subIndex(element, vIdx, lowDimDim);
for (unsigned int vIdx2 = vIdx; vIdx2 < element.subEntities(lowDimDim); ++vIdx2)
{
const auto globalJ = problem_().model().lowDimVertexMapper().subIndex(element, vIdx2, lowDimDim);
lowDimPattern.add(globalI, globalJ);
lowDimPattern.add(globalJ, globalI);
}
for (auto&& globalJ : couplingStencil)
lowDimCouplingPattern.add(globalI, globalJ);
//TODO: additionalDofDependencies
}
}
// // export occupation patterns to matrices
// lowDimPattern.exportIdx(m);
// lowDimCouplingPattern.exportIdx(cm);
// }
// export occupation patterns to matrices
lowDimPattern.exportIdx(m);
lowDimCouplingPattern.exportIdx(cm);
}
// Construct the multitype matrix for the global jacobian
void createMatrix_()
......@@ -318,18 +330,8 @@ protected:
auto A22 = LowDimMatrixBlock(lowDimSize, lowDimSize, LowDimMatrixBlock::random);
auto A21 = LowDimCouplingMatrixBlock(lowDimSize, bulkSize, LowDimCouplingMatrixBlock::random);
// cell-centered
if (!bulkIsBox)
buildBulkMatrixBlocksCC_(A11, A12);
else{
// buildBulkMatrixBlocksBox_(A11, A12);
}
if (!lowDimIsBox)
buildLowDimMatrixBlocksCC_(A22, A21);
else{
// buildLowDimMatrixBlocksBox_(A22, A21);
}
buildBulkMatrixBlocks_(A11, A12);
buildLowDimMatrixBlocks_(A22, A21);
(*matrix_)[bulkIdx][bulkIdx] = A11;
(*matrix_)[bulkIdx][lowDimIdx] = A12;
......
dumux/mixeddimension/subproblemlocaljacobian.hh
+ 188
− 6
View file @ f88374d6
......@@ -104,6 +104,8 @@ class SubProblemLocalJacobian : public GET_PROP_TYPE(SubProblemTypeTag, LocalJac
using Element = typename GridView::template Codim<0>::Entity;
using SolutionVector = typename GET_PROP_TYPE(SubProblemTypeTag, SolutionVector);
enum { isBox = GET_PROP_VALUE(SubProblemTypeTag, ImplicitIsBox) };
public:
// copying a local jacobian is not a good idea
......@@ -138,8 +140,6 @@ public:
JacobianMatrixCoupling& couplingMatrix,
SolutionVector& residual)
{
const bool isGhost = (element.partitionType() == Dune::GhostEntity);
// prepare the volvars/fvGeometries in case caching is disabled
auto fvGeometry = localView(this->model_().globalFvGeometry());
fvGeometry.bind(element);
......@@ -153,13 +153,32 @@ public:
auto elemFluxVarsCache = localView(this->model_().globalFluxVarsCache());
elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
// set the actual dof index
this->globalI_ = this->problem_().elementMapper().index(element);
// check for boundaries on the element
ElementBoundaryTypes elemBcTypes;
elemBcTypes.update(this->problem_(), element, fvGeometry);
assemble_(element, fvGeometry, prevElemVolVars, curElemVolVars, elemFluxVarsCache, elemBcTypes, matrix, couplingMatrix, residual);
}
protected:
// for cell-centered models
template<class JacobianMatrix, class JacobianMatrixCoupling, class T = SubProblemTypeTag, typename std::enable_if<((!GET_PROP_VALUE(T, ImplicitIsBox))), bool>::type = 0>
void assemble_(const Element& element,
const FVElementGeometry& fvGeometry,
ElementVolumeVariables& prevElemVolVars,
ElementVolumeVariables& curElemVolVars,
ElementFluxVariablesCache& elemFluxVarsCache,
const ElementBoundaryTypes& elemBcTypes,
JacobianMatrix& matrix,
JacobianMatrixCoupling& couplingMatrix,
SolutionVector& residual)
{
const bool isGhost = (element.partitionType() == Dune::GhostEntity);
// set the actual dof index
this->globalI_ = this->problem_().elementMapper().index(element);
// Evaluate the undeflected element local residual
this->localResidual().eval(element,
fvGeometry,
......@@ -204,7 +223,73 @@ public:
residual);
}
protected:
// for box models
template<class JacobianMatrix, class JacobianMatrixCoupling, class T = SubProblemTypeTag, typename std::enable_if<((GET_PROP_VALUE(T, ImplicitIsBox))), bool>::type = 0>
void assemble_(const Element& element,
const FVElementGeometry& fvGeometry,
ElementVolumeVariables& prevElemVolVars,
ElementVolumeVariables& curElemVolVars,
ElementFluxVariablesCache& elemFluxVarsCache,
const ElementBoundaryTypes& elemBcTypes,
JacobianMatrix& matrix,
JacobianMatrixCoupling& couplingMatrix,
SolutionVector& residual)
{
constexpr auto numEq = GET_PROP_VALUE(T, NumEq);
// the element solution
auto curElemSol = this->problem_().model().elementSolution(element, this->problem_().model().curSol());
// calculate the actual element residual
this->localResidual().eval(element, fvGeometry, prevElemVolVars, curElemVolVars, elemBcTypes, elemFluxVarsCache);
this->residual_ = this->localResidual().residual();
// residual[this->globalI_] = this->localResidual().residual(0);
this->problem_().model().updatePVWeights(fvGeometry);
// calculation of the derivatives
for (auto&& scv : scvs(fvGeometry))
{
// dof index and corresponding actual pri vars
const auto dofIdx = scv.dofIndex();
auto& curVolVars = this->getCurVolVars(curElemVolVars, scv);
VolumeVariables origVolVars(curVolVars);
// add precalculated residual for this scv into the global container
residual[dofIdx] += this->residual_[scv.index()];
// calculate derivatives w.r.t to the privars at the dof at hand
for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
{
this->evalPartialDerivative_(matrix,
element,
fvGeometry,
prevElemVolVars,
curElemVolVars,
curElemSol,
scv,
elemBcTypes,
elemFluxVarsCache,
pvIdx);
// restore the original state of the scv's volume variables
curVolVars = origVolVars;
// restore the original element solution
curElemSol[scv.index()][pvIdx] = this->problem_().model().curSol()[scv.dofIndex()][pvIdx];
}
// TODO: what if we have an extended source stencil????
}
//TODO: is otherElement in this method required? is otherResidual correct?
evalPartialDerivativeCouplingBox_(element,
fvGeometry,
curElemVolVars,
elemFluxVarsCache,
elemBcTypes,
couplingMatrix,
residual);
}
/*!
* \brief Returns a reference to the problem.
......@@ -234,6 +319,7 @@ protected:
decltype(auto) otherProblem_(typename std::enable_if<!std::is_same<typename GET_PROP_TYPE(T, LowDimProblemTypeTag), SubProblemTypeTag>::value, void>::type* x = nullptr)
{ return globalProblem_().lowDimProblem(); }
// cell-centered
template<class JacobianMatrixCoupling>
void evalPartialDerivativeCoupling_(const Element& element,
const FVElementGeometry& fvGeometry,
......@@ -329,6 +415,102 @@ protected:
}
}
// box
template<class JacobianMatrixCoupling>
void evalPartialDerivativeCouplingBox_(const Element& element,
const FVElementGeometry& fvGeometry,
ElementVolumeVariables& curElemVolVars,
ElementFluxVariablesCache& elemFluxVarsCache,
const ElementBoundaryTypes& elemBcTypes,
JacobianMatrixCoupling& couplingMatrix,
SolutionVector& residual)
{
const auto& couplingStencil = globalProblem_().couplingManager().couplingStencil(element);
constexpr auto numEq = GET_PROP_VALUE(SubProblemTypeTag, NumEq);
for (auto globalJ : couplingStencil)
{
const auto otherResidual = globalProblem_().couplingManager().evalCouplingResidual(element,
fvGeometry,
curElemVolVars,
elemBcTypes,
elemFluxVarsCache);
auto& otherPriVars = otherProblem_().model().curSol()[globalJ];
auto originalOtherPriVars = otherPriVars;
// derivatives in the neighbors with repect to the current elements
ElementSolutionVector partialDeriv(fvGeometry.numScv());
for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
{
const Scalar eps = this->numericEpsilon(otherPriVars[pvIdx]);
Scalar delta = 0;
if (numericDifferenceMethod_ >= 0)
{
// we are not using backward differences, i.e. we need to
// calculate f(x + \epsilon)
// deflect primary variables
otherPriVars[pvIdx] += eps;
delta += eps;
// calculate the residual with the deflected primary variables
partialDeriv = globalProblem_().couplingManager().evalCouplingResidual(element,
fvGeometry,
curElemVolVars,
elemBcTypes,
elemFluxVarsCache);
}
else
{
// we are using backward differences, i.e. we don't need
// to calculate f(x + \epsilon) and we can recycle the
// (already calculated) residual f(x)
partialDeriv = otherResidual;
}
if (numericDifferenceMethod_ <= 0)
{
// we are not using forward differences, i.e. we
// need to calculate f(x - \epsilon)
// deflect the primary variables
otherPriVars[pvIdx] -= 2*eps;
delta += eps;
// calculate the residual with the deflected primary variables
partialDeriv -= globalProblem_().couplingManager().evalCouplingResidual(element,
fvGeometry,
curElemVolVars,
elemBcTypes,
elemFluxVarsCache);
}
else
{
// we are using forward differences, i.e. we don't need to
// calculate f(x - \epsilon) and we can recycle the
// (already calculated) residual f(x)
partialDeriv -= otherResidual;
}
// divide difference in residuals by the magnitude of the
// deflections between the two function evaluation
partialDeriv /= delta;
// restore the original state of the element solution vector
otherPriVars = originalOtherPriVars;
// update the global jacobian matrix (coupling block)
// this->updateGlobalJacobian_(couplingMatrix, scv.dofIndex(), globalJ, pvIdx, partialDeriv[scv.dofIndex()]);
for (auto&& scv : scvs(fvGeometry))
this->updateGlobalJacobian_(couplingMatrix, scv.dofIndex(), globalJ, pvIdx, partialDeriv[scv.index()]);
}
}
}
// The problem we would like to solve
GlobalProblem *globalProblemPtr_;
// The type of the numeric difference method (forward, center, backward)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment