diff --git a/apps/benchmarks/FluidParticleCoupling/CMakeLists.txt b/apps/benchmarks/FluidParticleCoupling/CMakeLists.txt
index a22871bb3eb9db47c10b47ca3a08f7812d5d5068..c42c1e497a02f30875f566b62cfeac33bfc21338 100644
--- a/apps/benchmarks/FluidParticleCoupling/CMakeLists.txt
+++ b/apps/benchmarks/FluidParticleCoupling/CMakeLists.txt
@@ -8,10 +8,13 @@ waLBerla_python_file_generates(GeneratedLBMWithForce.py
GeneratedLBMWithForce.cpp GeneratedLBMWithForce.h
)
+
+if( WALBERLA_BUILD_WITH_CODEGEN )
+
waLBerla_add_executable ( NAME SphereWallCollision FILES SphereWallCollision.cpp GeneratedLBM.py
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
-waLBerla_add_executable ( NAME SettlingSphereTenCate FILES SettlingSphereInBox.cpp GeneratedLBM.py
+waLBerla_add_executable ( NAME SettlingSphereInBox FILES SettlingSphereInBox.cpp GeneratedLBM.py
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
waLBerla_add_executable ( NAME SphereMovingWithPrescribedVelocity FILES SphereMovingWithPrescribedVelocity.cpp GeneratedLBM.py
@@ -26,8 +29,35 @@ waLBerla_add_executable ( NAME DragForceSphere FILES DragForceSphere.cpp Generat
waLBerla_add_executable ( NAME ForcesOnSphereNearPlane FILES ForcesOnSphereNearPlane.cpp GeneratedLBM.py
DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
-waLBerla_add_executable ( NAME ObliqueDryCollision FILES ObliqueDryCollision.cpp
- DEPENDS blockforest core mesa_pd postprocessing )
-
waLBerla_add_executable ( NAME ObliqueWetCollision FILES ObliqueWetCollision.cpp GeneratedLBM.py
- DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
\ No newline at end of file
+ DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+
+else()
+
+waLBerla_add_executable ( NAME SphereWallCollision FILES SphereWallCollision.cpp
+ DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+
+waLBerla_add_executable ( NAME SettlingSphereInBox FILES SettlingSphereInBox.cpp
+ DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+
+waLBerla_add_executable ( NAME SphereMovingWithPrescribedVelocity FILES SphereMovingWithPrescribedVelocity.cpp
+ DEPENDS blockforest boundary core domain_decomposition field lbm mesa_pd lbm_mesapd_coupling postprocessing timeloop vtk )
+
+waLBerla_add_executable ( NAME LubricationForceEvaluation FILES LubricationForceEvaluation.cpp
+ DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+
+waLBerla_add_executable ( NAME DragForceSphere FILES DragForceSphere.cpp
+ DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+
+waLBerla_add_executable ( NAME ForcesOnSphereNearPlane FILES ForcesOnSphereNearPlane.cpp
+ DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+
+waLBerla_add_executable ( NAME ObliqueWetCollision FILES ObliqueWetCollision.cpp
+ DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+
+endif()
+
+
+
+waLBerla_add_executable ( NAME ObliqueDryCollision FILES ObliqueDryCollision.cpp
+ DEPENDS blockforest core mesa_pd postprocessing )
\ No newline at end of file
diff --git a/apps/benchmarks/FluidParticleCoupling/DragForceSphere.cpp b/apps/benchmarks/FluidParticleCoupling/DragForceSphere.cpp
index 2f269c5d944d021820fe9f63af86ff5c3671798d..a28651957a527509cc3575172614ff2242a6c869 100644
--- a/apps/benchmarks/FluidParticleCoupling/DragForceSphere.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/DragForceSphere.cpp
@@ -77,9 +77,9 @@
#include <iomanip>
#include <iostream>
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
#include "GeneratedLBMWithForce.h"
-
-#define USE_TRT_LIKE_LATTICE_MODEL
+#endif
namespace drag_force_sphere_mem
{
@@ -91,7 +91,13 @@ namespace drag_force_sphere_mem
using namespace walberla;
using walberla::uint_t;
+
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
using LatticeModel_T = lbm::GeneratedLBMWithForce;
+#else
+using ForceModel_T = lbm::force_model::LuoConstant;
+using LatticeModel_T = lbm::D3Q19< lbm::collision_model::D3Q19MRT, false, ForceModel_T>;
+#endif
using Stencil_T = LatticeModel_T::Stencil;
using PdfField_T = lbm::PdfField<LatticeModel_T>;
@@ -517,8 +523,7 @@ int main( int argc, char **argv )
///////////////////////
// ADD DATA TO BLOCKS //
////////////////////////
-
-
+
real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega );
real_t lambda_d = (useSRT) ? lambda_e : lbm::collision_model::TRT::lambda_d( omega, magicNumber );
real_t omegaBulk = (useSRT) ? lambda_e : lbm_mesapd_coupling::omegaBulkFromOmega(omega, bulkViscRateFactor);
@@ -527,30 +532,23 @@ int main( int argc, char **argv )
BlockDataID omegaBulkFieldID = field::addToStorage<ScalarField_T>( blocks, "omega bulk field", omegaBulk, field::fzyx, uint_t(0) );
// create the lattice model
-
- // generated MRT
-#ifdef USE_TRT_LIKE_LATTICE_MODEL
- WALBERLA_LOG_INFO_ON_ROOT("Using TRT-like lattice model!");
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
+ WALBERLA_LOG_INFO_ON_ROOT("Using generated TRT-like lattice model!");
WALBERLA_LOG_INFO_ON_ROOT(" - magic number " << magicNumber);
WALBERLA_LOG_INFO_ON_ROOT(" - omegaBulk = " << omegaBulk << ", bulk visc. = " << lbm_mesapd_coupling::bulkViscosityFromOmegaBulk(omegaBulk) << " (bvrf " << bulkViscRateFactor << ")");
WALBERLA_LOG_INFO_ON_ROOT(" - lambda_e " << lambda_e << ", lambda_d " << lambda_d << ", omegaBulk " << omegaBulk );
WALBERLA_LOG_INFO_ON_ROOT(" - use omega bulk adaption = " << useOmegaBulkAdaption << " (adaption layer size = " << adaptionLayerSize << ")");
LatticeModel_T latticeModel = LatticeModel_T(omegaBulkFieldID, setup.extForce, lambda_d, lambda_e);
- //LatticeModel_T latticeModel = LatticeModel_T(setup.extForce, omegaBulk, lambda_d, lambda_e);
- //LatticeModel_T latticeModel = LatticeModel_T(setup.extForce, lambda_e, lambda_d);
#else
- WALBERLA_LOG_INFO_ON_ROOT("Using other lattice model!");
- // generated KBC
- LatticeModel_T latticeModel = LatticeModel_T(setup.extForce, omega);
+ WALBERLA_LOG_INFO_ON_ROOT("Using waLBerla built-in MRT lattice model and ignoring omega bulk field since not supported!");
+ WALBERLA_LOG_INFO_ON_ROOT(" - magic number " << magicNumber);
+ WALBERLA_LOG_INFO_ON_ROOT(" - omegaBulk = " << omegaBulk << ", bulk visc. = " << lbm_mesapd_coupling::bulkViscosityFromOmegaBulk(omegaBulk) << " (bvrf " << bulkViscRateFactor << ")");
+ WALBERLA_LOG_INFO_ON_ROOT(" - lambda_e " << lambda_e << ", lambda_d " << lambda_d << ", omegaBulk " << omegaBulk );
+
+ LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::D3Q19MRT( omegaBulk, omegaBulk, lambda_d, lambda_e, lambda_e, lambda_d ), ForceModel_T( Vector3<real_t> ( setup.extForce, 0, 0 ) ));
#endif
- // set s_e and s_eps to omegaBulk
- // s_e allows to set a specific bulk viscosity that is different from the SRT/TRT default: nu_b = 2/3 nu_shear
- // Khirevich et al propose to set s_eps in the same way to avoid huge differences in s_e and s_eps (would be omega in SRT/TRT) that could lead to stability problems
- // to not break TRT properties, it is proposed to set Lambda_bulk = conste * Lambda_shear (Khirevich), where conste sets the ratio between bulk and shear viscosity-terms (1 in SRT/TRT -> nu_b = 2/3 nu_shear)
- // all in all, this allows to set the bulk viscosity but to retain the "TRT" properties, i.e. independent drag force of tau
- //LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::D3Q19MRT( omegaBulk, omegaBulk, lambda_d, lambda_e, lambda_e, lambda_d ), ForceModel_T( Vector3<real_t> ( setup.extForce, 0, 0 ) ));
// add PDF field
BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (fzyx)", latticeModel,
@@ -636,9 +634,9 @@ int main( int argc, char **argv )
using DragForceEval_T = DragForceEvaluator<ParticleAccessor_T>;
auto forceEval = make_shared< DragForceEval_T >( &timeloop, &setup, blocks, flagFieldID, pdfFieldID, accessor, sphereID );
- auto lbmSweep = LatticeModel_T::Sweep( pdfFieldID );
-
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
+ auto lbmSweep = LatticeModel_T::Sweep( pdfFieldID );
// splitting stream and collide is currently buggy in lbmpy with vectorization -> build without optimize for local host and vectorization
// collide LBM step
timeloop.add() << Sweep([&lbmSweep](IBlock * const block){lbmSweep.collide(block);}, "collide LB sweep" );
@@ -650,16 +648,21 @@ int main( int argc, char **argv )
// stream LBM step
timeloop.add() << Sweep([&lbmSweep](IBlock * const block){lbmSweep.stream(block);}, "stream LB sweep" )
<< AfterFunction( SharedFunctor< DragForceEval_T >( forceEval ), "drag force evaluation" );
+#else
+ auto lbmSweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag );
+
+ // collision sweep
+ timeloop.add() << Sweep( lbm::makeCollideSweep( lbmSweep ), "cell-wise LB sweep (collide)" );
-/*
// add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
timeloop.add() << BeforeFunction( optimizedPDFCommunicationScheme, "LBM Communication" )
<< Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingID ), "Boundary Handling" );
- // stream + collide LBM step
- timeloop.add() << Sweep( lbmSweep, "LB sweep" )
+ // streaming & force evaluation
+ timeloop.add() << Sweep( lbm::makeStreamSweep( lbmSweep ), "cell-wise LB sweep (stream)" )
<< AfterFunction( SharedFunctor< DragForceEval_T >( forceEval ), "drag force evaluation" );
-*/
+#endif
+
// resetting force
timeloop.addFuncAfterTimeStep( [ps,accessor](){ps->forEachParticle(false, mesa_pd::kernel::SelectAll(), *accessor, lbm_mesapd_coupling::ResetHydrodynamicForceTorqueKernel(),*accessor);}, "reset force on sphere");
diff --git a/apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp b/apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp
index eae674e9aa141bcee24ae4057b9da15a057a35f7..77524bf2da21e491336e7ada054fb3ac76483747 100644
--- a/apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/ForcesOnSphereNearPlane.cpp
@@ -78,10 +78,11 @@
#include "field/vtk/all.h"
#include "lbm/vtk/all.h"
-#include "GeneratedLBM.h"
#include <functional>
-#define USE_TRT_LIKE_LATTICE_MODEL
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
+#include "GeneratedLBM.h"
+#endif
namespace forces_on_sphere_near_plane
{
@@ -97,8 +98,11 @@ using walberla::uint_t;
// TYPEDEFS //
//////////////
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
using LatticeModel_T = lbm::GeneratedLBM;
-//using LatticeModel_T = lbm::D3Q19< lbm::collision_model::D3Q19MRT>;
+#else
+using LatticeModel_T = lbm::D3Q19< lbm::collision_model::D3Q19MRT>;
+#endif
using Stencil_T = LatticeModel_T::Stencil;
using PdfField_T = lbm::PdfField<LatticeModel_T>;
@@ -439,10 +443,6 @@ int main( int argc, char **argv )
const real_t domainWidth = real_t(16) * diameter; //y
const real_t domainHeight = real_t( 8) * diameter; //z
- //const real_t domainLength = real_t(3) * diameter; //x
- //const real_t domainWidth = real_t(3) * diameter; //y
- //const real_t domainHeight = real_t(8) * diameter; //z
-
Vector3<uint_t> domainSize( uint_c( std::ceil(domainLength)), uint_c( std::ceil(domainWidth)), uint_c( std::ceil(domainHeight)) );
@@ -559,15 +559,14 @@ int main( int argc, char **argv )
BlockDataID omegaBulkFieldID = field::addToStorage<ScalarField_T>( blocks, "omega bulk field", omegaBulk, field::fzyx, uint_t(0) );
// create the lattice model
-#ifdef USE_TRT_LIKE_LATTICE_MODEL
real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega );
real_t lambda_d = lbm::collision_model::TRT::lambda_d( omega, magicNumber );
- //LatticeModel_T latticeModel = LatticeModel_T(omegaBulk, lambda_d, lambda_e);
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
+ WALBERLA_LOG_INFO_ON_ROOT("Using generated TRT-like lattice model!");
LatticeModel_T latticeModel = LatticeModel_T(omegaBulkFieldID, lambda_d, lambda_e);
- //LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::D3Q19MRT( omegaBulk, omegaBulk, lambda_d, lambda_e, lambda_e, lambda_d ));
#else
- // generated KBC
- LatticeModel_T latticeModel = LatticeModel_T(omega);
+ WALBERLA_LOG_INFO_ON_ROOT("Using waLBerla built-in MRT lattice model and ignoring omega bulk field since not supported!");
+ LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::D3Q19MRT( omegaBulk, omegaBulk, lambda_d, lambda_e, lambda_e, lambda_d ));
#endif
// add PDF field
@@ -593,7 +592,7 @@ int main( int argc, char **argv )
ps->forEachParticle(false, lbm_mesapd_coupling::RegularParticlesSelector(), *accessor, movingParticleMappingKernel, *accessor, CLI_Flag);
// map planes
- ps->forEachParticle(false, lbm_mesapd_coupling::GlobalParticlesSelector(), *accessor, movingParticleMappingKernel, *accessor, CLI_Flag); //TODO: was SBB_Flag?
+ ps->forEachParticle(false, lbm_mesapd_coupling::GlobalParticlesSelector(), *accessor, movingParticleMappingKernel, *accessor, CLI_Flag);
// initialize Couette velocity profile in whole domain
if(initializeVelocityProfile)
@@ -628,10 +627,6 @@ int main( int argc, char **argv )
if( vtkIOFreq != uint_t(0) )
{
- // flag field (written only once in the first time step, ghost layers are also written)
- //auto flagFieldVTK = vtk::createVTKOutput_BlockData( blocks, "flag_field", vtkIOFreq, uint_t(0), false, baseFolderVTK );
- //flagFieldVTK->addCellDataWriter( make_shared< field::VTKWriter< FlagField_T > >( flagFieldID, "FlagField" ) );
- //timeloop.addFuncBeforeTimeStep( vtk::writeFiles( flagFieldVTK ), "VTK (flag field data)" );
auto pdfFieldVTK = vtk::createVTKOutput_BlockData( blocks, "fluid_field", vtkIOFreq, uint_t(0), false, baseFolderVTK );
@@ -639,16 +634,6 @@ int main( int argc, char **argv )
fluidFilter.addFlag( Fluid_Flag );
pdfFieldVTK->addCellInclusionFilter( fluidFilter );
- //AABB sliceAABB( real_t(0), real_c(domainSize[1])*real_t(0.5)-real_t(1), real_t(0),
- // real_c(domainSize[0]), real_c(domainSize[1])*real_t(0.5)+real_t(1), real_c(domainSize[2]) );
- //vtk::AABBCellFilter aabbSliceFilter( sliceAABB );
- //pdfFieldVTK->addCellInclusionFilter( aabbSliceFilter );
-
- //vtk::ChainedFilter combinedSliceFilter;
- //combinedSliceFilter.addFilter( fluidFilter );
- //combinedSliceFilter.addFilter( aabbSliceFilter );
- //pdfFieldVTK->addCellInclusionFilter( combinedSliceFilter );
-
pdfFieldVTK->addCellDataWriter( make_shared< lbm::VelocityVTKWriter< LatticeModel_T, float > >( pdfFieldID, "VelocityFromPDF" ) );
pdfFieldVTK->addCellDataWriter( make_shared< lbm::DensityVTKWriter < LatticeModel_T, float > >( pdfFieldID, "DensityFromPDF" ) );
@@ -656,14 +641,19 @@ int main( int argc, char **argv )
}
auto bhSweep = BoundaryHandling_T::getBlockSweep( boundaryHandlingID );
- auto lbmSweep = LatticeModel_T::Sweep( pdfFieldID );
// add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
timeloop.add() << BeforeFunction( optimizedPDFCommunicationScheme, "LBM Communication" )
<< Sweep(bhSweep, "Boundary Handling" );
// stream + collide LBM step
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
+ auto lbmSweep = LatticeModel_T::Sweep( pdfFieldID );
timeloop.add() << Sweep( lbmSweep, "LB sweep" );
+#else
+ auto lbmSweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag );
+ timeloop.add() << Sweep( makeSharedSweep( lbmSweep ), "cell-wise LB sweep" );
+#endif
// add force evaluation and logging
diff --git a/apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py b/apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py
index 37d1a35ee5e3c18f79764ec4308cd4f2a93bbbbc..43b08e8ec940daa31e4fa126d100a1076a59416c 100644
--- a/apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py
+++ b/apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py
@@ -7,182 +7,38 @@ from pystencils_walberla import CodeGeneration
from lbmpy.creationfunctions import create_lb_collision_rule
from lbmpy.moments import MOMENT_SYMBOLS, is_even, get_order
from lbmpy.stencils import get_stencil
-from lbmpy.maxwellian_equilibrium import get_moments_of_discrete_maxwellian_equilibrium
from lbmpy.forcemodels import *
-from collections import OrderedDict
-from lbmpy.methods import create_with_discrete_maxwellian_eq_moments
with CodeGeneration() as ctx:
-
- # methods: TRTlike, KBC, SRT, cumulant
- generatedMethod = "TRTlike"
-
- print("Generating " + generatedMethod + " LBM method")
-
- #
- # x,y,z = MOMENT_SYMBOLS
- #
- # e_sq = x**2 + y**2 + z**2
- #
- # moments = [0] * 19
- #
- # moments[0] = sp.sympify(1)
- # moments[1] = 19 * e_sq - 30
- # moments[2] = (21 * e_sq**2 - 53 * e_sq + 24) / 2
- #
- # moments[3] = x
- # moments[5] = y
- # moments[7] = z
- #
- # moments[4] = (5*e_sq - 9) * x
- # moments[6] = (5*e_sq - 9) * y
- # moments[8] = (5*e_sq - 9) * z
- #
- # moments[9] = 3 * x**2 - e_sq
- # moments[11] = y**2 - z**2
- #
- # moments[13] = x * y
- # moments[14] = y * z
- # moments[15] = x * z
- #
- # moments[10] = (3* e_sq - 5)*(3*x**2 - e_sq)
- # moments[12] = (3* e_sq - 5)*(y**2 - z**2)
- # moments[16] = (y**2 - z**2) * x
- # moments[17] = (z**2 - x**2) * y
- # moments[18] = (x**2 - y**2) * z
- #
- # m_eq = get_moments_of_discrete_maxwellian_equilibrium(stencil, moments, order=2, c_s_sq=sp.Rational(1,3))
- #
- # rr_dict = OrderedDict(zip(moments, omega))
- #
- # forcing = sp.symbols("forcing_:%d" % 3)
- # forcing=(sp.symbols("fx"),0,0)
- # forcemodel=Luo(forcing) #None
- # method = create_with_discrete_maxwellian_eq_moments(stencil, rr_dict, compressible=False, force_model=forcemodel)
- # at this stage, we have the MRT model from the LBM book!
-
- #this is the schiller / walberla MRT
-
- if generatedMethod == "TRTlike":
- stencil = get_stencil("D3Q19", 'walberla')
- omega = sp.symbols("omega_:%d" % len(stencil))
- method = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False)
-
- def modification_func(moment, eq, rate):
- omegaVisc = sp.Symbol("omega_visc")
- omegaBulk = ps.fields("omega_bulk: [3D]", layout='fzyx')# sp.Symbol("omega_bulk")
- omegaMagic = sp.Symbol("omega_magic")
- if get_order(moment) <= 1:
- return moment, eq, 0
- elif rate == omega[1]:
- return moment, eq, omegaBulk.center_vector
- elif rate == omega[2]:
- return moment, eq, omegaBulk.center_vector
- elif is_even(moment):
- return moment, eq, omegaVisc
- else:
- return moment, eq, omegaMagic
-
- my_method = create_lb_method_from_existing(method, modification_func)
-
- # optimizations
-
- collision_rule = create_lb_collision_rule(lb_method=my_method,
- optimization={"cse_global": True,
- "cse_pdfs": False,
- #"split": True
- } )
-
- generate_lattice_model(ctx, 'GeneratedLBM', collision_rule)
-
- elif generatedMethod == "KBC":
- methodName = 'trt-kbc-n4'
- #omega_field = ps.fields("omegaField(1): [3D]", layout='fzyx') omega_output_field=omega_field.center,
- collision_rule = create_lb_collision_rule(method=methodName,entropic=True, stencil=get_stencil("D3Q27", 'walberla'), compressible=True,
- optimization={"cse_global": True,
- "cse_pdfs": False,
- "split": True})
- generate_lattice_model(ctx, 'GeneratedLBM', collision_rule)
-
- elif generatedMethod == "SRT":
- methodName = 'srt'
- collision_rule = create_lb_collision_rule(method=methodName,stencil=get_stencil("D3Q19", 'walberla'),
+ stencil = get_stencil("D3Q19", 'walberla')
+ omega = sp.symbols("omega_:%d" % len(stencil))
+ method = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False)
+
+ def modification_func(moment, eq, rate):
+ omegaVisc = sp.Symbol("omega_visc")
+ omegaBulk = ps.fields("omega_bulk: [3D]", layout='fzyx')# sp.Symbol("omega_bulk")
+ omegaMagic = sp.Symbol("omega_magic")
+ if get_order(moment) <= 1:
+ return moment, eq, 0
+ elif rate == omega[1]:
+ return moment, eq, omegaBulk.center_vector
+ elif rate == omega[2]:
+ return moment, eq, omegaBulk.center_vector
+ elif is_even(moment):
+ return moment, eq, omegaVisc
+ else:
+ return moment, eq, omegaMagic
+
+ my_method = create_lb_method_from_existing(method, modification_func)
+
+ collision_rule = create_lb_collision_rule(lb_method=my_method,
optimization={"cse_global": True,
- "cse_pdfs": False,
- "split": True})
- generate_lattice_model(ctx, 'GeneratedLBM', collision_rule)
-
- elif generatedMethod == "cumulant":
-
- x,y,z = MOMENT_SYMBOLS
-
- cumulants = [0] * 27
-
- cumulants[0] = sp.sympify(1) #000
-
- cumulants[1] = x #100
- cumulants[2] = y #010
- cumulants[3] = z #001
-
- cumulants[4] = x*y #110
- cumulants[5] = x*z #101
- cumulants[6] = y*z #011
-
- cumulants[7] = x**2 - y**2 #200 - 020
- cumulants[8] = x**2 - z**2 #200 - 002
- cumulants[9] = x**2 + y**2 + z**2 #200 + 020 + 002
-
- cumulants[10] = x*y**2 + x*z**2 #120 + 102
- cumulants[11] = x**2*y + y*z**2 #210 + 012
- cumulants[12] = x**2*z + y**2*z #201 + 021
- cumulants[13] = x*y**2 - x*z**2 #120 - 102
- cumulants[14] = x**2*y - y*z**2 #210 - 012
- cumulants[15] = x**2*z - y**2*z #201 - 021
-
- cumulants[16] = x*y*z #111
-
- cumulants[17] = x**2*y**2 - 2*x**2*z**2 + y**2*z**2 # 220- 2*202 +022
- cumulants[18] = x**2*y**2 + x**2*z**2 - 2*y**2*z**2 # 220 + 202 - 2*002
- cumulants[19] = x**2*y**2 + x**2*z**2 + y**2*z**2 # 220 + 202 + 022
-
- cumulants[20] = x**2*y*z # 211
- cumulants[21] = x*y**2*z # 121
- cumulants[22] = x*y*z**2 # 112
-
- cumulants[23] = x**2*y**2*z # 221
- cumulants[24] = x**2*y*z**2 # 212
- cumulants[25] = x*y**2*z**2 # 122
-
- cumulants[26] = x**2*y**2*z**2 # 222
-
-
- from lbmpy.moments import get_order
- def get_relaxation_rate(cumulant, omega):
- if get_order(cumulant) <= 1:
- return 0
- elif get_order(cumulant) == 2 and cumulant != x**2 + y**2 + z**2:
- return omega
- else:
- return 1
-
- stencil = get_stencil("D3Q27", 'walberla')
-
- omega = sp.Symbol("omega")
- rr_dict = OrderedDict((c, get_relaxation_rate(c, omega))
- for c in cumulants)
-
- from lbmpy.methods import create_with_continuous_maxwellian_eq_moments
- my_method = create_with_continuous_maxwellian_eq_moments(stencil, rr_dict, cumulant=True, compressible=True)
-
- collision_rule = create_lb_collision_rule(lb_method=my_method,
- optimization={"cse_global": True,
- "cse_pdfs": False})
+ "cse_pdfs": False
+ } )
- generate_lattice_model(ctx, 'GeneratedLBM', collision_rule)
+ generate_lattice_model(ctx, 'GeneratedLBM', collision_rule)
- else:
- print("Invalid generated method string! " + generatedMethod)
diff --git a/apps/benchmarks/FluidParticleCoupling/GeneratedLBMWithForce.py b/apps/benchmarks/FluidParticleCoupling/GeneratedLBMWithForce.py
index 00200a679512ecfcce9dd57113cac2fcf116cbd0..665fc0226d3a30f66a01d5f46006dd83b481d7db 100644
--- a/apps/benchmarks/FluidParticleCoupling/GeneratedLBMWithForce.py
+++ b/apps/benchmarks/FluidParticleCoupling/GeneratedLBMWithForce.py
@@ -7,75 +7,35 @@ from pystencils_walberla import CodeGeneration
from lbmpy.creationfunctions import create_lb_collision_rule
from lbmpy.moments import MOMENT_SYMBOLS, is_even, get_order
from lbmpy.stencils import get_stencil
-from lbmpy.maxwellian_equilibrium import get_moments_of_discrete_maxwellian_equilibrium
from lbmpy.forcemodels import *
-from collections import OrderedDict
-from lbmpy.methods import create_with_discrete_maxwellian_eq_moments
with CodeGeneration() as ctx:
forcing=(sp.symbols("fx"),0,0)
- forcemodel=Luo(forcing) #None
-
- # methods: TRTlike, KBC, SRT
- generatedMethod = "TRTlike"
-
- if generatedMethod == "TRTlike":
- stencil = get_stencil("D3Q19", 'walberla')
- omega = sp.symbols("omega_:%d" % len(stencil))
-
- methodWithForce = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False,force_model=forcemodel)
-
- def modification_func(moment, eq, rate):
- omegaVisc = sp.Symbol("omega_visc")
- omegaBulk = ps.fields("omega_bulk: [3D]", layout='fzyx')# = sp.Symbol("omega_bulk")
- #omegaBulk = sp.Symbol("omega_bulk")
- omegaMagic = sp.Symbol("omega_magic")
- if get_order(moment) <= 1:
- return moment, eq, 0
- elif rate == omega[1]:
- return moment, eq, omegaBulk.center_vector
- elif rate == omega[2]:
- return moment, eq, omegaBulk.center_vector
- elif is_even(moment):
- return moment, eq, omegaVisc
- else:
- return moment, eq, omegaMagic
-
- my_methodWithForce = create_lb_method_from_existing(methodWithForce, modification_func)
-
- collision_rule = create_lb_collision_rule(lb_method=my_methodWithForce)
- # ,
- # optimization={"cse_global": True,
- # "cse_pdfs": False,
- # "split": True,
- # "vectorization":{'instruction_set': 'avx',
- # 'nontemporal': True,
- # 'assume_inner_stride_one': True}})
-
- generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule)
-
- elif generatedMethod == "SRT":
- collision_rule = create_lb_collision_rule(method='srt',stencil=get_stencil("D3Q19", 'walberla'), force_model=forcemodel,
- optimization={"cse_global": True,
- "cse_pdfs": False,
- #"split": True
- })
- generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule)
-
- elif generatedMethod == "TRT":
- collision_rule = create_lb_collision_rule(method='trt',stencil=get_stencil("D3Q19", 'walberla'), force_model=forcemodel, maxwellian_moments=False)
- generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule)
-
-
- elif generatedMethod == "KBC":
- collision_rule = create_lb_collision_rule(method='trt-kbc-n4',entropic=True, stencil=get_stencil("D3Q27", 'walberla'), compressible=True, force_model=forcemodel,
- optimization={"cse_global": True,
- "cse_pdfs": False,
- #"split": True
- })
- generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule)
-
- else:
- print("Invalid generated method string! " + generatedMethod)
+ forcemodel=Luo(forcing)
+
+ stencil = get_stencil("D3Q19", 'walberla')
+ omega = sp.symbols("omega_:%d" % len(stencil))
+
+ methodWithForce = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False,force_model=forcemodel)
+
+ def modification_func(moment, eq, rate):
+ omegaVisc = sp.Symbol("omega_visc")
+ omegaBulk = ps.fields("omega_bulk: [3D]", layout='fzyx')
+ omegaMagic = sp.Symbol("omega_magic")
+ if get_order(moment) <= 1:
+ return moment, eq, 0
+ elif rate == omega[1]:
+ return moment, eq, omegaBulk.center_vector
+ elif rate == omega[2]:
+ return moment, eq, omegaBulk.center_vector
+ elif is_even(moment):
+ return moment, eq, omegaVisc
+ else:
+ return moment, eq, omegaMagic
+
+ my_methodWithForce = create_lb_method_from_existing(methodWithForce, modification_func)
+
+ collision_rule = create_lb_collision_rule(lb_method=my_methodWithForce)
+ generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule)
diff --git a/apps/benchmarks/FluidParticleCoupling/LubricationForceEvaluation.cpp b/apps/benchmarks/FluidParticleCoupling/LubricationForceEvaluation.cpp
index 6182af2dcd78f41e8b5da294644560b7502e50e7..0be7fe660efa3adc5d4199bdb8f85e7fb9c00e0c 100644
--- a/apps/benchmarks/FluidParticleCoupling/LubricationForceEvaluation.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/LubricationForceEvaluation.cpp
@@ -84,11 +84,11 @@
#include "field/vtk/all.h"
#include "lbm/vtk/all.h"
-#include "GeneratedLBM.h"
-
#include <functional>
-#define USE_TRT_LIKE_LATTICE_MODEL
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
+#include "GeneratedLBM.h"
+#endif
namespace lubrication_force_evaluation
{
@@ -100,7 +100,11 @@ namespace lubrication_force_evaluation
using namespace walberla;
using walberla::uint_t;
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
using LatticeModel_T = lbm::GeneratedLBM;
+#else
+using LatticeModel_T = lbm::D3Q19< lbm::collision_model::D3Q19MRT>;
+#endif
using Stencil_T = LatticeModel_T::Stencil;
using PdfField_T = lbm::PdfField<LatticeModel_T>;
@@ -351,7 +355,7 @@ int main( int argc, char **argv )
uint_t randomSeed = uint_c(std::chrono::system_clock::now().time_since_epoch().count());
mpi::broadcastObject(randomSeed); // root process chooses seed and broadcasts it
- std::mt19937 randomNumberGenerator(randomSeed);
+ std::mt19937 randomNumberGenerator(static_cast<unsigned int>(randomSeed));
Vector3<real_t> domainCenter( real_c(xSize) * real_t(0.5), real_c(ySize) * real_t(0.5), real_c(zSize) * real_t(0.5) );
Vector3<real_t> offsetVector(math::realRandom<real_t>(real_t(0), real_t(1), randomNumberGenerator),
@@ -499,17 +503,14 @@ int main( int argc, char **argv )
BlockDataID omegaBulkFieldID = field::addToStorage<ScalarField_T>( blocks, "omega bulk field", omegaBulk, field::fzyx, uint_t(0) );
// create the lattice model
-#ifdef USE_TRT_LIKE_LATTICE_MODEL
- WALBERLA_LOG_INFO_ON_ROOT("Using TRTlike model!");
real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega );
real_t lambda_d = lbm::collision_model::TRT::lambda_d( omega, magicNumber );
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
+ WALBERLA_LOG_INFO_ON_ROOT("Using generated TRT-like lattice model!");
LatticeModel_T latticeModel = LatticeModel_T(omegaBulkFieldID, lambda_d, lambda_e);
- //LatticeModel_T latticeModel = LatticeModel_T(omegaBulk, lambda_d, lambda_e);
- //LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::D3Q19MRT( omegaBulk, omegaBulk, lambda_d, lambda_e, lambda_e, lambda_d ));
#else
- WALBERLA_LOG_INFO_ON_ROOT("Using KBC model!");
- // generated KBC
- LatticeModel_T latticeModel = LatticeModel_T(omega);
+ WALBERLA_LOG_INFO_ON_ROOT("Using waLBerla built-in MRT lattice model and ignoring omega bulk field since not supported!");
+ LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::D3Q19MRT( omegaBulk, omegaBulk, lambda_d, lambda_e, lambda_e, lambda_d ));
#endif
// add PDF field
@@ -621,7 +622,13 @@ int main( int argc, char **argv )
<< Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingID ), "Boundary Handling" );
// stream + collide LBM step
- timeloop.add() << Sweep( LatticeModel_T::Sweep( pdfFieldID ), "LB stream & collide" );
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
+ auto lbmSweep = LatticeModel_T::Sweep( pdfFieldID );
+ timeloop.add() << Sweep( lbmSweep, "LB sweep" );
+#else
+ auto lbmSweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag );
+ timeloop.add() << Sweep( makeSharedSweep( lbmSweep ), "cell-wise LB sweep" );
+#endif
////////////////////////
diff --git a/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp b/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
index c8601416f4e27d421b913ea52dd56dde3abee7e4..966a55925f85fea9a051687ae6836afbfc27ca7b 100644
--- a/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
@@ -97,11 +97,11 @@
#include "field/vtk/all.h"
#include "lbm/vtk/all.h"
-#include "GeneratedLBM.h"
-
#include <functional>
-#define USE_TRT_LIKE_LATTICE_MODEL
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
+#include "GeneratedLBM.h"
+#endif
namespace oblique_wet_collision
{
@@ -113,7 +113,11 @@ namespace oblique_wet_collision
using namespace walberla;
using walberla::uint_t;
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
using LatticeModel_T = lbm::GeneratedLBM;
+#else
+using LatticeModel_T = lbm::D3Q19< lbm::collision_model::D3Q19MRT>;
+#endif
using Stencil_T = LatticeModel_T::Stencil;
using PdfField_T = lbm::PdfField<LatticeModel_T>;
@@ -216,7 +220,7 @@ class SpherePropertyLogger
{
public:
SpherePropertyLogger( const shared_ptr< ParticleAccessor_T > & ac, walberla::id_t sphereUid,
- const std::string fileNameLogging, const std::string fileNameForceLogging, bool fileIO,
+ const std::string & fileNameLogging, const std::string & fileNameForceLogging, bool fileIO,
real_t diameter, real_t uIn, real_t impactRatio, uint_t numberOfSubSteps, real_t fGravX, real_t fGravZ) :
ac_( ac ), sphereUid_( sphereUid ), fileNameLogging_( fileNameLogging ), fileNameForceLogging_( fileNameForceLogging ),
fileIO_(fileIO),
@@ -599,7 +603,9 @@ int main( int argc, char **argv )
const real_t dx_SI = diameter_SI / diameter;
const real_t sphereVolume = math::pi / real_t(6) * diameter * diameter * diameter;
- real_t dt_SI, viscosity, omega;
+ real_t dt_SI;
+ real_t viscosity;
+ real_t omega;
real_t uIn = real_t(1);
real_t accelerationFactor = real_t(0);
bool applyArtificialGravityAfterAccelerating = false;
@@ -828,17 +834,14 @@ int main( int argc, char **argv )
//LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( real_t(1) / relaxationTime ) );
// generated MRT
-#ifdef USE_TRT_LIKE_LATTICE_MODEL
- WALBERLA_LOG_INFO_ON_ROOT("Using TRT-like lattice model!");
real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega );
real_t lambda_d = lbm::collision_model::TRT::lambda_d( omega, magicNumber );
- //LatticeModel_T latticeModel = LatticeModel_T(omegaBulk, lambda_d, lambda_e);
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
+ WALBERLA_LOG_INFO_ON_ROOT("Using generated TRT-like lattice model!");
LatticeModel_T latticeModel = LatticeModel_T(omegaBulkFieldID, lambda_d, lambda_e);
- //LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::D3Q19MRT( omegaBulk, omegaBulk, lambda_d, lambda_e, lambda_e, lambda_d ));
#else
- WALBERLA_LOG_INFO_ON_ROOT("Using different lattice model!");
- // generated model with a single omega
- LatticeModel_T latticeModel = LatticeModel_T(omega);
+ WALBERLA_LOG_INFO_ON_ROOT("Using waLBerla built-in MRT lattice model and ignoring omega bulk field since not supported!");
+ LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::D3Q19MRT( omegaBulk, omegaBulk, lambda_d, lambda_e, lambda_e, lambda_d ));
#endif
// add PDF field
@@ -931,13 +934,6 @@ int main( int argc, char **argv )
auto particleVtkWriter = vtk::createVTKOutput_PointData(particleVtkOutput, "Particles", vtkIOFreq, baseFolder, "simulation_step");
timeloop.addFuncBeforeTimeStep( vtk::writeFiles( particleVtkWriter ), "VTK (sphere data)" );
- // flag field (written only once in the first time step, ghost layers are also written)
- /*
- auto flagFieldVTK = vtk::createVTKOutput_BlockData( blocks, "flag_field", timesteps, FieldGhostLayers, false, baseFolder );
- flagFieldVTK->addCellDataWriter( make_shared< field::VTKWriter< FlagField_T > >( flagFieldID, "FlagField" ) );
- timeloop.addFuncBeforeTimeStep( vtk::writeFiles( flagFieldVTK ), "VTK (flag field data)" );
- */
-
// pdf field, as slice
auto pdfFieldVTK = vtk::createVTKOutput_BlockData( blocks, "fluid_field", vtkIOFreq, 0, false, baseFolder );
@@ -1022,10 +1018,13 @@ int main( int argc, char **argv )
<< Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingID ), "Boundary Handling" );
// stream + collide LBM step
- // generated LBM sweep
- timeloop.add() << Sweep( LatticeModel_T::Sweep( pdfFieldID ), "LB stream & collide" );
- // walberla sweeps:
- //timeloop.add() << Sweep( makeSharedSweep( lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag ) ), "cell-wise LB sweep" );
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
+ auto lbmSweep = LatticeModel_T::Sweep( pdfFieldID );
+ timeloop.add() << Sweep( lbmSweep, "LB sweep" );
+#else
+ auto lbmSweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag );
+ timeloop.add() << Sweep( makeSharedSweep( lbmSweep ), "cell-wise LB sweep" );
+#endif
// evaluation functionality
std::string loggingFileName( baseFolder + "/LoggingObliqueWetCollision.txt");
@@ -1048,8 +1047,9 @@ int main( int argc, char **argv )
uint_t numBounces = uint_t(0);
uint_t tImpact = uint_t(0);
- real_t curVel(real_t(0)), oldVel(real_t(0));
- real_t maxSettlingVel = real_t(0);
+ real_t curVel(real_t(0));
+ real_t oldVel(real_t(0));
+ real_t maxSettlingVel(real_t(0));
real_t minGapSize(real_t(0));
@@ -1194,7 +1194,7 @@ int main( int argc, char **argv )
if( artificiallyAccelerateSphere )
{
lbm_mesapd_coupling::RegularParticlesSelector sphereSelector;
- real_t newSphereVel = uIn * (exp(-accelerationFactor * real_t(i) / responseTime ) - real_t(1));
+ real_t newSphereVel = uIn * (std::exp(-accelerationFactor * real_t(i) / responseTime ) - real_t(1));
ps->forEachParticle(useOpenMP, sphereSelector, *accessor,
[newSphereVel,impactAngle](const size_t idx, ParticleAccessor_T& ac){
ac.setLinearVelocity(idx, Vector3<real_t>( -std::sin(impactAngle), real_t(0), std::cos(impactAngle) ) * newSphereVel);
diff --git a/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp b/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp
index 6fa7822310fd770d223ec3beeefd3279b7aca029..b8a2b2c1929518665c6dc6f84fbe2b6d78558afc 100644
--- a/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp
@@ -87,11 +87,11 @@
#include "field/vtk/all.h"
#include "lbm/vtk/all.h"
-#include "GeneratedLBM.h"
-
#include <functional>
-#define USE_TRT_LIKE_LATTICE_MODEL
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
+#include "GeneratedLBM.h"
+#endif
namespace settling_sphere_in_box
{
@@ -103,8 +103,11 @@ namespace settling_sphere_in_box
using namespace walberla;
using walberla::uint_t;
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
using LatticeModel_T = lbm::GeneratedLBM;
-//using LatticeModel_T = lbm::D3Q19< lbm::collision_model::D3Q19MRT>;
+#else
+using LatticeModel_T = lbm::D3Q19< lbm::collision_model::D3Q19MRT>;
+#endif
using Stencil_T = LatticeModel_T::Stencil;
using PdfField_T = lbm::PdfField<LatticeModel_T>;
@@ -426,7 +429,8 @@ int main( int argc, char **argv )
const real_t diameter_SI = real_t(15e-3);
const real_t densitySphere_SI = real_t(1120);
- real_t densityFluid_SI, dynamicViscosityFluid_SI;
+ real_t densityFluid_SI;
+ real_t dynamicViscosityFluid_SI;
real_t expectedSettlingVelocity_SI;
switch( fluidType )
{
@@ -600,17 +604,14 @@ int main( int argc, char **argv )
BlockDataID omegaBulkFieldID = field::addToStorage<ScalarField_T>( blocks, "omega bulk field", omegaBulk, field::fzyx, uint_t(0) );
// create the lattice model
-#ifdef USE_TRT_LIKE_LATTICE_MODEL
- WALBERLA_LOG_INFO_ON_ROOT("Using TRTlike model!");
real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega );
real_t lambda_d = lbm::collision_model::TRT::lambda_d( omega, magicNumber );
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
+ WALBERLA_LOG_INFO_ON_ROOT("Using generated TRT-like lattice model!");
LatticeModel_T latticeModel = LatticeModel_T(omegaBulkFieldID, lambda_d, lambda_e);
- //LatticeModel_T latticeModel = LatticeModel_T(omegaBulk, lambda_d, lambda_e);
- //LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::D3Q19MRT( omegaBulk, omegaBulk, lambda_d, lambda_e, lambda_e, lambda_d ));
#else
- WALBERLA_LOG_INFO_ON_ROOT("Using KBC model!");
- // generated KBC
- LatticeModel_T latticeModel = LatticeModel_T(omega);
+ WALBERLA_LOG_INFO_ON_ROOT("Using waLBerla built-in MRT lattice model and ignoring omega bulk field since not supported!");
+ LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::D3Q19MRT( omegaBulk, omegaBulk, lambda_d, lambda_e, lambda_e, lambda_d ));
#endif
@@ -759,32 +760,21 @@ int main( int argc, char **argv )
timeloop.add() << Sweep( makeSharedSweep(omegaBulkAdapter), "Omega Bulk Adapter");
}
- bool useStreamCollide = true;
+
+ // add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
auto bhSweep = BoundaryHandling_T::getBlockSweep( boundaryHandlingID );
+ timeloop.add() << BeforeFunction( optimizedPDFCommunicationScheme, "LBM Communication" )
+ << Sweep(bhSweep, "Boundary Handling" );
- // generated sweeps
+ // stream + collide LBM step
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
auto lbmSweep = LatticeModel_T::Sweep( pdfFieldID );
- if( useStreamCollide )
- {
- // add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
- timeloop.add() << BeforeFunction( optimizedPDFCommunicationScheme, "LBM Communication" )
- << Sweep(bhSweep, "Boundary Handling" );
-
- // stream + collide LBM step
- timeloop.add() << Sweep( lbmSweep, "LB sweep" );
- }
- else
- {
- // collide LBM step
- timeloop.add() << Sweep([&lbmSweep](IBlock * const block){lbmSweep.collide(block);}, "cell-wise collide LB sweep" );
-
- // add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
- timeloop.add() << BeforeFunction( optimizedPDFCommunicationScheme, "LBM Communication" )
- << Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingID ), "Boundary Handling" );
+ timeloop.add() << Sweep( lbmSweep, "LB sweep" );
+#else
+ auto lbmSweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag );
+ timeloop.add() << Sweep( makeSharedSweep( lbmSweep ), "cell-wise LB sweep" );
+#endif
- // stream LBM step
- timeloop.add() << Sweep([&lbmSweep](IBlock * const block){lbmSweep.stream(block);}, "cell-wise stream LB sweep" );
- }
// evaluation functionality
std::string loggingFileName( baseFolder + "/LoggingSettlingSphere_");
diff --git a/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp b/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp
index 9e1b866792a5e5088c5e6940f8b89aa9074ede81..644a887df3aeb44096f2f267934bb9729035a89d 100644
--- a/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp
@@ -91,11 +91,11 @@
#include "field/vtk/all.h"
#include "lbm/vtk/all.h"
-#include "GeneratedLBM.h"
-
#include <functional>
-#define USE_TRT_LIKE_LATTICE_MODEL
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
+#include "GeneratedLBM.h"
+#endif
namespace sphere_moving_with_prescribed_velocity
{
@@ -107,8 +107,11 @@ namespace sphere_moving_with_prescribed_velocity
using namespace walberla;
using walberla::uint_t;
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
using LatticeModel_T = lbm::GeneratedLBM;
-//using LatticeModel_T = lbm::D3Q19< lbm::collision_model::D3Q19MRT>;
+#else
+using LatticeModel_T = lbm::D3Q19< lbm::collision_model::D3Q19MRT>;
+#endif
using Stencil_T = LatticeModel_T::Stencil;
using PdfField_T = lbm::PdfField<LatticeModel_T>;
@@ -458,8 +461,6 @@ int main( int argc, char **argv )
filesystem::create_directory( tpath );
}
- //if(adaptionLayerSize > real_t(2)) WALBERLA_LOG_WARNING_ON_ROOT("Adaption layer size is larger than permitted by parallel setup!");
-
//////////////////////////
// NUMERICAL PARAMETERS //
//////////////////////////
@@ -557,15 +558,14 @@ int main( int argc, char **argv )
BlockDataID omegaBulkFieldID = field::addToStorage<ScalarField_T>( blocks, "omega bulk field", omegaBulk, field::fzyx, uint_t(0) );
// create the lattice model
-#ifdef USE_TRT_LIKE_LATTICE_MODEL
real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega );
real_t lambda_d = lbm::collision_model::TRT::lambda_d( omega, magicNumber );
- //LatticeModel_T latticeModel = LatticeModel_T(omegaBulk, lambda_d, lambda_e);
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
+ WALBERLA_LOG_INFO_ON_ROOT("Using generated TRT-like lattice model!");
LatticeModel_T latticeModel = LatticeModel_T(omegaBulkFieldID, lambda_d, lambda_e);
- //LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::D3Q19MRT( omegaBulk, omegaBulk, lambda_d, lambda_e, lambda_e, lambda_d ));
#else
- // generated KBC
- LatticeModel_T latticeModel = LatticeModel_T(omega);
+ WALBERLA_LOG_INFO_ON_ROOT("Using waLBerla built-in MRT lattice model and ignoring omega bulk field since not supported!");
+ LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::D3Q19MRT( omegaBulk, omegaBulk, lambda_d, lambda_e, lambda_e, lambda_d ));
#endif
// add PDF field
@@ -604,9 +604,6 @@ int main( int argc, char **argv )
mesa_pd::mpi::ReduceProperty reduceProperty;
// set up coupling functionality
- //Vector3<real_t> gravitationalForce( real_t(0), real_t(0), -(densityRatio - real_t(1)) * gravitationalAcceleration * sphereVolume );
- //lbm_mesapd_coupling::AddForceOnParticlesKernel addGravitationalForce(gravitationalForce);
- lbm_mesapd_coupling::AddHydrodynamicInteractionKernel addHydrodynamicInteraction;
lbm_mesapd_coupling::ResetHydrodynamicForceTorqueKernel resetHydrodynamicForceTorque;
lbm_mesapd_coupling::AverageHydrodynamicForceTorqueKernel averageHydrodynamicForceTorque;
@@ -653,11 +650,6 @@ int main( int argc, char **argv )
auto particleVtkWriter = vtk::createVTKOutput_PointData(particleVtkOutput, "Particles", vtkIOFreq, baseFolder, "simulation_step");
timeloop.addFuncBeforeTimeStep( vtk::writeFiles( particleVtkWriter ), "VTK (sphere data)" );
- // flag field (written only once in the first time step, ghost layers are also written)
- //auto flagFieldVTK = vtk::createVTKOutput_BlockData( blocks, "flag_field", timesteps, FieldGhostLayers, false, baseFolder );
- //flagFieldVTK->addCellDataWriter( make_shared< field::VTKWriter< FlagField_T > >( flagFieldID, "FlagField" ) );
- //timeloop.addFuncBeforeTimeStep( vtk::writeFiles( flagFieldVTK ), "VTK (flag field data)" );
-
// pdf field
auto pdfFieldVTK = vtk::createVTKOutput_BlockData( blocks, "fluid_field", vtkIOFreq, 0, false, baseFolder );
@@ -720,33 +712,19 @@ int main( int argc, char **argv )
timeloop.add() << Sweep( makeSharedSweep(omegaBulkAdapter), "Omega Bulk Adapter");
}
-
- bool useStreamCollide = true;
+ // add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
auto bhSweep = BoundaryHandling_T::getBlockSweep( boundaryHandlingID );
+ timeloop.add() << BeforeFunction( optimizedPDFCommunicationScheme, "LBM Communication" )
+ << Sweep(bhSweep, "Boundary Handling" );
- // generated sweeps
+ // stream + collide LBM step
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
auto lbmSweep = LatticeModel_T::Sweep( pdfFieldID );
- if( useStreamCollide )
- {
- // add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
- timeloop.add() << BeforeFunction( optimizedPDFCommunicationScheme, "LBM Communication" )
- << Sweep(bhSweep, "Boundary Handling" );
-
- // stream + collide LBM step
- timeloop.add() << Sweep( lbmSweep, "LB sweep" );
- }
- else
- {
- // collide LBM step
- timeloop.add() << Sweep([&lbmSweep](IBlock * const block){lbmSweep.collide(block);}, "cell-wise collide LB sweep" );
-
- // add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
- timeloop.add() << BeforeFunction( optimizedPDFCommunicationScheme, "LBM Communication" )
- << Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingID ), "Boundary Handling" );
-
- // stream LBM step
- timeloop.add() << Sweep([&lbmSweep](IBlock * const block){lbmSweep.stream(block);}, "cell-wise stream LB sweep" );
- }
+ timeloop.add() << Sweep( lbmSweep, "LB sweep" );
+#else
+ auto lbmSweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag );
+ timeloop.add() << Sweep( makeSharedSweep( lbmSweep ), "cell-wise LB sweep" );
+#endif
// evaluation functionality
std::string loggingFileName( baseFolder + "/LoggingPrescribedMovingSphere");
@@ -760,7 +738,6 @@ int main( int argc, char **argv )
if( conserveMomentum ) loggingFileName += "_conserveMomentum";
if( artificiallyAccelerateSphere ) loggingFileName += "_acc";
if( !fileNameEnding.empty()) loggingFileName += "_" + fileNameEnding;
- //loggingFileName += "_SBB";
loggingFileName += ".txt";
if( fileIO )
@@ -806,7 +783,7 @@ int main( int argc, char **argv )
if( artificiallyAccelerateSphere )
{
// accelerating after reading from a checkpoint file does not make sense, as current actual runtime is not known
- real_t newSphereVel = velocity * (exp(-accelerationFactor * real_t(i) ) - real_t(1));
+ real_t newSphereVel = velocity * (std::exp(-accelerationFactor * real_t(i) ) - real_t(1));
ps->forEachParticle(useOpenMP, sphereSelector, *accessor, [newSphereVel](const size_t idx, ParticleAccessor_T& ac){ ac.setLinearVelocity(idx, Vector3<real_t>(real_t(0), real_t(0), newSphereVel));}, *accessor);
}
diff --git a/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp b/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
index 61a3844f14a27fc0c6bfc89e953a419f6653ed11..8a339e6915f29dff6fcaea552b68923deb5540ce 100644
--- a/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
@@ -97,11 +97,11 @@
#include "field/vtk/all.h"
#include "lbm/vtk/all.h"
-#include "GeneratedLBM.h"
-
#include <functional>
-#define USE_TRT_LIKE_LATTICE_MODEL
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
+#include "GeneratedLBM.h"
+#endif
namespace sphere_wall_collision
{
@@ -113,7 +113,11 @@ namespace sphere_wall_collision
using namespace walberla;
using walberla::uint_t;
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
using LatticeModel_T = lbm::GeneratedLBM;
+#else
+using LatticeModel_T = lbm::D3Q19< lbm::collision_model::D3Q19MRT>;
+#endif
using Stencil_T = LatticeModel_T::Stencil;
using PdfField_T = lbm::PdfField<LatticeModel_T>;
@@ -292,7 +296,7 @@ public:
return position_;
}
- void writeResult(std::string filename, uint_t tImpact)
+ void writeResult(const std::string & filename, uint_t tImpact)
{
WALBERLA_ROOT_SECTION()
{
@@ -455,22 +459,16 @@ real_t getForce(walberla::id_t uid, ParticleAccessor_T & ac)
return force;
}
-
-real_t getAcceleratedSphereVelocity(real_t accelerationRate, uint_t timestep, real_t tref, real_t finalVelocity)
-{
- return finalVelocity * (std::exp(-real_c(timestep)/(tref*accelerationRate))-real_t(1));
-}
-
template< typename DensityInterpolator_T>
void logDensityToFile(const shared_ptr<StructuredBlockStorage> & blocks, BlockDataID densityInterpolatorID, real_t surfaceDistance,
- Vector3<real_t> spherePosition, real_t diameter, std::string fileName, uint_t timestep, real_t averageDensityInSystem )
+ Vector3<real_t> spherePosition, real_t diameter, const std::string & fileName, uint_t timestep, real_t averageDensityInSystem )
{
real_t evalRadius = real_t(0.5) * diameter + surfaceDistance;
std::vector< Vector3<real_t> > evaluationPositions;
evaluationPositions.push_back(spherePosition + Vector3<real_t>(real_t(0), real_t(0), evalRadius));
- evaluationPositions.push_back(spherePosition + Vector3<real_t>(evalRadius / sqrt(real_t(2)), real_t(0), evalRadius / sqrt(real_t(2))));
+ evaluationPositions.push_back(spherePosition + Vector3<real_t>(evalRadius / real_c(sqrt(real_t(2))), real_t(0), real_c(evalRadius / sqrt(real_t(2)))));
evaluationPositions.push_back(spherePosition + Vector3<real_t>(evalRadius, real_t(0), real_t(0)));
- evaluationPositions.push_back(spherePosition + Vector3<real_t>(evalRadius / sqrt(real_t(2)), real_t(0), -evalRadius / sqrt(real_t(2))));
+ evaluationPositions.push_back(spherePosition + Vector3<real_t>(evalRadius / real_c(sqrt(real_t(2))), real_t(0), -evalRadius / real_c(sqrt(real_t(2)))));
evaluationPositions.push_back(spherePosition + Vector3<real_t>(real_t(0), real_t(0), -evalRadius));
std::vector<real_t> densityAtPos(evaluationPositions.size(), real_t(0));
@@ -689,7 +687,7 @@ int main( int argc, char **argv )
{
uint_t seed1 = uint_c(std::chrono::system_clock::now().time_since_epoch().count());
mpi::broadcastObject(seed1); // root process chooses seed and broadcasts it
- std::mt19937 g1 (seed1);
+ std::mt19937 g1 (static_cast<unsigned int>(seed1));
real_t initialPositionOffset = real_t(0.5)*math::realRandom<real_t>(real_t(-1), real_t(1), g1);
initialSpherePosition += initialPositionOffset;
@@ -901,22 +899,14 @@ int main( int argc, char **argv )
BlockDataID omegaBulkFieldID = field::addToStorage<ScalarField_T>( blocks, "omega bulk field", omegaBulk, field::fzyx, uint_t(0) );
// create the lattice model
-
- // TRT
- //LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( real_t(1) / relaxationTime ) );
-
- // generated MRT
-#ifdef USE_TRT_LIKE_LATTICE_MODEL
- WALBERLA_LOG_INFO_ON_ROOT("Using TRT-like lattice model!");
real_t lambda_e = lbm::collision_model::TRT::lambda_e( omega );
real_t lambda_d = lbm::collision_model::TRT::lambda_d( omega, magicNumber );
- //LatticeModel_T latticeModel = LatticeModel_T(omegaBulk, lambda_d, lambda_e);
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
+ WALBERLA_LOG_INFO_ON_ROOT("Using generated TRT-like lattice model!");
LatticeModel_T latticeModel = LatticeModel_T(omegaBulkFieldID, lambda_d, lambda_e);
- //LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::D3Q19MRT( omegaBulk, omegaBulk, lambda_d, lambda_e, lambda_e, lambda_d ));
#else
- WALBERLA_LOG_INFO_ON_ROOT("Using different lattice model!");
- // generated model with a single omega
- LatticeModel_T latticeModel = LatticeModel_T(omega);
+ WALBERLA_LOG_INFO_ON_ROOT("Using waLBerla built-in MRT lattice model and ignoring omega bulk field since not supported!");
+ LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::D3Q19MRT( omegaBulk, omegaBulk, lambda_d, lambda_e, lambda_e, lambda_d ));
#endif
// add PDF field
@@ -1025,7 +1015,6 @@ int main( int argc, char **argv )
auto particleVtkOutput = make_shared<mesa_pd::vtk::ParticleVtkOutput>(ps);
particleVtkOutput->addOutput<mesa_pd::data::SelectParticleOwner>("owner");
particleVtkOutput->addOutput<mesa_pd::data::SelectParticleUid>("uid");
- particleVtkOutput->addOutput<mesa_pd::data::SelectParticleShapeID>("shapeID");
particleVtkOutput->addOutput<mesa_pd::data::SelectParticleLinearVelocity>("velocity");
particleVtkOutput->setParticleSelector( [sphereShape](const mesa_pd::data::ParticleStorage::iterator& pIt) {return pIt->getShapeID() == sphereShape;} ); //limit output to sphere
auto particleVtkWriter = vtk::createVTKOutput_PointData(particleVtkOutput, "Particles", vtkIOFreq, baseFolder, "simulation_step");
@@ -1089,10 +1078,13 @@ int main( int argc, char **argv )
<< Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingID ), "Boundary Handling" );
// stream + collide LBM step
- // generated LBM sweep
- timeloop.add() << Sweep( LatticeModel_T::Sweep( pdfFieldID ), "LB stream & collide" );
- // walberla sweeps:
- //timeloop.add() << Sweep( makeSharedSweep( lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag ) ), "cell-wise LB sweep" );
+#ifdef WALBERLA_BUILD_WITH_CODEGEN
+ auto lbmSweep = LatticeModel_T::Sweep( pdfFieldID );
+ timeloop.add() << Sweep( lbmSweep, "LB sweep" );
+#else
+ auto lbmSweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag );
+ timeloop.add() << Sweep( makeSharedSweep( lbmSweep ), "cell-wise LB sweep" );
+#endif
// this is carried out after the particle integration, it corrects the flag field and restores missing PDF information
// then, the checkpointing file can be written, as otherwise some cells are invalid and can not be recovered
@@ -1134,8 +1126,6 @@ int main( int argc, char **argv )
}
-
-
// evaluation functionality
std::string loggingFileName( baseFolder + "/LoggingSphereWallCollision.txt");
if( fileIO )
@@ -1162,8 +1152,9 @@ int main( int argc, char **argv )
uint_t tImpact = uint_t(0);
std::vector<uint_t> impactTimes;
- real_t curVel(real_t(0)), oldVel(real_t(0));
- real_t maxSettlingVel = real_t(0);
+ real_t curVel(real_t(0));
+ real_t oldVel(real_t(0));
+ real_t maxSettlingVel(real_t(0));
std::vector<real_t> maxSettlingVelBetweenBounces;
real_t minGapSize(real_t(0));
@@ -1328,7 +1319,7 @@ int main( int argc, char **argv )
{
// overwrite velocity of particle with prescribed one
lbm_mesapd_coupling::RegularParticlesSelector sphereSelector;
- real_t newSphereVel = uIn * (exp(-accelerationFactor * real_t(i) ) - real_t(1));
+ real_t newSphereVel = uIn * (std::exp(-accelerationFactor * real_t(i) ) - real_t(1));
ps->forEachParticle(useOpenMP, sphereSelector, *accessor, [newSphereVel](const size_t idx, ParticleAccessor_T& ac){ ac.setLinearVelocity(idx, Vector3<real_t>(real_t(0), real_t(0), newSphereVel));}, *accessor);
}