Skip to content

GitLab

Commit 9ca59c7f authored by Nils Kohl's avatar Nils Kohl 🌝
Browse files

Finished implementing vertex DoF ghost-layers required for facet integrals... 

Finished implementing vertex DoF ghost-layers required for facet integrals with coupling into neighboring macros. This feature is for example relevant for the P1CG x enriched Galerkin coupling (P1 -> P0).

Tested for vector Laplace problem on unit square in P1DGEConvergenceTests2D.cpp.
parent 9dbdf172
Pipeline #40631 failed with stages
in 55 minutes and 41 seconds
Hide whitespace changes
Inline Side-by-side
src/hyteg/dgfunctionspace/DGOperator.hpp
View file @ 9ca59c7f
......@@ -444,11 +444,12 @@ class DGOperator : public Operator< DGFunction< real_t >, DGFunction< real_t > >
// The neighboring micro-element coords have to be computed since they are now different as for an
// element on the same macro-volume.
std::vector< Index > neighborElementVertexIndices;
std::vector< Eigen::Matrix< real_t, 3, 1 > > neighborElementVertexCoords;
Eigen::Matrix< real_t, 3, 1 > neighborOppositeVertexCoords;
neighborInfo.macroBoundaryNeighborElementVertexCoords(
n, neighborElementVertexCoords, neighborOppositeVertexCoords );
n, neighborElementVertexIndices, neighborElementVertexCoords, neighborOppositeVertexCoords );
localMat.setZero();
form_->integrateFacetCoupling( dim,
......
src/hyteg/mixedoperators/P0ToP1Operator.hpp
View file @ 9ca59c7f
......@@ -158,7 +158,7 @@ class P0ToP1Operator : public Operator< P0Function< real_t >, P1Function< real_t
const auto face = storage->getFace( pid );
for ( const auto& [n, _] : face->getIndirectNeighborFaceIDsOverEdges() )
{
glMemory[n] = src.volumeDoFFunction()->glMemory( pid, level, n );
glMemory[n] = src.getDGFunction()->volumeDoFFunction()->glMemory( pid, level, n );
}
}
else
......@@ -167,7 +167,7 @@ class P0ToP1Operator : public Operator< P0Function< real_t >, P1Function< real_t
const auto cell = storage->getCell( pid );
for ( const auto& [n, _] : cell->getIndirectNeighborCellIDsOverFaces() )
{
glMemory[n] = src.volumeDoFFunction()->glMemory( pid, level, n );
glMemory[n] = src.getDGFunction()->volumeDoFFunction()->glMemory( pid, level, n );
}
}
......@@ -405,11 +405,12 @@ class P0ToP1Operator : public Operator< P0Function< real_t >, P1Function< real_t
// The neighboring micro-element coords have to be computed since they are now different as for an
// element on the same macro-volume.
std::vector< Index > neighborElementVertexIndices;
std::vector< Eigen::Matrix< real_t, 3, 1 > > neighborElementVertexCoords;
Eigen::Matrix< real_t, 3, 1 > neighborOppositeVertexCoords;
neighborInfo.macroBoundaryNeighborElementVertexCoords(
n, neighborElementVertexCoords, neighborOppositeVertexCoords );
n, neighborElementVertexIndices, neighborElementVertexCoords, neighborOppositeVertexCoords );
localMat.setZero();
form_->integrateFacetCoupling( dim,
......@@ -419,8 +420,8 @@ class P0ToP1Operator : public Operator< P0Function< real_t >, P1Function< real_t
neighborInfo.oppositeVertexCoords( n ),
neighborOppositeVertexCoords,
neighborInfo.outwardNormal( n ),
*src.basis(),
*dst.basis(),
*src.getDGFunction()->basis(),
dstBasis,
srcPolyDegree,
dstPolyDegree,
localMat );
......@@ -514,7 +515,7 @@ class P0ToP1Operator : public Operator< P0Function< real_t >, P1Function< real_t
const auto nSrcDof =
srcDofMemory[volumedofspace::indexing::index( neighborInfo.neighborElementIndices( n ).x(),
neighborInfo.neighborElementIndices( n ).y(),
faceType,
neighborInfo.neighborFaceType( n ),
0,
1,
level,
......
src/hyteg/mixedoperators/P1ToP0Operator.hpp
View file @ 9ca59c7f
......@@ -163,7 +163,7 @@ class P1ToP0Operator : public Operator< P1Function< real_t >, P0Function< real_t
const auto face = storage->getFace( pid );
for ( const auto& [n, _] : face->getIndirectNeighborFaceIDsOverEdges() )
{
glMemory[n] = face->template getData( src.getFaceGLDataID( n ) );
glMemory[n] = face->getData( src.getFaceGLDataID( n ) )->getPointer( level );
}
}
else
......@@ -172,7 +172,7 @@ class P1ToP0Operator : public Operator< P1Function< real_t >, P0Function< real_t
const auto cell = storage->getCell( pid );
for ( const auto& [n, _] : cell->getIndirectNeighborCellIDsOverFaces() )
{
glMemory[n] = cell->template getData( src.getCellGLDataID( n ) );
glMemory[n] = cell->getData( src.getCellGLDataID( n ) )->getPointer( level );
}
}
......@@ -440,20 +440,22 @@ class P1ToP0Operator : public Operator< P1Function< real_t >, P0Function< real_t
neighborInfo.oppositeVertexCoords( n ),
neighborOppositeVertexCoords,
neighborInfo.outwardNormal( n ),
*src.basis(),
*dst.basis(),
srcBasis,
*dst.getDGFunction()->basis(),
srcPolyDegree,
dstPolyDegree,
localMat );
// Now we need the DoFs from the neighboring element.
// --- START vertex DoF GL handling at macro-macro boundary -------------------------------------------
// Now we need the vertex DoFs (i.e. their indices) from the neighboring element.
//
// Those are partly on the current macro, and partly on the ghost-layer. This makes things a little
// more difficult for the CG element. We have to make sure that the order in which we plug in the
// element coords into the form is the same order we put the src vertex DoFs.
//
// What we have at hand are the local micro-vertex from the perspective of the neighboring macro in the
// correct order. Looping over those, we translate the indices to the indices of the micro-vertices
// What we have at hand are the local micro-vertices from the perspective of the neighboring macro in
// the correct order. Looping over those, we translate the indices to the indices of the micro-vertices
// local to the dst macro. One of those is _not_ located on the macro-macro-boundary. This is the
// micro-vertex that must be taken from the ghost-layer.
......@@ -462,55 +464,88 @@ class P1ToP0Operator : public Operator< P1Function< real_t >, P0Function< real_t
std::vector< uint_t > nSrcDoFArrIndices( numSrcDofs );
std::vector< bool > onGhostLayer( numSrcDofs );
for ( const auto& nElementVertexIdx : neighborElementVertexIndices )
if ( !storage->hasGlobalCells() )
{
// basis trafo to local volume macro
if ( vertexdof::macrocell::isOnCellFace( ... ) )
{
// get index on GL
}
else
// 2D
const auto face = storage->getFace( pid );
const auto edgePID = face->neighborEdges().at( neighborInfo.macroBoundaryID( n ) );
const auto neighborFace = storage->getFace(
face->getIndirectNeighborFaceIDsOverEdges().at( neighborInfo.macroBoundaryID( n ) ) );
const auto localEdgeIDNeighborFace = neighborFace->edge_index( edgePID );
for ( uint_t i = 0; i < neighborElementVertexIndices.size(); i++ )
{
const auto nElementVertexIdx = neighborElementVertexIndices[i];
// Check vertex DoF on interface or ghost-layer.
switch ( localEdgeIDNeighborFace )
{
case 0:
onGhostLayer[i] = nElementVertexIdx.y() != 0;
break;
case 1:
onGhostLayer[i] = nElementVertexIdx.x() != 0;
break;
case 2:
onGhostLayer[i] = nElementVertexIdx.x() + nElementVertexIdx.y() !=
levelinfo::num_microvertices_per_edge( level ) - 1;
break;
}
if ( !onGhostLayer[i] )
{
// If the DoF is not on the ghost-layer (i.e. it is on the interface) we need to obtain the
// logical index on the local macro volume. This is done via index "basis trafo".
// fetch DoFs from local mem
std::array< uint_t, 4 > srcBasis;
for ( uint_t ii = 0; ii < 3; ii++ )
{
if ( algorithms::contains( face->neighborVertices(),
neighborFace->neighborVertices().at( ii ) ) )
{
srcBasis[ii] = face->vertex_index( neighborFace->neighborVertices().at( ii ) );
}
else
{
srcBasis[ii] = face->vertex_index( face->get_vertex_opposite_to_edge(
face->neighborEdges().at( neighborInfo.macroBoundaryID( n ) ) ) );
}
}
srcBasis[3] = 3;
// Basis trafo to local macro.
const auto localIndex =
indexing::basisConversion( nElementVertexIdx,
srcBasis,
{ 0, 1, 2, 3 },
levelinfo::num_microvertices_per_edge( level ) );
nSrcDoFArrIndices[i] = vertexdof::macroface::index( level, localIndex.x(), localIndex.y() );
nSrcDofs[i] = srcDofMemory[nSrcDoFArrIndices[i]];
}
else
{
// Take DoF from GL memory.
nSrcDoFArrIndices[i] = volumedofspace::indexing::indexNeighborInGhostLayer(
neighborInfo.macroBoundaryID( n ),
elementIdx.x(),
elementIdx.y(),
faceType,
0,
1,
level,
volumedofspace::indexing::VolumeDoFMemoryLayout::AoS );
nSrcDofs[i] = glMemory[neighborInfo.macroBoundaryID( n )][nSrcDoFArrIndices[i]];
}
}
}
for ( uint_t srcDofIdx = 0; srcDofIdx < numSrcDofs; srcDofIdx++ )
else
{
if ( dim == 2 )
{
nSrcDoFArrIndices[srcDofIdx] =
volumedofspace::indexing::indexNeighborInGhostLayer( neighborInfo.macroBoundaryID( n ),
elementIdx.x(),
elementIdx.y(),
faceType,
srcDofIdx,
numSrcDofs,
level,
srcMemLayout );
}
else
{
nSrcDoFArrIndices[srcDofIdx] =
volumedofspace::indexing::indexNeighborInGhostLayer( neighborInfo.macroBoundaryID( n ),
elementIdx.x(),
elementIdx.y(),
elementIdx.z(),
cellType,
srcDofIdx,
numSrcDofs,
level,
srcMemLayout );
}
nSrcDofs( srcDofIdx ) = glMemory[neighborInfo.macroBoundaryID( n )][nSrcDoFArrIndices[srcDofIdx]];
// 3D
WALBERLA_ABORT( "Not implemented for 3D." );
}
// --- END vertex DoF GL handling at macro-macro boundary ---------------------------------------------
if ( mat == nullptr )
{
// Matrix-vector multiplication.
......@@ -519,38 +554,30 @@ class P1ToP0Operator : public Operator< P1Function< real_t >, P0Function< real_t
else
{
// Sparse assembly.
// TODO: maybe there is a nicer way to do the gl stuff ...
for ( uint_t dstDofIdx = 0; dstDofIdx < numDstDofs; dstDofIdx++ )
for ( uint_t srcDofIdx = 0; srcDofIdx < numSrcDofs; srcDofIdx++ )
{
for ( uint_t srcDofIdx = 0; srcDofIdx < numSrcDofs; srcDofIdx++ )
uint_t globalColIdx, globalRowIdx;
if ( dim == 2 )
{
uint_t globalRowIdx;
if ( dim == 2 )
{
globalRowIdx = dstDofMemory[volumedofspace::indexing::index( elementIdx.x(),
elementIdx.y(),
faceType,
dstDofIdx,
numDstDofs,
level,
dstMemLayout )];
}
else
{
globalRowIdx = dstDofMemory[volumedofspace::indexing::index( elementIdx.x(),
elementIdx.y(),
elementIdx.z(),
cellType,
dstDofIdx,
numDstDofs,
level,
dstMemLayout )];
}
const auto globalColIdx =
glMemory[neighborInfo.macroBoundaryID( n )][nSrcDoFArrIndices[srcDofIdx]];
globalRowIdx = dstDofMemory[volumedofspace::indexing::index(
elementIdx.x(), elementIdx.y(), faceType, 0, 1, level, dstMemLayout )];
}
else
{
globalRowIdx = dstDofMemory[volumedofspace::indexing::index(
elementIdx.x(), elementIdx.y(), elementIdx.z(), cellType, 0, 1, level, dstMemLayout )];
}
mat->addValue( globalRowIdx, globalColIdx, localMat( dstDofIdx, srcDofIdx ) );
if ( !onGhostLayer[srcDofIdx] )
{
globalColIdx = srcDofMemory[nSrcDoFArrIndices[srcDofIdx]];
}
else
{
globalColIdx = glMemory[neighborInfo.macroBoundaryID( n )][nSrcDoFArrIndices[srcDofIdx]];
}
mat->addValue( globalRowIdx, globalColIdx, localMat( 0, srcDofIdx ) );
}
}
}
......
src/hyteg/p1dgefunctionspace/P1DGEFunction.hpp
View file @ 9ca59c7f
......@@ -417,7 +417,7 @@ P1DGEFunction< ValueType >::P1DGEFunction( const std::string&
uint_t maxLevel,
BoundaryCondition bc )
: Function< P1DGEFunction< ValueType > >( name, storage, minLevel, maxLevel )
, u_conforming_{ std::make_shared< P1VectorFunction< ValueType > >( name + "_conforming", storage, minLevel, maxLevel, bc ) }
, u_conforming_{ std::make_shared< P1VectorFunction< ValueType > >( name + "_conforming", storage, minLevel, maxLevel, bc, true ) }
, u_discontinuous_{ std::make_shared< P0Function< ValueType > >( name + "_discontinuous", storage, minLevel, maxLevel, bc ) }
, basis_conforming_()
, basis_discontinuous_()
......
src/hyteg/p1functionspace/P1VectorFunction.hpp
View file @ 9ca59c7f
......@@ -53,15 +53,27 @@ class P1VectorFunction final : public CSFVectorFunction< P1VectorFunction< Value
size_t minLevel,
size_t maxLevel,
BoundaryCondition bc )
: P1VectorFunction( _name, storage, minLevel, maxLevel, bc, false )
{}
P1VectorFunction( const std::string& _name,
const std::shared_ptr< PrimitiveStorage >& storage,
size_t minLevel,
size_t maxLevel,
BoundaryCondition bc,
bool addVolumeGhostLayer )
: CSFVectorFunction< P1VectorFunction< ValueType > >( _name )
{
this->compFunc_.clear();
this->compFunc_.push_back( std::make_shared< VectorComponentType >( _name + "_u", storage, minLevel, maxLevel, bc ) );
this->compFunc_.push_back( std::make_shared< VectorComponentType >( _name + "_v", storage, minLevel, maxLevel, bc ) );
this->compFunc_.push_back(
std::make_shared< VectorComponentType >( _name + "_u", storage, minLevel, maxLevel, bc, addVolumeGhostLayer ) );
this->compFunc_.push_back(
std::make_shared< VectorComponentType >( _name + "_v", storage, minLevel, maxLevel, bc, addVolumeGhostLayer ) );
if ( storage->hasGlobalCells() )
{
this->compFunc_.push_back( std::make_shared< VectorComponentType >( _name + "_w", storage, minLevel, maxLevel, bc ) );
this->compFunc_.push_back(
std::make_shared< VectorComponentType >( _name + "_w", storage, minLevel, maxLevel, bc, addVolumeGhostLayer ) );
}
}
......
src/hyteg/p1functionspace/VertexDoFFunction.cpp
View file @ 9ca59c7f
......@@ -152,8 +152,14 @@ VertexDoFFunction< ValueType >::VertexDoFFunction( const std::string&
allocateMemory( level, *it.second );
}
communicators_[level]->addPackInfo( std::make_shared< VertexDoFPackInfo< ValueType > >(
level, vertexDataID_, edgeDataID_, faceDataID_, cellDataID_, this->getStorage() ) );
communicators_[level]->addPackInfo( std::make_shared< VertexDoFPackInfo< ValueType > >( level,
vertexDataID_,
edgeDataID_,
faceDataID_,
cellDataID_,
faceGhostLayerDataIDs_,
cellGhostLayerDataIDs_,
this->getStorage() ) );
additiveCommunicators_[level]->addPackInfo( std::make_shared< VertexDoFAdditivePackInfo< ValueType > >(
level, vertexDataID_, edgeDataID_, faceDataID_, cellDataID_, this->getStorage() ) );
}
......
src/hyteg/p1functionspace/VertexDoFFunction.hpp
View file @ 9ca59c7f
......@@ -107,13 +107,13 @@ class VertexDoFFunction final : public Function< VertexDoFFunction< ValueType >
{
WALBERLA_CHECK( hasVolumeGhostLayer() );
WALBERLA_CHECK( !this->getStorage()->hasGlobalCells() );
return faceGhostLayerDataIDs_[glID];
return faceGhostLayerDataIDs_.at( glID );
}
const PrimitiveDataID< FunctionMemory< ValueType >, Cell >& getCellGLDataID( uint_t glID ) const
{
WALBERLA_CHECK( hasVolumeGhostLayer() );
WALBERLA_CHECK( this->getStorage()->hasGlobalCells() );
return cellGhostLayerDataIDs_[glID];
return cellGhostLayerDataIDs_.at( glID );
}
void swap( const VertexDoFFunction< ValueType >& other, const uint_t& level, const DoFType& flag = All ) const;
......
src/hyteg/p1functionspace/VertexDoFPackInfo.hpp
View file @ 9ca59c7f
......@@ -888,7 +888,7 @@ void VertexDoFPackInfo< ValueType >::unpackFaceFromFace( Face*
template < typename ValueType >
void VertexDoFPackInfo< ValueType >::communicateLocalFaceToFace( const Face* sender, Face* receiver ) const
{
this->storage_->getTimingTree()->start( "VertexDoF - Face to Face" );
this->storage_.lock()->getTimingTree()->start( "VertexDoF - Face to Face" );
// Which local edges do we iterate on?
......@@ -946,7 +946,7 @@ void VertexDoFPackInfo< ValueType >::communicateLocalFaceToFace( const Face* sen
glMicroVolumeIdx++;
}
this->storage_->getTimingTree()->stop( "VertexDoF - Face to Face" );
this->storage_.lock()->getTimingTree()->stop( "VertexDoF - Face to Face" );
}
template < typename ValueType >
......
tests/hyteg/p1dgefunctionspace/P1DGEConvergenceTests2D.cpp
View file @ 9ca59c7f
......@@ -31,8 +31,8 @@
#include "hyteg/mesh/MeshInfo.hpp"
#include "hyteg/p1dgefunctionspace/P1DGEOperators.hpp"
#include "hyteg/petsc/PETScCGSolver.hpp"
#include "hyteg/petsc/PETScMinResSolver.hpp"
#include "hyteg/petsc/PETScManager.hpp"
#include "hyteg/petsc/PETScMinResSolver.hpp"
#include "hyteg/petsc/PETScSparseMatrix.hpp"
#include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp"
#include "hyteg/solvers/CGSolver.hpp"
......@@ -56,7 +56,7 @@ real_t testHomogeneousDirichlet( const std::string& meshFile, const uint_t& leve
P1DGEMassOperator M( storage, level, level );
numerator.enumerate( level );
#if 0
// solution as a lambda function
std::function< real_t( const Point3D& p ) > u_x_expr = []( const Point3D& p ) -> real_t {
const real_t x = p[0];
......@@ -99,7 +99,32 @@ real_t testHomogeneousDirichlet( const std::string& meshFile, const uint_t& leve
const real_t x8 = 8 * std::cos( x2 );
return -( x4 * x5 * x8 * std::cos( x6 ) - 16 * x4 * x7 * std::sin( x2 ) + x7 * x8 * std::cos( x3 ) );
};
#else
// solution as a lambda function
std::function< real_t( const Point3D& p ) > u_x_expr = []( const Point3D& p ) -> real_t {
const real_t x = p[0];
const real_t y = p[1];
return std::sin( M_PI * x ) * std::sin( M_PI * y );
};
std::function< real_t( const Point3D& p ) > u_y_expr = []( const Point3D& p ) -> real_t {
const real_t x = p[0];
const real_t y = p[1];
return std::sin( M_PI * x ) * std::sin( M_PI * y );
};
// rhs as a lambda function
std::function< real_t( const Point3D& p ) > f_x_expr = []( const Point3D& p ) -> real_t {
const real_t x = p[0];
const real_t y = p[1];
return 2 * M_PI * M_PI * std::sin( M_PI * x ) * std::sin( M_PI * y );
};
std::function< real_t( const Point3D& p ) > f_y_expr = []( const Point3D& p ) -> real_t {
const real_t x = p[0];
const real_t y = p[1];
return 2 * M_PI * M_PI * std::sin( M_PI * x ) * std::sin( M_PI * y );
};
#endif
P1DGEFunction< real_t > u( "u", storage, level, level );
P1DGEFunction< real_t > f( "f", storage, level, level );
P1DGEFunction< real_t > rhs( "rhs", storage, level, level );
......@@ -284,9 +309,9 @@ void runLaplace()
real_t currentRate = std::nan( "" );
for ( uint_t level = 3; level <= 6; level++ )
{
lastError = currentError;
lastError = currentError;
// currentError = hyteg::testHomogeneousDirichlet( "../../data/meshes/tri_1el.msh", level, false );
currentError = hyteg::testHomogeneousDirichlet( "../../data/meshes/tri_2el.msh", level, true );
currentError = hyteg::testHomogeneousDirichlet( "../../data/meshes/quad_16el.msh", level, true );
currentRate = lastError / currentError;
WALBERLA_LOG_INFO_ON_ROOT( walberla::format( "%6d|%15.2e|%15.2e", level, currentError, currentRate ) );
}
......
Markdown is supported
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