lbmpy_tests/full_scenarios/phasefield_allen_cahn/phasefield-capillary-wave.ipynb → tests/full_scenarios/phasefield_allen_cahn/phasefield-capillary-wave.ipynb

 %% Cell type:code id: tags:
 
 ``` python
 import os
 
 import time
 from tqdm.notebook import tqdm
 
 import numpy as np
 import math
 from scipy.special import erfc
 
 from pystencils.session import *
 from lbmpy.session import *
 
 from pystencils.simp import sympy_cse
 from pystencils.boundaries import BoundaryHandling
 
 from lbmpy.phasefield_allen_cahn.parameter_calculation import AllenCahnParameters
 from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle
 from lbmpy.phasefield_allen_cahn.kernel_equations import *
 
 from lbmpy.advanced_streaming import LBMPeriodicityHandling
 from lbmpy.boundaries import NoSlip, LatticeBoltzmannBoundaryHandling
 ```
 
 %% Cell type:markdown id: tags:
 
-If `pycuda` is installed the simulation automatically runs on GPU
+If `cupy` is installed the simulation automatically runs on GPU
 
 %% Cell type:code id: tags:
 
 ``` python
 try:
-    import pycuda
+    import cupy
 except ImportError:
-    pycuda = None
+    cupy = None
     gpu = False
     target = ps.Target.CPU
-    print('No pycuda installed')
+    print('No cupy installed')
 
-if pycuda:
+if cupy:
     gpu = True
     target = ps.Target.GPU
 ```
 
-%% Output
-
-    No pycuda installed
-
 %% Cell type:markdown id: tags:
 
 # Capillary wave simulated with a phase-field model for immiscible fluids
 
 %% Cell type:markdown id: tags:
 
 ## Geometry Setup
 
 First of all the stencils for the phase-field LB step as well as the stencil for the hydrodynamic LB step are defined. According to the stencils the simulation runs either in 2D or 3D
 
 %% Cell type:code id: tags:
 
 ``` python
 stencil_phase = LBStencil(Stencil.D2Q9)
 stencil_hydro = LBStencil(Stencil.D2Q9)
 assert(stencil_hydro.D == stencil_phase.D)
 
 dimensions = stencil_phase.D
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # user defined input
 
 Re = 10 # Reynolds number
 domain_width = 50
 fluid_depth = 0.5 * domain_width
 amplitude = 0.01 * domain_width
 relaxation_rate_heavy = 1.99
 mobility = 0.02 # phase field mobility
 interface_width = 5 # phase field interface width
 density_heavy = 1.0 # density of heavy phase
 density_ratio = 1000
 density_light = density_heavy / density_ratio # density of light phase
 kinematic_viscosity_ratio = 1
 
 kinematic_viscosity_heavy = 1 / 3 * (1 / relaxation_rate_heavy - 0.5)
 kinematic_viscosity_light = kinematic_viscosity_heavy / kinematic_viscosity_ratio
 wavelength = domain_width
 wavenumber = 2.0 * np.pi / domain_width
 wave_frequency = Re * kinematic_viscosity_heavy / domain_width / amplitude # angular wave frequency
 surface_tension = wave_frequency**2 * (density_heavy + density_light) / wavenumber**3
 gravitational_acceleration = 0
 Pe = domain_width * amplitude * wave_frequency / mobility # Peclet number
 Cn = interface_width / domain_width # Cahn number
 dynamic_viscosity_heavy = kinematic_viscosity_heavy * density_heavy
 relaxation_time_heavy = 3.0 * kinematic_viscosity_heavy
 kinematic_viscosity_light = kinematic_viscosity_heavy / kinematic_viscosity_ratio
 dynamic_viscosity_light = kinematic_viscosity_light * density_light
 relaxation_time_light = 3.0 * kinematic_viscosity_light
 
 timesteps = int(80 / wave_frequency)
 
 data_extract_frequency = int(0.1 / wave_frequency)
 vtk_output_frequency = int(1 / wave_frequency)
 vtk_output_path = "vtk_out/capillary-wave"
 vtk_base_directory = vtk_output_path.split("/")[0] # create directory for vtk-output if it does not yet exist
 if not os.path.exists(vtk_base_directory):
     os.mkdir(os.getcwd() + "/" + vtk_base_directory)
 
 domain_size = (domain_width, domain_width)
 filename = "pf-re-" + str(Re) + "-resolution-" + str(domain_width) + ".txt"
 
 print("timesteps =", timesteps)
 print("Re =", Re)
 print("Pe =", Pe)
 print("Cn =", Cn)
 print("domain_width =", domain_width)
 print("fluid_depth =", fluid_depth)
 print("amplitude =", amplitude)
 print("relaxation_rate_heavy =", relaxation_rate_heavy)
 print("mobility =", mobility)
 print("interface_width =", interface_width)
 print("density_heavy =", density_heavy)
 print("density_light =", density_light)
 print("density_ratio =", density_ratio)
 print("kinematic_viscosity_heavy =", kinematic_viscosity_heavy)
 print("kinematic_viscosity_light =", kinematic_viscosity_light)
 print("kinematic_viscosity_ratio =", kinematic_viscosity_ratio)
 print("dynamic_viscosity_heavy =", dynamic_viscosity_heavy)
 print("dynamic_viscosity_light =", dynamic_viscosity_light)
 print("dynamic_viscosity_ratio =", dynamic_viscosity_heavy/dynamic_viscosity_light)
 print("wavelength =", wavelength)
 print("wavenumber =", wavenumber)
 print("wave_frequency =", wave_frequency)
 print("surface_tension =", surface_tension)
 print("gravitational_acceleration =", gravitational_acceleration)
 ```
 
 %% Output
 
     timesteps = 238800
     Re = 10
     Pe = 0.41876046901171776
     Cn = 0.1
     domain_width = 50
     fluid_depth = 25.0
     amplitude = 0.5
     relaxation_rate_heavy = 1.99
     mobility = 0.02
     interface_width = 5
     density_heavy = 1.0
     density_light = 0.001
     density_ratio = 1000
     kinematic_viscosity_heavy = 0.0008375209380234356
     kinematic_viscosity_light = 0.0008375209380234356
     kinematic_viscosity_ratio = 1
     dynamic_viscosity_heavy = 0.0008375209380234356
     dynamic_viscosity_light = 8.375209380234357e-07
     dynamic_viscosity_ratio = 1000.0
     wavelength = 50
     wavenumber = 0.12566370614359174
     wave_frequency = 0.00033500837520937423
     surface_tension = 5.6612953740701554e-05
     gravitational_acceleration = 0
 
 %% Cell type:code id: tags:
 
 ``` python
 parameters = AllenCahnParameters(density_heavy=density_heavy, density_light=density_light,
                                  dynamic_viscosity_heavy=dynamic_viscosity_heavy,
                                  dynamic_viscosity_light=dynamic_viscosity_light,
                                  surface_tension=surface_tension, mobility=mobility,
                                  interface_thickness=interface_width)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 parameters
 ```
 
 %% Output
 
     <lbmpy.phasefield_allen_cahn.parameter_calculation.AllenCahnParameters at 0x146131580>
 
 %% Cell type:markdown id: tags:
 
 ## Fields
 
 As a next step all fields which are needed get defined. To do so we create a `datahandling` object. More details about it can be found in the third tutorial of the [pystencils framework]( http://pycodegen.pages.walberla.net/pystencils/). Basically it holds all fields and manages the kernel runs.
 
 %% Cell type:code id: tags:
 
 ``` python
 # create a datahandling object
 dh = ps.create_data_handling((domain_size), periodicity=(True, False), parallel=False, default_target=target)
 
 # fields
 g = dh.add_array("g", values_per_cell=len(stencil_hydro))
 dh.fill("g", 0.0, ghost_layers=True)
 h = dh.add_array("h",values_per_cell=len(stencil_phase))
 dh.fill("h", 0.0, ghost_layers=True)
 
 g_tmp = dh.add_array("g_tmp", values_per_cell=len(stencil_hydro))
 dh.fill("g_tmp", 0.0, ghost_layers=True)
 h_tmp = dh.add_array("h_tmp",values_per_cell=len(stencil_phase))
 dh.fill("h_tmp", 0.0, ghost_layers=True)
 
 u = dh.add_array("u", values_per_cell=dh.dim)
 dh.fill("u", 0.0, ghost_layers=True)
 
 rho = dh.add_array("rho", values_per_cell=1)
 dh.fill("rho", 1.0, ghost_layers=True)
 
 C = dh.add_array("C")
 dh.fill("C", 0.0, ghost_layers=True)
 C_tmp = dh.add_array("C_tmp")
 dh.fill("C_tmp", 0.0, ghost_layers=True)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # set the frequency of the file outputs
 vtk_writer = dh.create_vtk_writer(vtk_output_path, ["C", "u", "rho"], ghost_layers=True)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 rho_L = parameters.symbolic_density_light
 rho_H = parameters.symbolic_density_heavy
 density = rho_L + C.center * (rho_H - rho_L)
 
 body_force = [0, 0, 0]
 body_force[1] = parameters.symbolic_gravitational_acceleration * density
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Definition of the lattice Boltzmann methods
 
 %% Cell type:code id: tags:
 
 ``` python
 w_c = parameters.symbolic_omega_phi
 
 config_phase = LBMConfig(stencil=stencil_phase, method=Method.MRT, compressible=True,
                          delta_equilibrium=False,
                          force=sp.symbols("F_:2"), velocity_input=u,
                          weighted=True, relaxation_rates=[0, w_c, w_c, 1, 1, 1, 1, 1, 1],
                          output={'density': C_tmp}, kernel_type='stream_pull_collide')
 
 method_phase = create_lb_method(lbm_config=config_phase)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 omega = parameters.omega(C)
 config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.MRT, compressible=False,
                          weighted=True, relaxation_rates=[omega, 1, 1, 1],
                          force=sp.symbols("F_:2"),
                          output={'velocity': u}, kernel_type='collide_stream_push')
 
 method_hydro = create_lb_method(lbm_config=config_hydro)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Initialization
 
 %% Cell type:code id: tags:
 
 ``` python
 h_updates = initializer_kernel_phase_field_lb(method_phase, C, u, h, parameters)
 g_updates = initializer_kernel_hydro_lb(method_hydro, 1, u, g)
 
 h_init = ps.create_kernel(h_updates, target=dh.default_target, cpu_openmp=True).compile()
 g_init = ps.create_kernel(g_updates, target=dh.default_target, cpu_openmp=True).compile()
 ```
 
 %% Output
 
     [['spatialInner0'], ['spatialInner1']]
     [['spatialInner0'], ['spatialInner1']]
 
 %% Cell type:code id: tags:
 
 ``` python
 # initialize the phase field
 def initialize_phasefield():
     Nx = domain_size[0]
     Ny = domain_size[1]
 
     for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
         # get x as cell center coordinate, i.e., including shift with 0.5
         x = np.zeros_like(block.midpoint_arrays[0])
         x[:, :] = block.midpoint_arrays[0]
 
         # get y as cell center coordinate, i.e., including shift with 0.5
         y = np.zeros_like(block.midpoint_arrays[1])
         y[:, :] = block.midpoint_arrays[1]
 
         tmp = fluid_depth + amplitude * np.cos(x / (Nx - 1) * np.pi * 2 + np.pi)
 
         # initialize diffuse interface with tanh profile
         init_values = 0.5 - 0.5 * np.tanh((y - tmp) / (interface_width / 2))
 
         block["C"][:, :] = init_values
         block["C_tmp"][:, :] = init_values
 
     if gpu:
         dh.all_to_gpu()
 
     dh.run_kernel(h_init, **parameters.symbolic_to_numeric_map)
     dh.run_kernel(g_init)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # plot initial profile
 x = np.arange(0, domain_size[0], 1)   # start,stop,step
 y = fluid_depth + amplitude * np.cos(x / (domain_size[0] - 1) * np.pi * 2 + np.pi)
 
 plt.plot(x,y, marker='o')
 plt.show()
 ```
 
 %% Output
 
 
 %% Cell type:code id: tags:
 
 ``` python
 # check symmetry of sine-curve
 for i in range(0, domain_size[0] - 1):
     if (i >= domain_size[0] * 0.5):
         continue
     if (y[i] != y[domain_size[0] - 1 - i]):
         print("asymmetric:", y[i], y[domain_size[0] - 1 - i])
 ```
 
 %% Output
 
     asymmetric: 25.111260466978155 25.11126046697816
 
 %% Cell type:code id: tags:
 
 ``` python
 initialize_phasefield()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 plt.scalar_field(dh.cpu_arrays["C"])
 plt.show()
 ```
 
 %% Output
 
 
 %% Cell type:code id: tags:
 
 ``` python
 # plot initial interface profile
 x = np.arange(0, domain_size[1]+2, 1)   # start,stop,step
 y = dh.cpu_arrays["C"][int(domain_size[0] * 0.5), :]
 
 plt.plot(x,y, marker='o')
 plt.grid()
 plt.show()
 ```
 
 %% Output
 
 
 %% Cell type:code id: tags:
 
 ``` python
 force_h = interface_tracking_force(C, stencil_phase, parameters)
 hydro_force = hydrodynamic_force(g, C, method_hydro, parameters, body_force)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Definition of the LB update rules
 
 %% Cell type:code id: tags:
 
 ``` python
 lbm_optimisation = LBMOptimisation(symbolic_field=h, symbolic_temporary_field=h_tmp)
 allen_cahn_update_rule = create_lb_update_rule(lbm_config=config_phase,
                                                lbm_optimisation=lbm_optimisation)
 
 allen_cahn_update_rule = add_interface_tracking_force(allen_cahn_update_rule, force_h)
 
 ast_kernel = ps.create_kernel(allen_cahn_update_rule, target=dh.default_target, cpu_openmp=True)
 kernel_allen_cahn_lb = ast_kernel.compile()
 ```
 
 %% Output
 
     [['spatialInner0'], ['spatialInner1']]
 
 %% Cell type:code id: tags:
 
 ``` python
 force_Assignments = hydrodynamic_force_assignments(g, u, C, method_hydro, parameters, body_force)
 
 lbm_optimisation = LBMOptimisation(symbolic_field=g, symbolic_temporary_field=g_tmp)
 hydro_lb_update_rule = create_lb_update_rule(lbm_config=config_hydro,
                                              lbm_optimisation=lbm_optimisation)
 
 hydro_lb_update_rule = add_hydrodynamic_force(hydro_lb_update_rule, force_Assignments, C, g, parameters)
 
 ast_kernel = ps.create_kernel(hydro_lb_update_rule, target=dh.default_target, cpu_openmp=True)
 kernel_hydro_lb = ast_kernel.compile()
 ```
 
 %% Output
 
     [['spatialInner0'], ['spatialInner1']]
 
 %% Cell type:markdown id: tags:
 
 ## Boundary Conditions
 
 %% Cell type:code id: tags:
 
 ``` python
 # periodic Boundarys for g, h and C
 periodic_BC_C = dh.synchronization_function(C.name, target=dh.default_target, optimization = {"openmp": True})
 
 periodic_BC_g = LBMPeriodicityHandling(stencil=stencil_hydro, data_handling=dh, pdf_field_name=g.name,
                                        streaming_pattern='push')
 periodic_BC_h = LBMPeriodicityHandling(stencil=stencil_phase, data_handling=dh, pdf_field_name=h.name,
                                        streaming_pattern='pull')
 
 # No slip boundary for the phasefield lbm
 bh_allen_cahn = LatticeBoltzmannBoundaryHandling(method_phase, dh, 'h',
                                                  target=dh.default_target, name='boundary_handling_h',
                                                  streaming_pattern='pull')
 
 # No slip boundary for the velocityfield lbm
 bh_hydro = LatticeBoltzmannBoundaryHandling(method_hydro, dh, 'g' ,
                                             target=dh.default_target, name='boundary_handling_g',
                                             streaming_pattern='push')
 
 contact_angle = BoundaryHandling(dh, C.name, stencil_hydro, target=dh.default_target)
 contact = ContactAngle(90, interface_width)
 
 wall = NoSlip()
 if dimensions == 2:
     bh_allen_cahn.set_boundary(wall, make_slice[:, 0])
     bh_allen_cahn.set_boundary(wall, make_slice[:, -1])
 
     bh_hydro.set_boundary(wall, make_slice[:, 0])
     bh_hydro.set_boundary(wall, make_slice[:, -1])
 
     contact_angle.set_boundary(contact, make_slice[:, 0])
     contact_angle.set_boundary(contact, make_slice[:, -1])
 else:
     bh_allen_cahn.set_boundary(wall, make_slice[:, 0, :])
     bh_allen_cahn.set_boundary(wall, make_slice[:, -1, :])
 
     bh_hydro.set_boundary(wall, make_slice[:, 0, :])
     bh_hydro.set_boundary(wall, make_slice[:, -1, :])
 
     contact_angle.set_boundary(contact, make_slice[:, 0, :])
     contact_angle.set_boundary(contact, make_slice[:, -1, :])
 
 
 bh_allen_cahn.prepare()
 bh_hydro.prepare()
 contact_angle.prepare()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Full timestep
 
 %% Cell type:code id: tags:
 
 ``` python
 # definition of the timestep for the immiscible fluids model
 def timeloop():
     # Solve the interface tracking LB step with boundary conditions
     periodic_BC_h()
     bh_allen_cahn()
     dh.run_kernel(kernel_allen_cahn_lb, **parameters.symbolic_to_numeric_map)
     dh.swap("C", "C_tmp")
 
     # apply the three phase-phase contact angle
     contact_angle()
     # periodic BC of the phase-field
     periodic_BC_C()
 
     # solve the hydro LB step with boundary conditions
     dh.run_kernel(kernel_hydro_lb, **parameters.symbolic_to_numeric_map)
     periodic_BC_g()
     bh_hydro()
 
     # compute density (only for vtk output)
     # must be done BEFORE swapping fields to avoid having outdated values
     compute_density()
 
     # field swaps
     dh.swap("h", "h_tmp")
     dh.swap("g", "g_tmp")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def initialize_hydrostatic_pressure():
     for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
 
         # get y as cell center coordinate, i.e., including shift with 0.5
         y = np.zeros_like(block.midpoint_arrays[1])
         y[:, :] = block.midpoint_arrays[1]
 
         # compute hydrostatic density
         rho_hydrostatic = 1 + 3 * gravitational_acceleration * (y - fluid_depth)
 
         # subtract 1 because PDFs in incompressible LBM are centered around zero in lbmpy
         rho_hydrostatic -= 1
 
         # set equilibrium PDFs with velocity=0 and rho;
         for i in range(0, stencil_hydro.Q, 1):
             block["g"][:,:,i] = method_hydro.weights[i] * rho_hydrostatic[:,:]
             block["g_tmp"][:,:,i] = method_hydro.weights[i] * rho_hydrostatic[:,:]
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def compute_density():
     for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
         # PDFs in incompressible LBM are centered around zero in lbmpy
         # => add 1 to whole sum, i.e., initialize with 1
         block["rho"].fill(1);
 
         # compute density
         for i in range(block["g"].shape[-1]):
             block["rho"][:,:] += block["g"][:,:,i]
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if (gravitational_acceleration != 0):
     initialize_hydrostatic_pressure()
     compute_density()
     plt.scalar_field(dh.cpu_arrays["rho"])
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 print("================================= start of the simulation ===================================")
 
 start = time.time()
 
 pbar = tqdm(total=timesteps)
 
 timestep = []
 surface_position = []
 symmetry_norm = []
 
 for i in range(0, timesteps):
 
     sum_c_2 = 0.0
     sum_delta_c_2 = 0.0
 
     # write vtk output
     if(i % vtk_output_frequency == 0):
         if gpu:
             dh.to_cpu("C")
             dh.to_cpu("u")
             dh.to_cpu("rho")
         vtk_writer(i)
 
     # extract data (to be written to file)
     if(i % data_extract_frequency == 0):
         if gpu:
             dh.to_cpu("C")
             dh.to_cpu("u")
             dh.to_cpu("rho")
 
         pbar.update(data_extract_frequency)
 
         timestep.append(i)
 
         ny = domain_size[1]
 
         # get index containing phase field value < 0.5
         i1 = np.argmax(dh.cpu_arrays["C"][int(domain_size[0] * 0.5), :] < 0.5)
         i0 = i1 - 1 # index containing phase field value >= 0.5
 
         f0 = dh.cpu_arrays["C"][int(domain_size[0] * 0.5), i0] # phase field value >= 0.5
         f1 = dh.cpu_arrays["C"][int(domain_size[0] * 0.5), i1] # phase field value < 0.5
 
         # coordinate of cell center is index+0.5-1 (0.5 to get to cell center; -1 to remove ghost layer from index)
         y0 = i0 + 0.5 - 1
         y1 = i1 + 0.5 - 1
 
         #interpolate
         surface_position.append( y0 + (y1 - y0) / (f1 - f0) * (0.5 - f0) )
 
         # evaluate symmetry in x-direction
         for y in range(0, domain_size[1] - 1):
             for x in range(0, domain_size[0] - 1):
                 if (x >= domain_size[0] * 0.5):
                     continue
 
                 x_mirrored = domain_size[0] - 1 - x;
                 sum_c_2 += dh.cpu_arrays["C"][x, y]**2
                 sum_delta_c_2 += (dh.cpu_arrays["C"][x, y] - dh.cpu_arrays["C"][x_mirrored, y])**2
 
         symmetry_norm.append( (sum_delta_c_2 / sum_c_2)**2 )
 
     timeloop()
 
 pbar.close()
 end = time.time()
 sim_time = end - start
 print("\n")
 print("time needed for the calculation: %4.4f" % sim_time , "seconds")
 
 nrOfCells = np.prod(domain_size)
 mlups = nrOfCells * timesteps / sim_time * 1e-6
 
 print("MLUPS: %4.2f" % mlups)
 ```
 
 %% Output
 
     ================================= start of the simulation ===================================
 
 
     
     
     time needed for the calculation: 92.7764 seconds
     MLUPS: 6.43
 
 %% Cell type:code id: tags:
 
 ``` python
 if gpu:
     dh.to_cpu("C")
     dh.to_cpu("u")
     dh.to_cpu("rho")
 
 plt.scalar_field(dh.cpu_arrays["C"])
 ```
 
 %% Output
 
     <matplotlib.image.AxesImage at 0x146e329a0>
 
 
 %% Cell type:code id: tags:
 
 ``` python
 # non-dimensionalize time and surface position
 t_nd = [value * wave_frequency for value in timestep]
 a_nd = [(value - fluid_depth) / amplitude for value in surface_position]
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # plot surface position over time
 plt.xlabel('$t/-$', fontsize=18)
 plt.ylabel('$a/-$', fontsize=18)
 plt.grid(color='black', linestyle='-', linewidth=0.1)
 plt.plot(t_nd, a_nd, color=(0.121, 0.231, 0.4))
 plt.show()
 ```
 
 %% Output
 
 
 %% Cell type:code id: tags:
 
 ``` python
 # plot analytical solutions
 
 # analytical model from Prosperetti, Motion of two superposed viscous fluids, 1981, doi: 10.1063/1.863522
 nu = kinematic_viscosity_heavy
 k = wavenumber
 omega_0 = wave_frequency
 rho_u = density_heavy
 rho_l = density_light
 a0 = amplitude
 beta = rho_l * rho_u / (rho_l + rho_u)**2
 
 k2nu = k**2 * nu # helper variable
 
 # polynomials from equation (21)
 p0 =  (1 - 4 * beta) * k2nu**2 + omega_0**2
 p1 =  4 * (1 - 3 * beta) * k2nu**1.5
 p2 =  2 * (1 - 6 * beta) * k2nu
 p3 = -4 * beta * k2nu**(0.5)
 p4 =  1
 p = [p0, p1, p2, p3, p4]
 
 # get roots of equation (21)
 z = np.polynomial.polynomial.polyroots(p)
 
 Z = np.array([(z[1]  - z[0]) * (z[2] - z[0]) * (z[3] - z[0]),
               (z[0]  - z[1]) * (z[2] - z[1]) * (z[3] - z[1]),
               (z[0]  - z[2]) * (z[1] - z[2]) * (z[3] - z[2]),
               (z[0]  - z[3]) * (z[1] - z[3]) * (z[2] - z[3])])
 
 t = np.arange(0, timesteps, 1)
 
 # equation (20)
 tmp0 = (1 - 4 * beta) * k2nu**2 # helper variable
 term0 = 4 * tmp0 / (8 * tmp0 + omega_0**2) * a0 * erfc((k2nu * t)**0.5) # first term in equation (20)
 
 term1 = 0.0
 for i in range(0, 4, 1):
     tmp1 = z[i]**2 - k2nu # helper variable
     term1 += z[i] / Z[i] * omega_0**2 * a0 / tmp1 * np.exp(tmp1 * t) * erfc(z[i] * t**0.5)
 
 a = np.real(term0 + term1)
 
 # non-dimesionalize time and surface position
 pos_anal_lin_org = a / a0
 t_anal_org = t * omega_0
 
 # non-dimesionalize time and surface position
 a_anal_nd = a / amplitude
 t_anal_nd = t * wave_frequency
 
 plt.plot(t_anal_nd, a_anal_nd, label="analytical solution")
 plt.plot(t_nd, a_nd, label="simulation")
 
 # correction factor from Denner et al., Dispersion and viscous attenuation of capillary waves with finite amplitude, 2016, doi: 10.1140/epjst/e2016-60199-2
 # model is now valid for initial amplitudes of about 0.1*wavelength
 a0_hat = a0 / wavelength
 C = -4.5 * a0_hat**3 + 5.3 * a0_hat**2 + 0.18 * a0_hat + 1
 a_anal_nd_corr = a_anal_nd * C
 t_anal_nd_corr = t_anal_nd * C
 plt.plot(t_anal_nd_corr, a_anal_nd_corr, label="linear analytical solution with correction factor")
 
 plt.xlabel('$t/-$', fontsize=18)
 plt.ylabel('$a/-$', fontsize=18)
 plt.legend()
 plt.grid()
 plt.show()
 ```
 
 %% Output
 
 
 %% Cell type:code id: tags:
 
 ``` python
 # plot symmetry norm over time
 plt.xlabel('$t/-$', fontsize=18)
 plt.ylabel('$symmetry-norm/-$', fontsize=18)
 plt.grid(color='black', linestyle='-', linewidth=0.1)
 plt.plot(t_nd, symmetry_norm, color=(0.121, 0.231, 0.4))
 plt.show()
 ```
 
 %% Output
 
 
 %% Cell type:code id: tags:
 
 ``` python
 # store results to file
 import csv
 with open(filename, 'w') as f:
     writer = csv.writer(f, delimiter='\t')
     writer.writerows(zip(t_nd, a_nd, symmetry_norm))
 ```

lbmpy_tests/full_scenarios/phasefield_allen_cahn/phasefield-gravity-wave.ipynb → tests/full_scenarios/phasefield_allen_cahn/phasefield-gravity-wave.ipynb

 %% Cell type:code id: tags:
 
 ``` python
 import os
 
 import time
 from tqdm.notebook import tqdm
 
 import numpy as np
 import math
 
 from pystencils.session import *
 from lbmpy.session import *
 
 from pystencils.simp import sympy_cse
 from pystencils.boundaries import BoundaryHandling
 
 from lbmpy.phasefield_allen_cahn.parameter_calculation import AllenCahnParameters
 from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle
 from lbmpy.phasefield_allen_cahn.kernel_equations import *
 
 from lbmpy.advanced_streaming import LBMPeriodicityHandling
 from lbmpy.boundaries import NoSlip, LatticeBoltzmannBoundaryHandling
 ```
 
 %% Cell type:markdown id: tags:
 
-If `pycuda` is installed the simulation automatically runs on GPU
+If `cupy` is installed the simulation automatically runs on GPU
 
 %% Cell type:code id: tags:
 
 ``` python
 try:
-    import pycuda
+    import cupy
 except ImportError:
-    pycuda = None
+    cupy = None
     gpu = False
     target = ps.Target.CPU
-    print('No pycuda installed')
+    print('No cupy installed')
 
-if pycuda:
+if cupy:
     gpu = True
     target = ps.Target.GPU
 ```
 
-%% Output
-
-    No pycuda installed
-
 %% Cell type:markdown id: tags:
 
 # Gravity wave simulated with a phase-field model for immiscible fluids
 
 %% Cell type:markdown id: tags:
 
 ## Geometry Setup
 
 First of all the stencils for the phase-field LB step as well as the stencil for the hydrodynamic LB step are defined. According to the stencils the simulation runs either in 2D or 3D
 
 %% Cell type:code id: tags:
 
 ``` python
 stencil_phase = LBStencil(Stencil.D2Q9)
 stencil_hydro = LBStencil(Stencil.D2Q9)
 assert(stencil_hydro.D == stencil_phase.D)
 
 dimensions = stencil_phase.D
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # user defined input
 
 Re = 10 # Reynolds number
 domain_width = 50
 fluid_depth = 0.5 * domain_width
 amplitude = 0.01 * domain_width
 relaxation_rate_heavy = 1.99
 mobility = 0.02 # phase field mobility
 interface_width = 5 # phase field interface width
 density_heavy = 1.0 # density of heavy phase
 density_ratio = 1000
 density_light = density_heavy / density_ratio # density of light phase
 kinematic_viscosity_ratio = 1
 
 kinematic_viscosity_heavy = 1 / 3 * (1 / relaxation_rate_heavy - 0.5)
 kinematic_viscosity_light = kinematic_viscosity_heavy / kinematic_viscosity_ratio
 wavelength = domain_width
 wavenumber = 2.0 * np.pi / domain_width
 wave_frequency = Re * kinematic_viscosity_heavy / domain_width / amplitude # angular wave frequency
 surface_tension = 0
 gravitational_acceleration = - wave_frequency**2 / wavenumber / np.tanh(wavenumber * fluid_depth)
 Pe = domain_width * amplitude * wave_frequency / mobility # Peclet number
 Cn = interface_width / domain_width # Cahn number
 dynamic_viscosity_heavy = kinematic_viscosity_heavy * density_heavy
 relaxation_time_heavy = 3.0 * kinematic_viscosity_heavy
 kinematic_viscosity_light = kinematic_viscosity_heavy / kinematic_viscosity_ratio
 dynamic_viscosity_light = kinematic_viscosity_light * density_light
 relaxation_time_light = 3.0 * kinematic_viscosity_light
 
 timesteps = int(80 / wave_frequency)
 
 data_extract_frequency = int(0.1 / wave_frequency)
 vtk_output_frequency = int(1 / wave_frequency)
 vtk_output_path = "vtk_out/gravity-wave"
 vtk_base_directory = vtk_output_path.split("/")[0] # create directory for vtk-output if it does not yet exist
 if not os.path.exists(vtk_base_directory):
     os.mkdir(os.getcwd() + "/" + vtk_base_directory)
 
 domain_size = (domain_width, domain_width)
 filename = "pf-re-" + str(Re) + "-resolution-" + str(domain_width) + ".txt"
 
 print("timesteps =", timesteps)
 print("Re =", Re)
 print("Pe =", Pe)
 print("Cn =", Cn)
 print("domain_width =", domain_width)
 print("fluid_depth =", fluid_depth)
 print("amplitude =", amplitude)
 print("relaxation_rate_heavy =", relaxation_rate_heavy)
 print("mobility =", mobility)
 print("interface_width =", interface_width)
 print("density_heavy =", density_heavy)
 print("density_light =", density_light)
 print("density_ratio =", density_ratio)
 print("kinematic_viscosity_heavy =", kinematic_viscosity_heavy)
 print("kinematic_viscosity_light =", kinematic_viscosity_light)
 print("kinematic_viscosity_ratio =", kinematic_viscosity_ratio)
 print("dynamic_viscosity_heavy =", dynamic_viscosity_heavy)
 print("dynamic_viscosity_light =", dynamic_viscosity_light)
 print("dynamic_viscosity_ratio =", dynamic_viscosity_heavy/dynamic_viscosity_light)
 print("wavelength =", wavelength)
 print("wavenumber =", wavenumber)
 print("wave_frequency =", wave_frequency)
 print("surface_tension =", surface_tension)
 print("gravitational_acceleration =", gravitational_acceleration)
 print("data_extract_frequency =", data_extract_frequency)
 print("vtk_output_frequency =", vtk_output_frequency)
 ```
 
 %% Output
 
     timesteps = 238800
     Re = 10
     Pe = 0.41876046901171776
     Cn = 0.1
     domain_width = 50
     fluid_depth = 25.0
     amplitude = 0.5
     relaxation_rate_heavy = 1.99
     mobility = 0.02
     interface_width = 5
     density_heavy = 1.0
     density_light = 0.001
     density_ratio = 1000
     kinematic_viscosity_heavy = 0.0008375209380234356
     kinematic_viscosity_light = 0.0008375209380234356
     kinematic_viscosity_ratio = 1
     dynamic_viscosity_heavy = 0.0008375209380234356
     dynamic_viscosity_light = 8.375209380234357e-07
     dynamic_viscosity_ratio = 1000.0
     wavelength = 50
     wavenumber = 0.12566370614359174
     wave_frequency = 0.00033500837520937423
     surface_tension = 0
     gravitational_acceleration = -8.964447065459421e-07
     data_extract_frequency = 298
     vtk_output_frequency = 2985
 
 %% Cell type:code id: tags:
 
 ``` python
 parameters = AllenCahnParameters(density_heavy=density_heavy, density_light=density_light,
                                  dynamic_viscosity_heavy=dynamic_viscosity_heavy,
                                  dynamic_viscosity_light=dynamic_viscosity_light,
                                  gravitational_acceleration=gravitational_acceleration,
                                  surface_tension=surface_tension, mobility=mobility,
                                  interface_thickness=interface_width)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 parameters
 ```
 
 %% Output
 
     <lbmpy.phasefield_allen_cahn.parameter_calculation.AllenCahnParameters at 0x12d1bd550>
 
 %% Cell type:markdown id: tags:
 
 ## Fields
 
 As a next step all fields which are needed get defined. To do so we create a `datahandling` object. More details about it can be found in the third tutorial of the [pystencils framework]( http://pycodegen.pages.walberla.net/pystencils/). Basically it holds all fields and manages the kernel runs.
 
 %% Cell type:code id: tags:
 
 ``` python
 # create a datahandling object
 dh = ps.create_data_handling((domain_size), periodicity=(True, False), parallel=False, default_target=target)
 
 # fields
 g = dh.add_array("g", values_per_cell=len(stencil_hydro))
 dh.fill("g", 0.0, ghost_layers=True)
 h = dh.add_array("h",values_per_cell=len(stencil_phase))
 dh.fill("h", 0.0, ghost_layers=True)
 
 g_tmp = dh.add_array("g_tmp", values_per_cell=len(stencil_hydro))
 dh.fill("g_tmp", 0.0, ghost_layers=True)
 h_tmp = dh.add_array("h_tmp",values_per_cell=len(stencil_phase))
 dh.fill("h_tmp", 0.0, ghost_layers=True)
 
 u = dh.add_array("u", values_per_cell=dh.dim)
 dh.fill("u", 0.0, ghost_layers=True)
 
 rho = dh.add_array("rho", values_per_cell=1)
 dh.fill("rho", 1.0, ghost_layers=True)
 
 C = dh.add_array("C")
 dh.fill("C", 0.0, ghost_layers=True)
 C_tmp = dh.add_array("C_tmp")
 dh.fill("C_tmp", 0.0, ghost_layers=True)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # set the frequency of the file outputs
 vtk_writer = dh.create_vtk_writer(vtk_output_path, ["C", "u", "rho"], ghost_layers=True)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 rho_L = parameters.symbolic_density_light
 rho_H = parameters.symbolic_density_heavy
 density = rho_L + C.center * (rho_H - rho_L)
 
 body_force = [0, 0, 0]
 body_force[1] = parameters.symbolic_gravitational_acceleration * density
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Definition of the lattice Boltzmann methods
 
 %% Cell type:code id: tags:
 
 ``` python
 w_c = parameters.symbolic_omega_phi
 
 config_phase = LBMConfig(stencil=stencil_phase, method=Method.MRT, compressible=True,
                          delta_equilibrium=False,
                          force=sp.symbols("F_:2"), velocity_input=u,
                          weighted=True, relaxation_rates=[0, w_c, w_c, 1, 1, 1, 1, 1, 1],
                          output={'density': C_tmp}, kernel_type='stream_pull_collide')
 
 method_phase = create_lb_method(lbm_config=config_phase)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 omega = parameters.omega(C)
 config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.MRT, compressible=False,
                          weighted=True, relaxation_rates=[omega, 1, 1, 1],
                          force=sp.symbols("F_:2"),
                          output={'velocity': u}, kernel_type='collide_stream_push')
 
 method_hydro = create_lb_method(lbm_config=config_hydro)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 method_hydro
 ```
 
 %% Output
 
     <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x12d4a1e50>
 
 %% Cell type:markdown id: tags:
 
 ## Initialization
 
 %% Cell type:code id: tags:
 
 ``` python
 h_updates = initializer_kernel_phase_field_lb(method_phase, C, u, h, parameters)
 g_updates = initializer_kernel_hydro_lb(method_hydro, 1, u, g)
 
 h_init = ps.create_kernel(h_updates, target=dh.default_target, cpu_openmp=True).compile()
 g_init = ps.create_kernel(g_updates, target=dh.default_target, cpu_openmp=True).compile()
 ```
 
 %% Output
 
     [['spatialInner0'], ['spatialInner1']]
     [['spatialInner0'], ['spatialInner1']]
 
 %% Cell type:code id: tags:
 
 ``` python
 # initialize the phase field
 def initialize_phasefield():
     Nx = domain_size[0]
     Ny = domain_size[1]
 
     for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
         # get x as cell center coordinate, i.e., including shift with 0.5
         x = np.zeros_like(block.midpoint_arrays[0])
         x[:, :] = block.midpoint_arrays[0]
 
         # get y as cell center coordinate, i.e., including shift with 0.5
         y = np.zeros_like(block.midpoint_arrays[1])
         y[:, :] = block.midpoint_arrays[1]
 
         tmp = fluid_depth + amplitude * np.cos(x / (Nx - 1) * np.pi * 2 + np.pi)
 
         # initialize diffuse interface with tanh profile
         init_values = 0.5 - 0.5 * np.tanh((y - tmp) / (interface_width / 2))
 
         block["C"][:, :] = init_values
         block["C_tmp"][:, :] = init_values
 
     if gpu:
         dh.all_to_gpu()
 
     dh.run_kernel(h_init, **parameters.symbolic_to_numeric_map)
     dh.run_kernel(g_init)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # plot initial profile
 x = np.arange(0, domain_size[0], 1)   # start,stop,step
 y = fluid_depth + amplitude * np.cos(x / (domain_size[0] - 1) * np.pi * 2 + np.pi)
 
 plt.plot(x,y, marker='o')
 plt.show()
 ```
 
 %% Output
 
 
 %% Cell type:code id: tags:
 
 ``` python
 initialize_phasefield()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 plt.scalar_field(dh.cpu_arrays["C"])
 plt.show()
 ```
 
 %% Output
 
 
 %% Cell type:code id: tags:
 
 ``` python
 # plot initial interface profile
 x = np.arange(0, domain_size[1]+2, 1)   # start,stop,step
 y = dh.cpu_arrays["C"][int(domain_size[0] * 0.5), :]
 
 plt.plot(x,y, marker='o')
 plt.grid()
 plt.show()
 ```
 
 %% Output
 
 
 %% Cell type:code id: tags:
 
 ``` python
 force_h = interface_tracking_force(C, stencil_phase, parameters)
 hydro_force = hydrodynamic_force(g, C, method_hydro, parameters, body_force)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Definition of the LB update rules
 
 %% Cell type:code id: tags:
 
 ``` python
 lbm_optimisation = LBMOptimisation(symbolic_field=h, symbolic_temporary_field=h_tmp)
 allen_cahn_update_rule = create_lb_update_rule(lbm_config=config_phase,
                                                lbm_optimisation=lbm_optimisation)
 
 allen_cahn_update_rule = add_interface_tracking_force(allen_cahn_update_rule, force_h)
 
 ast_kernel = ps.create_kernel(allen_cahn_update_rule, target=dh.default_target, cpu_openmp=True)
 kernel_allen_cahn_lb = ast_kernel.compile()
 ```
 
 %% Output
 
     [['spatialInner0'], ['spatialInner1']]
 
 %% Cell type:code id: tags:
 
 ``` python
 force_Assignments = hydrodynamic_force_assignments(g, u, C, method_hydro, parameters, body_force)
 
 lbm_optimisation = LBMOptimisation(symbolic_field=g, symbolic_temporary_field=g_tmp)
 hydro_lb_update_rule = create_lb_update_rule(lbm_config=config_hydro,
                                              lbm_optimisation=lbm_optimisation)
 
 hydro_lb_update_rule = add_hydrodynamic_force(hydro_lb_update_rule, force_Assignments, C, g, parameters)
 
 ast_kernel = ps.create_kernel(hydro_lb_update_rule, target=dh.default_target, cpu_openmp=True)
 kernel_hydro_lb = ast_kernel.compile()
 ```
 
 %% Output
 
     [['spatialInner0'], ['spatialInner1']]
 
 %% Cell type:markdown id: tags:
 
 ## Boundary Conditions
 
 %% Cell type:code id: tags:
 
 ``` python
 # periodic Boundarys for g, h and C
 periodic_BC_C = dh.synchronization_function(C.name, target=dh.default_target, optimization = {"openmp": True})
 
 periodic_BC_g = LBMPeriodicityHandling(stencil=stencil_hydro, data_handling=dh, pdf_field_name=g.name,
                                        streaming_pattern='push')
 periodic_BC_h = LBMPeriodicityHandling(stencil=stencil_phase, data_handling=dh, pdf_field_name=h.name,
                                        streaming_pattern='pull')
 
 # No slip boundary for the phasefield lbm
 bh_allen_cahn = LatticeBoltzmannBoundaryHandling(method_phase, dh, 'h',
                                                  target=dh.default_target, name='boundary_handling_h',
                                                  streaming_pattern='pull')
 
 # No slip boundary for the velocityfield lbm
 bh_hydro = LatticeBoltzmannBoundaryHandling(method_hydro, dh, 'g' ,
                                             target=dh.default_target, name='boundary_handling_g',
                                             streaming_pattern='push')
 
 contact_angle = BoundaryHandling(dh, C.name, stencil_hydro, target=dh.default_target)
 contact = ContactAngle(90, interface_width)
 
 wall = NoSlip()
 if dimensions == 2:
     bh_allen_cahn.set_boundary(wall, make_slice[:, 0])
     bh_allen_cahn.set_boundary(wall, make_slice[:, -1])
 
     bh_hydro.set_boundary(wall, make_slice[:, 0])
     bh_hydro.set_boundary(wall, make_slice[:, -1])
 
     contact_angle.set_boundary(contact, make_slice[:, 0])
     contact_angle.set_boundary(contact, make_slice[:, -1])
 else:
     bh_allen_cahn.set_boundary(wall, make_slice[:, 0, :])
     bh_allen_cahn.set_boundary(wall, make_slice[:, -1, :])
 
     bh_hydro.set_boundary(wall, make_slice[:, 0, :])
     bh_hydro.set_boundary(wall, make_slice[:, -1, :])
 
     contact_angle.set_boundary(contact, make_slice[:, 0, :])
     contact_angle.set_boundary(contact, make_slice[:, -1, :])
 
 
 bh_allen_cahn.prepare()
 bh_hydro.prepare()
 contact_angle.prepare()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Full timestep
 
 %% Cell type:code id: tags:
 
 ``` python
 # definition of the timestep for the immiscible fluids model
 def timeloop():
     # Solve the interface tracking LB step with boundary conditions
     periodic_BC_h()
     bh_allen_cahn()
     dh.run_kernel(kernel_allen_cahn_lb, **parameters.symbolic_to_numeric_map)
     dh.swap("C", "C_tmp")
 
     # apply the three phase-phase contact angle
     contact_angle()
     # periodic BC of the phase-field
     periodic_BC_C()
 
     # solve the hydro LB step with boundary conditions
     dh.run_kernel(kernel_hydro_lb, **parameters.symbolic_to_numeric_map)
     periodic_BC_g()
     bh_hydro()
 
     # compute density (only for vtk output)
     # must be done BEFORE swapping fields to avoid having outdated values
     compute_density()
 
     # field swaps
     dh.swap("h", "h_tmp")
     dh.swap("g", "g_tmp")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def initialize_hydrostatic_pressure():
     for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
 
         # get y as cell center coordinate, i.e., including shift with 0.5
         y = np.zeros_like(block.midpoint_arrays[1])
         y[:, :] = block.midpoint_arrays[1]
 
         # compute hydrostatic density
         rho_hydrostatic = 1 + 3 * gravitational_acceleration * (y - fluid_depth)
 
         # subtract 1 because PDFs in incompressible LBM are centered around zero in lbmpy
         rho_hydrostatic -= 1
 
         # set equilibrium PDFs with velocity=0 and rho;
         for i in range(0, stencil_hydro.Q, 1):
             block["g"][:,:,i] = method_hydro.weights[i] * rho_hydrostatic[:,:]
             block["g_tmp"][:,:,i] = method_hydro.weights[i] * rho_hydrostatic[:,:]
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def compute_density():
     for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
         # PDFs in incompressible LBM are centered around zero in lbmpy
         # => add 1 to whole sum, i.e., initialize with 1
         block["rho"].fill(1);
 
         # compute density
         for i in range(block["g"].shape[-1]):
             block["rho"][:,:] += block["g"][:,:,i]
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if (gravitational_acceleration != 0):
     initialize_hydrostatic_pressure()
     compute_density()
     plt.scalar_field(dh.cpu_arrays["rho"])
 ```
 
 %% Output
 
 
 %% Cell type:code id: tags:
 
 ``` python
 print("================================= start of the simulation ===================================")
 
 start = time.time()
 
 pbar = tqdm(total=timesteps)
 
 timestep = []
 surface_position = []
 symmetry_norm = []
 mass = []
 
 for i in range(0, timesteps):
 
     sum_c_2 = 0.0
     sum_delta_c_2 = 0.0
 
     # write vtk output
     if(i % vtk_output_frequency == 0):
         if gpu:
             dh.to_cpu("C")
             dh.to_cpu("u")
             dh.to_cpu("rho")
         vtk_writer(i)
 
     # extract data (to be written to file)
     if(i % data_extract_frequency == 0):
         pbar.update(data_extract_frequency)
 
         timestep.append(i)
 
         ny = domain_size[1]
 
         # get index containing phase field value < 0.5
         i1 = np.argmax(dh.cpu_arrays["C"][int(domain_size[0] * 0.5), :] < 0.5)
         i0 = i1 - 1 # index containing phase field value >= 0.5
 
         f0 = dh.cpu_arrays["C"][int(domain_size[0] * 0.5), i0] # phase field value >= 0.5
         f1 = dh.cpu_arrays["C"][int(domain_size[0] * 0.5), i1] # phase field value < 0.5
 
         # coordinate of cell center is index+0.5-1 (0.5 to get to cell center; -1 to remove ghost layer from index)
         y0 = i0 + 0.5 - 1
         y1 = i1 + 0.5 - 1
 
         #interpolate
         surface_position.append( y0 + (y1 - y0) / (f1 - f0) * (0.5 - f0) )
 
         # evaluate symmetry in x-direction
         for y in range(0, domain_size[1] - 1):
             for x in range(0, domain_size[0] - 1):
                 if (x >= domain_size[0] * 0.5):
                     continue
 
                 x_mirrored = domain_size[0] - 1 - x;
                 sum_c_2 += dh.cpu_arrays["C"][x, y]**2
                 sum_delta_c_2 += (dh.cpu_arrays["C"][x, y] - dh.cpu_arrays["C"][x_mirrored, y])**2
 
         symmetry_norm.append( (sum_delta_c_2 / sum_c_2)**2 )
 
         mass.append(np.sum(dh.cpu_arrays["C"]))
 
     timeloop()
 
 pbar.close()
 end = time.time()
 sim_time = end - start
 print("\n")
 print("time needed for the calculation: %4.4f" % sim_time , "seconds")
 
 nrOfCells = np.prod(domain_size)
 mlups = nrOfCells * timesteps / sim_time * 1e-6
 
 print("MLUPS: %4.2f" % mlups)
 ```
 
 %% Output
 
     ================================= start of the simulation ===================================
 
 
     
     
     time needed for the calculation: 99.6110 seconds
     MLUPS: 5.99
 
 %% Cell type:code id: tags:
 
 ``` python
 if gpu:
     dh.to_cpu("C")
     dh.to_cpu("u")
     dh.to_cpu("rho")
 
 plt.scalar_field(dh.cpu_arrays["C"])
 ```
 
 %% Output
 
     <matplotlib.image.AxesImage at 0x12db854f0>
 
 
 %% Cell type:code id: tags:
 
 ``` python
 # non-dimensionalize time and surface position
 t_nd = [value * wave_frequency for value in timestep]
 a_nd = [(value - fluid_depth) / amplitude for value in surface_position]
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 plt.plot(t_nd, mass, color=(0.121, 0.231, 0.4))
 plt.show()
 ```
 
 %% Output
 
 
 %% Cell type:code id: tags:
 
 ``` python
 # plot surface position over time
 plt.xlabel('$t/-$', fontsize=18)
 plt.ylabel('$a/-$', fontsize=18)
 plt.grid(color='black', linestyle='-', linewidth=0.1)
 plt.plot(t_nd, a_nd, color=(0.121, 0.231, 0.4))
 plt.show()
 ```
 
 %% Output
 
 
 %% Cell type:code id: tags:
 
 ``` python
 # plot analytical solutions
 x = 0.5 # position where surface elevation is evaluated;
         # 0.5 means center of first cell;
         # amplitude of analytical solution is at 0 (not in domain center)
 
 t = np.arange(0, timesteps, 1)
 a = amplitude * np.exp(-2 * kinematic_viscosity_heavy * wavenumber**2 * t) # damping of wave amplitude
 
 # linear analytical solution
 a_anal_lin = a * np.cos(wavenumber * x - wave_frequency * t) + fluid_depth
 
 # non-linear analytical solution
 a_anal_nonlin = (a * np.cos(wavenumber * x - wave_frequency * t) +
                         wavenumber * a**2 * 0.5 *
                         (1 + 3 / (2 * np.sinh(wavenumber * fluid_depth)**2)) *
                         1 / np.tanh(wavenumber * fluid_depth) *
                         np.cos(2 * wavenumber * x - 2 * wave_frequency * t) + fluid_depth)
 
 # non-dimesionalize time and surface position
 a_anal_lin_nd = (a_anal_lin - fluid_depth) / amplitude
 a_anal_nonlin_nd = (a_anal_nonlin - fluid_depth) / amplitude
 t_anal_nd = t * wave_frequency
 
 plt.plot(t_anal_nd, a_anal_lin_nd, label="linear analytical solution")
 #plt.plot(t_anal_nd, a_anal_nonlin_nd, label="nonlinear analytical solution")
 plt.plot(t_nd, a_nd, label="simulation")
 plt.xlabel('$t/-$', fontsize=18)
 plt.ylabel('$a/-$', fontsize=18)
 plt.legend()
 plt.grid()
 plt.show()
 ```
 
 %% Output
 
 
 %% Cell type:code id: tags:
 
 ``` python
 # plot symmetry norm over time
 plt.xlabel('$t/-$', fontsize=18)
 plt.ylabel('$symmetry-norm/-$', fontsize=18)
 plt.grid(color='black', linestyle='-', linewidth=0.1)
 plt.plot(t_nd, symmetry_norm, color=(0.121, 0.231, 0.4))
 plt.show()
 ```
 
 %% Output
 
 
 %% Cell type:code id: tags:
 
 ``` python
 # store results to file
 import csv
 with open(filename, 'w') as f:
     writer = csv.writer(f, delimiter='\t')
     writer.writerows(zip(t_nd, a_nd, symmetry_norm))
 ```

lbmpy_tests/full_scenarios/shear_wave/scenario_shear_wave.py → tests/full_scenarios/shear_wave/scenario_shear_wave.py

@@ -2,18 +2,23 @@
     The cumulant lattice Boltzmann equation in three dimensions: Theory and validation
     by  Geier, Martin; Schönherr, Martin; Pasquali, Andrea; Krafczyk, Manfred (2015)
 
-    Chapter 5.1
+    :cite:`geier2015` Chapter 5.1
+
+    NOTE: for integration tests, the parameter study is greatly shortened, i.e., the runs are shortened in time and
+    not all resolutions and viscosities are considered. Nevertheless, all values used by Geier et al. are still in
+    the setup, only commented, and remain ready to be used (check for comments that start with `NOTE`).
 """
 import numpy as np
-import sympy as sp
-
 import pytest
-from pystencils import Target
+import sympy as sp
 
-from lbmpy.creationfunctions import update_with_default_parameters
+from lbmpy import LatticeBoltzmannStep, LBStencil
+from lbmpy.creationfunctions import LBMConfig, LBMOptimisation
+from lbmpy.db import LbmpyJsonSerializer
+from lbmpy.enums import Method, Stencil
 from lbmpy.relaxationrates import (
     relaxation_rate_from_lattice_viscosity, relaxation_rate_from_magic_number)
-from lbmpy.scenarios import create_fully_periodic_flow
+from pystencils import Target, create_data_handling, CreateKernelConfig
 
 
 def get_exponent_term(l, **kwargs):
@@ -53,14 +58,13 @@ def get_analytical_solution(l, l_0, u_0, v_0, nu, t):
 
 def plot_y_velocity(vel, **kwargs):
     import matplotlib.pyplot as plt
-    from matplotlib import cm
-    vel = vel[:, vel.shape[1]//2, :, 1]
+    vel = vel[:, vel.shape[1] // 2, :, 1]
     ranges = [np.arange(n, dtype=float) for n in vel.shape]
     x, y = np.meshgrid(*ranges, indexing='ij')
     fig = plt.gcf()
     ax = fig.gca(projection='3d')
 
-    ax.plot_surface(x, y, vel, cmap=cm.coolwarm, rstride=2, cstride=2,
+    ax.plot_surface(x, y, vel, cmap='coolwarm', rstride=2, cstride=2,
                     linewidth=0, antialiased=True, **kwargs)
 
 
@@ -71,9 +75,9 @@ def fit_and_get_slope(x_values, y_values):
 
 
 def get_phase_error(phases, evaluation_interval):
-    xValues = np.arange(len(phases) * evaluation_interval, step=evaluation_interval)
-    phaseError = fit_and_get_slope(xValues, phases)
-    return phaseError
+    x_values = np.arange(len(phases) * evaluation_interval, step=evaluation_interval)
+    phase_error = fit_and_get_slope(x_values, phases)
+    return phase_error
 
 
 def get_viscosity(abs_values, evaluation_interval, l):
@@ -92,33 +96,50 @@ def get_amplitude_and_phase(vel_slice):
     return fft_abs[max_index], fft_phase[max_index]
 
 
-def run(l, l_0, u_0, v_0, nu, y_size, periodicity_in_kernel, **kwargs):
+def run(l, l_0, u_0, v_0, nu, y_size, lbm_config, lbm_optimisation, config):
     inv_result = {
         'viscosity_error': np.nan,
         'phase_error': np.nan,
         'mlups': np.nan,
     }
-    if 'initial_velocity' in kwargs:
-        del kwargs['initial_velocity']
+    if lbm_config.initial_velocity:
+        # no manually specified initial velocity allowed in shear wave setup
+        lbm_config.initial_velocity = None
 
-    print("Running size l=%d nu=%f " % (l, nu), kwargs)
+    print(f"Running size l={l} nu={nu}")
     initial_vel_arr = get_initial_velocity_field(l, l_0, u_0, v_0, y_size)
     omega = relaxation_rate_from_lattice_viscosity(nu)
 
-    kwargs['relaxation_rates'] = [sp.sympify(r) for r in kwargs['relaxation_rates']]
-    if 'entropic' not in kwargs or not kwargs['entropic']:
-        kwargs['relaxation_rates'] = [r.subs(sp.Symbol("omega"), omega) for r in kwargs['relaxation_rates']]
+    if lbm_config.fourth_order_correction and omega < 1.75:
+        pytest.skip("Fourth-order correction only allowed for omega >= 1.75.")
+
+    lbm_config.relaxation_rates = [sp.sympify(r) for r in lbm_config.relaxation_rates]
+    lbm_config.relaxation_rates = [r.subs(sp.Symbol("omega"), omega) for r in lbm_config.relaxation_rates]
+
+    periodicity_in_kernel = (lbm_optimisation.builtin_periodicity != (False, False, False))
+    domain_size = initial_vel_arr.shape[:-1]
 
-    scenario = create_fully_periodic_flow(initial_vel_arr, periodicity_in_kernel=periodicity_in_kernel, **kwargs)
+    data_handling = create_data_handling(domain_size, periodicity=not periodicity_in_kernel,
+                                         default_ghost_layers=1, parallel=False)
 
-    if 'entropic' in kwargs and kwargs['entropic']:
-        scenario.kernelParams['omega'] = kwargs['relaxation_rates'][0].subs(sp.Symbol("omega"), omega)
+    scenario = LatticeBoltzmannStep(data_handling=data_handling, name="periodic_scenario", lbm_config=lbm_config,
+                                    lbm_optimisation=lbm_optimisation, config=config)
+    for b in scenario.data_handling.iterate(ghost_layers=False):
+        np.copyto(b[scenario.velocity_data_name], initial_vel_arr[b.global_slice])
+    scenario.set_pdf_fields_from_macroscopic_values()
+
+    # NOTE: use those values to limit the runtime in integration tests
+    total_time_steps = 2000 * (l // l_0) ** 2
+    initial_time_steps = 1100 * (l // l_0) ** 2
+    eval_interval = 100 * (l // l_0) ** 2
+    # NOTE: for simulating the real shear-wave scenario from Geier et al. use the following values
+    # total_time_steps = 20000 * (l // l_0) ** 2
+    # initial_time_steps = 11000 * (l // l_0) ** 2
+    # eval_interval = 1000 * (l // l_0) ** 2
 
-    total_time_steps = 20000 * (l // l_0) ** 2
-    initial_time_steps = 11000 * (l // l_0) ** 2
-    eval_interval = 1000 * (l // l_0) ** 2
     scenario.run(initial_time_steps)
     if np.isnan(scenario.velocity_slice()).any():
+        print("   Result", inv_result)
         return inv_result
 
     magnitudes = []
@@ -149,65 +170,81 @@ def run(l, l_0, u_0, v_0, nu, y_size, periodicity_in_kernel, **kwargs):
 def create_full_parameter_study():
     from pystencils.runhelper import ParameterStudy
 
-    params = {
+    setup_params = {
         'l_0': 32,
         'u_0': 0.096,
         'v_0': 0.1,
-        'ySize': 1,
-        'periodicityInKernel': True,
-        'stencil': 'D3Q27',
-        'compressible': True,
+        'y_size': 1
     }
-    ls = [32 * 2 ** i for i in range(0, 5)]
-    nus = [1e-2, 1e-3, 1e-4, 1e-5]
-
-    srt_and_trt_methods = [{'method': method,
-                            'stencil': stencil,
-                            'compressible': comp,
-                            'relaxation_rates': ["omega", str(relaxation_rate_from_magic_number(sp.Symbol("omega")))],
-                            'equilibrium_order': eqOrder,
-                            'continuous_equilibrium': mbEq,
-                            'optimization': {'target': Target.CPU, 'split': True, 'cse_pdfs': True}}
-                           for method in ('srt', 'trt')
-                           for stencil in ('D3Q19', 'D3Q27')
+
+    omega, omega_f = sp.symbols("omega, omega_f")
+
+    # NOTE: use those values to limit the runtime in integration tests
+    ls = [32]
+    nus = [1e-5]
+    # NOTE: for simulating the real shear-wave scenario from Geier et al. use the following values
+    # ls = [32 * 2 ** i for i in range(0, 5)]
+    # nus = [1e-2, 1e-3, 1e-4, 1e-5]
+
+    srt_and_trt_methods = [LBMConfig(method=method,
+                                     stencil=LBStencil(stencil),
+                                     compressible=comp,
+                                     relaxation_rates=[omega, str(relaxation_rate_from_magic_number(omega))],
+                                     equilibrium_order=eqOrder,
+                                     continuous_equilibrium=mbEq)
+                           for method in (Method.SRT, Method.TRT)
+                           for stencil in (Stencil.D3Q19, Stencil.D3Q27)
                            for comp in (True, False)
                            for eqOrder in (1, 2, 3)
                            for mbEq in (True, False)]
 
-    methods = [{'method': 'trt', 'relaxation_rates': ["omega", 1]},
-               {'method': 'mrt3', 'relaxation_rates': ["omega", 1, 1], 'equilibrium_order': 2},
-               {'method': 'mrt3', 'relaxation_rates': ["omega", 1, 1], 'equilibrium_order': 3},
-               {'method': 'mrt3', 'cumulant': True, 'relaxation_rates': ["omega", 1, 1],  # cumulant
-                'continuous_equilibrium': True, 'equilibrium_order': 3,
-                'optimization': {'cse_global': True}},
-               {'method': 'trt-kbc-n4', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True,  # entropic order 2
-                'continuous_equilibrium': True, 'equilibrium_order': 2},
-               {'method': 'trt-kbc-n4', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True,  # entropic order 3
-                'continuous_equilibrium': True, 'equilibrium_order': 3},
-
-               # entropic cumulant:
-               {'method': 'trt-kbc-n4', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True,
-                'cumulant': True, 'continuous_equilibrium': True, 'equilibrium_order': 3},
-               {'method': 'trt-kbc-n2', 'relaxation_rates': ["omega", "omega_f"], 'entropic': True,
-                'cumulant': True, 'continuous_equilibrium': True, 'equilibrium_order': 3},
-               {'method': 'mrt3', 'relaxation_rates': ["omega", "omega_f", "omega_f"], 'entropic': True,
-                'cumulant': True, 'continuous_equilibrium': True, 'equilibrium_order': 3},
+    optimization_srt_trt = LBMOptimisation(split=True, cse_pdfs=True)
+
+    config = CreateKernelConfig(target=Target.CPU)
+
+    methods = [LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.TRT, relaxation_rates=[omega, 1]),
+               LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.MRT, relaxation_rates=[omega],
+                         equilibrium_order=2),
+               LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.MRT, relaxation_rates=[omega],
+                         equilibrium_order=3),
+               LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.CUMULANT, relaxation_rates=[omega],  # cumulant
+                         compressible=True, continuous_equilibrium=True, equilibrium_order=3),
+               LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.CUMULANT, relaxation_rates=[omega],  # cumulant with fourth-order correction
+                         compressible=True, continuous_equilibrium=True, fourth_order_correction=0.1),
+               LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.TRT_KBC_N4, relaxation_rates=[omega, omega_f],
+                         entropic=True, zero_centered=False, continuous_equilibrium=True, equilibrium_order=2),
+               LBMConfig(stencil=LBStencil(Stencil.D3Q27), method=Method.TRT_KBC_N4, relaxation_rates=[omega, omega_f],
+                         entropic=True, zero_centered=False, continuous_equilibrium=True, equilibrium_order=3),
+
+               # entropic cumulant: not supported for the moment
+               # LBMConfig(method=Method.CUMULANT, relaxation_rates=["omega", "omega_f"], entropic=True, zero_centered=False,
+               #           compressible=True, continuous_equilibrium=True, equilibrium_order=3)
                ]
 
-    parameter_study = ParameterStudy(run, database_connector="shear_wave_db")
+    parameter_study = ParameterStudy(run, database_connector="shear_wave_db",
+                                     serializer_info=('lbmpy_serializer', LbmpyJsonSerializer))
     for L in ls:
         for nu in nus:
-            for methodParams in methods + srt_and_trt_methods:
-                simulation_params = params.copy()
-                simulation_params.update(methodParams)
-                if 'optimization' not in simulation_params:
-                    simulation_params['optimization'] = {}
-                new_params, new_opt_params = update_with_default_parameters(simulation_params,
-                                                                            simulation_params['optimization'], False)
-                simulation_params = new_params
-                simulation_params['optimization'] = new_opt_params
-
-                simulation_params['L'] = L
+            for methodParams in srt_and_trt_methods:
+                simulation_params = setup_params.copy()
+
+                simulation_params['lbm_config'] = methodParams
+                simulation_params['lbm_optimisation'] = optimization_srt_trt
+                simulation_params['config'] = config
+
+                simulation_params['l'] = L
+                simulation_params['nu'] = nu
+                l_0 = simulation_params['l_0']
+                parameter_study.add_run(simulation_params.copy(), weight=(L / l_0) ** 4)
+
+            for methodParams in methods:
+                simulation_params = setup_params.copy()
+
+                simulation_params['lbm_config'] = methodParams
+                simulation_params['lbm_optimisation'] = LBMOptimisation()
+                simulation_params['config'] = config
+
+                simulation_params['l'] = L
                 simulation_params['nu'] = nu
                 l_0 = simulation_params['l_0']
                 parameter_study.add_run(simulation_params.copy(), weight=(L / l_0) ** 4)
@@ -215,15 +252,25 @@ def create_full_parameter_study():
 
 
 def test_shear_wave():
-    pytest.importorskip('pycuda')
+    pytest.importorskip('cupy')
     params = {
+        'y_size': 1,
         'l_0': 32,
         'u_0': 0.096,
         'v_0': 0.1,
 
-        'stencil': 'D3Q19',
-        'compressible': True,
-        "optimization": {"target": Target.GPU}
+        'l': 32,
+        'nu': 1e-2,
     }
-    run(32, nu=1e-2, equilibrium_order=2, method='srt', y_size=1, periodicity_in_kernel=True,
-        relaxation_rates=[sp.Symbol("omega"), 5, 5], continuous_equilibrium=True, **params)
+
+    kernel_config = CreateKernelConfig(target=Target.GPU)
+    lbm_config = LBMConfig(method=Method.SRT, relaxation_rates=[sp.Symbol("omega")], stencil=LBStencil(Stencil.D3Q27),
+                           compressible=True, continuous_equilibrium=True, equilibrium_order=2)
+
+    run(**params, lbm_config=lbm_config, config=kernel_config, lbm_optimisation=LBMOptimisation())
+
+
+@pytest.mark.longrun
+def test_all_scenarios():
+    parameter_study = create_full_parameter_study()
+    parameter_study.run()

lbmpy_tests/phasefield/test_n_phase_boyer_analytical.ipynb → tests/phasefield/test_n_phase_boyer_analytical.ipynb

@@ -415,7 +415,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -429,7 +429,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.2"
+   "version": "3.11.0rc1"
   }
  },
  "nbformat": 4,

lbmpy_tests/phasefield/test_n_phase_boyer_noncoupled.ipynb → tests/phasefield/test_n_phase_boyer_noncoupled.ipynb

 %% Cell type:code id: tags:
 
 ``` python
 import pytest
-pytest.importorskip('pycuda')
+pytest.importorskip('cupy')
 ```
 
 %% Output
 
-    ---------------------------------------------------------------------------
-    Skipped                                   Traceback (most recent call last)
-    /var/folders/07/0d7kq8fd0sx24cs53zz90_qc0000gp/T/ipykernel_16968/622163826.py in <module>
-          1 import pytest
-    ----> 2 pytest.importorskip('pycuda')
-
-    /opt/local/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/_pytest/outcomes.py in importorskip(modname, minversion, reason)
-        210             if reason is None:
-        211                 reason = f"could not import {modname!r}: {exc}"
-    --> 212             raise Skipped(reason, allow_module_level=True) from None
-        213     mod = sys.modules[modname]
-        214     if minversion is None:
-    Skipped: could not import 'pycuda': No module named 'pycuda'
+    <module 'cupy' from '/home/markus/.local/lib/python3.11/site-packages/cupy/__init__.py'>
 
 %% Cell type:code id: tags:
 
 ``` python
 from lbmpy.session import *
 from lbmpy.phasefield.n_phase_boyer import *
 from lbmpy.phasefield.kerneleqs import *
 from lbmpy.phasefield.contact_angle_circle_fitting import *
-from scipy.ndimage.filters import gaussian_filter
+from scipy.ndimage import gaussian_filter
 from pystencils.simp import sympy_cse_on_assignment_list
 one = sp.sympify(1)
 
 import pyximport
 pyximport.install(language_level=3)
 from lbmpy.phasefield.simplex_projection import simplex_projection_2d  # NOQA
 ```
 
 %% Cell type:markdown id: tags:
 
 # Simulation arbitrary surface tension case
 
 %% Cell type:code id: tags:
 
 ``` python
 n = 4
 dx, dt = 1, 1
 mobility = 2e-3
 domain_size = (150, 150)
 ε = one * 4
 penalty_factor = 0
 stabilization_factor = 10
 
 κ = (one,  one/2, one/3, one/4)
 sigma_factor = one / 15
 σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 )
 #σ
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 dh = create_data_handling(domain_size, periodicity=True, default_target=ps.Target.GPU)
 c = dh.add_array('c', values_per_cell=n)
 c_tmp = dh.add_array_like('c_tmp', 'c')
 
 μ = dh.add_array('mu', values_per_cell=n)
 
 cvec = c.center_vector
 μvec = μ.center_vector
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 α, _ = diffusion_coefficients(σ)
 
 f = lambda c: c**2 * ( 1 - c ) **2
 a, b = compute_ab(f)
 
 capital_f = capital_f0(cvec, σ) + correction_g(cvec, σ) + stabilization_factor * stabilization_term(cvec, α)
 
 f_bulk = free_energy_bulk(capital_f, b, ε) + penalty_factor * (one - sum(cvec))
 f_if = free_energy_interfacial(cvec, σ, a, ε)
 f = f_bulk + f_if
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 #f_bulk
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 μ_assignments = mu_kernel(f, cvec, c, μ)
 μ_assignments = [Assignment(a.lhs, a.rhs.doit()) for a in μ_assignments]
 μ_assignments = sympy_cse_on_assignment_list(μ_assignments)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 discretize = fd.Discretization2ndOrder(dx=dx, dt=dt)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def lapl(e):
     return sum(ps.fd.diff(e, d, d) for d in range(dh.dim))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 rhs = α * μvec
 discretized_rhs = [discretize(fd.expand_diff_full( lapl(mobility * rhs_i) + fd.transient(cvec[i], idx=i), functions=μvec))
                    for i, rhs_i in enumerate(rhs)]
 c_assignments = [Assignment(lhs, rhs)
                  for lhs, rhs in zip(c_tmp.center_vector, discretized_rhs)]
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 #c_assignments
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 μ_sync = dh.synchronization_function(μ.name)
 c_sync = dh.synchronization_function(c.name)
 optimization = {'cpu_openmp': False, 'cpu_vectorize_info': None}
 
 config = ps.CreateKernelConfig(cpu_openmp=False, target=dh.default_target)
 
 μ_kernel = create_kernel(μ_assignments, config=config).compile()
 c_kernel = create_kernel(c_assignments, config=config).compile()
 
 def set_c(slice_obj, values):
     for block in dh.iterate(slice_obj):
         arr = block[c.name]
         arr[..., : ] = values
 
 def smooth():
     for block in dh.iterate(ghost_layers=True):
         c_arr = block[c.name]
         for i in range(n):
             gaussian_filter(c_arr[..., i], sigma=2, output=c_arr[..., i])
 
 def time_loop(steps):
     dh.all_to_gpu()
     for t in range(steps):
         c_sync()
         dh.run_kernel(μ_kernel)
         μ_sync()
         dh.run_kernel(c_kernel)
         dh.swap(c.name, c_tmp.name)
         #simplex_projection_2d(dh.cpu_arrays[c.name])
     dh.all_to_cpu()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 set_c(make_slice[:, :], [0, 0, 0, 0])
 set_c(make_slice[:, 0.5:], [1, 0, 0, 0])
 set_c(make_slice[:, :0.5], [0, 1, 0, 0])
 set_c(make_slice[0.3:0.7, 0.3:0.7], [0, 0, 1, 0])
 smooth()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 #dh.load_all('n_phases_state_size200_stab10.npz')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 plt.phase_plot(dh.gather_array(c.name))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 neumann_angles_from_surface_tensions(lambda i, j: float(σ[i, j]))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 import time
 for i in range(10):
     start = time.perf_counter()
     time_loop(1_000)
     end = time.perf_counter()
 
     try:
         print(i, end - start, liquid_lens_neumann_angles(dh.gather_array(c.name)))
     except Exception:
         print(i, end - start, "none found")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 plt.subplot(1,3,1)
 t = dh.gather_array(c.name, make_slice[25, :]).squeeze()
 plt.plot(t);
 
 plt.subplot(1,3,2)
 plt.phase_plot(dh.gather_array(c.name), linewidth=1)
 
 plt.subplot(1,3,3)
 plt.scalar_field(dh.gather_array(μ.name)[:, :, 2])
 plt.colorbar();
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 assert not np.isnan(dh.max(c.name))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 t = dh.gather_array(c.name, make_slice[25, 55:90]).squeeze()
 plt.hlines(0.5, 0, 30)
 plt.plot(t);
 ```