.gitignore

@@ -23,6 +23,7 @@ doc/bibtex.json
 
 # macOS
 **/.DS_Store
+*.uuid
 
 # benchmark database
 /lbmpy_tests/benchmark/db
\ No newline at end of file

.gitlab-ci.yml

@@ -22,7 +22,7 @@ tests-and-coverage:
     - pip install git+https://gitlab-ci-token:${CI_JOB_TOKEN}@i10git.cs.fau.de/pycodegen/pystencils.git@master#egg=pystencils
     - env
     - pip list
-    - py.test -v -n $NUM_CORES --cov-report html --cov-report term --cov=. -m "not longrun" --junitxml=report.xml
+    - py.test -v -n $NUM_CORES --cov-report html --cov-report xml --cov-report term --cov=. -m "not longrun" --junitxml=report.xml
     - python3 -m coverage xml
   tags:
     - docker
@@ -113,22 +113,22 @@ latest-python:
       junit: report.xml
 
 # Minimal tests in windows environment
-minimal-windows:
-  stage: test
-  except:
-    variables:
-      - $ENABLE_NIGHTLY_BUILDS
-  tags:
-    - win
-  script:
-    - export NUM_CORES=$(nproc --all)
-    - export MPLBACKEND=Agg
-    - source /cygdrive/c/Users/build/Miniconda3/Scripts/activate
-    - source activate pystencils
-    - pip install git+https://gitlab-ci-token:${CI_JOB_TOKEN}@i10git.cs.fau.de/pycodegen/pystencils.git@master#egg=pystencils
-    - python -c "import numpy"
-    - pip install sympy==1.9
-    - py.test -v -m "not (notebook or longrun)"
+#minimal-windows:
+#  stage: test
+#  except:
+#    variables:
+#      - $ENABLE_NIGHTLY_BUILDS
+#  tags:
+#    - win
+#  script:
+#    - export NUM_CORES=$(nproc --all)
+#    - export MPLBACKEND=Agg
+#    - source /cygdrive/c/Users/build/Miniconda3/Scripts/activate
+#    - source activate pystencils
+#    - pip install git+https://gitlab-ci-token:${CI_JOB_TOKEN}@i10git.cs.fau.de/pycodegen/pystencils.git@master#egg=pystencils
+#    - python -c "import numpy"
+#    - pip install sympy==1.9
+#    - py.test -v -m "not (notebook or longrun)"
 
 minimal-sympy-master:
   stage: test

README.md

@@ -71,6 +71,7 @@ Many thanks go to the [contributors](https://i10git.cs.fau.de/pycodegen/lbmpy/-/
 If you use lbmpy in a publication, please cite the following articles:
 
 Overview:
+  - F. Hennig et al, Advanced Automatic Code Generation for Multiple Relaxation-Time Lattice Boltzmann Methods. SIAM Journal on Scientific Computing, 2023. https://doi.org/10.1137/22M1531348 ([Preprint](https://arxiv.org/abs/2211.02435))
   - M. Bauer et al, lbmpy: Automatic code generation for efficient parallel lattice Boltzmann methods. Journal of Computational Science, 2021. https://doi.org/10.1016/j.jocs.2020.101269 ([Preprint](https://arxiv.org/abs/2001.11806))
 
 Multiphase:

binder/environment.yml

@@ -7,7 +7,7 @@
 #     conda env create -f conda_environment_user.yml
 #     . activate pystencils
 #
-# If you have CUDA installed and want to use your GPU, uncomment the last line to install pycuda
+# If you have CUDA or ROCm installed and want to use your GPU, uncomment the last line to install cupy
 #
 # ----------------------------------------------------------------------------------------------------------------------
 
@@ -33,4 +33,4 @@ dependencies:
       - pyevtk # VTK output for serial simulations
       - blitzdb # file-based No-SQL database to store simulation results
       - pystencils
-      #- pycuda # add this if you have CUDA installed
+      #- cupy # add this if you have CUDA or ROCm installed

conftest.py

@@ -46,7 +46,7 @@ add_path_to_ignore('pystencils_tests/benchmark')
 add_path_to_ignore('_local_tmp')
 
 try:
-    import pycuda
+    import cupy
 except ImportError:
     collect_ignore += [os.path.join(SCRIPT_FOLDER, "lbmpy_tests/test_cpu_gpu_equivalence.py")]
 
@@ -60,7 +60,10 @@ try:
 except ImportError:
     collect_ignore += [os.path.join(SCRIPT_FOLDER, "lbmpy_tests/benchmark"),
                        os.path.join(SCRIPT_FOLDER,
-                                    "lbmpy_tests/full_scenarios/kida_vortex_flow/scenario_kida_vortex_flow.py")]
+                                    "lbmpy_tests/full_scenarios/kida_vortex_flow/scenario_kida_vortex_flow.py"),
+                       os.path.join(SCRIPT_FOLDER, "lbmpy_tests/full_scenarios/shear_wave/scenario_shear_wave.py"),
+                       os.path.join(SCRIPT_FOLDER, "lbmpy_tests/test_json_serializer.py"),
+                       os.path.join(SCRIPT_FOLDER, "lbmpy/db.py")]
 
 if platform.system().lower() == 'windows':
     collect_ignore += [os.path.join(SCRIPT_FOLDER, "lbmpy_tests/test_quicktests.py")]

doc/conf.py

@@ -33,7 +33,7 @@ version = re.sub(r'(\d+\.\d+)\.\d+(.*)', r'\1\2', lbmpy.__version__)
 version = re.sub(r'(\.dev\d+).*?$', r'\1', version)
 # The full version, including alpha/beta/rc tags.
 release = lbmpy.__version__
-language = None
+language = 'en'
 exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store', '**.ipynb_checkpoints']
 default_role = 'any'
 pygments_style = 'sphinx'

doc/notebooks/01_tutorial_predefinedScenarios.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 lbmpy.session import *
 from pystencils import Target
 ```
 %% Cell type:markdown id: tags:
 # Tutorial 01: Running pre-defined scenarios
 *lbmpy* is a module to do Lattice Boltzmann simulations in Python.
 In this tutorial you will get a broad overview of *lbmpy*'s features. We will run some of the included scenarios that come with *lbmpy*, like a channel flow and a lid driven cavity. This tutorial uses the simple, high-level API of *lbmpy*, while the following tutorials go into the low-level details.
 The only prerequisite for this tutorial is basic Python and [numpy](http://www.numpy.org/) knowledge.
 > #### What's special about *lbmpy* ?
 > The LBM kernels (i.e. the functions that do all the computations) are not written in Python. Instead *lbmpy* generates optimized C or CUDA code for these kernels and compiles it using the *pystencils* module. In that way we get very fast LBM kernels, a lot faster than pure Python implementations and probably also faster than handwritten C kernels. This sounds complicated, but we don't have to care about all this background work, since all compiled kernels are available as Python functions again. Thus *lbmpy* can be used just like any other Python package.
 ## Lid Driven Cavity
 We start by simulating a fluid in a rectangular box, where one wall (the lid) is moving. This is called a 'lid driven cavity'. At the stationary walls *no-slip* boundary conditions are set, which enforce zero velocity at the wall. At the lid there is a *velocity bounce back (UBB)* boundary condition, which sets zero normal velocity and a prescribed tangential velocity.
 We don't have to set up all these boundary conditions manually since there is a function ``create_lid_driven_cavity``  that does all the work for us. This function takes the tangential velocity of the lid, which drives the flow. It is given in lattice units and to get a stable simulation it should be smaller than 0.1. The `relaxation_rate` determines the viscosity of the fluid: Small relaxation rates correspond to high viscosity. The `relaxation_rate` has to be between 0 and 2.
 %% Cell type:code id: tags:
 ``` python
 ldc_scenario = create_lid_driven_cavity(domain_size=(80,50), lid_velocity=0.01, relaxation_rate=1.95)
 ldc_scenario.method
 ```
 %% Output
     <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7fe7a8828d30>
 %% Cell type:markdown id: tags:
 The *run* method of the scenario runs the specified amount of time steps. When you run the next cell, 2000 time steps are executed and the velocity field is plotted. You can run the cell multiple times to see a time evolution.
 %% Cell type:code id: tags:
 ``` python
 ldc_scenario.run(2000)
 plt.figure(dpi=200)
 plt.vector_field(ldc_scenario.velocity_slice(), step=2);
 ```
 %% Output
 %% Cell type:markdown id: tags:
 ### Variations to experiment with:
 - simulate with a higher ``relaxation_rate`` (i.e. higher [Reynolds number](https://en.wikipedia.org/wiki/Reynolds_number)), keep in mind that the ``relaxation_rate`` has to be smaller than 2. You might have to increase (``domain_size``) to keep the simulation stable and run more time steps to get to the stationary solution. You also might want to increase the ``step`` parameter for the plot, to reduce the number of arrows.
 - run a 3D simulation by adding a third dimension size to ``domain_size``. The ``velocity`` property of the scenario is now a 3D field that has to be sliced before it can be plotted, e.g. ``ldc_scenario.velocity[:, :, 10, 0:2]`` generates a slice at ``z=10`` and plot the ``x`` and ``y`` component of the velocity.
 %% Cell type:markdown id: tags:
 ##  Fully periodic flow
 Another simple scenario is a box with periodic boundary conditions in all directions. We initialize a non-zero initial velocity field, which is decaying over time due to viscous effects and the absence of driving forces or boundary conditions. In this example we initialize a shear flow where in one stripe the fluid is moving to the left, and everywhere else to the right. We perturbe this initial velocity field with random noise to get an instable shear layer.
 %% Cell type:code id: tags:
 ``` python
 width, height = 200, 60
 velocity_magnitude = 0.05
 init_vel = np.zeros((width,height,2))
 # fluid moving to the right everywhere...
 init_vel[:, :, 0] = velocity_magnitude
  # ...except at a stripe in the middle, where it moves left
 init_vel[:, height//3 : height//3*2, 0] = -velocity_magnitude
 # small random y velocity component
 init_vel[:, :, 1] = 0.1 * velocity_magnitude * np.random.rand(width,height)
 plt.figure(dpi=200)
 plt.vector_field(init_vel, step=4);
 ```
 %% Output
 %% Cell type:markdown id: tags:
 With this initial velocity field we create a simulation scenario:
 %% Cell type:code id: tags:
 ``` python
 shear_flow_scenario = create_fully_periodic_flow(initial_velocity=init_vel, relaxation_rate=1.97)
 ```
 %% Cell type:code id: tags:
 ``` python
 shear_flow_scenario.run(500)
 plt.figure(dpi=200)
 plt.vector_field(shear_flow_scenario.velocity[:, :])
 ```
 %% Output
     <matplotlib.quiver.Quiver at 0x7fe7a80e2bb0>
 %% Cell type:markdown id: tags:
 Instead of plotting a single point in time we create an animation. For this we first have to
 define an update function that runs a few time steps and returns the field to plot.
 This function is called ``iterations`` times, then the animation stops.
 To cancel the animation while it is running, hit the stop button in the IPython menu bar.
 %% Cell type:code id: tags:
 ``` python
 # def next_frame():
 #     shear_flow_scenario.run(50)
 #     return shear_flow_scenario.velocity[:, :]
 # plt.figure(dpi=200)
 # display_animation(plt.vector_field_animation(next_frame, step=2), iterations=50)
 ```
 %% Cell type:markdown id: tags:
 Vortices are created between the two layers. This phenomenon is called [Kelvin-Helmholz Instability](https://en.wikipedia.org/wiki/Kelvin%E2%80%93Helmholtz_instability). For a better visualization of the vortices we can plot the [vorticity](https://en.wikipedia.org/wiki/Vorticity) of velocity field:
 %% Cell type:code id: tags:
 ``` python
 plt.figure(dpi=200)
 plt.scalar_field(vorticity_2d(shear_flow_scenario.velocity[:, :]));
 ```
 %% Output
 %% Cell type:markdown id: tags:
 ###  Variations to experiment with:
 - make an animation of the vorticity
 - increase the ``relaxation_rate``. What is the maximum relaxation rate you can get before the simulation gets unstable?
 - use an entropic method to get to higher relaxation rates:
 ```
 entropic_shear_flow_scenario = create_fully_periodic_flow(initial_velocity=init_vel, method='trt-kbc-n4',
                                                           entropic=True, compressible=True)
 entropic_shear_flow_scenario.kernel_params['omega_0'] = 1.999
 ```
 %% Cell type:markdown id: tags:
 ## Channel
 In the last part of this tutorial you learn how to modify the boundary handling of scenarios.
 Therefor we set up a channel flow and place some objects into it.
 The channel will be driven by a constant body force e.g. gravity which acts in x direction. Along the flow direction periodic boundary conditions are used whereas the walls are modeled with a *noslip* boundary condition.
 %% Cell type:code id: tags:
 ``` python
 channel_scenario = create_channel(domain_size=(300, 100), force=1e-7, initial_velocity=(0.025, 0),
                                   relaxation_rate=1.97,
                                   config=CreateKernelConfig(target=Target.GPU))
 ```
 %% Cell type:code id: tags:
 ``` python
 channel_scenario._lbmKernels[0].ast
 ```
 %% Output
     KernelFunction kernel([_data_force_driven_channel_pdfSrc, _data_force_driven_channel_pdfTmp])
 %% Cell type:markdown id: tags:
 As in the last scenario, we specify an initial velocity here. Instead of passing a velocity value for every cell we specify here a constant for the complete domain.
 %% Cell type:code id: tags:
 ``` python
 channel_scenario.run(10000)
 plt.figure(dpi=200)
 plt.vector_field(channel_scenario.velocity[:, :], step=4);
 ```
 %% Output
 %% Cell type:markdown id: tags:
 This is a 2D [Poiseuille flow](https://en.wikipedia.org/wiki/Hagen%E2%80%93Poiseuille_equation) where a parabolic profile of the x velocity is expected for the stationary case.
 %% Cell type:code id: tags:
 ``` python
 vel_profile = channel_scenario.velocity[0.5, :, 0]
 plt.figure(dpi=200)
 plt.plot(vel_profile);
 ```
 %% Output
 %% Cell type:markdown id: tags:
 The stationary state is not yet reached, you can run more time steps and see how the profile gets closer to a parabola.
 ### Modifying boundaries
 Lets first view the current boundary configuration:
 %% Cell type:code id: tags:
 ``` python
 def draw_boundary_setup():
     fig = plt.figure(figsize=(10.0, 3.0), dpi=200)
     plt.boundary_handling(channel_scenario.boundary_handling)
     plt.axis('off');
 draw_boundary_setup()
 ```
 %% Output
 %% Cell type:markdown id: tags:
 In above plot you can see that the *no-slip* boundaries at the top and bottom of the channel, otherwise the channel is empty.
 Since an empty channel is pretty boring, lets put an obstacle in it. We start with the simplest option: a rectangular, solid block. The rectangle is specified as a slice, similar to advanced *numpy* indexing. The ``make_slice`` function can also take ``float`` as indices for specifying the slice relative to the domain size. The following cell puts an obstacle into the domain that has one third of the channel height.
 Additionally we have to pass a function that defines what should happen at the boundary (here ``noSlip``).
 By default, the name of the function is also the boundary name.
 %% Cell type:code id: tags:
 ``` python
 from lbmpy.boundaries import NoSlip
 wall = NoSlip()
 channel_scenario.boundary_handling.set_boundary(wall, make_slice[0.2:0.25, 0:0.333])
 ```
 %% Output
     2
 %% Cell type:markdown id: tags:
 When plotting the boundary handling again, we see the rectangular obstacle:
 %% Cell type:code id: tags:
 ``` python
 draw_boundary_setup()
 ```
 %% Output
 %% Cell type:markdown id: tags:
 When setting and plotting boundaries the domain is actually 2 cells larger than originally specified. These so called 'ghost layer' slices are one cell thick and are automatically added at each domain boundary. They can be used to set boundaries and are also used for communication when running distributed memory parallel simulations.
 When specifying the slice, keep in mind that the domain is actually slightly larger. When plotting the simulation results the ghost layers are automatically removed.
 To convert a cell back to ``domain``, the same method can be used. To demonstrate this, we cut a piece out of the obstacle:
 %% Cell type:code id: tags:
 ``` python
 channel_scenario.boundary_handling.set_boundary('domain', make_slice[0.2:0.235, 0.0333:0.3])
 fig = plt.figure(figsize=(10.0, 3.0), dpi=200)
 plt.boundary_handling(channel_scenario.boundary_handling)
 plt.axis('off');
 ```
 %% Output
 %% Cell type:markdown id: tags:
 To add non-rectangular obstacles one can also pass a mask array, which is ``True`` for cells where the boundary should be set. This is demonstrated in the next cell, where a sphere is placed in the channel.
 %% Cell type:code id: tags:
 ``` python
 def set_sphere(x, y):
     shape = channel_scenario.domain_size
     mid = (0.5 * shape[0], 0.5 * shape[1])
     radius = 13
     return (x-mid[0])**2 + (y-mid[1])**2 < radius**2
 channel_scenario.boundary_handling.set_boundary(wall, mask_callback=set_sphere)
 draw_boundary_setup()
 ```
 %% Output
 %% Cell type:markdown id: tags:
 Now, we run the simulation to see the flow aroung the obstacles:
 %% Cell type:code id: tags:
 ``` python
 channel_scenario.run(10000)
 plt.figure(dpi=200)
 plt.vector_field_magnitude(channel_scenario.velocity[:,:]);
 ```
 %% Output
 %% Cell type:markdown id: tags:
 ###  Variations to experiment with:
 - increase the Reynolds number. You might also have to increase the resolution, and/or use a more advanced method like cumulant or entropic stabilization

doc/notebooks/10_tutorial_conservative_allen_cahn_two_phase.ipynb

 %% Cell type:markdown id: tags:
 
 # The conservative Allen-Cahn model for high Reynolds number, two phase flow with large-density and viscosity constrast
 
 %% 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.kernel_equations import *
 from lbmpy.phasefield_allen_cahn.parameter_calculation import calculate_parameters_rti, AllenCahnParameters
 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:
 
 The conservative Allen-Cahn model (CACM) for two-phase flow is based on the work of Fakhari et al. (2017) [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. (2018) [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 two phase flow problems with high-density ratios and high Reynolds numbers. In this tutorial, an overview is provided on how to derive the model with lbmpy. For this, the model is defined with two LBM populations. One for the interface tracking, which we call the 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 can be performed in either 2D- or 3D-space. For 2D simulations, only the D2Q9 stencil is supported. For 3D simulations, the D3Q15, D3Q19 and the D3Q27 stencil are supported. Note here that the cascaded LBM can not be derived for D3Q15 stencils.
 
 %% Cell type:code id: tags:
 
 ``` python
 stencil_phase = LBStencil(Stencil.D2Q9)
 stencil_hydro = LBStencil(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 describe the problem.
 
 %% Cell type:code id: tags:
 
 ``` python
 # time step
 timesteps = 8000
 
 # reference time
 reference_time = 4000
 
 # calculate the parameters for the RTI
 parameters = calculate_parameters_rti(reference_length=L0,
                                       reference_time=reference_time,
                                       density_heavy=1.0,
                                       capillary_number=0.44,
                                       reynolds_number=3000,
                                       atwood_number=0.998,
                                       peclet_number=1000,
                                       density_ratio=1000,
                                       viscosity_ratio=100)
 ```
 
 %% Cell type:markdown id: tags:
 
 This function returns a `AllenCahnParameters` class. It is struct like class holding all parameters for the conservative Allen Cahn model:
 
 %% Cell type:code id: tags:
 
 ``` python
 parameters
 ```
 
 %% Output
 
     <lbmpy.phasefield_allen_cahn.parameter_calculation.AllenCahnParameters at 0x126d30cd0>
 
 %% 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/). This object 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
 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: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 with `Method.CENTRAL_MOMENT`. Note here that the hydrodynamic LB step is formulated as an incompressible velocity-based LBM. Thus, the velocity terms can not be removed from the equilibrium in the central moment space.
 
 %% Cell type:code id: tags:
 
 ``` python
 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)
 
 method_phase
 ```
 
 %% Output
 
     <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x126d44ee0>
 
 %% 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)
 
 method_hydro
 ```
 
 %% Output
 
     <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x126d22eb0>
 
 %% Cell type:markdown id: tags:
 
 ## Initialization
 
 %% Cell type:markdown id: tags:
 
 The probability distribution functions (pdfs) are initialised with the equilibrium distribution for the LB methods.
 
 %% 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()
 ```
 
 %% Cell type:markdown id: tags:
 
 Following this, the phase field is initialised 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) / (parameters.interface_thickness / 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
 Initialize_distributions()
 plt.scalar_field(dh.gather_array(C.name))
 ```
 
 %% Output
 
     <matplotlib.image.AxesImage at 0x1417f86d0>
 
 
 %% 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 = 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: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. As a result of this, a second temporary phase-field is needed.
 
 %% 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()
 ```
 
 %% Cell type:markdown id: tags:
 
 The update rule for 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 required for the viscous force term. Furthermore, the velocity update is directly done in the kernel.
 
 %% 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()
 ```
 
 %% 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, parameters.interface_thickness)
 
 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()
 
 
     # 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!

doc/sphinx/lbmpy.bib

@@ -110,6 +110,24 @@ issn = {0898-1221},
 doi = {10.1016/j.camwa.2015.05.001},
 }
 
+@article{geier2017,
+author = {Geier, Martin and Pasquali, Andrea and Sch{\"{o}}nherr, Martin},
+title = {Parametrization of the Cumulant Lattice Boltzmann Method for Fourth Order Accurate Diffusion Part I},
+year = {2017},
+issue_date = {November 2017},
+publisher = {Academic Press Professional, Inc.},
+address = {USA},
+volume = {348},
+number = {C},
+issn = {0021-9991},
+url = {https://doi.org/10.1016/j.jcp.2017.05.040},
+doi = {10.1016/j.jcp.2017.05.040},
+journal = {J. Comput. Phys.},
+month = {nov},
+pages = {862–888},
+numpages = {27}
+}
+
 @Article{Coreixas2019,
   title = {Comprehensive comparison of collision models in the lattice Boltzmann framework: Theoretical investigations},
   author = {Coreixas, Christophe and Chopard, Bastien and Latt, Jonas},

lbmpy/advanced_streaming/communication.py

@@ -12,15 +12,15 @@ class LBMPeriodicityHandling:
 
     def __init__(self, stencil, data_handling, pdf_field_name,
                  streaming_pattern='pull', ghost_layers=1,
-                 pycuda_direct_copy=True):
+                 cupy_direct_copy=True):
         """
             Periodicity Handling for Lattice Boltzmann Streaming.
 
             **On the usage with cuda:**
-            - pycuda allows the copying of sliced arrays within device memory using the numpy syntax,
+            - cupy allows the copying of sliced arrays within device memory using the numpy syntax,
             e.g. `dst[:,0] = src[:,-1]`. In this implementation, this is the default for periodicity
-            handling. Alternatively, if you set `pycuda_direct_copy=False`, GPU kernels are generated and
-            compiled. The compiled kernels are almost twice as fast in execution as pycuda array copying,
+            handling. Alternatively, if you set `cupy_direct_copy=False`, GPU kernels are generated and
+            compiled. The compiled kernels are almost twice as fast in execution as cupy array copying,
             but especially for large stencils like D3Q27, their compilation can take up to 20 seconds.
             Choose your weapon depending on your use case.
         """
@@ -40,7 +40,7 @@ class LBMPeriodicityHandling:
         self.inplace_pattern = is_inplace(streaming_pattern)
 
         self.cpu = self.target == Target.CPU
-        self.pycuda_direct_copy = self.target == Target.GPU and pycuda_direct_copy
+        self.cupy_direct_copy = self.target == Target.GPU and cupy_direct_copy
 
         def is_copy_direction(direction):
             s = 0
@@ -63,7 +63,7 @@ class LBMPeriodicityHandling:
                                                            ghost_layers=ghost_layers)
             self.comm_slices.append(list(chain.from_iterable(v for k, v in slices_per_comm_dir.items())))
 
-        if self.target == Target.GPU and not pycuda_direct_copy:
+        if self.target == Target.GPU and not cupy_direct_copy:
             self.device_copy_kernels = list()
             for timestep in timesteps:
                 self.device_copy_kernels.append(self._compile_copy_kernels(timestep))
@@ -90,7 +90,7 @@ class LBMPeriodicityHandling:
 
     def _periodicity_handling_gpu(self, prev_timestep):
         arr = self.dh.gpu_arrays[self.pdf_field_name]
-        if self.pycuda_direct_copy:
+        if self.cupy_direct_copy:
             for src, dst in self.comm_slices[prev_timestep.idx]:
                 arr[dst] = arr[src]
         else:

lbmpy/chapman_enskog/chapman_enskog.py

@@ -370,7 +370,7 @@ def take_moments(eqn, pdf_to_moment_name=(('f', '\\Pi'), ('\\Omega f', '\\Upsilo
                 if new_f_index is None:
                     rest *= factor
                 else:
-                    assert not(new_f_index and f_index)
+                    assert not (new_f_index and f_index)
                     f_index = new_f_index
 
         moment_tuple = [0] * len(velocity_terms)

lbmpy/creationfunctions.py

@@ -87,7 +87,7 @@ from pystencils.simp import sympy_cse, SimplificationStrategy
 from lbmpy.methods.abstractlbmethod import LbmCollisionRule, AbstractLbMethod
 from lbmpy.methods.cumulantbased import CumulantBasedLbMethod
 
-# Filter out JobLib warnings. They are not usefull for use:
+# Filter out JobLib warnings. They are not useful for use:
 # https://github.com/joblib/joblib/issues/683
 filterwarnings("ignore", message="Persisting input arguments took")
 
@@ -99,7 +99,7 @@ class LBMConfig:
     """
     stencil: lbmpy.stencils.LBStencil = LBStencil(Stencil.D2Q9)
     """
-    All stencils are defined in :class:`lbmpy.enums.Stencil`. From that :class:`lbmpy.stencils.LBStenil` 
+    All stencils are defined in :class:`lbmpy.enums.Stencil`. From that :class:`lbmpy.stencils.LBStencil` 
     class will be created
     """
     method: Method = Method.SRT
@@ -197,6 +197,13 @@ class LBMConfig:
     Special correction for D3Q27 cumulant LBMs. For Details see
     :mod:`lbmpy.methods.cumulantbased.galilean_correction`
     """
+    fourth_order_correction: Union[float, bool] = False
+    """
+    Special correction for rendering D3Q27 cumulant LBMs fourth-order accurate in diffusion. For Details see
+    :mod:`lbmpy.methods.cumulantbased.fourth_order_correction`. If set to `True`, the fourth-order correction is 
+    employed without limiters (or more precisely with a very high limiter, practically disabling the limiters). If this 
+    variable is set to a number, the latter is used for the limiters (uniformly for omega_3, omega_4 and omega_5). 
+    """
     collision_space_info: CollisionSpaceInfo = None
     """
     Information about the LB method's collision space (see :class:`lbmpy.methods.creationfunctions.CollisionSpaceInfo`)
@@ -412,7 +419,10 @@ class LBMConfig:
                 force_not_zero = True
 
         if self.force_model is None and force_not_zero:
-            self.force_model = forcemodels.Guo(self.force[:self.stencil.D])
+            if self.method == Method.CUMULANT:
+                self.force_model = forcemodels.CentralMoment(self.force[:self.stencil.D])
+            else:
+                self.force_model = forcemodels.Guo(self.force[:self.stencil.D])
 
         force_model_dict = {
             'simple': forcemodels.Simple,
@@ -424,7 +434,8 @@ class LBMConfig:
             'edm': forcemodels.EDM,
             'kupershtokh': forcemodels.EDM,
             'he': forcemodels.He,
-            'shanchen': forcemodels.ShanChen
+            'shanchen': forcemodels.ShanChen,
+            'centralmoment': forcemodels.CentralMoment
         }
 
         if isinstance(self.force_model, str):
@@ -642,6 +653,19 @@ def create_lb_collision_rule(lb_method=None, lbm_config=None, lbm_optimisation=N
         from lbmpy.methods.cumulantbased import add_galilean_correction
         collision_rule = add_galilean_correction(collision_rule)
 
+    if lbm_config.fourth_order_correction:
+        from lbmpy.methods.cumulantbased import add_fourth_order_correction
+
+        # must provide a second relaxation rate in implementation; defaults to 1
+        if len(lbm_config.relaxation_rates) == 1:
+            lbm_config.relaxation_rates.append(1)
+
+        cumulant_limiter = 1e6 if lbm_config.fourth_order_correction is True else lbm_config.fourth_order_correction
+        collision_rule = add_fourth_order_correction(collision_rule=collision_rule,
+                                                     shear_relaxation_rate=lbm_config.relaxation_rates[0],
+                                                     bulk_relaxation_rate=lbm_config.relaxation_rates[1],
+                                                     limiter=cumulant_limiter)
+
     if lbm_config.entropic:
         if lbm_config.smagorinsky or lbm_config.cassons:
             raise ValueError("Choose either entropic, smagorinsky or cassons")
@@ -748,11 +772,25 @@ def create_lb_method(lbm_config=None, **params):
         method_nr = lbm_config.method.name[-1]
         method = create_trt_kbc(dim, relaxation_rates[0], relaxation_rates[1], 'KBC-N' + method_nr, **common_params)
     elif lbm_config.method == Method.CUMULANT:
+        if lbm_config.fourth_order_correction:
+            if lbm_config.stencil.D != 3 and lbm_config.stencil.Q != 27:
+                raise ValueError("Fourth-order correction can only be applied to D3Q27 cumulant methods.")
+
+            assert len(relaxation_rates) <= 2, "Optimal parametrisation for fourth-order cumulants needs either one " \
+                                               "or two relaxation rates, associated with the shear (and bulk) " \
+                                               "viscosity. All other relaxation rates are automatically chosen " \
+                                               "optimally"
+
+            # define method in terms of symbolic relaxation rates and assign optimal values later
+            from lbmpy.methods.cumulantbased.fourth_order_correction import FOURTH_ORDER_RELAXATION_RATE_SYMBOLS
+            relaxation_rates = FOURTH_ORDER_RELAXATION_RATE_SYMBOLS
+
         if lbm_config.nested_moments is not None:
             method = create_cumulant(
                 lbm_config.stencil, relaxation_rates, lbm_config.nested_moments, **cumulant_params)
         else:
             method = create_with_default_polynomial_cumulants(lbm_config.stencil, relaxation_rates, **cumulant_params)
+
     elif lbm_config.method == Method.MONOMIAL_CUMULANT:
         method = create_with_monomial_cumulants(lbm_config.stencil, relaxation_rates, **cumulant_params)
     else:

lbmpy/db.py

+import json
+import six
+import inspect
+
+from pystencils.runhelper.db import PystencilsJsonEncoder
+from pystencils.simp import SimplificationStrategy
+from lbmpy import LBStencil, Method, CollisionSpace
+from lbmpy.creationfunctions import LBMConfig, LBMOptimisation
+from lbmpy.methods import CollisionSpaceInfo
+from lbmpy.forcemodels import AbstractForceModel
+from lbmpy.non_newtonian_models import CassonsParameters
+
+
+class LbmpyJsonEncoder(PystencilsJsonEncoder):
+
+    def default(self, obj):
+        if isinstance(obj, (LBMConfig, LBMOptimisation, CollisionSpaceInfo, CassonsParameters)):
+            return obj.__dict__
+        if isinstance(obj, (LBStencil, Method, CollisionSpace)):
+            return obj.name
+        if isinstance(obj, AbstractForceModel):
+            return obj.__class__.__name__
+        if isinstance(obj, SimplificationStrategy):
+            return obj.__str__()
+        if inspect.isclass(obj):
+            return obj.__name__
+        return PystencilsJsonEncoder.default(self, obj)
+
+
+class LbmpyJsonSerializer(object):
+
+    @classmethod
+    def serialize(cls, data):
+        if six.PY3:
+            if isinstance(data, bytes):
+                return json.dumps(data.decode('utf-8'), cls=LbmpyJsonEncoder, ensure_ascii=False).encode('utf-8')
+            else:
+                return json.dumps(data, cls=LbmpyJsonEncoder, ensure_ascii=False).encode('utf-8')
+        else:
+            return json.dumps(data, cls=LbmpyJsonEncoder, ensure_ascii=False).encode('utf-8')
+
+    @classmethod
+    def deserialize(cls, data):
+        if six.PY3:
+            return json.loads(data.decode('utf-8'))
+        else:
+            return json.loads(data.decode('utf-8'))

lbmpy/enums.py

@@ -217,3 +217,7 @@ class ForceModel(Enum):
     """
     See :class:`lbmpy.forcemodels.ShanChen`
     """
+    CENTRALMOMENT = auto()
+    """
+    See :class:`lbmpy.forcemodels`
+    """

lbmpy/forcemodels.py

@@ -210,6 +210,28 @@ class Simple(AbstractForceModel):
         return before, after
 
 
+class CentralMoment(AbstractForceModel):
+    r"""
+    A force model that only shifts the macroscopic and equilibrium velocities and does not introduce a forcing term to
+    the collision process. Forcing is then applied through relaxation of the first central moments in the shifted
+    frame of reference (cf. https://doi.org/10.1016/j.camwa.2015.05.001).
+    """
+    def __call__(self, lb_method):
+        raise ValueError("This force model can only be combined with the Cumulant collision model")
+
+    def symmetric_central_moment_forcing(self, lb_method, *_):
+        before = sp.Matrix([0, ] * lb_method.stencil.Q)
+        after = sp.Matrix([0, ] * lb_method.stencil.Q)
+        for i, sf in enumerate(self.symbolic_force_vector):
+            before[i + 1] = sp.Rational(1, 2) * sf
+            after[i + 1] = sp.Rational(1, 2) * sf
+
+        return before, after
+
+    def equilibrium_velocity_shift(self, density):
+        return default_velocity_shift(density, self.symbolic_force_vector)
+
+
 class Luo(AbstractForceModel):
     r"""Force model by Luo :cite:`luo1993lattice`.
 

lbmpy/max_domain_size_info.py

@@ -27,16 +27,16 @@ import warnings
 
 import numpy as np
 
-# Optional packages cpuinfo, pycuda and psutil for hardware queries
+# Optional packages cpuinfo, cupy and psutil for hardware queries
 try:
     from cpuinfo import get_cpu_info
 except ImportError:
     get_cpu_info = None
 
 try:
-    from pycuda.autoinit import device
+    import cupy
 except ImportError:
-    device = None
+    cupy = None
 
 try:
     from psutil import virtual_memory
@@ -113,9 +113,11 @@ def memory_sizes_of_current_machine():
         if 'l3_cache_size' in cpu_info:
             result['L3'] = cpu_info['l3_cache_size']
 
-    if device:
-        size = device.total_memory() / (1024 * 1024)
-        result['GPU'] = "{0:.0f} MB".format(size)
+    if cupy:
+        for i in range(cupy.cuda.runtime.getDeviceCount()):
+            device = cupy.cuda.Device(i)
+            size = device.mem_info[1] / (1024 * 1024)
+            result[f'GPU{i}'] = "{0:.0f} MB".format(size)
 
     if virtual_memory:
         mem = virtual_memory()
@@ -124,7 +126,7 @@ def memory_sizes_of_current_machine():
 
     if not result:
         warnings.warn("Couldn't query for any local memory size."
-                      "Install py-cpuinfo to get cache sizes, psutil for RAM size and pycuda for GPU memory size.")
+                      "Install py-cpuinfo to get cache sizes, psutil for RAM size and cupy for GPU memory size.")
 
     return result
 

lbmpy/methods/abstractlbmethod.py

@@ -67,8 +67,8 @@ class AbstractLbMethod(abc.ABC):
         return d
 
     @property
-    def subs_dict_relxation_rate(self):
-        """returns a dictonary which maps the replaced numerical relaxation rates to its original numerical value"""
+    def subs_dict_relaxation_rate(self):
+        """returns a dictionary which maps the replaced numerical relaxation rates to its original numerical value"""
         result = dict()
         for i in range(self._stencil.Q):
             result[self.symbolic_relaxation_matrix[i, i]] = self.relaxation_matrix[i, i]

lbmpy/methods/cumulantbased/__init__.py

 from .cumulantbasedmethod import CumulantBasedLbMethod
 from .galilean_correction import add_galilean_correction
+from .fourth_order_correction import add_fourth_order_correction
 
-__all__ = ['CumulantBasedLbMethod', 'add_galilean_correction']
+__all__ = ['CumulantBasedLbMethod', 'add_galilean_correction', 'add_fourth_order_correction']

lbmpy/methods/cumulantbased/cumulant_simplifications.py

@@ -16,11 +16,8 @@ def insert_logs(ac, **kwargs):
 def insert_log_products(ac, **kwargs):
     def callback(asm):
         rhs = asm.rhs
-        if isinstance(rhs, sp.log):
+        if rhs.find(sp.log):
             return True
-        if isinstance(rhs, sp.Mul):
-            if any(isinstance(arg, sp.log) for arg in rhs.args):
-                return True
         return False
 
     return insert_subexpressions(ac, callback, **kwargs)