Bump minimum SymPy version and add Python 3.9 to CI

Same as pystencils!208 (merged)

lbmpy_tests/centeredcumulant/test_periodic_pipe_flow.ipynb

 %% Cell type:code id: tags:
 
 ``` python
+import pytest
+pytest.importorskip('pycuda')
 ```
 
+%% Output
 
+    <module 'pycuda' from '/home/markus/miniconda3/envs/pystencils/lib/python3.8/site-packages/pycuda/__init__.py'>
 
+%% Cell type:code id: tags:
 
+``` python
 %load_ext autoreload
 %autoreload 2
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 import sympy as sp
 import numpy as np
 from numpy.testing import assert_allclose
 
 from pystencils.session import *
 from lbmpy.session import *
 
 from lbmpy.stencils import get_stencil
 
 from lbmpy.methods.centeredcumulant import CenteredCumulantForceModel
 
 from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Pipe Flow Scenario
 
 %% Cell type:code id: tags:
 
 ``` python
 class PeriodicPipeFlow:
     def __init__(self, method_params, optimization=dict()):
 
         wall_boundary = NoSlip()
 
         self.stencil = method_params.get('stencil', get_stencil("D2Q9"))
         self.streaming_pattern = method_params.get('streaming_pattern', 'pull')
         self.target = optimization.get('target', 'cpu')
 
         self.gpu = self.target in ['gpu', 'opencl']
 
         #   Stencil
         self.q = len(self.stencil)
         self.dim = len(self.stencil[0])
 
         #   Streaming
         self.inplace = is_inplace(self.streaming_pattern)
         self.timesteps = get_timesteps(self.streaming_pattern)
         self.zeroth_timestep = self.timesteps[0]
 
         #   Domain, Data Handling and PDF fields
         self.pipe_length = 60
         self.pipe_radius = 15
         self.domain_size = (self.pipe_length, ) + (2 * self.pipe_radius,) * (self.dim - 1)
         self.periodicity = (True, ) + (False, ) * (self.dim - 1)
         force = (0.0001, ) + (0.0,) * (self.dim - 1)
         self.force = method_params.get('force', force)
 
         self.dh = create_data_handling(domain_size=self.domain_size,
                                        periodicity=self.periodicity, default_target=self.target)
 
         self.pdfs = self.dh.add_array('pdfs', self.q)
         if not self.inplace:
             self.pdfs_tmp = self.dh.add_array_like('pdfs_tmp', self.pdfs.name)
 
         method_params['force'] = self.force
         optimization['symbolic_field'] = self.pdfs
 
         if not self.inplace:
             optimization['symbolic_temporary_field'] = self.pdfs_tmp
 
         self.lb_collision = create_lb_collision_rule(optimization=optimization, **method_params)
         self.lb_method = self.lb_collision.method
 
         self.lb_kernels = []
         for t in self.timesteps:
             self.lb_kernels.append(create_lb_function(collision_rule=self.lb_collision,
                                                       optimization=optimization,
                                                       timestep=t,
                                                       **method_params))
 
         #   Macroscopic Values
         self.density = 1.0
         self.density_field = self.dh.add_array('rho', 1)
         u_x = 0.0
         self.velocity = (u_x,) * self.dim
         self.velocity_field = self.dh.add_array('u', self.dim)
 
         setter = macroscopic_values_setter(
             self.lb_method, self.density, self.velocity, self.pdfs,
             streaming_pattern=self.streaming_pattern, previous_timestep=self.zeroth_timestep)
         self.init_kernel = create_kernel(setter, ghost_layers=1, target=self.target).compile()
 
         self.getter_kernels = []
         for t in self.timesteps:
             getter = macroscopic_values_getter(
                 self.lb_method, self.density_field, self.velocity_field, self.pdfs,
                 streaming_pattern=self.streaming_pattern, previous_timestep=t)
             self.getter_kernels.append(create_kernel(getter, ghost_layers=1, target=self.target).compile())
 
         #   Periodicity
         self.periodicity_handler = LBMPeriodicityHandling(
             self.stencil, self.dh, self.pdfs.name, streaming_pattern=self.streaming_pattern)
 
         #   Boundary Handling
         self.wall = wall_boundary
         self.bh = LatticeBoltzmannBoundaryHandling(
             self.lb_method, self.dh, self.pdfs.name,
             streaming_pattern=self.streaming_pattern, target=self.target)
 
         self.bh.set_boundary(boundary_obj=self.wall, mask_callback=self.mask_callback)
 
         self.current_timestep = self.zeroth_timestep
 
     def mask_callback(self, x, y, z=None):
         y = y - self.pipe_radius
         z = z - self.pipe_radius if z is not None else 0
         return np.sqrt(y**2 + z**2) >= self.pipe_radius
 
     def init(self):
         self.current_timestep = self.zeroth_timestep
         self.dh.run_kernel(self.init_kernel)
 
     def step(self):
         #   Order matters! First communicate, then boundaries, otherwise
         #   periodicity handling overwrites reflected populations
         # Periodicty
         self.periodicity_handler(self.current_timestep)
 
         # Boundaries
         self.bh(prev_timestep=self.current_timestep)
 
         # Here, the next time step begins
         self.current_timestep = self.current_timestep.next()
 
         # LBM Step
         self.dh.run_kernel(self.lb_kernels[self.current_timestep.idx])
 
         # Field Swaps
         if not self.inplace:
             self.dh.swap(self.pdfs.name, self.pdfs_tmp.name)
 
         # Macroscopic Values
         self.dh.run_kernel(self.getter_kernels[self.current_timestep.idx])
 
     def run(self, iterations):
         for _ in range(iterations):
             self.step()
 
     @property
     def velocity_array(self):
         if self.gpu:
             self.dh.to_cpu(self.velocity_field.name)
         return self.dh.gather_array(self.velocity_field.name)
 
     def get_trimmed_velocity_array(self):
         if self.gpu:
             self.dh.to_cpu(self.velocity_field.name)
         u = np.copy(self.dh.gather_array(self.velocity_field.name))
         mask = self.bh.get_mask(None, self.wall)
         for idx in np.ndindex(u.shape[:-1]):
             if mask[idx] != 0:
                 u[idx] = np.full((self.dim, ), np.nan)
         return u
 ```
 
 %% Cell type:markdown id: tags:
 
 ## General Setup
 
 %% Cell type:code id: tags:
 
 ``` python
 stencil = get_stencil('D3Q19')
 dim = len(stencil[0])
 target = 'gpu'
 streaming_pattern = 'aa'
 
 optimization = { 'target': target }
 
 viscous_rr = 1.1
 force = (0.0001, ) + (0.0,) * (dim - 1)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## 1. Reference: SRT Method
 
 %% Cell type:code id: tags:
 
 ``` python
 srt_params = {
     'stencil': stencil,
     'method': 'srt',
     'relaxation_rate': viscous_rr,
     'force_model': 'guo',
     'force' : force,
     'streaming_pattern': streaming_pattern
 }
 
 srt_flow = PeriodicPipeFlow(srt_params, optimization)
 srt_flow.init()
 srt_flow.run(400)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 srt_u = srt_flow.get_trimmed_velocity_array()
 ps.plot.vector_field_magnitude(srt_u[30,:,:,:])
 ps.plot.colorbar()
 ```
 
 %% Output
 
-    <matplotlib.colorbar.Colorbar at 0x7f8678eccbe0>
+    <matplotlib.colorbar.Colorbar at 0x7f20b2f75bb0>
 
-
-
+
 
 %% Cell type:markdown id: tags:
 
 ## 2. Centered Cumulant Method with Implicit Forcing
 
 The default setup of the centered cumulant method does not use a force model to compute the force contributions on all populations, but instead applies the force in cumulant space simply by setting the momentum relaxation rate $\omega_1 = 2$. Due to the half-force shift of the equilibrium input velocity, the first central moments (which correspond to momentum relative to the moving frame of reference) are not zero, but equal to $-F/2$. The relaxation process causes the first central moments to simply change sign:
 $$
     \kappa_{100}^* = - \kappa_{100}.
 $$
 Thus, $\kappa_{100}^* = F_x /2$. In total, the entire acceleration given by the force is added onto the momentum.
 
 %% Cell type:code id: tags:
 
 ``` python
 cm_method_params = {
     'method' : 'monomial_cumulant',
     'stencil': stencil,
     'relaxation_rate': viscous_rr,  #   Specify viscous relaxation rate only
     'force_model': CenteredCumulantForceModel(force),
     'force' : force,
     'streaming_pattern' : streaming_pattern
 }
 
 optimization = {
     'target': target,
     'pre_simplification' : True
 }
 
 cm_impl_f_flow = PeriodicPipeFlow(cm_method_params, optimization)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 lb_method = cm_impl_f_flow.lb_method
 assert all(rr == 2 for rr in lb_method.relaxation_rates[1:4])
 assert all(rr == viscous_rr for rr in lb_method.relaxation_rates[4:9])
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 cm_impl_f_flow.init()
 cm_impl_f_flow.run(400)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 cm_impl_f_u = cm_impl_f_flow.get_trimmed_velocity_array()
 ps.plot.vector_field_magnitude(cm_impl_f_u[30,:,:,:])
 ps.plot.colorbar()
 ```
 
 %% Output
 
-    <matplotlib.colorbar.Colorbar at 0x7f8678e35130>
+    <matplotlib.colorbar.Colorbar at 0x7f20aa846670>
 
-
-
+
 
 %% Cell type:code id: tags:
 
 ``` python
 assert_allclose(cm_impl_f_u, srt_u, rtol=1, atol=1e-4)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## 3. Centered Cumulant Method with Explicit Forcing
 
 %% Cell type:code id: tags:
 
 ``` python
 from lbmpy.forcemodels import Schiller
 
 cm_method_params_expl_force = {
     'method' : 'cumulant',
     'stencil': stencil,
     'relaxation_rates': [viscous_rr],  #   Specify viscous relaxation rate only
     'force_model': Schiller(force),
     'force' : force,
     'streaming_pattern' : streaming_pattern
 }
 
 optimization = {
     'target': target,
     'pre_simplification' : True
 }
 
 cm_expl_f_flow = PeriodicPipeFlow(cm_method_params_expl_force, optimization)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 lb_method = cm_expl_f_flow.lb_method
 assert all(rr == 0 for rr in lb_method.relaxation_rates[1:4])
 assert all(rr == viscous_rr for rr in lb_method.relaxation_rates[4:9])
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 cm_expl_f_flow.init()
 cm_expl_f_flow.run(400)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 cm_expl_f_u = cm_expl_f_flow.get_trimmed_velocity_array()
 ps.plot.vector_field_magnitude(cm_expl_f_u[30,:,:,:])
 ps.plot.colorbar()
 ```
 
 %% Output
 
-    <matplotlib.colorbar.Colorbar at 0x7f8678b73c10>
+    <matplotlib.colorbar.Colorbar at 0x7f20aa641070>
 
-
-
+
 
 %% Cell type:code id: tags:
 
 ``` python
 assert_allclose(cm_expl_f_u, srt_u, rtol=1, atol=1e-4)
 assert_allclose(cm_expl_f_u, cm_impl_f_u, rtol=1, atol=1e-4)
 ```