Merge branch 'fix/subdomain-ghosts' into 'master'

Fix handling of ghost elements in subdomain local assemblers.

See merge request !2386

dumux/multidomain/subdomainboxlocalassembler.hh

@@ -118,21 +118,55 @@ public:
         this->asImp_().bindLocalViews();
         this->elemBcTypes().update(problem(), this->element(), this->fvGeometry());
 
-        // for the diagonal jacobian block
-        // forward to the internal implementation
-        const auto residual = this->asImp_().assembleJacobianAndResidualImpl(jacRow[domainId], *std::get<domainId>(gridVariables));
+        if (!this->elementIsGhost())
+        {
+            // for the diagonal jacobian block
+            // forward to the internal implementation
+            const auto residual = this->asImp_().assembleJacobianAndResidualImpl(jacRow[domainId], *std::get<domainId>(gridVariables));
 
-        // update the residual vector
-        for (const auto& scv : scvs(this->fvGeometry()))
-            res[scv.dofIndex()] += residual[scv.localDofIndex()];
+            // update the residual vector
+            for (const auto& scv : scvs(this->fvGeometry()))
+                res[scv.dofIndex()] += residual[scv.localDofIndex()];
 
-        // assemble the coupling blocks
-        using namespace Dune::Hybrid;
-        forEach(integralRange(Dune::Hybrid::size(jacRow)), [&](auto&& i)
+            // assemble the coupling blocks
+            using namespace Dune::Hybrid;
+            forEach(integralRange(Dune::Hybrid::size(jacRow)), [&](auto&& i)
+            {
+                if (i != id)
+                    this->assembleJacobianCoupling(i, jacRow, residual, gridVariables);
+            });
+        }
+        else
         {
-            if (i != id)
-                this->assembleJacobianCoupling(i, jacRow, residual, gridVariables);
-        });
+            using GridGeometry = typename GridVariables::GridGeometry;
+            using GridView = typename GridGeometry::GridView;
+            static constexpr auto dim = GridView::dimension;
+
+            int numVerticesLocal = this->element().subEntities(dim);
+
+            for (int i = 0; i < numVerticesLocal; ++i)
+            {
+                const auto vertex = this->element().template subEntity<dim>(i);
+
+                if (vertex.partitionType() == Dune::InteriorEntity ||
+                    vertex.partitionType() == Dune::BorderEntity)
+                {
+                    // do not change the non-ghost vertices
+                    continue;
+                }
+
+                // set main diagonal entries for the vertex
+                const auto vIdx = this->assembler().gridGeometry(domainId).vertexMapper().index(vertex);
+
+                typedef typename JacobianMatrix::block_type BlockType;
+                BlockType &J = jacRow[domainId][vIdx][vIdx];
+                for (int j = 0; j < BlockType::rows; ++j)
+                    J[j][j] = 1.0;
+
+                // set residual for the vertex
+                res[vIdx] = 0;
+            }
+        }
 
         // lambda for the incorporation of Dirichlet Bcs
         auto applyDirichlet = [&] (const auto& scvI,

dumux/multidomain/subdomaincclocalassembler.hh

@@ -337,6 +337,17 @@ public:
         Residuals origResiduals(numNeighbors + 1); origResiduals = 0.0;
         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
+        const auto evalNeighborFlux = [&] (const auto& neighbor, const auto& scvf)
+        {
+            if (neighbor.partitionType() == Dune::GhostEntity)
+                return LocalResidualValues(0.0);
+            else
+                return this->evalFluxResidual(neighbor, scvf);
+        };
+
         // get the elements in which we need to evaluate the fluxes
         // and calculate these in the undeflected state
         unsigned int j = 1;
@@ -344,7 +355,7 @@ public:
         {
             neighborElements[j-1] = gridGeometry.element(dataJ.globalJ);
             for (const auto scvfIdx : dataJ.scvfsJ)
-                origResiduals[j] += this->evalFluxResidual(neighborElements[j-1], fvGeometry.scvf(scvfIdx));
+                origResiduals[j] += evalNeighborFlux(neighborElements[j-1], fvGeometry.scvf(scvfIdx));
 
             ++j;
         }
@@ -368,12 +379,7 @@ public:
 
         for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
         {
-            // for ghost elements we assemble a 1.0 where the primary variable and zero everywhere else
-            // as we always solve for a delta of the solution with repect to the initial solution this
-            // results in a delta of zero for ghosts, we still need to do the neighbor derivatives though
-            // so we are not done yet here.
             partialDerivs = 0.0;
-            if (this->elementIsGhost()) partialDerivs[0][pvIdx] = 1.0;
 
             auto evalResiduals = [&](Scalar priVar)
             {
@@ -388,12 +394,12 @@ public:
                     gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
 
                 // calculate the residual with the deflected primary variables
-                if (!this->elementIsGhost()) partialDerivsTmp[0] = this->evalLocalResidual()[0];
+                partialDerivsTmp[0] = this->evalLocalResidual()[0];
 
                 // calculate the fluxes in the neighbors with the deflected primary variables
                 for (std::size_t k = 0; k < connectivityMap[globalI].size(); ++k)
                     for (auto scvfIdx : connectivityMap[globalI][k].scvfsJ)
-                        partialDerivsTmp[k+1] += this->evalFluxResidual(neighborElements[k], fvGeometry.scvf(scvfIdx));
+                        partialDerivsTmp[k+1] += evalNeighborFlux(neighborElements[k], fvGeometry.scvf(scvfIdx));
 
                 return partialDerivsTmp;
             };
@@ -404,6 +410,15 @@ public:
             NumericDifferentiation::partialDerivative(evalResiduals, elemSol[0][pvIdx], partialDerivs, origResiduals,
                                                       eps_(elemSol[0][pvIdx], pvIdx), numDiffMethod);
 
+            // Correct derivative for ghost elements, i.e. set a 1 for the derivative w.r.t. the
+            // current primary variable and a 0 elsewhere. As we always solve for a delta of the
+            // solution with repect to the initial one, this results in a delta of 0 for ghosts.
+            if (this->elementIsGhost())
+            {
+                partialDerivs[0] = 0.0;
+                partialDerivs[0][pvIdx] = 1.0;
+            }
+
             // add the current partial derivatives to the global jacobian matrix
             if constexpr (Problem::enableInternalDirichletConstraints())
             {
@@ -788,6 +803,17 @@ public:
     template<class JacobianMatrixDiagBlock, class GridVariables>
     LocalResidualValues assembleJacobianAndResidualImpl(JacobianMatrixDiagBlock& A, GridVariables& gridVariables)
     {
+        // treat ghost separately, we always want zero update for ghosts
+        if (this->elementIsGhost())
+        {
+            const auto globalI = this->assembler().gridGeometry(domainI).elementMapper().index(this->element());
+            for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
+                A[globalI][globalI][pvIdx][pvIdx] = 1.0;
+
+            // return zero residual
+            return LocalResidualValues(0.0);
+        }
+
         // get some aliases for convenience
         const auto& problem = this->problem();
         const auto& element = this->element();