From 05b337550bbe759432e25d49db7c81e6d2e049cc Mon Sep 17 00:00:00 2001
From: Markus Holzer <markus.holzer@fau.de>
Date: Mon, 12 Apr 2021 08:09:31 +0000
Subject: [PATCH] FIX: SP/DP conversion warnings
---
.../FluidParticleCoupling/GeneratedLBM.py | 104 +++++-----
.../GeneratedLBMWithForce.py | 191 +++++++++---------
.../InitializerFunctions.cpp | 19 +-
.../benchmark_multiphase.cpp | 2 +-
.../UniformGridGPU/InitShearVelocity.h | 2 +-
.../UniformGridGPU/UniformGridGPU.cpp | 2 +-
.../UniformGridGPU/UniformGridGPU.py | 2 +-
.../UniformGridGenerated/InitShearVelocity.h | 2 +-
.../UniformGridGenerated/ManualKernels.h | 16 +-
.../UniformGridGenerated.cpp | 10 +-
.../UniformGridGenerated.py | 20 +-
.../lbmpy_walberla/walberla_lbm_generation.py | 48 +++--
12 files changed, 221 insertions(+), 197 deletions(-)
diff --git a/apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py b/apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py
index 9772c5ba1..982aa6167 100644
--- a/apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py
+++ b/apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py
@@ -14,19 +14,21 @@ from collections import OrderedDict
with CodeGeneration() as ctx:
- generatedMethod = "TRTlike"
- #generatedMethod = "D3Q27TRTlike"
- #generatedMethod = "cumulant"
+ generatedMethod = 'TRTlike'
+ #generatedMethod = 'D3Q27TRTlike'
+ #generatedMethod = 'cumulant'
clear_cache()
+ dtype_string = 'float64' if ctx.double_accuracy else 'float32'
+
cpu_vectorize_info = {'instruction_set': get_vectorize_instruction_set(ctx)}
- if generatedMethod == "TRTlike":
- omegaVisc = sp.Symbol("omega_visc")
- omegaBulk = ps.fields("omega_bulk: [3D]", layout='fzyx')
- omegaMagic = sp.Symbol("omega_magic")
- stencil = get_stencil("D3Q19", 'walberla')
+ if generatedMethod == 'TRTlike':
+ 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')
x, y, z = MOMENT_SYMBOLS
one = sp.Rational(1, 1)
@@ -53,61 +55,62 @@ with CodeGeneration() as ctx:
generate_lattice_model(ctx, 'GeneratedLBM', collision_rule, field_layout='fzyx', cpu_vectorize_info=cpu_vectorize_info)
- if generatedMethod == "D3Q27TRTlike":
+ if generatedMethod == 'D3Q27TRTlike':
- omegaVisc = sp.Symbol("omega_visc")
- omegaBulk = ps.fields("omega_bulk: [3D]", layout='fzyx')
- omegaMagic = sp.Symbol("omega_magic")
- stencil = get_stencil("D3Q27", 'walberla')
+ 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')
- relaxation_rates=[omegaVisc, omegaBulk.center_vector, omegaMagic, omegaVisc, omegaMagic, omegaVisc]
+ 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)
+ methodWithForce = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False, weighted=True,
+ compressible=False, relaxation_rates=relaxation_rates)
collision_rule = create_lb_collision_rule(lb_method=methodWithForce, optimization={'cse_global': True})
- generate_lattice_model(ctx, 'GeneratedLBM', collision_rule, field_layout='fzyx', cpu_vectorize_info=cpu_vectorize_info)
+ generate_lattice_model(ctx, 'GeneratedLBM', collision_rule, field_layout='fzyx',
+ cpu_vectorize_info=cpu_vectorize_info)
- if generatedMethod == "cumulant":
+ if generatedMethod == 'cumulant':
- x,y,z = MOMENT_SYMBOLS
+ x, y, z = MOMENT_SYMBOLS
cumulants = [0] * 27
- cumulants[0] = sp.sympify(1) #000
+ cumulants[0] = sp.sympify(1) # 000
- cumulants[1] = x #100
- cumulants[2] = y #010
- cumulants[3] = z #001
+ cumulants[1] = x # 100
+ cumulants[2] = y # 010
+ cumulants[3] = z # 001
- cumulants[4] = x*y #110
- cumulants[5] = x*z #101
- cumulants[6] = y*z #011
+ cumulants[4] = x*y # 110
+ cumulants[5] = x*z # 101
+ cumulants[6] = y*z # 011
- cumulants[7] = x**2 - y**2 #200 - 020
- cumulants[8] = x**2 - z**2 #200 - 002
- cumulants[9] = x**2 + y**2 + z**2 #200 + 020 + 002
+ cumulants[7] = x**2 - y**2 # 200 - 020
+ cumulants[8] = x**2 - z**2 # 200 - 002
+ cumulants[9] = x**2 + y**2 + z**2 # 200 + 020 + 002
- cumulants[10] = x*y**2 + x*z**2 #120 + 102
- cumulants[11] = x**2*y + y*z**2 #210 + 012
- cumulants[12] = x**2*z + y**2*z #201 + 021
- cumulants[13] = x*y**2 - x*z**2 #120 - 102
- cumulants[14] = x**2*y - y*z**2 #210 - 012
- cumulants[15] = x**2*z - y**2*z #201 - 021
+ cumulants[10] = x*y**2 + x*z**2 # 120 + 102
+ cumulants[11] = x**2*y + y*z**2 # 210 + 012
+ cumulants[12] = x**2*z + y**2*z # 201 + 021
+ cumulants[13] = x*y**2 - x*z**2 # 120 - 102
+ cumulants[14] = x**2*y - y*z**2 # 210 - 012
+ cumulants[15] = x**2*z - y**2*z # 201 - 021
- cumulants[16] = x*y*z #111
+ cumulants[16] = x*y*z # 111
- cumulants[17] = x**2*y**2 - 2*x**2*z**2 + y**2*z**2 # 220- 2*202 +022
- cumulants[18] = x**2*y**2 + x**2*z**2 - 2*y**2*z**2 # 220 + 202 - 2*002
- cumulants[19] = x**2*y**2 + x**2*z**2 + y**2*z**2 # 220 + 202 + 022
+ cumulants[17] = x**2*y**2 - 2*x**2*z**2 + y**2*z**2 # 220- 2*202 +022
+ cumulants[18] = x**2*y**2 + x**2*z**2 - 2*y**2*z**2 # 220 + 202 - 2*002
+ cumulants[19] = x**2*y**2 + x**2*z**2 + y**2*z**2 # 220 + 202 + 022
- cumulants[20] = x**2*y*z # 211
- cumulants[21] = x*y**2*z # 121
- cumulants[22] = x*y*z**2 # 112
+ cumulants[20] = x**2*y*z # 211
+ cumulants[21] = x*y**2*z # 121
+ cumulants[22] = x*y*z**2 # 112
- cumulants[23] = x**2*y**2*z # 221
- cumulants[24] = x**2*y*z**2 # 212
- cumulants[25] = x*y**2*z**2 # 122
+ cumulants[23] = x**2*y**2*z # 221
+ cumulants[24] = x**2*y*z**2 # 212
+ cumulants[25] = x*y**2*z**2 # 122
cumulants[26] = x**2*y**2*z**2 # 222
@@ -119,9 +122,9 @@ with CodeGeneration() as ctx:
else:
return 1
- stencil = get_stencil("D3Q27", 'walberla')
+ stencil = get_stencil('D3Q27', 'walberla')
- omega = sp.Symbol("omega")
+ omega = sp.Symbol('omega')
rr_dict = OrderedDict((c, get_relaxation_rate(c, omega))
for c in cumulants)
@@ -129,12 +132,13 @@ with CodeGeneration() as ctx:
my_method = create_with_continuous_maxwellian_eq_moments(stencil, rr_dict, cumulant=True, compressible=True,)
collision_rule = create_lb_collision_rule(lb_method=my_method,
- optimization={"cse_global": True,
- "cse_pdfs": False})
+ optimization={'cse_global': True,
+ 'cse_pdfs': False})
print(my_method.relaxation_rates)
- generate_lattice_model(ctx, 'GeneratedLBM', collision_rule, field_layout='fzyx', cpu_vectorize_info=cpu_vectorize_info)
+ 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 6f7202525..9eb861b8c 100644
--- a/apps/benchmarks/FluidParticleCoupling/GeneratedLBMWithForce.py
+++ b/apps/benchmarks/FluidParticleCoupling/GeneratedLBMWithForce.py
@@ -15,25 +15,27 @@ from collections import OrderedDict
with CodeGeneration() as ctx:
- forcing=(sp.symbols("fx"),0,0)
- forcemodel=Luo(forcing)
+ forcing = (sp.symbols('fx'), 0, 0)
+ forcemodel = Luo(forcing)
- generatedMethod = "TRTlike"
- #generatedMethod = "D3Q27TRTlike"
- #generatedMethod = "cumulant"
- #generatedMethod = "cumulantTRT"
+ generatedMethod = 'TRTlike'
+ # generatedMethod = 'D3Q27TRTlike'
+ # generatedMethod = 'cumulant'
+ # generatedMethod = 'cumulantTRT'
- print("Generating " + generatedMethod + " LBM method")
+ print('Generating ' + generatedMethod + ' LBM method')
clear_cache()
+ dtype_string = 'float64' if ctx.double_accuracy else 'float32'
+
cpu_vectorize_info = {'instruction_set': get_vectorize_instruction_set(ctx)}
- if generatedMethod == "TRTlike":
- omegaVisc = sp.Symbol("omega_visc")
- omegaBulk = ps.fields("omega_bulk: [3D]", layout='fzyx')
- omegaMagic = sp.Symbol("omega_magic")
- stencil = get_stencil("D3Q19", 'walberla')
+ if generatedMethod == 'TRTlike':
+ 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')
x, y, z = MOMENT_SYMBOLS
one = sp.Rational(1, 1)
@@ -59,63 +61,64 @@ with CodeGeneration() as ctx:
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":
+ if generatedMethod == 'D3Q27TRTlike':
- omegaVisc = sp.Symbol("omega_visc")
- omegaBulk = ps.fields("omega_bulk: [3D]", layout='fzyx')
- omegaMagic = sp.Symbol("omega_magic")
- stencil = get_stencil("D3Q27", 'walberla')
+ 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')
- relaxation_rates=[omegaVisc, omegaBulk.center_vector, omegaMagic, omegaVisc, omegaMagic, omegaVisc]
+ 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)
+ methodWithForce = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False, weighted=True,
+ compressible=False, force_model=forcemodel, relaxation_rates=relaxation_rates)
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)
+ generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule, field_layout='fzyx',
+ cpu_vectorize_info=cpu_vectorize_info)
- if generatedMethod == "cumulant":
+ if generatedMethod == 'cumulant':
- x,y,z = MOMENT_SYMBOLS
+ x, y, z = MOMENT_SYMBOLS
cumulants = [0] * 27
- cumulants[0] = sp.sympify(1) #000
+ cumulants[0] = sp.sympify(1) # 000
- cumulants[1] = x #100
- cumulants[2] = y #010
- cumulants[3] = z #001
+ cumulants[1] = x # 100
+ cumulants[2] = y # 010
+ cumulants[3] = z # 001
- cumulants[4] = x*y #110
- cumulants[5] = x*z #101
- cumulants[6] = y*z #011
+ cumulants[4] = x*y # 110
+ cumulants[5] = x*z # 101
+ cumulants[6] = y*z # 011
- cumulants[7] = x**2 - y**2 #200 - 020
- cumulants[8] = x**2 - z**2 #200 - 002
- cumulants[9] = x**2 + y**2 + z**2 #200 + 020 + 002
+ cumulants[7] = x**2 - y**2 # 200 - 020
+ cumulants[8] = x**2 - z**2 # 200 - 002
+ cumulants[9] = x**2 + y**2 + z**2 # 200 + 020 + 002
- cumulants[10] = x*y**2 + x*z**2 #120 + 102
- cumulants[11] = x**2*y + y*z**2 #210 + 012
- cumulants[12] = x**2*z + y**2*z #201 + 021
- cumulants[13] = x*y**2 - x*z**2 #120 - 102
- cumulants[14] = x**2*y - y*z**2 #210 - 012
- cumulants[15] = x**2*z - y**2*z #201 - 021
+ cumulants[10] = x*y**2 + x*z**2 # 120 + 102
+ cumulants[11] = x**2*y + y*z**2 # 210 + 012
+ cumulants[12] = x**2*z + y**2*z # 201 + 021
+ cumulants[13] = x*y**2 - x*z**2 # 120 - 102
+ cumulants[14] = x**2*y - y*z**2 # 210 - 012
+ cumulants[15] = x**2*z - y**2*z # 201 - 021
- cumulants[16] = x*y*z #111
+ cumulants[16] = x*y*z # 111
- cumulants[17] = x**2*y**2 - 2*x**2*z**2 + y**2*z**2 # 220- 2*202 +022
- cumulants[18] = x**2*y**2 + x**2*z**2 - 2*y**2*z**2 # 220 + 202 - 2*002
- cumulants[19] = x**2*y**2 + x**2*z**2 + y**2*z**2 # 220 + 202 + 022
+ cumulants[17] = x**2*y**2 - 2*x**2*z**2 + y**2*z**2 # 220- 2*202 +022
+ cumulants[18] = x**2*y**2 + x**2*z**2 - 2*y**2*z**2 # 220 + 202 - 2*002
+ cumulants[19] = x**2*y**2 + x**2*z**2 + y**2*z**2 # 220 + 202 + 022
- cumulants[20] = x**2*y*z # 211
- cumulants[21] = x*y**2*z # 121
- cumulants[22] = x*y*z**2 # 112
+ cumulants[20] = x**2*y*z # 211
+ cumulants[21] = x*y**2*z # 121
+ cumulants[22] = x*y*z**2 # 112
- cumulants[23] = x**2*y**2*z # 221
- cumulants[24] = x**2*y*z**2 # 212
- cumulants[25] = x*y**2*z**2 # 122
+ cumulants[23] = x**2*y**2*z # 221
+ cumulants[24] = x**2*y*z**2 # 212
+ cumulants[25] = x*y**2*z**2 # 122
- cumulants[26] = x**2*y**2*z**2 # 222
+ cumulants[26] = x**2*y**2*z**2 # 222
def get_relaxation_rate(cumulant, omega):
if get_order(cumulant) <= 1:
@@ -125,66 +128,68 @@ with CodeGeneration() as ctx:
else:
return 1
- stencil = get_stencil("D3Q27", 'walberla')
+ stencil = get_stencil('D3Q27', 'walberla')
- omega = sp.Symbol("omega")
+ omega = sp.Symbol('omega')
rr_dict = OrderedDict((c, get_relaxation_rate(c, omega))
for c in cumulants)
from lbmpy.methods import create_with_continuous_maxwellian_eq_moments
- my_method = create_with_continuous_maxwellian_eq_moments(stencil, rr_dict, cumulant=True, compressible=True, force_model=forcemodel)
+ my_method = create_with_continuous_maxwellian_eq_moments(stencil, rr_dict, cumulant=True, compressible=True,
+ force_model=forcemodel)
collision_rule = create_lb_collision_rule(lb_method=my_method,
- optimization={"cse_global": True,
- "cse_pdfs": False})
+ optimization={'cse_global': True,
+ 'cse_pdfs': False})
print(my_method.relaxation_rates)
- generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule, field_layout='fzyx', cpu_vectorize_info=cpu_vectorize_info)
+ generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule, field_layout='fzyx',
+ cpu_vectorize_info=cpu_vectorize_info)
- if generatedMethod == "cumulantTRT":
+ if generatedMethod == 'cumulantTRT':
- x,y,z = MOMENT_SYMBOLS
+ x, y, z = MOMENT_SYMBOLS
cumulants = [0] * 27
- cumulants[0] = sp.sympify(1) #000
+ cumulants[0] = sp.sympify(1) # 000
- cumulants[1] = x #100
- cumulants[2] = y #010
- cumulants[3] = z #001
+ cumulants[1] = x # 100
+ cumulants[2] = y # 010
+ cumulants[3] = z # 001
- cumulants[4] = x*y #110
- cumulants[5] = x*z #101
- cumulants[6] = y*z #011
+ cumulants[4] = x*y # 110
+ cumulants[5] = x*z # 101
+ cumulants[6] = y*z # 011
- cumulants[7] = x**2 - y**2 #200 - 020
- cumulants[8] = x**2 - z**2 #200 - 002
- cumulants[9] = x**2 + y**2 + z**2 #200 + 020 + 002
+ cumulants[7] = x**2 - y**2 # 200 - 020
+ cumulants[8] = x**2 - z**2 # 200 - 002
+ cumulants[9] = x**2 + y**2 + z**2 # 200 + 020 + 002
- cumulants[10] = x*y**2 + x*z**2 #120 + 102
- cumulants[11] = x**2*y + y*z**2 #210 + 012
- cumulants[12] = x**2*z + y**2*z #201 + 021
- cumulants[13] = x*y**2 - x*z**2 #120 - 102
- cumulants[14] = x**2*y - y*z**2 #210 - 012
- cumulants[15] = x**2*z - y**2*z #201 - 021
+ cumulants[10] = x*y**2 + x*z**2 # 120 + 102
+ cumulants[11] = x**2*y + y*z**2 # 210 + 012
+ cumulants[12] = x**2*z + y**2*z # 201 + 021
+ cumulants[13] = x*y**2 - x*z**2 # 120 - 102
+ cumulants[14] = x**2*y - y*z**2 # 210 - 012
+ cumulants[15] = x**2*z - y**2*z # 201 - 021
- cumulants[16] = x*y*z #111
+ cumulants[16] = x*y*z # 111
- cumulants[17] = x**2*y**2 - 2*x**2*z**2 + y**2*z**2 # 220- 2*202 +022
- cumulants[18] = x**2*y**2 + x**2*z**2 - 2*y**2*z**2 # 220 + 202 - 2*002
- cumulants[19] = x**2*y**2 + x**2*z**2 + y**2*z**2 # 220 + 202 + 022
+ cumulants[17] = x**2*y**2 - 2*x**2*z**2 + y**2*z**2 # 220- 2*202 +022
+ cumulants[18] = x**2*y**2 + x**2*z**2 - 2*y**2*z**2 # 220 + 202 - 2*002
+ cumulants[19] = x**2*y**2 + x**2*z**2 + y**2*z**2 # 220 + 202 + 022
- cumulants[20] = x**2*y*z # 211
- cumulants[21] = x*y**2*z # 121
- cumulants[22] = x*y*z**2 # 112
+ cumulants[20] = x**2*y*z # 211
+ cumulants[21] = x*y**2*z # 121
+ cumulants[22] = x*y*z**2 # 112
- cumulants[23] = x**2*y**2*z # 221
- cumulants[24] = x**2*y*z**2 # 212
- cumulants[25] = x*y**2*z**2 # 122
+ cumulants[23] = x**2*y**2*z # 221
+ cumulants[24] = x**2*y*z**2 # 212
+ cumulants[25] = x*y**2*z**2 # 122
- cumulants[26] = x**2*y**2*z**2 # 222
+ cumulants[26] = x**2*y**2*z**2 # 222
def get_relaxation_rate(cumulant, omegaVisc, omegaMagic):
if get_order(cumulant) <= 1:
@@ -194,22 +199,24 @@ with CodeGeneration() as ctx:
else:
return omegaMagic
- stencil = get_stencil("D3Q27", 'walberla')
+ stencil = get_stencil('D3Q27', 'walberla')
- omegaVisc = sp.Symbol("omega_visc")
- omegaMagic = sp.Symbol("omega_magic")
+ omegaVisc = sp.Symbol('omega_visc')
+ omegaMagic = sp.Symbol('omega_magic')
rr_dict = OrderedDict((c, get_relaxation_rate(c, omegaVisc, omegaMagic))
for c in cumulants)
from lbmpy.methods import create_with_continuous_maxwellian_eq_moments
- my_method = create_with_continuous_maxwellian_eq_moments(stencil, rr_dict, cumulant=True, compressible=True, force_model=forcemodel)
+ my_method = create_with_continuous_maxwellian_eq_moments(stencil, rr_dict, cumulant=True, compressible=True,
+ force_model=forcemodel)
collision_rule = create_lb_collision_rule(lb_method=my_method,
- optimization={"cse_global": True,
- "cse_pdfs": False})
+ optimization={'cse_global': True,
+ 'cse_pdfs': False})
print(my_method.relaxation_rates)
- generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule, field_layout='fzyx', cpu_vectorize_info=cpu_vectorize_info)
+ generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule, field_layout='fzyx',
+ cpu_vectorize_info=cpu_vectorize_info)
diff --git a/apps/benchmarks/PhaseFieldAllenCahn/InitializerFunctions.cpp b/apps/benchmarks/PhaseFieldAllenCahn/InitializerFunctions.cpp
index 5f7e5b2f6..ae8c2e2e7 100644
--- a/apps/benchmarks/PhaseFieldAllenCahn/InitializerFunctions.cpp
+++ b/apps/benchmarks/PhaseFieldAllenCahn/InitializerFunctions.cpp
@@ -40,10 +40,10 @@ void initPhaseField_bubble(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]));
- phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(2.0 * (Ri - R) / W);
+ real_t Ri = std::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]));
+ phaseField->get(x, y, z) = real_t(0.5) + real_t(0.5) * std::tanh(real_t(2.0) * (Ri - R) / W);
)
// clang-format on
}
@@ -52,9 +52,9 @@ void initPhaseField_bubble(const shared_ptr< StructuredBlockStorage >& blocks, B
void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
const real_t W = 5)
{
- auto X = blocks->getDomainCellBB().xMax();
- auto halfY = (blocks->getDomainCellBB().yMax()) / 2.0;
- double perturbation = 0.05;
+ const real_t X = real_c(blocks->getDomainCellBB().xMax());
+ const real_t halfY = real_c(blocks->getDomainCellBB().yMax()) / real_t(2.0);
+ const real_t perturbation = real_t(0.05);
for (auto& block : *blocks)
{
@@ -62,8 +62,9 @@ void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, Bloc
// clang-format off
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_t tmp = perturbation * X * (std::cos((real_t(2.0) * math::pi * real_c(globalCell[0])) / X)
+ + std::cos((real_t(2.0) * math::pi * real_c(globalCell[2])) / X));
+ phaseField->get(x, y, z) = real_t(0.5) + real_t(0.5) * std::tanh(((real_t(globalCell[1]) - halfY) - tmp) / (W / real_t(2.0)));
)
// clang-format on
}
diff --git a/apps/benchmarks/PhaseFieldAllenCahn/benchmark_multiphase.cpp b/apps/benchmarks/PhaseFieldAllenCahn/benchmark_multiphase.cpp
index 8efb4756d..bbdf8c144 100644
--- a/apps/benchmarks/PhaseFieldAllenCahn/benchmark_multiphase.cpp
+++ b/apps/benchmarks/PhaseFieldAllenCahn/benchmark_multiphase.cpp
@@ -84,7 +84,7 @@ using flag_t = walberla::uint8_t;
using FlagField_T = FlagField< flag_t >;
#if defined(WALBERLA_BUILD_WITH_CUDA)
-typedef cuda::GPUField< double > GPUField;
+typedef cuda::GPUField< real_t > GPUField;
#endif
// using CommScheme_T = cuda::communication::UniformGPUScheme<stencil::D2Q9>;
diff --git a/apps/benchmarks/UniformGridGPU/InitShearVelocity.h b/apps/benchmarks/UniformGridGPU/InitShearVelocity.h
index 2aed66b1a..9a6c7d1db 100644
--- a/apps/benchmarks/UniformGridGPU/InitShearVelocity.h
+++ b/apps/benchmarks/UniformGridGPU/InitShearVelocity.h
@@ -6,7 +6,7 @@ namespace walberla {
inline void initShearVelocity(const shared_ptr<StructuredBlockStorage> & blocks, BlockDataID velFieldID,
- const real_t xMagnitude=0.005, const real_t fluctuationMagnitude=0.05 )
+ const real_t xMagnitude=real_t(0.005), const real_t fluctuationMagnitude=real_t(0.05) )
{
math::seedRandomGenerator(0);
auto halfZ = blocks->getDomainCellBB().zMax() / 2;
diff --git a/apps/benchmarks/UniformGridGPU/UniformGridGPU.cpp b/apps/benchmarks/UniformGridGPU/UniformGridGPU.cpp
index deaea0815..a123cbf34 100644
--- a/apps/benchmarks/UniformGridGPU/UniformGridGPU.cpp
+++ b/apps/benchmarks/UniformGridGPU/UniformGridGPU.cpp
@@ -55,7 +55,7 @@ using FlagField_T = FlagField<flag_t>;
void initShearVelocity(const shared_ptr<StructuredBlockStorage> & blocks, BlockDataID velFieldID,
- const real_t xMagnitude=0.1, const real_t fluctuationMagnitude=0.05 )
+ const real_t xMagnitude=real_t(0.1), const real_t fluctuationMagnitude=real_t(0.05) )
{
math::seedRandomGenerator(0);
auto halfZ = blocks->getDomainCellBB().zMax() / 2;
diff --git a/apps/benchmarks/UniformGridGPU/UniformGridGPU.py b/apps/benchmarks/UniformGridGPU/UniformGridGPU.py
index 20ea3f0f7..6f9670d1d 100644
--- a/apps/benchmarks/UniformGridGPU/UniformGridGPU.py
+++ b/apps/benchmarks/UniformGridGPU/UniformGridGPU.py
@@ -135,7 +135,7 @@ with CodeGeneration() as ctx:
]
lb_method = create_lb_method(**options)
update_rule = create_lb_update_rule(lb_method=lb_method, **options)
-
+
if not noopt:
update_rule = insert_fast_divisions(update_rule)
update_rule = insert_fast_sqrts(update_rule)
diff --git a/apps/benchmarks/UniformGridGenerated/InitShearVelocity.h b/apps/benchmarks/UniformGridGenerated/InitShearVelocity.h
index 2aed66b1a..9a6c7d1db 100644
--- a/apps/benchmarks/UniformGridGenerated/InitShearVelocity.h
+++ b/apps/benchmarks/UniformGridGenerated/InitShearVelocity.h
@@ -6,7 +6,7 @@ namespace walberla {
inline void initShearVelocity(const shared_ptr<StructuredBlockStorage> & blocks, BlockDataID velFieldID,
- const real_t xMagnitude=0.005, const real_t fluctuationMagnitude=0.05 )
+ const real_t xMagnitude=real_t(0.005), const real_t fluctuationMagnitude=real_t(0.05) )
{
math::seedRandomGenerator(0);
auto halfZ = blocks->getDomainCellBB().zMax() / 2;
diff --git a/apps/benchmarks/UniformGridGenerated/ManualKernels.h b/apps/benchmarks/UniformGridGenerated/ManualKernels.h
index efc7b59d7..2e7dc7ee6 100644
--- a/apps/benchmarks/UniformGridGenerated/ManualKernels.h
+++ b/apps/benchmarks/UniformGridGenerated/ManualKernels.h
@@ -30,20 +30,20 @@ public:
for( auto d = LM::Stencil::begin(); d != LM::Stencil::end(); ++d )
{
const real_t pdf = src->get(x - d.cx(), y - d.cy(), z - d.cz(), d.toIdx());
- u[0] += pdf * d.cx();
- u[1] += pdf * d.cy();
- u[2] += pdf * d.cz();
+ u[0] += pdf * real_c(d.cx());
+ u[1] += pdf * real_c(d.cy());
+ u[2] += pdf * real_c(d.cz());
rho += pdf;
}
// collide
- const real_t vel_sqr = u.sqrLength() * 1.5;
+ const real_t vel_sqr = u.sqrLength() * real_t(1.5);
for( auto d = LM::Stencil::begin(); d != LM::Stencil::end(); ++d )
{
const real_t pdf = src->get(x - d.cx(), y - d.cy(), z - d.cz(), d.toIdx());
- const real_t vel = d.cx()*u[0] + d.cy()*u[1] + d.cz()*u[2];
- dst->get(x, y, z, d.toIdx()) = ( 1.0 - omega_ ) * pdf +
- omega_ * LM::w[ d.toIdx() ] * ( rho - vel_sqr + 3.0*vel + 4.5*vel*vel );
+ const real_t vel = real_c(d.cx())*u[0] + real_c(d.cy())*u[1] + real_c(d.cz())*u[2];
+ dst->get(x, y, z, d.toIdx()) = ( real_t(1.0) - omega_ ) * pdf +
+ omega_ * LM::w[ d.toIdx() ] * ( rho - vel_sqr + real_t(3.0)*vel + real_t(4.5)*vel*vel );
}
});
src->swapDataPointers(*dst);
@@ -166,7 +166,7 @@ public:
dst->get(x,y,z,Stencil::idx[TN]) = omega_trm_ * dd_tmp_TN + omega_w2_ * ( vel_trm_TN_BS + velYpZ );
dst->get(x,y,z,Stencil::idx[BS]) = omega_trm_ * dd_tmp_BS + omega_w2_ * ( vel_trm_TN_BS - velYpZ );
- });
+ })
src->swapDataPointers(*dst);
}
diff --git a/apps/benchmarks/UniformGridGenerated/UniformGridGenerated.cpp b/apps/benchmarks/UniformGridGenerated/UniformGridGenerated.cpp
index b4275d395..844ae03d6 100644
--- a/apps/benchmarks/UniformGridGenerated/UniformGridGenerated.cpp
+++ b/apps/benchmarks/UniformGridGenerated/UniformGridGenerated.cpp
@@ -49,7 +49,7 @@ int main( int argc, char **argv )
for( auto cfg = python_coupling::configBegin( argc, argv ); cfg != python_coupling::configEnd(); ++cfg )
{
- WALBERLA_MPI_WORLD_BARRIER();
+ WALBERLA_MPI_WORLD_BARRIER()
auto config = *cfg;
logging::configureLogging( config );
@@ -61,7 +61,7 @@ int main( int argc, char **argv )
const std::string timeStepMode = parameters.getParameter<std::string>( "timeStepMode", "twoField");
const real_t omega = parameters.getParameter<real_t>( "omega", real_c( 1.4 ));
uint_t timesteps = parameters.getParameter<uint_t>( "timesteps", uint_c( 60 ));
- const real_t shearVelocityMagnitude = parameters.getParameter<real_t>("shearVelocityMagnitude", 0.08);
+ const real_t shearVelocityMagnitude = parameters.getParameter<real_t>("shearVelocityMagnitude", real_t(0.08));
const bool directComm = parameters.getParameter<bool>("directComm", false);
auto pdfFieldAdder = [](IBlock* const block, StructuredBlockStorage * const storage) {
@@ -120,7 +120,7 @@ int main( int argc, char **argv )
twoFieldKernel = StreamPullCollideD3Q19<MyLM>(pdfFieldId, tmpPdfFieldId, omega);
} else {
WALBERLA_ABORT_NO_DEBUG_INFO("Invalid option for \"twoFieldKernelType\", "
- "valid options are \"generated\", \"manualGeneric\", \"manualD3Q19\"");
+ "valid options are \"generated\", \"manualGeneric\", \"manualD3Q19\"")
}
using F = std::function<void()>;
@@ -144,7 +144,7 @@ int main( int argc, char **argv )
timeLoop.add() << Sweep( pystencils::GenLbKernelAAEven(pdfFieldId, omega), "AA Even" );
timeLoop.add() << Sweep( pystencils::GenLbKernelAAOdd(pdfFieldId, omega), "AA Odd");
} else {
- WALBERLA_ABORT("Invalid value for timeStepMode");
+ WALBERLA_ABORT("Invalid value for timeStepMode")
}
@@ -205,7 +205,7 @@ int main( int argc, char **argv )
" time steps: " << timesteps <<
setw(15) << " block size: " << cellsPerBlock <<
" mlups/core: " << int(mlupsPerProcess/ threads) <<
- " mlups: " << int(mlupsPerProcess) * MPIManager::instance()->numProcesses());
+ " mlups: " << int(mlupsPerProcess) * MPIManager::instance()->numProcesses())
WALBERLA_ROOT_SECTION()
{
diff --git a/apps/benchmarks/UniformGridGenerated/UniformGridGenerated.py b/apps/benchmarks/UniformGridGenerated/UniformGridGenerated.py
index d42452222..3b489ab1a 100644
--- a/apps/benchmarks/UniformGridGenerated/UniformGridGenerated.py
+++ b/apps/benchmarks/UniformGridGenerated/UniformGridGenerated.py
@@ -6,9 +6,9 @@ from pystencils_walberla import CodeGeneration, generate_pack_info_from_kernel,
from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
from lbmpy.fieldaccess import AAEvenTimeStepAccessor, AAOddTimeStepAccessor
-omega = sp.symbols("omega")
-omega_free = sp.Symbol("omega_free")
-omega_fill = sp.symbols("omega_:10")
+omega = sp.symbols('omega')
+omega_free = sp.Symbol('omega_free')
+omega_fill = sp.symbols('omega_:10')
options_dict = {
'srt': {
@@ -95,9 +95,9 @@ with CodeGeneration() as ctx:
config_name = ctx.config
noopt = False
d3q27 = False
- if config_name.endswith("_d3q27"):
+ if config_name.endswith('_d3q27'):
d3q27 = True
- config_name = config_name[:-len("_d3q27")]
+ config_name = config_name[:-len('_d3q27')]
if config_name == '':
config_name = 'trt'
@@ -108,9 +108,11 @@ with CodeGeneration() as ctx:
if d3q27:
options['stencil'] = 'D3Q27'
+ dtype_string = 'float64' if ctx.double_accuracy else 'float32'
+
stencil_str = options['stencil']
q = int(stencil_str[stencil_str.find('Q') + 1:])
- pdfs, velocity_field = ps.fields("pdfs({q}), velocity(3) : double[3D]".format(q=q), layout='fzyx')
+ pdfs, velocity_field = ps.fields(f'pdfs({q}), velocity(3) : {dtype_string}[3D]', layout='fzyx')
update_rule_two_field = create_lb_update_rule(optimization={'symbolic_field': pdfs,
'split': opts['two_field_split'],
@@ -129,10 +131,10 @@ with CodeGeneration() as ctx:
((0, 0, 1), UBB([0.05, 0, 0])),
((0, 0, -1), NoSlip()),
))
- cr_even = create_lb_collision_rule(stencil="D3Q19", compressible=False,
+ cr_even = create_lb_collision_rule(stencil='D3Q19', compressible=False,
optimization={'cse_global': opts['aa_even_cse_global'],
'cse_pdfs': opts['aa_even_cse_pdfs']})
- cr_odd = create_lb_collision_rule(stencil="D3Q19", compressible=False,
+ cr_odd = create_lb_collision_rule(stencil='D3Q19', compressible=False,
optimization={'cse_global': opts['aa_odd_cse_global'],
'cse_pdfs': opts['aa_odd_cse_pdfs']})
update_rule_aa_even = update_rule_with_push_boundaries(cr_even, pdfs, boundaries,
@@ -196,4 +198,4 @@ with CodeGeneration() as ctx:
'configName': ctx.config,
'optimizationDict': str(opts),
}
- ctx.write_file("GenDefines.h", info_header.format(**infoHeaderParams))
+ ctx.write_file('GenDefines.h', info_header.format(**infoHeaderParams))
diff --git a/python/lbmpy_walberla/walberla_lbm_generation.py b/python/lbmpy_walberla/walberla_lbm_generation.py
index 82867d8b8..698843f3b 100644
--- a/python/lbmpy_walberla/walberla_lbm_generation.py
+++ b/python/lbmpy_walberla/walberla_lbm_generation.py
@@ -13,7 +13,7 @@ from lbmpy.updatekernels import create_lbm_kernel, create_stream_pull_only_kerne
from pystencils import AssignmentCollection, create_kernel
from pystencils.astnodes import SympyAssignment
from pystencils.backends.cbackend import CBackend, CustomSympyPrinter, get_headers
-from pystencils.data_types import TypedSymbol, type_all_numbers
+from pystencils.data_types import TypedSymbol, type_all_numbers, cast_func
from pystencils.field import Field
from pystencils.stencil import have_same_entries, offset_to_direction_string
from pystencils.sympyextensions import get_symmetric_part
@@ -31,33 +31,36 @@ def __lattice_model(generation_context, class_name, lb_method, stream_collide_as
if not stencil_name:
raise ValueError("lb_method uses a stencil that is not supported in waLBerla")
is_float = not generation_context.double_accuracy
- dtype = np.float32 if is_float else np.float64
+ dtype_string = "float32" if is_float else "float64"
vel_symbols = lb_method.conserved_quantity_computation.first_order_moment_symbols
- rho_sym = sp.Symbol("rho")
- pdfs_sym = sp.symbols("f_:%d" % (len(lb_method.stencil),))
+ rho_sym = sp.Symbol('rho')
+ pdfs_sym = sp.symbols(f'f_:{len(lb_method.stencil)}')
vel_arr_symbols = [IndexedBase(sp.Symbol('u'), shape=(1,))[i] for i in range(len(vel_symbols))]
- momentum_density_symbols = [sp.Symbol("md_%d" % (i,)) for i in range(len(vel_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, "float32" if is_float else "float64", False)
+ _, _, equilibrium = add_types(equilibrium.main_assignments, dtype_string, False)
equilibrium = sp.Matrix([e.rhs for e in equilibrium])
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)
force_model = lb_method.force_model
macroscopic_velocity_shift = None
if force_model:
if hasattr(force_model, 'macroscopic_velocity_shift'):
- macroscopic_velocity_shift = [expression_to_code(e, "lm.", ["rho"], dtype=dtype)
+ 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(rho_sym)]
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,
+ 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,
variables_without_prefix=['rho_in', 'u'])
momentum_density_getter = cqc.output_equations_from_pdfs(pdfs_sym, {'density': rho_sym,
'momentum_density': momentum_density_symbols})
@@ -66,7 +69,7 @@ def __lattice_model(generation_context, class_name, lb_method, stream_collide_as
required_headers = get_headers(stream_collide_ast)
if refinement_scaling:
- refinement_scaling_info = [(e0, e1, expression_to_code(e2, '', dtype=dtype)) for e0, e1, e2 in
+ refinement_scaling_info = [(e0, e1, expression_to_code(e2, '', dtype=dtype_string)) 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)):
@@ -96,9 +99,9 @@ def __lattice_model(generation_context, class_name, lb_method, stream_collide_as
'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),
+ variables_without_prefix=[e.name for e in pdfs_sym], dtype=dtype_string),
'momentum_density_getter': equations_to_code(momentum_density_getter, variables_without_prefix=pdfs_sym,
- dtype=dtype),
+ dtype=dtype_string),
'density_velocity_setter_macroscopic_values': density_velocity_setter_macroscopic_values,
'refinement_scaling_info': refinement_scaling_info,
@@ -120,8 +123,8 @@ def __lattice_model(generation_context, class_name, lb_method, stream_collide_as
header = env.get_template('LatticeModel.tmpl.h').render(**jinja_context)
source = env.get_template('LatticeModel.tmpl.cpp').render(**jinja_context)
- generation_context.write_file("{}.h".format(class_name), header)
- generation_context.write_file("{}.cpp".format(class_name), source)
+ generation_context.write_file(f"{class_name}.h", header)
+ generation_context.write_file(f"{class_name}.cpp", source)
def generate_lattice_model(generation_context, class_name, collision_rule, field_layout='zyxf', refinement_scaling=None,
@@ -232,7 +235,9 @@ def stencil_switch_statement(stencil, values):
return template.render(dir_to_value_dict=dir_to_value_dict, undefined=StrictUndefined)
-def field_and_symbol_substitute(expr, variable_prefix="lm.", variables_without_prefix=[]):
+def field_and_symbol_substitute(expr, variable_prefix="lm.", variables_without_prefix=None):
+ if variables_without_prefix is None:
+ variables_without_prefix = []
variables_without_prefix = [v.name if isinstance(v, sp.Symbol) else v for v in variables_without_prefix]
substitutions = {}
# check for member access
@@ -249,26 +254,29 @@ def field_and_symbol_substitute(expr, variable_prefix="lm.", variables_without_p
else:
postfix2 = ''
if fa.field.index_dimensions == 0:
- substitutions[fa] = sp.Symbol("%s%s->get(x,y,z)" % (prefix, fa.field.name + postfix2))
+ substitutions[fa] = sp.Symbol(f"{prefix}{fa.field.name + postfix2}->get(x,y,z)")
else:
assert fa.field.index_dimensions == 1, "walberla supports only 0 or 1 index dimensions"
- substitutions[fa] = sp.Symbol("%s%s->get(x,y,z,%s)" % (prefix, fa.field.name + postfix2, fa.index[0]))
+ substitutions[fa] = sp.Symbol(f"{prefix}{fa.field.name + postfix2}->get(x,y,z,{fa.index[0]})")
else:
if sym.name not in variables_without_prefix:
substitutions[sym] = sp.Symbol(variable_prefix + sym.name + postfix)
return expr.subs(substitutions)
-def expression_to_code(expr, variable_prefix="lm.", variables_without_prefix=[], dtype="double"):
+def expression_to_code(expr, variable_prefix="lm.", variables_without_prefix=None, dtype="double"):
"""
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)
+ :param dtype: default data type used for numbers in the code
:param expr: sympy expression
:param variable_prefix: all variables (and field) are prefixed with this string
this is used for member variables mostly
:param variables_without_prefix: this variables are not prefixed
:return: code string
"""
+ if variables_without_prefix is None:
+ variables_without_prefix = []
return cpp_printer.doprint(
type_expr(field_and_symbol_substitute(expr, variable_prefix, variables_without_prefix), dtype=dtype))
@@ -285,7 +293,9 @@ def type_expr(eq, 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=[], dtype="double"):
+def equations_to_code(equations, variable_prefix="lm.", variables_without_prefix=None, dtype="double"):
+ if variables_without_prefix is None:
+ variables_without_prefix = []
if isinstance(equations, AssignmentCollection):
equations = equations.all_assignments
--
GitLab