Merge branch 'master' into master-p

.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-*

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)

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,

apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py

 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)
-
-
-

apps/benchmarks/FluidParticleCoupling/GeneratedLBMWithForce.py

 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)
-

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)

apps/benchmarks/PhaseFieldAllenCahn/multiphase_codegen.py

-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
+    # 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,
+