[box] Make periodic boundaries work in sequential

dumux/assembly/boxlocalassembler.hh

@@ -107,6 +107,16 @@ public:
                 if (scvI.localDofIndex() == scvJ.localDofIndex())
                     jac[scvI.dofIndex()][scvI.dofIndex()][eqIdx][pvIdx] = 1.0;
             }
+
+            // if a periodic dof has dirichlet values also apply the same dirichlet values to the other dof
+            if (this->assembler().fvGridGeometry().dofOnPeriodicBoundary(scvI.dofIndex()))
+            {
+                const auto periodicDof = this->assembler().fvGridGeometry().periodicallyMappedDof(scvI.dofIndex());
+                res[periodicDof][eqIdx] = this->curElemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx];
+                const auto end = jac[periodicDof].end();
+                for (auto it = jac[periodicDof].begin(); it != end; ++it)
+                    (*it) = periodicDof != it.index() ? 0.0 : 1.0;
+            }
         };
 
         this->asImp_().evalDirichletBoundaries(applyDirichlet);

dumux/assembly/fvassembler.hh

@@ -27,6 +27,7 @@
 #include <type_traits>
 
 #include <dune/istl/matrixindexset.hh>
+#include <dune/istl/io.hh>
 
 #include <dumux/common/properties.hh>
 #include <dumux/common/timeloop.hh>
@@ -121,6 +122,8 @@ public:
             LocalAssembler localAssembler(*this, element, curSol);
             localAssembler.assembleJacobianAndResidual(*jacobian_, *residual_, *gridVariables_, partialReassembler);
         });
+
+        enforcePeriodicConstraints_(*jacobian_, *residual_, *fvGridGeometry_);
     }
 
     /*!
@@ -395,6 +398,29 @@ private:
             DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system");
     }
 
+    template<class GG> std::enable_if_t<GG::discMethod == DiscretizationMethod::box, void>
+    enforcePeriodicConstraints_(JacobianMatrix& jac, SolutionVector& res, const GG& fvGridGeometry)
+    {
+        for (const auto& m : fvGridGeometry.periodicVertexMap())
+        {
+            if (m.first < m.second)
+            {
+                // add the second row to the first
+                res[m.first] += res[m.second];
+                const auto end = jac[m.second].end();
+                for (auto it = jac[m.second].begin(); it != end; ++it)
+                    jac[m.first][it.index()] += (*it);
+
+                // enforce constraint in second row
+                for (auto it = jac[m.second].begin(); it != end; ++it)
+                    (*it) = it.index() == m.second ? 1.0 : it.index() == m.first ? -1.0 : 0.0;
+            }
+        }
+    }
+
+    template<class GG> std::enable_if_t<GG::discMethod != DiscretizationMethod::box, void>
+    enforcePeriodicConstraints_(JacobianMatrix& jac, SolutionVector& res, const GG& fvGridGeometry) {}
+
     //! pointer to the problem to be solved
     std::shared_ptr<const Problem> problem_;
 

dumux/assembly/jacobianpattern.hh

@@ -51,11 +51,20 @@ Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
             for (unsigned int vIdx = 0; vIdx < element.subEntities(dim); ++vIdx)
             {
                 const auto globalI = gridGeometry.vertexMapper().subIndex(element, vIdx, dim);
-                for (unsigned int vIdx2 = vIdx; vIdx2 < element.subEntities(dim); ++vIdx2)
+                for (unsigned int vIdx2 = 0; vIdx2 < element.subEntities(dim); ++vIdx2)
                 {
                     const auto globalJ = gridGeometry.vertexMapper().subIndex(element, vIdx2, dim);
                     pattern.add(globalI, globalJ);
                     pattern.add(globalJ, globalI);
+
+                    if (gridGeometry.dofOnPeriodicBoundary(globalI) && globalI != globalJ)
+                    {
+                        const auto globalIP = gridGeometry.periodicallyMappedDof(globalI);
+                        pattern.add(globalIP, globalI);
+                        pattern.add(globalI, globalIP);
+                        if (globalI > globalIP)
+                            pattern.add(globalIP, globalJ);
+                    }
                 }
             }
         }

dumux/discretization/box/fvelementgeometry.hh

@@ -313,7 +313,7 @@ private:
         // construct the sub control volume faces on the domain boundary
         for (const auto& intersection : intersections(fvGridGeometry().gridView(), element))
         {
-            if (intersection.boundary())
+            if (intersection.boundary() && !intersection.neighbor())
             {
                 const auto isGeometry = intersection.geometry();
 

dumux/discretization/box/fvgridgeometry.hh

@@ -26,6 +26,8 @@
 #ifndef DUMUX_DISCRETIZATION_BOX_GRID_FVGEOMETRY_HH
 #define DUMUX_DISCRETIZATION_BOX_GRID_FVGEOMETRY_HH
 
+#include <unordered_map>
+
 #include <dune/geometry/referenceelements.hh>
 #include <dune/localfunctions/lagrange/pqkfactory.hh>
 
@@ -116,9 +118,9 @@ public:
     : ParentType(gridView)
     {
         // Check if the overlap size is what we expect
-        if (!CheckOverlapSize<DiscretizationMethod::box>::isValid(gridView))
-            DUNE_THROW(Dune::InvalidStateException, "The box discretization method only works with zero overlap for parallel computations. "
-                                                     << " Set the parameter \"Grid.Overlap\" in the input file.");
+        // if (!CheckOverlapSize<DiscretizationMethod::box>::isValid(gridView))
+        //     DUNE_THROW(Dune::InvalidStateException, "The box discretization method only works with zero overlap for parallel computations. "
+        //                                              << " Set the parameter \"Grid.Overlap\" in the input file.");
     }
 
     //! the vertex mapper is the dofMapper
@@ -141,7 +143,7 @@ public:
 
     //! The total number of degrees of freedom
     std::size_t numDofs() const
-    { return this->gridView().size(dim); }
+    { return this->vertexMapper().size(); }
 
     //! update all fvElementGeometries (do this again after grid adaption)
     void update()
@@ -210,7 +212,7 @@ public:
             // construct the sub control volume faces on the domain boundary
             for (const auto& intersection : intersections(this->gridView(), element))
             {
-                if (intersection.boundary())
+                if (intersection.boundary() && !intersection.neighbor())
                 {
                     const auto isGeometry = intersection.geometry();
                     // count
@@ -246,6 +248,43 @@ public:
                         boundaryDofIndices_[vIdxGlobal] = true;
                     }
                 }
+
+                // inform the grid geometry if we have periodic boundaries
+                else if (intersection.boundary() && intersection.neighbor())
+                {
+                    this->setPeriodic();
+
+                    // find the mapped periodic vertex of all vertices on periodic boundaries
+                    const auto fIdx = intersection.indexInInside();
+                    const auto numFaceVerts = referenceElement.size(fIdx, 1, dim);
+                    const auto eps = 1e-7*(elementGeometry.corner(1) - elementGeometry.corner(0)).two_norm();
+                    for (int localVIdx = 0; localVIdx < numFaceVerts; ++localVIdx)
+                    {
+                        const auto vIdx = referenceElement.subEntity(fIdx, 1, localVIdx, dim);
+                        const auto vIdxGlobal = this->vertexMapper().subIndex(element, vIdx, dim);
+                        const auto vPos = elementGeometry.corner(vIdx);
+
+                        const auto& outside = intersection.outside();
+                        const auto outsideGeometry = outside.geometry();
+                        for (const auto& isOutside : intersections(this->gridView(), outside))
+                        {
+                            // only check periodic vertices of the periodic neighbor
+                            if (isOutside.boundary() && isOutside.neighbor())
+                            {
+                                const auto fIdxOutside = isOutside.indexInInside();
+                                const auto numFaceVertsOutside = referenceElement.size(fIdxOutside, 1, dim);
+                                for (int localVIdxOutside = 0; localVIdxOutside < numFaceVertsOutside; ++localVIdxOutside)
+                                {
+                                    const auto vIdxOutside = referenceElement.subEntity(fIdxOutside, 1, localVIdxOutside, dim);
+                                    const auto vPosOutside = outsideGeometry.corner(vIdxOutside);
+                                    const auto shift = std::abs((this->bBoxMax()-this->bBoxMin())*intersection.centerUnitOuterNormal());
+                                    if (std::abs((vPosOutside-vPos).two_norm() - shift) < eps)
+                                        periodicVertexMap_[vIdxGlobal] = this->vertexMapper().subIndex(outside, vIdxOutside, dim);
+                                }
+                            }
+                        }
+                    }
+                }
             }
         }
     }
@@ -266,6 +305,14 @@ public:
     bool dofOnBoundary(unsigned int dofIdx) const
     { return boundaryDofIndices_[dofIdx]; }
 
+    //! If a vertex / d.o.f. is on a periodic boundary
+    bool dofOnPeriodicBoundary(std::size_t dofIdx) const
+    { return periodicVertexMap_.count(dofIdx); }
+
+    //! The index of the vertex / d.o.f. on the other side of the periodic boundary
+    std::size_t periodicallyMappedDof(std::size_t dofIdx) const
+    { return periodicVertexMap_.at(dofIdx); }
+
 private:
 
     const FeCache feCache_;
@@ -279,6 +326,9 @@ private:
 
     // vertices on the boudary
     std::vector<bool> boundaryDofIndices_;
+
+    // a map for periodic boundary vertices
+    std::unordered_map<std::size_t, std::size_t> periodicVertexMap_;
 };
 
 /*!
@@ -326,9 +376,9 @@ public:
     : ParentType(gridView)
     {
         // Check if the overlap size is what we expect
-        if (!CheckOverlapSize<DiscretizationMethod::box>::isValid(gridView))
-            DUNE_THROW(Dune::InvalidStateException, "The box discretization method only works with zero overlap for parallel computations. "
-                                                     << " Set the parameter \"Grid.Overlap\" in the input file.");
+        // if (!CheckOverlapSize<DiscretizationMethod::box>::isValid(gridView))
+        //     DUNE_THROW(Dune::InvalidStateException, "The box discretization method only works with zero overlap for parallel computations. "
+        //                                              << " Set the parameter \"Grid.Overlap\" in the input file.");
     }
 
     //! the vertex mapper is the dofMapper
@@ -351,7 +401,7 @@ public:
 
     //! The total number of degrees of freedom
     std::size_t numDofs() const
-    { return this->gridView().size(dim); }
+    { return this->vertexMapper().size(); }
 
     //! update all fvElementGeometries (do this again after grid adaption)
     void update()
@@ -376,7 +426,7 @@ public:
             // store the sub control volume face indices on the domain boundary
             for (const auto& intersection : intersections(this->gridView(), element))
             {
-                if (intersection.boundary())
+                if (intersection.boundary() && !intersection.neighbor())
                 {
                     const auto isGeometry = intersection.geometry();
                     numScvf_ += isGeometry.corners();
@@ -393,8 +443,49 @@ public:
                         boundaryDofIndices_[vIdxGlobal] = true;
                     }
                 }
+
+                // inform the grid geometry if we have periodic boundaries
+                else if (intersection.boundary() && intersection.neighbor())
+                {
+                    this->setPeriodic();
+
+                    // find the mapped periodic vertex of all vertices on periodic boundaries
+                    const auto fIdx = intersection.indexInInside();
+                    const auto numFaceVerts = referenceElement.size(fIdx, 1, dim);
+                    const auto eps = 1e-7*(elementGeometry.corner(1) - elementGeometry.corner(0)).two_norm();
+                    for (int localVIdx = 0; localVIdx < numFaceVerts; ++localVIdx)
+                    {
+                        const auto vIdx = referenceElement.subEntity(fIdx, 1, localVIdx, dim);
+                        const auto vIdxGlobal = this->vertexMapper().subIndex(element, vIdx, dim);
+                        const auto vPos = elementGeometry.corner(vIdx);
+
+                        const auto& outside = intersection.outside();
+                        const auto outsideGeometry = outside.geometry();
+                        for (const auto& isOutside : intersections(this->gridView(), outside))
+                        {
+                            // only check periodic vertices of the periodic neighbor
+                            if (isOutside.boundary() && isOutside.neighbor())
+                            {
+                                const auto fIdxOutside = isOutside.indexInInside();
+                                const auto numFaceVertsOutside = referenceElement.size(fIdxOutside, 1, dim);
+                                for (int localVIdxOutside = 0; localVIdxOutside < numFaceVertsOutside; ++localVIdxOutside)
+                                {
+                                    const auto vIdxOutside = referenceElement.subEntity(fIdxOutside, 1, localVIdxOutside, dim);
+                                    const auto vPosOutside = outsideGeometry.corner(vIdxOutside);
+                                    const auto shift = std::abs((this->bBoxMax()-this->bBoxMin())*intersection.centerUnitOuterNormal());
+                                    if (std::abs((vPosOutside-vPos).two_norm() - shift) < eps)
+                                        periodicVertexMap_[vIdxGlobal] = this->vertexMapper().subIndex(outside, vIdxOutside, dim);
+                                }
+                            }
+                        }
+                    }
+                }
             }
         }
+
+
+                for (const auto& pair : periodicVertexMap_)
+                    std::cout << pair.first << "  " << pair.second << std::endl;
     }
 
     //! The finite element cache for creating local FE bases
@@ -402,9 +493,21 @@ public:
     { return feCache_; }
 
     //! If a vertex / d.o.f. is on the boundary
-    bool dofOnBoundary(unsigned int dofIdx) const
+    bool dofOnBoundary(std::size_t dofIdx) const
     { return boundaryDofIndices_[dofIdx]; }
 
+    //! If a vertex / d.o.f. is on a periodic boundary
+    bool dofOnPeriodicBoundary(std::size_t dofIdx) const
+    { return periodicVertexMap_.count(dofIdx); }
+
+    //! The index of the vertex / d.o.f. on the other side of the periodic boundary
+    std::size_t periodicallyMappedDof(std::size_t dofIdx) const
+    { return periodicVertexMap_.at(dofIdx); }
+
+    //! The index of the vertex / d.o.f. on the other side of the periodic boundary
+    const std::unordered_map<std::size_t, std::size_t> periodicVertexMap() const
+    { return periodicVertexMap_; }
+
 private:
 
     const FeCache feCache_;
@@ -417,6 +520,9 @@ private:
 
     // vertices on the boudary
     std::vector<bool> boundaryDofIndices_;
+
+    // a map for periodic boundary vertices
+    std::unordered_map<std::size_t, std::size_t> periodicVertexMap_;
 };
 
 } // end namespace Dumux

test/porousmediumflow/1p/implicit/periodicbc/CMakeLists.txt

@@ -38,12 +38,12 @@ dune_add_test(NAME test_1p_periodic_tpfa_caching
                                 ${CMAKE_CURRENT_BINARY_DIR}/test_1p_tpfa_caching-00001.vtu
                         --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_periodic_tpfa test_1p.input -Problem.Name test_1p_tpfa_caching")
 
-#dune_add_test(NAME test_1p_periodic_box
-#              SOURCES test_1pfv.cc
-#              COMPILE_DEFINITIONS TYPETAG=OnePIncompressibleBox
-#              CMAKE_GUARD dune-spgrid_FOUND
-#              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
-#              CMD_ARGS  --script fuzzy
-#                        --files ${CMAKE_SOURCE_DIR}/test/references/1ptestcc-reference.vtu
-#                                ${CMAKE_CURRENT_BINARY_DIR}/test_1p_box-00001.vtu
-#                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_periodic_box test_1p.input -Problem.Name test_1p_box")
+dune_add_test(NAME test_1p_periodic_box
+              SOURCES test_1pfv.cc
+              COMPILE_DEFINITIONS TYPETAG=OnePIncompressibleBox
+              CMAKE_GUARD dune-spgrid_FOUND
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMD_ARGS  --script fuzzy
+                        --files ${CMAKE_SOURCE_DIR}/test/references/1ptestcc-reference.vtu
+                                ${CMAKE_CURRENT_BINARY_DIR}/test_1p_box-00001.vtu
+                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_periodic_box test_1p.input -Problem.Name test_1p_box")

test/porousmediumflow/1p/implicit/periodicbc/periodic.dgf

@@ -3,7 +3,7 @@ DGF
 INTERVAL
 0 0            % lower left corner
 1 1            % upper right corner
-100 100        % number of cells in each direction
+100 100            % number of cells in each direction
 #
 
 PERIODICFACETRANSFORMATION