Commit 35ca769e authored by Dennis Gläser's avatar Dennis Gläser
Browse files

Merge branch 'fix/cc-ghost-assembly' into 'master'

Fix/cc ghost assembly

Closes #817

See merge request !1862
parents a5cda132 1872459d
......@@ -189,13 +189,18 @@ public:
origResiduals[0] = this->evalLocalResidual()[0];
// lambda for convenient evaluation of the fluxes across scvfs in the neighbors
// if the neighbor is a ghost we don't want to add anything to their residual
// so we return 0 and omit computing the flux
auto evalNeighborFlux = [&] (const auto& neighbor, const auto& scvf)
{
return this->localResidual().evalFlux(this->problem(),
neighbor,
this->fvGeometry(),
this->curElemVolVars(),
this->elemFluxVarsCache(), scvf);
if (neighbor.partitionType() == Dune::GhostEntity)
return NumEqVector(0.0);
else
return this->localResidual().evalFlux(this->problem(),
neighbor,
this->fvGeometry(),
this->curElemVolVars(),
this->elemFluxVarsCache(), scvf);
};
// get the elements in which we need to evaluate the fluxes
......@@ -242,8 +247,8 @@ public:
if (enableGridFluxVarsCache)
gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
// calculate the residual with the deflected primary variables (except for ghosts)
if (!this->elementIsGhost()) partialDerivsTmp[0] = this->evalLocalResidual()[0];
// calculate the residual with the deflected primary variables
partialDerivsTmp[0] = this->evalLocalResidual()[0];
// calculate the fluxes in the neighbors with the deflected primary variables
for (std::size_t k = 0; k < numNeighbors; ++k)
......@@ -269,6 +274,7 @@ 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++)
{
// the diagonal entries
......@@ -430,6 +436,8 @@ class CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/tr
using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };
public:
using ParentType::ParentType;
......@@ -441,6 +449,17 @@ public:
*/
NumEqVector assembleJacobianAndResidualImpl(JacobianMatrix& A, const GridVariables& gridVariables)
{
// treat ghost separately, we always want zero update for ghosts
if (this->elementIsGhost())
{
const auto globalI = this->assembler().gridGeometry().elementMapper().index(this->element());
for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
A[globalI][globalI][pvIdx][pvIdx] = 1.0;
// return zero residual
return NumEqVector(0.0);
}
// assemble the undeflected residual
const auto residual = this->evalLocalResidual()[0];
......@@ -509,6 +528,8 @@ class CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/fa
using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };
public:
using ParentType::ParentType;
......@@ -520,6 +541,17 @@ public:
*/
NumEqVector assembleJacobianAndResidualImpl(JacobianMatrix& A, const GridVariables& gridVariables)
{
// treat ghost separately, we always want zero update for ghosts
if (this->elementIsGhost())
{
const auto globalI = this->assembler().gridGeometry().elementMapper().index(this->element());
for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
A[globalI][globalI][pvIdx][pvIdx] = 1.0;
// return zero residual
return NumEqVector(0.0);
}
// assemble the undeflected residual
const auto residual = this->evalLocalResidual()[0];
......
Supports Markdown
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