Commit ede05f5f authored by Nils Kohl's avatar Nils Kohl 🌝
Browse files

Enabling HyTeG apply in mixed P1/P0 operators.

parent 8bbf638f
Pipeline #39557 failed with stages
in 32 minutes and 28 seconds
......@@ -165,8 +165,6 @@ class DGOperator : public Operator< DGFunction< real_t >, DGFunction< real_t > >
using indexing::Index;
using volumedofspace::indexing::ElementNeighborInfo;
WALBERLA_CHECK( updateType == Replace );
src.communicate( level );
const auto storage = this->getStorage();
......@@ -652,20 +650,51 @@ class DGOperator : public Operator< DGFunction< real_t >, DGFunction< real_t > >
{
if ( dim == 2 )
{
dstDofMemory[volumedofspace::indexing::index(
elementIdx.x(), elementIdx.y(), faceType, dstDofIdx, numDstDofs, level, dstMemLayout )] =
dstDofs( dstDofIdx );
if ( updateType == Replace )
{
dstDofMemory[volumedofspace::indexing::index(
elementIdx.x(), elementIdx.y(), faceType, dstDofIdx, numDstDofs, level, dstMemLayout )] =
dstDofs( dstDofIdx );
}
else if ( updateType == Add )
{
dstDofMemory[volumedofspace::indexing::index(
elementIdx.x(), elementIdx.y(), faceType, dstDofIdx, numDstDofs, level, dstMemLayout )] +=
dstDofs( dstDofIdx );
}
else
{
WALBERLA_ABORT( "Invalid update type." );
}
}
else
{
dstDofMemory[volumedofspace::indexing::index( elementIdx.x(),
elementIdx.y(),
elementIdx.z(),
cellType,
dstDofIdx,
numDstDofs,
level,
dstMemLayout )] = dstDofs( dstDofIdx );
if ( updateType == Replace )
{
dstDofMemory[volumedofspace::indexing::index( elementIdx.x(),
elementIdx.y(),
elementIdx.z(),
cellType,
dstDofIdx,
numDstDofs,
level,
dstMemLayout )] = dstDofs( dstDofIdx );
}
else if ( updateType == Add )
{
dstDofMemory[volumedofspace::indexing::index( elementIdx.x(),
elementIdx.y(),
elementIdx.z(),
cellType,
dstDofIdx,
numDstDofs,
level,
dstMemLayout )] += dstDofs( dstDofIdx );
}
else
{
WALBERLA_ABORT( "Invalid update type." );
}
}
}
}
......
......@@ -257,27 +257,10 @@ class P0ToP1Operator : public Operator< P0Function< real_t >, P1Function< real_t
auto vertexDoFIndicesArray = celldof::macrocell::getMicroVerticesFromMicroCell( elementIdx, cellType );
vertexDoFIndices.insert( vertexDoFIndices.begin(), vertexDoFIndicesArray.begin(), vertexDoFIndicesArray.end() );
}
#if 0
for ( uint_t dstDofIdx = 0; dstDofIdx < numDstDofs; dstDofIdx++ )
{
if ( dim == 2 )
{
srcDofs( srcDofIdx ) = srcDofMemory[vertexdof::macroface::index(
level, vertexDoFIndices[srcDofIdx].x(), vertexDoFIndices[srcDofIdx].y() )];
}
else
{
srcDofs( srcDofIdx ) = srcDofMemory[vertexdof::macrocell::index( level,
vertexDoFIndices[srcDofIdx].x(),
vertexDoFIndices[srcDofIdx].y(),
vertexDoFIndices[srcDofIdx].z() )];
}
}
#endif
if ( mat == nullptr )
{
// Matrix-vector multiplication.
WALBERLA_ABORT( "Not implemented / tested." );
dstDofs += localMat * srcDoF;
}
else
......@@ -295,7 +278,7 @@ class P0ToP1Operator : public Operator< P0Function< real_t >, P1Function< real_t
}
else
{
WALBERLA_ABORT( "Not implemented." );
WALBERLA_ABORT( "P0 to P1 assembly not implemented in 3D." );
}
}
}
......@@ -351,7 +334,7 @@ class P0ToP1Operator : public Operator< P0Function< real_t >, P1Function< real_t
}
else
{
WALBERLA_ABORT( "Not implemented." );
WALBERLA_ABORT( "P0 to P1 assembly not implemented in 3D." );
}
}
}
......@@ -407,7 +390,7 @@ class P0ToP1Operator : public Operator< P0Function< real_t >, P1Function< real_t
}
else
{
WALBERLA_ABORT( "Not implemented." );
WALBERLA_ABORT( "P0 to P1 assembly not implemented in 3D." );
}
}
}
......@@ -639,7 +622,7 @@ class P0ToP1Operator : public Operator< P0Function< real_t >, P1Function< real_t
}
else
{
WALBERLA_ABORT( "Not implemented." );
WALBERLA_ABORT( "P0 to P1 assembly not implemented in 3D." );
}
}
}
......
......@@ -113,8 +113,6 @@ class P1ToP0Operator : public Operator< P1Function< real_t >, P0Function< real_t
using indexing::Index;
using volumedofspace::indexing::ElementNeighborInfo;
WALBERLA_CHECK( updateType == Replace );
communication::syncFunctionBetweenPrimitives( src, level );
const auto storage = this->getStorage();
......@@ -640,20 +638,51 @@ class P1ToP0Operator : public Operator< P1Function< real_t >, P0Function< real_t
{
if ( dim == 2 )
{
dstDofMemory[volumedofspace::indexing::index(
elementIdx.x(), elementIdx.y(), faceType, dstDofIdx, numDstDofs, level, dstMemLayout )] =
dstDofs( dstDofIdx );
if ( updateType == Replace )
{
dstDofMemory[volumedofspace::indexing::index(
elementIdx.x(), elementIdx.y(), faceType, dstDofIdx, numDstDofs, level, dstMemLayout )] =
dstDofs( dstDofIdx );
}
else if ( updateType == Add )
{
dstDofMemory[volumedofspace::indexing::index(
elementIdx.x(), elementIdx.y(), faceType, dstDofIdx, numDstDofs, level, dstMemLayout )] +=
dstDofs( dstDofIdx );
}
else
{
WALBERLA_ABORT( "Invalid update type." );
}
}
else
{
dstDofMemory[volumedofspace::indexing::index( elementIdx.x(),
elementIdx.y(),
elementIdx.z(),
cellType,
dstDofIdx,
numDstDofs,
level,
dstMemLayout )] = dstDofs( dstDofIdx );
if ( updateType == Replace )
{
dstDofMemory[volumedofspace::indexing::index( elementIdx.x(),
elementIdx.y(),
elementIdx.z(),
cellType,
dstDofIdx,
numDstDofs,
level,
dstMemLayout )] = dstDofs( dstDofIdx );
}
else if ( updateType == Add )
{
dstDofMemory[volumedofspace::indexing::index( elementIdx.x(),
elementIdx.y(),
elementIdx.z(),
cellType,
dstDofIdx,
numDstDofs,
level,
dstMemLayout )] += dstDofs( dstDofIdx );
}
else
{
WALBERLA_ABORT( "Invalid update type." );
}
}
}
}
......
......@@ -32,6 +32,7 @@
#include "hyteg/petsc/PETScManager.hpp"
#include "hyteg/petsc/PETScSparseMatrix.hpp"
#include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp"
#include "hyteg/solvers/CGSolver.hpp"
using walberla::real_t;
using walberla::uint_t;
......@@ -112,14 +113,7 @@ real_t testHomogeneousDirichlet( const std::string& meshFile, const uint_t& leve
// Why do I need this?
f.interpolate( { u_x_expr, u_y_expr }, level, DirichletBoundary );
// simulate a multiply operation:
// M.apply( f, rhs, level, All, Replace );
auto f_vec = std::make_shared< PETScVector< real_t, P1DGEFunction > >();
f_vec->createVectorFromFunction( f, numerator, level, All );
auto rhs_vec = std::make_shared< PETScVector< real_t, P1DGEFunction > >();
rhs_vec->createVectorFromFunction( f, numerator, level, All );
Mpetsc.multiply( *f_vec, *rhs_vec );
rhs_vec->createFunctionFromVector( rhs, numerator, level, All );
M.apply( f, rhs, level, All, Replace );
rhs.interpolate( 0, level, DirichletBoundary );
u.getConformingPart()->interpolate( 0, level, DirichletBoundary );
......@@ -128,14 +122,9 @@ real_t testHomogeneousDirichlet( const std::string& meshFile, const uint_t& leve
solver.solve( L, u, rhs, level );
err.assign( { 1.0, -1.0 }, { u, sol }, level );
// calculate the error in the L2 norm
auto err_vec = std::make_shared< PETScVector< real_t, P1DGEFunction > >();
err_vec->createVectorFromFunction( err, numerator, level, All );
auto Merr_vec = std::make_shared< PETScVector< real_t, P1DGEFunction > >();
Merr_vec->createVectorFromFunction( err, numerator, level, All );
Mpetsc.multiply( *err_vec, *Merr_vec );
Merr_vec->createFunctionFromVector( Merr, numerator, level, All );
// calculate the error in the L2 norm
M.apply( err, Merr, level, All, Replace );
auto discrL2 = sqrt( err.dotGlobal( Merr, level, Inner ) );
if ( writeVTK )
......
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