Central moments

This MR implements the central moment collision operator. Furthermore, some clean-up in the derivation pipeline is provided.

• derivation of the central moments from the continuous equilibrium
• derivation of the central moments via shift matrix for the discrete equilibrium
• implementation of a central moment class
• finish implementation of the collision
• provide test cases for the operator
• add documentation

doc/notebooks/10_tutorial_conservative_allen_cahn_two_phase.ipynb

 %% Cell type:markdown id: tags:
 
-# The conservative Allen-Cahn model for two phase flows with high-denstiy ratios and high Reynolds numbers
+# The conservative Allen-Cahn model for two phase flows with high-density ratios and high Reynolds numbers
 
 %% Cell type:code id: tags:
 
 ``` python
 from pystencils.session import *
 from lbmpy.session import *
 
 from pystencils.simp import sympy_cse
 from pystencils.boundaries import BoundaryHandling
 
 from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle
 from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel, CentralMomentMultiphaseForceModel
 from lbmpy.phasefield_allen_cahn.kernel_equations import *
 from lbmpy.phasefield_allen_cahn.parameter_calculation import calculate_parameters_rti
 
 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
 
 %% Cell type:code id: tags:
 
 ``` python
 try:
     import pycuda
 except ImportError:
     pycuda = None
     gpu = False
     target = 'cpu'
     print('No pycuda installed')
 
 if pycuda:
     gpu = True
     target = 'gpu'
 ```
 
 %% Cell type:markdown id: tags:
 
 The conservative Allen-Cahn model (CACM) for two-phase flows is based on the work of Fakhari et al. [Improved locality of the phase-field lattice-Boltzmann model for immiscible fluids at high density ratios](http://dx.doi.org/10.1103/PhysRevE.96.053301). The model can be created for two-dimensional problems as well as three-dimensional problems, which have been described by Mitchell et al. [Development of a three-dimensional
 phase-field lattice Boltzmann method for the study of immiscible fluids at high density ratios](http://dx.doi.org/10.1103/PhysRevE.96.053301). Furthermore, cascaded lattice Boltzmann methods can be combined with the model which was described in [A cascaded phase-field lattice Boltzmann model for the simulation of incompressible, immiscible fluids with high density contrast](http://dx.doi.org/10.1016/j.camwa.2019.08.018)
 
 
 The CACM is suitable for simulating highly complex tow phase flow problems with high-density rations and high Reynolds numbers. In this tutorial, an overview should be given on how to derive the model with lbmpy. Basically, the model consists of two LBM steps. One for the interface tracking, which we call phase-field LB step and one for recovering the hydrodynamic properties. The latter is called the hydrodynamic LB step.
 
 %% 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. For 2D simulations, only the D2Q9 stencil is supported. 3D simulations support the D3Q15, D3Q19 and the D3Q27 stencil. Note here that the cascaded LBM can not be derived for D3Q15 stencils.
 
 %% Cell type:code id: tags:
 
 ``` python
 stencil_phase = get_stencil("D2Q9")
 stencil_hydro = get_stencil("D2Q9")
 assert(len(stencil_phase[0]) == len(stencil_hydro[0]))
 
 dimensions = len(stencil_phase[0])
 ```
 
 %% Cell type:markdown id: tags:
 
 Definition of the domain size
 
 %% Cell type:code id: tags:
 
 ``` python
 # domain
 L0 = 256
 domain_size = (L0, 4 * L0)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Parameter definition
 
 The next step is to calculate all parameters which are needed for the simulation. In this example, a Rayleigh-Taylor instability test case is set up. The parameter calculation for this setup is already implemented in lbmpy and can be used with the dimensionless parameters which describes the problem.
 
 %% Cell type:code id: tags:
 
 ``` python
 # time step
 timesteps = 8000
 
 # reference time
 reference_time = 4000
 # density of the heavier fluid
 rho_H = 1.0
 
 # calculate the parameters for the RTI
 parameters = calculate_parameters_rti(reference_length=L0,
                                       reference_time=reference_time,
                                       density_heavy=rho_H,
                                       capillary_number=0.44,
                                       reynolds_number=3000,
                                       atwood_number=0.998,
                                       peclet_number=1000,
                                       density_ratio=1000,
                                       viscosity_ratio=100)
 # get the parameters
 rho_L = parameters.get("density_light")
 
 mu_H = parameters.get("dynamic_viscosity_heavy")
 mu_L = parameters.get("dynamic_viscosity_light")
 
 tau_H = parameters.get("relaxation_time_heavy")
 tau_L = parameters.get("relaxation_time_light")
 
 sigma = parameters.get("surface_tension")
 M = parameters.get("mobility")
 gravitational_acceleration = parameters.get("gravitational_acceleration")
 
 
 drho3 = (rho_H - rho_L)/3
 # interface thickness
 W = 5
 # coeffcient related to surface tension
 beta = 12.0 * (sigma/W)
 # coeffcient related to surface tension
 kappa = 1.5 * sigma*W
 # relaxation rate allen cahn (h)
 w_c = 1.0/(0.5 + (3.0 * M))
 ```
 
 %% 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)
 
 # pdf fields. g is used for hydrodynamics and h for the interface tracking
 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)
 
 # velocity field
 u = dh.add_array("u", values_per_cell=dh.dim)
 dh.fill("u", 0.0, ghost_layers=True)
 
 # phase-field
 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:markdown id: tags:
 
 As a next step the relaxation time is stated in a symbolic form. It is calculated via interpolation.
 
 %% Cell type:code id: tags:
 
 ``` python
 # relaxation time and rate
 tau = 0.5 + tau_L + (C.center) * (tau_H - tau_L)
 s8 = 1/(tau)
 
 # density for the whole domain
 rho = rho_L + (C.center) * (rho_H - rho_L)
 
 # body force
 body_force = [0, 0, 0]
 body_force[1] = gravitational_acceleration * rho
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Definition of the lattice Boltzmann methods
 
 %% Cell type:markdown id: tags:
 
 For both LB steps, a weighted orthogonal MRT (WMRT) method is used. It is also possible to change the method to a simpler SRT scheme or a more complicated CLBM scheme. The CLBM scheme can be obtained by commenting in the python snippets in the notebook cells below. Note here that the hydrodynamic LB step is formulated as an incompressible velocity bases LBM. Thus the velocity terms can not be removed from the equilibrium in the central moment space.
 
 %% Cell type:code id: tags:
 
 ``` python
 method_phase = create_lb_method(stencil=stencil_phase, method="mrt", compressible=True, weighted=True,
                                 relaxation_rates=[1, 1, 1, 1])
 
 method_phase.set_first_moment_relaxation_rate(w_c)
 
 # method_phase = create_lb_method(stencil=stencil_phase, method="central_moment", compressible=True,
 #                                 relaxation_rates=[0, w_c, 1, 1, 1], equilibrium_order=4)
 method_phase
 ```
 
 %% Output
 
     <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f8d7a8b3dc0>
 
 %% Cell type:code id: tags:
 
 ``` python
 method_hydro = create_lb_method(stencil=stencil_phase, method="mrt", compressible=False, weighted=True,
                                 relaxation_rates=[s8, 1, 1, 1])
 
 # method_hydro = create_lb_method(stencil=stencil_phase, method="central_moment", compressible=False,
 #                                 relaxation_rates=[s8, 1, 1], equilibrium_order=4)
 
 method_hydro
 ```
 
 %% Output
 
     <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f8d61b4a640>
 
 %% Cell type:markdown id: tags:
 
 ## Initialization
 
 %% Cell type:markdown id: tags:
 
 The initialisation of the pdfs is done with the equilibrium of the LB methods.
 
 %% Cell type:code id: tags:
 
 ``` python
 h_updates = initializer_kernel_phase_field_lb(h, C, u, method_phase, W)
 g_updates = initializer_kernel_hydro_lb(g, u, method_hydro)
 
 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()
 ```
 
 %% Cell type:markdown id: tags:
 
 Afterwards the phase-field gets initialized directly in python
 
 %% Cell type:code id: tags:
 
 ``` python
 # initialize the domain
 def Initialize_distributions():
     Nx = domain_size[0]
     Ny = domain_size[1]
 
     for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
         x = np.zeros_like(block.midpoint_arrays[0])
         x[:, :] = block.midpoint_arrays[0]
 
         y = np.zeros_like(block.midpoint_arrays[1])
         y[:, :] = block.midpoint_arrays[1]
 
         y -= 2 * L0
         tmp = 0.1 * Nx * np.cos((2 * np.pi * x) / Nx)
         init_values = 0.5 + 0.5 * np.tanh((y - tmp) / (W / 2))
         block["C"][:, :] = init_values
         block["C_tmp"][:, :] = init_values
 
     if gpu:
         dh.all_to_gpu()
 
     dh.run_kernel(h_init)
     dh.run_kernel(g_init)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 Initialize_distributions()
 plt.scalar_field(dh.gather_array(C.name))
 ```
 
 %% Output
 
     <matplotlib.image.AxesImage at 0x7f8d617272e0>
 
 
 %% Cell type:markdown id: tags:
 
 ## Source Terms
 
 %% Cell type:markdown id: tags:
 
 For the Allen Cahn LB step, the Allen Cahn equation needs to be applied as a source term. Here a simple forcing model is used which is directly applied in the moment space:
 $$
 F_i^\phi (\boldsymbol{x}, t) = \Delta t \frac{\left[1 - 4 \left(\phi - \phi_0\right)^2\right]}{\xi} w_i \boldsymbol{c}_i \cdot \frac{\nabla \phi}{|{\nabla \phi}|},
 $$
 where $\phi$ is the phase-field, $\phi_0$ is the interface location, $\Delta t$ it the timestep size $\xi$ is the interface width, $\boldsymbol{c}_i$ is the discrete direction from stencil_phase and $w_i$ are the weights. Furthermore, the equilibrium needs to be shifted:
 
 $$
 \bar{h}^{eq}_\alpha = h^{eq}_\alpha - \frac{1}{2} F^\phi_\alpha
 $$
 
 The hydrodynamic force is given by:
 $$
 F_i (\boldsymbol{x}, t) = \Delta t w_i \frac{\boldsymbol{c}_i \boldsymbol{F}}{\rho c_s^2},
 $$
 where $\rho$ is the interpolated density and $\boldsymbol{F}$ is the source term which consists of the pressure force
 $$
 \boldsymbol{F}_p = -p^* c_s^2 \nabla \rho,
 $$
 the surface tension force:
 $$
 \boldsymbol{F}_s = \mu_\phi \nabla \phi
 $$
 and the viscous force term:
 $$
 F_{\mu, i}^{\mathrm{MRT}} = - \frac{\nu}{c_s^2 \Delta t} \left[\sum_{\beta} c_{\beta i} c_{\beta j} \times \sum_{\alpha} \Omega_{\beta \alpha}(g_\alpha - g_\alpha^{\mathrm{eq}})\right] \frac{\partial \rho}{\partial x_j}.
 $$
 
 In the above equations $p^*$ is the normalised pressure which can be obtained from the zeroth order moment of the hydrodynamic distribution function $g$. The lattice speed of sound is given with $c_s$ and the chemical potential is $\mu_\phi$. Furthermore, the viscosity is $\nu$ and $\Omega$ is the moment based collision operator. Note here that the hydrodynamic equilibrium is also adjusted as shown above for the phase-field distribution functions.
 
 
 For CLBM methods the forcing is applied directly in the central moment space. This is done with the `CentralMomentMultiphaseForceModel`. Furthermore, the GUO force model is applied here to be consistent with [A cascaded phase-field lattice Boltzmann model for the simulation of incompressible, immiscible fluids with high density contrast](http://dx.doi.org/10.1016/j.camwa.2019.08.018). Here we refer to equation D.7 which can be derived for 3D stencils automatically with lbmpy.
 
 %% Cell type:code id: tags:
 
 ``` python
 force_h = [f / 3 for f in interface_tracking_force(C, stencil_phase, W)]
 force_g = hydrodynamic_force(g, C, method_hydro, tau, rho_H, rho_L, kappa, beta, body_force)
 
 
 if isinstance(method_phase, CentralMomentBasedLbMethod):
     force_model_h = CentralMomentMultiphaseForceModel(force=force_h)
 else:
     force_model_h = MultiphaseForceModel(force=force_h)
 
 
 if isinstance(method_hydro, CentralMomentBasedLbMethod):
     force_model_g = CentralMomentMultiphaseForceModel(force=force_g, rho=rho)
 else:
     force_model_g = MultiphaseForceModel(force=force_g, rho=rho)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Definition of the LB update rules
 
 %% Cell type:markdown id: tags:
 
 The update rule for the phase-field LB step is defined as:
 
 $$
 h_i (\boldsymbol{x} + \boldsymbol{c}_i \Delta t, t + \Delta t) = h_i(\boldsymbol{x}, t) + \Omega_{ij}^h(\bar{h_j}^{eq} - h_j)|_{(\boldsymbol{x}, t)} + F_i^\phi(\boldsymbol{x}, t).
 $$
 In our framework the pull scheme is applied as streaming step. Furthermore, the update of the phase-field is directly integrated into the kernel. Thus a second temporary phase-field is needed.
 
 %% Cell type:code id: tags:
 
 ``` python
 allen_cahn_lb = get_collision_assignments_phase(lb_method=method_phase,
                                                 velocity_input=u,
                                                 output={'density': C_tmp},
                                                 force_model=force_model_h,
                                                 symbolic_fields={"symbolic_field": h,
                                                                  "symbolic_temporary_field": h_tmp},
                                                 kernel_type='stream_pull_collide')
 
 # allen_cahn_lb = sympy_cse(allen_cahn_lb)
 
 ast_allen_cahn_lb = ps.create_kernel(allen_cahn_lb, target=dh.default_target, cpu_openmp=True)
 kernel_allen_cahn_lb = ast_allen_cahn_lb.compile()
 ```
 
 %% Cell type:markdown id: tags:
 
 The update rule of the hydrodynmaic LB step is defined as:
 
 $$
 g_i (\boldsymbol{x} + \boldsymbol{c}_i \Delta t, t + \Delta t) = g_i(\boldsymbol{x}, t) + \Omega_{ij}^g(\bar{g_j}^{eq} - g_j)|_{(\boldsymbol{x}, t)} + F_i(\boldsymbol{x}, t).
 $$
 
 Here the push scheme is applied which is easier due to the data access of the viscous force term. Furthermore, the velocity update is directly done in the kernel.
 
 %% Cell type:code id: tags:
 
 ``` python
 hydro_lb_update_rule = get_collision_assignments_hydro(lb_method=method_hydro,
                                                        density=rho,
                                                        velocity_input=u,
                                                        force_model=force_model_g,
                                                        sub_iterations=2,
                                                        symbolic_fields={"symbolic_field": g,
                                                                         "symbolic_temporary_field": g_tmp},
                                                        kernel_type='collide_stream_push')
 
 # hydro_lb_update_rule = sympy_cse(hydro_lb_update_rule)
 
 ast_hydro_lb = ps.create_kernel(hydro_lb_update_rule, target=dh.default_target, cpu_openmp=True)
 kernel_hydro_lb = ast_hydro_lb.compile()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Boundary Conditions
 
 %% Cell type:markdown id: tags:
 
 As a last step suitable boundary conditions are applied
 
 %% 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, W)
 
 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)
     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)
     periodic_BC_g()
     bh_hydro()
 
 
     # field swaps
     dh.swap("h", "h_tmp")
     dh.swap("g", "g_tmp")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 Initialize_distributions()
 
 frames = 300
 steps_per_frame = (timesteps//frames) + 1
 
 if 'is_test_run' not in globals():
     def run():
         for i in range(steps_per_frame):
             timeloop()
 
         if gpu:
             dh.to_cpu("C")
         return dh.gather_array(C.name)
 
     animation = plt.scalar_field_animation(run, frames=frames, rescale=True)
     set_display_mode('video')
     res = display_animation(animation)
 else:
     timeloop()
     res = None
 res
 ```
 
 %% Output
 
     <IPython.core.display.HTML object>
 
 %% Cell type:markdown id: tags:
 
 Note that the video is played for 10 seconds while the simulation time is only 2 seconds!