Commit d9b8ee59 authored by Marcus Mohr's avatar Marcus Mohr
Browse files

Update P[12]ElementwiseOperator

Changing PetscInt to matIdx_t allows to remove the preprocessor checks
for HYTEG_BUILD_WITH_PETSC and should fix issue with pipeline on the
new toMatrix() method.

Commit also adds P1ElementwiseOperator::toMatrix() and makes Petsc-
MatrixAssembly test use the new API for testing P[12]ElementwiseOperator.
parent bd1b3a13
......@@ -573,13 +573,11 @@ void P1ElementwiseOperator< P1Form >::computeLocalDiagonalContributions3D( const
}
}
#ifdef HYTEG_BUILD_WITH_PETSC
// Assemble operator as sparse matrix
template < class P1Form >
void P1ElementwiseOperator< P1Form >::assembleLocalMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P1Function< PetscInt >& src,
const P1Function< PetscInt >& dst,
const P1Function< matIdx_t >& src,
const P1Function< matIdx_t >& dst,
uint_t level,
DoFType flag ) const
{
......@@ -599,11 +597,11 @@ void P1ElementwiseOperator< P1Form >::assembleLocalMatrix( const std::shared_ptr
Cell& cell = *macroIter.second;
// get hold of the actual numerical data in the two indexing functions
PrimitiveDataID< FunctionMemory< PetscInt >, Cell > dstVertexDoFIdx = dst.getCellDataID();
PrimitiveDataID< FunctionMemory< PetscInt >, Cell > srcVertexDoFIdx = src.getCellDataID();
PrimitiveDataID< FunctionMemory< matIdx_t >, Cell > dstVertexDoFIdx = dst.getCellDataID();
PrimitiveDataID< FunctionMemory< matIdx_t >, Cell > srcVertexDoFIdx = src.getCellDataID();
PetscInt* srcIdx = cell.getData( srcVertexDoFIdx )->getPointer( level );
PetscInt* dstIdx = cell.getData( dstVertexDoFIdx )->getPointer( level );
matIdx_t* srcIdx = cell.getData( srcVertexDoFIdx )->getPointer( level );
matIdx_t* dstIdx = cell.getData( dstVertexDoFIdx )->getPointer( level );
// loop over micro-cells
for ( const auto& cType : celldof::allCellTypes )
......@@ -635,11 +633,11 @@ void P1ElementwiseOperator< P1Form >::assembleLocalMatrix( const std::shared_ptr
indexing::IndexIncrement offset;
// get hold of the actual numerical data in the two functions
PrimitiveDataID< FunctionMemory< PetscInt >, Face > dstVertexDoFIdx = dst.getFaceDataID();
PrimitiveDataID< FunctionMemory< PetscInt >, Face > srcVertexDoFIdx = src.getFaceDataID();
PrimitiveDataID< FunctionMemory< matIdx_t >, Face > dstVertexDoFIdx = dst.getFaceDataID();
PrimitiveDataID< FunctionMemory< matIdx_t >, Face > srcVertexDoFIdx = src.getFaceDataID();
PetscInt* srcIndices = face.getData( srcVertexDoFIdx )->getPointer( level );
PetscInt* dstIndices = face.getData( dstVertexDoFIdx )->getPointer( level );
matIdx_t* srcIndices = face.getData( srcVertexDoFIdx )->getPointer( level );
matIdx_t* dstIndices = face.getData( dstVertexDoFIdx )->getPointer( level );
// the explicit uint_c cast prevents a segfault in intel compiler 2018.4
// now loop over micro-faces of macro-face
......@@ -671,8 +669,8 @@ void P1ElementwiseOperator< P1Form >::localMatrixAssembly2D( const std::shared_p
const uint_t xIdx,
const uint_t yIdx,
const P1Elements::P1Elements2D::P1Element& element,
const PetscInt* const srcIdx,
const PetscInt* const dstIdx ) const
const matIdx_t* const srcIdx,
const matIdx_t* const dstIdx ) const
{
Matrix3r elMat;
indexing::Index nodeIdx;
......@@ -725,8 +723,8 @@ void P1ElementwiseOperator< P1Form >::localMatrixAssembly3D( const std::shared_p
const uint_t level,
const indexing::Index& microCell,
const celldof::CellType cType,
const PetscInt* const srcIdx,
const PetscInt* const dstIdx ) const
const matIdx_t* const srcIdx,
const matIdx_t* const dstIdx ) const
{
// determine coordinates of vertices of micro-element
std::array< indexing::Index, 4 > verts = celldof::macrocell::getMicroVerticesFromMicroCell( microCell, cType );
......@@ -765,8 +763,6 @@ void P1ElementwiseOperator< P1Form >::localMatrixAssembly3D( const std::shared_p
mat->addValues( rowIdx, colIdx, blockMatData );
}
#endif
// P1ElementwiseLaplaceOperator
template class P1ElementwiseOperator<
P1FenicsForm< p1_diffusion_cell_integral_0_otherwise, p1_tet_diffusion_cell_integral_0_otherwise > >;
......
......@@ -72,7 +72,6 @@ P1ElementwiseOperator( const std::shared_ptr< PrimitiveStorage >& storage, size_
size_t level,
DoFType flag ) const override;
#ifdef HYTEG_BUILD_WITH_PETSC
/// Assemble operator as sparse matrix
///
/// \param mat a sparse matrix proxy
......@@ -83,11 +82,29 @@ P1ElementwiseOperator( const std::shared_ptr< PrimitiveStorage >& storage, size_
///
/// \note src and dst are legal to and often will be the same function object
void assembleLocalMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P1Function< PetscInt >& src,
const P1Function< PetscInt >& dst,
const P1Function< matIdx_t >& src,
const P1Function< matIdx_t >& dst,
uint_t level,
DoFType flag ) const;
#endif
/// Assemble operator as sparse matrix.
///
/// \param mat a sparse matrix proxy
/// \param src P2Function for determining column indices
/// \param dst P2Function for determining row indices
/// \param level level in mesh hierarchy for which local operator is to be assembled
/// \param flag ignored
///
/// \note src and dst are legal to and often will be the same function object
/// \note we should rename assembleLocalMatrix() to toMatrix() and skip the delegation
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P1Function< matIdx_t >& src,
const P1Function< matIdx_t >& dst,
uint_t level,
DoFType flag ) const override
{
this->assembleLocalMatrix( mat, src, dst, level, flag );
}
/// \brief Pre-computes the local stiffness matrices for each (micro-)element and stores them all in memory.
///
......@@ -182,24 +199,22 @@ P1ElementwiseOperator( const std::shared_ptr< PrimitiveStorage >& storage, size_
const celldof::CellType cType,
real_t* const vertexData );
#ifdef HYTEG_BUILD_WITH_PETSC
void localMatrixAssembly2D( const std::shared_ptr< SparseMatrixProxy >& mat,
const Face& face,
const uint_t level,
const uint_t xIdx,
const uint_t yIdx,
const P1Elements::P1Elements2D::P1Element& element,
const PetscInt* const srcIdx,
const PetscInt* const dstIdx ) const;
const matIdx_t* const srcIdx,
const matIdx_t* const dstIdx ) const;
void localMatrixAssembly3D( const std::shared_ptr< SparseMatrixProxy >& mat,
const Cell& cell,
const uint_t level,
const indexing::Index& microCell,
const celldof::CellType cType,
const PetscInt* const srcIdx,
const PetscInt* const dstIdx ) const;
#endif
const matIdx_t* const srcIdx,
const matIdx_t* const dstIdx ) const;
/// Trigger (re)computation of diagonal matrix entries (central operator weights)
/// Allocates the required memory if the function was not yet allocated.
......
......@@ -320,11 +320,11 @@ void P2ElementwiseOperator< P2Form >::smooth_jac( const P2Function< real_t >& ds
// compute the current residual
this->apply( src, dst, level, flag );
dst.assign( { real_c( 1 ), real_c( -1 ) }, { rhs, dst }, level, flag );
dst.assign( {real_c( 1 ), real_c( -1 )}, {rhs, dst}, level, flag );
// perform Jacobi update step
dst.multElementwise( { *getInverseDiagonalValues(), dst }, level, flag );
dst.assign( { 1.0, omega }, { src, dst }, level, flag );
dst.multElementwise( {*getInverseDiagonalValues(), dst}, level, flag );
dst.assign( {1.0, omega}, {src, dst}, level, flag );
this->stopTiming( "smooth_jac" );
}
......@@ -596,7 +596,7 @@ void P2ElementwiseOperator< P2Form >::computeLocalDiagonalContributions2D( const
// assemble local element matrix
form.setGeometryMap( face.getGeometryMap() );
form.integrateAll( { v0, v1, v2 }, elMat );
form.integrateAll( {v0, v1, v2}, elMat );
// get global indices for local dofs
dofDataIdx[0] = vertexdof::macroface::indexFromVertex( level, xIdx, yIdx, element[0] );
......@@ -662,13 +662,11 @@ P2Form P2ElementwiseOperator< P2Form >::getForm() const
return form_;
}
#ifdef HYTEG_BUILD_WITH_PETSC
// Assemble operator as sparse matrix
template < class P2Form >
void P2ElementwiseOperator< P2Form >::assembleLocalMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P2Function< PetscInt >& src,
const P2Function< PetscInt >& dst,
const P2Function< matIdx_t >& src,
const P2Function< matIdx_t >& dst,
uint_t level,
DoFType flag ) const
{
......@@ -688,17 +686,17 @@ void P2ElementwiseOperator< P2Form >::assembleLocalMatrix( const std::shared_ptr
Cell& cell = *macroIter.second;
// get hold of the actual numerical data in the two indexing functions
PrimitiveDataID< FunctionMemory< PetscInt >, Cell > dstVertexDoFIdx = dst.getVertexDoFFunction().getCellDataID();
PrimitiveDataID< FunctionMemory< PetscInt >, Cell > srcVertexDoFIdx = src.getVertexDoFFunction().getCellDataID();
PrimitiveDataID< FunctionMemory< matIdx_t >, Cell > dstVertexDoFIdx = dst.getVertexDoFFunction().getCellDataID();
PrimitiveDataID< FunctionMemory< matIdx_t >, Cell > srcVertexDoFIdx = src.getVertexDoFFunction().getCellDataID();
PrimitiveDataID< FunctionMemory< PetscInt >, Cell > dstEdgeDoFIdx = dst.getEdgeDoFFunction().getCellDataID();
PrimitiveDataID< FunctionMemory< PetscInt >, Cell > srcEdgeDoFIdx = src.getEdgeDoFFunction().getCellDataID();
PrimitiveDataID< FunctionMemory< matIdx_t >, Cell > dstEdgeDoFIdx = dst.getEdgeDoFFunction().getCellDataID();
PrimitiveDataID< FunctionMemory< matIdx_t >, Cell > srcEdgeDoFIdx = src.getEdgeDoFFunction().getCellDataID();
PetscInt* srcVertexIndices = cell.getData( srcVertexDoFIdx )->getPointer( level );
PetscInt* dstVertexIndices = cell.getData( dstVertexDoFIdx )->getPointer( level );
matIdx_t* srcVertexIndices = cell.getData( srcVertexDoFIdx )->getPointer( level );
matIdx_t* dstVertexIndices = cell.getData( dstVertexDoFIdx )->getPointer( level );
PetscInt* srcEdgeIndices = cell.getData( srcEdgeDoFIdx )->getPointer( level );
PetscInt* dstEdgeIndices = cell.getData( dstEdgeDoFIdx )->getPointer( level );
matIdx_t* srcEdgeIndices = cell.getData( srcEdgeDoFIdx )->getPointer( level );
matIdx_t* dstEdgeIndices = cell.getData( dstEdgeDoFIdx )->getPointer( level );
// loop over micro-cells
for ( const auto& cType : celldof::allCellTypes )
......@@ -731,17 +729,17 @@ void P2ElementwiseOperator< P2Form >::assembleLocalMatrix( const std::shared_ptr
indexing::IndexIncrement offset;
// get hold of the actual numerical data in the two functions
PrimitiveDataID< FunctionMemory< PetscInt >, Face > dstVertexDoFIdx = dst.getVertexDoFFunction().getFaceDataID();
PrimitiveDataID< FunctionMemory< PetscInt >, Face > srcVertexDoFIdx = src.getVertexDoFFunction().getFaceDataID();
PrimitiveDataID< FunctionMemory< matIdx_t >, Face > dstVertexDoFIdx = dst.getVertexDoFFunction().getFaceDataID();
PrimitiveDataID< FunctionMemory< matIdx_t >, Face > srcVertexDoFIdx = src.getVertexDoFFunction().getFaceDataID();
PrimitiveDataID< FunctionMemory< PetscInt >, Face > dstEdgeDoFIdx = dst.getEdgeDoFFunction().getFaceDataID();
PrimitiveDataID< FunctionMemory< PetscInt >, Face > srcEdgeDoFIdx = src.getEdgeDoFFunction().getFaceDataID();
PrimitiveDataID< FunctionMemory< matIdx_t >, Face > dstEdgeDoFIdx = dst.getEdgeDoFFunction().getFaceDataID();
PrimitiveDataID< FunctionMemory< matIdx_t >, Face > srcEdgeDoFIdx = src.getEdgeDoFFunction().getFaceDataID();
PetscInt* srcVertexIndices = face.getData( srcVertexDoFIdx )->getPointer( level );
PetscInt* dstVertexIndices = face.getData( dstVertexDoFIdx )->getPointer( level );
matIdx_t* srcVertexIndices = face.getData( srcVertexDoFIdx )->getPointer( level );
matIdx_t* dstVertexIndices = face.getData( dstVertexDoFIdx )->getPointer( level );
PetscInt* srcEdgeIndices = face.getData( srcEdgeDoFIdx )->getPointer( level );
PetscInt* dstEdgeIndices = face.getData( dstEdgeDoFIdx )->getPointer( level );
matIdx_t* srcEdgeIndices = face.getData( srcEdgeDoFIdx )->getPointer( level );
matIdx_t* dstEdgeIndices = face.getData( dstEdgeDoFIdx )->getPointer( level );
// now loop over micro-faces of macro-face
for ( yIdx = 0; yIdx < rowsize - 2; ++yIdx )
......@@ -808,10 +806,10 @@ void P2ElementwiseOperator< P2Form >::localMatrixAssembly2D( const std::shared_p
const uint_t xIdx,
const uint_t yIdx,
const P2Elements::P2Element& element,
const PetscInt* const srcVertexIdx,
const PetscInt* const srcEdgeIdx,
const PetscInt* const dstVertexIdx,
const PetscInt* const dstEdgeIdx ) const
const matIdx_t* const srcVertexIdx,
const matIdx_t* const srcEdgeIdx,
const matIdx_t* const dstVertexIdx,
const matIdx_t* const dstEdgeIdx ) const
{
Matrix6r elMat;
......@@ -831,7 +829,7 @@ void P2ElementwiseOperator< P2Form >::localMatrixAssembly2D( const std::shared_p
// assemble local element matrix
form.setGeometryMap( face.getGeometryMap() );
form.integrateAll( { v0, v1, v2 }, elMat );
form.integrateAll( {v0, v1, v2}, elMat );
// determine global indices of our local DoFs (note the tweaked ordering to go along with FEniCS indexing)
dofDataIdx[0] = vertexdof::macroface::indexFromVertex( level, xIdx, yIdx, element[0] );
......@@ -877,10 +875,10 @@ void P2ElementwiseOperator< P2Form >::localMatrixAssembly3D( const std::shared_p
const uint_t level,
const indexing::Index& microCell,
const celldof::CellType cType,
const PetscInt* const srcVertexIdx,
const PetscInt* const srcEdgeIdx,
const PetscInt* const dstVertexIdx,
const PetscInt* const dstEdgeIdx ) const
const matIdx_t* const srcVertexIdx,
const matIdx_t* const srcEdgeIdx,
const matIdx_t* const dstVertexIdx,
const matIdx_t* const dstEdgeIdx ) const
{
// determine coordinates of vertices of micro-element
std::array< indexing::Index, 4 > verts = celldof::macrocell::getMicroVerticesFromMicroCell( microCell, cType );
......@@ -928,8 +926,6 @@ void P2ElementwiseOperator< P2Form >::localMatrixAssembly3D( const std::shared_p
mat->addValues( rowIdx, colIdx, blockMatData );
}
#endif
// P2ElementwiseLaplaceOperator
template class P2ElementwiseOperator<
P2FenicsForm< p2_diffusion_cell_integral_0_otherwise, p2_tet_diffusion_cell_integral_0_otherwise > >;
......
......@@ -113,7 +113,6 @@ class P2ElementwiseOperator : public Operator< P2Function< real_t >, P2Function<
WALBERLA_ABORT( "SOR not implemented for P2ElementwiseOperator." )
}
#ifdef HYTEG_BUILD_WITH_PETSC
/// Assemble operator as sparse matrix.
///
/// \param mat a sparse matrix proxy
......@@ -124,11 +123,10 @@ class P2ElementwiseOperator : public Operator< P2Function< real_t >, P2Function<
///
/// \note src and dst are legal to and often will be the same function object
void assembleLocalMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P2Function< PetscInt >& src,
const P2Function< PetscInt >& dst,
const P2Function< matIdx_t >& src,
const P2Function< matIdx_t >& dst,
uint_t level,
DoFType flag ) const;
#endif
/// Assemble operator as sparse matrix.
///
......@@ -139,6 +137,7 @@ class P2ElementwiseOperator : public Operator< P2Function< real_t >, P2Function<
/// \param flag ignored
///
/// \note src and dst are legal to and often will be the same function object
/// \note we should rename assembleLocalMatrix() to toMatrix() and skip the delegation
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P2Function< matIdx_t >& src,
const P2Function< matIdx_t >& dst,
......@@ -204,28 +203,26 @@ class P2ElementwiseOperator : public Operator< P2Function< real_t >, P2Function<
real_t* const vertexData,
real_t* const edgeData );
#ifdef HYTEG_BUILD_WITH_PETSC
void localMatrixAssembly2D( const std::shared_ptr< SparseMatrixProxy >& mat,
const Face& face,
const uint_t level,
const uint_t xIdx,
const uint_t yIdx,
const P2Elements::P2Element& element,
const PetscInt* const srcVertexIdx,
const PetscInt* const srcEdgeIdx,
const PetscInt* const dstVertexIdx,
const PetscInt* const dstEdgeIdx ) const;
const matIdx_t* const srcVertexIdx,
const matIdx_t* const srcEdgeIdx,
const matIdx_t* const dstVertexIdx,
const matIdx_t* const dstEdgeIdx ) const;
void localMatrixAssembly3D( const std::shared_ptr< SparseMatrixProxy >& mat,
const Cell& cell,
const uint_t level,
const indexing::Index& microCell,
const celldof::CellType cType,
const PetscInt* const srcVertexIdx,
const PetscInt* const srcEdgeIdx,
const PetscInt* const dstVertexIdx,
const PetscInt* const dstEdgeIdx ) const;
#endif
const matIdx_t* const srcVertexIdx,
const matIdx_t* const srcEdgeIdx,
const matIdx_t* const dstVertexIdx,
const matIdx_t* const dstEdgeIdx ) const;
/// Trigger (re)computation of diagonal matrix entries (central operator weights)
/// Allocates the required memory if the function was not yet allocated.
......
......@@ -274,8 +274,8 @@ int main( int argc, char* argv[] )
// -----------------------
WALBERLA_LOG_INFO_ON_ROOT( "" << rule << "\n ELEMENTWISE OPERATORS\n" << rule );
testAssembly< P1ElementwiseMassOperator >( storage, level, "P1ElementwiseOperator" );
testAssembly< P2ElementwiseMassOperator >( storage, level, "P2ElementwiseOperator" );
testAssembly_newAPI< P1ElementwiseMassOperator >( storage, level, "P1ElementwiseOperator" );
testAssembly_newAPI< P2ElementwiseMassOperator >( storage, level, "P2ElementwiseOperator" );
auto p2MassFormHyTeG = std::make_shared< P2Form_mass >();
std::shared_ptr< P2RowSumForm > lumpedMassFormP2HyTeG = std::make_shared< P2RowSumForm >( p2MassFormHyTeG );
......@@ -357,6 +357,7 @@ int main( int argc, char* argv[] )
// ------------------
// Block Operators
// ------------------
// #define FUTURE_OR_WIP
#ifdef FUTURE_OR_WIP
WALBERLA_LOG_INFO_ON_ROOT( "" << rule << "\n BLOCK OPERATORS\n" << rule );
......
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