diff --git a/tests/lbm/CMakeLists.txt b/tests/lbm/CMakeLists.txt index 87d27692a08414ffb67ab69f9c0d1f9ddea37d00..b67913654b775abb7fe5b9d246bf6be5bb52d7d3 100644 --- a/tests/lbm/CMakeLists.txt +++ b/tests/lbm/CMakeLists.txt @@ -106,6 +106,7 @@ waLBerla_compile_test( FILES codegen/LbCodeGenerationExample.cpp DEPENDS LbCodeG waLBerla_generate_target_from_python(NAME FluctuatingMRTGenerated FILE codegen/FluctuatingMRT.py OUT_FILES FluctuatingMRT_LatticeModel.cpp FluctuatingMRT_LatticeModel.h ) waLBerla_compile_test( FILES codegen/FluctuatingMRT.cpp DEPENDS FluctuatingMRTGenerated) +waLBerla_execute_test( NAME FluctuatingMRT COMMAND $<TARGET_FILE:FluctuatingMRT> ${CMAKE_CURRENT_SOURCE_DIR}/codegen/FluctuatingMRT.prm PROCESSES 1) waLBerla_generate_target_from_python(NAME FieldLayoutAndVectorizationTestGenerated FILE codegen/FieldLayoutAndVectorizationTest.py OUT_FILES FieldLayoutAndVectorizationTest_FZYX_Vec_LatticeModel.cpp FieldLayoutAndVectorizationTest_FZYX_Vec_LatticeModel.h diff --git a/tests/lbm/codegen/FluctuatingMRT.cpp b/tests/lbm/codegen/FluctuatingMRT.cpp index 6707c3b383e66c43da1ef3892c590762135e342f..a8a425d9f8dc9f52041e1f13fdebde0c41cd3ef6 100644 --- a/tests/lbm/codegen/FluctuatingMRT.cpp +++ b/tests/lbm/codegen/FluctuatingMRT.cpp @@ -18,19 +18,21 @@ // //====================================================================================================================== +/* This tests momentum balance for a fluctuating MRT LB with a constant + force applied +*/ + #include "blockforest/all.h" #include "core/all.h" #include "domain_decomposition/all.h" #include "field/all.h" #include "geometry/all.h" -#include "gui/all.h" #include "timeloop/all.h" -#include "lbm/field/PdfField.h" -#include "lbm/field/AddToStorage.h" #include "lbm/communication/PdfFieldPackInfo.h" -#include "lbm/gui/Connection.h" -#include "lbm/vtk/VTKOutput.h" +#include "lbm/field/Adaptors.h" +#include "lbm/field/AddToStorage.h" +#include "lbm/field/PdfField.h" #include "FluctuatingMRT_LatticeModel.h" @@ -59,18 +61,21 @@ int main( int argc, char ** argv ) // read parameters auto parameters = walberlaEnv.config()->getOneBlock( "Parameters" ); - const real_t omega = parameters.getParameter< real_t > ( "omega", real_c( 1.4 ) ); + const real_t omega = parameters.getParameter< real_t > ( "omega", real_c(1.4) ); + const real_t magic_number = 3. / 16; + const real_t omega_2 = (4 - 2 * omega) / (4 * magic_number * omega + 2 - omega); const Vector3<real_t> initialVelocity = parameters.getParameter< Vector3<real_t> >( "initialVelocity", Vector3<real_t>() ); - const uint_t timesteps = parameters.getParameter< uint_t > ( "timesteps", uint_c( 10 ) ); + const uint_t timesteps = parameters.getParameter< uint_t > ( "timesteps", uint_c(10) ); const real_t temperature = real_t(0.01); - const uint_t seed = uint_t(0); + const uint_t seed = uint_t(0); - const double remainingTimeLoggerFrequency = parameters.getParameter< double >( "remainingTimeLoggerFrequency", 3.0 ); // in seconds + const real_t remainingTimeLoggerFrequency = parameters.getParameter< double >( "remainingTimeLoggerFrequency", 3.0 ); // in seconds // create fields - BlockDataID forceFieldId = field::addToStorage<VectorField_T>( blocks, "Force", real_t( 0.0 ), field::fzyx); + real_t force = 2E-4; // Force to apply on each node on each axis + BlockDataID forceFieldId = field::addToStorage<VectorField_T>( blocks, "Force", force, field::fzyx ); - LatticeModel_T latticeModel = LatticeModel_T( forceFieldId, omega, 0, omega, omega, seed, temperature, uint_t(0) ); + LatticeModel_T latticeModel = LatticeModel_T( forceFieldId, omega, omega, omega_2, omega, seed, temperature, uint_t(0) ); BlockDataID pdfFieldId = lbm::addPdfFieldToStorage( blocks, "pdf field", latticeModel, initialVelocity, real_t(1), uint_t(1), field::fzyx ); BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >( blocks, "flag field" ); @@ -86,8 +91,8 @@ int main( int argc, char ** argv ) communication.addPackInfo( make_shared< lbm::PdfFieldPackInfo< LatticeModel_T > >( pdfFieldId ) ); // set the RNG counter to match the time step and propagate it to the fields' copies of the lattice model - timeloop.add() << BeforeFunction( [&](){ latticeModel.time_step_ = uint32_c(timeloop.getCurrentTimeStep()); }, "set RNG counter" ) - << Sweep( [&]( IBlock * block ){ + timeloop.add() << BeforeFunction( [&]() { latticeModel.time_step_ = uint32_c(timeloop.getCurrentTimeStep()); }, "set RNG counter" ) + << Sweep( [&]( IBlock * block ) { auto field = block->getData< PdfField_T >( pdfFieldId ); field->latticeModel().time_step_ = latticeModel.time_step_; }, "set RNG counter" ); @@ -101,25 +106,32 @@ int main( int argc, char ** argv ) "LBM stability check" ); // log remaining time - timeloop.addFuncAfterTimeStep( timing::RemainingTimeLogger( timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency ), "remaining time logger" ); + timeloop.addFuncAfterTimeStep( timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency), "remaining time logger" ); - // add VTK output to time loop - lbm::VTKOutput< LatticeModel_T, FlagField_T >::addToTimeloop( timeloop, blocks, walberlaEnv.config(), pdfFieldId, flagFieldId, fluidFlagUID ); + auto densityAdaptorId = field::addFieldAdaptor<lbm::Adaptor<LatticeModel_T>::Density> ( blocks, pdfFieldId, "DensityAdaptor" ); + auto velocityAdaptorId = field::addFieldAdaptor<lbm::Adaptor<LatticeModel_T>::VelocityVector>( blocks, pdfFieldId, "VelocityAdaptor" ); - // create adaptors, so that the GUI also displays density and velocity - // adaptors are like fields with the difference that they do not store values - // but calculate the values based on other fields ( here the PdfField ) - field::addFieldAdaptor<lbm::Adaptor<LatticeModel_T>::Density> ( blocks, pdfFieldId, "DensityAdaptor" ); - field::addFieldAdaptor<lbm::Adaptor<LatticeModel_T>::VelocityVector>( blocks, pdfFieldId, "VelocityAdaptor" ); + timeloop.run(); - if( parameters.getParameter<bool>( "useGui", false ) ) + // Calculate momentum + Vector3<real_t> momentum; // observed momentum + int count = 0; // count of lb nodes traversed + for (auto block = blocks->begin(); block != blocks->end(); ++block) { - GUI gui ( timeloop, blocks, argc, argv ); - lbm::connectToGui<LatticeModel_T> ( gui ); - gui.run(); + auto v = block->getData< lbm::Adaptor< LatticeModel_T >::VelocityVector >(velocityAdaptorId); + auto rho = block->getData< lbm::Adaptor< LatticeModel_T >::Density >(densityAdaptorId); + WALBERLA_FOR_ALL_CELLS_XYZ_OMP(v, omp critical, { + momentum += rho->get(x, y, z) * v->get(x, y, z); + count++; + }); } - else - timeloop.run(); + + // check + real_t expected_momentum = real_c(count) * force * real_c(timesteps); + printf("%g %g %g | %g\n", momentum[0], momentum[1], momentum[2], expected_momentum); + WALBERLA_CHECK_FLOAT_EQUAL(momentum[0], expected_momentum); + WALBERLA_CHECK_FLOAT_EQUAL(momentum[1], expected_momentum); + WALBERLA_CHECK_FLOAT_EQUAL(momentum[2], expected_momentum); return EXIT_SUCCESS; } diff --git a/tests/lbm/codegen/FluctuatingMRT.prm b/tests/lbm/codegen/FluctuatingMRT.prm new file mode 100644 index 0000000000000000000000000000000000000000..fbacd24099274a9edc5adc2d486e55298158ece4 --- /dev/null +++ b/tests/lbm/codegen/FluctuatingMRT.prm @@ -0,0 +1,18 @@ + +Parameters +{ + omega 0.5; + timesteps 2; + + useGui 0; + remainingTimeLoggerFrequency 3; // in seconds +} + +DomainSetup +{ + blocks < 1, 1, 1 >; + cellsPerBlock < 30, 30, 10 >; + periodic < 1, 1, 1 >; +} + +