diff --git a/hyteg_integration_tests/src/VektorAddition.cpp b/hyteg_integration_tests/src/VektorAddition.cpp new file mode 100644 index 0000000000000000000000000000000000000000..7f6e08b0ea4bdd9cdffb9373e74b34f78a39d2dd --- /dev/null +++ b/hyteg_integration_tests/src/VektorAddition.cpp @@ -0,0 +1,144 @@ +/* + * HyTeG Operator Generator + * Copyright (C) 2017-2024 Nils Kohl, Fabian Böhm, Daniel Bauer + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <https://www.gnu.org/licenses/>. + */ + +#include "core/DataTypes.h" + +#include "hyteg/elementwiseoperators/P1ElementwiseOperator.hpp" + +#include "Diffusion/TestOpDiffusion.hpp" +#include "OperatorGenerationTest.hpp" +#include "mixed_operator/VectorLaplaceOperator.hpp" +//#include "hyteg/p1functionspace/P1Function.hpp" + +using namespace hyteg; +using walberla::real_t; + +template < typename RefOpType, typename OpType> +void compareVecAddition( OperatorFactory< RefOpType > refOpFactory, + OperatorFactory< OpType > testOpFactory, + const uint_t level, + const StorageSetup& storageSetup, + const std::string& message//, + // const real_t thresholdOverMachineEps, + // const bool writeVTK = false ) + ) +{ + //not sure about all this + using RefType = typename RefOpType::srcType::valueType; + using TestType = typename TestOpType::srcType::valueType; + + using RefSrcFncType = typename RefOpType::srcType::template FunctionType< RefType >; + using RefDstFncType = typename RefOpType::dstType::template FunctionType< RefType >; + using TestSrcFncType = typename OpType::srcType::template FunctionType< TestType >; + using TestDstFncType = typename OpType::dstType::template FunctionType< TestType >; + + WALBERLA_LOG_INFO_ON_ROOT( message ) + + std::shared_ptr< hyteg::PrimitiveStorage > storage = storageSetup.createStorage(); + + TestSrcFncType src( "src", storage, level, level ); + TestDstFncType dst( "dst", storage, level, level ); + RefSrcFncType refSrc( "refSrc", storage, level, level ); + RefDstFncType refDst( "refDst", storage, level, level ); + + auto srcFunction = []( const hyteg::Point3D& x ) -> RefType { + return walberla::numeric_cast< RefType >( x[0] * x[0] * x[0] * x[0] * std::sinh( x[1] ) * std::cos( x[2] ) ); + }; + auto srcFunctionTest = []( const hyteg::Point3D& x ) -> TestType { + return walberla::numeric_cast< TestType >( x[0] * x[0] * x[0] * x[0] * std::sinh( x[1] ) * std::cos( x[2] ) ); + }; + + auto rand = []( const hyteg::Point3D& ) -> RefType { + return numeric_cast< RefType >( walberla::math::realRandom< RefType >() ); + }; + auto randTest = []( const hyteg::Point3D& ) -> TestType { + return numeric_cast< TestType >( walberla::math::realRandom< RefType >() ); + }; + // zusammenhang -> src[0] in VectorAddition Meshdaten + // init operand and fill destinations with random noise + src.interpolate( srcFunctionTest, level, All ); + dst.interpolate( randTest, level, All ); + refSrc.interpolate( srcFunction, level, All ); + refDst.interpolate( rand, level, All ); + + + RefOpType opRef = refOpFactory( storage, level, level ); + opRef.VectorAddition( refSrc, refDst, level, All, Replace );//???? + + OpType op = testOpFactory( storage, level, level ); + op.VectorAddition( src, dst, level, All, Replace ); + // ergebnis vergleichen + // wie vorgehen? in apply/invdiag werden auch nicht einträge miteinander verglichen? +} + + + + +template < typename RefOpType, typename OpType > +void compareVecAdd( const uint_t level, + const StorageSetup& storageSetup, + const std::string& message//, + // const real_t thresholdOverMachineEps, + // const bool writeVTK = false + ) +{ + OperatorFactory< RefOpType > makeRefOp = []( std::shared_ptr< PrimitiveStorage > storage, uint_t minLevel, uint_t maxLevel ) { + return RefOpType( storage, minLevel, maxLevel ); + }; + OperatorFactory< OpType > makeTestOp = []( std::shared_ptr< PrimitiveStorage > storage, uint_t minLevel, uint_t maxLevel ) { + return OpType( storage, minLevel, maxLevel ); + }; + compareVecAddition( makeRefOp, makeTestOp, level, storageSetup, message/*, thresholdOverMachineEps, writeVTK */); +} + +template < class Diagonal, class Reference, class... Operators > +void testOperator( const uint_t level, + const StorageSetup& storageSetup//, + // const real_t thresholdOverMachineEpsApply, + // const real_t thresholdOverMachineEpsInvDiag, + // const bool writeVTK = false + ) +{ + WALBERLA_LOG_INFO_ON_ROOT( "\n" << storageSetup.description() << " VectorAddition\n" ) + (compareVecAdd< Reference, Operators >( + level, storageSetup, typeid( Operators ).name()/*, thresholdOverMachineEpsApply, writeVTK */), + ... ); + +} + + +int main( int argc, char* argv[] ) +{ + walberla::MPIManager::instance()->initializeMPI( &argc, &argv ); + walberla::MPIManager::instance()->useWorldComm(); + + const uint_t level = 3; + StorageSetup storageSetup; + storageSetup = StorageSetup( + "quad_4el", MeshInfo::fromGmshFile( prependHyTeGMeshDir( "2D/quad_4el.msh" ) ), GeometryMap::Type::IDENTITY ); // for now just 2-dim + //TODO: - Operator erstellen (src vecs problem) + // - vektoraddition "durchführen" + // - referenz vektoraddition + // - ergebnis vergleichen + + testOperator< P1Function< real_t >, P1ElementwiseLaplaceOperator, operatorgeneration::TestOpDiffusion >( + level, storageSetup/*, 225, 9.0e6 */); + + + return EXIT_SUCCESS; +} \ No newline at end of file