lbmpy/methods/cumulantbased/cumulantbasedmethod.py

@@ -23,7 +23,7 @@ class CumulantBasedLbMethod(AbstractLbMethod):
     This class implements cumulant-based lattice boltzmann methods which relax all the non-conserved quantities
     as either monomial or polynomial cumulants. It is mostly inspired by the work presented in :cite:`geier2015`.
 
-    This method is implemented modularily as the transformation from populations to central moments to cumulants
+    This method is implemented modularly as the transformation from populations to central moments to cumulants
     is governed by subclasses of :class:`lbmpy.moment_transforms.AbstractMomentTransform`
     which can be specified by constructor argument. This allows the selection of the most efficient transformation
     for a given setup.
@@ -124,7 +124,7 @@ class CumulantBasedLbMethod(AbstractLbMethod):
     @property
     def zeroth_order_equilibrium_moment_symbol(self, ):
         """Returns a symbol referring to the zeroth-order moment of this method's equilibrium distribution,
-        which is the area under it's curve
+        which is the area under its curve
         (see :attr:`lbmpy.equilibrium.AbstractEquilibrium.zeroth_order_moment_symbol`)."""
         return self._equilibrium.zeroth_order_moment_symbol
 

lbmpy/methods/cumulantbased/fourth_order_correction.py

+import sympy as sp
+
+from lbmpy.moment_transforms import PRE_COLLISION_MONOMIAL_CUMULANT, POST_COLLISION_CUMULANT
+from lbmpy.moments import MOMENT_SYMBOLS, statistical_quantity_symbol
+from lbmpy.stencils import Stencil, LBStencil
+from pystencils import Assignment
+from pystencils.simp.assignment_collection import AssignmentCollection
+from .cumulantbasedmethod import CumulantBasedLbMethod
+
+X, Y, Z = MOMENT_SYMBOLS
+CORRECTED_FOURTH_ORDER_POLYNOMIALS = [X ** 2 * Y ** 2 - 2 * X ** 2 * Z ** 2 + Y ** 2 * Z ** 2,
+                                      X ** 2 * Y ** 2 + X ** 2 * Z ** 2 - 2 * Y ** 2 * Z ** 2,
+                                      X ** 2 * Y ** 2 + X ** 2 * Z ** 2 + Y ** 2 * Z ** 2,
+                                      X ** 2 * Y * Z,
+                                      X * Y ** 2 * Z,
+                                      X * Y * Z ** 2]
+FOURTH_ORDER_CORRECTION_SYMBOLS = sp.symbols("corr_fourth_:6")
+FOURTH_ORDER_RELAXATION_RATE_SYMBOLS = sp.symbols("corr_rr_:10")
+
+
+def add_fourth_order_correction(collision_rule, shear_relaxation_rate, bulk_relaxation_rate, limiter):
+    """Adds the fourth order correction terms (:cite:`geier2017`, eq. 44-49) to a given polynomial D3Q27
+    cumulant collision rule."""
+    method = collision_rule.method
+
+    if not isinstance(method, CumulantBasedLbMethod) or method.stencil != LBStencil(Stencil.D3Q27):
+        raise ValueError("Fourth-order correction is only defined for D3Q27 cumulant methods.")
+
+    polynomials = method.cumulants
+    rho = method.zeroth_order_equilibrium_moment_symbol
+
+    if not (set(CORRECTED_FOURTH_ORDER_POLYNOMIALS) < set(polynomials)):
+        raise ValueError("Fourth order correction requires polynomial cumulants "
+                         f"{', '.join(CORRECTED_FOURTH_ORDER_POLYNOMIALS)} to be present")
+
+    #   Call
+    #   PC1 = (x^2 * y^2 - 2 * x^2 * z^2 + y^2 * z^2),
+    #   PC2 = (x^2 * y^2 + x^2 * z^2 - 2 * y^2 * z^2),
+    #   PC3 = (x^2 * y^2 + x^2 * z^2 + y^2 * z^2),
+    #   PC4 = (x^2 * y * z),
+    #   PC5 = (x * y^2 * z),
+    #   PC6 = (x * y * z^2)
+    poly_symbols = [sp.Symbol(f'{POST_COLLISION_CUMULANT}_{polynomials.index(poly)}')
+                    for poly in CORRECTED_FOURTH_ORDER_POLYNOMIALS]
+
+    a_symb, b_symb = sp.symbols("a_corr, b_corr")
+    a, b = get_optimal_additional_parameters(shear_relaxation_rate, bulk_relaxation_rate)
+    correction_terms = get_fourth_order_correction_terms(method.relaxation_rate_dict, rho, a_symb, b_symb)
+    optimal_parametrisation = get_optimal_parametrisation_with_limiters(collision_rule, shear_relaxation_rate,
+                                                                        bulk_relaxation_rate, limiter)
+
+    subexp_dict = collision_rule.subexpressions_dict
+    subexp_dict = {**subexp_dict,
+                   a_symb: a,
+                   b_symb: b,
+                   **correction_terms.subexpressions_dict,
+                   **optimal_parametrisation.subexpressions_dict,
+                   **correction_terms.main_assignments_dict,
+                   **optimal_parametrisation.main_assignments_dict}
+    for sym, cor in zip(poly_symbols, FOURTH_ORDER_CORRECTION_SYMBOLS):
+        subexp_dict[sym] += cor
+
+    collision_rule.set_sub_expressions_from_dict(subexp_dict)
+    collision_rule.topological_sort()
+
+    return collision_rule
+
+
+def get_optimal_additional_parameters(shear_relaxation_rate, bulk_relaxation_rate):
+    """
+    Calculates the optimal parametrization for the additional parameters provided in :cite:`geier2017`
+    Equations 115-116.
+    """
+
+    omega_1 = shear_relaxation_rate
+    omega_2 = bulk_relaxation_rate
+
+    omega_11 = omega_1 * omega_1
+    omega_12 = omega_1 * omega_2
+    omega_22 = omega_2 * omega_2
+
+    two = sp.Float(2)
+    three = sp.Float(3)
+    four = sp.Float(4)
+    nine = sp.Float(9)
+
+    denom_ab = (omega_1 - omega_2) * (omega_2 * (two + three * omega_1) - sp.Float(8) * omega_1)
+
+    a = (four * omega_11 + two * omega_12 * (omega_1 - sp.Float(6)) + omega_22 * (
+        omega_1 * (sp.Float(10) - three * omega_1) - four)) / denom_ab
+    b = (four * omega_12 * (nine * omega_1 - sp.Float(16)) - four * omega_11 - two * omega_22 * (
+        two + nine * omega_1 * (omega_1 - two))) / (three * denom_ab)
+
+    return a, b
+
+
+def get_fourth_order_correction_terms(rrate_dict, rho, a, b):
+    pc1, pc2, pc3, pc4, pc5, pc6 = CORRECTED_FOURTH_ORDER_POLYNOMIALS
+
+    omega_1 = rrate_dict[X * Y]
+    omega_2 = rrate_dict[X ** 2 + Y ** 2 + Z ** 2]
+
+    try:
+        omega_6 = rrate_dict[pc1]
+        assert omega_6 == rrate_dict[pc2], \
+            "Cumulants (x^2 - y^2) and (x^2 - z^2) must have the same relaxation rate"
+        omega_7 = rrate_dict[pc3]
+        omega_8 = rrate_dict[pc4]
+        assert omega_8 == rrate_dict[pc5] == rrate_dict[pc6]
+    except IndexError:
+        raise ValueError("For the fourth order correction, all six polynomial cumulants"
+                         "(x^2 * y^2 - 2 * x^2 * z^2 + y^2 * z^2),"
+                         "(x^2 * y^2 + x^2 * z^2 - 2 * y^2 * z^2),"
+                         "(x^2 * y^2 + x^2 * z^2 + y^2 * z^2),"
+                         "(x^2 * y * z), (x * y^2 * z) and (x * y * z^2) must be present!")
+
+    dxu, dyv, dzw, dxvdyu, dxwdzu, dywdzv = sp.symbols('Dx, Dy, Dz, DxvDyu, DxwDzu, DywDzv')
+    c_xy = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (1, 1, 0))
+    c_xz = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (1, 0, 1))
+    c_yz = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (0, 1, 1))
+    c_xx = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (2, 0, 0))
+    c_yy = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (0, 2, 0))
+    c_zz = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (0, 0, 2))
+
+    cor1, cor2, cor3, cor4, cor5, cor6 = FOURTH_ORDER_CORRECTION_SYMBOLS
+
+    one = sp.Float(1)
+    two = sp.Float(2)
+    three = sp.Float(3)
+
+    #   Derivative Approximations
+    subexpressions = [
+        Assignment(dxu, - omega_1 / (two * rho) * (two * c_xx - c_yy - c_zz)
+                   - omega_2 / (two * rho) * (c_xx + c_yy + c_zz - rho)),
+        Assignment(dyv, dxu + (three * omega_1) / (two * rho) * (c_xx - c_yy)),
+        Assignment(dzw, dxu + (three * omega_1) / (two * rho) * (c_xx - c_zz)),
+        Assignment(dxvdyu, - three * omega_1 / rho * c_xy),
+        Assignment(dxwdzu, - three * omega_1 / rho * c_xz),
+        Assignment(dywdzv, - three * omega_1 / rho * c_yz)]
+
+    one_half = sp.Rational(1, 2)
+    one_over_three = sp.Rational(1, 3)
+    two_over_three = sp.Rational(2, 3)
+    four_over_three = sp.Rational(4, 3)
+
+    #   Correction Terms
+    main_assignments = [
+        Assignment(cor1, two_over_three * (one / omega_1 - one_half) * omega_6 * a * rho * (dxu - two * dyv + dzw)),
+        Assignment(cor2, two_over_three * (one / omega_1 - one_half) * omega_6 * a * rho * (dxu + dyv - two * dzw)),
+        Assignment(cor3, - four_over_three * (one / omega_1 - one_half) * omega_7 * a * rho * (dxu + dyv + dzw)),
+        Assignment(cor4, - one_over_three * (one / omega_1 - one_half) * omega_8 * b * rho * dywdzv),
+        Assignment(cor5, - one_over_three * (one / omega_1 - one_half) * omega_8 * b * rho * dxwdzu),
+        Assignment(cor6, - one_over_three * (one / omega_1 - one_half) * omega_8 * b * rho * dxvdyu)]
+
+    return AssignmentCollection(main_assignments=main_assignments, subexpressions=subexpressions)
+
+
+def get_optimal_parametrisation_with_limiters(collision_rule, shear_relaxation_rate, bulk_relaxation_rate, limiter):
+    """
+    Calculates the optimal parametrization for the relaxation rates provided in :cite:`geier2017`
+    Equations 112-114.
+    """
+
+    # if omega numbers
+    # assert omega_1 != omega_2, "Relaxation rates associated with shear and bulk viscosity must not be identical."
+    # assert omega_1 > 7/4
+
+    omega_1 = FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[0]
+    omega_2 = FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[1]
+
+    non_limited_omegas = sp.symbols("non_limited_omega_3:6")
+    limited_omegas = sp.symbols("limited_omega_3:6")
+
+    omega_11 = omega_1 * omega_1
+    omega_12 = omega_1 * omega_2
+    omega_22 = omega_2 * omega_2
+
+    one = sp.Float(1)
+    two = sp.Float(2)
+    three = sp.Float(3)
+    five = sp.Float(5)
+    six = sp.Float(6)
+    seven = sp.Float(7)
+    eight = sp.Float(8)
+    nine = sp.Float(9)
+
+    omega_3 = (eight * (omega_1 - two) * (omega_2 * (three * omega_1 - one) - five * omega_1)) / (
+        eight * (five - two * omega_1) * omega_1 + omega_2 * (eight + omega_1 * (nine * omega_1 - sp.Float(26))))
+
+    omega_4 = (eight * (omega_1 - two) * (omega_1 + omega_2 * (three * omega_1 - seven))) / (
+        omega_2 * (sp.Float(56) - sp.Float(42) * omega_1 + nine * omega_11) - eight * omega_1)
+
+    omega_5 = (sp.Float(24) * (omega_1 - two) * (sp.Float(4) * omega_11 + omega_12 * (
+        sp.Float(18) - sp.Float(13) * omega_1) + omega_22 * (two + omega_1 * (
+            six * omega_1 - sp.Float(11))))) / (sp.Float(16) * omega_11 * (omega_1 - six) - two * omega_12 * (
+                sp.Float(216) + five * omega_1 * (nine * omega_1 - sp.Float(46))) + omega_22 * (omega_1 * (
+                    three * omega_1 - sp.Float(10)) * (sp.Float(15) * omega_1 - sp.Float(28)) - sp.Float(48)))
+
+    rho = collision_rule.method.zeroth_order_equilibrium_moment_symbol
+
+    # add limiters to improve stability
+    c_xyy = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (1, 2, 0))
+    c_xzz = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (1, 0, 2))
+    c_xyz = statistical_quantity_symbol(PRE_COLLISION_MONOMIAL_CUMULANT, (1, 1, 1))
+    abs_c_xyy_plus_xzz = sp.Abs(c_xyy + c_xzz)
+    abs_c_xyy_minus_xzz = sp.Abs(c_xyy - c_xzz)
+    abs_c_xyz = sp.Abs(c_xyz)
+
+    limited_omega_3 = non_limited_omegas[0] + ((one - non_limited_omegas[0]) * abs_c_xyy_plus_xzz) / \
+        (rho * limiter + abs_c_xyy_plus_xzz)
+    limited_omega_4 = non_limited_omegas[1] + ((one - non_limited_omegas[1]) * abs_c_xyy_minus_xzz) / \
+        (rho * limiter + abs_c_xyy_minus_xzz)
+    limited_omega_5 = non_limited_omegas[2] + ((one - non_limited_omegas[2]) * abs_c_xyz) / (rho * limiter + abs_c_xyz)
+
+    subexpressions = [
+        Assignment(non_limited_omegas[0], omega_3),
+        Assignment(non_limited_omegas[1], omega_4),
+        Assignment(non_limited_omegas[2], omega_5),
+        Assignment(limited_omegas[0], limited_omega_3),
+        Assignment(limited_omegas[1], limited_omega_4),
+        Assignment(limited_omegas[2], limited_omega_5)]
+
+    #   Correction Terms
+    main_assignments = [
+        Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[0], shear_relaxation_rate),
+        Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[1], bulk_relaxation_rate),
+        Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[2], limited_omegas[0]),
+        Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[3], limited_omegas[1]),
+        Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[4], limited_omegas[2]),
+        Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[5], one),
+        Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[6], one),
+        Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[7], one),
+        Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[8], one),
+        Assignment(FOURTH_ORDER_RELAXATION_RATE_SYMBOLS[9], one),
+    ]
+
+    return AssignmentCollection(main_assignments=main_assignments, subexpressions=subexpressions)

lbmpy/methods/momentbased/entropic.py

@@ -286,7 +286,7 @@ def _get_relaxation_rates(collision_rule):
     omega_s = get_shear_relaxation_rate(method)
 
     # if the shear relaxation rate is not specified as a symbol look for its symbolic counter part in the subs dict
-    for symbolic_rr, rr in method.subs_dict_relxation_rate.items():
+    for symbolic_rr, rr in method.subs_dict_relaxation_rate.items():
         if omega_s == rr:
             omega_s = symbolic_rr
 

lbmpy/methods/momentbased/momentbasedmethod.py

@@ -16,8 +16,8 @@ from lbmpy.moment_transforms import PdfsToMomentsByChimeraTransform
 class MomentBasedLbMethod(AbstractLbMethod):
     """
     Moment based LBM is a class to represent the single (SRT), two (TRT) and multi relaxation time (MRT) methods.
-    These methods work by transforming the pdfs into moment space using a linear transformation. In the moment
-    space each component (moment) is relaxed to an equilibrium moment by a certain relaxation rate. These
+    These methods work by transforming the pdfs into moment space using a linear transformation. In the moment-space
+    each component (moment) is relaxed to an equilibrium moment by a certain relaxation rate. These
     equilibrium moments can e.g. be determined by taking the equilibrium moments of the continuous Maxwellian.
 
     Parameters:
@@ -378,7 +378,7 @@ class MomentBasedLbMethod(AbstractLbMethod):
             subexpressions += m_post_to_f_post_eqs.subexpressions
             main_assignments = m_post_to_f_post_eqs.main_assignments
         else:
-            #   For SRT, TRT by default, and whenever customly required, derive equations entirely in
+            #   For SRT, TRT by default, and whenever customarily required, derive equations entirely in
             #   population space
             if self._zero_centered and not self._equilibrium.deviation_only:
                 raise Exception("Can only derive population-space equations for zero-centered storage"

lbmpy/moment_transforms/cumulanttransforms.py

@@ -67,10 +67,10 @@ class CentralMomentsToCumulantsByGeneratingFunc(AbstractMomentTransform):
         self.cumulant_exponents = self.moment_exponents
         self.cumulant_polynomials = self.moment_polynomials
 
-        if(len(self.cumulant_exponents) != stencil.Q):
+        if len(self.cumulant_exponents) != stencil.Q:
             raise ValueError("Number of cumulant exponent tuples must match stencil size.")
 
-        if(len(self.cumulant_polynomials) != stencil.Q):
+        if len(self.cumulant_polynomials) != stencil.Q:
             raise ValueError("Number of cumulant polynomials must match stencil size.")
 
         self.central_moment_exponents = self.compute_required_central_moments()

lbmpy/moments.py

@@ -86,8 +86,8 @@ def moment_permutations(exponent_tuple):
 def moments_of_order(order, dim=3, include_permutations=True):
     """All tuples of length 'dim' which sum equals 'order'"""
     for item in __fixed_sum_tuples(dim, order, ordered=not include_permutations):
-        assert(len(item) == dim)
-        assert(sum(item) == order)
+        assert len(item) == dim
+        assert sum(item) == order
         yield item
 
 

lbmpy/scenarios.py

@@ -42,10 +42,10 @@ def create_fully_periodic_flow(initial_velocity, periodicity_in_kernel=False, lb
         initial_velocity: numpy array that defines an initial velocity for each cell. The shape of this
                          array determines the domain size.
         periodicity_in_kernel: don't use boundary handling for periodicity, but directly generate the kernel periodic
-        lbm_kernel: a LBM function, which would otherwise automatically created
+        lbm_kernel: an LBM function, which would otherwise automatically created
         data_handling: data handling instance that is used to create the necessary fields. If a data handling is
                        passed the periodicity and parallel arguments are ignored.
-        parallel: True for distributed memory parallelization with walberla
+        parallel: True for distributed memory parallelization with waLBerla
         kwargs: other parameters are passed on to the method, see :mod:`lbmpy.creationfunctions`
 
     Returns:
@@ -76,9 +76,9 @@ def create_lid_driven_cavity(domain_size=None, lid_velocity=0.005, lbm_kernel=No
     Args:
         domain_size: tuple specifying the number of cells in each dimension
         lid_velocity: x velocity of lid in lattice coordinates.
-        lbm_kernel: a LBM function, which would otherwise automatically created
+        lbm_kernel: an LBM function, which would otherwise automatically created
         kwargs: other parameters are passed on to the method, see :mod:`lbmpy.creationfunctions`
-        parallel: True for distributed memory parallelization with walberla
+        parallel: True for distributed memory parallelization with waLBerla
         data_handling: see documentation of :func:`create_fully_periodic_flow`
     Returns:
         instance of :class:`Scenario`

lbmpy_tests/advanced_streaming/test_fully_periodic_flow.py

@@ -20,7 +20,7 @@ all_results = dict()
 targets = [Target.CPU]
 
 try:
-    import pycuda.autoinit
+    import cupy
     targets += [Target.GPU]
 except Exception:
     pass

lbmpy_tests/advanced_streaming/test_periodic_pipe_with_force.py

@@ -21,7 +21,7 @@ all_results = dict()
 targets = [Target.CPU]
 
 try:
-    import pycuda.autoinit
+    import cupy
     targets += [Target.GPU]
 except Exception:
     pass

lbmpy_tests/cumulantmethod/test_flow_around_sphere.py

@@ -16,7 +16,7 @@ from pystencils import create_kernel, create_data_handling, Assignment, Target,
 from pystencils.slicing import slice_from_direction, get_slice_before_ghost_layer
 
 
-def flow_around_sphere(stencil, galilean_correction, L_LU, total_steps):
+def flow_around_sphere(stencil, galilean_correction, fourth_order_correction, L_LU, total_steps):
     if galilean_correction and stencil.Q != 27:
         pytest.skip()
 
@@ -38,7 +38,7 @@ def flow_around_sphere(stencil, galilean_correction, L_LU, total_steps):
 
     lbm_config = LBMConfig(stencil=stencil, method=Method.CUMULANT, relaxation_rate=omega_v,
                            compressible=True,
-                           galilean_correction=galilean_correction)
+                           galilean_correction=galilean_correction, fourth_order_correction=fourth_order_correction)
 
     lbm_opt = LBMOptimisation(pre_simplification=True)
     config = CreateKernelConfig(target=target)
@@ -135,14 +135,22 @@ def flow_around_sphere(stencil, galilean_correction, L_LU, total_steps):
 
 @pytest.mark.parametrize('stencil', [Stencil.D2Q9, Stencil.D3Q19, Stencil.D3Q27])
 @pytest.mark.parametrize('galilean_correction', [False, True])
-def test_flow_around_sphere_short(stencil, galilean_correction):
-    pytest.importorskip('pycuda')
-    flow_around_sphere(LBStencil(stencil), galilean_correction, 5, 200)
+@pytest.mark.parametrize('fourth_order_correction', [False, True])
+def test_flow_around_sphere_short(stencil, galilean_correction, fourth_order_correction):
+    pytest.importorskip('cupy')
+    stencil = LBStencil(stencil)
+    if fourth_order_correction and stencil.Q != 27:
+        pytest.skip("Fourth-order correction only defined for D3Q27 stencil.")
+    flow_around_sphere(stencil, galilean_correction, fourth_order_correction, 5, 200)
 
 
 @pytest.mark.parametrize('stencil', [Stencil.D2Q9, Stencil.D3Q19, Stencil.D3Q27])
 @pytest.mark.parametrize('galilean_correction', [False, True])
+@pytest.mark.parametrize('fourth_order_correction', [False, 0.01])
 @pytest.mark.longrun
-def test_flow_around_sphere_long(stencil, galilean_correction):
-    pytest.importorskip('pycuda')
-    flow_around_sphere(LBStencil(stencil), galilean_correction, 20, 3000)
+def test_flow_around_sphere_long(stencil, galilean_correction, fourth_order_correction):
+    pytest.importorskip('cupy')
+    stencil = LBStencil(stencil)
+    if fourth_order_correction and stencil.Q != 27:
+        pytest.skip("Fourth-order correction only defined for D3Q27 stencil.")
+    flow_around_sphere(stencil, galilean_correction, fourth_order_correction, 20, 3000)

lbmpy_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

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

@@ -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(method=Method.TRT, relaxation_rates=[omega, 1]),
+               LBMConfig(method=Method.MRT, relaxation_rates=[omega], equilibrium_order=2),
+               LBMConfig(method=Method.MRT, relaxation_rates=[omega], equilibrium_order=3),
+               LBMConfig(method=Method.CUMULANT, relaxation_rates=[omega],  # cumulant
+                         compressible=True, continuous_equilibrium=True, equilibrium_order=3),
+               LBMConfig(method=Method.CUMULANT, relaxation_rates=[omega],  # cumulant with fourth-order correction
+                         compressible=True, continuous_equilibrium=True, fourth_order_correction=0.1),
+               LBMConfig(method=Method.TRT_KBC_N4, relaxation_rates=[omega, omega_f], entropic=True,
+                         zero_centered=False,  # entropic order 2
+                         continuous_equilibrium=True, equilibrium_order=2),
+               LBMConfig(method=Method.TRT_KBC_N4, relaxation_rates=[omega, omega_f], entropic=True,
+                         zero_centered=False,  # entropic order 3
+                         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

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

 %% 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 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);
 ```

lbmpy_tests/phasefield/test_numerical_1D_nphase_model.ipynb

 %% Cell type:code id: tags:
 ``` python
 import pytest
-pytest.importorskip('pycuda')
+pytest.importorskip('cupy')
 ```
 %% Output
-    <module 'pycuda' from '/home/markus/miniconda3/envs/pystencils/lib/python3.8/site-packages/pycuda/__init__.py'>
+    <module 'cupy' from '/home/markus/.local/lib/python3.11/site-packages/cupy/__init__.py'>
 %% Cell type:code id: tags:
 ``` python
 from pystencils.datahandling import SerialDataHandling
 from pystencils import Target
 from lbmpy.session import *
 from lbmpy.phasefield.scenarios import create_n_phase_model
 from lbmpy.phasefield.analytical import *
 from lbmpy.phasefield.experiments1D import *
 from lbmpy.phasefield.analytical import analytic_interface_profile
 from functools import partial
 ```
 %% Cell type:markdown id: tags:
 ## Tanh test
 %% Cell type:code id: tags:
 ``` python
 width = 100
 f2 = partial(n_phases_correction_function_sign_switch, beta=1)
 sc1 = create_n_phase_model(data_handling=SerialDataHandling((width, 1),
                                                             periodicity=True),
                            f2=f2,
                            num_phases=4, alpha=1)
 phaseIdx = 1
 init_sharp_interface(sc1, phase_idx=1, inverse=False)
 #init_sharp_interface(sc1, phase_idx=0, inverse=True)
 sc1.set_pdf_fields_from_macroscopic_values()
 f2(sp.Symbol("c"))
 ```
 %% Output
     $\displaystyle \begin{cases} - c^{2} \left(1 - c\right)^{2} & \text{for}\: c > 1 \vee c < 0 \\c^{2} \left(1 - c\right)^{2} & \text{otherwise} \end{cases}$
     ⎧  2        2
     ⎪-c ⋅(1 - c)   for c > 1 ∨ c < 0
     ⎨
     ⎪ 2        2
     ⎩c ⋅(1 - c)        otherwise
 %% Cell type:code id: tags:
 ``` python
 if 'is_test_run' in globals():
     sc1.run(1000)
 else:
     sc1.run(100_000)
 plot_status(sc1)
 sc1.time_steps_run
 ```
 %% Output
     $\displaystyle 100000$
     100000
 %% Cell type:code id: tags:
 ``` python
 alpha = 1
 visWidth = 25
 x = np.arange(2 * visWidth) - visWidth
 analytic = np.array([analytic_interface_profile(x_i - 0.5, alpha) for x_i in x], dtype=np.float64)
 visSlice = make_slice[25 - visWidth:25 + visWidth, 0, phaseIdx]
 simulated = sc1.phi_slice(visSlice).copy()
 plt.plot(analytic, label='analytic', marker='o')
 plt.plot(simulated, label='simulated', marker='x')
 plt.legend();
 ```
 %% Output
 %% Cell type:code id: tags:
 ``` python
 assert not np.isnan(np.max(simulated))
 ```
 %% Cell type:markdown id: tags:
 # Phase separation
 %% Cell type:code id: tags:
 ``` python
 num_phases = 3
 f2 = partial(n_phases_correction_function, power=2, beta=0.001)
 ll = create_n_phase_model(data_handling=SerialDataHandling((200, 200),
                                                             periodicity=True),
                           f2=f2,
                           surface_tensions=lambda i, j: 0.0025/2 if i != j else 0,
                           num_phases=num_phases, alpha=1,
                           cahn_hilliard_relaxation_rates=1.2,
                           optimization={'target': Target.GPU})
 ```
 %% Cell type:code id: tags:
 ``` python
 ll.set_concentration(make_slice[:, :], [1/num_phases] * (ll.num_order_parameters + 1))
 φ_arr = ll.data_handling.cpu_arrays['pf_phi']
 φ_arr += (np.random.rand(*φ_arr.shape)-0.5) * 0.2
 ll.set_pdf_fields_from_macroscopic_values()
 #plt.phase_plot_for_step(ll)
 ```
 %% Cell type:code id: tags:
 ``` python
 f_bulk, f_if = separate_into_bulk_and_interface(ll.free_energy)
 φ = sp.symbols("phi_:{}".format(num_phases))
 hessian = sp.hessian(f_bulk, φ)
 ```
 %% Cell type:code id: tags:
 ``` python
 evaluated_hessian = hessian.subs({ f: 1/num_phases for f in φ})
 eigenvalues = [a.evalf()  for a,b in evaluated_hessian.eigenvals().items()]
 assert any(a < 0 for a in eigenvalues)
 evaluated_hessian, eigenvalues
 ```
 %% Output
     $\displaystyle \left( \left[\begin{matrix}-0.0025 & -0.00125 & 0\\-0.00125 & -0.0025 & 0\\0 & 0 & 0\end{matrix}\right], \  \left[ -0.00375, \  -0.00125, \  0\right]\right)$
     ⎛⎡-0.0025   -0.00125  0⎤                         ⎞
     ⎜⎢                     ⎥                         ⎟
     ⎜⎢-0.00125  -0.0025   0⎥, [-0.00375, -0.00125, 0]⎟
     ⎜⎢                     ⎥                         ⎟
     ⎝⎣   0         0      0⎦                         ⎠
 %% Cell type:code id: tags:
 ``` python
 eigenvalues
 ```
 %% Output
     $\displaystyle \left[ -0.00375, \  -0.00125, \  0\right]$
     [-0.00375, -0.00125, 0]
 %% Cell type:code id: tags:
 ``` python
 if 'is_test_run' in globals():
     ll.run(1000)
 else:
     for i in range(24):
         ll.run(2_000)
         min_φ, max_φ = np.min(ll.phi[:, :]), np.max(ll.phi[:, :])
         print(i, min_φ, max_φ)
         if max_φ > 0.6:
             break
     assert max_φ > 0.6
     ll.time_steps_run
 ```
 %% Output
     0 0.2844762021511006 0.3838781484847188
     1 0.2812488166654865 0.3845137789370807
     2 0.27659193269243015 0.3875453524100341
     3 0.27075461978147 0.3923292814826181
     4 0.2635704413996355 0.3983998843331933
     5 0.25151664507776433 0.4061626772130177
     6 0.23663369143081694 0.41555059289162133
     7 0.21821188439371514 0.42652011369339826
     8 0.196755994839884 0.4387150193246893
     9 0.17308577769384706 0.4518920343703341
     10 0.14901980890265892 0.4657525526156948
     11 0.12446112645939521 0.4796481639615766
     12 0.10319874740737675 0.49550774462815844
     13 0.08587012318951315 0.5129226962978193
     14 0.07292322237902125 0.5296966903033001
     15 0.0616852126807945 0.5451477887017905
     16 0.05204360948282743 0.5584496646561611
     17 0.04419330583874016 0.5767485456188547
     18 0.037619319353642565 0.5937927941163773
     19 0.03367764166025847 0.6103474655281609
 %% Cell type:code id: tags:
 ``` python
 ll.run(4_000)
 np.min(ll.phi[:, :]), np.max(ll.phi[:, :])
 ```
 %% Output
     $\displaystyle \left( 0.0292654136493928, \  0.645196441767244\right)$
     (0.029265413649392794, 0.6451964417672443)
 %% Cell type:code id: tags:
 ``` python
 if 'is_test_run' not in globals():
     plt.phase_plot_for_step(ll)
 ```
 %% Output
 %% Cell type:markdown id: tags:
 ## 2) Advection Test
 %% Cell type:code id: tags:
 ``` python
 sc2 = create_n_phase_model(data_handling=SerialDataHandling((width, 1), periodicity=True),
                            num_phases=3, alpha=alpha)
 ux = 0.05
 phase_idx =1
 sc2.data_handling.fill(sc2.vel_field_name, ux, value_idx=0)
 init_sharp_interface(sc2, phase_idx=phase_idx, inverse=False)
 sc2.set_pdf_fields_from_macroscopic_values()
 plot_status(sc2)
 ```
 %% Output
 %% Cell type:code id: tags:
 ``` python
 if 'is_test_run' in globals():
     sc2.run(1000)
 else:
     sc2.run(10_000)
     assert abs(np.average(sc2.velocity[:, :, 0]) - 0.05) < 1e-10
 plot_status(sc2)
 ```
 %% Output

lbmpy_tests/test_boundary_handling.py

@@ -26,7 +26,7 @@ def mirror_stencil(direction, mirror_axis):
 def test_simple(target):
     if target == Target.GPU:
         import pytest
-        pytest.importorskip('pycuda')
+        pytest.importorskip('cupy')
 
     dh = create_data_handling((4, 4), parallel=False, default_target=target)
     dh.add_array('pdfs', values_per_cell=9, cpu=True, gpu=target != Target.CPU)

lbmpy_tests/test_cpu_gpu_equivalence.py

@@ -31,7 +31,7 @@ def run_equivalence_test(domain_size, lbm_config, lbm_opt, config, time_steps=13
                                       ((18, 20), Method.MRT, True, (4, 2), 'zyxf'),
                                       ((7, 11, 18), Method.TRT, False, False, 'numpy')])
 def test_force_driven_channel_short(scenario):
-    pytest.importorskip("pycuda")
+    pytest.importorskip("cupy")
     ds = scenario[0]
     method = scenario[1]
     compressible = scenario[2]

lbmpy_tests/test_diffusion.py

@@ -77,7 +77,7 @@ def test_diffusion():
 
       The hydrodynamic field is not simulated, instead a constant velocity is assumed.
     """
-    pytest.importorskip("pycuda")
+    pytest.importorskip("cupy")
     # Parameters
     domain_size = (1600, 160)
     omega = 1.38