Commit 32c633e8 authored by Timo Koch's avatar Timo Koch
Browse files

[box] Make periodic boundaries work in sequential

parent 10718be6
......@@ -137,6 +137,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);
......@@ -206,7 +216,7 @@ public:
{
const auto dirichletValues = this->problem().dirichlet(this->element(), scvI);
// set the dirichlet conditions in residual and jacobian
// set the Dirichlet conditions in residual and jacobian
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
{
if (bcTypes.isDirichlet(eqIdx))
......
......@@ -121,6 +121,8 @@ public:
LocalAssembler localAssembler(*this, element, curSol);
localAssembler.assembleJacobianAndResidual(*jacobian_, *residual_, *gridVariables_, partialReassembler);
});
enforcePeriodicConstraints_(*jacobian_, *residual_, *fvGridGeometry_);
}
/*!
......@@ -395,6 +397,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_;
......
......@@ -51,11 +51,19 @@ 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);
}
}
}
}
......
......@@ -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();
......
......@@ -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>
......@@ -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,6 +443,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);
}
}
}
}
}
}
}
}
......@@ -402,9 +489,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 +516,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
......
......@@ -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/test-1p-box-periodic-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")
......@@ -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
......
This diff is collapsed.
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