From 2f7b33ed781f03fdecf2d7bcef09ef88f6466e92 Mon Sep 17 00:00:00 2001
From: Markus Holzer <markus.holzer@fau.de>
Date: Tue, 24 May 2022 16:08:47 +0200
Subject: [PATCH] Adaption for CodeGen pipeline
---
.gitignore | 4 +-
.gitlab-ci.yml | 121 +++++++++-----
.../FlowAroundSphereCodeGen.cpp | 24 +--
.../FlowAroundSphereCodeGen.py | 7 +-
.../FluidParticleCoupling/GeneratedLBM.py | 47 +++---
.../GeneratedLBMWithForce.py | 50 +++---
.../benchmark_multiphase.cpp | 18 +--
.../PhaseFieldAllenCahn/multiphase_codegen.py | 117 +++++++-------
.../UniformGridCPU/UniformGridCPU.cpp | 8 +-
.../UniformGridCPU/UniformGridCPU.py | 3 +-
.../UniformGridGPU/UniformGridGPU.cpp | 8 +-
.../UniformGridGPU/UniformGridGPU.py | 5 +-
.../old_ideas/UniformGridGPU.py | 2 +-
.../CPU/InitializerFunctions.cpp | 142 ++++++++---------
.../PhaseFieldAllenCahn/CPU/multiphase.cpp | 14 +-
.../CPU/multiphase_codegen.py | 125 ++++++---------
.../GPU/InitializerFunctions.cpp | 148 +++++++++---------
.../PhaseFieldAllenCahn/GPU/multiphase.cpp | 18 +--
.../GPU/multiphase_codegen.py | 140 +++++++----------
.../PhaseFieldAllenCahn/GPU/util.cpp | 30 ++--
.../PorousMediaCumulantLBMKernel.py | 50 +++---
.../codegen/02_LBMLatticeModelGeneration.cpp | 4 +-
.../codegen/02_LBMLatticeModelGeneration.py | 3 +-
.../codegen/03_AdvancedLBMCodegen.cpp | 6 +-
.../codegen/03_AdvancedLBMCodegen.py | 93 ++++++-----
apps/tutorials/codegen/HeatEquationKernel.py | 30 ++--
.../lbmpy_walberla/additional_data_handler.py | 2 +-
python/lbmpy_walberla/packinfo.py | 4 +-
python/lbmpy_walberla/sparse.py | 4 +-
.../templates/LatticeModel.tmpl.cpp | 4 +-
.../templates/LatticeModel.tmpl.h | 135 ++++++++--------
.../tests/test_walberla_codegen.py | 3 +-
.../lbmpy_walberla/walberla_lbm_generation.py | 103 +++++++-----
python/pystencils_walberla/boundary.py | 12 +-
python/pystencils_walberla/codegen.py | 9 +-
python/pystencils_walberla/jinja_filters.py | 28 +++-
.../pystencils_walberla/kernel_selection.py | 2 +-
.../templates/Boundary.tmpl.cpp | 2 +-
.../templates/Boundary.tmpl.h | 6 +-
python/pystencils_walberla/utility.py | 2 +-
python/waLBerla_tests/test_blockforest.py | 9 +-
python/waLBerla_tests/test_field.py | 10 +-
src/core/DataTypes.h | 2 +-
src/python_coupling/export/BasicExport.cpp | 5 +
tests/cuda/SimpleKernelTest.cpp | 38 ++---
tests/cuda/codegen/CodegenJacobiGPU.cpp | 20 +--
tests/cuda/codegen/CodegenPoissonGPU.cpp | 14 +-
tests/cuda/codegen/CudaJacobiKernel.py | 5 +-
tests/cuda/codegen/CudaPoisson.py | 3 +-
.../codegen/GeneratedFieldPackInfoTestGPU.cpp | 2 +-
.../codegen/GeneratedFieldPackInfoTestGPU.py | 6 +-
tests/cuda/codegen/MicroBenchmarkGpuLbm.cpp | 8 +-
tests/cuda/codegen/MicroBenchmarkGpuLbm.py | 10 +-
tests/field/codegen/CodegenJacobiCPU.cpp | 20 +--
tests/field/codegen/CodegenPoissonCPU.cpp | 16 +-
tests/field/codegen/EK.cpp | 20 +--
tests/field/codegen/EK.py | 2 +
.../codegen/GeneratedFieldPackInfoTest.cpp | 66 ++++----
.../codegen/GeneratedFieldPackInfoTest.py | 6 +-
tests/field/codegen/JacobiKernel.py | 5 +-
tests/field/codegen/MultipleFieldSwaps.cpp | 16 +-
tests/field/codegen/MultipleFieldSwaps.py | 7 +-
tests/field/codegen/Poisson.py | 3 +-
.../FieldLayoutAndVectorizationTest.cpp | 14 +-
.../FieldLayoutAndVectorizationTest.py | 3 +-
tests/lbm/codegen/FluctuatingMRT.cpp | 18 +--
tests/lbm/codegen/FluctuatingMRT.py | 28 ++--
tests/lbm/codegen/GeneratedFreeSlip.cpp | 13 +-
tests/lbm/codegen/GeneratedFreeSlip.py | 65 ++++----
tests/lbm/codegen/GeneratedOutflowBC.cpp | 16 +-
tests/lbm/codegen/GeneratedOutflowBC.py | 57 +++----
tests/lbm/codegen/InplaceStreamingCodegen.py | 68 ++++----
.../codegen/InplaceStreamingCodegenTest.cpp | 3 +-
tests/lbm/codegen/LbCodeGenerationExample.cpp | 10 +-
tests/lbm/codegen/LbCodeGenerationExample.py | 7 +-
.../lbm/codegen/LbmPackInfoGenerationTest.py | 2 +-
.../codegen/StreamInCellIntervalCodegen.py | 15 +-
.../StreamInCellIntervalCodegenTest.cpp | 8 +-
78 files changed, 1115 insertions(+), 1038 deletions(-)
diff --git a/.gitignore b/.gitignore
index 4430e54af..c755c9f0b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -72,6 +72,4 @@ cmake_install.cmake
CMakeDefs.h
/moduleStatistics.json
/walberla-config.cmake
-/cmake-build-debug/
-/cmake-build-release/
-/cmake-build-debug-remote/
+cmake-build-*
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index c47ea2f9c..af80860b1 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -503,7 +503,7 @@ gcc_9_serial:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:9
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -526,7 +526,7 @@ gcc_9_mpionly:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:9
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -547,7 +547,7 @@ gcc_9_hybrid:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:9
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -567,7 +567,7 @@ gcc_9_serial_dbg:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:9
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -591,7 +591,7 @@ gcc_9_mpionly_dbg:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:9
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -613,7 +613,7 @@ gcc_9_hybrid_dbg:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:9
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -633,12 +633,20 @@ gcc_9_hybrid_dbg:
gcc_9_hybrid_dbg_sp:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:9
+ before_script:
+ - pip3 install lbmpy==1.0 jinja2 pytest
+ - cd python
+ - python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
+ - cd ..
+ - CC=gcc CXX=g++ pip3 install pycuda
variables:
WALBERLA_BUILD_WITH_CUDA: "ON"
CMAKE_BUILD_TYPE: "DebugOptimized"
WALBERLA_DOUBLE_ACCURACY: "OFF"
WALBERLA_BUILD_WITH_PARMETIS: "OFF"
WALBERLA_BUILD_WITH_METIS: "OFF"
+ WALBERLA_BUILD_WITH_CODEGEN: "ON"
+ WALBERLA_BUILD_WITH_PYTHON: "ON"
only:
variables:
- $ENABLE_NIGHTLY_BUILDS
@@ -650,7 +658,7 @@ gcc_10_serial:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:10
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -673,7 +681,7 @@ gcc_10_mpionly:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:10
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -694,7 +702,7 @@ gcc_10_hybrid:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:10
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -714,7 +722,7 @@ gcc_10_serial_dbg:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:10
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -738,7 +746,7 @@ gcc_10_mpionly_dbg:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:10
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -760,7 +768,7 @@ gcc_10_hybrid_dbg:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:10
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -780,12 +788,20 @@ gcc_10_hybrid_dbg:
gcc_10_hybrid_dbg_sp:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:10
+ before_script:
+ - pip3 install lbmpy==1.0 jinja2 pytest
+ - cd python
+ - python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
+ - cd ..
+ - CC=gcc CXX=g++ pip3 install pycuda
variables:
WALBERLA_BUILD_WITH_CUDA: "ON"
CMAKE_BUILD_TYPE: "DebugOptimized"
WALBERLA_DOUBLE_ACCURACY: "OFF"
WALBERLA_BUILD_WITH_PARMETIS: "OFF"
WALBERLA_BUILD_WITH_METIS: "OFF"
+ WALBERLA_BUILD_WITH_CODEGEN: "ON"
+ WALBERLA_BUILD_WITH_PYTHON: "ON"
only:
variables:
- $ENABLE_NIGHTLY_BUILDS
@@ -797,7 +813,7 @@ gcc_11_serial:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:11
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -820,7 +836,7 @@ gcc_11_mpionly:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:11
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -842,11 +858,12 @@ gcc_11_hybrid:
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:11
stage: pretest
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
- CC=gcc CXX=g++ pip3 install pycuda
+ - pip3 list
variables:
WALBERLA_BUILD_WITH_CUDA: "ON"
WALBERLA_BUILD_WITH_CODEGEN: "ON"
@@ -859,7 +876,7 @@ gcc_11_serial_dbg:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:11
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -880,7 +897,7 @@ gcc_11_mpionly_dbg:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:11
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -899,7 +916,7 @@ gcc_11_hybrid_dbg:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:11
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -916,12 +933,20 @@ gcc_11_hybrid_dbg:
gcc_11_hybrid_dbg_sp:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/gcc:11
+ before_script:
+ - pip3 install lbmpy==1.0 jinja2 pytest
+ - cd python
+ - python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
+ - cd ..
+ - CC=gcc CXX=g++ pip3 install pycuda
variables:
WALBERLA_BUILD_WITH_CUDA: "ON"
CMAKE_BUILD_TYPE: "DebugOptimized"
WALBERLA_DOUBLE_ACCURACY: "OFF"
WALBERLA_BUILD_WITH_PARMETIS: "OFF"
WALBERLA_BUILD_WITH_METIS: "OFF"
+ WALBERLA_BUILD_WITH_CODEGEN: "ON"
+ WALBERLA_BUILD_WITH_PYTHON: "ON"
tags:
- cuda11
- docker
@@ -1422,7 +1447,7 @@ clang_11.0_serial:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:11.0
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -1445,7 +1470,7 @@ clang_11.0_mpionly:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:11.0
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -1466,7 +1491,7 @@ clang_11.0_hybrid:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:11.0
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -1486,7 +1511,7 @@ clang_11.0_serial_dbg:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:11.0
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -1510,7 +1535,7 @@ clang_11.0_mpionly_dbg:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:11.0
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -1532,7 +1557,7 @@ clang_11.0_hybrid_dbg:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:11.0
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -1552,12 +1577,20 @@ clang_11.0_hybrid_dbg:
clang_11.0_hybrid_dbg_sp:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:11.0
+ before_script:
+ - pip3 install lbmpy==1.0 jinja2 pytest
+ - cd python
+ - python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
+ - cd ..
+ - CC=gcc CXX=g++ pip3 install pycuda
variables:
WALBERLA_BUILD_WITH_CUDA: "ON"
CMAKE_BUILD_TYPE: "DebugOptimized"
WALBERLA_DOUBLE_ACCURACY: "OFF"
WALBERLA_BUILD_WITH_PARMETIS: "OFF"
WALBERLA_BUILD_WITH_METIS: "OFF"
+ WALBERLA_BUILD_WITH_CODEGEN: "ON"
+ WALBERLA_BUILD_WITH_PYTHON: "ON"
only:
variables:
- $ENABLE_NIGHTLY_BUILDS
@@ -1569,7 +1602,7 @@ clang_12.0_serial:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:12.0
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -1592,7 +1625,7 @@ clang_12.0_mpionly:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:12.0
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -1613,7 +1646,7 @@ clang_12.0_hybrid:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:12.0
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -1633,7 +1666,7 @@ clang_12.0_serial_dbg:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:12.0
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -1657,7 +1690,7 @@ clang_12.0_mpionly_dbg:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:12.0
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -1679,7 +1712,7 @@ clang_12.0_hybrid_dbg:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:12.0
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -1699,12 +1732,20 @@ clang_12.0_hybrid_dbg:
clang_12.0_hybrid_dbg_sp:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:12.0
+ before_script:
+ - pip3 install lbmpy==1.0 jinja2 pytest
+ - cd python
+ - python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
+ - cd ..
+ - CC=gcc CXX=g++ pip3 install pycuda
variables:
WALBERLA_BUILD_WITH_CUDA: "ON"
CMAKE_BUILD_TYPE: "DebugOptimized"
WALBERLA_DOUBLE_ACCURACY: "OFF"
WALBERLA_BUILD_WITH_PARMETIS: "OFF"
WALBERLA_BUILD_WITH_METIS: "OFF"
+ WALBERLA_BUILD_WITH_CODEGEN: "ON"
+ WALBERLA_BUILD_WITH_PYTHON: "ON"
only:
variables:
- $ENABLE_NIGHTLY_BUILDS
@@ -1716,7 +1757,7 @@ clang_13.0_serial:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:13.0
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -1739,7 +1780,7 @@ clang_13.0_mpionly:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:13.0
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -1760,7 +1801,7 @@ clang_13.0_hybrid:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:13.0
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -1777,7 +1818,7 @@ clang_13.0_serial_dbg:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:13.0
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -1798,7 +1839,7 @@ clang_13.0_mpionly_dbg:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:13.0
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -1817,7 +1858,7 @@ clang_13.0_hybrid_dbg:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:13.0
before_script:
- - pip3 install lbmpy==0.4.1 jinja2 pytest
+ - pip3 install lbmpy==1.0 jinja2 pytest
- cd python
- python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
- cd ..
@@ -1835,12 +1876,20 @@ clang_13.0_hybrid_dbg_sp:
extends: .build_template
image: i10git.cs.fau.de:5005/walberla/buildenvs/clang:13.0
stage: pretest
+ before_script:
+ - pip3 install lbmpy==1.0 jinja2 pytest
+ - cd python
+ - python3 -m pytest --junitxml=report.xml pystencils_walberla lbmpy_walberla
+ - cd ..
+ - CC=gcc CXX=g++ pip3 install pycuda
variables:
WALBERLA_BUILD_WITH_CUDA: "ON"
CMAKE_BUILD_TYPE: "DebugOptimized"
WALBERLA_DOUBLE_ACCURACY: "OFF"
WALBERLA_BUILD_WITH_PARMETIS: "OFF"
WALBERLA_BUILD_WITH_METIS: "OFF"
+ WALBERLA_BUILD_WITH_CODEGEN: "ON"
+ WALBERLA_BUILD_WITH_PYTHON: "ON"
tags:
- cuda11
- docker
diff --git a/apps/benchmarks/FlowAroundSphereCodeGen/FlowAroundSphereCodeGen.cpp b/apps/benchmarks/FlowAroundSphereCodeGen/FlowAroundSphereCodeGen.cpp
index 48f029dee..bdd1ccbe1 100644
--- a/apps/benchmarks/FlowAroundSphereCodeGen/FlowAroundSphereCodeGen.cpp
+++ b/apps/benchmarks/FlowAroundSphereCodeGen/FlowAroundSphereCodeGen.cpp
@@ -73,21 +73,21 @@ auto VelocityCallback = [](const Cell& pos, const shared_ptr< StructuredBlockFor
real_t inflow_velocity, const bool constant_inflow = true) {
if (constant_inflow)
{
- Vector3< real_t > result(inflow_velocity, 0.0, 0.0);
+ Vector3< real_t > result(inflow_velocity, real_c(0.0), real_c(0.0));
return result;
}
else
{
Cell globalCell;
CellInterval domain = SbF->getDomainCellBB();
- real_t h_y = real_c(domain.ySize());
- real_t h_z = real_c(domain.zSize());
+ auto h_y = real_c(domain.ySize());
+ auto h_z = real_c(domain.zSize());
SbF->transformBlockLocalToGlobalCell(globalCell, block, pos);
- real_t y1 = real_c(globalCell[1] - (h_y / 2.0 - 0.5));
- real_t z1 = real_c(globalCell[2] - (h_z / 2.0 - 0.5));
+ auto y1 = real_c(globalCell[1] - (h_y / 2.0 - 0.5));
+ auto z1 = real_c(globalCell[2] - (h_z / 2.0 - 0.5));
- real_t u = (inflow_velocity * real_c(16)) / (h_y * h_y * h_z * h_z) * (h_y / real_c(2.0) - y1) *
+ real_t u = (inflow_velocity * real_c(16.0)) / (h_y * h_y * h_z * h_z) * (h_y / real_c(2.0) - y1) *
(h_y / real_c(2.0) + y1) * (h_z / real_c(2.0) - z1) * (h_z / real_c(2.0) + z1);
Vector3< real_t > result(u, 0.0, 0.0);
@@ -151,9 +151,9 @@ int main(int argc, char** argv)
auto parameters = config->getOneBlock("Parameters");
const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(10));
- const real_t omega = parameters.getParameter< real_t >("omega", real_t(1.9));
- const real_t u_max = parameters.getParameter< real_t >("u_max", real_t(0.05));
- const real_t reynolds_number = parameters.getParameter< real_t >("reynolds_number", real_t(1000));
+ const real_t omega = parameters.getParameter< real_t >("omega", real_c(1.9));
+ const real_t u_max = parameters.getParameter< real_t >("u_max", real_c(0.05));
+ const real_t reynolds_number = parameters.getParameter< real_t >("reynolds_number", real_c(1000.0));
const uint_t diameter_sphere = parameters.getParameter< uint_t >("diameter_sphere", uint_t(5));
const bool constant_inflow = parameters.getParameter< bool >("constant_inflow", true);
@@ -162,8 +162,8 @@ int main(int argc, char** argv)
// create fields
BlockDataID pdfFieldID = blocks->addStructuredBlockData< PdfField_T >(pdfFieldAdder, "PDFs");
- BlockDataID velFieldID = field::addToStorage< VelocityField_T >(blocks, "velocity", real_t(0), field::fzyx);
- BlockDataID densityFieldID = field::addToStorage< ScalarField_T >(blocks, "density", real_t(0), field::fzyx);
+ BlockDataID velFieldID = field::addToStorage< VelocityField_T >(blocks, "velocity", real_c(0.0), field::fzyx);
+ BlockDataID densityFieldID = field::addToStorage< ScalarField_T >(blocks, "density", real_c(0.0), field::fzyx);
#if defined(WALBERLA_BUILD_WITH_CUDA)
BlockDataID pdfFieldIDGPU = cuda::addGPUFieldToStorage< PdfField_T >(blocks, pdfFieldID, "PDFs on GPU", true);
@@ -300,7 +300,7 @@ int main(int argc, char** argv)
timeloop.run();
simTimer.end();
WALBERLA_LOG_INFO_ON_ROOT("Simulation finished")
- auto time = simTimer.last();
+ auto time = real_c(simTimer.last());
auto nrOfCells = real_c(cellsPerBlock[0] * cellsPerBlock[1] * cellsPerBlock[2]);
auto mlupsPerProcess = nrOfCells * real_c(timesteps) / time * 1e-6;
WALBERLA_LOG_RESULT_ON_ROOT("MLUPS per process " << mlupsPerProcess)
diff --git a/apps/benchmarks/FlowAroundSphereCodeGen/FlowAroundSphereCodeGen.py b/apps/benchmarks/FlowAroundSphereCodeGen/FlowAroundSphereCodeGen.py
index 8d3b556bb..c170a8101 100644
--- a/apps/benchmarks/FlowAroundSphereCodeGen/FlowAroundSphereCodeGen.py
+++ b/apps/benchmarks/FlowAroundSphereCodeGen/FlowAroundSphereCodeGen.py
@@ -33,7 +33,8 @@ with CodeGeneration() as ctx:
'velocity': velocity_field
}
- lbm_config = LBMConfig(stencil=stencil, method=Method.CUMULANT, relaxation_rate=omega, galilean_correction=True,
+ lbm_config = LBMConfig(stencil=stencil, method=Method.CUMULANT, compressible=True,
+ relaxation_rate=omega, galilean_correction=True,
field_name='pdfs', streaming_pattern=streaming_pattern, output=output)
lbm_optimisation = LBMOptimisation(symbolic_field=pdfs, cse_global=False, cse_pdfs=False)
@@ -66,9 +67,9 @@ with CodeGeneration() as ctx:
generate_sweep(ctx, 'FlowAroundSphereCodeGen_MacroSetter', setter_assignments, target=target)
# boundaries
- ubb = UBB(lambda *args: None, dim=dim)
+ ubb = UBB(lambda *args: None, dim=dim, data_type=data_type)
ubb_data_handler = UBBAdditionalDataHandler(stencil, ubb)
- outflow = ExtrapolationOutflow(stencil[4], lb_method, streaming_pattern=streaming_pattern)
+ outflow = ExtrapolationOutflow(stencil[4], lb_method, streaming_pattern=streaming_pattern, data_type=data_type)
outflow_data_handler = OutflowAdditionalDataHandler(stencil, outflow, target=target)
generate_alternating_lbm_boundary(ctx, 'FlowAroundSphereCodeGen_UBB', ubb, lb_method,
diff --git a/apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py b/apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py
index 982aa6167..8274c20a5 100644
--- a/apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py
+++ b/apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py
@@ -1,22 +1,22 @@
import sympy as sp
from sympy.core.cache import clear_cache
import pystencils as ps
-from lbmpy.creationfunctions import create_lb_method_from_existing, create_lb_method
+from lbmpy.creationfunctions import LBMConfig, LBMOptimisation, Method, Stencil, create_lb_method
from lbmpy_walberla import generate_lattice_model
from pystencils_walberla import CodeGeneration
from pystencils_walberla import get_vectorize_instruction_set
from lbmpy.creationfunctions import create_lb_collision_rule
from lbmpy.moments import is_even, get_order, MOMENT_SYMBOLS
-from lbmpy.stencils import get_stencil
+from lbmpy.stencils import LBStencil
from collections import OrderedDict
with CodeGeneration() as ctx:
generatedMethod = 'TRTlike'
- #generatedMethod = 'D3Q27TRTlike'
- #generatedMethod = 'cumulant'
+ # generatedMethod = 'D3Q27TRTlike'
+ # generatedMethod = 'cumulant'
clear_cache()
@@ -28,7 +28,7 @@ with CodeGeneration() as ctx:
omegaVisc = sp.Symbol('omega_visc')
omegaBulk = ps.fields(f'omega_bulk: {dtype_string}[3D]', layout='fzyx')
omegaMagic = sp.Symbol('omega_magic')
- stencil = get_stencil('D3Q19', 'walberla')
+ stencil = LBStencil(Stencil.D3Q19)
x, y, z = MOMENT_SYMBOLS
one = sp.Rational(1, 1)
@@ -44,33 +44,43 @@ with CodeGeneration() as ctx:
]
# relaxation rate for first group of moments (1,x,y,z) is set to zero internally
- relaxation_rates=[omegaBulk.center_vector, omegaBulk.center_vector, omegaMagic, omegaVisc, omegaVisc, omegaMagic]
+ relaxation_rates = [omegaBulk.center_vector, omegaBulk.center_vector,
+ omegaMagic, omegaVisc, omegaVisc, omegaMagic]
- methodWithForce = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False,
- nested_moments=moments, relaxation_rates=relaxation_rates)
+ lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, continuous_equilibrium=False,
+ zero_centered=False, delta_equilibrium=False,
+ nested_moments=moments, relaxation_rates=relaxation_rates)
- #print(methodWithForce.relaxation_rates)
- #print(methodWithForce.moment_matrix)
- collision_rule = create_lb_collision_rule(lb_method=methodWithForce, optimization={'cse_global': True})
+ lbm_opt = LBMOptimisation(cse_global=True)
+
+ methodWithForce = create_lb_method(lbm_config=lbm_config)
+
+ # print(methodWithForce.relaxation_rates)
+ # print(methodWithForce.moment_matrix)
+ collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
generate_lattice_model(ctx, 'GeneratedLBM', collision_rule, field_layout='fzyx', cpu_vectorize_info=cpu_vectorize_info)
-
if generatedMethod == 'D3Q27TRTlike':
omegaVisc = sp.Symbol('omega_visc')
omegaBulk = ps.fields(f'omega_bulk: {dtype_string}[3D]', layout='fzyx')
omegaMagic = sp.Symbol('omega_magic')
- stencil = get_stencil('D3Q27', 'walberla')
+ stencil = LBStencil(Stencil.D3Q27)
relaxation_rates = [omegaVisc, omegaBulk.center_vector, omegaMagic, omegaVisc, omegaMagic, omegaVisc]
- methodWithForce = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False, weighted=True,
- compressible=False, relaxation_rates=relaxation_rates)
+ lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, maxwellian_moments=False, weighted=True,
+ compressible=False, relaxation_rates=relaxation_rates)
+
+ lbm_opt = LBMOptimisation(cse_global=True)
+
+ methodWithForce = create_lb_method(lbm_config=lbm_config)
- collision_rule = create_lb_collision_rule(lb_method=methodWithForce, optimization={'cse_global': True})
+ collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
generate_lattice_model(ctx, 'GeneratedLBM', collision_rule, field_layout='fzyx',
cpu_vectorize_info=cpu_vectorize_info)
+# using create_with_continuous_maxwellian_eq_moments won't work probably
if generatedMethod == 'cumulant':
x, y, z = MOMENT_SYMBOLS
@@ -122,7 +132,7 @@ with CodeGeneration() as ctx:
else:
return 1
- stencil = get_stencil('D3Q27', 'walberla')
+ stencil = LBStencil(Stencil.D3Q27)
omega = sp.Symbol('omega')
rr_dict = OrderedDict((c, get_relaxation_rate(c, omega))
@@ -139,6 +149,3 @@ with CodeGeneration() as ctx:
generate_lattice_model(ctx, 'GeneratedLBM', collision_rule, field_layout='fzyx',
cpu_vectorize_info=cpu_vectorize_info)
-
-
-
diff --git a/apps/benchmarks/FluidParticleCoupling/GeneratedLBMWithForce.py b/apps/benchmarks/FluidParticleCoupling/GeneratedLBMWithForce.py
index 9eb861b8c..4968f0670 100644
--- a/apps/benchmarks/FluidParticleCoupling/GeneratedLBMWithForce.py
+++ b/apps/benchmarks/FluidParticleCoupling/GeneratedLBMWithForce.py
@@ -1,22 +1,20 @@
import sympy as sp
from sympy.core.cache import clear_cache
import pystencils as ps
-from lbmpy.creationfunctions import create_lb_method_from_existing, create_lb_ast, create_lb_method
+from lbmpy.creationfunctions import LBMConfig, LBMOptimisation, ForceModel, Method, Stencil, create_lb_method
from lbmpy_walberla import generate_lattice_model
from pystencils_walberla import CodeGeneration
from pystencils_walberla import get_vectorize_instruction_set
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.forcemodels import Luo
+from lbmpy.stencils import LBStencil
from collections import OrderedDict
with CodeGeneration() as ctx:
forcing = (sp.symbols('fx'), 0, 0)
- forcemodel = Luo(forcing)
generatedMethod = 'TRTlike'
# generatedMethod = 'D3Q27TRTlike'
@@ -35,7 +33,7 @@ with CodeGeneration() as ctx:
omegaVisc = sp.Symbol('omega_visc')
omegaBulk = ps.fields(f'omega_bulk: {dtype_string}[3D]', layout='fzyx')
omegaMagic = sp.Symbol('omega_magic')
- stencil = get_stencil('D3Q19', 'walberla')
+ stencil = LBStencil(Stencil.D3Q19)
x, y, z = MOMENT_SYMBOLS
one = sp.Rational(1, 1)
@@ -51,29 +49,42 @@ with CodeGeneration() as ctx:
]
# relaxation rate for first group of moments (1,x,y,z) is set to zero internally
- relaxation_rates=[omegaBulk.center_vector, omegaBulk.center_vector, omegaMagic, omegaVisc, omegaVisc, omegaMagic]
+ relaxation_rates = [omegaBulk.center_vector, omegaBulk.center_vector,
+ omegaMagic, omegaVisc, omegaVisc, omegaMagic]
- methodWithForce = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False,
- force_model=forcemodel, nested_moments=moments, relaxation_rates=relaxation_rates)
+ lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, continuous_equilibrium=False,
+ zero_centered=False, delta_equilibrium=False,
+ force=forcing, force_model=ForceModel.LUO,
+ nested_moments=moments, relaxation_rates=relaxation_rates)
+
+ lbm_opt = LBMOptimisation(cse_global=True)
+
+ methodWithForce = create_lb_method(lbm_config=lbm_config)
+
+ # print(methodWithForce.relaxation_rates)
+ # print(methodWithForce.moment_matrix)
+ collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+ generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule, field_layout='fzyx',
+ cpu_vectorize_info=cpu_vectorize_info)
- #print(methodWithForce.relaxation_rates)
- #print(methodWithForce.moment_matrix)
- collision_rule = create_lb_collision_rule(lb_method=methodWithForce, optimization={'cse_global': True})
- generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule, field_layout='fzyx', cpu_vectorize_info=cpu_vectorize_info)
-
if generatedMethod == 'D3Q27TRTlike':
omegaVisc = sp.Symbol('omega_visc')
omegaBulk = ps.fields(f'omega_bulk: {dtype_string}[3D]', layout='fzyx')
omegaMagic = sp.Symbol('omega_magic')
- stencil = get_stencil('D3Q27', 'walberla')
+ stencil = LBStencil(Stencil.D3Q27)
relaxation_rates = [omegaVisc, omegaBulk.center_vector, omegaMagic, omegaVisc, omegaMagic, omegaVisc]
- methodWithForce = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False, weighted=True,
- compressible=False, force_model=forcemodel, relaxation_rates=relaxation_rates)
+ lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, maxwellian_moments=False, weighted=True,
+ force=forcing, force_model=ForceModel.LUO,
+ compressible=False, relaxation_rates=relaxation_rates)
+
+ lbm_opt = LBMOptimisation(cse_global=True)
- collision_rule = create_lb_collision_rule(lb_method=methodWithForce, optimization={'cse_global': True})
+ methodWithForce = create_lb_method(lbm_config=lbm_config)
+
+ collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule, field_layout='fzyx',
cpu_vectorize_info=cpu_vectorize_info)
@@ -128,7 +139,7 @@ with CodeGeneration() as ctx:
else:
return 1
- stencil = get_stencil('D3Q27', 'walberla')
+ stencil = LBStencil(Stencil.D3Q27)
omega = sp.Symbol('omega')
rr_dict = OrderedDict((c, get_relaxation_rate(c, omega))
@@ -199,7 +210,7 @@ with CodeGeneration() as ctx:
else:
return omegaMagic
- stencil = get_stencil('D3Q27', 'walberla')
+ stencil = LBStencil(Stencil.D3Q27)
omegaVisc = sp.Symbol('omega_visc')
omegaMagic = sp.Symbol('omega_magic')
@@ -219,4 +230,3 @@ with CodeGeneration() as ctx:
generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule, field_layout='fzyx',
cpu_vectorize_info=cpu_vectorize_info)
-
diff --git a/apps/benchmarks/PhaseFieldAllenCahn/benchmark_multiphase.cpp b/apps/benchmarks/PhaseFieldAllenCahn/benchmark_multiphase.cpp
index 5af04a355..b757f3a5b 100644
--- a/apps/benchmarks/PhaseFieldAllenCahn/benchmark_multiphase.cpp
+++ b/apps/benchmarks/PhaseFieldAllenCahn/benchmark_multiphase.cpp
@@ -86,14 +86,14 @@ int main(int argc, char** argv)
const std::string timeStepStrategy = parameters.getParameter< std::string >("timeStepStrategy", "normal");
const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(50));
const real_t remainingTimeLoggerFrequency =
- parameters.getParameter< real_t >("remainingTimeLoggerFrequency", 3.0);
+ parameters.getParameter< real_t >("remainingTimeLoggerFrequency", real_c(3.0));
const uint_t scenario = parameters.getParameter< uint_t >("scenario", uint_c(1));
const uint_t warmupSteps = parameters.getParameter< uint_t >("warmupSteps", uint_t(2));
#if defined(WALBERLA_BUILD_WITH_CUDA)
// CPU fields
- BlockDataID vel_field = field::addToStorage< VelocityField_T >(blocks, "vel", real_t(0), field::fzyx);
- BlockDataID phase_field = field::addToStorage< PhaseField_T >(blocks, "phase", real_t(0), field::fzyx);
+ BlockDataID vel_field = field::addToStorage< VelocityField_T >(blocks, "vel", real_c(0.0), field::fzyx);
+ BlockDataID phase_field = field::addToStorage< PhaseField_T >(blocks, "phase", real_c(0.0), field::fzyx);
// GPU fields
BlockDataID lb_phase_field_gpu = cuda::addGPUFieldToStorage< cuda::GPUField< real_t > >(
blocks, "lb phase field on GPU", Stencil_phase_T::Size, field::fzyx, 1);
@@ -105,11 +105,11 @@ int main(int argc, char** argv)
cuda::addGPUFieldToStorage< PhaseField_T >(blocks, phase_field, "phase field on GPU", true);
#else
BlockDataID lb_phase_field =
- field::addToStorage< PdfField_phase_T >(blocks, "lb phase field", real_t(0), field::fzyx);
+ field::addToStorage< PdfField_phase_T >(blocks, "lb phase field", real_c(0.0), field::fzyx);
BlockDataID lb_velocity_field =
- field::addToStorage< PdfField_hydro_T >(blocks, "lb velocity field", real_t(0), field::fzyx);
- BlockDataID vel_field = field::addToStorage< VelocityField_T >(blocks, "vel", real_t(0), field::fzyx);
- BlockDataID phase_field = field::addToStorage< PhaseField_T >(blocks, "phase", real_t(0), field::fzyx);
+ field::addToStorage< PdfField_hydro_T >(blocks, "lb velocity field", real_c(0.0), field::fzyx);
+ BlockDataID vel_field = field::addToStorage< VelocityField_T >(blocks, "vel", real_c(0.0), field::fzyx);
+ BlockDataID phase_field = field::addToStorage< PhaseField_T >(blocks, "phase", real_c(0.0), field::fzyx);
#endif
if (timeStepStrategy != "phase_only" && timeStepStrategy != "hydro_only" && timeStepStrategy != "kernel_only")
@@ -120,7 +120,7 @@ int main(int argc, char** argv)
auto bubbleParameters = config->getOneBlock("Bubble");
const Vector3< real_t > bubbleMidPoint =
bubbleParameters.getParameter< Vector3< real_t > >("bubbleMidPoint");
- const real_t bubbleRadius = bubbleParameters.getParameter< real_t >("bubbleRadius", 20.0);
+ const real_t bubbleRadius = bubbleParameters.getParameter< real_t >("bubbleRadius", real_c(20.0));
initPhaseField_bubble(blocks, phase_field, bubbleRadius, bubbleMidPoint);
}
else if (scenario == 2)
@@ -327,7 +327,7 @@ int main(int argc, char** argv)
#endif
simTimer.end();
WALBERLA_LOG_INFO_ON_ROOT("Simulation finished")
- auto time = simTimer.last();
+ auto time = real_c(simTimer.last());
auto nrOfCells = real_c(cellsPerBlock[0] * cellsPerBlock[1] * cellsPerBlock[2]);
auto mlupsPerProcess = nrOfCells * real_c(timesteps) / time * 1e-6;
WALBERLA_LOG_RESULT_ON_ROOT("MLUPS per process: " << mlupsPerProcess)
diff --git a/apps/benchmarks/PhaseFieldAllenCahn/multiphase_codegen.py b/apps/benchmarks/PhaseFieldAllenCahn/multiphase_codegen.py
index f74db05e9..e9bfa1dae 100644
--- a/apps/benchmarks/PhaseFieldAllenCahn/multiphase_codegen.py
+++ b/apps/benchmarks/PhaseFieldAllenCahn/multiphase_codegen.py
@@ -1,20 +1,20 @@
-from pystencils import fields, Target, TypedSymbol
+import numpy as np
+import sympy as sp
+
+from pystencils import fields, Target
+from pystencils.typing import TypedSymbol
from pystencils.simp import sympy_cse
from lbmpy import LBMConfig, LBStencil, Method, Stencil
-from lbmpy.creationfunctions import create_lb_method
+from lbmpy.creationfunctions import LBMOptimisation, create_lb_method, create_lb_update_rule
+from lbmpy.phasefield_allen_cahn.kernel_equations import initializer_kernel_phase_field_lb, \
+ initializer_kernel_hydro_lb, interface_tracking_force, hydrodynamic_force, add_interface_tracking_force, \
+ add_hydrodynamic_force, hydrodynamic_force_assignments
+from lbmpy.phasefield_allen_cahn.parameter_calculation import AllenCahnParameters
from pystencils_walberla import CodeGeneration, generate_sweep, generate_pack_info_for_field, generate_info_header
from lbmpy_walberla import generate_lb_pack_info
-from lbmpy.phasefield_allen_cahn.kernel_equations import initializer_kernel_phase_field_lb, \
- initializer_kernel_hydro_lb, interface_tracking_force, \
- hydrodynamic_force, get_collision_assignments_hydro, get_collision_assignments_phase
-
-from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel
-
-import numpy as np
-
with CodeGeneration() as ctx:
field_type = "float64" if ctx.double_accuracy else "float32"
@@ -25,22 +25,11 @@ with CodeGeneration() as ctx:
########################
# PARAMETER DEFINITION #
########################
-
- density_liquid = 1.0
- density_gas = 0.001
- surface_tension = 0.0001
- M = 0.02
-
- # phase-field parameter
- drho3 = (density_liquid - density_gas) / 3
- # interface thickness
- W = 5
- # coefficient related to surface tension
- beta = 12.0 * (surface_tension / W)
- # coefficient related to surface tension
- kappa = 1.5 * surface_tension * W
- # relaxation rate allen cahn (h)
- w_c = 1.0 / (0.5 + (3.0 * M))
+ # In the codegneration skript we only need the symbolic parameters
+ parameters = AllenCahnParameters(density_heavy=1.0, density_light=0.001,
+ dynamic_viscosity_heavy=0.03, dynamic_viscosity_light=0.03,
+ surface_tension=0.0001, mobility=0.02, interface_thickness=5,
+ gravitational_acceleration=1e-6)
########################
# FIELDS #
@@ -63,65 +52,73 @@ with CodeGeneration() as ctx:
# RELAXATION RATES AND EXTERNAL FORCES #
########################################
- # calculate the relaxation rate for the hydro lb as well as the body force
- density = density_gas + C.center * (density_liquid - density_gas)
-
+ # relaxation rate for interface tracking LB step
+ relaxation_rate_allen_cahn = 1.0 / (0.5 + (3.0 * parameters.symbolic_mobility))
+ # calculate the relaxation rate for hydrodynamic LB step
+ rho_L = parameters.symbolic_density_light
+ rho_H = parameters.symbolic_density_heavy
+ density = rho_L + C.center * (rho_H - rho_L)
# force acting on all phases of the model
- body_force = np.array([0, 1e-6, 0])
+ body_force = [0, 0, 0]
+ body_force[1] = parameters.symbolic_gravitational_acceleration * density
- relaxation_time = 0.03 + 0.5
- relaxation_rate = 1.0 / relaxation_time
+ omega = parameters.omega(C)
###############
# LBM METHODS #
###############
- lbm_config_phase = LBMConfig(stencil=stencil_phase, method=Method.SRT, relaxation_rate=w_c, compressible=True)
+ rates = [0.0]
+ rates += [relaxation_rate_allen_cahn] * stencil_phase.D
+ rates += [1.0] * (stencil_phase.Q - stencil_phase.D - 1)
+
+ lbm_config_phase = LBMConfig(stencil=stencil_phase, method=Method.MRT, compressible=True,
+ delta_equilibrium=False,
+ force=sp.symbols(f"F_:{stencil_phase.D}"), velocity_input=u,
+ weighted=True, relaxation_rates=rates,
+ output={'density': C_tmp}, kernel_type='stream_pull_collide')
method_phase = create_lb_method(lbm_config=lbm_config_phase)
- lbm_config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.MRT, weighted=True,
- relaxation_rates=[relaxation_rate, 1, 1, 1, 1, 1])
+ lbm_config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.MRT, compressible=False,
+ weighted=True, relaxation_rate=omega,
+ force=sp.symbols(f"F_:{stencil_hydro.D}"),
+ output={'velocity': u}, kernel_type='collide_stream_push')
method_hydro = create_lb_method(lbm_config=lbm_config_hydro)
# create the kernels for the initialization of the g and h field
- h_updates = initializer_kernel_phase_field_lb(h, C, u, method_phase, W)
- g_updates = initializer_kernel_hydro_lb(g, u, method_hydro)
-
- force_h = [f / 3 for f in interface_tracking_force(C, stencil_phase, W, fd_stencil=LBStencil(Stencil.D3Q27))]
- force_model_h = MultiphaseForceModel(force=force_h)
+ h_updates = initializer_kernel_phase_field_lb(method_phase, C, u, h, parameters)
+ h_updates = h_updates.new_with_substitutions(parameters.parameter_map())
- force_g = hydrodynamic_force(g, C, method_hydro, relaxation_time, density_liquid, density_gas, kappa, beta,
- body_force,
- fd_stencil=LBStencil(Stencil.D3Q27))
+ g_updates = initializer_kernel_hydro_lb(method_hydro, 1, u, g)
+ g_updates = g_updates.new_with_substitutions(parameters.parameter_map())
- force_model_g = MultiphaseForceModel(force=force_g, rho=density)
+ force_h = interface_tracking_force(C, stencil_phase, parameters)
+ hydro_force = hydrodynamic_force(g, C, method_hydro, parameters, body_force)
####################
# LBM UPDATE RULES #
####################
- phase_field_LB_step = get_collision_assignments_phase(lb_method=method_phase,
- velocity_input=u,
- output={'density': C_tmp},
- force_model=force_model_h,
- symbolic_fields={"symbolic_field": h,
- "symbolic_temporary_field": h_tmp},
- kernel_type='stream_pull_collide')
+ lbm_optimisation = LBMOptimisation(symbolic_field=h, symbolic_temporary_field=h_tmp)
+ phase_field_LB_step = create_lb_update_rule(lbm_config=lbm_config_phase,
+ lbm_optimisation=lbm_optimisation)
- phase_field_LB_step = sympy_cse(phase_field_LB_step)
+ phase_field_LB_step = add_interface_tracking_force(phase_field_LB_step, force_h)
+ phase_field_LB_step = phase_field_LB_step.new_with_substitutions(parameters.parameter_map())
+ phase_field_LB_step = sympy_cse(phase_field_LB_step)
# ---------------------------------------------------------------------------------------------------------
+ force_Assignments = hydrodynamic_force_assignments(g, u, C, method_hydro, parameters, body_force)
- hydro_LB_step = get_collision_assignments_hydro(lb_method=method_hydro,
- density=density,
- velocity_input=u,
- force_model=force_model_g,
- sub_iterations=2,
- symbolic_fields={"symbolic_field": g,
- "symbolic_temporary_field": g_tmp},
- kernel_type='collide_stream_push')
+ lbm_optimisation = LBMOptimisation(symbolic_field=g, symbolic_temporary_field=g_tmp)
+ hydro_LB_step = create_lb_update_rule(lbm_config=lbm_config_hydro,
+ lbm_optimisation=lbm_optimisation)
+
+ hydro_LB_step = add_hydrodynamic_force(hydro_LB_step, force_Assignments, C, g, parameters)
+ hydro_LB_step = hydro_LB_step.new_with_substitutions(parameters.parameter_map())
hydro_LB_step = sympy_cse(hydro_LB_step)
+
###################
# GENERATE SWEEPS #
###################
diff --git a/apps/benchmarks/UniformGridCPU/UniformGridCPU.cpp b/apps/benchmarks/UniformGridCPU/UniformGridCPU.cpp
index c5fef7d89..9916a71e1 100644
--- a/apps/benchmarks/UniformGridCPU/UniformGridCPU.cpp
+++ b/apps/benchmarks/UniformGridCPU/UniformGridCPU.cpp
@@ -84,8 +84,8 @@ int main(int argc, char** argv)
// Creating fields
BlockDataID pdfFieldId = blocks->addStructuredBlockData< PdfField_T >(pdfFieldAdder, "pdfs");
- BlockDataID velFieldId = field::addToStorage< VelocityField_T >(blocks, "vel", real_t(0), field::fzyx);
- BlockDataID densityFieldId = field::addToStorage< ScalarField_T >(blocks, "density", real_t(1.0), field::fzyx);
+ BlockDataID velFieldId = field::addToStorage< VelocityField_T >(blocks, "vel", real_c(0.0), field::fzyx);
+ BlockDataID densityFieldId = field::addToStorage< ScalarField_T >(blocks, "density", real_c(1.0), field::fzyx);
// Initialize velocity on cpu
if (initShearFlow)
@@ -270,7 +270,7 @@ int main(int argc, char** argv)
timeLoop.singleStep();
real_t remainingTimeLoggerFrequency =
- parameters.getParameter< real_t >("remainingTimeLoggerFrequency", -1.0); // in seconds
+ parameters.getParameter< real_t >("remainingTimeLoggerFrequency", real_c(-1.0)); // in seconds
if (remainingTimeLoggerFrequency > 0)
{
auto logger = timing::RemainingTimeLogger(timeLoop.getNrOfTimeSteps() * uint_c(outerIterations),
@@ -287,7 +287,7 @@ int main(int argc, char** argv)
timeLoop.run();
simTimer.end();
WALBERLA_LOG_INFO_ON_ROOT("Simulation finished")
- real_t time = simTimer.last();
+ auto time = real_c(simTimer.last());
WALBERLA_MPI_SECTION()
{
walberla::mpi::reduceInplace(time, walberla::mpi::MAX);
diff --git a/apps/benchmarks/UniformGridCPU/UniformGridCPU.py b/apps/benchmarks/UniformGridCPU/UniformGridCPU.py
index 2ea96523a..c8f500174 100644
--- a/apps/benchmarks/UniformGridCPU/UniformGridCPU.py
+++ b/apps/benchmarks/UniformGridCPU/UniformGridCPU.py
@@ -69,6 +69,7 @@ options_dict = {
'entropic': {
'method': Method.TRT_KBC_N4,
'compressible': True,
+ 'zero_centered': False,
'relaxation_rates': [omega, omega_free],
'entropic': True,
'entropic_newton_iterations': False
@@ -148,8 +149,6 @@ with CodeGeneration() as ctx:
collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
if optimize:
- collision_rule = insert_fast_divisions(collision_rule)
- collision_rule = insert_fast_sqrts(collision_rule)
collision_rule = insert_constants(collision_rule)
collision_rule = insert_zeros(collision_rule)
collision_rule = insert_aliases(collision_rule)
diff --git a/apps/benchmarks/UniformGridGPU/UniformGridGPU.cpp b/apps/benchmarks/UniformGridGPU/UniformGridGPU.cpp
index b12c534a4..4d2ee1afa 100644
--- a/apps/benchmarks/UniformGridGPU/UniformGridGPU.cpp
+++ b/apps/benchmarks/UniformGridGPU/UniformGridGPU.cpp
@@ -86,8 +86,8 @@ int main(int argc, char** argv)
// Creating fields
BlockDataID pdfFieldCpuID =
- field::addToStorage< PdfField_T >(blocks, "pdfs cpu", real_t(std::nan("")), field::fzyx);
- BlockDataID velFieldCpuID = field::addToStorage< VelocityField_T >(blocks, "vel", real_t(0), field::fzyx);
+ field::addToStorage< PdfField_T >(blocks, "pdfs cpu", real_c(std::nan("")), field::fzyx);
+ BlockDataID velFieldCpuID = field::addToStorage< VelocityField_T >(blocks, "vel", real_c(0.0), field::fzyx);
// Initialize velocity on cpu
if (initShearFlow)
@@ -294,7 +294,7 @@ int main(int argc, char** argv)
timeLoop.singleStep();
double remainingTimeLoggerFrequency =
- parameters.getParameter< double >("remainingTimeLoggerFrequency", -1.0); // in seconds
+ parameters.getParameter< real_t >("remainingTimeLoggerFrequency", real_c(-1.0)); // in seconds
if (remainingTimeLoggerFrequency > 0)
{
auto logger = timing::RemainingTimeLogger(timeLoop.getNrOfTimeSteps() * uint_c(outerIterations),
@@ -316,7 +316,7 @@ int main(int argc, char** argv)
cudaDeviceSynchronize();
simTimer.end();
WALBERLA_LOG_INFO_ON_ROOT("Simulation finished")
- auto time = simTimer.last();
+ auto time = real_c(simTimer.last());
auto nrOfCells = real_c(cellsPerBlock[0] * cellsPerBlock[1] * cellsPerBlock[2]);
auto mlupsPerProcess = nrOfCells * real_c(timesteps) / time * 1e-6;
diff --git a/apps/benchmarks/UniformGridGPU/UniformGridGPU.py b/apps/benchmarks/UniformGridGPU/UniformGridGPU.py
index ce61f5f0e..e8fa9906a 100644
--- a/apps/benchmarks/UniformGridGPU/UniformGridGPU.py
+++ b/apps/benchmarks/UniformGridGPU/UniformGridGPU.py
@@ -4,7 +4,7 @@ import pystencils as ps
from dataclasses import replace
-from pystencils.data_types import TypedSymbol
+from pystencils.typing import TypedSymbol
from pystencils.fast_approximation import insert_fast_sqrts, insert_fast_divisions
from lbmpy import LBMConfig, LBMOptimisation, LBStencil, Method, Stencil
@@ -78,6 +78,7 @@ options_dict = {
'entropic': {
'method': Method.TRT_KBC_N4,
'compressible': True,
+ 'zero_centered': False,
'relaxation_rates': [omega, omega_free],
'entropic': True,
'entropic_newton_iterations': False
@@ -184,7 +185,7 @@ with CodeGeneration() as ctx:
# Boundaries
noslip = NoSlip()
- ubb = UBB((0.05, 0, 0))
+ ubb = UBB((0.05, 0, 0), data_type=field_type)
generate_alternating_lbm_boundary(ctx, 'UniformGridGPU_NoSlip', noslip, lb_method, field_name=pdfs.name,
streaming_pattern=streaming_pattern, target=ps.Target.GPU)
diff --git a/apps/benchmarks/UniformGridGPU/old_ideas/UniformGridGPU.py b/apps/benchmarks/UniformGridGPU/old_ideas/UniformGridGPU.py
index c861cb155..619c9fbfc 100644
--- a/apps/benchmarks/UniformGridGPU/old_ideas/UniformGridGPU.py
+++ b/apps/benchmarks/UniformGridGPU/old_ideas/UniformGridGPU.py
@@ -7,7 +7,7 @@ from lbmpy.fieldaccess import StreamPullTwoFieldsAccessor
from pystencils_walberla import generate_pack_info_from_kernel
from lbmpy_walberla import generate_lattice_model, generate_boundary
from pystencils_walberla import CodeGeneration, generate_sweep
-from pystencils.data_types import TypedSymbol
+from pystencils.typing import TypedSymbol
from pystencils.fast_approximation import insert_fast_sqrts, insert_fast_divisions
from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp
index 979c87ff3..d7cb146a0 100644
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp
@@ -33,9 +33,9 @@ using PhaseField_T = GhostLayerField< real_t, 1 >;
using FlagField_T = FlagField< uint8_t >;
void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
- const real_t R = 10,
+ const real_t R = real_c(10.0),
const Vector3< real_t > bubbleMidPoint = Vector3< real_t >(0.0, 0.0, 0.0),
- const bool bubble = true, const real_t W = 5)
+ const bool bubble = true, const real_t W = real_c(5.0))
{
for (auto& block : *blocks)
{
@@ -43,21 +43,21 @@ void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, B
// clang-format off
WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, Cell globalCell;
blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
- real_t Ri = sqrt((globalCell[0] - bubbleMidPoint[0]) * (globalCell[0] - bubbleMidPoint[0]) +
- (globalCell[1] - bubbleMidPoint[1]) * (globalCell[1] - bubbleMidPoint[1]) +
- (globalCell[2] - bubbleMidPoint[2]) * (globalCell[2] - bubbleMidPoint[2]));
- if (bubble) phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(2.0 * (Ri - R) / W);
- else phaseField->get(x, y, z) = 0.5 - 0.5 * tanh(2.0 * (Ri - R) / W);
+ real_t Ri = sqrt((real_c(globalCell[0]) - bubbleMidPoint[0]) * (real_c(globalCell[0]) - bubbleMidPoint[0]) +
+ (real_c(globalCell[1]) - bubbleMidPoint[1]) * (real_c(globalCell[1]) - bubbleMidPoint[1]) +
+ (real_c(globalCell[2]) - bubbleMidPoint[2]) * (real_c(globalCell[2]) - bubbleMidPoint[2]));
+ if (bubble) phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * real_c(tanh(2.0 * (Ri - R) / W));
+ else phaseField->get(x, y, z) = real_c(0.5) - real_c(0.5) * real_c(tanh(2.0 * (Ri - R) / W));
)
// clang-format on
}
}
void init_Taylor_bubble(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
- const real_t D = 5, const real_t H = 2, const real_t DT = 20, const real_t Donut_x0 = 40)
+ const real_t D = real_c(5.0), const real_t H = real_c(2.0), const real_t DT = real_c(20.0), const real_t Donut_x0 = real_c(40.0))
{
- auto Mx = blocks->getDomainCellBB().xMax() / 2.0;
- auto Mz = blocks->getDomainCellBB().zMax() / 2.0;
+ auto Mx = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0);
+ auto Mz = real_c(blocks->getDomainCellBB().zMax()) / real_c(2.0);
for (auto& block : *blocks)
{
@@ -65,42 +65,42 @@ void init_Taylor_bubble(const shared_ptr< StructuredBlockStorage >& blocks, Bloc
WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(
phaseField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
- real_t Ri = D * sqrt(pow(H, 2) - pow(DT - sqrt(pow(globalCell[0] - Mx, 2) + pow(globalCell[2] - Mz, 2)), 2));
+ real_t Ri = D * real_c(sqrt(pow(H, 2) - pow(DT - sqrt(pow(real_c(globalCell[0]) - Mx, 2) + pow(real_c(globalCell[2]) - Mz, 2)), 2)));
- real_t shifter = atan2((globalCell[2] - Mz), (globalCell[0] - Mx));
+ real_t shifter = atan2((real_c(globalCell[2]) - Mz), (real_c(globalCell[0]) - Mx));
if (shifter < 0) shifter = shifter + 2 * math::pi;
- if ((globalCell[1] < Donut_x0 + Ri * sin(shifter / 2.0)) && (globalCell[1] > Donut_x0 - Ri)) {
+ if ((real_c(globalCell[1]) < Donut_x0 + Ri * sin(shifter / 2.0)) && (real_c(globalCell[1]) > Donut_x0 - Ri)) {
phaseField->get(x, y, z) = 0.0;
- } else { phaseField->get(x, y, z) = 1.0; })
+ } else { phaseField->get(x, y, z) = real_c(1.0); })
}
}
void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t R,
- real_t W = 5)
+ real_t W = real_c(5.0))
{
Vector3< real_t > bubbleMidPoint;
- auto X = blocks->getDomainCellBB().xMax();
- auto Y = blocks->getDomainCellBB().yMax();
- auto Z = blocks->getDomainCellBB().zMax();
+ auto X = real_c(blocks->getDomainCellBB().xMax());
+ auto Y = real_c(blocks->getDomainCellBB().yMax());
+ auto Z = real_c(blocks->getDomainCellBB().zMax());
// 20 percent from the top are filled with the gas phase
- real_t gas_top = Y - Y / 5.0;
+ real_t gas_top = Y - Y / real_c(5.0);
// Diameter of the bubble
- real_t D = R * 2;
+ real_t D = R * real_c(2.0);
// distance in between the bubbles
- int dist = 4;
- auto nx = static_cast< unsigned int >(floor(X / (D + dist * W)));
- auto nz = static_cast< unsigned int >(floor(Z / (D + dist * W)));
+ real_t dist = real_c(4.0);
+ auto nx = uint_c(floor(X / (D + dist * W)));
+ auto nz = uint_c(floor(Z / (D + dist * W)));
// fluctuation of the bubble radii
std::vector< std::vector< real_t > > fluctuation_radius(nx, std::vector< real_t >(nz, 0.0));
std::vector< std::vector< real_t > > fluctuation_pos(nx, std::vector< real_t >(nz, 0.0));
- real_t max_fluctuation_radius = R / 5;
- real_t max_fluctuation_pos = (dist * W) / 3.0;
+ real_t max_fluctuation_radius = R / real_c(5.0);
+ real_t max_fluctuation_pos = (dist * W) / real_c(3.0);
for (unsigned int i = 0; i < nx; ++i)
{
for (unsigned int j = 0; j < nz; ++j)
@@ -118,35 +118,35 @@ void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, Block
for (unsigned int i = 0; i < nx; ++i) {
for (unsigned int j = 0; j < nz; ++j)
{
- bubbleMidPoint[0] = (i + 1) * (D + (dist * W)) - (D + (dist * W)) / 2.0 + fluctuation_pos[i][j];
+ bubbleMidPoint[0] = real_c(i + 1) * (D + (dist * W)) - (D + (dist * W)) / real_c(2.0) + fluctuation_pos[i][j];
bubbleMidPoint[1] = R + W + 4;
- bubbleMidPoint[2] = (j + 1) * (D + (dist * W)) - (D + (dist * W)) / 2.0 + fluctuation_pos[i][j];
+ bubbleMidPoint[2] = real_c(j + 1) * (D + (dist * W)) - (D + (dist * W)) / real_c(2.0) + fluctuation_pos[i][j];
- real_t Ri = sqrt((globalCell[0] - bubbleMidPoint[0]) * (globalCell[0] - bubbleMidPoint[0]) +
- (globalCell[1] - bubbleMidPoint[1]) * (globalCell[1] - bubbleMidPoint[1]) +
- (globalCell[2] - bubbleMidPoint[2]) * (globalCell[2] - bubbleMidPoint[2]));
- if (globalCell[0] >= i * (D + dist * W) && globalCell[0] <= (i + 1) * (D + dist * W) &&
- globalCell[2] >= j * (D + dist * W) && globalCell[2] <= (j + 1) * (D + dist * W))
- phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(2.0 * (Ri - (R - fluctuation_radius[i][j])) / W);
+ real_t Ri = sqrt((real_c(globalCell[0]) - bubbleMidPoint[0]) * (real_c(globalCell[0]) - bubbleMidPoint[0]) +
+ (real_c(globalCell[1]) - bubbleMidPoint[1]) * (real_c(globalCell[1]) - bubbleMidPoint[1]) +
+ (real_c(globalCell[2]) - bubbleMidPoint[2]) * (real_c(globalCell[2]) - bubbleMidPoint[2]));
+ if (real_c(globalCell[0]) >= real_c(i) * (D + dist * W) && real_c(globalCell[0]) <= real_c(i + 1) * (D + dist * W) &&
+ real_c(globalCell[2]) >= real_c(j) * (D + dist * W) && real_c(globalCell[2]) <= real_c(j + 1) * (D + dist * W))
+ phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * tanh(real_c(2.0) * (Ri - (R - fluctuation_radius[i][j])) / W);
- if (globalCell[0] > nx * (D + dist * W)) phaseField->get(x, y, z) = 1.0;
- if (globalCell[2] > nz * (D + dist * W)) phaseField->get(x, y, z) = 1.0;
+ if (real_c(globalCell[0]) > real_c(nx) * (D + dist * W)) phaseField->get(x, y, z) = real_c(1.0);
+ if (real_c(globalCell[2]) > real_c(nz) * (D + dist * W)) phaseField->get(x, y, z) = real_c(1.0);
}
}
- if (globalCell[1] > gas_top) {
- phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(2.0 * (gas_top + 10 - globalCell[1]) / W);
+ if (real_c(globalCell[1]) > gas_top) {
+ phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * tanh(real_c(2.0) * (gas_top + 10 - real_c(globalCell[1])) / W);
})
}
}
void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
- const real_t W = 5, const bool pipe = true)
+ const real_t W = real_c(5.0), const bool pipe = true)
{
- auto X = blocks->getDomainCellBB().xMax();
- auto Z = blocks->getDomainCellBB().zMax();
- auto halfY = (blocks->getDomainCellBB().yMax()) / 2.0;
- real_t perturbation = 0.05;
+ auto X = real_c(blocks->getDomainCellBB().xMax());
+ auto Z = real_c(blocks->getDomainCellBB().zMax());
+ auto halfY = real_c(blocks->getDomainCellBB().yMax()) / real_c(2.0);
+ auto perturbation = real_c(0.05);
if (pipe)
{
@@ -155,10 +155,10 @@ void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, Bloc
auto phaseField = block.getData< PhaseField_T >(phaseFieldID);
WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(
phaseField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
- real_t R = sqrt((globalCell[0] - X / 2) * (globalCell[0] - X / 2) +
- (globalCell[2] - Z / 2) * (globalCell[2] - Z / 2));
- if (R > X) R = X; real_t tmp = perturbation * X * cos((2.0 * math::pi * R) / X);
- phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(((globalCell[1] - halfY) + tmp) / (W / 2.0));)
+ real_t R = sqrt((real_c(globalCell[0]) - X / 2) * (real_c(globalCell[0]) - X / 2) +
+ (real_c(globalCell[2]) - Z / 2) * (real_c(globalCell[2]) - Z / 2));
+ if (R > X) R = X; real_t tmp = perturbation * X * cos((real_c(2.0) * math::pi * R) / X);
+ phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * tanh(((real_c(globalCell[1]) - halfY) + tmp) / (W / real_c(2.0)));)
}
}
else
@@ -169,8 +169,8 @@ void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, Bloc
WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(
phaseField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
real_t tmp = perturbation * X *
- (cos((2.0 * math::pi * globalCell[0]) / X) + cos((2.0 * math::pi * globalCell[2]) / X));
- phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(((globalCell[1] - halfY) - tmp) / (W / 2.0));)
+ (cos((real_c(2.0) * math::pi * real_c(globalCell[0])) / X) + cos((real_c(2.0) * math::pi * real_c(globalCell[2])) / X));
+ phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * tanh(((real_c(globalCell[1]) - halfY) - tmp) / (W / real_c(2.0)));)
}
}
}
@@ -182,12 +182,12 @@ void initTubeWithCylinder(const shared_ptr< StructuredBlockStorage >& blocks, Bl
{
if (eccentricity_or_pipe_ratio)
{
- auto Mx = blocks->getDomainCellBB().xMax() / 2.0;
- auto Mz = blocks->getDomainCellBB().zMax() / 2.0;
+ auto Mx = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0);
+ auto Mz = real_c(blocks->getDomainCellBB().zMax()) / real_c(2.0);
- auto R_outer = blocks->getDomainCellBB().xMax() / 2.0 + 1.0;
+ auto R_outer = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0) + real_c(1.0);
- real_t const shift = eccentricity * Mx / 2;
+ auto const shift = real_c(eccentricity * Mx / real_c(2.0));
for (auto& block : *blocks)
{
@@ -196,31 +196,31 @@ void initTubeWithCylinder(const shared_ptr< StructuredBlockStorage >& blocks, Bl
WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, Cell globalCell;
blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
real_t R1;
- if (globalCell[1] <= start_transition) {
- R1 = sqrt((globalCell[0] - Mx) * (globalCell[0] - Mx) + (globalCell[2] - Mz) * (globalCell[2] - Mz));
- } else if (globalCell[1] > start_transition && globalCell[1] < start_transition + length_transition) {
- real_t tmp = math::pi * (globalCell[1] - start_transition) / (length_transition);
- real_t shift_tmp = shift * 0.5 * (1 - cos(tmp));
- R1 = sqrt((globalCell[0] - Mx - shift_tmp) * (globalCell[0] - Mx - shift_tmp) +
- (globalCell[2] - Mz) * (globalCell[2] - Mz));
+ if (real_c(globalCell[1]) <= start_transition) {
+ R1 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz));
+ } else if (real_c(globalCell[1]) > start_transition && real_c(globalCell[1]) < start_transition + length_transition) {
+ real_t tmp = math::pi * (real_c(globalCell[1]) - start_transition) / (length_transition);
+ real_t shift_tmp = shift * real_c(0.5) * (1 - cos(tmp));
+ R1 = sqrt((real_c(globalCell[0]) - Mx - shift_tmp) * (real_c(globalCell[0]) - Mx - shift_tmp) +
+ (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz));
} else {
- R1 = sqrt((globalCell[0] - Mx - shift) * (globalCell[0] - Mx - shift) +
- (globalCell[2] - Mz) * (globalCell[2] - Mz));
+ R1 = sqrt((real_c(globalCell[0]) - Mx - shift) * (real_c(globalCell[0]) - Mx - shift) +
+ (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz));
}
- real_t R2 = sqrt((globalCell[0] - Mx) * (globalCell[0] - Mx) + (globalCell[2] - Mz) * (globalCell[2] - Mz));
+ real_t R2 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz));
if (R1 < R_in) addFlag(flagField->get(x, y, z), boundaryFlag);
if (R2 > R_outer) addFlag(flagField->get(x, y, z), boundaryFlag);)
}
}
else
{
- auto Mx = blocks->getDomainCellBB().xMax() / 2.0;
- auto Mz = blocks->getDomainCellBB().zMax() / 2.0;
+ auto Mx = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0);
+ auto Mz = real_c(blocks->getDomainCellBB().zMax()) / real_c(2.0);
- auto R_outer = blocks->getDomainCellBB().xMax() / 2.0 + 1.0;
+ auto R_outer = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0) + real_c(1.0);
- real_t const shift = eccentricity * R_in;
+ auto const shift = real_c(eccentricity * R_in);
real_t R_tmp;
for (auto& block : *blocks)
@@ -229,17 +229,17 @@ void initTubeWithCylinder(const shared_ptr< StructuredBlockStorage >& blocks, Bl
auto boundaryFlag = flagField->getOrRegisterFlag(boundaryFlagUID);
WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(
flagField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
- if (globalCell[1] <= start_transition) {
+ if (real_c(globalCell[1]) <= start_transition) {
R_tmp = R_in;
- } else if (globalCell[1] > start_transition && globalCell[1] < start_transition + length_transition) {
- real_t tmp = math::pi * (globalCell[1] - start_transition) / (length_transition);
- real_t shift_tmp = shift * 0.5 * (1 - cos(tmp));
+ } else if (real_c(globalCell[1]) > start_transition && real_c(globalCell[1]) < start_transition + length_transition) {
+ real_t tmp = math::pi * (real_c(globalCell[1]) - start_transition) / (length_transition);
+ real_t shift_tmp = shift * real_c(0.5) * (1 - cos(tmp));
R_tmp = R_in + shift_tmp;
} else {
R_tmp = R_in + shift;
}
- real_t R2 = sqrt((globalCell[0] - Mx) * (globalCell[0] - Mx) + (globalCell[2] - Mz) * (globalCell[2] - Mz));
+ real_t R2 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz));
if (R2 < R_tmp) addFlag(flagField->get(x, y, z), boundaryFlag);
if (R2 > R_outer) addFlag(flagField->get(x, y, z), boundaryFlag);)
}
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp
index e9f21f0f0..184672199 100644
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp
@@ -79,7 +79,7 @@ int main(int argc, char** argv)
const std::string timeStepStrategy = parameters.getParameter< std::string >("timeStepStrategy", "normal");
const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(50));
const real_t remainingTimeLoggerFrequency =
- parameters.getParameter< real_t >("remainingTimeLoggerFrequency", 3.0);
+ parameters.getParameter< real_t >("remainingTimeLoggerFrequency", real_c(3.0));
const uint_t scenario = parameters.getParameter< uint_t >("scenario", uint_c(1));
Vector3< int > overlappingWidth =
@@ -90,11 +90,11 @@ int main(int argc, char** argv)
///////////////////////
BlockDataID lb_phase_field =
- field::addToStorage< PdfField_phase_T >(blocks, "lb phase field", real_t(0), field::fzyx);
+ field::addToStorage< PdfField_phase_T >(blocks, "lb phase field", real_c(0.0), field::fzyx);
BlockDataID lb_velocity_field =
- field::addToStorage< PdfField_hydro_T >(blocks, "lb velocity field", real_t(0), field::fzyx);
- BlockDataID vel_field = field::addToStorage< VelocityField_T >(blocks, "vel", real_t(0), field::fzyx);
- BlockDataID phase_field = field::addToStorage< PhaseField_T >(blocks, "phase", real_t(0), field::fzyx);
+ field::addToStorage< PdfField_hydro_T >(blocks, "lb velocity field", real_c(0.0), field::fzyx);
+ BlockDataID vel_field = field::addToStorage< VelocityField_T >(blocks, "vel", real_c(0.0), field::fzyx);
+ BlockDataID phase_field = field::addToStorage< PhaseField_T >(blocks, "phase", real_c(0.0), field::fzyx);
BlockDataID flagFieldID = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
@@ -114,7 +114,7 @@ int main(int argc, char** argv)
{
auto bubbleParameters = config->getOneBlock("Bubble");
const Vector3< real_t > bubbleMidPoint = bubbleParameters.getParameter< Vector3< real_t > >("bubbleMidPoint");
- const real_t bubbleRadius = bubbleParameters.getParameter< real_t >("bubbleRadius", 20.0);
+ const real_t bubbleRadius = bubbleParameters.getParameter< real_t >("bubbleRadius", real_c(20.0));
const bool bubble = bubbleParameters.getParameter< bool >("bubble", true);
initPhaseField_sphere(blocks, phase_field, bubbleRadius, bubbleMidPoint, bubble);
}
@@ -300,7 +300,7 @@ int main(int argc, char** argv)
timeLoop->run();
simTimer.end();
WALBERLA_LOG_INFO_ON_ROOT("Simulation finished")
- auto time = simTimer.last();
+ auto time = real_c(simTimer.last());
auto nrOfCells = real_c(cellsPerBlock[0] * cellsPerBlock[1] * cellsPerBlock[2]);
auto mlupsPerProcess = nrOfCells * real_c(timesteps) / time * 1e-6;
WALBERLA_LOG_RESULT_ON_ROOT("MLUPS per process: " << mlupsPerProcess)
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py
index df7771178..a2dba9cb5 100644
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py
@@ -1,19 +1,18 @@
import numpy as np
import sympy as sp
-from pystencils import Assignment, fields, TypedSymbol, Target
-from pystencils.astnodes import Block, Conditional
-from pystencils.simp import sympy_cse
+from pystencils import fields, Target
+from pystencils.typing import TypedSymbol
-from lbmpy import LBMConfig, LBStencil, Method, Stencil
+from lbmpy import LBMConfig, LBMOptimisation, LBStencil, Method, Stencil, create_lb_update_rule
from lbmpy.creationfunctions import create_lb_method
from lbmpy.boundaries import NoSlip
from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle
-from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel
from lbmpy.phasefield_allen_cahn.kernel_equations import initializer_kernel_phase_field_lb, \
- initializer_kernel_hydro_lb, interface_tracking_force, hydrodynamic_force, get_collision_assignments_hydro, \
- get_collision_assignments_phase
+ initializer_kernel_hydro_lb, interface_tracking_force, hydrodynamic_force, add_interface_tracking_force, \
+ add_hydrodynamic_force, hydrodynamic_force_assignments
+from lbmpy.phasefield_allen_cahn.parameter_calculation import AllenCahnParameters
import pystencils_walberla
from pystencils_walberla import CodeGeneration, generate_sweep, generate_pack_info_for_field, generate_info_header
@@ -29,28 +28,8 @@ with CodeGeneration() as ctx:
stencil_hydro = LBStencil(Stencil.D3Q27)
assert (stencil_phase.D == stencil_hydro.D)
- ########################
- # PARAMETER DEFINITION #
- ########################
-
- density_liquid = sp.Symbol("rho_H")
- density_gas = sp.Symbol("rho_L")
-
- surface_tension = sp.Symbol("sigma")
- mobility = sp.Symbol("mobility")
-
- gravitational_acceleration = sp.Symbol("gravity")
-
- relaxation_time_liquid = sp.Symbol("tau_H")
- relaxation_time_gas = sp.Symbol("tau_L")
-
- # phase-field parameter
- drho3 = (density_liquid - density_gas) / 3
- # interface thickness
- W = sp.Symbol("interface_thickness")
- # coefficients related to surface tension
- beta = 12.0 * (surface_tension / W)
- kappa = 1.5 * surface_tension * W
+ # In the codegneration skript we only need the symbolic parameters
+ parameters = AllenCahnParameters(0, 0, 0, 0, 0)
########################
# FIELDS #
@@ -74,73 +53,65 @@ with CodeGeneration() as ctx:
########################################
# relaxation rate for interface tracking LB step
- relaxation_rate_allen_cahn = 1.0 / (0.5 + (3.0 * mobility))
+ relaxation_rate_allen_cahn = 1.0 / (0.5 + (3.0 * parameters.symbolic_mobility))
# calculate the relaxation rate for hydrodynamic LB step
- density = density_gas + C.center * (density_liquid - density_gas)
+ rho_L = parameters.symbolic_density_light
+ rho_H = parameters.symbolic_density_heavy
+ density = rho_L + C.center * (rho_H - rho_L)
# force acting on all phases of the model
- body_force = np.array([0, gravitational_acceleration * density, 0])
- # calculation of the relaxation time via viscosities
- # viscosity = viscosity_gas * viscosity_liquid + C.center\
- # * (density_liquid*viscosity_liquid - viscosity_liquid*viscosity_gas)
- # relaxation_time = 3 * viscosity / density + 0.5
-
- relaxation_time = 0.5 + relaxation_time_gas + C.center * (relaxation_time_liquid - relaxation_time_gas)
- # calculate the relaxation time if the phase-field has values over one
- # relaxation_rate = 1.0 / relaxation_time
- relaxation_rate = sp.Symbol("s8")
- relaxation_rate_cutoff = sp.Piecewise((1 / (0.5 + relaxation_time_liquid), C.center > 0.999), # True value
- (1 / relaxation_time, True)) # Else value
+ body_force = [0, 0, 0]
+ body_force[1] = parameters.symbolic_gravitational_acceleration * density
+
+ omega = parameters.omega(C)
###############
# LBM METHODS #
###############
- lbm_config_phase = LBMConfig(stencil=stencil_phase, method=Method.MRT, compressible=True, weighted=True,
- relaxation_rates=[1, 1, 1, 1, 1, 1])
+
+ rates = [0.0]
+ rates += [relaxation_rate_allen_cahn] * stencil_phase.D
+ rates += [1.0] * (stencil_phase.Q - stencil_phase.D - 1)
+
+ lbm_config_phase = LBMConfig(stencil=stencil_phase, method=Method.MRT, compressible=True,
+ delta_equilibrium=False,
+ force=sp.symbols(f"F_:{stencil_phase.D}"), velocity_input=u,
+ weighted=True, relaxation_rates=rates,
+ output={'density': C_tmp}, kernel_type='stream_pull_collide')
method_phase = create_lb_method(lbm_config=lbm_config_phase)
- method_phase.set_first_moment_relaxation_rate(relaxation_rate_allen_cahn)
- lbm_config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.MRT, weighted=True,
- relaxation_rates=[relaxation_rate, 1, 1, 1, 1, 1])
+ lbm_config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.MRT, compressible=False,
+ weighted=True, relaxation_rate=omega,
+ force=sp.symbols(f"F_:{stencil_hydro.D}"),
+ output={'velocity': u}, kernel_type='collide_stream_push')
method_hydro = create_lb_method(lbm_config=lbm_config_hydro)
# create the kernels for the initialization of the g and h field
- h_updates = initializer_kernel_phase_field_lb(h, C, u, method_phase, W)
- g_updates = initializer_kernel_hydro_lb(g, u, method_hydro)
-
- force_h = [f / 3 for f in interface_tracking_force(C, stencil_phase, W)]
- force_model_h = MultiphaseForceModel(force=force_h)
+ h_updates = initializer_kernel_phase_field_lb(method_phase, C, u, h, parameters)
+ g_updates = initializer_kernel_hydro_lb(method_hydro, 1, u, g)
- force_g = hydrodynamic_force(g, C, method_hydro, relaxation_time,
- density_liquid, density_gas, kappa, beta, body_force)
- force_model_g = MultiphaseForceModel(force=force_g, rho=density)
+ force_h = interface_tracking_force(C, stencil_phase, parameters)
+ hydro_force = hydrodynamic_force(g, C, method_hydro, parameters, body_force)
####################
# LBM UPDATE RULES #
####################
- phase_field_LB_step = get_collision_assignments_phase(lb_method=method_phase,
- velocity_input=u,
- output={'density': C_tmp},
- force_model=force_model_h,
- symbolic_fields={"symbolic_field": h,
- "symbolic_temporary_field": h_tmp},
- kernel_type='stream_pull_collide')
+ lbm_optimisation = LBMOptimisation(symbolic_field=h, symbolic_temporary_field=h_tmp)
+ allen_cahn_update_rule = create_lb_update_rule(lbm_config=lbm_config_phase,
+ lbm_optimisation=lbm_optimisation)
- phase_field_LB_step = sympy_cse(phase_field_LB_step)
+ allen_cahn_update_rule = add_interface_tracking_force(allen_cahn_update_rule, force_h)
# ---------------------------------------------------------------------------------------------------------
- hydro_LB_step = get_collision_assignments_hydro(lb_method=method_hydro,
- density=density,
- velocity_input=u,
- force_model=force_model_g,
- sub_iterations=2,
- symbolic_fields={"symbolic_field": g,
- "symbolic_temporary_field": g_tmp},
- kernel_type='collide_stream_push')
+ force_Assignments = hydrodynamic_force_assignments(g, u, C, method_hydro, parameters, body_force)
+
+ lbm_optimisation = LBMOptimisation(symbolic_field=g, symbolic_temporary_field=g_tmp)
+ hydro_lb_update_rule = create_lb_update_rule(lbm_config=lbm_config_hydro,
+ lbm_optimisation=lbm_optimisation)
- hydro_LB_step.set_sub_expressions_from_dict({**{relaxation_rate: relaxation_rate_cutoff},
- **hydro_LB_step.subexpressions_dict})
+ hydro_lb_update_rule = add_hydrodynamic_force(hydro_lb_update_rule, force_Assignments, C, g, parameters)
- contact_angle = ContactAngle(contact_angle_in_degrees, W)
+ contact_angle = ContactAngle(contact_angle_in_degrees, parameters.symbolic_interface_thickness,
+ data_type=field_type)
###################
# GENERATE SWEEPS #
@@ -171,13 +142,13 @@ with CodeGeneration() as ctx:
generate_sweep(ctx, 'initialize_phase_field_distributions', h_updates, target=Target.CPU)
generate_sweep(ctx, 'initialize_velocity_based_distributions', g_updates, target=Target.CPU)
- generate_sweep(ctx, 'phase_field_LB_step', phase_field_LB_step,
+ generate_sweep(ctx, 'phase_field_LB_step', allen_cahn_update_rule,
field_swaps=[(h, h_tmp), (C, C_tmp)],
inner_outer_split=True,
target=Target.CPU)
generate_boundary(ctx, 'phase_field_LB_NoSlip', NoSlip(), method_phase, target=Target.CPU, streaming_pattern='pull')
- generate_sweep(ctx, 'hydro_LB_step', hydro_LB_step,
+ generate_sweep(ctx, 'hydro_LB_step', hydro_lb_update_rule,
field_swaps=[(g, g_tmp)],
inner_outer_split=True,
target=Target.CPU)
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/InitializerFunctions.cpp b/apps/showcases/PhaseFieldAllenCahn/GPU/InitializerFunctions.cpp
index 979c87ff3..1d15f6a80 100644
--- a/apps/showcases/PhaseFieldAllenCahn/GPU/InitializerFunctions.cpp
+++ b/apps/showcases/PhaseFieldAllenCahn/GPU/InitializerFunctions.cpp
@@ -33,9 +33,9 @@ using PhaseField_T = GhostLayerField< real_t, 1 >;
using FlagField_T = FlagField< uint8_t >;
void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
- const real_t R = 10,
+ const real_t R = real_c(10.0),
const Vector3< real_t > bubbleMidPoint = Vector3< real_t >(0.0, 0.0, 0.0),
- const bool bubble = true, const real_t W = 5)
+ const bool bubble = true, const real_t W = real_c(5.0))
{
for (auto& block : *blocks)
{
@@ -43,21 +43,21 @@ void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, B
// clang-format off
WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, Cell globalCell;
blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
- real_t Ri = sqrt((globalCell[0] - bubbleMidPoint[0]) * (globalCell[0] - bubbleMidPoint[0]) +
- (globalCell[1] - bubbleMidPoint[1]) * (globalCell[1] - bubbleMidPoint[1]) +
- (globalCell[2] - bubbleMidPoint[2]) * (globalCell[2] - bubbleMidPoint[2]));
- if (bubble) phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(2.0 * (Ri - R) / W);
- else phaseField->get(x, y, z) = 0.5 - 0.5 * tanh(2.0 * (Ri - R) / W);
+ real_t Ri = sqrt((real_c(real_c(globalCell[0])) - bubbleMidPoint[0]) * (real_c(real_c(globalCell[0])) - bubbleMidPoint[0]) +
+ (real_c(globalCell[1]) - bubbleMidPoint[1]) * (real_c(globalCell[1]) - bubbleMidPoint[1]) +
+ (real_c(globalCell[2]) - bubbleMidPoint[2]) * (real_c(globalCell[2]) - bubbleMidPoint[2]));
+ if (bubble) phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * real_c(tanh(real_c(2.0) * (Ri - R) / W));
+ else phaseField->get(x, y, z) = real_c(0.5) - real_c(0.5) * tanh(real_c(2.0) * (Ri - R) / W);
)
// clang-format on
}
}
void init_Taylor_bubble(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
- const real_t D = 5, const real_t H = 2, const real_t DT = 20, const real_t Donut_x0 = 40)
+ const real_t D = 5.0, const real_t H = 2.0, const real_t DT = 20.0, const real_t Donut_x0 = 40.0)
{
- auto Mx = blocks->getDomainCellBB().xMax() / 2.0;
- auto Mz = blocks->getDomainCellBB().zMax() / 2.0;
+ auto Mx = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0);
+ auto Mz = real_c(blocks->getDomainCellBB().zMax()) / real_c(2.0);
for (auto& block : *blocks)
{
@@ -65,42 +65,42 @@ void init_Taylor_bubble(const shared_ptr< StructuredBlockStorage >& blocks, Bloc
WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(
phaseField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
- real_t Ri = D * sqrt(pow(H, 2) - pow(DT - sqrt(pow(globalCell[0] - Mx, 2) + pow(globalCell[2] - Mz, 2)), 2));
+ real_t Ri = D * real_c(sqrt(pow(H, 2) - pow(DT - sqrt(pow(real_c(real_c(globalCell[0])) - Mx, 2) + pow(real_c(globalCell[2]) - Mz, 2)), 2)));
- real_t shifter = atan2((globalCell[2] - Mz), (globalCell[0] - Mx));
+ real_t shifter = atan2((real_c(globalCell[2]) - Mz), (real_c(globalCell[0]) - Mx));
if (shifter < 0) shifter = shifter + 2 * math::pi;
- if ((globalCell[1] < Donut_x0 + Ri * sin(shifter / 2.0)) && (globalCell[1] > Donut_x0 - Ri)) {
+ if ((real_c(globalCell[1]) < Donut_x0 + Ri * sin(shifter / 2.0)) && (real_c(globalCell[1]) > Donut_x0 - Ri)) {
phaseField->get(x, y, z) = 0.0;
} else { phaseField->get(x, y, z) = 1.0; })
}
}
void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t R,
- real_t W = 5)
+ real_t W = 5.0)
{
Vector3< real_t > bubbleMidPoint;
- auto X = blocks->getDomainCellBB().xMax();
- auto Y = blocks->getDomainCellBB().yMax();
- auto Z = blocks->getDomainCellBB().zMax();
+ auto X = real_c(blocks->getDomainCellBB().xMax());
+ auto Y = real_c(blocks->getDomainCellBB().yMax());
+ auto Z = real_c(blocks->getDomainCellBB().zMax());
// 20 percent from the top are filled with the gas phase
- real_t gas_top = Y - Y / 5.0;
+ real_t gas_top = Y - Y / real_c(5.0);
// Diameter of the bubble
- real_t D = R * 2;
+ real_t D = R * real_c(2.0);
// distance in between the bubbles
- int dist = 4;
- auto nx = static_cast< unsigned int >(floor(X / (D + dist * W)));
- auto nz = static_cast< unsigned int >(floor(Z / (D + dist * W)));
+ auto dist = real_c(4.0);
+ auto nx = uint_c(floor(X / (D + dist * W)));
+ auto nz = uint_c(floor(Z / (D + dist * W)));
// fluctuation of the bubble radii
- std::vector< std::vector< real_t > > fluctuation_radius(nx, std::vector< real_t >(nz, 0.0));
- std::vector< std::vector< real_t > > fluctuation_pos(nx, std::vector< real_t >(nz, 0.0));
+ std::vector< std::vector< real_t > > fluctuation_radius(nx, std::vector< real_t >(nz, real_c(0.0)));
+ std::vector< std::vector< real_t > > fluctuation_pos(nx, std::vector< real_t >(nz, real_c(0.0)));
- real_t max_fluctuation_radius = R / 5;
- real_t max_fluctuation_pos = (dist * W) / 3.0;
+ real_t max_fluctuation_radius = R / real_c(5.0);
+ real_t max_fluctuation_pos = (dist * W) / real_c(3.0);
for (unsigned int i = 0; i < nx; ++i)
{
for (unsigned int j = 0; j < nz; ++j)
@@ -118,35 +118,35 @@ void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, Block
for (unsigned int i = 0; i < nx; ++i) {
for (unsigned int j = 0; j < nz; ++j)
{
- bubbleMidPoint[0] = (i + 1) * (D + (dist * W)) - (D + (dist * W)) / 2.0 + fluctuation_pos[i][j];
+ bubbleMidPoint[0] = real_c(i + 1) * (D + (dist * W)) - (D + (dist * W)) / real_c(2.0) + fluctuation_pos[i][j];
bubbleMidPoint[1] = R + W + 4;
- bubbleMidPoint[2] = (j + 1) * (D + (dist * W)) - (D + (dist * W)) / 2.0 + fluctuation_pos[i][j];
+ bubbleMidPoint[2] = real_c(j + 1) * (D + (dist * W)) - (D + (dist * W)) / real_c(2.0) + fluctuation_pos[i][j];
- real_t Ri = sqrt((globalCell[0] - bubbleMidPoint[0]) * (globalCell[0] - bubbleMidPoint[0]) +
- (globalCell[1] - bubbleMidPoint[1]) * (globalCell[1] - bubbleMidPoint[1]) +
- (globalCell[2] - bubbleMidPoint[2]) * (globalCell[2] - bubbleMidPoint[2]));
- if (globalCell[0] >= i * (D + dist * W) && globalCell[0] <= (i + 1) * (D + dist * W) &&
- globalCell[2] >= j * (D + dist * W) && globalCell[2] <= (j + 1) * (D + dist * W))
- phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(2.0 * (Ri - (R - fluctuation_radius[i][j])) / W);
+ real_t Ri = sqrt((real_c(globalCell[0]) - bubbleMidPoint[0]) * (real_c(globalCell[0]) - bubbleMidPoint[0]) +
+ (real_c(globalCell[1]) - bubbleMidPoint[1]) * (real_c(globalCell[1]) - bubbleMidPoint[1]) +
+ (real_c(globalCell[2]) - bubbleMidPoint[2]) * (real_c(globalCell[2]) - bubbleMidPoint[2]));
+ if (real_c(globalCell[0]) >= real_c(i) * (D + dist * W) && real_c(globalCell[0]) <= real_c(i + 1) * (D + dist * W) &&
+ real_c(globalCell[2]) >= real_c(j) * (D + dist * W) && real_c(globalCell[2]) <= real_c(j + 1) * (D + dist * W))
+ phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * real_c(tanh(2.0 * (Ri - (R - fluctuation_radius[i][j])) / W));
- if (globalCell[0] > nx * (D + dist * W)) phaseField->get(x, y, z) = 1.0;
- if (globalCell[2] > nz * (D + dist * W)) phaseField->get(x, y, z) = 1.0;
+ if (real_c(globalCell[0]) > real_c(nx) * (D + dist * W)) phaseField->get(x, y, z) = 1.0;
+ if (real_c(globalCell[2]) > real_c(nz) * (D + dist * W)) phaseField->get(x, y, z) = 1.0;
}
}
- if (globalCell[1] > gas_top) {
- phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(2.0 * (gas_top + 10 - globalCell[1]) / W);
+ if (real_c(globalCell[1]) > gas_top) {
+ phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * real_c(tanh(2.0 * (gas_top + 10 - real_c(globalCell[1])) / W));
})
}
}
void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
- const real_t W = 5, const bool pipe = true)
+ const real_t W = 5.0, const bool pipe = true)
{
- auto X = blocks->getDomainCellBB().xMax();
- auto Z = blocks->getDomainCellBB().zMax();
- auto halfY = (blocks->getDomainCellBB().yMax()) / 2.0;
- real_t perturbation = 0.05;
+ auto X = real_c(blocks->getDomainCellBB().xMax());
+ auto Z = real_c(blocks->getDomainCellBB().zMax());
+ auto halfY = real_c((blocks->getDomainCellBB().yMax())) / real_c(2.0);
+ auto perturbation = real_c(0.05);
if (pipe)
{
@@ -155,10 +155,10 @@ void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, Bloc
auto phaseField = block.getData< PhaseField_T >(phaseFieldID);
WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(
phaseField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
- real_t R = sqrt((globalCell[0] - X / 2) * (globalCell[0] - X / 2) +
- (globalCell[2] - Z / 2) * (globalCell[2] - Z / 2));
- if (R > X) R = X; real_t tmp = perturbation * X * cos((2.0 * math::pi * R) / X);
- phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(((globalCell[1] - halfY) + tmp) / (W / 2.0));)
+ real_t R = sqrt((real_c(globalCell[0]) - X / 2) * (real_c(globalCell[0]) - X / 2) +
+ (real_c(globalCell[2]) - Z / 2) * (real_c(globalCell[2]) - Z / 2));
+ if (R > X) R = X; real_t tmp = perturbation * X * real_c(cos((2.0 * math::pi * R) / X));
+ phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * real_c(tanh(((real_c(globalCell[1]) - halfY) + tmp)) / (W / 2.0));)
}
}
else
@@ -169,8 +169,8 @@ void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, Bloc
WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(
phaseField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
real_t tmp = perturbation * X *
- (cos((2.0 * math::pi * globalCell[0]) / X) + cos((2.0 * math::pi * globalCell[2]) / X));
- phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(((globalCell[1] - halfY) - tmp) / (W / 2.0));)
+ real_c(cos((real_c(2.0) * math::pi * real_c(globalCell[0]))) / X + cos((2.0 * math::pi * real_c(globalCell[2])) / X));
+ phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * real_c(tanh(((real_c(globalCell[1]) - halfY) - tmp)) / (W / 2.0));)
}
}
}
@@ -182,45 +182,45 @@ void initTubeWithCylinder(const shared_ptr< StructuredBlockStorage >& blocks, Bl
{
if (eccentricity_or_pipe_ratio)
{
- auto Mx = blocks->getDomainCellBB().xMax() / 2.0;
- auto Mz = blocks->getDomainCellBB().zMax() / 2.0;
+ auto Mx = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0);
+ auto Mz = real_c(blocks->getDomainCellBB().zMax()) / real_c(2.0);
- auto R_outer = blocks->getDomainCellBB().xMax() / 2.0 + 1.0;
+ auto R_outer = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0) + real_c(1.0);
- real_t const shift = eccentricity * Mx / 2;
+ real_t const shift = eccentricity * Mx / real_c(2.0);
for (auto& block : *blocks)
{
auto flagField = block.template getData< FlagField_T >(flagFieldID);
auto boundaryFlag = flagField->getOrRegisterFlag(boundaryFlagUID);
- WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, Cell globalCell;
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, Cell globalCell;
blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
real_t R1;
- if (globalCell[1] <= start_transition) {
- R1 = sqrt((globalCell[0] - Mx) * (globalCell[0] - Mx) + (globalCell[2] - Mz) * (globalCell[2] - Mz));
- } else if (globalCell[1] > start_transition && globalCell[1] < start_transition + length_transition) {
- real_t tmp = math::pi * (globalCell[1] - start_transition) / (length_transition);
- real_t shift_tmp = shift * 0.5 * (1 - cos(tmp));
- R1 = sqrt((globalCell[0] - Mx - shift_tmp) * (globalCell[0] - Mx - shift_tmp) +
- (globalCell[2] - Mz) * (globalCell[2] - Mz));
+ if (real_c(globalCell[1]) <= start_transition) {
+ R1 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz));
+ } else if (real_c(globalCell[1]) > start_transition && real_c(globalCell[1]) < start_transition + length_transition) {
+ real_t tmp = math::pi * (real_c(globalCell[1]) - start_transition) / (length_transition);
+ real_t shift_tmp = shift * real_c(0.5) * (real_c(1.0) - real_c(cos(tmp)));
+ R1 = sqrt((real_c(globalCell[0]) - Mx - shift_tmp) * (real_c(globalCell[0]) - Mx - shift_tmp) +
+ (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz));
} else {
- R1 = sqrt((globalCell[0] - Mx - shift) * (globalCell[0] - Mx - shift) +
- (globalCell[2] - Mz) * (globalCell[2] - Mz));
+ R1 = sqrt((real_c(globalCell[0]) - Mx - shift) * (real_c(globalCell[0]) - Mx - shift) +
+ (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz));
}
- real_t R2 = sqrt((globalCell[0] - Mx) * (globalCell[0] - Mx) + (globalCell[2] - Mz) * (globalCell[2] - Mz));
+ real_t R2 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz));
if (R1 < R_in) addFlag(flagField->get(x, y, z), boundaryFlag);
if (R2 > R_outer) addFlag(flagField->get(x, y, z), boundaryFlag);)
}
}
else
{
- auto Mx = blocks->getDomainCellBB().xMax() / 2.0;
- auto Mz = blocks->getDomainCellBB().zMax() / 2.0;
-
- auto R_outer = blocks->getDomainCellBB().xMax() / 2.0 + 1.0;
+ auto Mx = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0);
+ auto Mz = real_c(blocks->getDomainCellBB().zMax()) / real_c(2.0);
- real_t const shift = eccentricity * R_in;
+ auto R_outer = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0) + real_c(1.0);
+
+ auto const shift = real_c(eccentricity * R_in);
real_t R_tmp;
for (auto& block : *blocks)
@@ -229,17 +229,17 @@ void initTubeWithCylinder(const shared_ptr< StructuredBlockStorage >& blocks, Bl
auto boundaryFlag = flagField->getOrRegisterFlag(boundaryFlagUID);
WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(
flagField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
- if (globalCell[1] <= start_transition) {
+ if (real_c(globalCell[1]) <= start_transition) {
R_tmp = R_in;
- } else if (globalCell[1] > start_transition && globalCell[1] < start_transition + length_transition) {
- real_t tmp = math::pi * (globalCell[1] - start_transition) / (length_transition);
- real_t shift_tmp = shift * 0.5 * (1 - cos(tmp));
+ } else if (real_c(globalCell[1]) > start_transition && real_c(globalCell[1]) < start_transition + length_transition) {
+ real_t tmp = math::pi * (real_c(globalCell[1]) - start_transition) / (length_transition);
+ real_t shift_tmp = shift * real_c(0.5) * (real_c(1.0) - real_c(cos(tmp)));
R_tmp = R_in + shift_tmp;
} else {
R_tmp = R_in + shift;
}
- real_t R2 = sqrt((globalCell[0] - Mx) * (globalCell[0] - Mx) + (globalCell[2] - Mz) * (globalCell[2] - Mz));
+ real_t R2 = real_c(sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)));
if (R2 < R_tmp) addFlag(flagField->get(x, y, z), boundaryFlag);
if (R2 > R_outer) addFlag(flagField->get(x, y, z), boundaryFlag);)
}
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase.cpp b/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase.cpp
index f2c8f0087..852bd26e5 100644
--- a/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase.cpp
+++ b/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase.cpp
@@ -100,7 +100,7 @@ int main(int argc, char** argv)
const std::string timeStepStrategy = parameters.getParameter< std::string >("timeStepStrategy", "normal");
const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(50));
const real_t remainingTimeLoggerFrequency =
- parameters.getParameter< real_t >("remainingTimeLoggerFrequency", 3.0);
+ parameters.getParameter< real_t >("remainingTimeLoggerFrequency", real_c(3.0));
const uint_t scenario = parameters.getParameter< uint_t >("scenario", uint_c(1));
Vector3< int > gpuBlockSize =
parameters.getParameter< Vector3< int > >("gpuBlockSize", Vector3< int >(128, 1, 1));
@@ -111,8 +111,8 @@ int main(int argc, char** argv)
///////////////////////
// CPU fields
- BlockDataID vel_field = field::addToStorage< VelocityField_T >(blocks, "vel", real_t(0), field::fzyx);
- BlockDataID phase_field = field::addToStorage< PhaseField_T >(blocks, "phase", real_t(0), field::fzyx);
+ BlockDataID vel_field = field::addToStorage< VelocityField_T >(blocks, "vel", real_c(0.0), field::fzyx);
+ BlockDataID phase_field = field::addToStorage< PhaseField_T >(blocks, "phase", real_c(0.0), field::fzyx);
// GPU fields
BlockDataID lb_phase_field_gpu = cuda::addGPUFieldToStorage< cuda::GPUField< real_t > >(
blocks, "lb phase field on GPU", Stencil_phase_T::Size, field::fzyx, 1);
@@ -144,7 +144,7 @@ int main(int argc, char** argv)
{
auto bubbleParameters = config->getOneBlock("Bubble");
const Vector3< real_t > bubbleMidPoint = bubbleParameters.getParameter< Vector3< real_t > >("bubbleMidPoint");
- const real_t bubbleRadius = bubbleParameters.getParameter< real_t >("bubbleRadius", 20.0);
+ const real_t bubbleRadius = bubbleParameters.getParameter< real_t >("bubbleRadius", real_c(20.0));
const bool bubble = bubbleParameters.getParameter< bool >("bubble", true);
initPhaseField_sphere(blocks, phase_field, bubbleRadius, bubbleMidPoint, bubble, interface_thickness);
}
@@ -155,7 +155,7 @@ int main(int argc, char** argv)
else if (scenario == 3)
{
auto bubbleParameters = config->getOneBlock("Bubble");
- const real_t bubbleRadius = bubbleParameters.getParameter< real_t >("bubbleRadius", 20.0);
+ const real_t bubbleRadius = bubbleParameters.getParameter< real_t >("bubbleRadius", real_c(20.0));
init_bubble_field(blocks, phase_field, bubbleRadius);
}
else if (scenario == 4)
@@ -166,9 +166,9 @@ int main(int argc, char** argv)
const real_t diameter = TorusParameters.getParameter< real_t >("Donut_D");
const real_t donutTime = TorusParameters.getParameter< real_t >("DonutTime");
init_Taylor_bubble(blocks, phase_field, diameter, height, donutTime, midpoint);
- center_of_mass[0] = real_t(cellsPerBlock[0]);
- center_of_mass[1] = real_t(midpoint);
- center_of_mass[2] = real_t(cellsPerBlock[2]);
+ center_of_mass[0] = real_c(cellsPerBlock[0]);
+ center_of_mass[1] = real_c(midpoint);
+ center_of_mass[2] = real_c(cellsPerBlock[2]);
}
WALBERLA_LOG_INFO_ON_ROOT("initialization of the phase field done")
@@ -444,7 +444,7 @@ int main(int argc, char** argv)
simTimer.end();
WALBERLA_LOG_INFO_ON_ROOT("Simulation finished")
- auto time = simTimer.last();
+ auto time = real_c(simTimer.last());
auto nrOfCells = real_c(cellsPerBlock[0] * cellsPerBlock[1] * cellsPerBlock[2]);
auto mlupsPerProcess = nrOfCells * real_c(timesteps) / time * 1e-6;
WALBERLA_LOG_RESULT_ON_ROOT("MLUPS per process: " << mlupsPerProcess)
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_codegen.py b/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_codegen.py
index a41a6939b..471a3e9f4 100644
--- a/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_codegen.py
+++ b/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_codegen.py
@@ -1,25 +1,26 @@
import numpy as np
import sympy as sp
-from pystencils import Assignment, fields, TypedSymbol, Target
-from pystencils.astnodes import Block, Conditional
-from pystencils.simp import sympy_cse
+from pystencils import fields, Target
+from pystencils.astnodes import Conditional, Block, SympyAssignment
+from pystencils.node_collection import NodeCollection
+from pystencils.typing import TypedSymbol
+from pystencils.simp.simplifications import sympy_cse
-from lbmpy import LBMConfig, LBStencil, Method, Stencil
+from lbmpy import LBMConfig, LBMOptimisation, LBStencil, Method, Stencil, create_lb_update_rule
from lbmpy.creationfunctions import create_lb_method
from lbmpy.boundaries import NoSlip
from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle
-from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel
from lbmpy.phasefield_allen_cahn.kernel_equations import initializer_kernel_phase_field_lb, \
- initializer_kernel_hydro_lb, interface_tracking_force, hydrodynamic_force, get_collision_assignments_hydro, \
- get_collision_assignments_phase
+ initializer_kernel_hydro_lb, interface_tracking_force, hydrodynamic_force, add_interface_tracking_force, \
+ add_hydrodynamic_force, hydrodynamic_force_assignments
+from lbmpy.phasefield_allen_cahn.parameter_calculation import AllenCahnParameters
import pystencils_walberla
from pystencils_walberla import CodeGeneration, generate_sweep, generate_pack_info_for_field, generate_info_header
from lbmpy_walberla import generate_boundary, generate_lb_pack_info
-
with CodeGeneration() as ctx:
field_type = "float64" if ctx.double_accuracy else "float32"
@@ -29,28 +30,8 @@ with CodeGeneration() as ctx:
stencil_hydro = LBStencil(Stencil.D3Q27)
assert (stencil_phase.D == stencil_hydro.D)
- ########################
- # PARAMETER DEFINITION #
- ########################
-
- density_liquid = sp.Symbol("rho_H")
- density_gas = sp.Symbol("rho_L")
-
- surface_tension = sp.Symbol("sigma")
- mobility = sp.Symbol("mobility")
-
- gravitational_acceleration = sp.Symbol("gravity")
-
- relaxation_time_liquid = sp.Symbol("tau_H")
- relaxation_time_gas = sp.Symbol("tau_L")
-
- # phase-field parameter
- drho3 = (density_liquid - density_gas) / 3
- # interface thickness
- W = sp.Symbol("interface_thickness")
- # coefficients related to surface tension
- beta = 12.0 * (surface_tension / W)
- kappa = 1.5 * surface_tension * W
+ # In the codegeneration skript we only need the symbolic parameters
+ parameters = AllenCahnParameters(0, 0, 0, 0, 0)
########################
# FIELDS #
@@ -75,79 +56,74 @@ with CodeGeneration() as ctx:
########################################
# relaxation rate for interface tracking LB step
- relaxation_rate_allen_cahn = 1.0 / (0.5 + (3.0 * mobility))
+ relaxation_rate_allen_cahn = 1.0 / (0.5 + (3.0 * parameters.symbolic_mobility))
# calculate the relaxation rate for hydrodynamic LB step
- density = density_gas + C.center * (density_liquid - density_gas)
+ rho_L = parameters.symbolic_density_light
+ rho_H = parameters.symbolic_density_heavy
+ density = rho_L + C.center * (rho_H - rho_L)
# force acting on all phases of the model
- body_force = np.array([0, gravitational_acceleration * density, 0])
- # calculation of the relaxation time via viscosities
- # viscosity = viscosity_gas * viscosity_liquid + C.center\
- # * (density_liquid*viscosity_liquid - viscosity_liquid*viscosity_gas)
- # relaxation_time = 3 * viscosity / density + 0.5
-
- relaxation_time = 0.5 + relaxation_time_gas + C.center * (relaxation_time_liquid - relaxation_time_gas)
- # calculate the relaxation time if the phase-field has values over one
- # relaxation_rate = 1.0 / relaxation_time
- relaxation_rate = sp.Symbol("s8")
- relaxation_rate_cutoff = sp.Piecewise((1 / (0.5 + relaxation_time_liquid), C.center > 0.999), # True value
- (1 / relaxation_time, True)) # Else value
+ body_force = [0, 0, 0]
+ body_force[1] = parameters.symbolic_gravitational_acceleration * density
+
+ omega = parameters.omega(C)
###############
# LBM METHODS #
###############
- lbm_config_phase = LBMConfig(stencil=stencil_phase, method=Method.MRT, compressible=True, weighted=True,
- relaxation_rates=[1, 1, 1, 1, 1, 1])
+
+ rates = [0.0]
+ rates += [relaxation_rate_allen_cahn] * stencil_phase.D
+ rates += [1.0] * (stencil_phase.Q - stencil_phase.D - 1)
+
+ lbm_config_phase = LBMConfig(stencil=stencil_phase, method=Method.MRT, compressible=True,
+ delta_equilibrium=False,
+ force=sp.symbols(f"F_:{stencil_phase.D}"), velocity_input=u,
+ weighted=True, relaxation_rates=rates,
+ output={'density': C_tmp}, kernel_type='stream_pull_collide')
method_phase = create_lb_method(lbm_config=lbm_config_phase)
- method_phase.set_first_moment_relaxation_rate(relaxation_rate_allen_cahn)
- lbm_config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.MRT, weighted=True,
- relaxation_rates=[relaxation_rate, 1, 1, 1, 1, 1])
+ lbm_config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.MRT, compressible=False,
+ weighted=True, relaxation_rate=omega,
+ force=sp.symbols(f"F_:{stencil_hydro.D}"),
+ output={'velocity': u}, kernel_type='collide_stream_push')
method_hydro = create_lb_method(lbm_config=lbm_config_hydro)
# create the kernels for the initialization of the g and h field
- h_updates = initializer_kernel_phase_field_lb(h, C, u, method_phase, W)
- g_updates = initializer_kernel_hydro_lb(g, u, method_hydro)
+ h_updates = initializer_kernel_phase_field_lb(method_phase, C, u, h, parameters)
+ g_updates = initializer_kernel_hydro_lb(method_hydro, 1, u, g)
- force_h = [f / 3 for f in interface_tracking_force(C, stencil_phase, W)]
- force_model_h = MultiphaseForceModel(force=force_h)
-
- force_g = hydrodynamic_force(g, C, method_hydro, relaxation_time,
- density_liquid, density_gas, kappa, beta, body_force)
- force_model_g = MultiphaseForceModel(force=force_g, rho=density)
+ force_h = interface_tracking_force(C, stencil_phase, parameters)
+ hydro_force = hydrodynamic_force(g, C, method_hydro, parameters, body_force)
####################
# LBM UPDATE RULES #
####################
- phase_field_LB_step = get_collision_assignments_phase(lb_method=method_phase,
- velocity_input=u,
- output={'density': C_tmp},
- force_model=force_model_h,
- symbolic_fields={"symbolic_field": h,
- "symbolic_temporary_field": h_tmp},
- kernel_type='stream_pull_collide')
+ lbm_optimisation = LBMOptimisation(symbolic_field=h, symbolic_temporary_field=h_tmp)
+ allen_cahn_update_rule = create_lb_update_rule(lbm_config=lbm_config_phase,
+ lbm_optimisation=lbm_optimisation)
+
+ allen_cahn_update_rule = add_interface_tracking_force(allen_cahn_update_rule, force_h)
- phase_field_LB_step = sympy_cse(phase_field_LB_step)
+ allen_cahn_update_rule = sympy_cse(allen_cahn_update_rule)
- phase_field_LB_step = [Conditional(sp.Eq(flag.center(), 2),
- Block(phase_field_LB_step),
- Block([Assignment(C_tmp.center, C.center)]))]
+ allen_cahn_update_rule = NodeCollection([Conditional(sp.Eq(flag.center(), 2),
+ Block(allen_cahn_update_rule.all_assignments),
+ Block([SympyAssignment(C_tmp.center, C.center)]))])
# ---------------------------------------------------------------------------------------------------------
- hydro_LB_step = get_collision_assignments_hydro(lb_method=method_hydro,
- density=density,
- velocity_input=u,
- force_model=force_model_g,
- sub_iterations=2,
- symbolic_fields={"symbolic_field": g,
- "symbolic_temporary_field": g_tmp},
- kernel_type='collide_stream_push')
+ force_Assignments = hydrodynamic_force_assignments(g, u, C, method_hydro, parameters, body_force)
+
+ lbm_optimisation = LBMOptimisation(symbolic_field=g, symbolic_temporary_field=g_tmp)
+ hydro_lb_update_rule = create_lb_update_rule(lbm_config=lbm_config_hydro,
+ lbm_optimisation=lbm_optimisation)
- hydro_LB_step.set_sub_expressions_from_dict({**{relaxation_rate: relaxation_rate_cutoff},
- **hydro_LB_step.subexpressions_dict})
+ hydro_lb_update_rule = add_hydrodynamic_force(hydro_lb_update_rule, force_Assignments, C, g, parameters)
- hydro_LB_step = [Conditional(sp.Eq(flag.center(), 2), Block(hydro_LB_step))]
+ hydro_lb_update_rule = NodeCollection([Conditional(sp.Eq(flag.center(), 2),
+ Block(hydro_lb_update_rule.all_assignments))])
- contact_angle = ContactAngle(contact_angle_in_degrees, W)
+ contact_angle = ContactAngle(contact_angle_in_degrees, parameters.symbolic_interface_thickness,
+ data_type=field_type)
###################
# GENERATE SWEEPS #
@@ -178,14 +154,14 @@ with CodeGeneration() as ctx:
generate_sweep(ctx, 'initialize_phase_field_distributions', h_updates, target=Target.GPU)
generate_sweep(ctx, 'initialize_velocity_based_distributions', g_updates, target=Target.GPU)
- generate_sweep(ctx, 'phase_field_LB_step', phase_field_LB_step,
+ generate_sweep(ctx, 'phase_field_LB_step', allen_cahn_update_rule,
field_swaps=[(h, h_tmp), (C, C_tmp)],
target=Target.GPU,
gpu_indexing_params=sweep_params,
varying_parameters=vp)
generate_boundary(ctx, 'phase_field_LB_NoSlip', NoSlip(), method_phase, target=Target.GPU, streaming_pattern='pull')
- generate_sweep(ctx, 'hydro_LB_step', hydro_LB_step,
+ generate_sweep(ctx, 'hydro_LB_step', hydro_lb_update_rule,
field_swaps=[(g, g_tmp)],
target=Target.GPU,
gpu_indexing_params=sweep_params,
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/util.cpp b/apps/showcases/PhaseFieldAllenCahn/GPU/util.cpp
index 9912edb76..f6694910c 100644
--- a/apps/showcases/PhaseFieldAllenCahn/GPU/util.cpp
+++ b/apps/showcases/PhaseFieldAllenCahn/GPU/util.cpp
@@ -46,13 +46,13 @@ void calc_total_velocity(const shared_ptr< StructuredBlockStorage >& blocks, std
phaseField,
auto fluidFlag = flagField->getFlag(fluidFlagUID); if (flagField->get(x, y, z) == fluidFlag) {
- if (phaseField->get(x, y, z) < 0.5)
+ if (phaseField->get(x, y, z) < real_c(0.5))
{
real_t invC = 1 - phaseField->get(x, y, z);
real_t U = velField->get(x, y, z, 0);
real_t V = velField->get(x, y, z, 1);
real_t W = velField->get(x, y, z, 2);
- total_velocity[0] += sqrt((U * invC) * (U * invC) + (V * invC) * (V * invC) + (W * invC) * (W * invC));
+ total_velocity[0] += real_c(sqrt((U * invC) * (U * invC) + (V * invC) * (V * invC) + (W * invC) * (W * invC)));
total_velocity[1] += U * invC;
total_velocity[2] += V * invC;
total_velocity[3] += W * invC;
@@ -69,11 +69,11 @@ void flood_fill(PhaseField_T& phaseField, VelocityField_T& velocityField, CellIn
Field< bool, 1 > visit(phaseField.xSize(), phaseField.ySize(), phaseField.zSize(), false, field::fzyx);
using namespace stencil;
- volume = 0;
- nrOfCells = 0;
- center_of_mass[0] = 0.0;
- center_of_mass[1] = 0.0;
- center_of_mass[2] = 0.0;
+ volume = real_c(0.0);
+ nrOfCells = uint_c(0);
+ center_of_mass[0] = real_c( 0.0);
+ center_of_mass[1] = real_c(0.0);
+ center_of_mass[2] = real_c(0.0);
while (phaseField.get(startCell) > 0.5 && startCell.x() > 0)
--startCell.x();
@@ -92,14 +92,14 @@ void flood_fill(PhaseField_T& phaseField, VelocityField_T& velocityField, CellIn
nrOfCells++;
volume += invC;
- total_velocity[0] += sqrt((v_U * invC) * (v_U * invC) + (v_V * invC) * (v_V * invC) + (v_W * invC) * (v_W * invC));
+ total_velocity[0] += real_c(sqrt((v_U * invC) * (v_U * invC) + (v_V * invC) * (v_V * invC) + (v_W * invC) * (v_W * invC)));
total_velocity[1] += v_U * invC;
total_velocity[2] += v_V * invC;
total_velocity[3] += v_W * invC;
- center_of_mass[0] += (startCell.x() + boundingBox.xMin());
- center_of_mass[1] += (startCell.y() + boundingBox.yMin());
- center_of_mass[2] += (startCell.z() + boundingBox.xMin());
+ center_of_mass[0] += real_c(startCell.x() + boundingBox.xMin());
+ center_of_mass[1] += real_c(startCell.y() + boundingBox.yMin());
+ center_of_mass[2] += real_c(startCell.z() + boundingBox.xMin());
const int DIRS[6] = { N, S, E, W, T, B };
@@ -126,14 +126,14 @@ void flood_fill(PhaseField_T& phaseField, VelocityField_T& velocityField, CellIn
volume += invC;
total_velocity[0] +=
- sqrt((v_U * invC) * (v_U * invC) + (v_V * invC) * (v_V * invC) + (v_W * invC) * (v_W * invC));
+ real_c(sqrt((v_U * invC) * (v_U * invC) + (v_V * invC) * (v_V * invC) + (v_W * invC) * (v_W * invC)));
total_velocity[1] += v_U * invC;
total_velocity[2] += v_V * invC;
total_velocity[3] += v_W * invC;
- center_of_mass[0] += (neighborCell.x() + boundingBox.xMin());
- center_of_mass[1] += (neighborCell.y() + boundingBox.yMin());
- center_of_mass[2] += (neighborCell.z() + boundingBox.xMin());
+ center_of_mass[0] += real_c(neighborCell.x() + boundingBox.xMin());
+ center_of_mass[1] += real_c(neighborCell.y() + boundingBox.yMin());
+ center_of_mass[2] += real_c(neighborCell.z() + boundingBox.xMin());
visit.get(neighborCell) = true;
cellQueue.push(neighborCell);
diff --git a/apps/showcases/PorousMedia/PorousMediaCumulantLBMKernel.py b/apps/showcases/PorousMedia/PorousMediaCumulantLBMKernel.py
index 96cc2f0aa..8bee23b8e 100644
--- a/apps/showcases/PorousMedia/PorousMediaCumulantLBMKernel.py
+++ b/apps/showcases/PorousMedia/PorousMediaCumulantLBMKernel.py
@@ -1,39 +1,29 @@
import sympy as sp
import pystencils as ps
-from lbmpy.stencils import get_stencil
-from lbmpy.creationfunctions import create_lb_collision_rule, create_lb_update_rule
+from lbmpy import Stencil, LBStencil
+from lbmpy.creationfunctions import LBMConfig, LBMOptimisation, Method, create_lb_collision_rule
from lbmpy_walberla import RefinementScaling, generate_lattice_model
from pystencils_walberla import CodeGeneration, generate_sweep, generate_pack_info_from_kernel
-stencil = get_stencil("D3Q19", ordering="walberla")
-omega = sp.Symbol('omega')
-layout = 'fzyx'
-pdf_field = ps.fields("pdfs(19): [3D]", layout=layout)
-
-optimizations = {'cse_global': True,
- 'cse_pdfs': False,
- 'split': True,
- 'symbolic_field': pdf_field,
- 'field_layout': layout,
- 'pre_simplification': True
- }
-
-params = {'stencil': stencil,
- 'method': 'cumulant',
- 'relaxation_rate': omega,
- 'equilibrium_order': 4,
- 'maxwellian_moments': True,
- 'compressible': True,
- }
-
-collision_rule = create_lb_collision_rule(optimization=optimizations, **params)
-update_rule = create_lb_update_rule(collision_rule=collision_rule, optimization=optimizations)
-
-refinement_scaling = RefinementScaling()
-refinement_scaling.add_standard_relaxation_rate_scaling(omega)
with CodeGeneration() as ctx:
- generate_lattice_model(ctx, "LbCodeGen_LatticeModel", collision_rule, refinement_scaling=refinement_scaling,
- field_layout=layout)
+ data_type = "float64" if ctx.double_accuracy else "float32"
+ stencil = LBStencil(Stencil.D3Q19)
+ omega = sp.Symbol('omega')
+ layout = 'fzyx'
+ pdf_field = ps.fields(f"pdfs({stencil.Q}): {data_type}[3D]", layout=layout)
+
+ lbm_config = LBMConfig(stencil=stencil, method=Method.CUMULANT, relaxation_rate=omega, compressible=True)
+ lbm_opt = LBMOptimisation(cse_global=True, cse_pdfs=False, split=False,
+ symbolic_field=pdf_field, field_layout=layout, pre_simplification=True)
+
+ collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+
+ refinement_scaling = RefinementScaling()
+ refinement_scaling.add_standard_relaxation_rate_scaling(omega)
+
+ generate_lattice_model(ctx, "LbCodeGen_LatticeModel", collision_rule, refinement_scaling=refinement_scaling,
+ field_layout=layout)
+
diff --git a/apps/tutorials/codegen/02_LBMLatticeModelGeneration.cpp b/apps/tutorials/codegen/02_LBMLatticeModelGeneration.cpp
index 934ab64ff..2a83abba2 100644
--- a/apps/tutorials/codegen/02_LBMLatticeModelGeneration.cpp
+++ b/apps/tutorials/codegen/02_LBMLatticeModelGeneration.cpp
@@ -86,7 +86,7 @@ struct ShearFlowInit
std::vector< real_t > operator()(const Cell& globalCell)
{
std::vector< real_t > densityAndVelocity = exprInitFunc_(globalCell);
- real_t yPerturbation = noiseMagnitude_ * rng_();
+ auto yPerturbation = real_c(noiseMagnitude_ * rng_());
densityAndVelocity[2] += yPerturbation;
return densityAndVelocity;
@@ -115,7 +115,7 @@ int main(int argc, char** argv)
const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(10));
const real_t omega = parameters.getParameter< real_t >("omega", real_c(1.8));
const double remainingTimeLoggerFrequency =
- parameters.getParameter< double >("remainingTimeLoggerFrequency", 3.0); // in seconds
+ parameters.getParameter< double >("remainingTimeLoggerFrequency", real_c(3.0)); // in seconds
///////////////////
/// Field Setup ///
diff --git a/apps/tutorials/codegen/02_LBMLatticeModelGeneration.py b/apps/tutorials/codegen/02_LBMLatticeModelGeneration.py
index 224ab742a..494c47bf2 100644
--- a/apps/tutorials/codegen/02_LBMLatticeModelGeneration.py
+++ b/apps/tutorials/codegen/02_LBMLatticeModelGeneration.py
@@ -21,7 +21,8 @@ lbm_opt = LBMOptimisation(cse_global=True, field_layout=layout)
# SRT Method Definition
# ===========================
-lbm_config = LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=omega)
+lbm_config = LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=omega,
+ zero_centered=False, delta_equilibrium=False)
srt_collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
srt_update_rule = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
diff --git a/apps/tutorials/codegen/03_AdvancedLBMCodegen.cpp b/apps/tutorials/codegen/03_AdvancedLBMCodegen.cpp
index 492d40b8b..5b7790c7d 100644
--- a/apps/tutorials/codegen/03_AdvancedLBMCodegen.cpp
+++ b/apps/tutorials/codegen/03_AdvancedLBMCodegen.cpp
@@ -87,7 +87,7 @@ void initShearFlowVelocityField(const shared_ptr< StructuredBlockForest >& block
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());
+ auto n_y = real_c(blocks->getNumberOfYCells());
for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
{
@@ -98,7 +98,7 @@ void initShearFlowVelocityField(const shared_ptr< StructuredBlockForest >& block
Cell globalCell(cellIt.cell());
blocks->transformBlockLocalToGlobalCell(globalCell, *blockIt);
- real_t relative_y = real_c(globalCell.y()) / n_y;
+ auto relative_y = real_c(globalCell.y()) / n_y;
u->get(cellIt.cell(), 0) = relative_y < 0.3 || relative_y > 0.7 ? velocityMagnitude : -velocityMagnitude;
@@ -129,7 +129,7 @@ int main(int argc, char** argv)
const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(10));
const real_t omega = parameters.getParameter< real_t >("omega", real_c(1.8));
const double remainingTimeLoggerFrequency =
- parameters.getParameter< double >("remainingTimeLoggerFrequency", 3.0); // in seconds
+ parameters.getParameter< double >("remainingTimeLoggerFrequency", real_c(3.0)); // in seconds
const uint_t VTKwriteFrequency = parameters.getParameter< uint_t >("VTKwriteFrequency", 1000);
////////////////////////////////////
diff --git a/apps/tutorials/codegen/03_AdvancedLBMCodegen.py b/apps/tutorials/codegen/03_AdvancedLBMCodegen.py
index c94eb81ba..b139c9999 100644
--- a/apps/tutorials/codegen/03_AdvancedLBMCodegen.py
+++ b/apps/tutorials/codegen/03_AdvancedLBMCodegen.py
@@ -10,59 +10,56 @@ from lbmpy.boundaries import NoSlip
from pystencils_walberla import CodeGeneration, generate_sweep, generate_pack_info_from_kernel
from lbmpy_walberla import generate_boundary
-
-# ========================
-# General Parameters
-# ========================
-
-stencil = LBStencil(Stencil.D2Q9)
-omega = sp.Symbol('omega')
-layout = 'fzyx'
-
-# PDF Fields
-pdfs, pdfs_tmp = ps.fields(f'pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): [2D]', layout=layout)
-
-# Velocity Output Field
-velocity = ps.fields(f"velocity({stencil.D}): [2D]", layout=layout)
-output = {'velocity': velocity}
-
-# LBM Optimisation
-lbm_opt = LBMOptimisation(cse_global=True,
- symbolic_field=pdfs,
- symbolic_temporary_field=pdfs_tmp,
- field_layout=layout)
-
-
-# ==================
-# Method Setup
-# ==================
-
-lbm_config = LBMConfig(stencil=stencil,
- method=Method.CUMULANT,
- relaxation_rate=omega,
- compressible=True,
- output=output)
-
-lbm_update_rule = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
-
-lbm_method = lbm_update_rule.method
-
-# ========================
-# PDF Initialization
-# ========================
-
-initial_rho = sp.Symbol('rho_0')
-
-pdfs_setter = macroscopic_values_setter(lbm_method,
- initial_rho,
- velocity.center_vector,
- pdfs.center_vector)
-
# =====================
# Code Generation
# =====================
with CodeGeneration() as ctx:
+ # ========================
+ # General Parameters
+ # ========================
+ data_type = "float64" if ctx.double_accuracy else "float32"
+ stencil = LBStencil(Stencil.D2Q9)
+ omega = sp.Symbol('omega')
+ layout = 'fzyx'
+
+ # PDF Fields
+ pdfs, pdfs_tmp = ps.fields(f'pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): {data_type}[2D]', layout=layout)
+
+ # Velocity Output Field
+ velocity = ps.fields(f"velocity({stencil.D}): {data_type}[2D]", layout=layout)
+ output = {'velocity': velocity}
+
+ # LBM Optimisation
+ lbm_opt = LBMOptimisation(cse_global=True,
+ symbolic_field=pdfs,
+ symbolic_temporary_field=pdfs_tmp,
+ field_layout=layout)
+
+ # ==================
+ # Method Setup
+ # ==================
+
+ lbm_config = LBMConfig(stencil=stencil,
+ method=Method.CUMULANT,
+ relaxation_rate=omega,
+ compressible=True,
+ output=output)
+
+ lbm_update_rule = create_lb_update_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+ lbm_method = lbm_update_rule.method
+
+ # ========================
+ # PDF Initialization
+ # ========================
+
+ initial_rho = sp.Symbol('rho_0')
+
+ pdfs_setter = macroscopic_values_setter(lbm_method,
+ initial_rho,
+ velocity.center_vector,
+ pdfs.center_vector)
+
if ctx.cuda:
target = ps.Target.GPU
else:
diff --git a/apps/tutorials/codegen/HeatEquationKernel.py b/apps/tutorials/codegen/HeatEquationKernel.py
index 2aefdba6d..024940a38 100644
--- a/apps/tutorials/codegen/HeatEquationKernel.py
+++ b/apps/tutorials/codegen/HeatEquationKernel.py
@@ -2,24 +2,26 @@ import sympy as sp
import pystencils as ps
from pystencils_walberla import CodeGeneration, generate_sweep
-u, u_tmp = ps.fields("u, u_tmp: [2D]", layout='fzyx')
-kappa = sp.Symbol("kappa")
-dx = sp.Symbol("dx")
-dt = sp.Symbol("dt")
-heat_pde = ps.fd.transient(u) - kappa * (ps.fd.diff(u, 0, 0) + ps.fd.diff(u, 1, 1))
+with CodeGeneration() as ctx:
+ data_type = "float64" if ctx.double_accuracy else "float32"
-discretize = ps.fd.Discretization2ndOrder(dx=dx, dt=dt)
-heat_pde_discretized = discretize(heat_pde)
-heat_pde_discretized = heat_pde_discretized.args[1] + heat_pde_discretized.args[0].simplify()
+ u, u_tmp = ps.fields(f"u, u_tmp: {data_type}[2D]", layout='fzyx')
+ kappa = sp.Symbol("kappa")
+ dx = sp.Symbol("dx")
+ dt = sp.Symbol("dt")
+ heat_pde = ps.fd.transient(u) - kappa * (ps.fd.diff(u, 0, 0) + ps.fd.diff(u, 1, 1))
+ discretize = ps.fd.Discretization2ndOrder(dx=dx, dt=dt)
+ heat_pde_discretized = discretize(heat_pde)
+ heat_pde_discretized = heat_pde_discretized.args[1] + heat_pde_discretized.args[0].simplify()
-@ps.kernel
-def update():
- u_tmp.center @= heat_pde_discretized
+ @ps.kernel
+ def update():
+ u_tmp.center @= heat_pde_discretized
-ac = ps.AssignmentCollection(update)
-ac = ps.simp.simplifications.add_subexpressions_for_divisions(ac)
-with CodeGeneration() as ctx:
+ ac = ps.AssignmentCollection(update)
+ ac = ps.simp.simplifications.add_subexpressions_for_divisions(ac)
+
generate_sweep(ctx, 'HeatEquationKernel', ac)
diff --git a/python/lbmpy_walberla/additional_data_handler.py b/python/lbmpy_walberla/additional_data_handler.py
index 214347451..32f8c0ee2 100644
--- a/python/lbmpy_walberla/additional_data_handler.py
+++ b/python/lbmpy_walberla/additional_data_handler.py
@@ -1,5 +1,5 @@
from pystencils import Target
-from pystencils.stencil import inverse_direction, offset_to_direction_string
+from pystencils.stencil import inverse_direction
from lbmpy.advanced_streaming import AccessPdfValues, numeric_offsets, numeric_index
diff --git a/python/lbmpy_walberla/packinfo.py b/python/lbmpy_walberla/packinfo.py
index 27f26b05b..796ccfd98 100644
--- a/python/lbmpy_walberla/packinfo.py
+++ b/python/lbmpy_walberla/packinfo.py
@@ -2,9 +2,11 @@ from collections import defaultdict
from lbmpy import LBStencil, Stencil
from lbmpy.advanced_streaming.utility import Timestep, get_accessor, get_timesteps
from lbmpy.advanced_streaming.communication import _extend_dir
+
+from pystencils import Assignment, Field, Target
from pystencils.stencil import inverse_direction
+
from pystencils_walberla.codegen import comm_directions, generate_pack_info
-from pystencils import Assignment, Field, Target
def generate_lb_pack_info(generation_context,
diff --git a/python/lbmpy_walberla/sparse.py b/python/lbmpy_walberla/sparse.py
index 95e45f64e..fc1482b39 100644
--- a/python/lbmpy_walberla/sparse.py
+++ b/python/lbmpy_walberla/sparse.py
@@ -3,7 +3,9 @@ import numpy as np
from lbmpy.sparse import (
create_lb_update_rule_sparse, create_macroscopic_value_getter_sparse,
create_macroscopic_value_setter_sparse, create_symbolic_list)
-from pystencils import TypedSymbol, create_kernel
+
+from pystencils import create_kernel
+from pystencils.typing import TypedSymbol
class ListLbGenerator:
diff --git a/python/lbmpy_walberla/templates/LatticeModel.tmpl.cpp b/python/lbmpy_walberla/templates/LatticeModel.tmpl.cpp
index 52528b172..dd50337e1 100644
--- a/python/lbmpy_walberla/templates/LatticeModel.tmpl.cpp
+++ b/python/lbmpy_walberla/templates/LatticeModel.tmpl.cpp
@@ -58,8 +58,8 @@ namespace {{namespace}} {
{{stream_kernel|generate_definition(target)}}
-const real_t {{class_name}}::w[{{Q}}] = { {{weights}} };
-const real_t {{class_name}}::wInv[{{Q}}] = { {{inverse_weights}} };
+const {{dtype}} {{class_name}}::w[{{Q}}] = { {{weights}} };
+const {{dtype}} {{class_name}}::wInv[{{Q}}] = { {{inverse_weights}} };
void {{class_name}}::Sweep::streamCollide( IBlock * block, const uint_t numberOfGhostLayersToInclude )
{
diff --git a/python/lbmpy_walberla/templates/LatticeModel.tmpl.h b/python/lbmpy_walberla/templates/LatticeModel.tmpl.h
index 39ae1afd9..d1e299578 100644
--- a/python/lbmpy_walberla/templates/LatticeModel.tmpl.h
+++ b/python/lbmpy_walberla/templates/LatticeModel.tmpl.h
@@ -20,6 +20,7 @@
#pragma once
#include "core/DataTypes.h"
#include "core/logging/Logging.h"
+#include "core/math/Matrix3.h"
#include "field/GhostLayerField.h"
#include "field/SwapableCompare.h"
@@ -98,8 +99,8 @@ class {{class_name}}
public:
typedef stencil::{{stencil_name}} Stencil;
typedef stencil::{{communication_stencil_name}} CommunicationStencil;
- static const real_t w[{{Q}}];
- static const real_t wInv[{{Q}}];
+ static const {{dtype}} w[{{Q}}];
+ static const {{dtype}} wInv[{{Q}}];
static const bool compressible = {% if compressible %}true{% else %}false{% endif %};
@@ -158,11 +159,11 @@ private:
if( targetLevel != currentLevel )
{
- real_t level_scale_factor(1);
+ {{dtype}} level_scale_factor(1);
if( currentLevel < targetLevel )
- level_scale_factor = real_c( uint_t(1) << ( targetLevel - currentLevel ) );
+ level_scale_factor = {{dtype}}( uint_t(1) << ( targetLevel - currentLevel ) );
else // currentLevel > targetLevel
- level_scale_factor = real_t(1) / real_c( uint_t(1) << ( currentLevel - targetLevel ) );
+ level_scale_factor = {{dtype}}(1) / {{dtype}}( uint_t(1) << ( currentLevel - targetLevel ) );
{% for scalingType, name, expression in refinement_scaling_info -%}
{% if scalingType == 'normal' %}
@@ -225,44 +226,44 @@ class EquilibriumDistribution< {{class_name}}, void>
public:
typedef typename {{class_name}}::Stencil Stencil;
- static real_t get( const stencil::Direction direction,
- const Vector3< real_t > & u = Vector3< real_t >( real_t(0.0) ),
- real_t rho = real_t(1.0) )
+ static {{dtype}} get( const stencil::Direction direction,
+ const Vector3< {{dtype}} > & u = Vector3< {{dtype}} >( {{dtype}}(0.0) ),
+ {{dtype}} rho = {{dtype}}(1.0) )
{
{% if not compressible %}
- rho -= real_t(1.0);
+ rho -= {{dtype}}(1.0);
{% endif %}
{{equilibrium_from_direction}}
}
- static real_t getSymmetricPart( const stencil::Direction direction,
- const Vector3<real_t> & u = Vector3< real_t >(real_t(0.0)),
- real_t rho = real_t(1.0) )
+ static {{dtype}} getSymmetricPart( const stencil::Direction direction,
+ const Vector3<{{dtype}}> & u = Vector3< {{dtype}} >({{dtype}}(0.0)),
+ {{dtype}} rho = {{dtype}}(1.0) )
{
{% if not compressible %}
- rho -= real_t(1.0);
+ rho -= {{dtype}}(1.0);
{% endif %}
{{symmetric_equilibrium_from_direction}}
}
- static real_t getAsymmetricPart( const stencil::Direction direction,
- const Vector3< real_t > & u = Vector3<real_t>( real_t(0.0) ),
- real_t rho = real_t(1.0) )
+ static {{dtype}} getAsymmetricPart( const stencil::Direction direction,
+ const Vector3< {{dtype}} > & u = Vector3<{{dtype}}>( {{dtype}}(0.0) ),
+ {{dtype}} rho = {{dtype}}(1.0) )
{
{% if not compressible %}
- rho -= real_t(1.0);
+ rho -= {{dtype}}(1.0);
{% endif %}
{{asymmetric_equilibrium_from_direction}}
}
- static std::vector< real_t > get( const Vector3< real_t > & u = Vector3<real_t>( real_t(0.0) ),
- real_t rho = real_t(1.0) )
+ static std::vector< {{dtype}} > get( const Vector3< {{dtype}} > & u = Vector3<{{dtype}}>( {{dtype}}(0.0) ),
+ {{dtype}} rho = {{dtype}}(1.0) )
{
{% if not compressible %}
- rho -= real_t(1.0);
+ rho -= {{dtype}}(1.0);
{% endif %}
- std::vector< real_t > equilibrium( Stencil::Size );
+ std::vector< {{dtype}} > equilibrium( Stencil::Size );
for( auto d = Stencil::begin(); d != Stencil::end(); ++d )
{
equilibrium[d.toIdx()] = get(*d, u, rho);
@@ -278,25 +279,25 @@ template<>
struct AdaptVelocityToForce<{{class_name}}, void>
{
template< typename FieldPtrOrIterator >
- static Vector3<real_t> get( FieldPtrOrIterator & it, const {{class_name}} & lm,
- const Vector3< real_t > & velocity, const real_t rho )
+ static Vector3<{{dtype}}> get( FieldPtrOrIterator & it, const {{class_name}} & lm,
+ const Vector3< {{dtype}} > & velocity, const {{dtype}} rho )
{
auto x = it.x();
auto y = it.y();
auto z = it.z();
{% if macroscopic_velocity_shift %}
- return velocity - Vector3<real_t>({{macroscopic_velocity_shift | join(",") }} {% if D == 2 %}, real_t(0.0) {%endif %} );
+ return velocity - Vector3<{{dtype}}>({{macroscopic_velocity_shift | join(",") }} {% if D == 2 %}, {{dtype}}(0.0) {%endif %} );
{% else %}
return velocity;
{% endif %}
}
- static Vector3<real_t> get( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const {{class_name}} & lm,
- const Vector3< real_t > & velocity, const real_t rho )
+ static Vector3<{{dtype}}> get( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const {{class_name}} & lm,
+ const Vector3< {{dtype}} > & velocity, const {{dtype}} rho )
{
{% if macroscopic_velocity_shift %}
- return velocity - Vector3<real_t>({{macroscopic_velocity_shift | join(",") }} {% if D == 2 %}, real_t(0.0) {%endif %} );
+ return velocity - Vector3<{{dtype}}>({{macroscopic_velocity_shift | join(",") }} {% if D == 2 %}, {{dtype}}(0.0) {%endif %} );
{% else %}
return velocity;
{% endif %}
@@ -312,10 +313,10 @@ struct Equilibrium< {{class_name}}, void >
template< typename FieldPtrOrIterator >
static void set( FieldPtrOrIterator & it,
- const Vector3< real_t > & u = Vector3< real_t >( real_t(0.0) ), real_t rho = real_t(1.0) )
+ const Vector3< {{dtype}} > & u = Vector3< {{dtype}} >( {{dtype}}(0.0) ), {{dtype}} rho = {{dtype}}(1.0) )
{
{%if not compressible %}
- rho -= real_t(1.0);
+ rho -= {{dtype}}(1.0);
{%endif %}
{% for eqTerm in equilibrium -%}
@@ -325,13 +326,13 @@ struct Equilibrium< {{class_name}}, void >
template< typename PdfField_T >
static void set( PdfField_T & pdf, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
- const Vector3< real_t > & u = Vector3< real_t >( real_t(0.0) ), real_t rho = real_t(1.0) )
+ const Vector3< {{dtype}} > & u = Vector3< {{dtype}} >( {{dtype}}(0.0) ), {{dtype}} rho = {{dtype}}(1.0) )
{
{%if not compressible %}
- rho -= real_t(1.0);
+ rho -= {{dtype}}(1.0);
{%endif %}
- real_t & xyz0 = pdf(x,y,z,0);
+ {{dtype}} & xyz0 = pdf(x,y,z,0);
{% for eqTerm in equilibrium -%}
pdf.getF( &xyz0, {{loop.index0 }})= {{eqTerm}};
{% endfor -%}
@@ -343,22 +344,22 @@ template<>
struct Density<{{class_name}}, void>
{
template< typename FieldPtrOrIterator >
- static inline real_t get( const {{class_name}} & , const FieldPtrOrIterator & it )
+ static inline {{dtype}} get( const {{class_name}} & , const FieldPtrOrIterator & it )
{
{% for i in range(Q) -%}
- const real_t f_{{i}} = it[{{i}}];
+ const {{dtype}} f_{{i}} = it[{{i}}];
{% endfor -%}
{{density_getters | indent(8)}}
return rho;
}
template< typename PdfField_T >
- static inline real_t get( const {{class_name}} & ,
+ static inline {{dtype}} get( const {{class_name}} & ,
const PdfField_T & pdf, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
{
- const real_t & xyz0 = pdf(x,y,z,0);
+ const {{dtype}} & xyz0 = pdf(x,y,z,0);
{% for i in range(Q) -%}
- const real_t f_{{i}} = pdf.getF( &xyz0, {{i}});
+ const {{dtype}} f_{{i}} = pdf.getF( &xyz0, {{i}});
{% endfor -%}
{{density_getters | indent(8)}}
return rho;
@@ -371,7 +372,7 @@ struct DensityAndVelocity<{{class_name}}>
{
template< typename FieldPtrOrIterator >
static void set( FieldPtrOrIterator & it, const {{class_name}} & lm,
- const Vector3< real_t > & u = Vector3< real_t >( real_t(0.0) ), const real_t rho_in = real_t(1.0) )
+ const Vector3< {{dtype}} > & u = Vector3< {{dtype}} >( {{dtype}}(0.0) ), const {{dtype}} rho_in = {{dtype}}(1.0) )
{
auto x = it.x();
auto y = it.y();
@@ -379,22 +380,22 @@ struct DensityAndVelocity<{{class_name}}>
{{density_velocity_setter_macroscopic_values | indent(8)}}
{% if D == 2 -%}
- const real_t u_2(0.0);
+ const {{dtype}} u_2(0.0);
{% endif %}
- Equilibrium<{{class_name}}>::set(it, Vector3<real_t>(u_0, u_1, u_2), rho{%if not compressible %} + real_t(1) {%endif%});
+ Equilibrium<{{class_name}}>::set(it, Vector3<{{dtype}}>(u_0, u_1, u_2), rho{%if not compressible %} + {{dtype}}(1) {%endif%});
}
template< typename PdfField_T >
static void set( PdfField_T & pdf, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const {{class_name}} & lm,
- const Vector3< real_t > & u = Vector3< real_t >( real_t(0.0) ), const real_t rho_in = real_t(1.0) )
+ const Vector3< {{dtype}} > & u = Vector3< {{dtype}} >( {{dtype}}(0.0) ), const {{dtype}} rho_in = {{dtype}}(1.0) )
{
{{density_velocity_setter_macroscopic_values | indent(8)}}
{% if D == 2 -%}
- const real_t u_2(0.0);
+ const {{dtype}} u_2(0.0);
{% endif %}
- Equilibrium<{{class_name}}>::set(pdf, x, y, z, Vector3<real_t>(u_0, u_1, u_2), rho {%if not compressible %} + real_t(1) {%endif%});
+ Equilibrium<{{class_name}}>::set(pdf, x, y, z, Vector3<{{dtype}}>(u_0, u_1, u_2), rho {%if not compressible %} + {{dtype}}(1) {%endif%});
}
};
@@ -404,7 +405,7 @@ struct DensityAndVelocityRange<{{class_name}}, FieldIteratorXYZ>
{
static void set( FieldIteratorXYZ & begin, const FieldIteratorXYZ & end, const {{class_name}} & lm,
- const Vector3< real_t > & u = Vector3< real_t >( real_t(0.0) ), const real_t rho_in = real_t(1.0) )
+ const Vector3< {{dtype}} > & u = Vector3< {{dtype}} >( {{dtype}}(0.0) ), const {{dtype}} rho_in = {{dtype}}(1.0) )
{
for( auto cellIt = begin; cellIt != end; ++cellIt )
{
@@ -413,10 +414,10 @@ struct DensityAndVelocityRange<{{class_name}}, FieldIteratorXYZ>
const auto z = cellIt.z();
{{density_velocity_setter_macroscopic_values | indent(12)}}
{% if D == 2 -%}
- const real_t u_2(0.0);
+ const {{dtype}} u_2(0.0);
{% endif %}
- Equilibrium<{{class_name}}>::set(cellIt, Vector3<real_t>(u_0, u_1, u_2), rho{%if not compressible %} + real_t(1) {%endif%});
+ Equilibrium<{{class_name}}>::set(cellIt, Vector3<{{dtype}}>(u_0, u_1, u_2), rho{%if not compressible %} + {{dtype}}(1) {%endif%});
}
}
};
@@ -427,7 +428,7 @@ template<>
struct DensityAndMomentumDensity<{{class_name}}>
{
template< typename FieldPtrOrIterator >
- static real_t get( Vector3< real_t > & momentumDensity, const {{class_name}} & lm,
+ static {{dtype}} get( Vector3< {{dtype}} > & momentumDensity, const {{class_name}} & lm,
const FieldPtrOrIterator & it )
{
const auto x = it.x();
@@ -435,7 +436,7 @@ struct DensityAndMomentumDensity<{{class_name}}>
const auto z = it.z();
{% for i in range(Q) -%}
- const real_t f_{{i}} = it[{{i}}];
+ const {{dtype}} f_{{i}} = it[{{i}}];
{% endfor -%}
{{momentum_density_getter | indent(8) }}
@@ -446,12 +447,12 @@ struct DensityAndMomentumDensity<{{class_name}}>
}
template< typename PdfField_T >
- static real_t get( Vector3< real_t > & momentumDensity, const {{class_name}} & lm, const PdfField_T & pdf,
+ static {{dtype}} get( Vector3< {{dtype}} > & momentumDensity, const {{class_name}} & lm, const PdfField_T & pdf,
const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
{
- const real_t & xyz0 = pdf(x,y,z,0);
+ const {{dtype}} & xyz0 = pdf(x,y,z,0);
{% for i in range(Q) -%}
- const real_t f_{{i}} = pdf.getF( &xyz0, {{i}});
+ const {{dtype}} f_{{i}} = pdf.getF( &xyz0, {{i}});
{% endfor -%}
{{momentum_density_getter | indent(8) }}
@@ -467,14 +468,14 @@ template<>
struct MomentumDensity< {{class_name}}>
{
template< typename FieldPtrOrIterator >
- static void get( Vector3< real_t > & momentumDensity, const {{class_name}} & lm, const FieldPtrOrIterator & it )
+ static void get( Vector3< {{dtype}} > & momentumDensity, const {{class_name}} & lm, const FieldPtrOrIterator & it )
{
const auto x = it.x();
const auto y = it.y();
const auto z = it.z();
{% for i in range(Q) -%}
- const real_t f_{{i}} = it[{{i}}];
+ const {{dtype}} f_{{i}} = it[{{i}}];
{% endfor -%}
{{momentum_density_getter | indent(8) }}
@@ -484,12 +485,12 @@ struct MomentumDensity< {{class_name}}>
}
template< typename PdfField_T >
- static void get( Vector3< real_t > & momentumDensity, const {{class_name}} & lm, const PdfField_T & pdf,
+ static void get( Vector3< {{dtype}} > & momentumDensity, const {{class_name}} & lm, const PdfField_T & pdf,
const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
{
- const real_t & xyz0 = pdf(x,y,z,0);
+ const {{dtype}} & xyz0 = pdf(x,y,z,0);
{% for i in range(Q) -%}
- const real_t f_{{i}} = pdf.getF( &xyz0, {{i}});
+ const {{dtype}} f_{{i}} = pdf.getF( &xyz0, {{i}});
{% endfor -%}
{{momentum_density_getter | indent(8) }}
@@ -504,13 +505,13 @@ template<>
struct PressureTensor<{{class_name}}>
{
template< typename FieldPtrOrIterator >
- static void get( Matrix3< real_t > & /* pressureTensor */, const {{class_name}} & /* latticeModel */, const FieldPtrOrIterator & /* it */ )
+ static void get( Matrix3< {{dtype}} > & /* pressureTensor */, const {{class_name}} & /* latticeModel */, const FieldPtrOrIterator & /* it */ )
{
WALBERLA_ABORT("Not implemented");
}
template< typename PdfField_T >
- static void get( Matrix3< real_t > & /* pressureTensor */, const {{class_name}} & /* latticeModel */, const PdfField_T & /* pdf */,
+ static void get( Matrix3< {{dtype}} > & /* pressureTensor */, const {{class_name}} & /* latticeModel */, const PdfField_T & /* pdf */,
const cell_idx_t /* x */, const cell_idx_t /* y */, const cell_idx_t /* z */ )
{
WALBERLA_ABORT("Not implemented");
@@ -522,27 +523,27 @@ template<>
struct ShearRate<{{class_name}}>
{
template< typename FieldPtrOrIterator >
- static inline real_t get( const {{class_name}} & /* latticeModel */, const FieldPtrOrIterator & /* it */,
- const Vector3< real_t > & /* velocity */, const real_t /* rho */)
+ static inline {{dtype}} get( const {{class_name}} & /* latticeModel */, const FieldPtrOrIterator & /* it */,
+ const Vector3< {{dtype}} > & /* velocity */, const {{dtype}} /* rho */)
{
WALBERLA_ABORT("Not implemented");
- return real_t(0.0);
+ return {{dtype}}(0.0);
}
template< typename PdfField_T >
- static inline real_t get( const {{class_name}} & latticeModel,
+ static inline {{dtype}} get( const {{class_name}} & latticeModel,
const PdfField_T & /* pdf */, const cell_idx_t /* x */, const cell_idx_t /* y */, const cell_idx_t /* z */,
- const Vector3< real_t > & /* velocity */, const real_t /* rho */ )
+ const Vector3< {{dtype}} > & /* velocity */, const {{dtype}} /* rho */ )
{
WALBERLA_ABORT("Not implemented");
- return real_t(0.0);
+ return {{dtype}}(0.0);
}
- static inline real_t get( const std::vector< real_t > & /* nonEquilibrium */, const real_t /* relaxationParam */,
- const real_t /* rho */ = real_t(1) )
+ static inline {{dtype}} get( const std::vector< {{dtype}} > & /* nonEquilibrium */, const {{dtype}} /* relaxationParam */,
+ const {{dtype}} /* rho */ = {{dtype}}(1) )
{
WALBERLA_ABORT("Not implemented");
- return real_t(0.0);
+ return {{dtype}}(0.0);
}
};
diff --git a/python/lbmpy_walberla/tests/test_walberla_codegen.py b/python/lbmpy_walberla/tests/test_walberla_codegen.py
index c36bf659f..4aae0a04e 100644
--- a/python/lbmpy_walberla/tests/test_walberla_codegen.py
+++ b/python/lbmpy_walberla/tests/test_walberla_codegen.py
@@ -78,6 +78,7 @@ class WalberlaLbmpyCodegenTest(unittest.TestCase):
stencil = LBStencil(stencil=Stencil.D3Q27)
lbm_config = LBMConfig(stencil=stencil, method=Method.TRT_KBC_N1,
+ zero_centered=False,
relaxation_rates=[1.5, sp.Symbol('omega_free')],
compressible=True, entropic=True,
omega_output_field=omega_field)
@@ -92,7 +93,7 @@ class WalberlaLbmpyCodegenTest(unittest.TestCase):
force_field, vel_field = ps.fields("force(3), velocity(3): [3D]", layout='fzyx')
stencil = LBStencil(stencil=Stencil.D3Q19)
- lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, compressible=True,
+ lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, compressible=True, zero_centered=False,
fluctuating={'seed': 0, 'temperature': 1e-6}, force_model=ForceModel.GUO,
relaxation_rates=[omega_shear] * 19, force=force_field.center_vector)
diff --git a/python/lbmpy_walberla/walberla_lbm_generation.py b/python/lbmpy_walberla/walberla_lbm_generation.py
index cb3d56f41..8eda3e76f 100644
--- a/python/lbmpy_walberla/walberla_lbm_generation.py
+++ b/python/lbmpy_walberla/walberla_lbm_generation.py
@@ -12,11 +12,12 @@ from lbmpy.updatekernels import create_lbm_kernel, create_stream_only_kernel
from pystencils import AssignmentCollection, create_kernel, Target
from pystencils.astnodes import SympyAssignment
from pystencils.backends.cbackend import CBackend, CustomSympyPrinter, get_headers
-from pystencils.data_types import TypedSymbol, type_all_numbers, cast_func
+from pystencils.typing import BasicType, CastFunc, TypedSymbol
from pystencils.field import Field
+from pystencils.node_collection import NodeCollection
from pystencils.stencil import offset_to_direction_string
from pystencils.sympyextensions import get_symmetric_part
-from pystencils.transformations import add_types
+from pystencils.typing.transformations import add_types
from pystencils_walberla.codegen import KernelInfo, config_from_context
from pystencils_walberla.jinja_filters import add_pystencils_filters_to_jinja_env
@@ -24,55 +25,75 @@ cpp_printer = CustomSympyPrinter()
REFINEMENT_SCALE_FACTOR = sp.Symbol("level_scale_factor")
-def __lattice_model(generation_context, class_name, lb_method, stream_collide_ast, collide_ast, stream_ast,
+def __type_equilibrium_assignments(assignments, config, subs_dict):
+ result = assignments.new_with_substitutions(subs_dict)
+ result = NodeCollection(result.main_assignments)
+ result.evaluate_terms()
+ result = add_types(result.all_assignments, config)
+ return result
+
+
+def __lattice_model(generation_context, class_name, config, lb_method, stream_collide_ast, collide_ast, stream_ast,
refinement_scaling):
stencil_name = lb_method.stencil.name
if not stencil_name:
raise ValueError("lb_method uses a stencil that is not supported in waLBerla")
communication_stencil_name = stencil_name if stencil_name != "D3Q15" else "D3Q27"
- is_float = not generation_context.double_accuracy
- dtype_string = "float32" if is_float else "float64"
+ default_dtype = config.data_type.default_factory()
+
+ cqc = lb_method.conserved_quantity_computation
- vel_symbols = lb_method.conserved_quantity_computation.first_order_moment_symbols
+ vel_symbols = cqc.velocity_symbols
rho_sym = sp.Symbol('rho')
+ reference_density = rho_sym if cqc.compressible else cqc.background_density
pdfs_sym = sp.symbols(f'f_:{lb_method.stencil.Q}')
- vel_arr_symbols = [IndexedBase(sp.Symbol('u'), shape=(1,))[i] for i in range(len(vel_symbols))]
+
+ vel_arr_symbols = [IndexedBase(TypedSymbol('u', default_dtype), shape=(1,))[i] for i in range(len(vel_symbols))]
+ subs_dict = {a: b for a, b in zip(vel_symbols, vel_arr_symbols)}
+
momentum_density_symbols = sp.symbols(f'md_:{len(vel_symbols)}')
equilibrium = lb_method.get_equilibrium()
- equilibrium = equilibrium.new_with_substitutions({a: b for a, b in zip(vel_symbols, vel_arr_symbols)})
- _, _, equilibrium = add_types(equilibrium.main_assignments, dtype_string, False)
- equilibrium = sp.Matrix([e.rhs for e in equilibrium])
+ lhs_list = [a.lhs for a in equilibrium.main_assignments]
+
+ equilibrium_matrix = sp.Matrix([e.rhs for e in equilibrium.main_assignments])
+ symmetric_equilibrium_matrix = get_symmetric_part(equilibrium_matrix, vel_symbols)
+ asymmetric_equilibrium_matrix = sp.expand(equilibrium_matrix - symmetric_equilibrium_matrix)
+
+ equilibrium = AssignmentCollection([ps.Assignment(lhs, rhs)
+ for lhs, rhs in zip(lhs_list, equilibrium_matrix)])
+ symmetric_equilibrium = AssignmentCollection([ps.Assignment(lhs, rhs)
+ for lhs, rhs in zip(lhs_list, symmetric_equilibrium_matrix)])
+ asymmetric_equilibrium = AssignmentCollection([ps.Assignment(lhs, rhs)
+ for lhs, rhs in zip(lhs_list, asymmetric_equilibrium_matrix)])
- symmetric_equilibrium = get_symmetric_part(equilibrium, vel_arr_symbols)
- symmetric_equilibrium = symmetric_equilibrium.subs(sp.Rational(1, 2), cast_func(sp.Rational(1, 2), dtype_string))
- asymmetric_equilibrium = sp.expand(equilibrium - symmetric_equilibrium)
+ equilibrium = __type_equilibrium_assignments(equilibrium, config, subs_dict)
+ symmetric_equilibrium = __type_equilibrium_assignments(symmetric_equilibrium, config, subs_dict)
+ asymmetric_equilibrium = __type_equilibrium_assignments(asymmetric_equilibrium, config, subs_dict)
force_model = lb_method.force_model
macroscopic_velocity_shift = None
if force_model:
if hasattr(force_model, 'macroscopic_velocity_shift'):
macroscopic_velocity_shift = [e.subs(force_model.subs_dict_force)
- for e in force_model.macroscopic_velocity_shift(rho_sym)]
- macroscopic_velocity_shift = [expression_to_code(e.subs(sp.Rational(1, 2), cast_func(sp.Rational(1, 2),
- dtype_string)),
- "lm.", ['rho'], dtype=dtype_string)
+ for e in force_model.macroscopic_velocity_shift(reference_density)]
+ macroscopic_velocity_shift = [expression_to_code(e, "lm.", ['rho'], dtype=default_dtype)
for e in macroscopic_velocity_shift]
- cqc = lb_method.conserved_quantity_computation
-
eq_input_from_input_eqs = cqc.equilibrium_input_equations_from_init_values(sp.Symbol('rho_in'), vel_arr_symbols)
- density_velocity_setter_macroscopic_values = equations_to_code(eq_input_from_input_eqs, dtype=dtype_string,
+ density_velocity_setter_macroscopic_values = equations_to_code(eq_input_from_input_eqs, dtype=default_dtype,
variables_without_prefix=['rho_in', 'u'])
momentum_density_getter = cqc.output_equations_from_pdfs(pdfs_sym, {'density': rho_sym,
'momentum_density': momentum_density_symbols})
+
+ is_float = True if issubclass(default_dtype.numpy_dtype.type, np.float32) else False
constant_suffix = "f" if is_float else ""
required_headers = get_headers(stream_collide_ast)
if refinement_scaling:
- refinement_scaling_info = [(e0, e1, expression_to_code(e2, '', dtype=dtype_string)) for e0, e1, e2 in
+ refinement_scaling_info = [(e0, e1, expression_to_code(e2, '', dtype=default_dtype)) for e0, e1, e2 in
refinement_scaling.scaling_info]
# append '_' to entries since they are used as members
for i in range(len(refinement_scaling_info)):
@@ -95,17 +116,18 @@ def __lattice_model(generation_context, class_name, lb_method, stream_collide_as
'compressible': lb_method.conserved_quantity_computation.compressible,
'weights': ",".join(str(w.evalf()) + constant_suffix for w in lb_method.weights),
'inverse_weights': ",".join(str((1 / w).evalf()) + constant_suffix for w in lb_method.weights),
+ 'dtype': "float" if is_float else "double",
'equilibrium_from_direction': stencil_switch_statement(lb_method.stencil, equilibrium),
'symmetric_equilibrium_from_direction': stencil_switch_statement(lb_method.stencil, symmetric_equilibrium),
'asymmetric_equilibrium_from_direction': stencil_switch_statement(lb_method.stencil, asymmetric_equilibrium),
- 'equilibrium': [cpp_printer.doprint(e) for e in equilibrium],
+ 'equilibrium': [cpp_printer.doprint(e.rhs) for e in equilibrium],
'macroscopic_velocity_shift': macroscopic_velocity_shift,
'density_getters': equations_to_code(cqc.output_equations_from_pdfs(pdfs_sym, {"density": rho_sym}),
- variables_without_prefix=[e.name for e in pdfs_sym], dtype=dtype_string),
+ variables_without_prefix=[e.name for e in pdfs_sym], dtype=default_dtype),
'momentum_density_getter': equations_to_code(momentum_density_getter, variables_without_prefix=pdfs_sym,
- dtype=dtype_string),
+ dtype=default_dtype),
'density_velocity_setter_macroscopic_values': density_velocity_setter_macroscopic_values,
'refinement_scaling_info': refinement_scaling_info,
@@ -140,7 +162,6 @@ def generate_lattice_model(generation_context, class_name, collision_rule, field
# usually a numpy layout is chosen by default i.e. xyzf - which is bad for waLBerla where at least the spatial
# coordinates should be ordered in reverse direction i.e. zyx
- dtype = np.float64 if config.data_type == "float64" else np.float32
lb_method = collision_rule.method
q = lb_method.stencil.Q
@@ -154,8 +175,10 @@ def generate_lattice_model(generation_context, class_name, collision_rule, field
elif field_layout == 'zyxf':
config.cpu_vectorize_info['assume_inner_stride_one'] = False
- src_field = ps.Field.create_generic('pdfs', dim, dtype, index_dimensions=1, layout=field_layout, index_shape=(q,))
- dst_field = ps.Field.create_generic('pdfs_tmp', dim, dtype, index_dimensions=1, layout=field_layout,
+ src_field = ps.Field.create_generic('pdfs', dim, config.data_type['pdfs'].numpy_dtype,
+ index_dimensions=1, layout=field_layout, index_shape=(q,))
+ dst_field = ps.Field.create_generic('pdfs_tmp', dim, config.data_type['pdfs_tmp'].numpy_dtype,
+ index_dimensions=1, layout=field_layout,
index_shape=(q,))
stream_collide_update_rule = create_lbm_kernel(collision_rule, src_field, dst_field, StreamPullTwoFieldsAccessor())
@@ -173,7 +196,7 @@ def generate_lattice_model(generation_context, class_name, collision_rule, field
stream_ast = create_kernel(stream_update_rule, config=config)
stream_ast.function_name = 'kernel_stream'
stream_ast.assumed_inner_stride_one = config.cpu_vectorize_info['assume_inner_stride_one']
- __lattice_model(generation_context, class_name, lb_method, stream_collide_ast, collide_ast, stream_ast,
+ __lattice_model(generation_context, class_name, config, lb_method, stream_collide_ast, collide_ast, stream_ast,
refinement_scaling)
@@ -234,11 +257,11 @@ def stencil_switch_statement(stencil, values):
case {{direction_name}}: return {{value}};
{% endfor -%}
default:
- WALBERLA_ABORT("Invalid Direction");
+ WALBERLA_ABORT("Invalid Direction")
}
""")
- dir_to_value_dict = {offset_to_direction_string(d): cpp_printer.doprint(v) for d, v in zip(stencil, values)}
+ dir_to_value_dict = {offset_to_direction_string(d): cpp_printer.doprint(v.rhs) for d, v in zip(stencil, values)}
return template.render(dir_to_value_dict=dir_to_value_dict, undefined=StrictUndefined)
@@ -271,7 +294,7 @@ def field_and_symbol_substitute(expr, variable_prefix="lm.", variables_without_p
return expr.subs(substitutions)
-def expression_to_code(expr, variable_prefix="lm.", variables_without_prefix=None, dtype="double"):
+def expression_to_code(expr, variable_prefix="lm.", variables_without_prefix=None, dtype=None):
"""
Takes a sympy expression and creates a C code string from it. Replaces field accesses by
walberla field accesses i.e. field_W^1 -> field->get(-1, 0, 0, 1)
@@ -282,6 +305,9 @@ def expression_to_code(expr, variable_prefix="lm.", variables_without_prefix=Non
:param variables_without_prefix: this variables are not prefixed
:return: code string
"""
+ if dtype is None:
+ dtype = BasicType("float64")
+
if variables_without_prefix is None:
variables_without_prefix = []
return cpp_printer.doprint(
@@ -289,18 +315,15 @@ def expression_to_code(expr, variable_prefix="lm.", variables_without_prefix=Non
def type_expr(eq, dtype):
- def recurse(expr):
- for i in range(len(expr.args)):
- if expr.args[i] == sp.Rational or expr.args[i] == sp.Float:
- expr.args[i] = type_all_numbers(expr.args[i], dtype=dtype)
- else:
- recurse(expr.args[i])
-
- recurse(eq)
+ # manually cast 0.5 to dtype since this is somehow not done automatically
+ eq = eq.subs(sp.Rational(1, 2), CastFunc(sp.Rational(1, 2), dtype))
return eq.subs({s: TypedSymbol(s.name, dtype) for s in eq.atoms(sp.Symbol)})
-def equations_to_code(equations, variable_prefix="lm.", variables_without_prefix=None, dtype="double"):
+def equations_to_code(equations, variable_prefix="lm.", variables_without_prefix=None, dtype=None):
+ if dtype is None:
+ dtype = BasicType("float64")
+
if variables_without_prefix is None:
variables_without_prefix = []
if isinstance(equations, AssignmentCollection):
diff --git a/python/pystencils_walberla/boundary.py b/python/pystencils_walberla/boundary.py
index 5cf8f1fe5..329b68059 100644
--- a/python/pystencils_walberla/boundary.py
+++ b/python/pystencils_walberla/boundary.py
@@ -5,7 +5,7 @@ from pystencils.boundaries.boundaryhandling import create_boundary_kernel
from pystencils.boundaries.createindexlist import (
boundary_index_array_coordinate_names, direction_member_name,
numpy_data_type_for_boundary_object)
-from pystencils.data_types import TypedSymbol, create_type
+from pystencils.typing import TypedSymbol, create_type
from pystencils.stencil import inverse_direction
from pystencils_walberla.codegen import config_from_context
@@ -46,8 +46,10 @@ def generate_boundary(generation_context,
create_kernel_params = config.__dict__
del create_kernel_params['target']
del create_kernel_params['index_fields']
+ del create_kernel_params['default_number_int']
+ del create_kernel_params['skip_independence_check']
- field_data_type = np.float64 if config.data_type == "float64" else np.float32
+ field_data_type = config.data_type[field_name].numpy_dtype
index_struct_dtype = numpy_data_type_for_boundary_object(boundary_object, dim)
@@ -57,7 +59,7 @@ def generate_boundary(generation_context,
field_type=field_type)
index_field = Field('indexVector', FieldType.INDEXED, index_struct_dtype, layout=[0],
- shape=(TypedSymbol("indexVectorSize", create_type(np.int64)), 1), strides=(1, 1))
+ shape=(TypedSymbol("indexVectorSize", create_type("int32")), 1), strides=(1, 1))
if not kernel_creation_function:
kernel_creation_function = create_boundary_kernel
@@ -79,6 +81,9 @@ def generate_boundary(generation_context,
if additional_data_handler is None:
additional_data_handler = AdditionalDataHandler(stencil=neighbor_stencil)
+ default_dtype = config.data_type.default_factory()
+ is_float = True if issubclass(default_dtype.numpy_dtype.type, np.float32) else False
+
context = {
'kernel': kernel_family,
'class_name': boundary_object.name,
@@ -92,6 +97,7 @@ def generate_boundary(generation_context,
'inner_or_boundary': boundary_object.inner_or_boundary,
'single_link': boundary_object.single_link,
'additional_data_handler': additional_data_handler,
+ 'dtype': "double" if is_float else "float",
'layout': layout
}
diff --git a/python/pystencils_walberla/codegen.py b/python/pystencils_walberla/codegen.py
index de792b8fc..c5cefc06e 100644
--- a/python/pystencils_walberla/codegen.py
+++ b/python/pystencils_walberla/codegen.py
@@ -14,8 +14,8 @@ from pystencils.backends.simd_instruction_sets import get_supported_instruction_
from pystencils.stencil import inverse_direction, offset_to_direction_string
from pystencils.backends.cuda_backend import CudaSympyPrinter
-from pystencils.kernelparameters import SHAPE_DTYPE
-from pystencils.data_types import TypedSymbol
+from pystencils.typing.typed_sympy import SHAPE_DTYPE
+from pystencils.typing import TypedSymbol
from pystencils_walberla.jinja_filters import add_pystencils_filters_to_jinja_env
from pystencils_walberla.kernel_selection import KernelCallNode, KernelFamily, HighLevelInterfaceSpec
@@ -275,6 +275,9 @@ def generate_pack_info(generation_context, class_name: str,
config = replace(config, cpu_vectorize_info=None)
config_zero_gl = replace(config_zero_gl, cpu_vectorize_info=None)
+ config = replace(config, allow_double_writes=True)
+ config_zero_gl = replace(config_zero_gl, allow_double_writes=True)
+
template_name = "CpuPackInfo.tmpl" if config.target == Target.CPU else 'GpuPackInfo.tmpl'
fields_accessed = set()
@@ -500,7 +503,7 @@ def config_from_context(generation_context, target=Target.CPU, data_type=None,
cpu_vectorize_info['assume_aligned'] = cpu_vectorize_info.get('assume_aligned', False)
cpu_vectorize_info['nontemporal'] = cpu_vectorize_info.get('nontemporal', False)
- config = CreateKernelConfig(target=target, data_type=data_type,
+ config = CreateKernelConfig(target=target, data_type=data_type, default_number_float=data_type,
cpu_openmp=cpu_openmp, cpu_vectorize_info=cpu_vectorize_info,
**kwargs)
diff --git a/python/pystencils_walberla/jinja_filters.py b/python/pystencils_walberla/jinja_filters.py
index b9cfd32d0..0ee3f3cd1 100644
--- a/python/pystencils_walberla/jinja_filters.py
+++ b/python/pystencils_walberla/jinja_filters.py
@@ -10,7 +10,7 @@ import sympy as sp
from pystencils import Target, Backend
from pystencils.backends.cbackend import generate_c
-from pystencils.data_types import TypedSymbol, get_base_type
+from pystencils.typing import TypedSymbol, get_base_type
from pystencils.field import FieldType
from pystencils.sympyextensions import prod
@@ -78,7 +78,7 @@ def get_field_fsize(field):
def get_field_stride(param):
field = param.fields[0]
- type_str = get_base_type(param.symbol.dtype).base_name
+ type_str = get_base_type(param.symbol.dtype).c_name
stride_names = ['xStride()', 'yStride()', 'zStride()', 'fStride()']
stride_names = [f"{type_str}({param.field_name}->{e})" for e in stride_names]
strides = stride_names[:field.spatial_dimensions]
@@ -322,12 +322,12 @@ def generate_call(ctx, kernel, ghost_layers_to_include=0, cell_interval=None, st
f"{instruction_set['bytes']}, 0);")
elif param.is_field_stride:
casted_stride = get_field_stride(param)
- type_str = param.symbol.dtype.base_name
+ type_str = param.symbol.dtype.c_name
kernel_call_lines.append(f"const {type_str} {param.symbol.name} = {casted_stride};")
elif param.is_field_shape:
coord = param.symbol.coordinate
field = param.fields[0]
- type_str = param.symbol.dtype.base_name
+ type_str = param.symbol.dtype.c_name
shape = f"{type_str}({get_end_coordinates(field)[coord]})"
assert coord < 3
max_value = f"{field.name}->{('x', 'y', 'z')[coord]}SizeWithGhostLayer()"
@@ -364,11 +364,16 @@ def generate_constructor_initializer_list(kernel_info, parameters_to_ignore=None
parameters_to_ignore += kernel_info.temporary_fields
parameter_initializer_list = []
+ # First field pointer
for param in kernel_info.parameters:
if param.is_field_pointer and param.field_name not in parameters_to_ignore:
parameter_initializer_list.append(f"{param.field_name}ID({param.field_name}ID_)")
- elif not param.is_field_parameter and param.symbol.name not in parameters_to_ignore:
+
+ # Then free parameters
+ for param in kernel_info.parameters:
+ if not param.is_field_parameter and param.symbol.name not in parameters_to_ignore:
parameter_initializer_list.append(f"{param.symbol.name}_({param.symbol.name})")
+
return ", ".join(parameter_initializer_list)
@@ -383,11 +388,16 @@ def generate_constructor_parameters(kernel_info, parameters_to_ignore=None):
parameters_to_ignore += kernel_info.temporary_fields + varying_parameter_names
parameter_list = []
+ # First field pointer
for param in kernel_info.parameters:
if param.is_field_pointer and param.field_name not in parameters_to_ignore:
parameter_list.append(f"BlockDataID {param.field_name}ID_")
- elif not param.is_field_parameter and param.symbol.name not in parameters_to_ignore:
+
+ # Then free parameters
+ for param in kernel_info.parameters:
+ if not param.is_field_parameter and param.symbol.name not in parameters_to_ignore:
parameter_list.append(f"{param.symbol.dtype} {param.symbol.name}")
+
varying_parameters = ["%s %s" % e for e in varying_parameters]
return ", ".join(parameter_list + varying_parameters)
@@ -427,7 +437,11 @@ def generate_members(ctx, kernel_info, parameters_to_ignore=(), only_fields=Fals
continue
if param.is_field_pointer and param.field_name not in params_to_skip:
result.append(f"BlockDataID {param.field_name}ID;")
- elif not param.is_field_parameter and param.symbol.name not in params_to_skip:
+
+ for param in kernel_info.parameters:
+ if only_fields and not param.is_field_parameter:
+ continue
+ if not param.is_field_parameter and param.symbol.name not in params_to_skip:
result.append(f"{param.symbol.dtype} {param.symbol.name}_;")
for field_name in kernel_info.temporary_fields:
diff --git a/python/pystencils_walberla/kernel_selection.py b/python/pystencils_walberla/kernel_selection.py
index 5215e599e..b2e831cc8 100644
--- a/python/pystencils_walberla/kernel_selection.py
+++ b/python/pystencils_walberla/kernel_selection.py
@@ -6,7 +6,7 @@ from jinja2.filters import do_indent
from pystencils import Target, TypedSymbol
from pystencils.backends.cbackend import get_headers
from pystencils.backends.cuda_backend import CudaSympyPrinter
-from pystencils.kernelparameters import SHAPE_DTYPE
+from pystencils.typing.typed_sympy import SHAPE_DTYPE
"""
diff --git a/python/pystencils_walberla/templates/Boundary.tmpl.cpp b/python/pystencils_walberla/templates/Boundary.tmpl.cpp
index 9c41e7965..a7b1c064e 100644
--- a/python/pystencils_walberla/templates/Boundary.tmpl.cpp
+++ b/python/pystencils_walberla/templates/Boundary.tmpl.cpp
@@ -72,7 +72,7 @@ void {{class_name}}::run_impl(
)
{
auto * indexVectors = block->getData<IndexVectors>(indexVectorID);
- int64_t indexVectorSize = int64_c( indexVectors->indexVector(type).size() );
+ int32_t indexVectorSize = int32_c( indexVectors->indexVector(type).size() );
if( indexVectorSize == 0)
return;
diff --git a/python/pystencils_walberla/templates/Boundary.tmpl.h b/python/pystencils_walberla/templates/Boundary.tmpl.h
index 8a673d32e..03b403aea 100644
--- a/python/pystencils_walberla/templates/Boundary.tmpl.h
+++ b/python/pystencils_walberla/templates/Boundary.tmpl.h
@@ -80,7 +80,7 @@ public:
{% endif -%}
CpuIndexVector & indexVector(Type t) { return cpuVectors_[t]; }
- {{StructName}} * pointerCpu(Type t) { return &(cpuVectors_[t][0]); }
+ {{StructName}} * pointerCpu(Type t) { return cpuVectors_[t].data(); }
{% if target == 'gpu' -%}
{{StructName}} * pointerGpu(Type t) { return gpuVectors_[t]; }
@@ -238,7 +238,7 @@ public:
{%endif%}
{%else%}
auto flagWithGLayers = flagField->xyzSizeWithGhostLayer();
- real_t dot = 0.0; real_t maxn = 0.0;
+ {{dtype}} dot = 0.0; {{dtype}} maxn = 0.0;
cell_idx_t calculated_idx = 0;
cell_idx_t dx = 0; cell_idx_t dy = 0; {%if dim == 3%} cell_idx_t dz = 0; {% endif %}
cell_idx_t sum_x = 0; cell_idx_t sum_y = 0; {%if dim == 3%} cell_idx_t sum_z = 0; {%endif %}
@@ -272,7 +272,7 @@ public:
{
{%- for dirIdx, dirVec, offset in additional_data_handler.stencil_info %}
dx = {{dirVec[0]}}; dy = {{dirVec[1]}}; {%if dim == 3%} dz = {{dirVec[2]}}; {% endif %}
- dot = real_c( dx*sum_x + dy*sum_y {%if dim == 3%} + dz*sum_z {% endif %});
+ dot = {{dtype}}( dx*sum_x + dy*sum_y {%if dim == 3%} + dz*sum_z {% endif %});
if (dot > maxn)
{
maxn = dot;
diff --git a/python/pystencils_walberla/utility.py b/python/pystencils_walberla/utility.py
index 8800130bc..c109265ef 100644
--- a/python/pystencils_walberla/utility.py
+++ b/python/pystencils_walberla/utility.py
@@ -1,5 +1,5 @@
from os import path
-from pystencils.data_types import get_base_type
+from pystencils.typing import get_base_type
from pystencils_walberla.cmake_integration import CodeGenerationContext
from lbmpy import LBStencil
diff --git a/python/waLBerla_tests/test_blockforest.py b/python/waLBerla_tests/test_blockforest.py
index 09a158151..296cb570b 100644
--- a/python/waLBerla_tests/test_blockforest.py
+++ b/python/waLBerla_tests/test_blockforest.py
@@ -1,14 +1,17 @@
import unittest
import numpy as np
+import waLBerla as wlb
from waLBerla import field, createUniformBlockGrid, AABB
+field_data_type = np.float64 if wlb.build_info.build_with_double_accuracy else np.float32
+
class BlockforestModuleTest(unittest.TestCase):
def testMemoryManagement1(self):
"""Testing correct reference counting of block data"""
blocks = createUniformBlockGrid(blocks=(1, 1, 1), cellsPerBlock=(2, 2, 2))
- field.addToStorage(blocks, "TestField", np.float64)
+ field.addToStorage(blocks, "TestField", field_data_type)
f = blocks[0]["TestField"]
strides_before = f.strides
del blocks
@@ -24,7 +27,7 @@ class BlockforestModuleTest(unittest.TestCase):
"""Testing correct reference counting of block data
Holding only a numpy array pointing to a waLBerla field should still hold the blockstructure alive"""
blocks = createUniformBlockGrid(blocks=(1, 1, 1), cellsPerBlock=(2, 2, 2))
- field.addToStorage(blocks, "TestField", np.float64)
+ field.addToStorage(blocks, "TestField", field_data_type)
npf = field.toArray(blocks[0]["TestField"])
npf[:, :, :] = 42.0
del blocks
@@ -36,7 +39,7 @@ class BlockforestModuleTest(unittest.TestCase):
def testMemoryManagement3(self):
"""Same as testMemoryManagement2, but with iterators"""
blocks = createUniformBlockGrid(blocks=(1, 1, 1), cellsPerBlock=(2, 2, 2))
- field.addToStorage(blocks, "TestField", np.float64)
+ field.addToStorage(blocks, "TestField", field_data_type)
for block in blocks:
for name in block.fieldNames:
if name == "TestField":
diff --git a/python/waLBerla_tests/test_field.py b/python/waLBerla_tests/test_field.py
index 04ab5afb4..e97fdd26d 100644
--- a/python/waLBerla_tests/test_field.py
+++ b/python/waLBerla_tests/test_field.py
@@ -3,12 +3,14 @@ import numpy as np
import waLBerla as wlb
from waLBerla import field, createUniformBlockGrid
+field_data_type = np.float64 if wlb.build_info.build_with_double_accuracy else np.float32
+
class FieldModuleTest(unittest.TestCase):
def testFieldAsBlockData(self):
blocks = createUniformBlockGrid(blocks=(1, 1, 1), cellsPerBlock=(3, 2, 2), periodic=(True, False, False))
- field.addToStorage(blocks, 'myField', np.float64, fSize=3, ghostLayers=0, initValue=0.0)
+ field.addToStorage(blocks, 'myField', field_data_type, fSize=3, ghostLayers=0, initValue=0.0)
my_field = wlb.field.toArray(blocks[0]['myField'])
self.assertEqual(my_field[0, 0, 0, 0], 0)
my_field[0, 0, 0, 0] = 42.0
@@ -18,8 +20,8 @@ class FieldModuleTest(unittest.TestCase):
def testNumpyConversionWithoutGhostLayers(self):
blocks = createUniformBlockGrid(blocks=(1, 1, 1), cellsPerBlock=(1, 2, 3), periodic=(True, False, False))
- field.addToStorage(blocks, 'f1', np.float64, fSize=2, ghostLayers=0, initValue=2.0)
- field.addToStorage(blocks, 'f2', np.float64, fSize=3, ghostLayers=0, initValue=2.0)
+ field.addToStorage(blocks, 'f1', field_data_type, fSize=2, ghostLayers=0, initValue=2.0)
+ field.addToStorage(blocks, 'f2', field_data_type, fSize=3, ghostLayers=0, initValue=2.0)
f1np = field.toArray(blocks[0]['f1'])
f2np = field.toArray(blocks[0]['f2'])
@@ -31,7 +33,7 @@ class FieldModuleTest(unittest.TestCase):
size = (10, 5, 4)
gl = 3
blocks = createUniformBlockGrid(blocks=(1, 1, 1), cellsPerBlock=size, periodic=(True, False, False))
- field.addToStorage(blocks, 'f', np.float64, fSize=3, ghostLayers=gl, initValue=0.0)
+ field.addToStorage(blocks, 'f', field_data_type, fSize=3, ghostLayers=gl, initValue=0.0)
f = blocks[0]['f']
diff --git a/src/core/DataTypes.h b/src/core/DataTypes.h
index 397f2c8d6..2f8687197 100644
--- a/src/core/DataTypes.h
+++ b/src/core/DataTypes.h
@@ -164,7 +164,7 @@ template< typename T > inline cell_idx_t cell_idx_c( T t ) { return numeric_cast
#ifdef WALBERLA_DOUBLE_ACCURACY
using real_t = double;
#else
-typedef float real_t;
+using real_t = float;
#endif
inline constexpr real_t operator"" _r( long double t ) { return static_cast< real_t >(t); }
diff --git a/src/python_coupling/export/BasicExport.cpp b/src/python_coupling/export/BasicExport.cpp
index 22b444d1b..efbb0226c 100644
--- a/src/python_coupling/export/BasicExport.cpp
+++ b/src/python_coupling/export/BasicExport.cpp
@@ -629,6 +629,11 @@ void exportBuildInfo(py::module_ &m)
m2.attr("build_machine" ) = WALBERLA_BUILD_MACHINE;
m2.attr("source_dir") = WALBERLA_SOURCE_DIR;
m2.attr("build_dir") = WALBERLA_BUILD_DIR;
+#ifdef WALBERLA_DOUBLE_ACCURACY
+ m2.attr("build_with_double_accuracy") = true;
+#else
+ m2.attr("build_with_double_accuracy") = false;
+#endif
}
diff --git a/tests/cuda/SimpleKernelTest.cpp b/tests/cuda/SimpleKernelTest.cpp
index 6b79fa8db..f2f9a2a8b 100644
--- a/tests/cuda/SimpleKernelTest.cpp
+++ b/tests/cuda/SimpleKernelTest.cpp
@@ -42,23 +42,23 @@ void kernel_double( cuda::FieldAccessor<double> f );
GhostLayerField<double,1> * createCPUField( IBlock* const block, StructuredBlockStorage* const storage )
{
return new GhostLayerField<double,1> (
- storage->getNumberOfXCells( *block ), // number of cells in x direction
- storage->getNumberOfYCells( *block ), // number of cells in y direction
- storage->getNumberOfZCells( *block ), // number of cells in z direction
- 1, // number of ghost layers
- double(1), // initial value
- field::fzyx);
+ storage->getNumberOfXCells( *block ), // number of cells in x direction
+ storage->getNumberOfYCells( *block ), // number of cells in y direction
+ storage->getNumberOfZCells( *block ), // number of cells in z direction
+ 1, // number of ghost layers
+ double(1), // initial value
+ field::fzyx);
}
cuda::GPUField<double> * createGPUField( IBlock* const block, StructuredBlockStorage* const storage )
{
return new cuda::GPUField<double> (
- storage->getNumberOfXCells( *block ), // number of cells in x direction
- storage->getNumberOfYCells( *block ), // number of cells in y direction
- storage->getNumberOfZCells( *block ), // number of cells in z direction
- 1, // fSize
- 1, // number of ghost layers
- field::fzyx );
+ storage->getNumberOfXCells( *block ), // number of cells in x direction
+ storage->getNumberOfYCells( *block ), // number of cells in y direction
+ storage->getNumberOfZCells( *block ), // number of cells in z direction
+ 1, // fSize
+ 1, // number of ghost layers
+ field::fzyx );
}
@@ -68,11 +68,11 @@ int main( int argc, char ** argv )
debug::enterTestMode();
shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid (
- uint_t(1), uint_t(1), uint_t(1), // number of blocks in x,y,z direction
- uint_t(14), uint_t(14), uint_t(14), // how many cells per block (x,y,z)
- real_c(0.5), // dx: length of one cell in physical coordinates
- false, // one block per process - "false" means all blocks to one process
- false, false, false ); // no periodicity
+ uint_t(1), uint_t(1), uint_t(1), // number of blocks in x,y,z direction
+ uint_t(14), uint_t(14), uint_t(14), // how many cells per block (x,y,z)
+ real_c(0.5), // dx: length of one cell in physical coordinates
+ false, // one block per process - "false" means all blocks to one process
+ false, false, false ); // no periodicity
@@ -98,7 +98,7 @@ int main( int argc, char ** argv )
cuda::fieldCpy( *cpuField, *gpuField );
- WALBERLA_ASSERT_FLOAT_EQUAL( cpuField->get(0,0,0), real_t(2) );
+ WALBERLA_ASSERT_FLOAT_EQUAL( cpuField->get(0,0,0), real_t(2) )
}
@@ -108,5 +108,5 @@ int main( int argc, char ** argv )
//gui.run();
- return 0;
+ return EXIT_SUCCESS;
}
diff --git a/tests/cuda/codegen/CodegenJacobiGPU.cpp b/tests/cuda/codegen/CodegenJacobiGPU.cpp
index 7842220f6..93814e0a5 100644
--- a/tests/cuda/codegen/CodegenJacobiGPU.cpp
+++ b/tests/cuda/codegen/CodegenJacobiGPU.cpp
@@ -49,8 +49,8 @@
using namespace walberla;
-typedef GhostLayerField<double,1> ScalarField;
-typedef cuda::GPUField<double> GPUField;
+typedef GhostLayerField<real_t, 1> ScalarField;
+typedef cuda::GPUField<real_t> GPUField;
ScalarField * createField( IBlock* const block, StructuredBlockStorage* const storage )
@@ -62,7 +62,7 @@ ScalarField * createField( IBlock* const block, StructuredBlockStorage* const st
1, // one ghost layer
double(0), // initial value
field::fzyx, // layout
- make_shared<cuda::HostFieldAllocator<double> >() // allocator for host pinned memory
+ make_shared<cuda::HostFieldAllocator<real_t> >() // allocator for host pinned memory
);
}
@@ -75,7 +75,7 @@ void testJacobi2D()
shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid (
uint_t(1) , uint_t(1), uint_t(1), // number of blocks in x,y,z direction
xSize, ySize, uint_t(1), // how many cells per block (x,y,z)
- real_t(1), // dx: length of one cell in physical coordinates
+ real_c(1.0), // dx: length of one cell in physical coordinates
false, // one block per process - "false" means all blocks to one process
true, true, true ); // no periodicity
@@ -90,7 +90,7 @@ void testJacobi2D()
auto f = blockIt->getData<ScalarField>( cpuFieldID );
for( cell_idx_t y = 0; y < cell_idx_c( f->ySize() / 2 ); ++y )
for( cell_idx_t x = 0; x < cell_idx_c( f->xSize() / 2 ); ++x )
- f->get( x, y, 0 ) = 1.0;
+ f->get( x, y, 0 ) = real_c(1.0);
}
typedef blockforest::communication::UniformBufferedScheme<stencil::D2Q9> CommScheme;
@@ -114,7 +114,7 @@ void testJacobi2D()
auto firstBlock = blocks->begin();
auto f = firstBlock->getData<ScalarField>( cpuFieldID );
- WALBERLA_CHECK_FLOAT_EQUAL(f->get(0,0,0), real_t(1.0 / 4.0));
+ WALBERLA_CHECK_FLOAT_EQUAL(f->get(0,0,0), real_c(1.0 / 4.0))
}
@@ -128,7 +128,7 @@ void testJacobi3D()
shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid (
uint_t(1) , uint_t(1), uint_t(1), // number of blocks in x,y,z direction
xSize, ySize, zSize, // how many cells per block (x,y,z)
- real_t(1), // dx: length of one cell in physical coordinates
+ real_c(1.0), // dx: length of one cell in physical coordinates
false, // one block per process - "false" means all blocks to one process
true, true, true ); // no periodicity
@@ -144,7 +144,7 @@ void testJacobi3D()
for( cell_idx_t z = 0; z < cell_idx_c( f->zSize() / 2 ); ++z )
for( cell_idx_t y = 0; y < cell_idx_c( f->ySize() / 2 ); ++y )
for( cell_idx_t x = 0; x < cell_idx_c( f->xSize() / 2 ); ++x )
- f->get( x, y, z ) = 1.0;
+ f->get( x, y, z ) = real_c(1.0);
}
typedef blockforest::communication::UniformBufferedScheme<stencil::D3Q7> CommScheme;
@@ -168,7 +168,7 @@ void testJacobi3D()
auto firstBlock = blocks->begin();
auto f = firstBlock->getData<ScalarField>( cpuFieldID );
- WALBERLA_CHECK_FLOAT_EQUAL(f->get(0,0,0), real_t(1.0 / 8.0));
+ WALBERLA_CHECK_FLOAT_EQUAL(f->get(0,0,0), real_c(1.0 / 8.0))
}
int main( int argc, char ** argv )
@@ -179,5 +179,5 @@ int main( int argc, char ** argv )
testJacobi2D();
testJacobi3D();
- return 0;
+ return EXIT_SUCCESS;
}
diff --git a/tests/cuda/codegen/CodegenPoissonGPU.cpp b/tests/cuda/codegen/CodegenPoissonGPU.cpp
index d160f7999..ef5ae96c0 100644
--- a/tests/cuda/codegen/CodegenPoissonGPU.cpp
+++ b/tests/cuda/codegen/CodegenPoissonGPU.cpp
@@ -88,26 +88,26 @@ void testPoisson()
{
const uint_t xCells = uint_t(200);
const uint_t yCells = uint_t(100);
- const real_t xSize = real_t(2);
- const real_t ySize = real_t(1);
+ const real_t xSize = real_c(2.0);
+ const real_t ySize = real_c(1.0);
const real_t dx = xSize / real_c( xCells + uint_t(1) );
const real_t dy = ySize / real_c( yCells + uint_t(1) );
// Create blocks
shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid (
- math::AABB( real_t(0.5) * dx, real_t(0.5) * dy, real_t(0),
- xSize - real_t(0.5) * dx, ySize - real_t(0.5) * dy, dx ),
+ math::AABB( real_c(0.5) * dx, real_c(0.5) * dy, real_c(0.0),
+ xSize - real_c(0.5) * dx, ySize - real_c(0.5) * dy, dx ),
uint_t(1) , uint_t(1), uint_t(1), // number of blocks in x,y,z direction
xCells, yCells, uint_t(1), // how many cells per block (x,y,z)
false, // one block per process - "false" means all blocks to one process
false, false, false ); // no periodicity
- BlockDataID cpuFieldID = field::addToStorage< ScalarField_T >( blocks, "CPU Field src", real_t(0.0) );
+ BlockDataID cpuFieldID = field::addToStorage< ScalarField_T >( blocks, "CPU Field src", real_c(0.0) );
BlockDataID gpuField = cuda::addGPUFieldToStorage<ScalarField_T>( blocks, cpuFieldID, "GPU Field src" );
initU( blocks, cpuFieldID );
- BlockDataID cpufId = field::addToStorage< ScalarField_T >( blocks, "CPU Field f", real_t(0.0));
+ BlockDataID cpufId = field::addToStorage< ScalarField_T >( blocks, "CPU Field f", real_c(0.0));
BlockDataID gpufId = cuda::addGPUFieldToStorage<ScalarField_T>( blocks, cpufId, "GPU Field f" );
initF( blocks, cpufId );
@@ -145,5 +145,5 @@ int main( int argc, char ** argv )
testPoisson();
- return 0;
+ return EXIT_SUCCESS;
}
diff --git a/tests/cuda/codegen/CudaJacobiKernel.py b/tests/cuda/codegen/CudaJacobiKernel.py
index 79c2d767c..60fd64c1d 100644
--- a/tests/cuda/codegen/CudaJacobiKernel.py
+++ b/tests/cuda/codegen/CudaJacobiKernel.py
@@ -4,8 +4,9 @@ from pystencils_walberla import CodeGeneration, generate_sweep
with CodeGeneration() as ctx:
+ field_type = "float64" if ctx.double_accuracy else "float32"
# ----- Stencil 2D - created by specifying weights in nested list --------------------------
- src, dst = ps.fields("src, src_tmp: [2D]", layout='fzyx')
+ src, dst = ps.fields(f"src, src_tmp: {field_type}[2D]", layout='fzyx')
stencil = [[1.11, 2.22, 3.33],
[4.44, 5.55, 6.66],
[7.77, 8.88, 9.99]]
@@ -13,7 +14,7 @@ with CodeGeneration() as ctx:
generate_sweep(ctx, 'CudaJacobiKernel2D', assignments, field_swaps=[(src, dst)], target=ps.Target.GPU)
# ----- Stencil 3D - created by using kernel_decorator with assignments in '@=' format -----
- src, dst = ps.fields("src, src_tmp: [3D]")
+ src, dst = ps.fields(f"src, src_tmp: {field_type}[3D]")
@ps.kernel
def kernel_func():
diff --git a/tests/cuda/codegen/CudaPoisson.py b/tests/cuda/codegen/CudaPoisson.py
index eedff3609..24e471a8d 100644
--- a/tests/cuda/codegen/CudaPoisson.py
+++ b/tests/cuda/codegen/CudaPoisson.py
@@ -4,10 +4,11 @@ from pystencils_walberla import CodeGeneration, generate_sweep
with CodeGeneration() as ctx:
+ field_type = "float64" if ctx.double_accuracy else "float32"
# ----- Solving the 2D Poisson equation with rhs --------------------------
dx = sp.Symbol("dx")
dy = sp.Symbol("dy")
- src, dst, rhs = ps.fields("src, src_tmp, rhs: [2D]", layout='fzyx')
+ src, dst, rhs = ps.fields(f"src, src_tmp, rhs: {field_type}[2D]", layout='fzyx')
@ps.kernel
def kernel_func():
diff --git a/tests/cuda/codegen/GeneratedFieldPackInfoTestGPU.cpp b/tests/cuda/codegen/GeneratedFieldPackInfoTestGPU.cpp
index 4c4bed8b4..47c7d6741 100644
--- a/tests/cuda/codegen/GeneratedFieldPackInfoTestGPU.cpp
+++ b/tests/cuda/codegen/GeneratedFieldPackInfoTestGPU.cpp
@@ -202,7 +202,7 @@ int main(int argc, char **argv) {
testScalarFieldPullReduction( blocks, scalarGPUFieldId );
- return 0;
+ return EXIT_SUCCESS;
}
diff --git a/tests/cuda/codegen/GeneratedFieldPackInfoTestGPU.py b/tests/cuda/codegen/GeneratedFieldPackInfoTestGPU.py
index bac0b89da..bd0635e2c 100644
--- a/tests/cuda/codegen/GeneratedFieldPackInfoTestGPU.py
+++ b/tests/cuda/codegen/GeneratedFieldPackInfoTestGPU.py
@@ -3,12 +3,10 @@ import pystencils as ps
from pystencils_walberla import CodeGeneration, generate_pack_info_for_field
with CodeGeneration() as ctx:
-
layout = 'fzyx'
-
field = ps.fields("field: int32[3D]", layout=layout)
# communication
generate_pack_info_for_field(ctx, 'ScalarFieldCommunicationGPU', field, target=ps.Target.GPU)
- generate_pack_info_for_field(ctx, 'ScalarFieldPullReductionGPU', field, target=ps.Target.GPU, operator=op.add,
- gl_to_inner=True)
+ generate_pack_info_for_field(ctx, 'ScalarFieldPullReductionGPU', field, target=ps.Target.GPU,
+ operator=op.add, gl_to_inner=True)
diff --git a/tests/cuda/codegen/MicroBenchmarkGpuLbm.cpp b/tests/cuda/codegen/MicroBenchmarkGpuLbm.cpp
index 8516c0933..8b03e0277 100644
--- a/tests/cuda/codegen/MicroBenchmarkGpuLbm.cpp
+++ b/tests/cuda/codegen/MicroBenchmarkGpuLbm.cpp
@@ -44,8 +44,8 @@ int main( int argc, char **argv )
shared_ptr<StructuredBlockForest> blocks = blockforest::createUniformBlockGrid(1u, 1u, 1u,
128u, 128u, 128u, 1.0, false, false, false, false);
- BlockDataID srcID = cuda::addGPUFieldToStorage<cuda::GPUField<double> >(blocks, "src", 19, field::fzyx, 1);
- BlockDataID dstID = cuda::addGPUFieldToStorage<cuda::GPUField<double> >(blocks, "dst", 19, field::fzyx, 1);
+ BlockDataID srcID = cuda::addGPUFieldToStorage<cuda::GPUField<real_t> >(blocks, "src", 19, field::fzyx, 1);
+ BlockDataID dstID = cuda::addGPUFieldToStorage<cuda::GPUField<real_t> >(blocks, "dst", 19, field::fzyx, 1);
int iterations = 3;
@@ -60,7 +60,7 @@ int main( int argc, char **argv )
for( auto &block: *blocks )
stream( &block );
- WALBERLA_CUDA_CHECK(cudaDeviceSynchronize());
+ WALBERLA_CUDA_CHECK(cudaDeviceSynchronize())
- return 0;
+ return EXIT_SUCCESS;
}
diff --git a/tests/cuda/codegen/MicroBenchmarkGpuLbm.py b/tests/cuda/codegen/MicroBenchmarkGpuLbm.py
index c461df018..09cefc055 100644
--- a/tests/cuda/codegen/MicroBenchmarkGpuLbm.py
+++ b/tests/cuda/codegen/MicroBenchmarkGpuLbm.py
@@ -4,18 +4,18 @@ from lbmpy import LBStencil, Stencil
from pystencils_walberla import CodeGeneration, generate_sweep
with CodeGeneration() as ctx:
- f_size = 19
- dtype = 'float64' if ctx.double_accuracy else 'float32'
+ stencil = LBStencil(Stencil.D3Q19)
+ field_type = 'float64' if ctx.double_accuracy else 'float32'
# Copy sweep
- src, dst = ps.fields(f"src({f_size}), dst({f_size}) : {dtype}[3D]", layout='fzyx')
+ src, dst = ps.fields(f"src({stencil.Q}), dst({stencil.Q}) : {field_type}[3D]", layout='fzyx')
- copy_only = [ps.Assignment(dst(i), src(i)) for i in range(f_size)]
+ copy_only = [ps.Assignment(dst(i), src(i)) for i in range(stencil.Q)]
generate_sweep(ctx, 'MicroBenchmarkCopyKernel', copy_only,
target=ps.Target.GPU, gpu_indexing_params={'block_size': (128, 1, 1)})
# Stream-only sweep
- stencil = LBStencil(Stencil.D3Q19)
+
stream_only = create_stream_only_kernel(stencil, src_field=src, dst_field=dst)
generate_sweep(ctx, 'MicroBenchmarkStreamKernel', stream_only,
target=ps.Target.GPU, gpu_indexing_params={'block_size': (128, 1, 1)})
diff --git a/tests/field/codegen/CodegenJacobiCPU.cpp b/tests/field/codegen/CodegenJacobiCPU.cpp
index d0a9693ae..3bba9623e 100644
--- a/tests/field/codegen/CodegenJacobiCPU.cpp
+++ b/tests/field/codegen/CodegenJacobiCPU.cpp
@@ -38,7 +38,7 @@
using namespace walberla;
-typedef GhostLayerField<double,1> ScalarField;
+typedef GhostLayerField<real_t, 1> ScalarField;
void testJacobi2D()
@@ -49,12 +49,12 @@ void testJacobi2D()
shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid (
uint_t(1) , uint_t(1), uint_t(1), // number of blocks in x,y,z direction
xSize, ySize, uint_t(1), // how many cells per block (x,y,z)
- real_t(1), // dx: length of one cell in physical coordinates
+ real_c(1.0), // dx: length of one cell in physical coordinates
false, // one block per process - "false" means all blocks to one process
true, true, true ); // full periodicity
- BlockDataID fieldID = field::addToStorage<ScalarField>(blocks, "Field", real_t(0.0));
+ BlockDataID fieldID = field::addToStorage<ScalarField>(blocks, "Field", real_c(0.0));
// Initialize a quarter of the field with ones, the rest remains 0
// Jacobi averages the domain -> every cell should be at 0.25 at sufficiently many timesteps
@@ -63,7 +63,7 @@ void testJacobi2D()
auto f = blockIt->getData<ScalarField>( fieldID );
for( cell_idx_t y = 0; y < cell_idx_c( f->ySize() / 2 ); ++y )
for( cell_idx_t x = 0; x < cell_idx_c( f->xSize() / 2 ); ++x )
- f->get( x, y, 0 ) = 1.0;
+ f->get( x, y, 0 ) = real_c(1.0);
}
typedef blockforest::communication::UniformBufferedScheme<stencil::D2Q9> CommScheme;
@@ -84,7 +84,7 @@ void testJacobi2D()
auto firstBlock = blocks->begin();
auto f = firstBlock->getData<ScalarField>( fieldID );
- WALBERLA_CHECK_FLOAT_EQUAL(f->get(0,0,0), real_t(1.0 / 4.0));
+ WALBERLA_CHECK_FLOAT_EQUAL(f->get(0,0,0), real_c(1.0 / 4.0));
}
@@ -97,12 +97,12 @@ void testJacobi3D()
shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid (
uint_t(1) , uint_t(1), uint_t(1), // number of blocks in x,y,z direction
xSize, ySize, zSize, // how many cells per block (x,y,z)
- real_t(1), // dx: length of one cell in physical coordinates
+ real_c(1.0), // dx: length of one cell in physical coordinates
false, // one block per process - "false" means all blocks to one process
true, true, true ); // no periodicity
- BlockDataID fieldID = field::addToStorage<ScalarField>(blocks, "Field", real_t(0.0));
+ BlockDataID fieldID = field::addToStorage<ScalarField>(blocks, "Field", real_c(0.0));
// Initialize a quarter of the field with ones, the rest remains 0
// Jacobi averages the domain -> every cell should be at 0.25 at sufficiently many timesteps
@@ -112,7 +112,7 @@ void testJacobi3D()
for( cell_idx_t z = 0; z < cell_idx_c( f->zSize() / 2); ++z )
for( cell_idx_t y = 0; y < cell_idx_c( f->ySize() / 2 ); ++y )
for( cell_idx_t x = 0; x < cell_idx_c( f->xSize() / 2 ); ++x )
- f->get( x, y, z ) = 1.0;
+ f->get( x, y, z ) = real_c(1.0);
}
typedef blockforest::communication::UniformBufferedScheme<stencil::D3Q7> CommScheme;
@@ -132,7 +132,7 @@ void testJacobi3D()
auto firstBlock = blocks->begin();
auto f = firstBlock->getData<ScalarField>( fieldID );
- WALBERLA_CHECK_FLOAT_EQUAL(f->get(0,0,0), real_t(1.0 / 8.0));
+ WALBERLA_CHECK_FLOAT_EQUAL(f->get(0,0,0), real_c(1.0 / 8.0));
}
@@ -144,5 +144,5 @@ int main( int argc, char ** argv )
testJacobi2D();
testJacobi3D();
- return 0;
+ return EXIT_SUCCESS;
}
diff --git a/tests/field/codegen/CodegenPoissonCPU.cpp b/tests/field/codegen/CodegenPoissonCPU.cpp
index 4ecbab91e..6c5696d40 100644
--- a/tests/field/codegen/CodegenPoissonCPU.cpp
+++ b/tests/field/codegen/CodegenPoissonCPU.cpp
@@ -67,7 +67,7 @@ void initF( const shared_ptr< StructuredBlockStorage > & blocks, const BlockData
CellInterval xyz = f->xyzSize();
for( auto cell = xyz.begin(); cell != xyz.end(); ++cell )
{
- f->get( *cell ) = 0.0;
+ f->get( *cell ) = real_c(0.0);
}
}
}
@@ -77,25 +77,25 @@ void testPoisson()
const uint_t xCells = uint_t(100);
const uint_t yCells = uint_t(100);
- const real_t xSize = real_t(1);
- const real_t ySize = real_t(1);
+ const real_t xSize = real_c(1.0);
+ const real_t ySize = real_c(1.0);
const real_t dx = xSize / real_c( xCells + uint_t(1) );
const real_t dy = ySize / real_c( yCells + uint_t(1) );
// Create blocks
shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid (
- math::AABB( real_t(0.5) * dx, real_t(0.5) * dy, real_t(0),
- xSize - real_t(0.5) * dx, ySize - real_t(0.5) * dy, dx ),
+ math::AABB( real_c(0.5) * dx, real_c(0.5) * dy, real_c(0.0),
+ xSize - real_c(0.5) * dx, ySize - real_c(0.5) * dy, dx ),
uint_t(1) , uint_t(1), uint_t(1), // number of blocks in x,y,z direction
xCells, yCells, uint_t(1), // how many cells per block (x,y,z)
false, // one block per process - "false" means all blocks to one process
false, false, false ); // no periodicity
- BlockDataID fieldID = field::addToStorage<ScalarField_T>(blocks, "Field", real_t(0.0));
+ BlockDataID fieldID = field::addToStorage<ScalarField_T>(blocks, "Field", real_c(0.0));
initU( blocks, fieldID );
- BlockDataID fId = field::addToStorage< ScalarField_T >( blocks, "f", real_t(0.0));
+ BlockDataID fId = field::addToStorage< ScalarField_T >( blocks, "f", real_c(0.0));
initF( blocks, fId );
typedef blockforest::communication::UniformBufferedScheme<stencil::D2Q9> CommScheme;
@@ -115,7 +115,7 @@ void testPoisson()
auto firstBlock = blocks->begin();
auto f = firstBlock->getData<ScalarField_T>( fieldID );
- WALBERLA_CHECK_LESS(f->get(50,99,0) - std::sin( math::pi * 0.5 ) * std::sinh( math::pi * 0.99 ), 0.01);
+ WALBERLA_CHECK_LESS(f->get(50,99,0) - std::sin( math::pi * 0.5 ) * std::sinh( math::pi * 0.99 ), 0.01)
}
diff --git a/tests/field/codegen/EK.cpp b/tests/field/codegen/EK.cpp
index dbfdcb32d..1f795b873 100644
--- a/tests/field/codegen/EK.cpp
+++ b/tests/field/codegen/EK.cpp
@@ -56,7 +56,7 @@ void initC(const shared_ptr< StructuredBlockStorage > & blocks, BlockDataID cID)
void checkC(const shared_ptr< StructuredBlockStorage > & blocks, BlockDataID cID, real_t D, uint_t t)
{
- real_t sum(0);
+ real_t sum(0.0);
uint_t size(0);
for(auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
{
@@ -76,12 +76,12 @@ void checkC(const shared_ptr< StructuredBlockStorage > & blocks, BlockDataID cID
const real_t val = c->get(x, y, 0);
sum += val;
- const real_t x0 = real_c(x) - real_c(c->xSize()/2) + real_t(0.5);
- const real_t y0 = real_c(y) - real_c(c->ySize()/2) + real_t(0.5);
+ const real_t x0 = real_c(x) - real_c(c->xSize()/2) + real_c(0.5);
+ const real_t y0 = real_c(y) - real_c(c->ySize()/2) + real_c(0.5);
const real_t r2 = x0*x0 + y0*y0;
// solution to the diffusion equation in 2D
- const real_t ref = real_t(size)/(real_t(4)*math::pi*D*real_c(t)) * std::exp(-r2/real_c(4*D*real_c(t)));
+ const real_t ref = real_c(size)/(real_c(4.0)*math::pi*D*real_c(t)) * std::exp(-r2/real_c(4*D*real_c(t)));
real_t rel = std::abs(real_c(1) - val/ref);
if(ref >= 1) // we only expect good agreement where the values are larger than sum/size
@@ -91,12 +91,12 @@ void checkC(const shared_ptr< StructuredBlockStorage > & blocks, BlockDataID cID
std::cout << y << " " << y0 << " " << val << " " << ref << " " << (rel*100) << "%" << std::endl;
}
- WALBERLA_CHECK_LESS(rel, 0.03); // 3% deviation is okay
+ WALBERLA_CHECK_LESS(rel, 0.03) // 3% deviation is okay
}
}
}
}
- WALBERLA_CHECK_FLOAT_EQUAL(sum, real_c(size));
+ WALBERLA_CHECK_FLOAT_EQUAL(sum, real_c(size))
}
void testEK()
@@ -109,13 +109,13 @@ void testEK()
shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid (
uint_t(1) , uint_t(1), uint_t(1), // number of blocks in x,y,z direction
xSize, ySize, uint_t(1), // how many cells per block (x,y,z)
- real_t(1), // dx: length of one cell in physical coordinates
+ real_c(1.0), // dx: length of one cell in physical coordinates
false, // one block per process - "false" means all blocks to one process
true, true, true ); // full periodicity
- BlockDataID c = field::addToStorage<DensityField_T>(blocks, "c", real_t(0.0), field::fzyx);
+ BlockDataID c = field::addToStorage<DensityField_T>(blocks, "c", real_c(0.0), field::fzyx);
initC(blocks, c);
- BlockDataID j = field::addToStorage<FluxField_T>(blocks, "j", real_t(0.0), field::fzyx);
+ BlockDataID j = field::addToStorage<FluxField_T>(blocks, "j", real_c(0.0), field::fzyx);
typedef blockforest::communication::UniformBufferedScheme<stencil::D2Q9> CommScheme;
typedef field::communication::PackInfo<DensityField_T> Packing;
@@ -127,7 +127,7 @@ void testEK()
// Registering the sweep
timeloop.add() << BeforeFunction( commScheme, "Communication" )
- << Sweep( pystencils::EKFlux(D, c, j), "EK flux Kernel" );
+ << Sweep( pystencils::EKFlux(c, j, D), "EK flux Kernel" );
timeloop.add() << Sweep( pystencils::EKContinuity(c, j), "EK continuity Kernel" );
timeloop.run();
diff --git a/tests/field/codegen/EK.py b/tests/field/codegen/EK.py
index 7cf53c035..79985ece0 100644
--- a/tests/field/codegen/EK.py
+++ b/tests/field/codegen/EK.py
@@ -2,9 +2,11 @@ import sympy as sp
import pystencils as ps
from pystencils_walberla import CodeGeneration, generate_sweep
+
def grad(f):
return sp.Matrix([ps.fd.diff(f, i) for i in range(f.spatial_dimensions)])
+
D = sp.Symbol("D")
with CodeGeneration() as ctx:
diff --git a/tests/field/codegen/GeneratedFieldPackInfoTest.cpp b/tests/field/codegen/GeneratedFieldPackInfoTest.cpp
index 2f8031064..a4835531b 100644
--- a/tests/field/codegen/GeneratedFieldPackInfoTest.cpp
+++ b/tests/field/codegen/GeneratedFieldPackInfoTest.cpp
@@ -47,9 +47,9 @@ void testScalarField( IBlock * block, BlockDataID fieldId )
GhostLayerField<int,1> & field = *(block->getData<GhostLayerField<int,1> > (fieldId));
field.setWithGhostLayer( 0 );
- WALBERLA_CHECK_EQUAL(field.xSize(), 2);
- WALBERLA_CHECK_EQUAL(field.ySize(), 2);
- WALBERLA_CHECK_EQUAL(field.zSize(), 2);
+ WALBERLA_CHECK_EQUAL(field.xSize(), 2)
+ WALBERLA_CHECK_EQUAL(field.ySize(), 2)
+ WALBERLA_CHECK_EQUAL(field.zSize(), 2)
// initialize the bottom boundary
field(0,0,0) = 1;
@@ -63,10 +63,10 @@ void testScalarField( IBlock * block, BlockDataID fieldId )
pystencils::ScalarFieldCommunication pi(fieldId);
pi.communicateLocal( block, block, stencil::B );
- WALBERLA_CHECK_EQUAL ( field(0,0,+2), 1 );
- WALBERLA_CHECK_EQUAL ( field(0,1,+2), 2 );
- WALBERLA_CHECK_EQUAL ( field(1,0,+2), 3 );
- WALBERLA_CHECK_EQUAL ( field(1,1,+2), 4 );
+ WALBERLA_CHECK_EQUAL ( field(0,0,+2), 1 )
+ WALBERLA_CHECK_EQUAL ( field(0,1,+2), 2 )
+ WALBERLA_CHECK_EQUAL ( field(1,0,+2), 3 )
+ WALBERLA_CHECK_EQUAL ( field(1,1,+2), 4 )
// -------------- Buffer Communication Test ---------------------
@@ -86,10 +86,10 @@ void testScalarField( IBlock * block, BlockDataID fieldId )
pi.unpackData( block, stencil::T, recvBuf );
- WALBERLA_CHECK_EQUAL ( field(0,0,+2), 1 );
- WALBERLA_CHECK_EQUAL ( field(0,1,+2), 2 );
- WALBERLA_CHECK_EQUAL ( field(1,0,+2), 3 );
- WALBERLA_CHECK_EQUAL ( field(1,1,+2), 4 );
+ WALBERLA_CHECK_EQUAL ( field(0,0,+2), 1 )
+ WALBERLA_CHECK_EQUAL ( field(0,1,+2), 2 )
+ WALBERLA_CHECK_EQUAL ( field(1,0,+2), 3 )
+ WALBERLA_CHECK_EQUAL ( field(1,1,+2), 4 )
}
void testScalarFieldPullReduction( IBlock * block, BlockDataID fieldId )
@@ -97,9 +97,9 @@ void testScalarFieldPullReduction( IBlock * block, BlockDataID fieldId )
GhostLayerField<int,1> & field = *(block->getData<GhostLayerField<int,1> > (fieldId));
field.setWithGhostLayer( 0 );
- WALBERLA_CHECK_EQUAL(field.xSize(), 2);
- WALBERLA_CHECK_EQUAL(field.ySize(), 2);
- WALBERLA_CHECK_EQUAL(field.zSize(), 2);
+ WALBERLA_CHECK_EQUAL(field.xSize(), 2)
+ WALBERLA_CHECK_EQUAL(field.ySize(), 2)
+ WALBERLA_CHECK_EQUAL(field.zSize(), 2)
// initialize the bottom ghost layer cells
field(0,0,-1) = 1;
@@ -118,32 +118,32 @@ void testScalarFieldPullReduction( IBlock * block, BlockDataID fieldId )
pi1.communicateLocal( block, block, stencil::B );
// check values in top ghost layer
- WALBERLA_CHECK_EQUAL ( field(0,0,2), 0 );
- WALBERLA_CHECK_EQUAL ( field(0,1,2), 0 );
- WALBERLA_CHECK_EQUAL ( field(1,0,2), 0 );
- WALBERLA_CHECK_EQUAL ( field(1,1,2), 0 );
+ WALBERLA_CHECK_EQUAL ( field(0,0,2), 0 )
+ WALBERLA_CHECK_EQUAL ( field(0,1,2), 0 )
+ WALBERLA_CHECK_EQUAL ( field(1,0,2), 0 )
+ WALBERLA_CHECK_EQUAL ( field(1,1,2), 0 )
// check values in top interior cells
- WALBERLA_CHECK_EQUAL ( field(0,0,1), 2 );
- WALBERLA_CHECK_EQUAL ( field(0,1,1), 3 );
- WALBERLA_CHECK_EQUAL ( field(1,0,1), 4 );
- WALBERLA_CHECK_EQUAL ( field(1,1,1), 5 );
+ WALBERLA_CHECK_EQUAL ( field(0,0,1), 2 )
+ WALBERLA_CHECK_EQUAL ( field(0,1,1), 3 )
+ WALBERLA_CHECK_EQUAL ( field(1,0,1), 4 )
+ WALBERLA_CHECK_EQUAL ( field(1,1,1), 5 )
// communicate periodic from top to bottom with standard form to sync ghost layers
pystencils::ScalarFieldCommunication pi2 (fieldId);
pi2.communicateLocal( block, block, stencil::T );
// check values in bottom ghost layer
- WALBERLA_CHECK_EQUAL ( field(0,0,-1), 2 );
- WALBERLA_CHECK_EQUAL ( field(0,1,-1), 3 );
- WALBERLA_CHECK_EQUAL ( field(1,0,-1), 4 );
- WALBERLA_CHECK_EQUAL ( field(1,1,-1), 5 );
+ WALBERLA_CHECK_EQUAL ( field(0,0,-1), 2 )
+ WALBERLA_CHECK_EQUAL ( field(0,1,-1), 3 )
+ WALBERLA_CHECK_EQUAL ( field(1,0,-1), 4 )
+ WALBERLA_CHECK_EQUAL ( field(1,1,-1), 5 )
// check values in top interior cells
- WALBERLA_CHECK_EQUAL ( field(0,0,1), 2 );
- WALBERLA_CHECK_EQUAL ( field(0,1,1), 3 );
- WALBERLA_CHECK_EQUAL ( field(1,0,1), 4 );
- WALBERLA_CHECK_EQUAL ( field(1,1,1), 5 );
+ WALBERLA_CHECK_EQUAL ( field(0,0,1), 2 )
+ WALBERLA_CHECK_EQUAL ( field(0,1,1), 3 )
+ WALBERLA_CHECK_EQUAL ( field(1,0,1), 4 )
+ WALBERLA_CHECK_EQUAL ( field(1,1,1), 5 )
}
@@ -159,7 +159,7 @@ int main(int argc, char **argv)
uint_t processes = uint_c( MPIManager::instance()->numProcesses() );
auto blocks = createUniformBlockGrid(processes,1 ,1, //blocks
2,2,2, //cells
- 1, //dx
+ real_c(1.0), //dx
false, //one block per process
true,true,true);//periodicity
@@ -172,7 +172,7 @@ int main(int argc, char **argv)
// Create a BlockForest with 8x8x8 cells per block
blocks = createUniformBlockGrid(processes,1 ,1, //blocks
8,8,8, //cells
- 1, //dx
+ real_c(1.0), //dx
false, //one block per process
true,true,true);//periodicity
@@ -188,7 +188,7 @@ int main(int argc, char **argv)
for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) // block loop
testScalarFieldPullReduction( &(*blockIt), scalarFieldId );
- return 0;
+ return EXIT_SUCCESS;
}
} // namespace walberla
diff --git a/tests/field/codegen/GeneratedFieldPackInfoTest.py b/tests/field/codegen/GeneratedFieldPackInfoTest.py
index 8124060c3..91ef712d8 100644
--- a/tests/field/codegen/GeneratedFieldPackInfoTest.py
+++ b/tests/field/codegen/GeneratedFieldPackInfoTest.py
@@ -3,12 +3,10 @@ import pystencils as ps
from pystencils_walberla import CodeGeneration, generate_pack_info_for_field
with CodeGeneration() as ctx:
-
layout = 'fzyx'
-
field = ps.fields("field: int32[3D]", layout=layout)
# communication
generate_pack_info_for_field(ctx, 'ScalarFieldCommunication', field, target=ps.Target.CPU)
- generate_pack_info_for_field(ctx, 'ScalarFieldPullReduction', field, target=ps.Target.CPU, operator=op.add,
- gl_to_inner=True)
+ generate_pack_info_for_field(ctx, 'ScalarFieldPullReduction', field, target=ps.Target.CPU,
+ operator=op.add, gl_to_inner=True)
diff --git a/tests/field/codegen/JacobiKernel.py b/tests/field/codegen/JacobiKernel.py
index 29e677a72..ba39b1b56 100644
--- a/tests/field/codegen/JacobiKernel.py
+++ b/tests/field/codegen/JacobiKernel.py
@@ -4,8 +4,9 @@ from pystencils_walberla import CodeGeneration, generate_sweep
with CodeGeneration() as ctx:
+ field_type = "float64" if ctx.double_accuracy else "float32"
# ----- Stencil 2D - created by specifying weights in nested list --------------------------
- src, dst = ps.fields("src, src_tmp: [2D]", layout='fzyx')
+ src, dst = ps.fields(f"src, src_tmp: {field_type}[2D]", layout='fzyx')
stencil = [[1.11, 2.22, 3.33],
[4.44, 5.55, 6.66],
[7.77, 8.88, 9.99]]
@@ -13,7 +14,7 @@ with CodeGeneration() as ctx:
generate_sweep(ctx, 'JacobiKernel2D', assignments, field_swaps=[(src, dst)])
# ----- Stencil 3D - created by using kernel_decorator with assignments in '@=' format -----
- src, dst = ps.fields("src, src_tmp: [3D]", layout='fzyx')
+ src, dst = ps.fields(f"src, src_tmp: {field_type}[3D]", layout='fzyx')
@ps.kernel
def kernel_func():
diff --git a/tests/field/codegen/MultipleFieldSwaps.cpp b/tests/field/codegen/MultipleFieldSwaps.cpp
index 965bf15a3..5fa7eec87 100644
--- a/tests/field/codegen/MultipleFieldSwaps.cpp
+++ b/tests/field/codegen/MultipleFieldSwaps.cpp
@@ -32,7 +32,7 @@
using namespace walberla;
-typedef GhostLayerField<double,1> ScalarField;
+typedef GhostLayerField<real_t ,1> ScalarField;
void testMultipleFieldSwaps()
{
uint_t xSize = 5;
@@ -41,14 +41,14 @@ void testMultipleFieldSwaps()
shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid (
uint_t(1) , uint_t(1), uint_t(1), // number of blocks in x,y,z direction
xSize, ySize, uint_t(1), // how many cells per block (x,y,z)
- real_t(1), // dx: length of one cell in physical coordinates
+ real_c(1.0), // dx: length of one cell in physical coordinates
false, // one block per process - "false" means all blocks to one process
true, true, true ); // full periodicity
- BlockDataID fieldID_1 = field::addToStorage<ScalarField>(blocks, "Field_1", real_t(1.0));
- BlockDataID fieldID_2 = field::addToStorage<ScalarField>(blocks, "Field_2", real_t(1.0));
- BlockDataID fieldID_3 = field::addToStorage<ScalarField>(blocks, "Field_3", real_t(1.0));
+ BlockDataID fieldID_1 = field::addToStorage<ScalarField>(blocks, "Field_1", real_c(1.0));
+ BlockDataID fieldID_2 = field::addToStorage<ScalarField>(blocks, "Field_2", real_c(1.0));
+ BlockDataID fieldID_3 = field::addToStorage<ScalarField>(blocks, "Field_3", real_c(1.0));
pystencils::MultipleFieldSwaps kernel(fieldID_1, fieldID_2, fieldID_3);
@@ -63,9 +63,9 @@ void testMultipleFieldSwaps()
// clang-format off
WALBERLA_FOR_ALL_CELLS_XYZ(field_1, Cell globalCell;
blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
- WALBERLA_CHECK_FLOAT_EQUAL(field_1->get(x, y, z), real_t(2.0))
- WALBERLA_CHECK_FLOAT_EQUAL(field_2->get(x, y, z), real_t(2.0))
- WALBERLA_CHECK_FLOAT_EQUAL(field_3->get(x, y, z), real_t(2.0))
+ WALBERLA_CHECK_FLOAT_EQUAL(field_1->get(x, y, z), real_c(2.0))
+ WALBERLA_CHECK_FLOAT_EQUAL(field_2->get(x, y, z), real_c(2.0))
+ WALBERLA_CHECK_FLOAT_EQUAL(field_3->get(x, y, z), real_c(2.0))
)
// clang-format on
}
diff --git a/tests/field/codegen/MultipleFieldSwaps.py b/tests/field/codegen/MultipleFieldSwaps.py
index d4268d8ea..c8abbfc9b 100644
--- a/tests/field/codegen/MultipleFieldSwaps.py
+++ b/tests/field/codegen/MultipleFieldSwaps.py
@@ -2,9 +2,10 @@ import pystencils as ps
from pystencils_walberla import CodeGeneration, generate_sweep
with CodeGeneration() as ctx:
- src_1, dst_1 = ps.fields("src1, src1_tmp: [2D]", layout='fzyx')
- src_2, dst_2 = ps.fields("src2, src2_tmp: [2D]", layout='fzyx')
- src_3, dst_3 = ps.fields("src3, src3_tmp: [2D]", layout='fzyx')
+ field_type = "float64" if ctx.double_accuracy else "float32"
+ src_1, dst_1 = ps.fields(f"src1, src1_tmp: {field_type}[2D]", layout='fzyx')
+ src_2, dst_2 = ps.fields(f"src2, src2_tmp: {field_type}[2D]", layout='fzyx')
+ src_3, dst_3 = ps.fields(f"src3, src3_tmp: {field_type}[2D]", layout='fzyx')
assignments = ps.AssignmentCollection([ps.Assignment(dst_1.center, 2 * src_1.center),
ps.Assignment(dst_2.center, 2 * src_2.center),
diff --git a/tests/field/codegen/Poisson.py b/tests/field/codegen/Poisson.py
index a4da698dd..a4202ad5a 100644
--- a/tests/field/codegen/Poisson.py
+++ b/tests/field/codegen/Poisson.py
@@ -4,10 +4,11 @@ from pystencils_walberla import CodeGeneration, generate_sweep
with CodeGeneration() as ctx:
+ field_type = "float64" if ctx.double_accuracy else "float32"
# ----- Solving the 2D Poisson equation with rhs --------------------------
dx = sp.Symbol("dx")
dy = sp.Symbol("dy")
- src, dst, rhs = ps.fields("src, src_tmp, rhs: [2D]", layout='fzyx')
+ src, dst, rhs = ps.fields(f"src, src_tmp, rhs: {field_type}[2D]", layout='fzyx')
@ps.kernel
def kernel_func():
diff --git a/tests/lbm/codegen/FieldLayoutAndVectorizationTest.cpp b/tests/lbm/codegen/FieldLayoutAndVectorizationTest.cpp
index db86a8d31..b1a435546 100644
--- a/tests/lbm/codegen/FieldLayoutAndVectorizationTest.cpp
+++ b/tests/lbm/codegen/FieldLayoutAndVectorizationTest.cpp
@@ -67,18 +67,18 @@ int main(int argc, char **argv) {
walberla::Environment walberlaEnv(argc, argv);
auto blocks = blockforest::createUniformBlockGrid( 1, 1, 1,
- 32, 32, 32, real_t(1),
+ 32, 32, 32, real_c(1.0),
0, false, false,
true, true, true, //periodicity
false );
- real_t omega = real_t(1.6);
- Vector3<real_t> initialVel(real_t(0.1), real_t(0), real_t(0));
+ real_t omega = real_c(1.6);
+ Vector3<real_t> initialVel(real_c(0.1), real_c(0.0), real_c(0.0));
using LatticeModel1_T = lbm::FieldLayoutAndVectorizationTest_FZYX_Vec_LatticeModel;
LatticeModel1_T latticeModel1(omega);
- BlockDataID pdfField1ID = lbm::addPdfFieldToStorage(blocks, "pdf field1", latticeModel1, initialVel, real_t(1), field::fzyx);
+ BlockDataID pdfField1ID = lbm::addPdfFieldToStorage(blocks, "pdf field1", latticeModel1, initialVel, real_c(1.0), field::fzyx);
for(auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt) {
auto sweep = LatticeModel1_T::Sweep(pdfField1ID);
sweep(&(*blockIt));
@@ -86,7 +86,7 @@ int main(int argc, char **argv) {
using LatticeModel2_T = lbm::FieldLayoutAndVectorizationTest_FZYX_NoVec_LatticeModel;
LatticeModel2_T latticeModel2(omega);
- BlockDataID pdfField2ID = lbm::addPdfFieldToStorage(blocks, "pdf field2", latticeModel2, initialVel, real_t(1), field::fzyx);
+ BlockDataID pdfField2ID = lbm::addPdfFieldToStorage(blocks, "pdf field2", latticeModel2, initialVel, real_c(1.0), field::fzyx);
for(auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt) {
auto sweep = LatticeModel2_T::Sweep(pdfField2ID);
sweep(&(*blockIt));
@@ -96,7 +96,7 @@ int main(int argc, char **argv) {
using LatticeModel3_T = lbm::FieldLayoutAndVectorizationTest_ZYXF_Vec_LatticeModel;
LatticeModel3_T latticeModel3(omega);
- BlockDataID pdfField3ID = lbm::addPdfFieldToStorage(blocks, "pdf field3", latticeModel3, initialVel, real_t(1), field::zyxf);
+ BlockDataID pdfField3ID = lbm::addPdfFieldToStorage(blocks, "pdf field3", latticeModel3, initialVel, real_c(1.0), field::zyxf);
for(auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt) {
auto sweep = LatticeModel3_T::Sweep(pdfField3ID);
sweep(&(*blockIt));
@@ -106,7 +106,7 @@ int main(int argc, char **argv) {
using LatticeModel4_T = lbm::FieldLayoutAndVectorizationTest_ZYXF_NoVec_LatticeModel;
LatticeModel4_T latticeModel4(omega);
- BlockDataID pdfField4ID = lbm::addPdfFieldToStorage(blocks, "pdf field4", latticeModel4, initialVel, real_t(1), field::zyxf);
+ BlockDataID pdfField4ID = lbm::addPdfFieldToStorage(blocks, "pdf field4", latticeModel4, initialVel, real_c(1.0), field::zyxf);
for(auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt) {
auto sweep = LatticeModel4_T::Sweep(pdfField4ID);
sweep(&(*blockIt));
diff --git a/tests/lbm/codegen/FieldLayoutAndVectorizationTest.py b/tests/lbm/codegen/FieldLayoutAndVectorizationTest.py
index f04e357bc..e639e0a60 100644
--- a/tests/lbm/codegen/FieldLayoutAndVectorizationTest.py
+++ b/tests/lbm/codegen/FieldLayoutAndVectorizationTest.py
@@ -10,7 +10,8 @@ from collections import namedtuple
with CodeGeneration() as ctx:
omega_shear = sp.symbols("omega")
- lbm_config = LBMConfig(stencil=LBStencil(Stencil.D2Q9), compressible=False, method=Method.SRT)
+ lbm_config = LBMConfig(stencil=LBStencil(Stencil.D2Q9), compressible=False, method=Method.SRT,
+ zero_centered=False, delta_equilibrium=False)
collision_rule = create_lb_collision_rule(lbm_config=lbm_config)
SetupDefinition = namedtuple('SetupDefinition', ['name', 'field_layout', 'vectorization_dict'])
diff --git a/tests/lbm/codegen/FluctuatingMRT.cpp b/tests/lbm/codegen/FluctuatingMRT.cpp
index a8a425d9f..62e81c943 100644
--- a/tests/lbm/codegen/FluctuatingMRT.cpp
+++ b/tests/lbm/codegen/FluctuatingMRT.cpp
@@ -65,18 +65,18 @@ int main( int argc, char ** argv )
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 real_t temperature = real_t(0.01);
+ const uint_t timesteps = parameters.getParameter< uint_t > ( "timesteps", uint_c(10.0) );
+ const real_t temperature = real_c(0.01);
const uint_t seed = uint_t(0);
- const real_t remainingTimeLoggerFrequency = parameters.getParameter< double >( "remainingTimeLoggerFrequency", 3.0 ); // in seconds
+ const real_t remainingTimeLoggerFrequency = parameters.getParameter< real_t >( "remainingTimeLoggerFrequency", real_c(3.0) ); // in seconds
// create fields
- real_t force = 2E-4; // Force to apply on each node on each axis
+ real_t force = real_c(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, 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 pdfFieldId = lbm::addPdfFieldToStorage( blocks, "pdf field", latticeModel, initialVelocity, real_c(1.0), uint_t(1), field::fzyx );
BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >( blocks, "flag field" );
// create and initialize flag field
@@ -123,15 +123,15 @@ int main( int argc, char ** argv )
WALBERLA_FOR_ALL_CELLS_XYZ_OMP(v, omp critical, {
momentum += rho->get(x, y, z) * v->get(x, y, z);
count++;
- });
+ })
}
// 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);
+ 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.py b/tests/lbm/codegen/FluctuatingMRT.py
index eed48db71..ff9ad6da9 100644
--- a/tests/lbm/codegen/FluctuatingMRT.py
+++ b/tests/lbm/codegen/FluctuatingMRT.py
@@ -1,16 +1,18 @@
import sympy as sp
import pystencils as ps
-from lbmpy.creationfunctions import create_lb_collision_rule, create_mrt_orthogonal
+from lbmpy.creationfunctions import create_lb_collision_rule
from lbmpy.moments import is_bulk_moment, is_shear_moment, get_order
from lbmpy.forcemodels import Guo
-from lbmpy import LBMConfig, LBMOptimisation, LBStencil, Stencil
+from lbmpy import LBMConfig, LBMOptimisation, LBStencil, Stencil, Method
from pystencils_walberla import CodeGeneration
from lbmpy_walberla import generate_lattice_model
with CodeGeneration() as ctx:
+ data_type = "float64" if ctx.double_accuracy else "float32"
+
omega_shear = sp.symbols("omega_shear")
temperature = sp.symbols("temperature")
- force_field, vel_field = ps.fields("force(3), velocity(3): [3D]", layout='fzyx')
+ force_field, vel_field = ps.fields(f"force(3), velocity(3): {data_type}[3D]", layout='fzyx')
def rr_getter(moment_group):
"""Maps a group of moments to a relaxation rate (shear, bulk, even, odd)
@@ -23,7 +25,7 @@ with CodeGeneration() as ctx:
order = order[0]
if order < 2:
- return [0] * len(moment_group)
+ return [0.0] * len(moment_group)
elif any(is_bulk):
assert all(is_bulk)
return [sp.Symbol("omega_bulk")] * len(moment_group)
@@ -36,24 +38,22 @@ with CodeGeneration() as ctx:
else:
return [sp.Symbol("omega_odd")] * len(moment_group)
- method = create_mrt_orthogonal(
- stencil=LBStencil(Stencil.D3Q19),
- compressible=True,
- weighted=True,
- relaxation_rates=rr_getter,
- force_model=Guo(force_field.center_vector)
- )
fluctuating = {'temperature': temperature,
'block_offsets': 'walberla',
'rng_node': ps.rng.PhiloxTwoDoubles if ctx.double_accuracy else ps.rng.PhiloxFourFloats}
- lbm_config = LBMConfig(fluctuating=fluctuating)
+ lbm_config = LBMConfig(stencil=LBStencil(Stencil.D3Q19), method=Method.MRT, compressible=True,
+ weighted=True, zero_centered=False, relaxation_rates=rr_getter,
+ force_model=Guo(force=force_field.center_vector),
+ fluctuating=fluctuating)
+
lbm_opt = LBMOptimisation(cse_global=True)
- collision_rule = create_lb_collision_rule(lb_method=method, lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+ collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
params = {}
if ctx.optimize_for_localhost:
params['cpu_vectorize_info'] = {'assume_inner_stride_one': True, 'assume_aligned': True}
- generate_lattice_model(ctx, 'FluctuatingMRT_LatticeModel', collision_rule, field_layout='fzyx', **params)
+ generate_lattice_model(ctx, 'FluctuatingMRT_LatticeModel', collision_rule, field_layout='fzyx',
+ data_type=data_type, **params)
diff --git a/tests/lbm/codegen/GeneratedFreeSlip.cpp b/tests/lbm/codegen/GeneratedFreeSlip.cpp
index 24cb15ef5..6220ef4cb 100644
--- a/tests/lbm/codegen/GeneratedFreeSlip.cpp
+++ b/tests/lbm/codegen/GeneratedFreeSlip.cpp
@@ -83,13 +83,13 @@ int main(int argc, char** argv)
const real_t omega = parameters.getParameter< real_t >("omega", real_c(1.4));
const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(10));
- const double remainingTimeLoggerFrequency =
- parameters.getParameter< double >("remainingTimeLoggerFrequency", 3.0); // in seconds
+ const real_t remainingTimeLoggerFrequency =
+ parameters.getParameter< real_t >("remainingTimeLoggerFrequency", real_c(3.0)); // in seconds
// create fields
BlockDataID pdfFieldID = blocks->addStructuredBlockData< PdfField_T >(pdfFieldAdder, "PDFs");
- BlockDataID velFieldID = field::addToStorage< VelocityField_T >(blocks, "velocity", real_t(0), field::fzyx);
- BlockDataID densityFieldID = field::addToStorage< ScalarField_T >(blocks, "density", real_t(1.0), field::fzyx);
+ BlockDataID velFieldID = field::addToStorage< VelocityField_T >(blocks, "velocity", real_c(0.0), field::fzyx);
+ BlockDataID densityFieldID = field::addToStorage< ScalarField_T >(blocks, "density", real_c(1.0), field::fzyx);
BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
@@ -129,7 +129,7 @@ int main(int argc, char** argv)
if (tubeSetup)
{
- real_t R = sqrt((globalCell[1] - My) * (globalCell[1] - My) + (globalCell[2] - Mz) * (globalCell[2] - Mz));
+ real_t R = real_c(sqrt((globalCell[1] - My) * (globalCell[1] - My) + (globalCell[2] - Mz) * (globalCell[2] - Mz)));
if (R > real_c(cellsPerBlock[1]) * real_c(0.5)) addFlag(flagField->get(x, y, z), wallFlag);
}
if (slopedWall)
@@ -206,7 +206,8 @@ int main(int argc, char** argv)
real_t v2 = velField->get(cell_idx_c(cellsPerBlock[0] / 2), i + 1, cell_idx_c(cellsPerBlock[2] / 2), 0);
real_t grad = v2 - v1;
// WALBERLA_LOG_DEVEL_VAR(grad)
- WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(grad, 0.0, 1e-16)
+ WALBERLA_CHECK_FLOAT_EQUAL(grad, 0.0)
+
}
}
}
diff --git a/tests/lbm/codegen/GeneratedFreeSlip.py b/tests/lbm/codegen/GeneratedFreeSlip.py
index d0e252ddd..e079cbe75 100644
--- a/tests/lbm/codegen/GeneratedFreeSlip.py
+++ b/tests/lbm/codegen/GeneratedFreeSlip.py
@@ -12,48 +12,51 @@ from lbmpy_walberla import generate_boundary, generate_lb_pack_info, generate_al
import sympy as sp
-stencil = LBStencil(Stencil.D3Q19)
+with CodeGeneration() as ctx:
+ data_type = "float64" if ctx.double_accuracy else "float32"
+ stencil = LBStencil(Stencil.D3Q19)
-pdfs, pdfs_tmp = fields(f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): double[{stencil.D}D]", layout='fzyx')
-velocity_field, density_field = fields(f"velocity({stencil.D}), density(1) : double[{stencil.D}D]", layout='fzyx')
-omega = sp.Symbol("omega")
+ pdfs, pdfs_tmp = fields(f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): {data_type}[{stencil.D}D]",
+ layout='fzyx')
+ velocity_field, density_field = fields(f"velocity({stencil.D}), density(1) : {data_type}[{stencil.D}D]",
+ layout='fzyx')
+ omega = sp.Symbol("omega")
-output = {
- 'density': density_field,
- 'velocity': velocity_field
-}
+ output = {
+ 'density': density_field,
+ 'velocity': velocity_field
+ }
-streaming_pattern = 'esotwist'
-timesteps = get_timesteps(streaming_pattern)
+ streaming_pattern = 'esotwist'
+ timesteps = get_timesteps(streaming_pattern)
-lbm_config = LBMConfig(method=Method.SRT, stencil=stencil, relaxation_rate=omega, force=(1e-5, 0, 0),
- output=output, streaming_pattern=streaming_pattern)
+ lbm_config = LBMConfig(method=Method.SRT, stencil=stencil, relaxation_rate=omega, force=(1e-5, 0, 0),
+ output=output, streaming_pattern=streaming_pattern)
-lbm_opt = LBMOptimisation(symbolic_field=pdfs,
- cse_global=False, cse_pdfs=False)
+ lbm_opt = LBMOptimisation(symbolic_field=pdfs,
+ cse_global=False, cse_pdfs=False)
-method = create_lb_method(lbm_config=lbm_config)
+ method = create_lb_method(lbm_config=lbm_config)
-# getter & setter
-setter_assignments = macroscopic_values_setter(method, density=density_field.center,
- velocity=velocity_field.center_vector,
- pdfs=pdfs, streaming_pattern=streaming_pattern,
- previous_timestep=timesteps[0])
+ # getter & setter
+ setter_assignments = macroscopic_values_setter(method, density=density_field.center,
+ velocity=velocity_field.center_vector,
+ pdfs=pdfs, streaming_pattern=streaming_pattern,
+ previous_timestep=timesteps[0])
-if not is_inplace(streaming_pattern):
- lbm_opt = replace(lbm_opt, symbolic_temporary_field=pdfs_tmp)
- field_swaps = [(pdfs, pdfs_tmp)]
-else:
- field_swaps = []
+ if not is_inplace(streaming_pattern):
+ lbm_opt = replace(lbm_opt, symbolic_temporary_field=pdfs_tmp)
+ field_swaps = [(pdfs, pdfs_tmp)]
+ else:
+ field_swaps = []
-collision_rule = create_lb_collision_rule(lb_method=method, lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+ collision_rule = create_lb_collision_rule(lb_method=method, lbm_config=lbm_config, lbm_optimisation=lbm_opt)
-stencil_typedefs = {'Stencil_T': stencil}
-field_typedefs = {'PdfField_T': pdfs,
- 'VelocityField_T': velocity_field,
- 'ScalarField_T': density_field}
+ stencil_typedefs = {'Stencil_T': stencil}
+ field_typedefs = {'PdfField_T': pdfs,
+ 'VelocityField_T': velocity_field,
+ 'ScalarField_T': density_field}
-with CodeGeneration() as ctx:
# sweeps
generate_alternating_lbm_sweep(ctx, 'GeneratedFreeSlip_Sweep', collision_rule, lbm_config=lbm_config,
lbm_optimisation=lbm_opt, field_swaps=field_swaps)
diff --git a/tests/lbm/codegen/GeneratedOutflowBC.cpp b/tests/lbm/codegen/GeneratedOutflowBC.cpp
index dba8bdbf8..453a1e7cc 100644
--- a/tests/lbm/codegen/GeneratedOutflowBC.cpp
+++ b/tests/lbm/codegen/GeneratedOutflowBC.cpp
@@ -70,10 +70,10 @@ Vector3< real_t > ShearProfile::operator()( const Cell& pos, const shared_ptr< S
{
Cell globalCell;
CellInterval domain = SbF->getDomainCellBB();
- real_t h_y = domain.yMax() - domain.yMin();
+ real_t h_y = real_c(domain.yMax()) - real_c(domain.yMin());
SbF->transformBlockLocalToGlobalCell(globalCell, block, pos);
- real_t u = inflow_velocity_ * (globalCell[1] / h_y);
+ real_t u = inflow_velocity_ * (real_c(globalCell[1]) / h_y);
Vector3< real_t > result(u, 0.0, 0.0);
return result;
@@ -93,7 +93,7 @@ int main(int argc, char** argv)
auto parameters = walberlaEnv.config()->getOneBlock("Parameters");
const real_t omega = parameters.getParameter< real_t >("omega", real_c(1.4));
- const real_t u_max = parameters.getParameter< real_t >("u_max", real_t(0.05));
+ const real_t u_max = parameters.getParameter< real_t >("u_max", real_c(0.05));
const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(10));
const double remainingTimeLoggerFrequency =
@@ -101,8 +101,8 @@ int main(int argc, char** argv)
// create fields
BlockDataID pdfFieldID = blocks->addStructuredBlockData< PdfField_T >(pdfFieldAdder, "PDFs");
- BlockDataID velFieldID = field::addToStorage< VelocityField_T >(blocks, "velocity", real_t(0), field::fzyx);
- BlockDataID densityFieldID = field::addToStorage< ScalarField_T >(blocks, "density", real_t(0), field::fzyx);
+ BlockDataID velFieldID = field::addToStorage< VelocityField_T >(blocks, "velocity", real_c(0.0), field::fzyx);
+ BlockDataID densityFieldID = field::addToStorage< ScalarField_T >(blocks, "density", real_c(0.0), field::fzyx);
BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
@@ -170,16 +170,16 @@ int main(int argc, char** argv)
timeloop.run();
CellInterval domain = blocks->getDomainCellBB();
- real_t h_y = domain.yMax() - domain.yMin();
+ real_t h_y = real_c(domain.yMax()) - real_c(domain.yMin());
for (auto& block : *blocks)
{
auto velField = block.getData<VelocityField_T>(velFieldID);
WALBERLA_FOR_ALL_CELLS_XYZ
(
velField,
- Cell globalCell;
+ Cell globalCell;
blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
- WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(velField->get(x, y, z, 0), u_max * (globalCell[1] / h_y), 0.01)
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(velField->get(x, y, z, 0), u_max * (real_c(globalCell[1]) / h_y), real_c(0.01))
)
}
diff --git a/tests/lbm/codegen/GeneratedOutflowBC.py b/tests/lbm/codegen/GeneratedOutflowBC.py
index 3fae265e1..a765eee4c 100644
--- a/tests/lbm/codegen/GeneratedOutflowBC.py
+++ b/tests/lbm/codegen/GeneratedOutflowBC.py
@@ -9,52 +9,55 @@ from lbmpy_walberla import generate_boundary, generate_lb_pack_info
import sympy as sp
-stencil = LBStencil(Stencil.D2Q9)
+with CodeGeneration() as ctx:
+ data_type = "float64" if ctx.double_accuracy else "float32"
+ stencil = LBStencil(Stencil.D2Q9)
-pdfs, pdfs_tmp = fields(f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): double[{stencil.D}D]", layout='fzyx')
-velocity_field, density_field = fields(f"velocity({stencil.D}), density(1) : double[{stencil.D}D]", layout='fzyx')
-omega = sp.Symbol("omega")
-u_max = sp.Symbol("u_max")
+ pdfs, pdfs_tmp = fields(f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): {data_type}[{stencil.D}D]",
+ layout='fzyx')
+ velocity_field, density_field = fields(f"velocity({stencil.D}), density(1) : {data_type}[{stencil.D}D]",
+ layout='fzyx')
+ omega = sp.Symbol("omega")
+ u_max = sp.Symbol("u_max")
-output = {
- 'density': density_field,
- 'velocity': velocity_field
-}
+ output = {
+ 'density': density_field,
+ 'velocity': velocity_field
+ }
-lbm_config = LBMConfig(method=Method.CUMULANT, stencil=stencil, relaxation_rate=omega,
- galilean_correction=stencil.Q == 27, field_name='pdfs', output=output)
+ lbm_config = LBMConfig(method=Method.CUMULANT, stencil=stencil, relaxation_rate=omega, compressible=True,
+ galilean_correction=stencil.Q == 27, field_name='pdfs', output=output)
-lbm_opt = LBMOptimisation(symbolic_field=pdfs, symbolic_temporary_field=pdfs_tmp,
- cse_global=False, cse_pdfs=False)
+ lbm_opt = LBMOptimisation(symbolic_field=pdfs, symbolic_temporary_field=pdfs_tmp,
+ cse_global=False, cse_pdfs=False)
-method = create_lb_method(lbm_config=lbm_config)
+ method = create_lb_method(lbm_config=lbm_config)
-# getter & setter
-setter_assignments = macroscopic_values_setter(method, velocity=velocity_field.center_vector,
- pdfs=pdfs, density=1.0)
+ # getter & setter
+ setter_assignments = macroscopic_values_setter(method, velocity=velocity_field.center_vector,
+ pdfs=pdfs, density=1.0)
-update_rule = create_lb_update_rule(lb_method=method, lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+ update_rule = create_lb_update_rule(lb_method=method, lbm_config=lbm_config, lbm_optimisation=lbm_opt)
-stencil_typedefs = {'Stencil_T': stencil}
-field_typedefs = {'PdfField_T': pdfs,
- 'VelocityField_T': velocity_field,
- 'ScalarField_T': density_field}
+ stencil_typedefs = {'Stencil_T': stencil}
+ field_typedefs = {'PdfField_T': pdfs,
+ 'VelocityField_T': velocity_field,
+ 'ScalarField_T': density_field}
-with CodeGeneration() as ctx:
# sweeps
generate_sweep(ctx, 'GeneratedOutflowBC_Sweep', update_rule, field_swaps=[(pdfs, pdfs_tmp)])
generate_sweep(ctx, 'GeneratedOutflowBC_MacroSetter', setter_assignments)
# boundaries
- ubb_dynamic = UBB(lambda *args: None, dim=stencil.D)
+ ubb_dynamic = UBB(lambda *args: None, dim=stencil.D, data_type=data_type)
ubb_data_handler = UBBAdditionalDataHandler(stencil, ubb_dynamic)
if stencil.D == 2:
- ubb_static = UBB([sp.Symbol("u_max"), 0])
+ ubb_static = UBB([sp.Symbol("u_max"), 0], data_type=data_type)
else:
- ubb_static = UBB([sp.Symbol("u_max"), 0, 0])
+ ubb_static = UBB([sp.Symbol("u_max"), 0, 0], data_type=data_type)
- outflow = ExtrapolationOutflow(stencil[4], method)
+ outflow = ExtrapolationOutflow(stencil[4], method, data_type=data_type)
outflow_data_handler = OutflowAdditionalDataHandler(stencil, outflow)
# Dynamic UBB which is used to produce a specific velocity profile at the inflow.
diff --git a/tests/lbm/codegen/InplaceStreamingCodegen.py b/tests/lbm/codegen/InplaceStreamingCodegen.py
index 099ad4617..95f30ff59 100644
--- a/tests/lbm/codegen/InplaceStreamingCodegen.py
+++ b/tests/lbm/codegen/InplaceStreamingCodegen.py
@@ -12,51 +12,53 @@ from lbmpy.advanced_streaming import Timestep
from pystencils import Field
-# Common Setup
-
-stencil = LBStencil(Stencil.D3Q27)
-target = Target.CPU
-inplace_pattern = 'aa'
-two_fields_pattern = 'pull'
-namespace = 'lbmpy'
-
-f_field = Field.create_generic('f', stencil.D, index_shape=(stencil.Q,), layout='fzyx')
-f_field_tmp = Field.create_generic('f_tmp', stencil.D, index_shape=(stencil.Q,), layout='fzyx')
-u_field = Field.create_generic('u', stencil.D, index_shape=(stencil.D,), layout='fzyx')
+with CodeGeneration() as ctx:
+ # Common Setup
+ data_type = "float64" if ctx.double_accuracy else "float32"
+ stencil = LBStencil(Stencil.D3Q27)
+ target = Target.CPU
+ inplace_pattern = 'aa'
+ two_fields_pattern = 'pull'
+ namespace = 'lbmpy'
-output = {
- 'velocity': u_field
-}
+ f_field = Field.create_generic('f', stencil.D, dtype=data_type, index_shape=(stencil.Q,), layout='fzyx')
+ f_field_tmp = Field.create_generic('f_tmp', stencil.D, dtype=data_type, index_shape=(stencil.Q,), layout='fzyx')
+ u_field = Field.create_generic('u', stencil.D, dtype=data_type, index_shape=(stencil.D,), layout='fzyx')
-lbm_config = LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=1.5, output=output)
-lbm_opt = LBMOptimisation(symbolic_field=f_field,
- symbolic_temporary_field=f_field_tmp)
+ output = {
+ 'velocity': u_field
+ }
-config = CreateKernelConfig(target=target)
+ lbm_config = LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=1.5, output=output)
+ lbm_opt = LBMOptimisation(symbolic_field=f_field,
+ symbolic_temporary_field=f_field_tmp)
-collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt, config=config)
-lb_method = collision_rule.method
-noslip = NoSlip()
-ubb = UBB((0.05,) + (0,) * (stencil.D - 1))
+ collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+ lb_method = collision_rule.method
+ noslip = NoSlip()
+ ubb = UBB((0.05,) + (0.0,) * (stencil.D - 1), data_type=data_type)
-outflow_normal = (1,) + (0,) * (stencil.D - 1)
-outflow_pull = ExtrapolationOutflow(outflow_normal, lb_method, streaming_pattern=two_fields_pattern)
+ outflow_normal = (1,) + (0,) * (stencil.D - 1)
+ outflow_pull = ExtrapolationOutflow(outflow_normal, lb_method, streaming_pattern=two_fields_pattern,
+ data_type=data_type)
-outflow_inplace = ExtrapolationOutflow(outflow_normal, lb_method, streaming_pattern=inplace_pattern)
+ outflow_inplace = ExtrapolationOutflow(outflow_normal, lb_method, streaming_pattern=inplace_pattern,
+ data_type=data_type)
-init_velocity = (0,) * stencil.D
+ init_velocity = (0.0,) * stencil.D
-init_kernel_pull = macroscopic_values_setter(lb_method, 1, init_velocity, f_field, streaming_pattern=two_fields_pattern)
-init_kernel_inplace = macroscopic_values_setter(
- lb_method, 1, init_velocity, f_field, streaming_pattern=inplace_pattern, previous_timestep=Timestep.ODD)
+ init_kernel_pull = macroscopic_values_setter(lb_method, 1.0, init_velocity, f_field,
+ streaming_pattern=two_fields_pattern)
+ init_kernel_inplace = macroscopic_values_setter(
+ lb_method, 1.0, init_velocity, f_field, streaming_pattern=inplace_pattern, previous_timestep=Timestep.ODD)
-stencil_typedefs = {'Stencil_T': stencil}
-field_typedefs = {'PdfField_T': f_field, 'VelocityField_T': u_field}
+ stencil_typedefs = {'Stencil_T': stencil}
+ field_typedefs = {'PdfField_T': f_field, 'VelocityField_T': u_field}
-with CodeGeneration() as ctx:
# Pull-Pattern classes
+ config = CreateKernelConfig(data_type=data_type, default_number_float=data_type)
ast_pull = create_lb_ast(collision_rule=collision_rule,
- streaming_pattern=two_fields_pattern, lbm_optimisation=lbm_opt)
+ streaming_pattern=two_fields_pattern, lbm_optimisation=lbm_opt, config=config)
generate_sweep(ctx, 'PullSweep', ast_pull, field_swaps=[(f_field, f_field_tmp)], namespace=namespace)
generate_boundary(ctx, 'PullNoSlip', noslip, lb_method,
diff --git a/tests/lbm/codegen/InplaceStreamingCodegenTest.cpp b/tests/lbm/codegen/InplaceStreamingCodegenTest.cpp
index 718042efa..f43f5dbd7 100644
--- a/tests/lbm/codegen/InplaceStreamingCodegenTest.cpp
+++ b/tests/lbm/codegen/InplaceStreamingCodegenTest.cpp
@@ -32,7 +32,6 @@
#include "timeloop/SweepTimeloop.h"
-#include <string>
using namespace walberla::lbmpy;
@@ -151,7 +150,7 @@ int main(int argc, char** argv)
{
for (uint_t i = 0; i < pullVelocityField->fSize(); i++)
{
- WALBERLA_CHECK_FLOAT_EQUAL(pullVelocityField->get(c, i), inplaceVelocityField->get(c, i));
+ WALBERLA_CHECK_FLOAT_EQUAL(pullVelocityField->get(c, i), inplaceVelocityField->get(c, i))
}
}
}
diff --git a/tests/lbm/codegen/LbCodeGenerationExample.cpp b/tests/lbm/codegen/LbCodeGenerationExample.cpp
index 8bb6f96f7..4711fb1b9 100644
--- a/tests/lbm/codegen/LbCodeGenerationExample.cpp
+++ b/tests/lbm/codegen/LbCodeGenerationExample.cpp
@@ -69,16 +69,16 @@ int main(int argc, char** argv)
const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(10));
const double remainingTimeLoggerFrequency =
- parameters.getParameter< double >("remainingTimeLoggerFrequency", 3.0); // in seconds
+ parameters.getParameter< double >("remainingTimeLoggerFrequency", real_c(3.0)); // in seconds
// create fields
- BlockDataID forceFieldId = field::addToStorage< VectorField_T >(blocks, "Force", real_t(0.0), field::fzyx);
- BlockDataID velFieldId = field::addToStorage< VectorField_T >(blocks, "Velocity", real_t(0.0), field::fzyx);
- BlockDataID omegaFieldId = field::addToStorage< ScalarField_T >(blocks, "Omega", real_t(0.0), field::fzyx);
+ BlockDataID forceFieldId = field::addToStorage< VectorField_T >(blocks, "Force", real_c(0.0), field::fzyx);
+ BlockDataID velFieldId = field::addToStorage< VectorField_T >(blocks, "Velocity", real_c(0.0), field::fzyx);
+ BlockDataID omegaFieldId = field::addToStorage< ScalarField_T >(blocks, "Omega", real_c(0.0), field::fzyx);
LatticeModel_T latticeModel = LatticeModel_T(forceFieldId, omegaFieldId, velFieldId, omega);
BlockDataID pdfFieldId =
- lbm::addPdfFieldToStorage(blocks, "pdf field", latticeModel, initialVelocity, real_t(1), field::fzyx);
+ lbm::addPdfFieldToStorage(blocks, "pdf field", latticeModel, initialVelocity, real_c(log10(1)), field::fzyx);
BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
// create and initialize boundary handling
diff --git a/tests/lbm/codegen/LbCodeGenerationExample.py b/tests/lbm/codegen/LbCodeGenerationExample.py
index 182a501d4..870c3bf81 100644
--- a/tests/lbm/codegen/LbCodeGenerationExample.py
+++ b/tests/lbm/codegen/LbCodeGenerationExample.py
@@ -7,11 +7,12 @@ from pystencils_walberla import CodeGeneration, generate_info_header
from lbmpy_walberla import RefinementScaling, generate_boundary, generate_lattice_model
with CodeGeneration() as ctx:
+ data_type = "float64" if ctx.double_accuracy else "float32"
omega, omega_free = sp.symbols("omega, omega_free")
- force_field, vel_field, omega_out = ps.fields("force(3), velocity(3), omega_out: [3D]", layout='fzyx')
+ force_field, vel_field, omega_out = ps.fields(f"force(3), velocity(3), omega_out: {data_type}[3D]", layout='fzyx')
stencil = LBStencil(Stencil.D3Q19)
- lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, entropic=True,
+ lbm_config = LBMConfig(stencil=stencil, method=Method.MRT, entropic=True, zero_centered=False,
compressible=True, omega_output_field=omega_out,
force=force_field.center_vector, output={'velocity': vel_field},
relaxation_rates=[omega, omega, omega_free, omega_free, omega_free, omega_free])
@@ -33,6 +34,6 @@ with CodeGeneration() as ctx:
# for vectorisation on AVX512 due to scatter and gather intrinsics
generate_lattice_model(ctx, 'LbCodeGenerationExample_LatticeModel', collision_rule,
field_layout='fzyx', refinement_scaling=scaling)
- generate_boundary(ctx, 'LbCodeGenerationExample_UBB', UBB([0.05, 0, 0]), collision_rule.method)
+ generate_boundary(ctx, 'LbCodeGenerationExample_UBB', UBB([0.05, 0, 0], data_type=data_type), collision_rule.method)
generate_boundary(ctx, 'LbCodeGenerationExample_NoSlip', NoSlip(), collision_rule.method)
generate_info_header(ctx, 'LbCodeGenerationExample')
diff --git a/tests/lbm/codegen/LbmPackInfoGenerationTest.py b/tests/lbm/codegen/LbmPackInfoGenerationTest.py
index 9e24470bb..87af2b6dc 100644
--- a/tests/lbm/codegen/LbmPackInfoGenerationTest.py
+++ b/tests/lbm/codegen/LbmPackInfoGenerationTest.py
@@ -14,7 +14,7 @@ with CodeGeneration() as ctx:
target = ps.Target.CPU
stencil = LBStencil(Stencil.D3Q19)
- pdf_field = Field.create_generic('pdfs', stencil.D, index_shape=(stencil.Q,), layout='fzyx')
+ pdf_field = Field.create_generic('pdfs', stencil.D, index_shape=(stencil.Q,), layout='fzyx', dtype=data_type)
lbm_config = LBMConfig(method=Method.SRT, stencil=stencil, streaming_pattern=streaming_pattern,
timestep=Timestep.ODD)
diff --git a/tests/lbm/codegen/StreamInCellIntervalCodegen.py b/tests/lbm/codegen/StreamInCellIntervalCodegen.py
index 4205d0c4d..7391a0d66 100644
--- a/tests/lbm/codegen/StreamInCellIntervalCodegen.py
+++ b/tests/lbm/codegen/StreamInCellIntervalCodegen.py
@@ -1,13 +1,14 @@
import sympy as sp
-from lbmpy.creationfunctions import create_lb_collision_rule, create_lb_update_rule
-from lbmpy import LBMConfig
+from lbmpy.creationfunctions import create_lb_collision_rule
+from lbmpy import LBMConfig, LBMOptimisation, Method, Stencil
+from lbmpy.stencils import LBStencil
from pystencils_walberla import CodeGeneration
from lbmpy_walberla import generate_lattice_model
# general parameters
-stencil = 'D3Q19'
+stencil = LBStencil(Stencil.D3Q19)
omega = sp.Symbol('omega')
layout = 'fzyx'
@@ -20,10 +21,12 @@ method_params = {'stencil': stencil,
'relaxation_rate': omega,
'compressible': True}
-lbm_config = LBMConfig(streaming_pattern='pull')
+lbm_config = LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=omega, compressible=True,
+ zero_centered=False, delta_equilibrium=False,
+ streaming_pattern='pull')
+lbm_opt = LBMOptimisation(cse_global=True, field_layout=layout)
-collision_rule = create_lb_collision_rule(optimization=optimizations, **method_params)
-update_rule = create_lb_update_rule(collision_rule=collision_rule, lbm_config=lbm_config, optimization=optimizations)
+collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
with CodeGeneration() as ctx:
generate_lattice_model(ctx, "GeneratedLatticeModel", collision_rule, field_layout=layout)
diff --git a/tests/lbm/codegen/StreamInCellIntervalCodegenTest.cpp b/tests/lbm/codegen/StreamInCellIntervalCodegenTest.cpp
index 0698ccc13..7981f5b2b 100644
--- a/tests/lbm/codegen/StreamInCellIntervalCodegenTest.cpp
+++ b/tests/lbm/codegen/StreamInCellIntervalCodegenTest.cpp
@@ -64,12 +64,12 @@ int main(int argc, char** argv)
const std::shared_ptr< StructuredBlockForest > blockForest =
blockforest::createUniformBlockGrid(numBlocks[0], numBlocks[1], numBlocks[2], // blocks
cellsPerBlock[0], cellsPerBlock[1], cellsPerBlock[2], // cells
- real_t(1.0), // dx
+ real_c(1.0), // dx
true, // one block per process
true, true, true); // periodicity
// relaxation rate
- real_t omega = real_t(1.8);
+ real_t omega = real_c(1.8);
// create lattice model
LatticeModel_T latticeModel = LatticeModel_T(omega);
@@ -87,8 +87,8 @@ int main(int argc, char** argv)
{
PdfField_T* pdfField = blockIt->getData< PdfField_T >(pdfFieldID);
- pdfField->setDensityAndVelocity(cell0, Vector3< real_t >(real_t(0.1), real_t(0.1), real_t(0.1)), real_t(1.1));
- pdfField->setDensityAndVelocity(cell1, Vector3< real_t >(real_t(-0.1), real_t(-0.1), real_t(-0.1)), real_t(0.9));
+ pdfField->setDensityAndVelocity(cell0, Vector3< real_t >(real_c(0.1), real_c(0.1), real_c(0.1)), real_c(1.1));
+ pdfField->setDensityAndVelocity(cell1, Vector3< real_t >(real_c(-0.1), real_c(-0.1), real_c(-0.1)), real_c(0.9));
}
// create communication for PDF fields
--
GitLab