Adapted documentation

7c1a5ae9 · Markus Holzer · d38d207b · 7c1a5ae9 · 7c1a5ae9 · 7c1a5ae9
Commit 7c1a5ae9 authored Sep 20, 2021 by Markus Holzer
--- a/doc/notebooks/00_tutorial_lbmpy_walberla_overview.ipynb
+++ b/doc/notebooks/00_tutorial_lbmpy_walberla_overview.ipynb
--- a/doc/notebooks/01_tutorial_predefinedScenarios.ipynb
+++ b/doc/notebooks/01_tutorial_predefinedScenarios.ipynb
--- a/doc/notebooks/02_tutorial_boundary_setup.ipynb
+++ b/doc/notebooks/02_tutorial_boundary_setup.ipynb
@@ -24,12 +24,17 @@
 ```

 %% Cell type:code id: tags:

 ``` python
-sc1 = LatticeBoltzmannStep(domain_size=domain_size, method='srt', relaxation_rate=1.9, force=(1e-6, 0, 0),
-                           periodicity=(True, False, False), optimization={'target': Target.CPU})
+lbm_config = LBMConfig(stencil=Stencil.D3Q19, method=Method.SRT, relaxation_rate=1.9, force=(1e-6, 0, 0))
+config = CreateKernelConfig(target=Target.CPU)
+
+sc1 = LatticeBoltzmannStep(domain_size=domain_size,
+                           periodicity=(True, False, False),
+                           lbm_config=lbm_config,
+                           config=config)

 ```

 %% Cell type:markdown id: tags:

@@ -95,14 +100,16 @@
 Again, we set up a scenario, but this time without external force and without periodicity.

 %% Cell type:code id: tags:

 ``` python
-sc2 = LatticeBoltzmannStep(domain_size=domain_size, method='srt', relaxation_rate=1.9,
-                           optimization={'target': Target.CPU},
-                           # optimization={'target': Target.GPU, 'gpu_indexing_params': {'block_size': (16, 8, 2)}}
-                           )
+config = CreateKernelConfig(target=Target.CPU)
+# config = CreateKernelConfig(target=Target.GPU, gpu_indexing_params={'block_size': (128, 1, 1)})
+
+sc2 = LatticeBoltzmannStep(domain_size=domain_size,
+                           lbm_config=LBMConfig(method=Method.SRT, relaxation_rate=1.9),
+                           config=config)

 ```

 %% Cell type:markdown id: tags:


--- a/doc/notebooks/03_tutorial_lbm_formulation.ipynb
+++ b/doc/notebooks/03_tutorial_lbm_formulation.ipynb
@@ -60,18 +60,20 @@

 %% Cell type:code id: tags:

 ``` python
 from lbmpy.creationfunctions import create_lb_method
-method = create_lb_method(stencil=Stencil.D2Q9, method='mrt_raw')
+lbm_config = LBMConfig(stencil=Stencil.D2Q9, method=Method.MRT_RAW)
+
+method = create_lb_method(lbm_config=lbm_config)
 # check also method='srt', 'trt', 'mrt'
 method
 ```

 %%%% Output: execute_result

-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7ff41c19b850>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f0d442c0fa0>

 %% Cell type:markdown id: tags:

 The first column labeled "Moment" defines the collision space and thus the transformation matrix $C$.
 The remaining columns specify the equilibrium vector in moment space $c^{(eq)}$ and the corresponding relaxation rate.
@@ -129,11 +131,11 @@
 method.stencil
 ```

 %%%% Output: execute_result

-    <lbmpy.stencils.LBStencil at 0x7ff400af3910>
+    <lbmpy.stencils.LBStencil at 0x7f0d442f2a90>

 %% Cell type:markdown id: tags:

 ### Orthogonal MRTs

@@ -144,28 +146,30 @@


 %% Cell type:code id: tags:

 ``` python
-weighted_ortho_mrt = create_lb_method(stencil=Stencil.D2Q9, method="mrt", weighted=True)
+lbm_config = LBMConfig(stencil=Stencil.D2Q9, method=Method.MRT, weighted=True)
+weighted_ortho_mrt = create_lb_method(lbm_config=lbm_config)
 weighted_ortho_mrt
 ```

 %%%% Output: execute_result

-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7ff400771250>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f0d39ecaee0>

 %% Cell type:code id: tags:

 ``` python
-ortho_mrt = create_lb_method(stencil=Stencil.D2Q9, method="mrt", weighted=False)
+lbm_config = LBMConfig(stencil=Stencil.D2Q9, method=Method.MRT, weighted=False)
+ortho_mrt = create_lb_method(lbm_config=lbm_config)
 ortho_mrt
 ```

 %%%% Output: execute_result

-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7ff4007f0d00>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f0d39f23760>

 %% Cell type:markdown id: tags:

 One can check if a method is orthogonalized:

@@ -189,18 +193,20 @@
 Another popular method is the cascaded lattice Boltzmann method. The cascaded LBM increases the numerical stability by shifting the collision step to the central moment space. Thus it is applied in the non-moving frame and achieves a better Galilean invariance. Typically the central moment collision operator is derived for compressible LB methods, and a higher-order equilibrium is used. Although incompressible LB methods with a second-order equilibrium can be derived with lbmpy, it violates the Galilean invariance and thus reduces the advantages of the method.

 %% Cell type:code id: tags:

 ``` python
-central_moment_method = create_lb_method(stencil=Stencil.D2Q9, method="central_moment",
-                                         equilibrium_order=4, compressible=True)
+lbm_config = LBMConfig(stencil=Stencil.D2Q9, method=Method.CENTRAL_MOMENT, equilibrium_order=4,
+                       compressible=True)
+
+central_moment_method = create_lb_method(lbm_config)
 central_moment_method
 ```

 %%%% Output: execute_result

-    <lbmpy.methods.momentbased.centralmomentbasedmethod.CentralMomentBasedLbMethod at 0x7ff400771910>
+    <lbmpy.methods.momentbased.centralmomentbasedmethod.CentralMomentBasedLbMethod at 0x7f0d39c84d90>

 %% Cell type:markdown id: tags:

 The shift to the central moment space is done by applying a so-called shift matrix. Usually, this introduces a high numerical overhead. This problem is solved with lbmpy because each transformation stage can be specifically optimised individually. Therefore, it is possible to derive a CLBM with only a little numerical overhead.

@@ -307,16 +313,16 @@
 This nested moment list can be passed to `create_lb_method`:

 %% Cell type:code id: tags:

 ``` python
-create_lb_method(stencil=Stencil.D2Q9, method="mrt", nested_moments=moments)
+create_lb_method(LBMConfig(stencil=Stencil.D2Q9, method=Method.MRT, nested_moments=moments))
 ```

 %%%% Output: execute_result

-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7ff4005d7160>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f0d39cceac0>

 %% Cell type:markdown id: tags:

 If one needs to also specify custom equilibrium moments the following approach can be used

@@ -347,11 +353,11 @@
 method
 ```

 %%%% Output: execute_result

-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7ff400465f10>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f0d39bf1c10>

 %% Cell type:markdown id: tags:

 Instead of manually defining all entries in the method table, *lbmpy* has functions to fill the table according to a specific pattern. For example:
 - for a full stencil (D2Q9, D3Q27) there exist exactly 9 or 27 linearly independent moments. These can either be taken as they are, or orthogonalized using Gram-Schmidt, weighted Gram-Schmidt or a Hermite approach
@@ -373,11 +379,11 @@
 method
 ```

 %%%% Output: execute_result

-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7ff4004bfac0>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f0d39bfbc70>

 %% Cell type:markdown id: tags:

 Our customized method can be directly passed into one of the scenarios. We can for example set up a channel flow with it. Since we used symbols as relaxation rates, we have to pass them in as `kernel_params`.


--- a/doc/notebooks/04_tutorial_cumulant_LBM.ipynb
+++ b/doc/notebooks/04_tutorial_cumulant_LBM.ipynb
@@ -91,34 +91,36 @@
 When creating a monomial cumulant method, one option is to specify only a single relaxation rate which will then be assigned to all cumulants related to the shear viscosity. In this case, all other non-conserved cumulants will be relaxed to equilibrium. Alternatively, individual relaxation rates for all non-conserved cumulants can be specified. The conserved cumulants are set to zero by default to save computational cost. They can be adjusted with `set_zeroth_moment_relaxation_rate`, `set_first_moment_relaxation_rate` and `set_conserved_moments_relaxation_rate`.
  
 %% Cell type:code id: tags:
  
 ``` python
-method_monomial = create_lb_method(method='monomial_cumulant', relaxation_rate=sp.Symbol('omega_v'))
+lbm_config = LBMConfig(method=Method.MONOMIAL_CUMULANT, relaxation_rate=sp.Symbol('omega_v'))
+method_monomial = create_lb_method(lbm_config=lbm_config)
 method_monomial
 ```
  
 %%%% Output: execute_result
  
-    <lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x7f52c30150a0>
+    <lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x7f1b2d7d6070>
  
 %% Cell type:markdown id: tags:
  
 ### Relaxation of Polynomial Cumulants
  
 By setting `method="cumulant"`, a set of default polynomial cumulants is chosen to be relaxed. Those cumulants are taken from literature and assembled into groups selected to enforce rotational invariance (see: `lbmpy.methods.centeredcumulant.centered_cumulants`). Default polynomial groups are available for the D2Q9, D3Q19 and D3Q27 stencils. As before it is possible to specify only a single relaxation rate assigned to the moments governing the shear viscosity. All other relaxation rates are then automatically set to one.
  
 %% Cell type:code id: tags:
  
 ``` python
-method_polynomial = create_lb_method(method='cumulant', relaxation_rate=sp.Symbol('omega_v'))
+lbm_config = LBMConfig(method=Method.CUMULANT, relaxation_rate=sp.Symbol('omega_v'))
+method_polynomial = create_lb_method(lbm_config=lbm_config)
 method_polynomial
 ```
  
 %%%% Output: execute_result
  
-    <lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x7f52d37f5940>
+    <lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x7f1b339eeb20>
  
 %% Cell type:markdown id: tags:
  
 ### Central Moments and Forcing
  
@@ -139,11 +141,11 @@
 method_with_force
 ```
  
 %%%% Output: execute_result
  
-    <lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x7f52c2fd6880>
+    <lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x7f1b33d73460>
  
 %% Cell type:markdown id: tags:
  
 ## C) Exemplary simulation: flow around sphere
  
@@ -209,17 +211,20 @@
 We choose a cumulant lattice Boltzmann method, as described above. Here the second-order cumulants are relaxed with the relaxation rate calculated above. All higher-order cumulants are relaxed with one, which means we set them to the equilibrium.
  
 %% Cell type:code id: tags:
  
 ``` python
-method = create_lb_method(stencil=Stencil.D2Q9, method='cumulant', relaxation_rate=omega)
+lbm_config = LBMConfig(stencil=Stencil.D2Q9, method=Method.CUMULANT, relaxation_rate=omega,
+                       output={'velocity': velField}, kernel_type='stream_pull_collide')
+
+method = create_lb_method(lbm_config=lbm_config)
 method
 ```
  
 %%%% Output: execute_result
  
-    <lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x7f52c2fa8fd0>
+    <lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x7f1b2d717280>
  
 %% Cell type:markdown id: tags:
  
 ### Initialisation with equilibrium
  
@@ -239,15 +244,14 @@
 ### Definition of the update rule
  
 %% Cell type:code id: tags:
  
 ``` python
+lbm_optimisation = LBMOptimisation(symbolic_field=src, symbolic_temporary_field=dst)
 update = create_lb_update_rule(lb_method=method,
-                               output={'velocity': velField},
-                               optimization={"symbolic_field": src,
-                                             "symbolic_temporary_field": dst},
-                               kernel_type='stream_pull_collide')
+                               lbm_config=lbm_config,
+                               lbm_optimisation=lbm_optimisation)
  
 ast_kernel = ps.create_kernel(update, target=dh.default_target, cpu_openmp=True)
 kernel = ast_kernel.compile()
 ```
  
--- a/doc/notebooks/06_tutorial_modifying_method_smagorinsky.ipynb
+++ b/doc/notebooks/06_tutorial_modifying_method_smagorinsky.ipynb
@@ -241,18 +241,20 @@
 We create a MRT method, where we fix all relaxation rates except the relaxation rate that controls the viscosity.
  
 %% Cell type:code id: tags:
  
 ``` python
-method = create_lb_method(method='mrt', weighted=True, stencil=Stencil.D2Q9, force=(1e-6, 0),
-                          force_model='luo', relaxation_rates=[0, 0, ω, 1.9, 1.9])
+lbm_conifg = LBMConfig(stencil=Stencil.D2Q9, method=Method.MRT, force=(1e-6, 0),
+                       force_model=ForceModel.LUO, relaxation_rates=[0, 0, ω, 1.9, 1.9])
+
+method = create_lb_method(lbm_config=lbm_conifg)
 method
 ```
  
 %%%% Output: execute_result
  
-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f0b1feea940>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7fc79c0cb7f0>
  
 %% Cell type:markdown id: tags:
  
 Only the collision rule has to be changed. Thus we first construct the collision rule, add the Smagorinsky equations and create a normal scenario from the modified collision rule. To avoid that the macroscopic quantity symbols in the Smagorinsky equations fall prey to optimization, we must disable simplification:
  
@@ -268,11 +270,11 @@
 collision_rule
 ```
  
 %%%% Output: execute_result
  
-    AssignmentCollection: d_2, d_7, d_4, d_6, d_0, d_8, d_3, d_1, d_5 <- f(f_2, f_3, f_1, f_7, f_4, f_0, omega_total, f_8, f_5, f_6)
+    AssignmentCollection: d_1, d_4, d_6, d_0, d_3, d_2, d_7, d_8, d_5 <- f(f_7, f_2, f_0, f_5, f_3, f_8, omega_total, f_1, f_4, f_6)
  
 %% Cell type:markdown id: tags:
  
 In the next cell the collision rule is simplified by extracting common subexpressions
  

--- a/doc/notebooks/08_tutorial_shanchen_twophase.ipynb
+++ b/doc/notebooks/08_tutorial_shanchen_twophase.ipynb
@@ -81,16 +81,14 @@
 ## Kernels

 %% Cell type:code id: tags:

 ``` python
-collision = create_lb_update_rule(stencil=stencil,
-                                  relaxation_rate=omega_a, 
-                                  compressible=True,
-                                  force_model='guo', 
-                                  force=force,
-                                  kernel_type='collide_only',
+lbm_config = LBMConfig(stencil=stencil, relaxation_rate=omega_a, compressible=True,
+                       force_model=ForceModel.GUO, force=force, kernel_type='collide_only')
+
+collision = create_lb_update_rule(lbm_config=lbm_config,
                                  optimization={'symbolic_field': src})

 stream = create_stream_pull_with_output_kernel(collision.method, src, dst, {'density': ρ})


@@ -106,11 +104,11 @@
 ## Initialization

 %% Cell type:code id: tags:

 ``` python
-method_without_force = create_lb_method(stencil=stencil, relaxation_rate=omega_a, compressible=True)
+method_without_force = create_lb_method(LBMConfig(stencil=stencil, relaxation_rate=omega_a, compressible=True))
 init_assignments = macroscopic_values_setter(method_without_force, velocity=(0, 0), 
                                             pdfs=src.center_vector, density=ρ.center)


 init_kernel = ps.create_kernel(init_assignments, ghost_layers=0).compile()

--- a/doc/notebooks/09_tutorial_shanchen_twocomponent.ipynb
+++ b/doc/notebooks/09_tutorial_shanchen_twocomponent.ipynb
@@ -122,26 +122,24 @@
 ## Kernels

 %% Cell type:code id: tags:

 ``` python
-collision_a = create_lb_update_rule(stencil=stencil,
-                                    relaxation_rate=omega_a, 
-                                    compressible=True,
-                                    velocity_input=u_full, density_input=ρ_a,
-                                    force_model='guo', 
-                                    force=force_a,
-                                    kernel_type='collide_only',
+lbm_config_a = LBMConfig(stencil=stencil, relaxation_rate=omega_a, compressible=True,
+                         velocity_input=u_full, density_input=ρ_a, force_model=ForceModel.GUO, 
+                         force=force_a, kernel_type='collide_only')
+
+lbm_config_b = LBMConfig(stencil=stencil, relaxation_rate=omega_b, compressible=True,
+                         velocity_input=u_full, density_input=ρ_b, force_model=ForceModel.GUO, 
+                         force=force_b, kernel_type='collide_only')
+
+
+
+collision_a = create_lb_update_rule(lbm_config=lbm_config_a,
                                    optimization={'symbolic_field': src_a})

-collision_b = create_lb_update_rule(stencil=stencil,
-                                    relaxation_rate=omega_b,
-                                    compressible=True,
-                                    velocity_input=u_full, density_input=ρ_b,
-                                    force_model='guo', 
-                                    force=force_b,
-                                    kernel_type='collide_only',
+collision_b = create_lb_update_rule(lbm_config=lbm_config_b,
                                    optimization={'symbolic_field': src_b})
 ```

 %% Cell type:code id: tags:


--- a/doc/notebooks/10_tutorial_conservative_allen_cahn_two_phase.ipynb
+++ b/doc/notebooks/10_tutorial_conservative_allen_cahn_two_phase.ipynb
@@ -190,39 +190,47 @@
 For both LB steps, a weighted orthogonal MRT (WMRT) method is used. It is also possible to change the method to a simpler SRT scheme or a more complicated CLBM scheme. The CLBM scheme can be obtained by commenting in the python snippets in the notebook cells below. Note here that the hydrodynamic LB step is formulated as an incompressible velocity-based LBM. Thus, the velocity terms can not be removed from the equilibrium in the central moment space.

 %% Cell type:code id: tags:

 ``` python
-method_phase = create_lb_method(stencil=stencil_phase, method="mrt", compressible=True, weighted=True,
-                                relaxation_rates=[1, 1, 1, 1])
+config_phase = LBMConfig(stencil=stencil_phase, method=Method.MRT, compressible=True,
+                         weighted=True, relaxation_rates=[1, 1, 1, 1, 1])

+method_phase = create_lb_method(lbm_config=config_phase)
 method_phase.set_first_moment_relaxation_rate(w_c)

-# method_phase = create_lb_method(stencil=stencil_phase, method="central_moment", compressible=True,
-#                                 relaxation_rates=[0, w_c, 1, 1, 1], equilibrium_order=4)
+# config_phase = LBMConfig(stencil=stencil_phase, method=Method.CENTRAL_MOMENT, compressible=True,
+#                          weighted=True, relaxation_rates=[0, w_c, 1, 1, 1], equilibrium_order=4)
+
+# method_phase = create_lb_method(lbm_config=config_phase)
 method_phase
 ```

 %%%% Output: execute_result

-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7efec8cbdee0>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f5d7c26fdc0>

 %% Cell type:code id: tags:

 ``` python
-method_hydro = create_lb_method(stencil=stencil_phase, method="mrt", compressible=False, weighted=True,
-                                relaxation_rates=[s8, 1, 1, 1])
+config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.MRT, compressible=False,
+                         weighted=True, relaxation_rates=[s8, 1, 1, 1])
+
+method_hydro = create_lb_method(lbm_config=config_hydro)
+
+
+# config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.CENTRAL_MOMENT, compressible=False,
+#                          weighted=True, relaxation_rates=[s8, 1, 1], equilibrium_order=4)

-# method_hydro = create_lb_method(stencil=stencil_phase, method="central_moment", compressible=False,
-#                                 relaxation_rates=[s8, 1, 1], equilibrium_order=4)
+# method_hydro = create_lb_method(lbm_config=config_hydro)

 method_hydro
 ```

 %%%% Output: execute_result

-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7efec677b760>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f5d7c249fa0>

 %% Cell type:markdown id: tags:

 ## Initialization

@@ -279,11 +287,11 @@
 plt.scalar_field(dh.gather_array(C.name))
 ```

 %%%% Output: execute_result

-    <matplotlib.image.AxesImage at 0x7efec6339160>
+    <matplotlib.image.AxesImage at 0x7f5d758dcac0>

 %%%% Output: display_data

    ![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAA64AAAFlCAYAAADrtrUsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbkUlEQVR4nO3df7DddX3n8df75ibhl5IAMQLBBiVV6Q+pjZSq46i0CnSn0K51dXaUseyws6Pd7tqZLd3dGXfa3f7Y6dbWbtcdtlhwp611aa2MtVqKWmstSLQUREQigiRCEkISfuTnveezf9xv7BUJkNybez+59/GYuXO/53O+55zPzXzm3vvM93u+t1prAQAAgF6NzfcEAAAA4OkIVwAAALomXAEAAOiacAUAAKBrwhUAAICuCVcAAAC6Nj7fE3g6p512Wlu7du18TwOAheLAl+d7BgvD0u+f7xkAsEB88YtffLi1tuqZ9us6XNeuXZsNGzbM9zQAWCBGD62b7yksCGPP97MZgNlRVfc/m/2cKgwAAEDXhCsAAABde8ZwraoPVNXWqvrytLFTqurGqrpn+LxyGK+qel9Vbayq26vq5dMec/mw/z1VdfnR+XIAAABYaJ7NEddrk1z0pLGrktzUWluX5KbhdpJcnGTd8HFlkvcnU6Gb5D1JfiTJ+UneczB2AQAA4Ok8Y7i21j6b5JEnDV+a5Lph+7okl00b/2CbcnOSFVV1epI3JrmxtfZIa21Hkhvz3TEMAAAA3+VI3+O6urX24LD9UJLVw/aZSR6Ytt+mYexQ49+lqq6sqg1VtWHbtm1HOD0AAAAWihlfnKm11pK0WZjLwee7urW2vrW2ftWqZ/xzPgAAACxwRxquW4ZTgDN83jqMb05y1rT91gxjhxoHAACAp3Wk4XpDkoNXBr48yUenjb99uLrwBUl2DacUfzLJG6pq5XBRpjcMYwAAAPC0xp9ph6r64ySvTXJaVW3K1NWBfz3Jh6vqiiT3J3nzsPvHk1ySZGOS3UnekSSttUeq6leS3Drs98uttSdf8AkAAAC+yzOGa2vtrYe468Kn2LcleechnucDST5wWLMDAABg0ZvxxZkAAADgaBKuAAAAdE24AgAA0DXhCgAAQNeEKwAAAF0TrgAAAHRNuAIAANA14QoAAEDXhCsAAABdE64AAAB0TbgCAADQNeEKAABA14QrAAAAXROuAAAAdE24AgAA0DXhCgAAQNeEKwAAAF0TrgAAAHRNuAIAANA14QoAAEDXhCsAAABdE64AAAB0TbgCAADQNeEKAABA14QrAAAAXROuAAAAdE24AgAA0DXhCgAAQNeEKwAAAF0TrgAAAHRNuAIAANA14QoAAEDXhCsAAABdE64AAAB0TbgCAADQNeEKAABA14QrAAAAXROuAAAAdE24AgAA0DXhCgAAQNeEKwAAAF0TrgAAAHRNuAIAANA14QoAAEDXhCsAAABdE64AAAB0TbgCAADQNeEKAABA14QrAAAAXROuAAAAdE24AgAA0DXhCgAAQNeEKwAAAF0TrgAAAHRtRuFaVf++qu6sqi9X1R9X1XFVdXZV3VJVG6vqT6pq2bDv8uH2xuH+tbPyFQAAALCgHXG4VtWZSf5tkvWtte9PsiTJW5L8RpL3ttbOSbIjyRXDQ65IsmMYf++wHwAAADytmZ4qPJ7k+KoaT3JCkgeTvD7J9cP91yW5bNi+dLid4f4Lq6pm+PoAAAAscEccrq21zUl+M8k3MxWsu5J8McnO1trEsNumJGcO22cmeWB47MSw/6lH+voAAAAsDjM5VXhlpo6inp3kjCQnJrlophOqqiurakNVbdi2bdtMnw4AAIBj3ExOFf6xJN9orW1rrR1I8mdJXpVkxXDqcJKsSbJ52N6c5KwkGe4/Ocn2Jz9pa+3q1tr61tr6VatWzWB6AAAALAQzCddvJrmgqk4Y3qt6YZKvJPl0kjcN+1ye5KPD9g3D7Qz3f6q11mbw+gAAACwCM3mP6y2ZusjSl5LcMTzX1Ul+Mcm7q2pjpt7Des3wkGuSnDqMvzvJVTOYNwAAAIvE+DPvcmittfckec+Thu9Ncv5T7Ls3yc/M5PUAAABYfGb653AAAADgqBKuAAAAdE24AgAA0DXhCgAAQNeEKwAAAF0TrgAAAHRNuAIAANA14QoAAEDXhCsAAABdE64AAAB0TbgCAADQNeEKAABA14QrAAAAXROuAAAAdE24AgAA0DXhCgAAQNeEKwAAAF0TrgAAAHRNuAIAANA14QoAAEDXhCsAAABdE64AAAB0TbgCAADQNeEKAABA14QrAAAAXROuAAAAdE24AgAA0DXhCgAAQNeEKwAAAF0TrgAAAHRNuAIAANA14QoAAEDXhCsAAABdE64AAAB0TbgCAADQNeEKAABA14QrAAAAXROuAAAAdE24AgAA0DXhCgAAQNeEKwAAAF0TrgAAAHRNuAIAANA14QoAAEDXhCsAAABdE64AAAB0TbgCAADQNeEKAABA14QrAAAAXROuAAAAdE24AgAA0DXhCgAAQNeEKwAAAF0TrgAAAHRNuAIAANC1GYVrVa2oquur6qtVdVdV/WhVnVJVN1bVPcPnlcO+VVXvq6qNVXV7Vb18dr4EAAAAFrKZHnH9nSSfaK29JMnLktyV5KokN7XW1iW5abidJBcnWTd8XJnk/TN8bQAAABaBIw7Xqjo5yWuSXJMkrbX9rbWdSS5Nct2w23VJLhu2L03ywTbl5iQrqur0I319AAAAFoeZHHE9O8m2JH9QVf9QVb9fVScmWd1ae3DY56Ekq4ftM5M8MO3xm4YxAAAAOKSZhOt4kpcneX9r7YeSPJF/Oi04SdJaa0na4TxpVV1ZVRuqasO2bdtmMD0AAAAWgpmE66Ykm1prtwy3r89UyG45eArw8HnrcP/mJGdNe/yaYew7tNaubq2tb62tX7Vq1QymBwAAwEJwxOHaWnsoyQNV9eJh6MIkX0lyQ5LLh7HLk3x02L4hyduHqwtfkGTXtFOKAQAA4CmNz/DxP5fkD6tqWZJ7k7wjUzH84aq6Isn9Sd487PvxJJck2Zhk97AvAAAAPK0ZhWtr7bYk65/irgufYt+W5J0zeT0AAAAWn5n+HVcAAAA4qoQrAAAAXROuAAAAdE24AgAA0DXhCgAAQNeEKwAAAF0TrgAAAHRNuAIAANA14QoAAEDXhCsAAABdE64AAAB0TbgCAADQtfH5ngAAzJV97UCSZCxjGUslSZaU/8N9OpNtlCQZpWWUqe3j53NCACxKwhWAReP6x5+f4+pAVix5IqeO7c7JYweyYmwsJ4wtzXiWiNjBZBtlIpPZPTqQnaNRdo2WZvvohOycPDF729K8bb4nCMCiI1wBWDR+9dp/kdHy5MCJLZOnHshpz3s0563anFedfE/OW/5A1oxP5DljyxZlxB6M1cdG+7NpYjy37Tsrf7drXW7bdmYe3vLcLHlkaZY+Xhnbn7ztv833bAFYbIQrAIvG2j/4ejI+niwdz+jE4zOx8sTcecYP5PMvfFn2vGRvXrnu3vzUaV/Ky5d/K6uXLMvyGl/wATvZRtnXJrJlcn++tO+M/Om2H87N97wwx9+9PCffO8qKb+3LaTsez9gTe5L9B5LJyUS4AjDHhCsAi8bElq1TGzWWGqssWbIkJy9blpUnnpDR81bmvnUvzn962UuzYv22/OzZn8+FJ3wtpy/QgD0YrA9O7s+NT7w41973o9m5YVVO+8dRXnzProxt3ZTRE7vT9u9Pm5zMxKglw/tdAWCuCVcAFo/Whs+TUw02MZG2b19Gjz+eenh7nvP15Tn5Cyuz79PPy/vOvyz/+1Xb8851f5MfP3FjVi9ZnuW1dF6nP1v2tQPZMrkvNz5xTv7n116bib87Jatv3ZtTNz6Q0SM70vbty8Tk5D/9ewHAPBOuANBa2sRE2sRERnv2ZnzrtnzPXSuy//Nn5Lde/dP5v6/7Vt599l/llcdty8ljx2VpLZnvGR+RA20yu0Z78/m9q/Kb974x2z9zes743J4su+frmdyxMxMHJpLR5HxPEwC+i3AFgOlGk2n7JjOxZWuW7NiZtRtXZveX1uQXXnd5Xv3aL+ddq2/Ki5ceyPG17Jg5fXiyjbKn7c/dB8byOw9dnJv/5vuy5lMHsvaOezP5yI5M7N/v6CoAXROuAPBUWps6ZfahLTlux85879dW5+7bvi9veuP35ude8em86bm35/Qlx3d/9PVAm8yDk3ty/aM/mN+99XU5/RNLs+4L38rowS2Z2LdPsAJwTBCuAPB0Wsto796M7t+UFdt35OSvrsm1r78of/7Gl+U/v+hjecXyXXnu2HHdHX2dbKM8OtqbW/atzK9+/U3Z9YnT85JP70zde08mHn/CKcEAHFOEKwA8G6PJTD76aOrOjVnz4IrsufMF+TcXXZHLXvOF/KtTP5cXLl3azcWbdo/25/6Jifz+9lfnz//m/LzgkxM587aNGW1/JKOJifmeHgAcNuEKAIehHdifya3bctzfPp4Xf+P0fPaOH8lfXnJurvr+T+THT7gvp83j6cMH2mQentyTG3evza/dcXFO+suT8uK/3ZK2+aFM7tnjtGAAjlnCFQAOV2sZ7d6dbLwvq7Y9kpVffUF+7cfenD+78Ot591mfzMuW7clJtXzOTh+ebKM83vblH/cfn9964J/na3/9opx10xMZv+urmdz1qNOCATjmCVcAOFKjyUzu3JmxL+7OC795WrbfsTaX//iVedMFt+Ztp/x9Xjieo3r14YNXC753Irlu+6vykZtfkTU3tpy94f5Mbns4k64WDMACIVwBYCYOXn34Ww/lpL/amZd+9Yx8ZsMF+ciF5+Vnf+DzufS5t+V7xivLa+msnUJ8oE1mXzuQ+ydaPvroebnm9lfmlJuOy0tufjjtm9/KxJ69jrICsKAIVwCYDaPJjJ54InXPN7LqoYdz6j+ekf93wYW59tUX5K3nbshPPPe2vHB8b54ztizjWXLYR2En2ygTmcxjo/25d2JZ/uLRH84f3fmKnPS5E7Lull0Zu/e+jB57LM3FlwBYgIQrAMyiNjGRyZ07U7c/kdPvf05GXzgjH/uh1+SDr3hlfvTcjfnJ027LDyz/VlYtGeWEWpKltSRjGctY6jueZ5SWUUY50Cazu01m2+RY7th3Rm54+Lz8/VfOySm3judFX3osY/fdndGuxzI5ccBpwQAsWMIVAGZba1NXH97+SGrXY1l970lZ/Xer8s2XfG9++dyXZv9L9uQHz9qU9Su+mXOOeyjPH9+V59a+LK1RkmRvW5In2rI8NHFKNu59fjbsfEFuf2BNlt11fE79ymReevfO5KFtGT36uGAFYFEQrgBwtBwM2Ed2pHY9mufctzwn33xyJp+3MtvXnJ3rzzwne55f2XfKZNpJkxlbNvW+1NG+JaknlmT59iU5fkvLczZP5pxNT2TJ1i0Z7dyV0b59aZOTghWARUO4AsDR1lraxETaxERGu3entmzNCXcvy4nHH5c67ri045cny5eljU+977UmRsm+/ak9+9L27k3bszdt//5MiFUAFinhCgBzaVrEZs+eqbEaS43903tcW5I2akkbffsxALCYCVcAmC8Hg7RNfrtRAYDvdnT+IjoAAADMEuEKAABA14QrAAAAXROuAAAAdE24AgAA0DXhCgAAQNeEKwAAAF0TrgAAAHRNuAIAANA14QoAAEDXhCsAAABdE64AAAB0TbgCAADQNeEKAABA14QrAAAAXROuAAAAdE24AgAA0DXhCgAAQNeEKwAAAF0TrgAAAHRNuAIAANC1GYdrVS2pqn+oqo8Nt8+uqluqamNV/UlVLRvGlw+3Nw73r53pawMAALDwzcYR159Pcte027+R5L2ttXOS7EhyxTB+RZIdw/h7h/0AAADgac0oXKtqTZKfSPL7w+1K8vok1w+7XJfksmH70uF2hvsvHPYHAACAQ5rpEdffTvIfkoyG26cm2dlamxhub0py5rB9ZpIHkmS4f9ew/3eoqiurakNVbdi2bdsMpwcAAMCx7ojDtar+WZKtrbUvzuJ80lq7urW2vrW2ftWqVbP51AAAAByDxmfw2Fcl+cmquiTJcUmem+R3kqyoqvHhqOqaJJuH/TcnOSvJpqoaT3Jyku0zeH0AAAAWgSM+4tpa+6XW2prW2tokb0nyqdbav0zy6SRvGna7PMlHh+0bhtsZ7v9Ua60d6esDAACwOByNv+P6i0neXVUbM/Ue1muG8WuSnDqMvzvJVUfhtQEAAFhgZnKq8Le11j6T5DPD9r1Jzn+KffYm+ZnZeD0AAAAWj6NxxBUAAABmjXAFAACga8IVAACArglXAAAAuiZcAQAA6JpwBQAAoGvCFQAAgK4JVwAAALomXAEAAOiacAUAAKBrwhUAAICuCVcAAAC6JlwBAADomnAFAACga8IVAACArglXAAAAuiZcAQAA6JpwBQAAoGvCFQAAgK4JVwAAALomXAEAAOiacAUAAKBrwhUAAICuCVcAAAC6JlwBAADomnAFAACga8IVAACArglXAAAAuiZcAQAA6JpwBQAAoGvCFQAAgK4JVwAAALomXAEAAOiacAUAAKBrwhUAAICuCVcAAAC6JlwBAADomnAFAACga8IVAACArglXAAAAuiZcAQAA6JpwBQAAoGvCFQAAgK4JVwAAALomXAEAAOiacAUAAKBrwhUAAICuCVcAAAC6JlwBAADomnAFAACga8IVAACArglXAAAAuiZcAQAA6JpwBQAAoGtHHK5VdVZVfbqqvlJVd1bVzw/jp1TVjVV1z/B55TBeVfW+qtpYVbdX1ctn64sAAABg4ZrJEdeJJL/QWjs3yQVJ3llV5ya5KslNrbV1SW4abifJxUnWDR9XJnn/DF4bAACAReKIw7W19mBr7UvD9mNJ7kpyZpJLk1w37HZdksuG7UuTfLBNuTnJiqo6/UhfHwAAgMVhVt7jWlVrk/xQkluSrG6tPTjc9VCS1cP2mUkemPawTcMYAAAAHNKMw7WqTkryp0n+XWvt0en3tdZaknaYz3dlVW2oqg3btm2b6fQAAAA4xs0oXKtqaaai9Q9ba382DG85eArw8HnrML45yVnTHr5mGPsOrbWrW2vrW2vrV61aNZPpAQAAsADM5KrCleSaJHe11n5r2l03JLl82L48yUenjb99uLrwBUl2TTulGAAAAJ7S+Awe+6okb0tyR1XdNoz9xyS/nuTDVXVFkvuTvHm47+NJLkmyMcnuJO+YwWsDAACwSBxxuLbWPpekDnH3hU+xf0vyziN9PQAAABanWbmqMAAAABwtwhUAAICuCVcAAAC6JlwBAADomnAFAACga8IVAACArglXAAAAuiZcAQAA6JpwBQAAoGvCFQAAgK4JVwAAALomXAEAAOiacAUAAKBrwhUAAICuCVcAAAC6JlwBAADomnAFAACga8IVAACArglXAAAAuiZcAQAA6JpwBQAAoGvCFQAAgK4JVwAAALomXAEAAOiacAUAAKBrwhUAAICuCVcAAAC6JlwBAADomnAFAACga8IVAACArglXAAAAuiZcAQAA6JpwBQAAoGvCFQAAgK4JVwAAALomXAEAAOiacAUAAKBrwhUAAICuCVcAAAC6JlwBAADomnAFAACga8IVAACArglXAAAAuiZcAQAA6JpwBQAAoGvCFQAAgK4JVwAAALomXAEAAOiacAUAAKBrwhUAAICuCVcAAAC6JlwBAADomnAFAACga8IVAACArglXAAAAujbn4VpVF1XV3VW1saqumuvXBwAA4Ngyp+FaVUuS/F6Si5Ocm+StVXXuXM4BAACAY8tcH3E9P8nG1tq9rbX9ST6U5NI5ngMAAADHkLkO1zOTPDDt9qZhDAAAAJ7S+HxP4Mmq6sokVw43H6+qu+dzPgvcaUkenu9JsChZe8wXa28WVNV8T+FYZO0xn6w/5suzWXvf82yeaK7DdXOSs6bdXjOMfVtr7eokV8/lpBarqtrQWls/3/Ng8bH2mC/WHvPF2mM+WX/Ml9lce3N9qvCtSdZV1dlVtSzJW5LcMMdzAAAA4Bgyp0dcW2sTVfWuJJ9MsiTJB1prd87lHAAAADi2zPl7XFtrH0/y8bl+XZ6SU7KZL9Ye88XaY75Ye8wn64/5Mmtrr1prs/VcAAAAMOvm+j2uAAAAcFiE6wJVVT9TVXdW1aiq1j/pvl+qqo1VdXdVvXHa+EXD2Maqumra+NlVdcsw/ifDhbXgGVXVf6mqzVV12/BxybT7DmsdwkxZWxxtVXVfVd0xfL/bMIydUlU3VtU9w+eVw3hV1fuG9Xh7Vb18fmfPsaSqPlBVW6vqy9PGDnutVdXlw/73VNXl8/G1cGw5xNqbk9/3hOvC9eUkP53ks9MHq+rcTF3N+fuSXJTkf1XVkqpakuT3klyc5Nwkbx32TZLfSPLe1to5SXYkuWJuvgQWiPe21s4bPj6eHPE6hCNmbTGHXjd8vzv4n8ZXJbmptbYuyU3D7WRqLa4bPq5M8v45nynHsmsz9fNzusNaa1V1SpL3JPmRJOcnec/B2IWncW2+e+0lc/D7nnBdoFprd7XW7n6Kuy5N8qHW2r7W2jeSbMzUN6vzk2xsrd3bWtuf5ENJLq2pvzL/+iTXD4+/LsllR/0LYKE7rHU4j/Nk4bC2mC+XZupnZ/KdP0MvTfLBNuXmJCuq6vR5mB/HoNbaZ5M88qThw11rb0xyY2vtkdbajiQ35qmDBL7tEGvvUGb19z3huvicmeSBabc3DWOHGj81yc7W2sSTxuHZetdwatIHpv1P7uGuQ5gpa4u50JL8VVV9saquHMZWt9YeHLYfSrJ62LYmmW2Hu9asQWbTUf99T7gew6rqr6vqy0/x4SgCc+YZ1uH7k7woyXlJHkzyP+ZzrgBH2atbay/P1Olv76yq10y/s039KQd/zoGjzlpjjs3J73tz/ndcmT2ttR87godtTnLWtNtrhrEcYnx7pk4pGR+Ouk7fH571Oqyq/5PkY8PNw12HMFNPt+ZgVrTWNg+ft1bVRzJ1OtyWqjq9tfbgcHrm1mF3a5LZdrhrbXOS1z5p/DNzME8WmNbaloPbR/P3PUdcF58bkrylqpZX1dmZeqP+F5LcmmRdTV1BeFmm3kh9w/A/dp9O8qbh8Zcn+eg8zJtj0JPer/VTmbpoWHKY63Au58yCZW1xVFXViVX1nIPbSd6Qqe95N2TqZ2fynT9Db0jy9uGKrxck2TXtNE84Eoe71j6Z5A1VtXI4tfMNwxgclrn6fc8R1wWqqn4qye8mWZXkL6rqttbaG1trd1bVh5N8JclEkne21iaHx7wrU9+wliT5QGvtzuHpfjHJh6rqvyb5hyTXzPGXw7Hrv1fVeZk6Xem+JP86SY5wHcIRa61NWFscZauTfGTqmoYZT/JHrbVPVNWtST5cVVckuT/Jm4f9P57kkkxdrGR3knfM/ZQ5VlXVH2fqaOlpVbUpU1cH/vUcxlprrT1SVb+SqYhIkl9urT3bi+6wSB1i7b12Ln7fq6kDagAAANAnpwoDAADQNeEKAABA14QrAAAAXROuAAAAdE24AgAA0DXhCgAAQNeEKwAAAF0TrgAAAHTt/wO4DKC00SnD+AAAAABJRU5ErkJggg==)


--- a/doc/notebooks/demo_automatic_chapman_enskog_analysis.ipynb
+++ b/doc/notebooks/demo_automatic_chapman_enskog_analysis.ipynb
 %% Cell type:code id: tags:

 ``` python
-import sympy as sp
+from lbmpy.session import *
 from lbmpy.chapman_enskog import ChapmanEnskogAnalysis
-from lbmpy.creationfunctions import create_lb_method
-from lbmpy.enums import Stencil
 sp.init_printing()
 ```

 %% Cell type:markdown id: tags:

@@ -16,17 +14,17 @@
 First, we create a SRT lattice Boltzmann method. It is defined as the set of moments, together with one relaxation rate per moment.

 %% Cell type:code id: tags:

 ``` python
-method = create_lb_method(method='trt', stencil=Stencil.D3Q19, compressible=False)
+method = create_lb_method(LBMConfig(method=Method.TRT, stencil=Stencil.D3Q19, compressible=False))
 method
 ```

 %%%% Output: execute_result

-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f66bb69d6a0>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f5ea617b190>

 %% Cell type:markdown id: tags:

 Next, the Chapman Enskog analysis object is created. This may take a while...


--- a/doc/notebooks/demo_thermalized_lbm.ipynb
+++ b/doc/notebooks/demo_thermalized_lbm.ipynb
@@ -20,17 +20,17 @@
  
 %% Cell type:code id: tags:
  
 ``` python
 random_number_symbols = sp.symbols("rand_:3")
-method = create_lb_method(method='srt', relaxation_rate=1.8)
+method = create_lb_method(LBMConfig(method=Method.SRT, relaxation_rate=1.8))
 method
 ```
  
 %%%% Output: execute_result
  
-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f6e6c2a6280>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7fb1641ba820>
  
 %% Cell type:markdown id: tags:
  
 ...and modify its collision table. The `create_lb_method_from_existing` function provides a convenient way to do this. 
 We pass a custom function that receives a row of the collision table and returns a modified version of it.
@@ -52,11 +52,11 @@
 thermalized_method
 ```
  
 %%%% Output: execute_result
  
-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f6e45ef77c0>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7fb15feace80>
  
 %% Cell type:markdown id: tags:
  
 ## 2) Creating the kernel equations
  
@@ -123,11 +123,11 @@
 collision_rule
 ```
  
 %%%% Output: execute_result
  
-    AssignmentCollection: d_7, d_8, d_2, d_6, d_3, d_1, d_5, d_0, d_4 <- f(f_8, f_1, rand_0, f_3, f_5, f_7, f_0, rand_1, f_6, f_2, f_4)
+    AssignmentCollection: d_0, d_8, d_4, d_3, d_7, d_6, d_5, d_2, d_1 <- f(rand_1, f_4, f_5, f_0, f_8, f_2, rand_0, f_6, f_1, f_3, f_7)
  
 %% Cell type:markdown id: tags:
  
 Finally, lets test our method by running a lid-driven-cavity with it.
  

--- a/lbmpy/creationfunctions.py
+++ b/lbmpy/creationfunctions.py
@@ -191,7 +191,7 @@ For example, to modify the AST one can run::

 """
 from collections import OrderedDict
-from dataclasses import dataclass, field
+from dataclasses import dataclass, field, replace
 from typing import Union, List, Tuple, Any, Type
 import warnings

@@ -534,7 +534,7 @@ def update_with_default_parameters(params, opt_params=None, lbm_config=None, lbm
        for k, v in config_params.items():
            if not hasattr(config, k):
                raise KeyError(f'{v} is not a valid kwarg. Please look in CreateKernelConfig for valid settings')
-            setattr(config, k, v)
+        config = replace(config, **config_params)

    lbm_opt_params = ['cse_pdfs', 'cse_global', 'simplification', 'pre_simplification', 'split', 'field_size',
                      'field_layout', 'symbolic_field', 'symbolic_temporary_field', 'builtin_periodicity']
@@ -550,7 +550,7 @@ def update_with_default_parameters(params, opt_params=None, lbm_config=None, lbm
        for k, v in opt_params_dict.items():
            if not hasattr(lbm_optimisation, k):
                raise KeyError(f'{v} is not a valid kwarg. Please look in LBMOptimisation for valid settings')
-            setattr(lbm_optimisation, k, v)
+        lbm_optimisation = replace(lbm_optimisation, **opt_params_dict)

    if params is None:
        params = {}
@@ -561,7 +561,7 @@ def update_with_default_parameters(params, opt_params=None, lbm_config=None, lbm
        for k, v in params.items():
            if not hasattr(lbm_config, k):
                raise KeyError(f'{v} is not a valid kwarg. Please look in LBMConfig for valid settings')
-            setattr(lbm_config, k, v)
+        lbm_config = replace(lbm_config, **params)

    return lbm_config, lbm_optimisation, config

@@ -732,7 +732,7 @@ class LBMConfig:
            force_model_class = force_model_dict[new_force_model.name.lower()]
            self.force_model = force_model_class(force=self.force[:self.stencil.D])
        elif isinstance(self.force_model, ForceModel):
-            force_model_class = force_model_dict[new_force_model.name.lower()]
+            force_model_class = force_model_dict[self.force_model.name.lower()]
            self.force_model = force_model_class(force=self.force[:self.stencil.D])



--- a/lbmpy/lbstep.py
+++ b/lbmpy/lbstep.py
@@ -5,9 +5,11 @@ import numpy as np
 from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
 from lbmpy.creationfunctions import (
    create_lb_function, update_with_default_parameters)
+from lbmpy.enums import Stencil
 from lbmpy.macroscopic_value_kernels import (
    create_advanced_velocity_setter_collision_rule, pdf_initialization_assignments)
 from lbmpy.simplificationfactory import create_simplification_strategy
+from lbmpy.stencils import LBStencil
 from pystencils import create_data_handling, create_kernel, make_slice, Target, Backend
 from pystencils.slicing import SlicedGetter
 from pystencils.timeloop import TimeLoop
@@ -21,7 +23,8 @@ class LatticeBoltzmannStep:
                 compute_velocity_in_every_step=False, compute_density_in_every_step=False,
                 velocity_input_array_name=None, time_step_order='stream_collide', flag_interface=None,
                 alignment_if_vectorized=64, fixed_loop_sizes=True, fixed_relaxation_rates=True,
-                 timeloop_creation_function=TimeLoop, **method_parameters):
+                 timeloop_creation_function=TimeLoop,
+                 lbm_config=None, lbm_optimisation=None, config=None, **method_parameters):

        self._timeloop_creation_function = timeloop_creation_function

@@ -30,7 +33,11 @@ class LatticeBoltzmannStep: