diff --git a/apps/showcases/CMakeLists.txt b/apps/showcases/CMakeLists.txt
index 6330a83bdde5027a2d0fd8e030862b7bf521dcd4..64becac3b5b95193aa8ada372255ad92e9915b0b 100644
--- a/apps/showcases/CMakeLists.txt
+++ b/apps/showcases/CMakeLists.txt
@@ -8,6 +8,7 @@ add_subdirectory( ParticlePacking )
add_subdirectory( PegIntoSphereBed )
if ( WALBERLA_BUILD_WITH_CODEGEN AND WALBERLA_BUILD_WITH_PYTHON )
add_subdirectory( PhaseFieldAllenCahn )
+add_subdirectory( LDC )
endif()
if ( WALBERLA_BUILD_WITH_CODEGEN AND NOT WALBERLA_BUILD_WITH_OPENMP)
add_subdirectory( PorousMedia )
diff --git a/apps/showcases/LDC/CMakeLists.txt b/apps/showcases/LDC/CMakeLists.txt
index 339f648197b4d716c9357f19b5f03ad4fdb43ffd..c254106e37d2cec1c65e08b332219ea2fbed1ca7 100644
--- a/apps/showcases/LDC/CMakeLists.txt
+++ b/apps/showcases/LDC/CMakeLists.txt
@@ -1,51 +1,17 @@
# Code Generation Tutorials
-
+waLBerla_link_files_to_builddir(*.prm)
if( WALBERLA_BUILD_WITH_CODEGEN )
- # Tutorial 1: Heat Equation with pystencils
- walberla_generate_target_from_python( NAME CodegenHeatEquationKernel
- FILE HeatEquationKernel.py
- OUT_FILES HeatEquationKernel.cpp HeatEquationKernel.h )
-
- walberla_add_executable ( NAME 01_CodegenHeatEquation
- FILES 01_CodegenHeatEquation.cpp
- DEPENDS blockforest core field stencil timeloop vtk pde CodegenHeatEquationKernel )
-
- # Tutorial 2: lbmpy Lattice Model Generation
- waLBerla_link_files_to_builddir( *.prm )
-
- walberla_generate_target_from_python( NAME 02_LBMLatticeModelGenerationPython
- FILE 02_LBMLatticeModelGeneration.py
- OUT_FILES SRTLatticeModel.cpp SRTLatticeModel.h
- SRTPackInfo.cpp SRTPackInfo.h )
-
- walberla_add_executable ( NAME 02_LBMLatticeModelGenerationApp
- FILES 02_LBMLatticeModelGeneration.cpp
- DEPENDS blockforest core domain_decomposition field geometry timeloop lbm stencil vtk 02_LBMLatticeModelGenerationPython )
-
- # Tutorial 3: Advanced lbmpy Code Generation
- if(WALBERLA_BUILD_WITH_CUDA)
- walberla_generate_target_from_python( NAME 03_AdvancedLBMCodegenPython
- FILE 03_AdvancedLBMCodegen.py
- OUT_FILES CumulantMRTSweep.cu CumulantMRTSweep.h
- CumulantMRTPackInfo.cu CumulantMRTPackInfo.h
- InitialPDFsSetter.cu InitialPDFsSetter.h
- CumulantMRTNoSlip.cu CumulantMRTNoSlip.h)
-
- walberla_add_executable ( NAME 03_AdvancedLBMCodegenApp
- FILES 03_AdvancedLBMCodegen.cpp
- DEPENDS blockforest cuda core domain_decomposition field geometry timeloop lbm stencil vtk 03_AdvancedLBMCodegenPython )
- else()
- walberla_generate_target_from_python( NAME 03_AdvancedLBMCodegenPython
- FILE 03_AdvancedLBMCodegen.py
- OUT_FILES CumulantMRTSweep.cpp CumulantMRTSweep.h
- CumulantMRTPackInfo.cpp CumulantMRTPackInfo.h
- InitialPDFsSetter.cpp InitialPDFsSetter.h
- CumulantMRTNoSlip.cpp CumulantMRTNoSlip.h)
-
- walberla_add_executable ( NAME 03_AdvancedLBMCodegenApp
- FILES 03_AdvancedLBMCodegen.cpp
- DEPENDS blockforest core domain_decomposition field geometry timeloop lbm stencil vtk 03_AdvancedLBMCodegenPython )
- endif()
+ walberla_generate_target_from_python( NAME LDCCodegen
+ FILE LDC.py
+ OUT_FILES LDCSweep.cpp LDCSweep.h
+ LDCPackInfo.cpp LDCPackInfo.h
+ LDCSetter.cpp LDCSetter.h
+ LDCNoSlip.cpp LDCNoSlip.h
+ LDCUBB.cpp LDCUBB.h)
+
+ walberla_add_executable ( NAME LDC
+ FILES LDC.cpp
+ DEPENDS blockforest core domain_decomposition field geometry timeloop lbm stencil vtk LDCCodegen )
endif()
\ No newline at end of file
diff --git a/apps/showcases/LDC/LDC.cpp b/apps/showcases/LDC/LDC.cpp
index 492d40b8bf07458e2b8b847d36dff28e9356d72a..c76ea38518371cd2b112af51b0a565014839f539 100644
--- a/apps/showcases/LDC/LDC.cpp
+++ b/apps/showcases/LDC/LDC.cpp
@@ -19,34 +19,20 @@
//======================================================================================================================
#include "blockforest/all.h"
-
#include "core/all.h"
-
-#if defined(WALBERLA_BUILD_WITH_CUDA)
-# include "cuda/AddGPUFieldToStorage.h"
-# include "cuda/DeviceSelectMPI.h"
-# include "cuda/HostFieldAllocator.h"
-# include "cuda/ParallelStreams.h"
-# include "cuda/communication/GPUPackInfo.h"
-# include "cuda/communication/UniformGPUScheme.h"
-#endif
-
#include "domain_decomposition/all.h"
-
#include "field/all.h"
#include "field/vtk/VTKWriter.h"
-
#include "geometry/all.h"
-
-#include "stencil/D2Q9.h"
-
+#include "stencil/D3Q19.h"
#include "timeloop/all.h"
// Codegen Includes
-#include "CumulantMRTNoSlip.h"
-#include "CumulantMRTPackInfo.h"
-#include "CumulantMRTSweep.h"
-#include "InitialPDFsSetter.h"
+#include "LDCSweep.h"
+#include "LDCPackInfo.h"
+#include "LDCSetter.h"
+#include "LDCNoSlip.h"
+#include "LDCUBB.h"
namespace walberla
{
@@ -55,57 +41,23 @@ namespace walberla
///////////////////////
// Communication Pack Info
-typedef pystencils::CumulantMRTPackInfo PackInfo_T;
+typedef pystencils::LDCPackInfo PackInfo_T;
// LB Method Stencil
-typedef stencil::D2Q9 Stencil_T;
+typedef stencil::D3Q19 Stencil_T;
// PDF field type
typedef field::GhostLayerField< real_t, Stencil_T::Size > PdfField_T;
// Velocity Field Type
typedef field::GhostLayerField< real_t, Stencil_T::D > VectorField_T;
+typedef field::GhostLayerField< real_t, 1 > ScalarField_T;
// Boundary Handling
typedef walberla::uint8_t flag_t;
typedef FlagField< flag_t > FlagField_T;
-typedef lbm::CumulantMRTNoSlip NoSlip_T;
-
-#if defined(WALBERLA_BUILD_WITH_CUDA)
-typedef cuda::GPUField< real_t > GPUField;
-#endif
-
-//////////////////////////////////////////
-/// Shear Flow Velocity Initialization ///
-//////////////////////////////////////////
-
-void initShearFlowVelocityField(const shared_ptr< StructuredBlockForest >& blocks, const BlockDataID& velocityFieldId,
- const Config::BlockHandle& config)
-{
- math::RealRandom< real_t > rng(config.getParameter< std::mt19937::result_type >("noiseSeed", 42));
-
- real_t velocityMagnitude = config.getParameter< real_t >("velocityMagnitude", real_c(0.08));
- real_t noiseMagnitude = config.getParameter< real_t >("noiseMagnitude", real_c(0.1) * velocityMagnitude);
-
- real_t n_y = real_c(blocks->getNumberOfYCells());
-
- for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
- {
- auto u = (*blockIt).getData< VectorField_T >(velocityFieldId);
-
- for (auto cellIt = u->beginWithGhostLayerXYZ(); cellIt != u->end(); ++cellIt)
- {
- Cell globalCell(cellIt.cell());
- blocks->transformBlockLocalToGlobalCell(globalCell, *blockIt);
-
- real_t relative_y = real_c(globalCell.y()) / n_y;
-
- u->get(cellIt.cell(), 0) = relative_y < 0.3 || relative_y > 0.7 ? velocityMagnitude : -velocityMagnitude;
-
- u->get(cellIt.cell(), 1) = noiseMagnitude * rng();
- }
- }
-}
+typedef lbm::LDCNoSlip NoSlip_T;
+typedef lbm::LDCUBB UBB_T;
/////////////////////
/// Main Function ///
@@ -115,7 +67,7 @@ int main(int argc, char** argv)
{
walberla::Environment walberlaEnv(argc, argv);
- if (!walberlaEnv.config()) { WALBERLA_ABORT("No configuration file specified!"); }
+ if (!walberlaEnv.config()) { WALBERLA_ABORT("No configuration file specified!") }
///////////////////////////////////////////////////////
/// Block Storage Creation and Simulation Parameter ///
@@ -137,35 +89,45 @@ int main(int argc, char** argv)
////////////////////////////////////
// Common Fields
+ BlockDataID densityFieldId = field::addToStorage< ScalarField_T >(blocks, "density", real_c(1.0), field::fzyx);
BlockDataID velocityFieldId = field::addToStorage< VectorField_T >(blocks, "velocity", real_c(0.0), field::fzyx);
BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
-#if defined(WALBERLA_BUILD_WITH_CUDA)
- // GPU Field for PDFs
- BlockDataID pdfFieldId = cuda::addGPUFieldToStorage< cuda::GPUField< real_t > >(
- blocks, "pdf field on GPU", Stencil_T::Size, field::fzyx, uint_t(1));
+ for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt) {
+ auto densityField = blockIt->getData<ScalarField_T>(densityFieldId);
+ std::ifstream infile("density.txt");
+
+ WALBERLA_FOR_ALL_CELLS_XYZ(densityField,
+ real_t dens;
+ infile >> dens;
+ densityField->get(x, y, z) = dens;
+ )
+ infile.close();
+ }
+
+ for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt) {
+ auto velocityField = blockIt->getData<VectorField_T>(velocityFieldId);
+ std::ifstream infile("velocity.txt");
+
+ WALBERLA_FOR_ALL_CELLS_XYZ(velocityField,
+ for (uint_t i = 0; i < 3; ++i)
+ {
+ real_t velo;
+ infile >> velo;
+ velocityField->get(x, y, z, i) = velo;
+ }
+ )
+ infile.close();
+ }
+
+ WALBERLA_LOG_INFO_ON_ROOT("Finished parameter reading")
- // GPU Velocity Field
- BlockDataID velocityFieldIdGPU =
- cuda::addGPUFieldToStorage< VectorField_T >(blocks, velocityFieldId, "velocity on GPU", true);
-#else
// CPU Field for PDFs
BlockDataID pdfFieldId = field::addToStorage< PdfField_T >(blocks, "pdf field", real_c(0.0), field::fzyx);
-#endif
-
- // Velocity field setup
- auto shearFlowSetup = walberlaEnv.config()->getOneBlock("ShearFlowSetup");
- initShearFlowVelocityField(blocks, velocityFieldId, shearFlowSetup);
- real_t rho = shearFlowSetup.getParameter("rho", real_c(1.0));
// pdfs setup
-#if defined(WALBERLA_BUILD_WITH_CUDA)
- cuda::fieldCpy< GPUField, VectorField_T >(blocks, velocityFieldIdGPU, velocityFieldId);
- pystencils::InitialPDFsSetter pdfSetter(pdfFieldId, velocityFieldIdGPU, rho);
-#else
- pystencils::InitialPDFsSetter pdfSetter(pdfFieldId, velocityFieldId, rho);
-#endif
+ pystencils::LDCSetter pdfSetter(densityFieldId, pdfFieldId, velocityFieldId);
for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
{
@@ -176,11 +138,7 @@ int main(int argc, char** argv)
/// Sweep ///
/////////////
-#if defined(WALBERLA_BUILD_WITH_CUDA)
- pystencils::CumulantMRTSweep CumulantMRTSweep(pdfFieldId, velocityFieldIdGPU, omega);
-#else
- pystencils::CumulantMRTSweep CumulantMRTSweep(pdfFieldId, velocityFieldId, omega);
-#endif
+ pystencils::LDCSweep LBMKernel(densityFieldId, pdfFieldId, velocityFieldId, omega);
/////////////////////////
/// Boundary Handling ///
@@ -191,10 +149,12 @@ int main(int argc, char** argv)
auto boundariesConfig = walberlaEnv.config()->getOneBlock("Boundaries");
NoSlip_T noSlip(blocks, pdfFieldId);
+ UBB_T ubb(blocks, pdfFieldId);
geometry::initBoundaryHandling< FlagField_T >(*blocks, flagFieldId, boundariesConfig);
geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldId, fluidFlagUID);
+ ubb.fillFromFlagField< FlagField_T >(blocks, flagFieldId, FlagUID("UBB"), fluidFlagUID);
noSlip.fillFromFlagField< FlagField_T >(blocks, flagFieldId, FlagUID("NoSlip"), fluidFlagUID);
/////////////////
@@ -204,18 +164,13 @@ int main(int argc, char** argv)
SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);
// Communication
-#if defined(WALBERLA_BUILD_WITH_CUDA)
- cuda::communication::UniformGPUScheme< Stencil_T > com(blocks, 0);
- com.addPackInfo(make_shared< PackInfo_T >(pdfFieldId));
- auto communication = std::function< void() >([&]() { com.communicate(nullptr); });
-#else
blockforest::communication::UniformBufferedScheme< Stencil_T > communication(blocks);
communication.addPackInfo(make_shared< PackInfo_T >(pdfFieldId));
-#endif
// Timeloop
timeloop.add() << BeforeFunction(communication, "communication") << Sweep(noSlip);
- timeloop.add() << Sweep(CumulantMRTSweep);
+ timeloop.add() << Sweep(ubb);
+ timeloop.add() << Sweep(LBMKernel);
// Time logger
timeloop.addFuncAfterTimeStep(timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency),
@@ -223,16 +178,12 @@ int main(int argc, char** argv)
if (VTKwriteFrequency > 0)
{
- const std::string path = "vtk_out/tut03";
- auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "cumulant_mrt_velocity_field", VTKwriteFrequency, 0,
+ const std::string path = "vtk_out";
+ auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "LDC", VTKwriteFrequency, 0,
false, path, "simulation_step", false, true, true, false, 0);
-#if defined(WALBERLA_BUILD_WITH_CUDA)
- // Copy velocity data to CPU before output
- vtkOutput->addBeforeFunction(
- [&]() { cuda::fieldCpy< VectorField_T, GPUField >(blocks, velocityFieldId, velocityFieldIdGPU); });
-#endif
-
+ auto densWriter = make_shared< field::VTKWriter< ScalarField_T > >(densityFieldId, "Density");
+ vtkOutput->addCellDataWriter(densWriter);
auto velWriter = make_shared< field::VTKWriter< VectorField_T > >(velocityFieldId, "Velocity");
vtkOutput->addCellDataWriter(velWriter);
@@ -246,4 +197,7 @@ int main(int argc, char** argv)
} // namespace walberla
-int main(int argc, char** argv) { return walberla::main(argc, argv); }
+int main(int argc, char** argv)
+{
+ return walberla::main(argc, argv);
+}
diff --git a/apps/showcases/LDC/LDC.prm b/apps/showcases/LDC/LDC.prm
index d292c02d13450a5ae947f3284ded79baf2348538..46debbcba77c76c0d5597ebc4ee461ed9ce254df 100644
--- a/apps/showcases/LDC/LDC.prm
+++ b/apps/showcases/LDC/LDC.prm
@@ -1,28 +1,18 @@
Parameters
{
- omega 1.8;
- timesteps 10001;
+ omega 1.81818;
+ timesteps 201;
remainingTimeLoggerFrequency 3; // in seconds
- VTKwriteFrequency 1000;
-}
-
-ShearFlowSetup
-{
- rho 1.0;
-
- velocityMagnitude 0.08;
- noiseMagnitude 0.005;
-
- noiseSeed 42;
+ VTKwriteFrequency 100;
}
DomainSetup
{
blocks < 1, 1, 1 >;
- cellsPerBlock < 300, 80, 1 >;
- periodic < 1, 0, 1 >;
+ cellsPerBlock < 128, 128, 128 >;
+ periodic < 0, 0, 0 >;
}
StabilityChecker
@@ -33,6 +23,12 @@ StabilityChecker
}
Boundaries
-{
- Border { direction S,N; walldistance -1; flag NoSlip; }
+{
+
+ Border { direction E; walldistance -1; flag NoSlip; }
+ Border { direction W; walldistance -1; flag NoSlip; }
+ Border { direction S; walldistance -1; flag NoSlip; }
+ Border { direction T; walldistance -1; flag UBB; }
+ Border { direction N; walldistance -1; flag NoSlip; }
+ Border { direction B; walldistance -1; flag NoSlip; }
}
diff --git a/apps/showcases/LDC/LDC.py b/apps/showcases/LDC/LDC.py
index c94eb81bab1f584ccfaf2ec4b9ebf2545c4434f9..e7b99a3035f56a06f56fd0055c68c25fc4d353af 100644
--- a/apps/showcases/LDC/LDC.py
+++ b/apps/showcases/LDC/LDC.py
@@ -5,7 +5,7 @@ from lbmpy import LBMConfig, LBMOptimisation, LBStencil, Method, Stencil
from lbmpy.creationfunctions import create_lb_update_rule
from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
-from lbmpy.boundaries import NoSlip
+from lbmpy.boundaries import NoSlip, UBB
from pystencils_walberla import CodeGeneration, generate_sweep, generate_pack_info_from_kernel
from lbmpy_walberla import generate_boundary
@@ -15,16 +15,17 @@ from lbmpy_walberla import generate_boundary
# General Parameters
# ========================
-stencil = LBStencil(Stencil.D2Q9)
+stencil = LBStencil(Stencil.D3Q19)
omega = sp.Symbol('omega')
layout = 'fzyx'
# PDF Fields
-pdfs, pdfs_tmp = ps.fields(f'pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): [2D]', layout=layout)
+pdfs, pdfs_tmp = ps.fields(f'pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): [{stencil.D}D]', layout=layout)
# Velocity Output Field
-velocity = ps.fields(f"velocity({stencil.D}): [2D]", layout=layout)
-output = {'velocity': velocity}
+velocity = ps.fields(f"velocity({stencil.D}): [{stencil.D}D]", layout=layout)
+density = ps.fields(f"density: [{stencil.D}D]", layout=layout)
+output = {'density': density, 'velocity': velocity}
# LBM Optimisation
lbm_opt = LBMOptimisation(cse_global=True,
@@ -38,9 +39,9 @@ lbm_opt = LBMOptimisation(cse_global=True,
# ==================
lbm_config = LBMConfig(stencil=stencil,
- method=Method.CUMULANT,
+ method=Method.TRT,
relaxation_rate=omega,
- compressible=True,
+ compressible=False,
output=output)
lbm_update_rule = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
@@ -51,10 +52,8 @@ lbm_method = lbm_update_rule.method
# PDF Initialization
# ========================
-initial_rho = sp.Symbol('rho_0')
-
pdfs_setter = macroscopic_values_setter(lbm_method,
- initial_rho,
+ density,
velocity.center_vector,
pdfs.center_vector)
@@ -63,19 +62,17 @@ pdfs_setter = macroscopic_values_setter(lbm_method,
# =====================
with CodeGeneration() as ctx:
- if ctx.cuda:
- target = ps.Target.GPU
- else:
- target = ps.Target.CPU
+ target = ps.Target.CPU
# LBM Sweep
- generate_sweep(ctx, "CumulantMRTSweep", lbm_update_rule, field_swaps=[(pdfs, pdfs_tmp)], target=target)
+ generate_sweep(ctx, "LDCSweep", lbm_update_rule, field_swaps=[(pdfs, pdfs_tmp)], target=target)
# Pack Info
- generate_pack_info_from_kernel(ctx, "CumulantMRTPackInfo", lbm_update_rule, target=target)
+ generate_pack_info_from_kernel(ctx, "LDCPackInfo", lbm_update_rule, target=target)
# Macroscopic Values Setter
- generate_sweep(ctx, "InitialPDFsSetter", pdfs_setter, target=target)
+ generate_sweep(ctx, "LDCSetter", pdfs_setter, target=target)
# NoSlip Boundary
- generate_boundary(ctx, "CumulantMRTNoSlip", NoSlip(), lbm_method, target=target)
+ generate_boundary(ctx, "LDCNoSlip", NoSlip(), lbm_method, target=target)
+ generate_boundary(ctx, "LDCUBB", UBB((0.00014587, 0, 0)), lbm_method, target=target)