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

Enabled all possible test cases in PetscMatrixAssemblyTest.

All operators that are either currently refactored (P*__Operator_new), including
polynomial/surrogate, variable ops, and composites that are built upon those
are not yet equipped with the ::toMatrix() method.
parent fff77f0b
Pipeline #34476 failed with stages
in 16 minutes and 12 seconds
......@@ -43,6 +43,7 @@ class P1EpsilonStokesOperator : public Operator< P1StokesFunction< real_t >, P1S
void apply( const P1StokesFunction< real_t >& src, const P1StokesFunction< real_t >& dst, size_t level, DoFType flag ) const
{
WALBERLA_CHECK( !src.uvw[0].getStorage()->hasGlobalCells(), "P1EpsilonStokesOperator not implemented for 3D." );
WALBERLA_ASSERT_NOT_IDENTICAL( std::addressof( src ), std::addressof( dst ) );
A_uu.apply( src.uvw[0], dst.uvw[0], level, flag, Replace );
......@@ -58,6 +59,28 @@ class P1EpsilonStokesOperator : public Operator< P1StokesFunction< real_t >, P1S
pspg.apply( src.p, dst.p, level, flag | DirichletBoundary, Add );
}
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P1StokesFunction< matIdx_t >& src,
const P1StokesFunction< matIdx_t >& dst,
size_t level,
DoFType flag ) const
{
WALBERLA_CHECK( !src.uvw[0].getStorage()->hasGlobalCells(), "P1EpsilonStokesOperator not implemented for 3D." );
A_uu.toMatrix( mat, src.uvw[0], dst.uvw[0], level, flag );
A_uv.toMatrix( mat, src.uvw[1], dst.uvw[0], level, flag );
divT_x.toMatrix( mat, src.p, dst.uvw[0], level, flag );
A_vu.toMatrix( mat, src.uvw[0], dst.uvw[1], level, flag );
A_vv.toMatrix( mat, src.uvw[1], dst.uvw[1], level, flag );
divT_y.toMatrix( mat, src.p, dst.uvw[1], level, flag );
div_x.toMatrix( mat, src.uvw[0], dst.p, level, flag | DirichletBoundary );
div_y.toMatrix( mat, src.uvw[1], dst.p, level, flag | DirichletBoundary );
pspg.toMatrix( mat, src.p, dst.p, level, flag | DirichletBoundary );
}
P1ConstantEpsilonOperator_11 A_uu;
P1ConstantEpsilonOperator_12 A_uv;
P1ConstantEpsilonOperator_21 A_vu;
......
......@@ -54,6 +54,18 @@ class P1ScalarToP2VectorOperator : public Operator< P1Function< real_t >, P2Vect
operZ.apply( src, dst[2], level, flag, updateType );
}
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P1Function< matIdx_t >& src,
const P2VectorFunction< matIdx_t >& dst,
size_t level,
DoFType flag ) const
{
operX.toMatrix( mat, src, dst[0], level, flag );
operY.toMatrix( mat, src, dst[1], level, flag );
if ( dst.getDimension() == 3 )
operZ.toMatrix( mat, src, dst[2], level, flag );
}
private:
operX_t operX;
operY_t operY;
......
......@@ -19,18 +19,18 @@
*/
#pragma once
#include "hyteg/p2functionspace/P2Function.hpp"
#include "hyteg/p2functionspace/P2Elements.hpp"
#include "hyteg/p2functionspace/P2Function.hpp"
#ifdef _MSC_VER
# pragma warning(push, 0)
#pragma warning( push, 0 )
#endif
#include "hyteg/forms/form_fenics_base/P2FenicsForm.hpp"
#include "hyteg/forms/form_fenics_base/P1ToP2FenicsForm.hpp"
#include "hyteg/forms/form_fenics_base/P2FenicsForm.hpp"
#ifdef _MSC_VER
# pragma warning(pop)
#pragma warning( pop )
#endif
#include "hyteg/mixedoperators/VertexDoFToEdgeDoFOperator/VertexDoFToEdgeDoFOperator.hpp"
......@@ -40,24 +40,19 @@ namespace hyteg {
using walberla::real_t;
template< class P1ToP2Form >
class P1ToP2ConstantOperator : public Operator<P1Function < real_t>, P2Function<real_t> > {
template < class P1ToP2Form >
class P1ToP2ConstantOperator : public Operator< P1Function< real_t >, P2Function< real_t > >
{
public:
P1ToP2ConstantOperator( const std::shared_ptr< PrimitiveStorage >& storage, size_t minLevel, size_t maxLevel )
: Operator( storage, minLevel, maxLevel )
, vertexToVertex( storage, minLevel, maxLevel )
, vertexToEdge( storage, minLevel, maxLevel )
{}
P1ToP2ConstantOperator(const std::shared_ptr< PrimitiveStorage > & storage, size_t minLevel, size_t maxLevel)
: Operator(storage, minLevel, maxLevel), vertexToVertex(storage, minLevel, maxLevel),
vertexToEdge(storage, minLevel, maxLevel)
{
}
P1ConstantOperator< P1ToP2Form > const& getVertexToVertexOpr() const { return vertexToVertex; }
P1ConstantOperator< P1ToP2Form > const & getVertexToVertexOpr() const {
return vertexToVertex;
}
VertexDoFToEdgeDoFOperator< P1ToP2Form > const & getVertexToEdgeOpr() const {
return vertexToEdge;
}
VertexDoFToEdgeDoFOperator< P1ToP2Form > const& getVertexToEdgeOpr() const { return vertexToEdge; }
void apply( const P1Function< real_t >& src,
const P2Function< real_t >& dst,
......@@ -67,22 +62,35 @@ class P1ToP2ConstantOperator : public Operator<P1Function < real_t>, P2Function<
{
vertexToVertex.apply( src, dst.getVertexDoFFunction(), level, flag, updateType );
vertexToEdge.apply( src, dst.getEdgeDoFFunction(), level, flag, updateType );
}
void smooth_gs(P1Function< real_t > & dst, P2Function< real_t > & rhs, size_t level, DoFType flag)
{
WALBERLA_ABORT("not implemented");
}
}
private:
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P1Function< matIdx_t >& src,
const P2Function< matIdx_t >& dst,
size_t level,
DoFType flag ) const
{
vertexToVertex.toMatrix( mat, src, dst.getVertexDoFFunction(), level, flag );
vertexToEdge.toMatrix( mat, src, dst.getEdgeDoFFunction(), level, flag );
}
P1ConstantOperator< P1ToP2Form > vertexToVertex;
VertexDoFToEdgeDoFOperator< P1ToP2Form > vertexToEdge;
void smooth_gs( P1Function< real_t >& dst, P2Function< real_t >& rhs, size_t level, DoFType flag )
{
WALBERLA_ABORT( "not implemented" );
}
private:
P1ConstantOperator< P1ToP2Form > vertexToVertex;
VertexDoFToEdgeDoFOperator< P1ToP2Form > vertexToEdge;
};
typedef P1ToP2ConstantOperator< P1ToP2FenicsForm< p1_to_p2_divt_cell_integral_0_otherwise, p1_to_p2_tet_divt_tet_cell_integral_0_otherwise > > P1ToP2ConstantDivTxOperator;
typedef P1ToP2ConstantOperator< P1ToP2FenicsForm< p1_to_p2_divt_cell_integral_1_otherwise, p1_to_p2_tet_divt_tet_cell_integral_1_otherwise > > P1ToP2ConstantDivTyOperator;
typedef P1ToP2ConstantOperator< P1ToP2FenicsForm< fenics::NoAssemble, p1_to_p2_tet_divt_tet_cell_integral_2_otherwise > > P1ToP2ConstantDivTzOperator;
typedef P1ToP2ConstantOperator<
P1ToP2FenicsForm< p1_to_p2_divt_cell_integral_0_otherwise, p1_to_p2_tet_divt_tet_cell_integral_0_otherwise > >
P1ToP2ConstantDivTxOperator;
typedef P1ToP2ConstantOperator<
P1ToP2FenicsForm< p1_to_p2_divt_cell_integral_1_otherwise, p1_to_p2_tet_divt_tet_cell_integral_1_otherwise > >
P1ToP2ConstantDivTyOperator;
typedef P1ToP2ConstantOperator< P1ToP2FenicsForm< fenics::NoAssemble, p1_to_p2_tet_divt_tet_cell_integral_2_otherwise > >
P1ToP2ConstantDivTzOperator;
}
} // namespace hyteg
......@@ -19,18 +19,18 @@
*/
#pragma once
#include "hyteg/p2functionspace/P2Function.hpp"
#include "hyteg/p2functionspace/P2Elements.hpp"
#include "hyteg/p2functionspace/P2Function.hpp"
#ifdef _MSC_VER
# pragma warning(push, 0)
#pragma warning( push, 0 )
#endif
#include "hyteg/forms/form_fenics_base/P2FenicsForm.hpp"
#include "hyteg/forms/form_fenics_base/P2ToP1FenicsForm.hpp"
#ifdef _MSC_VER
# pragma warning(pop)
#pragma warning( pop )
#endif
#include "hyteg/mixedoperators/EdgeDoFToVertexDoFOperator/EdgeDoFToVertexDoFOperator.hpp"
......@@ -40,24 +40,19 @@ namespace hyteg {
using walberla::real_t;
template< class P2ToP1Form >
class P2ToP1ConstantOperator : public Operator<P2Function < real_t>, P1Function<real_t> > {
template < class P2ToP1Form >
class P2ToP1ConstantOperator : public Operator< P2Function< real_t >, P1Function< real_t > >
{
public:
P2ToP1ConstantOperator( const std::shared_ptr< PrimitiveStorage >& storage, size_t minLevel, size_t maxLevel )
: Operator( storage, minLevel, maxLevel )
, vertexToVertex( storage, minLevel, maxLevel )
, edgeToVertex( storage, minLevel, maxLevel )
{}
P2ToP1ConstantOperator(const std::shared_ptr< PrimitiveStorage > & storage, size_t minLevel, size_t maxLevel)
: Operator(storage, minLevel, maxLevel),
vertexToVertex(storage, minLevel, maxLevel),
edgeToVertex(storage, minLevel, maxLevel)
{
}
P1ConstantOperator< P2ToP1Form > const& getVertexToVertexOpr() const { return vertexToVertex; }
P1ConstantOperator< P2ToP1Form > const & getVertexToVertexOpr() const {
return vertexToVertex;
}
EdgeDoFToVertexDoFOperator< P2ToP1Form > const & getEdgeToVertexOpr() const {
return edgeToVertex;
}
EdgeDoFToVertexDoFOperator< P2ToP1Form > const& getEdgeToVertexOpr() const { return edgeToVertex; }
void apply( const P2Function< real_t >& src,
const P1Function< real_t >& dst,
......@@ -69,20 +64,33 @@ class P2ToP1ConstantOperator : public Operator<P2Function < real_t>, P1Function<
edgeToVertex.apply( src.getEdgeDoFFunction(), dst, level, flag, Add );
}
void smooth_gs(P2Function< real_t > & dst, P1Function< real_t > & rhs, size_t level, DoFType flag)
{
WALBERLA_ABORT("not implemented");
}
private:
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P2Function< matIdx_t >& src,
const P1Function< matIdx_t >& dst,
size_t level,
DoFType flag ) const
{
vertexToVertex.toMatrix( mat, src.getVertexDoFFunction(), dst, level, flag );
edgeToVertex.toMatrix( mat, src.getEdgeDoFFunction(), dst, level, flag );
}
P1ConstantOperator< P2ToP1Form > vertexToVertex;
EdgeDoFToVertexDoFOperator< P2ToP1Form > edgeToVertex;
void smooth_gs( P2Function< real_t >& dst, P1Function< real_t >& rhs, size_t level, DoFType flag )
{
WALBERLA_ABORT( "not implemented" );
}
private:
P1ConstantOperator< P2ToP1Form > vertexToVertex;
EdgeDoFToVertexDoFOperator< P2ToP1Form > edgeToVertex;
};
typedef P2ToP1ConstantOperator< P2ToP1FenicsForm< p2_to_p1_div_cell_integral_0_otherwise, p2_to_p1_tet_div_tet_cell_integral_0_otherwise > > P2ToP1ConstantDivxOperator;
typedef P2ToP1ConstantOperator< P2ToP1FenicsForm< p2_to_p1_div_cell_integral_1_otherwise, p2_to_p1_tet_div_tet_cell_integral_1_otherwise > > P2ToP1ConstantDivyOperator;
typedef P2ToP1ConstantOperator< P2ToP1FenicsForm< fenics::NoAssemble, p2_to_p1_tet_div_tet_cell_integral_2_otherwise > > P2ToP1ConstantDivzOperator;
typedef P2ToP1ConstantOperator<
P2ToP1FenicsForm< p2_to_p1_div_cell_integral_0_otherwise, p2_to_p1_tet_div_tet_cell_integral_0_otherwise > >
P2ToP1ConstantDivxOperator;
typedef P2ToP1ConstantOperator<
P2ToP1FenicsForm< p2_to_p1_div_cell_integral_1_otherwise, p2_to_p1_tet_div_tet_cell_integral_1_otherwise > >
P2ToP1ConstantDivyOperator;
typedef P2ToP1ConstantOperator< P2ToP1FenicsForm< fenics::NoAssemble, p2_to_p1_tet_div_tet_cell_integral_2_otherwise > >
P2ToP1ConstantDivzOperator;
}
} // namespace hyteg
......@@ -48,13 +48,25 @@ class P2VectorToP1ScalarOperator : public Operator< P2VectorFunction< real_t >,
DoFType flag,
UpdateType updateType = Replace ) const
{
std::array< UpdateType, 3 > ut = {updateType, Add, Add};
std::array< UpdateType, 3 > ut = { updateType, Add, Add };
operX.apply( src[0], dst, level, flag, ut[0] );
operY.apply( src[1], dst, level, flag, ut[1] );
if ( src.getDimension() == 3 )
operZ.apply( src[2], dst, level, flag, ut[2] );
}
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P2VectorFunction< matIdx_t >& src,
const P1Function< matIdx_t >& dst,
size_t level,
DoFType flag ) const
{
operX.toMatrix( mat, src[0], dst, level, flag );
operY.toMatrix( mat, src[1], dst, level, flag );
if ( dst.getDimension() == 3 )
operZ.toMatrix( mat, src[2], dst, level, flag );
}
private:
operX_t operX;
operY_t operY;
......
......@@ -317,10 +317,10 @@ int main( int argc, char* argv[] )
typedef EdgeDoFToVertexDoFOperator< P2FenicsForm< p2_mass_cell_integral_0_otherwise, p2_tet_mass_cell_integral_0_otherwise > >
E2VMassOperator;
testAssembly< E2VMassOperator >( storage, level, "EdgeDoFToVertexDoFOperator", false );
testAssembly< E2VMassOperator >( storage, level, "EdgeDoFToVertexDoFOperator" );
typedef VertexDoFToEdgeDoFOperator< P2FenicsForm< p2_mass_cell_integral_0_otherwise, p2_tet_mass_cell_integral_0_otherwise > >
V2EMassOperator;
testAssembly< V2EMassOperator >( storage, level, "VertexDoFToEdgeDoFOperator", false );
testAssembly< V2EMassOperator >( storage, level, "VertexDoFToEdgeDoFOperator" );
// -----------------------
// Elementwise Operators
......@@ -358,7 +358,7 @@ int main( int argc, char* argv[] )
testAssembly< P1StokesBlockLaplaceOperator >( storage, level, "P1StokesBlockLaplaceOperator" );
testAssembly< P1BlendingStokesOperator >( storage, level, "P1BlendingStokesOperator", false );
testAssembly< P1EpsilonStokesOperator >( storage, level, "P1EpsilonStokesOperator", false );
testAssembly< P1EpsilonStokesOperator >( storage, level, "P1EpsilonStokesOperator" );
testAssembly< P2P2UnstableStokesOperator >( storage, level, "P2P2UnstableStokesOperator" );
testAssembly< P2P1BlendingTaylorHoodStokesOperator >( storage, level, "P2P1BlendingTaylorHoodStokesOperator", false );
testAssembly< P2P1SurrogateTaylorHoodStokesOperator >( storage, level, level, "P2P1SurrogateTaylorHoodStokesOperator", false );
......@@ -368,24 +368,24 @@ int main( int argc, char* argv[] )
// ----------------------------
WALBERLA_LOG_INFO_ON_ROOT( "" << rule << "\n SCALAR TO VECTOR OPERATORS\n" << rule );
testAssembly< P1ToP2ElementwiseBlendingDivTOperator >( storage, level, "P1ScalarToP2VectorOperator", false );
testAssembly< P1ToP2ConstantDivTxOperator >( storage, level, "P1ToP2ConstantOperator", false );
testAssembly< P1ToP2ElementwiseBlendingDivTOperator >( storage, level, "P1ScalarToP2VectorOperator" );
testAssembly< P1ToP2ConstantDivTxOperator >( storage, level, "P1ToP2ConstantOperator" );
// ----------------------------
// Vector To Scalar Operators
// ----------------------------
WALBERLA_LOG_INFO_ON_ROOT( "" << rule << "\n VECTOR TO SCALAR OPERATORS\n" << rule );
testAssembly< P2ToP1ElementwiseBlendingDivOperator >( storage, level, "P2ScalarToP1VectorOperator", false );
testAssembly< P2ToP1ConstantDivxOperator >( storage, level, "P2ToP1ConstantOperator", false );
testAssembly< P2ToP1ElementwiseBlendingDivOperator >( storage, level, "P2ScalarToP1VectorOperator" );
testAssembly< P2ToP1ConstantDivxOperator >( storage, level, "P2ToP1ConstantOperator" );
// ------------------
// Vector Operators
// ------------------
WALBERLA_LOG_INFO_ON_ROOT( "" << rule << "\n VECTOR TO VECTOR OPERATORS\n" << rule );
testAssembly< P1ConstantVectorMassOperator >( storage, level, "P1ConstantVectorMassOperator", false );
testAssembly< P2ConstantVectorMassOperator >( storage, level, "P2ConstantVectorMassOperator", false );
testAssembly< P2ElementwiseBlendingVectorMassOperator >( storage, level, "P2ElementwiseBlendingVectorMassOperator", false );
testAssembly< P1ConstantVectorMassOperator >( storage, level, "P1ConstantVectorMassOperator" );
testAssembly< P2ConstantVectorMassOperator >( storage, level, "P2ConstantVectorMassOperator" );
testAssembly< P2ElementwiseBlendingVectorMassOperator >( storage, level, "P2ElementwiseBlendingVectorMassOperator" );
// WALBERLA_LOG_INFO_ON_ROOT( "Skipping: (since assembly doesn't compile)" );
// WALBERLA_LOG_INFO_ON_ROOT( " * [all of this kind]" );
......
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