Commit b04caa75 authored by Markus Holzer's avatar Markus Holzer Committed by Jan Hönig
Browse files

Dataclasses for method creation

parent 4e47c1eb
......@@ -2,7 +2,7 @@ import numpy as np
import sympy as sp
from pystencils.datahandling import create_data_handling
from pystencils import create_kernel
from pystencils import create_kernel, Target
from pystencils.slicing import make_slice
from lbmpy.creationfunctions import create_lb_collision_rule, create_lb_function
......@@ -18,11 +18,11 @@ from numpy.testing import assert_allclose
all_results = dict()
targets = ['cpu']
targets = [Target.CPU]
try:
import pycuda.autoinit
targets += ['gpu']
targets += [Target.GPU]
except Exception:
pass
......@@ -30,19 +30,19 @@ try:
import pystencils.opencl.autoinit
from pystencils.opencl.opencljit import get_global_cl_queue
if get_global_cl_queue() is not None:
targets += ['opencl']
targets += [Target.OPENCL]
except Exception:
pass
class PeriodicPipeFlow:
def __init__(self, stencil, streaming_pattern, wall_boundary=None, target='cpu'):
def __init__(self, stencil, streaming_pattern, wall_boundary=None, target=Target.CPU):
if wall_boundary is None:
wall_boundary = NoSlip()
self.target = target
self.gpu = target in ['gpu', 'opencl']
self.gpu = target in [Target.GPU, Target.OPENCL]
# Stencil
self.stencil = stencil
......
......@@ -9,7 +9,7 @@ from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
from lbmpy.stencils import get_stencil
from pystencils import create_kernel, create_data_handling, Assignment
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):
......@@ -21,7 +21,7 @@ def flow_around_sphere(stencil, galilean_correction, L_LU, total_steps):
dim = len(stencil[0])
Q = len(stencil)
target = 'gpu'
target = Target.GPU
streaming_pattern = 'aa'
timesteps = get_timesteps(streaming_pattern)
......@@ -49,7 +49,7 @@ def flow_around_sphere(stencil, galilean_correction, L_LU, total_steps):
'pre_simplification': True
}
lb_method = create_lb_method(optimization=optimization, **method_params)
lb_method = create_lb_method(**method_params)
def get_extrapolation_kernel(timestep):
boundary_assignments = []
......
%% Cell type:code id: tags:
``` python
import pytest
pytest.importorskip('pycuda')
```
%%%% Output: execute_result
<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 pystencils import Target
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.target = optimization.get('target', Target.CPU)
self.gpu = self.target in ['gpu', 'opencl']
self.gpu = self.target == Target.GPU
# 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'
target = 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': 'simple',
'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: execute_result
<matplotlib.colorbar.Colorbar at 0x7f4f5b6cbcd0>
<matplotlib.colorbar.Colorbar at 0x7f2eaa33eca0>
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAA0wAAAFoCAYAAABkP3u+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAArS0lEQVR4nO3df4zcd33n8ddrZ2d/eO3E8Y8EE4c6NEa9wLVUWCmn8gct5+D2aA1XUpLmINJFDScaCXTtndKeSHuISiBdS3tXDik0KSEKBALlsK6miSFUHFWbxoG05EfTmjQoNg7+HXvX3t3Zmff9Md/AZrOzfn+/u+ud2Xk+pJFnvvPaz352vjPf9Xs/n+/n64gQAAAAAODlBla6AwAAAADQrSiYAAAAAKADCiYAAAAA6ICCCQAAAAA6oGACAAAAgA4omAAAAACgAwomAAAAAF3J9i7bT9s+YPu2eZ4ftv254vmHbW8rtl9j+7Hi9ve235Ft82Xfg+swAQAAAOg2tmuS/knSTkkHJT0i6YaIeHJW5n2SfjIi/pPt6yW9IyLeZXuNpOmImLG9RdLfS3qlpDhfm3MxwgQAAACgG10j6UBEPBMR05Luk7R7Tma3pLuL+1+Q9BbbjoizETFTbB9Ru1DKtvkSg0vwg6Rt2rQptm3bdiG/JQBcUP/06DPprAeW6W9WXoY2l2kyQrRa6exr3vDq5ekEAHSBRx999FhEbF7pfpTx1p8bi+MnmpW//tF/mHpC0uSsTXdExB2zHl8u6blZjw9K+pk5zfwwU4wmvSBpo6Rjtn9G0l2SfkzSu4vnM22+xAUtmLZt26b9+/dfyG8JABfUzoHr0tmB0TX5hmu1dNRe+oqp1PTtZv6XZ+vs2XR23/77830AgB5j+3sr3Yeyjp1o6uEHtlb++vqW705GxI4l7NJLRMTDkl5r+19Jutv2V6q0c0ELJgAAAACrRagZ+ZkCFRySdMWsx1uLbfNlDtoelHSxpOMv6WXEU7bHJb0u2eZLcA4TAAAAgNJCUktR+ZbwiKTttq+0PSTpekl75mT2SLqpuP9OSQ9FRBRfMyhJtn9M0k9IejbZ5kswwgQAAACg6xTnHN0q6QFJNUl3RcQTtj8kaX9E7JF0p6R7bB+QdELtAkiS3iTpNtsNSS1J74uIY5I0X5sL9YOCCQAAAEAlLS3rlDxFxF5Je+dsu33W/UlJLzuBOCLukXRPts2FUDABAAAAKC0UavbBNV0pmAAAAABUkjwXqaex6AMAAAAAdMAIEwAAAIDSQlKzD0aYKJgArBqlLho7NpbOek3+ArO1q1+TzjbXDqezraH8hWtVW/oL16qZ/4U4MJ2/cG1tfCqd3XXZ+9LZKHFB3NbERDq7r8XFcwFgtn6YkkfBBAAAAKC0kFj0AQAAAAA6Wd5FxbsDiz4AAAAAQAeMMAEAAAAoLRQs+gAAAAAA84pSawL1LAomAAAAAKWF+uMcJgomAAAAABVYTS3DpSy6DIs+AAAAAEAHjDABAAAAKC0ktTiHCQAAAADm1w9T8iiYAHS1a4dvTGdr21+dzjY3rk1nG2vr+ey6/GF1ZjT/S6ZZz2ejlo6muZnP1hr5PzcOnluTztZfuS6fHW+ks7Xj4+lsmffjg1P3prMA0ItC/VEwcQ4TAAAAAHTACBMAAACASlqx+keYKJgAAAAAlNYvU/IomAAAAACUFrKafXCGDwUTAAAAgEr6YUre6i8JAQAAAKAiRpgAAAAAlMY5TAAAAADQkdWM1T9hjYIJAAAAQGkhqdUHZ/hQMAEAAACohCl5AJC0c+C6dLa2cUM669delc6ee8VYOjt5SS2dbazN/zJojOWzzeF0VK0SR+vI/2hpbuazAzP516A2lW+3PpH/K2Z9vJ7OjqzP74iRsZF0dtfm96azzeMnUrl9rfvTbQIAlgYFEwAAAIDSIjiHCQAAAAA6ajElDwAAAABerr2s+OofYVr9PyEAAAAAVMQIEwAAAIAKOIcJAAAAAObFdZgAAAAAYAHNYNEHAAAAAHiZkFn0AQAAAAD6GSNMAAAAACppsegDgH537dANqVztqivTbTa2rE9nz106lM5ObsgftKcuzs+5nlmbjmpmJNLZVols1EtkB/LZLLfyr5cb+ezAZD47WCY7ns821pbIjq1LZ0fXDaez9cMXp3LZz6MkPTj92XQWAKrol+swUTABAAAAKC3kvlj04bwloe0rbH/d9pO2n7D9/mL779k+ZPux4vaLy99dAAAAALhwMiNMM5J+MyK+ZXudpEdt7yue+1hE/I/l6x4AAACAbsV1mCRFxGFJh4v7Z2w/Jeny5e4YAAAAgO4VITX7YNGHUj+h7W2SflrSw8WmW23/g+27bF+y1J0DAAAA0K2s1iJuvSJdMNleK+mLkj4QEaclfULSj0t6vdojUH/Q4etusb3f9v6jR48uvscAAAAAVlyoPcJU9dYrUj21XVe7WLo3Iv5ckiLiBxHRjIiWpE9Kuma+r42IOyJiR0Ts2Lx581L1GwAAAACW3XnPYbJtSXdKeioi/nDW9i3F+U2S9A5Jjy9PFwEAAAB0I67D1Pazkt4t6Tu2Hyu2/Y6kG2y/Xu3RuGclvXcZ+gcAAACgC4WsVh9chymzSt43pXnPytq79N0BAAAA0CsYYQKwKu0cuC6drV11ZSo39aoN6TbPvqKezp7bmD8QT69PRzV9UaSzzbXNfMMj+WxtOJ+tD+azAwP5ny2r1cr/BXFmppbPTpXITuaztbF8tjla4mcbKdHu8HA6u6ae+/zkWyz3Od/Xur9EywDQFpJay7x4g+1dkv5YUk3Sn0bER+Y8Pyzp05LeIOm4pHdFxLO2d0r6iKQhSdOS/ktEPFR8zV9J2iLpXNHMtRFxpFMfKJgAAAAAdB3bNUkfl7RT0kFJj9jeExFPzordLOlkRFxl+3pJH5X0LknHJP1SRHzf9uskPaCXXkv2xojYn+nH6h9DAwAAALAMrOYibgnXSDoQEc9ExLSk+yTtnpPZLenu4v4XJL3FtiPi2xHx/WL7E5JGi9Go0iiYAAAAAJT24pS8qjdJm168Xmtxu2XOt7hc0nOzHh/US0eJXpKJiBlJL0jaOCfzK5K+FRFTs7b9me3HbH+wWBW8I6bkAQAAAKgkOVLUybGI2LFUfZmP7deqPU3v2lmbb4yIQ7bXqX2t2XerfR7UvBhhAgAAANCNDkm6YtbjrcW2eTO2ByVdrPbiD7K9VdKXJL0nIr774hdExKHi3zOSPqP21L+OKJgAAAAAlBbhxU7JO59HJG23faXtIUnXS9ozJ7NH0k3F/XdKeigiwvZ6SX8h6baI+OsXw7YHbW8q7tclvU3S4wt1gil5AAAAACppLuOy4hExY/tWtVe4q0m6KyKesP0hSfsjYo+kOyXdY/uApBNqF1WSdKukqyTdbvv2Ytu1kiYkPVAUSzVJX5X0yYX6QcEEAAAAoLSQ1FrcOUzn/x4ReyXtnbPt9ln3JyW97MJzEfFhSR/u0OwbyvSBggkAAABABV7WEaZuQcEE9KHaxg3pbGPL+lTu7Cvq6TbPXpo/uE5tiHR2en0rnY11M+ns8Nh0OrtmpER2qJHOjtbz2Zrzr0NWmV+I5xol3gvTJbKTQ+ns1Eg+Ozmc/1XYqudfh6iV+U9E7nUYaKzPt3gy/zkHAHRGwQQAAACgtPZ1mJZ3Sl43oGACAAAAUEmzDxbdpmACAAAAUFrIjDABAAAAQCetPhhhWv0/IQAAAABUxAgTAAAAgNIipCZT8gAAAABgfpzDBAAAAADzaC/6sPrP8Fn9PyEAAAAAVMQIEwAAAIBKmmJKHoAece3wjemsX3tVOnvu0qFcbmN+wHpqQ6Sz0xta6azXT6ez69ZOprOXrDmXzq4fzmcvGspnR2uNdLbu/GuW1Sgx5eJcs57Onp4eTWdPjeSzJ4fy2fH6SDo7PZD7PLTlXzM3c9naVP77D555RTpb5vjx4NS96SyA1S3EOUwAAAAA0EF/nMNEwQQAAACgklYfTMlb/SUhAAAAAFTECBMAAACA0rhwLQAAAAAsgHOYAAAAAGAe7QvXMsIEAAAAAPNi0QcAAAAA6GOMMAEAAAAojQvXAgAAAMACWPQBwIraOXBdOlvb/up09twrxtLZyQ25A+H0+nSTml7fSme9fjqdvfiis+nspWvH89nRM+ns5qF8uxcNTqaza2pT6WzdzXQ2qxG1dPZsczidPT00ks4eHVqbzo4Mrktnjwzk348vpJPSdGsonR1o5F7fycn8f0zqE/nP+ejE1nS2zHFpX+v+dBZAD4r+WPRh9ZeEAAAAAFARI0wAAAAASgv1xyp5FEwAAAAAKumHKXkUTAAAAABKY5U8AAAAAFhAPxRMLPoAAAAAAB0wwgQAAACgtFB/LCtOwQQAAACgElbJAwAAAID5RH+cw0TBBAAAAKA0VskDsOIGxsbS2ebGtens5CW1dHbq4tyBcPqiSLcZ62bS2XVrJ9PZS9eOp7OXr3khnd0yks9uqp9JZzfU8v0dG5hOZ+vOv75Zjcj/uphoDaWzJwbz79vRWiOdrbuVzpbRbOXXSjrdyGenp3LZ2rn8f0zq4/nP+VCJ48fA9/PHJQBYDSiYAAAAAFTCCBMAAAAAzKNfVsk77xwA21fY/rrtJ20/Yfv9xfYNtvfZ/ufi30uWv7sAAAAAukWEK996RWbS9Iyk34yIqyW9UdJv2L5a0m2SvhYR2yV9rXgMAAAAAKvGeQumiDgcEd8q7p+R9JSkyyXtlnR3Ebtb0tuXqY8AAAAAulBLrnzrFaXOYbK9TdJPS3pY0mURcbh46nlJl3X4mlsk3SJJr3rVqyp3FAAAAED3iD65DlN6zVPbayV9UdIHIuL07OciItReiv1lIuKOiNgRETs2b968qM4CAAAA6B79cA5TaoTJdl3tYuneiPjzYvMPbG+JiMO2t0g6slydBAAAANBtWCVPkmTbku6U9FRE/OGsp/ZIuqm4f5OkLy999wAAAABg5WRGmH5W0rslfcf2Y8W235H0EUmft32zpO9J+tVl6SEAAACArtRLU+uqOm/BFBHflDouY/GWpe0OgNm8Zk0621hbL5HNH9xm1uZyzbXNdJvDY9Pp7CVrzqWzl46eSWe3jLyQzm4dOpHObh7M92F9bSKdHXP+Nas7vy+yGlFLZydiKJ0dG8j/XCNupLNlNCJ9Oq8mZ/Kfs6np/LpKU5O513dmIt9mmc95mePHSInjEoDVLdQfiz6UWiUPAAAAACRJ0V4pb7XL/1kNAAAAAGZZ7usw2d5l+2nbB2zfNs/zw7Y/Vzz/cHEZJNneaftR298p/v35WV/zhmL7Adv/s1izoSMKJgAAAABdx3ZN0scl/YKkqyXdYPvqObGbJZ2MiKskfUzSR4vtxyT9UkT8a7UXqLtn1td8QtKvS9pe3HYt1A8KJgAAAAClhZb9OkzXSDoQEc9ExLSk+yTtnpPZLenu4v4XJL3FtiPi2xHx/WL7E5JGi9GoLZIuioi/La4l+2lJb1+oE5zDBAAAAKCCRV+HaZPt/bMe3xERd8x6fLmk52Y9PijpZ+a08cNMRMzYfkHSRrVHmF70K5K+FRFTti8v2pnd5uULdZKCCQAAAEAli1z04VhE7FiirszL9mvVnqZ3bdU2mJIHAAAAoBsdknTFrMdbi23zZmwPSrpY0vHi8VZJX5L0noj47qz81vO0+RIUTAAAAAAqWeZzmB6RtN32lbaHJF0vac+czB61F3WQpHdKeigiwvZ6SX8h6baI+Osf9TcOSzpt+43F6njvkfTlhTpBwQQAAACgtIjlLZgiYkbSrZIekPSUpM9HxBO2P2T7l4vYnZI22j4g6T9LenHp8VslXSXpdtuPFbdLi+feJ+lPJR2Q9F1JX1moH5zDBAAAAKCSRS76cF4RsVfS3jnbbp91f1LSdfN83YclfbhDm/slvS7bBwom4ALbOfCyz3RHtatfk8421uU/zo2x/MFtZiR5NudIM93mmpHpdHb98Ll0dvPQeDq7qX4m3+5gmezpdHb9wGQ6u84z6Wx9GX53NUqc1HsmGuls3fn3TRmTUU9nzzXz2dPDo/nsyHA6OzUylMrNjNTSbZb5nJc5foxsviSdLXO829e6P50F0D0WuehDT2BKHgAAAAB0wAgTAAAAgEqSizf0NAomAAAAAKWF0qvd9TQKJgAAAACV9MEpTBRMAAAAACqI/piSx6IPAAAAANABI0wAAAAAqumDOXkUTAAAAAAq6YcpeRRMAAAAACrhwrUAAAAA0McYYQIusIE1a9LZ5trhdHZmND8k3sw3q9ZI7k9HteFmus01Q4109qKhc/ns4GQ6u6E2ns6ur03kswMl+jAwk86OOf/3rXqJbFYjWvnvH/mfS8q/Xo1aLZ2daA2lsy8M5j+TZd6Pa4bG0tnTyc9P9vMoSc3h/DGh1PGjxHGpzPEOQO8JMSUPAAAAAOYXkiiYAAAAAGB+/XAOEwUTAAAAgGr6oGBi0QcAAAAA6IARJgAAAAAVmEUfAAAAAKCjPpiSR8EEAAAAoLzoj2XFOYcJAAAAADpghAkAAABANUzJAwAAAIBOVv+UPAom4EKr1dLR1lA+26znD1itEp/8qOf+dFQfbKbbHK038tlaPrumNpXOjg1M57POZ9d5pkS7+VnRawbq6eyg8u+brLrz+1et/D5rlHi9zpbYD2X2b5n3TZn3Y5n3+WDy8zOd/DxK5T7npY4fJY5LAyWOdwB6FCNMAAAAANBBHxRMLPoAAAAAAB0wwgQAAACgvJDUB8uKUzABAAAAqCT6YEoeBRMAAACAaiiYAAAAAKCDPpiSx6IPAAAAANABI0wAAAAAKjFT8gAAAABgHiHOYQIAAACA+bkvzmGiYAIuMLvEgaWWz0Yt32yp7EDuT0cDyZwk1dxKZ+ulss0S2ZllajcdVd3500gHld9ptRLtppX4C2K516vX9m++v2Xe59nPT/bzKC3fMaHMcanU8Q4AuhQFEwAAAIBq+mBK3nn/DGn7LttHbD8+a9vv2T5k+7Hi9ovL200AAAAAXScWcesRmXkbn5K0a57tH4uI1xe3vUvbLQAAAABdrw8KpvNOyYuIb9jedgH6AgAAAKBXhPpi0YfFnBl8q+1/KKbsXdIpZPsW2/tt7z969Ogivh0AAAAAXFhVC6ZPSPpxSa+XdFjSH3QKRsQdEbEjInZs3ry54rcDAAAA0G0c1W+9olLBFBE/iIhmRLQkfVLSNUvbLQAAAABdrw/OYapUMNneMuvhOyQ93ikLAAAAAL3qvIs+2P6spDdL2mT7oKTflfRm269XuzZ8VtJ7l6+LAAAAALpRL02tqyqzSt4N82y+cxn6AgAAAABd5bwFEwAAAADMqw+WFadgAgAAAFBejy3eUBUFEwAAAIBq+qBgWsyFawEAAABgVaNgAgAAAFDJcl+41vYu20/bPmD7tnmeH7b9ueL5h21vK7ZvtP112+O2/2TO1/xV0eZjxe3ShfrAlDwAAAAA1SzjlDzbNUkfl7RT0kFJj9jeExFPzordLOlkRFxl+3pJH5X0LkmTkj4o6XXFba4bI2J/ph+MMAEAAACoJhZxO79rJB2IiGciYlrSfZJ2z8nslnR3cf8Lkt5i2xExERHfVLtwWhQKJgAAAAClLWY6XjElb5Pt/bNut8z5FpdLem7W44PFtnkzETEj6QVJGxPd/7NiOt4HbS+4NjpT8gAAAACshGMRsWMFvu+NEXHI9jpJX5T0bkmf7hRmhAkAAABANeHqt/M7JOmKWY+3FtvmzdgelHSxpOMLdjniUPHvGUmfUXvqX0eMMAEXWESJsyOb+ayb+WZLZVu5K3i3kjlJakb+bzWNUtlaiWz+8Feu3XRUjWils/UyO20ZTsCdUf77l/m5yr1e3bB/8+/HMu/z7Ocn+3mUlu+YUOa4VOp4B6A3Le/H/BFJ221fqXZhdL2kX5uT2SPpJkl/I+mdkh6KBQ4+RVG1PiKO2a5Lepukry7UCQomAAAAAJVklwevIiJmbN8q6QFJNUl3RcQTtj8kaX9E7JF0p6R7bB+QdELtoqrdN/tZSRdJGrL9dknXSvqepAeKYqmmdrH0yYX6QcEEAAAAoJplHkiOiL2S9s7Zdvus+5OSruvwtds6NPuGMn3gHCYAAAAA6IARJgAAAADlxfJOyesWFEwAAAAAqqFgAgAAAIAO+qBg4hwmAAAAAOiAESYAAAAAlfTDOUyMMAEAAABAB4wwAQAAAKimD0aYKJiAC63ZTEcHpvPZWiN/xBqYcTrrRi47M1NLt3muUc9nm/ns2eZwOjvRGspnI589E410th4z6axaJdp1/n2T1YhWOjtRInsm8r+GyuyHMvu3zPumzPuxzPs8+/nJfh4laaDE26vU8aPEcanM8Q5AD+qTZcWZkgcAAAAAHTDCBAAAAKCaPhhhomACAAAAUA0FEwAAAAC8nNUf5zBRMAEAAACopg8KJhZ9AAAAAIAOGGECAAAAUF6fLCtOwQQAAACgGgomAAAAAOigDwomzmECAAAAgA4YYQIusNbZs+lsbXwqnR08tybfbr5ZDUw6lZuZqqXbPDtdT2dPT4/ms0Mj6eyJwbXp7NjAdDpbdzOdlSbTyYZnSvShVaIPye9f4i+IZyL/q+VUK7/PTjXH0tkTzfz+PT2T70OZ92OZ93kz+fkZTH4epXKf88Fz+R1c5rjULHG8A9CbOIcJAAAAADqhYAIAAACAeYQomAAAAACgk36YkseiDwAAAADQASNMAAAAAKrpgxEmCiYAAAAAlfTDlDwKJgAAAADVUDABAAAAwDz6ZJU8Fn0AAAAAgA4YYQIAAABQmovbakfBBFxg+1r3p7O7LntfOlt/5bp8diI/uDw4mTsUzkzW0m2enRxKZ0+NjKazR4fWprOjtUY6O+J8toxGrcRr5ul0tu5mle4sqBH5vk5Eif3bHEtnj87k3+PHGvns0en8++bUVP79WOZ9ruTnJ/t5lKT6RH6eTP3MTDqroyfT0TLHOwA9qg+m5FEwAQAAAKikH1bJ4xwmAAAAAOjgvAWT7btsH7H9+KxtG2zvs/3Pxb+XLG83AQAAAHSdWMStR2RGmD4ladecbbdJ+lpEbJf0teIxAAAAgH5CwSRFxDcknZizebeku4v7d0t6+9J2CwAAAEBXi/Y5TFVvvaLqOUyXRcTh4v7zki7rFLR9i+39tvcfPXq04rcDAAAA0HUYYTq/iFjwR46IOyJiR0Ts2Lx582K/HQAAAABcMFULph/Y3iJJxb9Hlq5LAAAAAHoBU/I62yPppuL+TZK+vDTdAQAAANAz+mBK3nkvXGv7s5LeLGmT7YOSflfSRyR93vbNkr4n6VeXs5MAAAAAuk8vjRRVdd6CKSJu6PDUW5a4LwDmiLNn09n6eKNEtp7ODo47lauN1dJtTo0MpbMnh0bT2ZHBdels3a10tozJyL+2E6386zA2MJ3O1j2TzmY14ry/Ln6ozM91ork2nT3WyO/fw5MXp7NHzuXbPXk2/36cmsi/DrXx3OdncDzdpOrj+f/FlDl+lDkuAcBqkP8NCAAAAAAv6rGpdVVRMAEAAACohoIJAAAAAF7O4hwmAAAAAOisDwqmRV+4FgAAAABWKwomAAAAAJU4ovIt1b69y/bTtg/Yvm2e54dtf654/mHb24rtG21/3fa47T+Z8zVvsP2d4mv+p+0FlwSmYAIAAABQ3mIuWpuol2zXJH1c0i9IulrSDbavnhO7WdLJiLhK0sckfbTYPinpg5J+a56mPyHp1yVtL267FuoHBRMAAACAShzVbwnXSDoQEc9ExLSk+yTtnpPZLenu4v4XJL3FtiNiIiK+qXbh9KP+2lskXRQRfxsRIenTkt6+UCcomAAAAABUs7gRpk2298+63TKn9cslPTfr8cFi27yZiJiR9IKkjQv0+PKinYXafAlWyQMAAACwEo5FxI6V7sT5UDABXaw1MZHO1o6Pp7Mj64fT2cbaBc+D/KHmaC4nSZPD+UPPeH0knT0y0Epny2hEfjD+XLOezr4wuCadXVObSmfrbqazWY2opbNnm/n31+mZ/P49Or02nT1ybl0+O55vd3w831+fyb/Ph07nPj/DL+TX7x05mX8flDl+NEsclwCsfst8HaZDkq6Y9XhrsW2+zEHbg5IulnT8PG1uPU+bL8GUPAAAAADVLOOiD5IekbTd9pW2hyRdL2nPnMweSTcV998p6aHi3KT5uxtxWNJp228sVsd7j6QvL9QJRpgAAAAAlJdfvKFa8xEztm+V9ICkmqS7IuIJ2x+StD8i9ki6U9I9tg9IOqF2USVJsv2spIskDdl+u6RrI+JJSe+T9ClJo5K+Utw6omACAAAA0JUiYq+kvXO23T7r/qSk6zp87bYO2/dLel22DxRMAAAAAKpZ3nOYugIFEwAAAIDSrGVf9KErUDABAAAAqKbz+gqrBgUTAAAAgEr6YYSJZcUBAAAAoANGmAAAAACUl7+eUk+jYAIAAABQiVsr3YPlR8EEdLF9rfvT2WuHb0xnR8ZG0tnG2LpUbmaklm6zVc/PBp4eGEpnX0gnpWYr34fJmXo6e3p4NJ29aOhcOjtaa6Sz9WX47dWI/Ot1rlni9ZrOv16npvLZk2fz2fHx/OchTuXfj0On8q/Z0KlcbuREft+OPD+RzraePZjOljkuAegDjDABAAAAwPxY9AEAAAAA+hgjTAAAAADKC3EdJgAAAADopB+m5FEwAQAAAKimDwomzmECAAAAgA4YYQIAAABQmsWUPAAAAACYXwSLPgAAAABAJ4wwAQAAAEAnFEwAesWDU/ems7s2vzedHV03nMo1h3M5SYpamfVm8tnp1lA6e7qRb3dqOn+oPD2Sfx3WDI2ls6P1RjpbcyudzWpG/vU616ins2enS2Qn8/t3aiKf9Zn8/h06lX8dhk84nR09nttno0em0236uefT2QdKHD8AoN9QMAEAAACohCl5AAAAADCfkNRa/RUTBRMAAACAalZ/vUTBBAAAAKCafpiSV+bMawAAAADoK4wwAQAAAKiGC9cCAAAAwPz6YUoeBRMAAACA8kJ9segD5zABAAAAQAeMMAEAAAAozZLMOUwAVqPm8RPpbP3wxancmvqGEj2op5Nu5gfCBxq1dHZ6Kt/u1GS+3amRoXT29HAznR0czGcHBpb+l1er5XR2Zib/ejWn8lmV2A+18Xx26HT+Zxs6lY5q9HgrnV3zfCOVqx/Od6DM5xwAKssf6noWBRMAAACAShhhOg/bz0o6I6kpaSYidixFpwAAAAB0uT5Z9GEpRph+LiKOLUE7AAAAANBVmJIHAAAAoILoiwvXLnZZ8ZD0oO1Hbd8yX8D2Lbb3295/9OjRRX47AAAAAN3CUf3WKxY7wvSmiDhk+1JJ+2z/Y0R8Y3YgIu6QdIck7dixo4deGgAAAAALYoRpYRFxqPj3iKQvSbpmKToFAAAAoMuF5Fb1W6+oXDDZHrO97sX7kq6V9PhSdQwAAAAAVtpipuRdJulLtl9s5zMR8ZdL0isAAAAA3a8PpuRVLpgi4hlJP7WEfQEAAADQS1Z/vcSy4kA/2te6P529duiGVG64xPcfaKxPZ2tTQ+ns5GR+lnHtnNPZmYn8oXJmpJbOtkbyv2Wm6/lsDCz9by+38q+XG/ns4OQyZcfTUQ2/kH+9Rk7kJ92PHplOZ+uHT6VyrWefS7dZ5nMOAFW5D0aYFrusOAAAAACsWowwAQAAAKimD0aYKJgAAAAAlBeSemh58KoomAAAAACUZkVfnMNEwQQAAACgmj4omFj0AQAAAAA6YIQJAAAAQDWMMAEAAADAPF5c9KHqLcH2LttP2z5g+7Z5nh+2/bni+Ydtb5v13G8X25+2/dZZ25+1/R3bj9nef74+MMIEAAAAoJLlXPTBdk3SxyXtlHRQ0iO290TEk7NiN0s6GRFX2b5e0kclvcv21ZKul/RaSa+U9FXbr4mIZvF1PxcRxzL9oGACsKAHpz+byu0cuC7dZv3khnR28Mwr8u1OjOWz47V0trHW+exYPtsczmdbJY7Wkf/R0tw8f+ZFAzP5bG0qn61P5H8p18fz2ZGT+R9u5PmJdNbPPZ/ONo+fSOX2te5PtwkAF8TyTsm7RtKBiHhGkmzfJ2m3pNkF025Jv1fc/4KkP7HtYvt9ETEl6V9sHyja+5uynWBKHgAAAICVsMn2/lm3W+Y8f7mk52Y9PlhsmzcTETOSXpC08TxfG5IetP3oPN/zZRhhAgAAAFBBLHaE6VhE7Fiq3pTwpog4ZPtSSfts/2NEfKNTmBEmAAAAAOWF2gVT1dv5HZJ0xazHW4tt82ZsD0q6WNLxhb42Il7894ikL6k9Va8jCiYAAAAA1SzvKnmPSNpu+0rbQ2ov4rBnTmaPpJuK+++U9FBERLH9+mIVvSslbZf0d7bHbK+TJNtjkq6V9PhCnWBKHgAAAICuExEztm+V9ICkmqS7IuIJ2x+StD8i9ki6U9I9xaIOJ9QuqlTkPq/2AhEzkn4jIpq2L5P0pfa6EBqU9JmI+MuF+kHBBAAAAKCS5VxWXJIiYq+kvXO23T7r/qSkeZfqjYjfl/T7c7Y9I+mnyvSBggkAAABANctcMHUDCiYAAAAA5YWkFgUTAAAAAMxj0cuK9wRWyQMAAACADhhhArAk9rXuX5Z2rx2+MZ0dndiazg5tXJvONtbW89l1+cPqzKjT2WY9n41aOprmZj5ba+T/2jh4Lp+tn5nJZ8cb6Wzt+Hg623r2YDr7wNS96SwA9Kw+GGGiYAIAAABQDQUTAAAAAMyDRR8AAAAAoJOQorXSnVh2LPoAAAAAAB0wwgQAAACgGs5hAgAAAIB5cA4TAAAAACygD0aYOIcJAAAAADpghAkAAABANX0wwkTBBAAAAKCCoGACgJX24NS96ezOgevS2YHvj6WzI2vW5LObL0lnm2uH09nWUC2dVc35bFYz/wtxYLqZztbGp/J9OHoyHY2zZ9PZ5sREOruvdX86CwCrXkhqrf7rMFEwAQAAAKimD0aYWPQBAAAAADpghAkAAABANX0wwkTBBAAAAKCC4MK1AAAAADCvkCJW/6IPnMMEAAAAAB0wwgQAAACgGqbkAQAAAEAHLPoAAAAAAPOI4MK1AAAAANARI0wA0Dv2te5f6S5o58B16ezAmjX5bK2WztpOZ7OizC/EZjMfPXs2ne2G/QsA6D8UTAAAAAAqiT6YkreoZcVt77L9tO0Dtm9bqk4BAAAA6HbRnpJX9dYjKo8w2a5J+riknZIOSnrE9p6IeHKpOgcAAACgS4VYVvw8rpF0ICKekSTb90naLYmCCQAAAOgHwZS8hVwu6blZjw8W217C9i2299vef/To0UV8OwAAAAC4sBZ1DlNGRNwRETsiYsfmzZuX+9sBAAAAuABCUrSi8q1XLGZK3iFJV8x6vLXYBgAAAGC1i+iLKXmLKZgekbTd9pVqF0rXS/q1JekVAAAAgK7XSyNFVVUumCJixvatkh6QVJN0V0Q8sWQ9AwAAAIAVtqgL10bEXkl7l6gvAAAAAHpJH0zJc1zAi0bZPirpexfsG/a3TZKOrXQnkMb+6j3ss97DPus97LPewv5anB+LiJ5aIc32X6q936s6FhG7lqo/y+WCFky4cGzvj4gdK90P5LC/eg/7rPewz3oP+6y3sL+wWi37suIAAAAA0KsomAAAAACgAwqm1euOle4ASmF/9R72We9hn/Ue9llvYX9hVeIcJgAAAADogBEmAAAAAOiAggkAAAAAOqBgWkVsX2f7Cdst2zvmPPfbtg/Yftr2W1eqj3g527uK/XLA9m0r3R+8nO27bB+x/fisbRts77P9z8W/l6xkH/Ejtq+w/XXbTxbHxPcX29lnXcr2iO2/s/33xT7778X2K20/XBwfP2d7aKX7ih+xXbP9bdv/t3jM/sKqRMG0ujwu6d9L+sbsjbavlnS9pNdK2iXpf9uuXfjuYa5iP3xc0i9IulrSDcX+Qnf5lNqfndluk/S1iNgu6WvFY3SHGUm/GRFXS3qjpN8oPlfss+41JennI+KnJL1e0i7bb5T0UUkfi4irJJ2UdPPKdRHzeL+kp2Y9Zn9hVaJgWkUi4qmIeHqep3ZLui8ipiLiXyQdkHTNhe0dOrhG0oGIeCYipiXdp/b+QheJiG9IOjFn825Jdxf375b09gvZJ3QWEYcj4lvF/TNq/4fucrHPula0jRcP68UtJP28pC8U29lnXcT2Vkn/TtKfFo8t9hdWKQqm/nC5pOdmPT5YbMPKY9/0rssi4nBx/3lJl61kZzA/29sk/bSkh8U+62rF9K7HJB2RtE/SdyWdioiZIsLxsbv8kaT/KqlVPN4o9hdWKQqmHmP7q7Yfn+fGqASwQqJ9fQau0dBlbK+V9EVJH4iI07OfY591n4hoRsTrJW1Ve/T9J1a2R+jE9tskHYmIR1e6L8CFMLjSHUA5EfFvK3zZIUlXzHq8tdiGlce+6V0/sL0lIg7b3qL2X8XRJWzX1S6W7o2IPy82s896QEScsv11Sf9G0nrbg8WoBcfH7vGzkn7Z9i9KGpF0kaQ/FvsLqxQjTP1hj6TrbQ/bvlLSdkl/t8J9QtsjkrYXKwsNqb04x54V7hNy9ki6qbh/k6Qvr2BfMEtxLsWdkp6KiD+c9RT7rEvZ3mx7fXF/VNJOtc89+7qkdxYx9lmXiIjfjoitEbFN7d9bD0XEjWJ/YZVye1YCVgPb75D0vyRtlnRK0mMR8dbiuf8m6T+qvXrUByLiKyvVT7xU8Re6P5JUk3RXRPz+yvYIc9n+rKQ3S9ok6QeSflfS/5H0eUmvkvQ9Sb8aEXMXhsAKsP0mSf9P0nf0o/Mrfkft85jYZ13I9k+qvUhATe0/5n4+Ij5k+9VqL4azQdK3Jf2HiJhauZ5iLttvlvRbEfE29hdWKwomAAAAAOiAKXkAAAAA0AEFEwAAAAB0QMEEAAAAAB1QMAEAAABABxRMAAAAANABBRMAAAAAdEDBBAAAAAAd/H+KhTErP8zvRQAAAABJRU5ErkJggg==)
%% 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: execute_result
<matplotlib.colorbar.Colorbar at 0x7f4f5b34bf40>
<matplotlib.colorbar.Colorbar at 0x7f2ea7a6e910>
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAA0wAAAFoCAYAAABkP3u+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAArdklEQVR4nO3de5BcZ33m8eeZnp6LZiTLlmVjLBOZWGzWsAkpVA5b4Q8SIqNkSQQbHOx4wVXritkiroLaZLecbOFkKVIFVZuQG0uViB2My2AwhEW1EbEFJsWSShzL4ARf4iAcs5bwRRfb0ow0Mz3dv/2jj5fxeHr0O2du3dPfT1WXuk8/ffqdPt2n9ev3Pe9xRAgAAAAA8HIDa90AAAAAAOhWFEwAAAAA0AEFEwAAAAB0QMEEAAAAAB1QMAEAAABABxRMAAAAANABBRMAAACArmR7t+3HbB+yfdMC9w/b/lxx/322txfLr7D9YHH5B9vvyK7zZc/BeZgAAAAAdBvbNUn/LGmXpMOS7pd0TUQ8MifzPkk/HhH/yfbVkt4REe+yvUHSTETM2r5I0j9IeqWkONs656OHCQAAAEA3ukLSoYh4PCJmJN0pac+8zB5JtxXXvyDpLbYdEacjYrZYPqJ2oZRd50sMLsMfknb++efH9u3bV/MpAWBV/fMDj6ezHlih36y8AutcocEI0Wqls695w6tXphEA0AUeeOCBYxGxda3bUcZbf2Ysjp9oVn78A/84fXdE7F4kcrGkJ+fcPizppzplit6kFyRtkXTM9k9JulXSj0h6d3F/Zp0vsaoF0/bt23Xw4MHVfEoAWFW7Bq5KZwdGN+RXXKulo/byV0ylhm8381+erdOn09kDB+/KtwEAeozt7691G8o6dqKp++7eVvnx9Yu+92O25xYHeyNi79Jb1hYR90l6re1/Lek221+psp5VLZgAAAAArBehZuRHCizgWETsXOT+I5IumXN7W7Fsocxh24OSzpF0/CWtjHjU9oSk1yXX+RIcwwQAAACgtJDUUlS+JNwvaYftS20PSbpa0r55mX2Sriuuv1PSvRERxWMGJcn2j0j6MUlPJNf5EvQwAQAAAOg6xTFHN0q6W1JN0q0R8bDtD0k6GBH7JN0i6XbbhySdULsAkqQ3SbrJdkNSS9L7IuKYJC20zsXaQcEEAAAAoJKWljQk76wiYr+k/fOW3Tzn+pSklx1AHBG3S7o9u87FUDABAAAAKC0UavbBOV0pmAAAAABUkjwWqacx6QMAAAAAdEAPEwAAAIDSQlKzD3qYKJgArBulTho7NpbOekP+BLO1y1+TzjbHh9PZ1lD+xLWqLf+Ja9XMfyEOzORPXFubmE5nd1/4vnQ2SpwQtzU5mc4eaHHyXACYqx+G5FEwAQAAACgtJCZ9AAAAAIBOVnZS8e7ApA8AAAAA0AE9TAAAAABKCwWTPgAAAADAgqLUnEA9i4IJAAAAQGmh/jiGiYIJAAAAQAVWUytwKosuw6QPAAAAANABPUwAAAAASgtJLY5hAgAAAICF9cOQPAomAF3tyuFr09najlens80t4+lsY7yez27M71ZnR/NfMs16Phu1dDTNzXy21sj/3Dh4ZkM6W3/lxnx2opHO1o5PpLNl3o/3TN+RzgJALwr1R8HEMUwAAAAA0AE9TAAAAAAqacX672GiYAIAAABQWr8MyaNgAgAAAFBayGr2wRE+FEwAAAAAKumHIXnrvyQEAAAAgIroYQIAAABQGscwAQAAAEBHVjPW/4A1CiYAAAAApYWkVh8c4UPBBAAAAKAShuQBQNKugavS2dqW89JZv/aydPbMK8bS2alza+lsYzz/ZdAYy2ebw+moWvV8diVGR7iVzw408q9BbTq/3vpk/g+rT+RfsJHN+Q0xMjaSzu7e+t50tnn8RCp3oHVXep0AgOVBwQQAAACgtAiOYQIAAACAjloMyQMAAACAl2tPK77+e5jW/18IAAAAABXRwwQAAACgAo5hAgAAAIAFcR4mAAAAAFhEM5j0AQAAAABeJmQmfQAAAACAfkYPEwAAAIBKWkz6AKDfXTl0TSpX2/Hq9DobF25KZ89cOJzOTp2X32lPn5Mfcz07no5qdjTS2dZwPhv1EtmBfDbLrfzr5UY+OzCdzw6eKZGdyGcb4yWyYxvT2dHxoXS2/szmVC77eZSke2Y+m84CQBX9ch4mCiYAAAAApYXcF5M+nLUktH2J7a/bfsT2w7bfXyz/XdtHbD9YXH5h5ZsLAAAAAKsn08M0K+k3IuJbtjdKesD2geK+j0XE/1i55gEAAADoVpyHSVJEPCXpqeL6KduPSrp4pRsGAAAAoHtFSM0+mPSh1F9oe7ukn5R0X7HoRtv/aPtW2+cud+MAAAAAdCurtYRLr0gXTLbHJX1R0gci4qSkT0j6UUmvV7sH6vc7PO4G2wdtHzx69OjSWwwAAABgzYXaPUxVL70i1VLbdbWLpTsi4i8kKSKeiYhmRLQkfVLSFQs9NiL2RsTOiNi5devW5Wo3AAAAAKy4sx7DZNuSbpH0aET8wZzlFxXHN0nSOyQ9tDJNBAAAANCNOA9T209Lerek79h+sFj225Kusf16tXvjnpD03hVoHwAAAIAuFLJafXAepswsed+UFjwqa//yNwcAAABAr6CHCcC6tGvgqnS2tuPVqdz0JfmJMk+/op7OntmS3xHPbE5HNbMp0tnmeDO/4pF8tjacz9YH89mBgfzfltVq5X9BnJ2t5bPTJbJT+WxtLJ9tjpb420ZKrHd4JJ3dMJRb73B6jeU+5wdad5VYMwC0haTWCk/eYHu3pD+SVJP0ZxHxkXn3D0v6tKQ3SDou6V0R8YTtXZI+ImlI0oyk/xIR9xaP+WtJF0k6U6zmyoh4tlMbKJgAAAAAdB3bNUkfl7RL0mFJ99veFxGPzIldL+m5iLjM9tWSPirpXZKOSfrFiPiB7ddJulsvPZfstRFxMNOO9d+HBgAAAGAFWM0lXBKukHQoIh6PiBlJd0raMy+zR9JtxfUvSHqLbUfEtyPiB8XyhyWNFr1RpVEwAQAAACjtxSF5VS+Szn/xfK3F5YZ5T3GxpCfn3D6sl/YSvSQTEbOSXpC0ZV7mlyV9KyKm5yz7c9sP2v5gMSt4RwzJAwAAAFBJsqeok2MRsXO52rIQ269Ve5jelXMWXxsRR2xvVPtcs+9W+zioBdHDBAAAAKAbHZF0yZzb24plC2ZsD0o6R+3JH2R7m6QvSXpPRHzvxQdExJHi31OSPqP20L+OKJgAAAAAlBbhpQ7JO5v7Je2wfantIUlXS9o3L7NP0nXF9XdKujciwvZmSX8p6aaI+JsXw7YHbZ9fXK9LepukhxZrBEPyAAAAAFTSXMFpxSNi1vaNas9wV5N0a0Q8bPtDkg5GxD5Jt0i63fYhSSfULqok6UZJl0m62fbNxbIrJU1KursolmqSvirpk4u1g4IJAAAAQGkhqbW0Y5jO/hwR+yXtn7fs5jnXpyS97MRzEfFhSR/usNo3lGkDBRMAAACACryiPUzdgoIJ6EO1Leels40LN6Vyp19RT6/z9AX5nev0eZHOzmxupbOxcTadHR6bSWc3jJTIDjXS2dF6Pltz/nXIKvOFeKZR4r0wUyI7NZTOTo/ks1PD+a/CVj3/OkStzH8icq/DwEzu8yhJ9RP5zzkAoDMKJgAAAACltc/DtLJD8roBBRMAAACASpp9MOk2BRMAAACA0kKmhwkAAAAAOmn1QQ/T+v8LAQAAAKAiepgAAAAAlBYhNRmSBwAAAAAL4xgmAAAAAFhAe9KH9X+Ez/r/CwEAAACgInqYAAAAAFTSFEPyAPSIK4evTWf92svS2TMXDudyW/Id1tPnRTo7c14rnfXmmXR24/hUOnvuhjPp7ObhfHbTUD47Wmuks3XnX7OsRokhF2ea9XT25MxoOvv8SD773FA+O1EfSWdnBobS2TKDONzMZWvTuc+jJA1OvCKdLbP/uGf6jnQWwPoW4hgmAAAAAOigP45homACAAAAUEmrD4bkrf+SEAAAAAAqoocJAAAAQGmcuBYAAAAAFsExTAAAAACwgPaJa+lhAgAAAIAFMekDAAAAAPQxepgAAAAAlMaJawEAAABgEUz6AGBN7Rq4Kp2t7Xh1OnvmFWPp7NR5uR3hzOb0KjWzuZXOevNMOnvOptPp7AXjE/ns6Kl0dutQfr2bBqfS2Q216XR2xLPpbNZU5L8uTjeH09mTQyPp7NGh8XR2ZHBjOvvsQP79+EI6Kc20htLZgUYtlZuayv/HpD6Z/5yPTm5LZ8vslw607kpnAfSg6I9JH9Z/SQgAAAAAFdHDBAAAAKC0UH/MkkfBBAAAAKCSfhiSR8EEAAAAoDRmyQMAAACARfRDwcSkDwAAAADQAT1MAAAAAEoL9ce04hRMAAAAACphljwAAAAAWEj0xzFMFEwAAAAASmOWPABrbmBsLJ1tbhlPZ6fOraWz0+fkdoQzmyK9ztg4m85uHJ9KZy8Yn0hnL97wQjp70Ug+e379VDp7Xi3f3rGBmXS27vzrm9WI/NfFZGsonT0xmH/fjtYa6WzdrXS2jGYrP1fSyUY+OzOdy9bO5P9jUp/If86HSuw/Bn6Q3y8BwHpAwQQAAACgkn7oYTrrT1q2L7H9dduP2H7Y9vuL5efZPmD7u8W/5658cwEAAAB0gxdnyat66RWZMQCzkn4jIi6X9EZJv277ckk3SfpaROyQ9LXiNgAAAIA+EeHKl15x1oIpIp6KiG8V109JelTSxZL2SLqtiN0m6e0r1EYAAAAAWBOljmGyvV3ST0q6T9KFEfFUcdfTki5c3qYBAAAA6Gach2kO2+OSvijpAxFx0v7hixMRYXvBKbJs3yDpBkl61atetbTWAgAAAOgK0SfnYUrNY2q7rnaxdEdE/EWx+BnbFxX3XyTp2YUeGxF7I2JnROzcunXrcrQZAAAAQBfgGCZJbncl3SLp0Yj4gzl37ZN0XXH9OklfXv7mAQAAAOhO/TFLXmZI3k9Lerek79h+sFj225I+Iunztq+X9H1Jv7IiLQQAAACANXLWgikivil1PJrrLcvbHAAAAAC9opeG1lVVapY8AKvLGzaks43xeolsfuc2O57LNceb6XUOj82ks+duOJPOXjB6Kp29aOSFdHbb0Il0dutgvg2ba5Pp7Jjzr1nd+W2R1YhaOjsZQ+ns2ED+7xpxI50toxGpw3klSVOz+c/Z9Ez+K3Z6Kvf6zk7m11nmc15m/zFSYr8EYH0L9cekDxRMAAAAAMqL9kx5613+ZzUAAAAAmKMlV75k2N5t+zHbh2zftMD9w7Y/V9x/X3HeWNneZfsB298p/v3ZOY95Q7H8kO0/9tzzJS2AggkAAABA17Fdk/RxST8v6XJJ19i+fF7seknPRcRlkj4m6aPF8mOSfjEi/o3aM3rfPucxn5D0a5J2FJfdi7WDggkAAABAaaEVPw/TFZIORcTjETEj6U5Je+Zl9ki6rbj+BUlvse2I+HZE/KBY/rCk0aI36iJJmyLi7yIiJH1a0tsXawTHMAEAAACoYMXPp3SxpCfn3D4s6ac6ZSJi1vYLkrao3cP0ol+W9K2ImLZ9cbGeueu8eLFGUDABAAAAqGSJkz6cb/vgnNt7I2Lv0lr0UrZfq/YwvSurroOCCQAAAMBaOBYROxe5/4ikS+bc3lYsWyhz2PagpHMkHZck29skfUnSeyLie3Py286yzpfgGCYAAAAAlazwMUz3S9ph+1LbQ5KulrRvXmaf2pM6SNI7Jd0bEWF7s6S/lHRTRPzND9sbT0k6afuNxex475H05cUaQcEEAAAAoLSIlS2YImJW0o2S7pb0qKTPR8TDtj9k+5eK2C2Sttg+JOk/S3px6vEbJV0m6WbbDxaXC4r73ifpzyQdkvQ9SV9ZrB0MyQMAAABQyQpP+qCI2C9p/7xlN8+5PiXpqgUe92FJH+6wzoOSXpdtAwUTsMp2DbzsM91R7fLXpLONjfmPc2Msv3ObHU0ezTnSTK9zw8hMOrt5+Ew6u3VoIp09v34qv97BMtmT6ezmgal0dqNn09n6Cnx3NUoc1HsqGuls3fn3TRlTUU9nzzTz2ZPDo/nsyHA6Oz0ylMrNjtbS6yzzOS+z/xjZem46W2Z/d6B1VzoLoHsscdKHnsCQPAAAAADogB4mAAAAAJUkJ2/oaRRMAAAAAEoLpWe762kUTAAAAAAq6YNDmCiYAAAAAFQQ/TEkj0kfAAAAAKADepgAAAAAVNMHY/IomAAAAABU0g9D8iiYAAAAAFTCiWsBAAAAoI/RwwSssoENG9LZ5vhwOjs7mu8Sb+ZXq9Zw7qej2nAzvc4NQ410dtPQmXx2cCqdPa82kc5urk3mswMl2jAwm86OOf/7Vr1ENqsRrfzzR/7vkvKvV6NWS2cnW0Pp7AuD+c9kmffjhqGxdPZk8vOT/TxKUnM4v08otf8osV8qs78D0HtCDMkDAAAAgIWFJAomAAAAAFhYPxzDRMEEAAAAoJo+KJiY9AEAAAAAOqCHCQAAAEAFZtIHAAAAAOioD4bkUTABAAAAKC/6Y1pxjmECAAAAgA7oYQIAAABQDUPyAAAAAKCT9T8kj4IJWG21WjraGspnm/X8DqtVT0cV9dxPR/XBZnqdo/VGPlvLZzfUptPZsYGZfNb57EbPllhvflT0hoH8RhtU/n2TVXd++6qV32aNEq/X6RLbocz2LfO+KfN+LPM+H0x+fmaSn0ep3Oe81P6jxH5poMT+DkCPoocJAAAAADrog4KJSR8AAAAAoAN6mAAAAACUF5L6YFpxCiYAAAAAlUQfDMmjYAIAAABQDQUTAAAAAHTQB0PymPQBAAAAADqghwkAAABAJWZIHgAAAAAsIMQxTAAAAACwMPfFMUwUTMAqs0vsWGr5bNTyq40SRy/GQO6no4FkTpJqbqWz9RLZEc+WWG+ZbLNENh1V3fkNMaj8Bq6VWG9aiV8Qy71eZd4LZdab377l3jf59pZ5n2c/P9nPYzubjpbaf5TZL5Xa3wFAl6JgAgAAAFBNHwzJO+vvT7Zvtf2s7YfmLPtd20dsP1hcfmFlmwkAAACg68QSLj0i02H/KUm7F1j+sYh4fXHZv7zNAgAAAND1+qBgOuuQvIj4hu3tq9AWAAAAAL0i1BeTPizlyOAbbf9jMWTv3E4h2zfYPmj74NGjR5fwdAAAAACwuqoWTJ+Q9KOSXi/pKUm/3ykYEXsjYmdE7Ny6dWvFpwMAAADQbRzVL72iUsEUEc9ERDMiWpI+KemK5W0WAAAAgK7XB8cwVSqYbF805+Y7JD3UKQsAAAAAveqskz7Y/qykN0s63/ZhSb8j6c22X692bfiEpPeuXBMBAAAAdKNeGlpXVWaWvGsWWHzLCrQFAAAAALrKWQsmAAAAAFhQH0wrTsEEAAAAoLwem7yhKgomAAAAANX0QcG0lBPXAgAAAMC6RsEEAAAAoJKVPnGt7d22H7N9yPZNC9w/bPtzxf332d5eLN9i++u2J2z/6bzH/HWxzgeLywWLtYEheQAAAACqWcEhebZrkj4uaZekw5Lut70vIh6ZE7te0nMRcZntqyV9VNK7JE1J+qCk1xWX+a6NiIOZdtDDBAAAAKCaWMLl7K6QdCgiHo+IGUl3StozL7NH0m3F9S9IeottR8RkRHxT7cJpSSiYAAAAAJS2lOF4xZC8820fnHO5Yd5TXCzpyTm3DxfLFsxExKykFyRtSTT/z4vheB+0vejc6AzJAwAAALAWjkXEzjV43msj4ojtjZK+KOndkj7dKUwPEwAAAIBqwtUvZ3dE0iVzbm8rli2YsT0o6RxJxxdtcsSR4t9Tkj6j9tC/juhhAlZZRImjI5v5rJv51bpVJps7g3crmZOkZuR/q2mUyE5FfpfWKJWtlcimo2pEfkPUy2zgFTgAd1b55y/zd5V7vcpsh/z2Lfe+yb8fy7zPs5+f7OexnU1HS+0/yuyXSu3vAPSmlf2Y3y9ph+1L1S6Mrpb0q/My+yRdJ+lvJb1T0r2xyM6nKKo2R8Qx23VJb5P01cUaQcEEAAAAoJLs9OBVRMSs7Rsl3S2pJunWiHjY9ockHYyIfZJukXS77UOSTqhdVLXbZj8haZOkIdtvl3SlpO9LursolmpqF0ufXKwdFEwAAAAAqlnhjuSI2C9p/7xlN8+5PiXpqg6P3d5htW8o0waOYQIAAACADuhhAgAAAFBerOyQvG5BwQQAAACgGgomAAAAAOigDwomjmECAAAAgA7oYQIAAABQST8cw0QPEwAAAAB0QA8TAAAAgGr6oIeJgglYbc1mOjowk8/WGvk91kDD6ayT2dnZWnqdZxr1fLaZz55uDqezk62hfDby2VPRSGfrMZvOqlVivc6/b7Ia0UpnJ0tkT0X+a6jMdiizfcu8b8q8H8u8z7Ofn+znUZIG8m+ZcvuPEvulMvs7AD2oT6YVZ0geAAAAAHRADxMAAACAavqgh4mCCQAAAEA1FEwAAAAA8HJWfxzDRMEEAAAAoJo+KJiY9AEAAAAAOqCHCQAAAEB5fTKtOAUTAAAAgGoomAAAAACggz4omDiGCQAAAAA6oIcJWGWt06fT2drEdDo7eGZDfr351Wpg2qnc7HQtvc7TM/V09uTMaD47NJLOnhgcT2fHBmbS2bqb6aw0lU42PFuiDa0SbUg+f4lfEE9F/qvl+VZ+mz3fHEtnTzTz2/fkbL4NZd6PZd7nzeTnZzD5eZTKfc4Hz+Q3cJn9UrPE/g5Ab+IYJgAAAADohIIJAAAAABYQomACAAAAgE76YUgekz4AAAAAQAf0MAEAAACopg96mCiYAAAAAFTSD0PyKJgAAAAAVEPBBAAAAAAL6JNZ8pj0AQAAAAA6oIcJAAAAQGkuLusdBROwyg607kpnd1/4vnS2/sqN+exkvnN58ExuVzg7VUuv8/TUUDr7/MhoOnt0aDydHa010tkR57NlNGolXjPPpLN1N6s0Z1GNyLd1Mkps3+ZYOnt0Nv8eP9bIZ4/O5N83z0/n349l3udKfn6yn0dJqk/mx8nUT82mszr6XDpaZn8HoEf1wZA8CiYAAAAAlfTDLHln/ZnZ9q22n7X90Jxl59k+YPu7xb/nrmwzAQAAAGD1ZcblfErS7nnLbpL0tYjYIelrxW0AAAAA/SSWcOkRZy2YIuIbkk7MW7xH0m3F9dskvX15mwUAAACg6/VBwVT1GKYLI+Kp4vrTki5cpvYAAAAA6AXBMUwpEbFojWj7BtsHbR88evToUp8OAAAAQLfogx6mqgXTM7YvkqTi32c7BSNib0TsjIidW7durfh0AAAAALD6qhZM+yRdV1y/TtKXl6c5AAAAAHqFo/qlV2SmFf+spL+V9K9sH7Z9vaSPSNpl+7uSfq64DQAAAKCf9MGQvLNO+hAR13S46y3L3BYAAAAAPaSXeoqqqjpLHoBVEKdPp7P1iUaJbD2dHZxwKlcbq6XXOT0ylM4+NzSazo4Mbkxn626ls2VMRf61nWzlX4exgZl0tu7ZdDarEfmvizJ/14nmeDp7rJHfvk9NnZPOPnsmv97nTuffj9OT+dehNpH7/AxOpFep+kT+fzFl9h9l9ksAsB5QMAEAAAAor8eG1lVFwQQAAACgGgomAAAAAHg5i2OYAAAAAKCzPiiYqp6HCQAAAADWPQomAAAAAJU4ovIltX57t+3HbB+yfdMC9w/b/lxx/322txfLt9j+uu0J23867zFvsP2d4jF/bHvRKYEpmAAAAACUt5ST1ibqJds1SR+X9POSLpd0je3L58Wul/RcRFwm6WOSPlosn5L0QUm/ucCqPyHp1yTtKC67F2sHBRMAAACAShzVLwlXSDoUEY9HxIykOyXtmZfZI+m24voXJL3FtiNiMiK+qXbh9MP22hdJ2hQRfxcRIenTkt6+WCMomAAAAABUs4I9TJIulvTknNuHi2ULZiJiVtILkracZZ2Hz7LOl2CWPAAAAABr4XzbB+fc3hsRe9esNR1QMAFdrDU5mc7Wjk+ksyObh9PZxviix0H+f83RXE6Spobzu56J+kg6++xAK50toxH5zvgzzXo6+8LghnR2Q206nR3xbDqbNRX5bXa6mX9/nZzNb9+jM+Pp7LNnNuazE/n1Tkzk2+tT+dds6GTu8zP8Qn7+3pHnmulsmf1Hs8R+CcD6t8TzMB2LiJ2L3H9E0iVzbm8rli2UOWx7UNI5ko6fZZ3bzrLOl2BIHgAAAIBqVnZI3v2Sdti+1PaQpKsl7ZuX2SfpuuL6OyXdWxybtHBzI56SdNL2G4vZ8d4j6cuLNYIeJgAAAADl5SdvqLb6iFnbN0q6W1JN0q0R8bDtD0k6GBH7JN0i6XbbhySdULuokiTZfkLSJklDtt8u6cqIeETS+yR9StKopK8Ul44omAAAAAB0pYjYL2n/vGU3z7k+JemqDo/d3mH5QUmvy7aBggkAAABANSvYw9QtKJgAAAAAlGat7JC8bkHBBAAAAKCazvMrrBsUTAAAAAAq6YceJqYVBwAAAIAO6GECAAAAUF7+fEo9jYIJAAAAQCVurXULVh4FE9DFDrTuSmevHL42nR0ZG0lnG2MbU7nZkVp6na16fjTwzMBQOvtCOik1W/k2TM3W09mTw6Pp7KahM+nsaK2RztZX4NurEfnX60yzxOs1k3+9np/OZ587nc9OTOQ/D/F8/v049Hz+NRt6PpcbOZHftiNPT6azrScOp7Nl9ksA+gA9TAAAAACwMCZ9AAAAAIA+Rg8TAAAAgPJCnIcJAAAAADrphyF5FEwAAAAAqumDgoljmAAAAACgA3qYAAAAAJRmMSQPAAAAABYWwaQPAAAAANAJPUwAAAAA0AkFE4Becc/0Hens7q3vTWdHx4dSuebwSHqdUSsz30w+O9PKtVWSTjby652eye8qT44Mp7MbhsbS2dF6I52tuZXOZjUj/3qdadTT2dMzJbJT+e07PZnP+lR++w49n38dhk84nR09nttmo89Mp9fpJ59OZ+8usf8AgH5DwQQAAACgEobkAQAAAMBCQlJr/VdMFEwAAAAAqln/9RIFEwAAAIBq+mFIXpkjrwEAAACgr9DDBAAAAKAaTlwLAAAAAAvrhyF5FEwAAAAAygv1xaQPHMMEAAAAAB3QwwQAAACgNEsyxzABWI+ax0+ks/VnNqdyG4ZqJVpQTyfdzHeEDzTybZiZzq93eiq/3umRoXT25HAznR0czGcHBpb/y6vVcjo7O5t/vZrTJd43JbZDbSKfHTqZ/9uGnk9HNXq8lc5ueLqRytWfOZleZ5nPOQBUlt/V9awlFUy2n5B0SlJT0mxE7FyORgEAAADofvQw5fxMRBxbhvUAAAAA6BVM+gAAAAAA/W2pBVNIusf2A7ZvWI4GAQAAAOgF0T5xbdVLj1jqkLw3RcQR2xdIOmD7nyLiG3MDRSF1gyS96lWvWuLTAQAAAOgW/XDi2iX1MEXEkeLfZyV9SdIVC2T2RsTOiNi5devWpTwdAAAAgG7SBz1MlQsm22O2N754XdKVkh5aroYBAAAA6GIhuVX90iuWMiTvQklfsv3iej4TEX+1LK0CAAAAgC5QuWCKiMcl/cQytgUAAABAL+mhoXVVLcd5mAAAAAD0o/VfL1EwAf3oQOuudPbKoWtSueESzz8wsymdrU3n1zw1lT8ss3bG6ezsZH5XOTtaS2dbw/lvmZl6PhsDy//t5Vb+9XIjnx2cLpEtsc0GJ9JRDb+Qf71GTuQH3Y8+M53O1p85mcq1/uX/ptdZ5nMOAFW5D3qYOHEtAAAAAHRADxMAAACAavqgh4mCCQAAAEB5IamHpgevioIJAAAAQGlW9MUxTBRMAAAAAKrpg4KJSR8AAAAAoAN6mAAAAABUQw8TAAAAACzgxUkfql4SbO+2/ZjtQ7ZvWuD+YdufK+6/z/b2Off9VrH8MdtvnbP8Cdvfsf2g7YNnawM9TAAAAAAqWclJH2zXJH1c0i5JhyXdb3tfRDwyJ3a9pOci4jLbV0v6qKR32b5c0tWSXivplZK+avs1EdEsHvczEXEs0w4KJgCLumfms6ncroGr0uusnzgvnR2ceEV+vZNj+exELZ1tjDufHctnm8P5bKuejipWYOyAS0wbO9DIZ2vT+Wx9Mv+lXJ/IZ0eea5499GL26cl01k8+nc42j59I5Q607kqvEwBWxcoOybtC0qGIeFySbN8paY+kuQXTHkm/W1z/gqQ/te1i+Z0RMS3pX2wfKtb3t2UbwZA8AAAAAN3oYklPzrl9uFi2YCYiZiW9IGnLWR4bku6x/YDtG87WCHqYAAAAAFQQS+1hOn/eMUR7I2LvEhuV8aaIOGL7AkkHbP9TRHyjU5iCCQAAAEB5oaUWTMciYuci9x+RdMmc29uKZQtlDtselHSOpOOLPTYiXvz3WdtfUnuoXseCiSF5AAAAAKpZ2Vny7pe0w/altofUnsRh37zMPknXFdffKeneiIhi+dXFLHqXStoh6e9tj9neKEm2xyRdKemhxRpBDxMAAACArhMRs7ZvlHS3pJqkWyPiYdsfknQwIvZJukXS7cWkDifULqpU5D6v9gQRs5J+PSKati+U9KX2vBAalPSZiPirxdpBwQQAAACgkpWcVlySImK/pP3zlt085/qUpAWn6o2I35P0e/OWPS7pJ8q0gYIJAAAAQDUrXDB1AwomAAAAAOWFpBYFEwAAAAAsYMnTivcEZskDAAAAgA7oYQKwLA607lqR9V45fG06Ozq5LZ0d2jKezjbG6/nsxvxudXbU6Wyzns9GLR1NczOfrTXyvzYOnsln66dm89mJRjpbOz6RzraeOJzO3j19RzoLAD2rD3qYKJgAAAAAVEPBBAAAAAALYNIHAAAAAOgkpGitdSNWHJM+AAAAAEAH9DABAAAAqIZjmAAAAABgARzDBAAAAACL6IMeJo5hAgAAAIAO6GECAAAAUE0f9DBRMAEAAACoICiYAGCt3TN9Rzq7a+CqdHbgB2Pp7MiGDfns1nPT2eb4cDrbGqqls6o5n81q5r8QB2aa6WxtYjrfhqPPpaNx+nQ625ycTGcPtO5KZwFg3QtJrfV/HiYKJgAAAADV9EEPE5M+AAAAAEAH9DABAAAAqKYPepgomAAAAABUEJy4FgAAAAAWFFLE+p/0gWOYAAAAAKADepgAAAAAVMOQPAAAAADogEkfAAAAAGABEZy4FgAAAAA6oocJAHrHgdZda90E7Rq4Kp0d2LAhn63V0lnb6WxWlPlCbDbz0dOn09lu2L4AgP5DwQQAAACgkuiDIXlLmlbc9m7bj9k+ZPum5WoUAAAAgG4X7SF5VS89onIPk+2apI9L2iXpsKT7be+LiEeWq3EAAAAAulSIacXP4gpJhyLicUmyfaekPZIomAAAAIB+EAzJW8zFkp6cc/twsewlbN9g+6Dtg0ePHl3C0wEAAADA6lrSMUwZEbE3InZGxM6tW7eu9NMBAAAAWAUhKVpR+dIrljIk74ikS+bc3lYsAwAAALDeRfTFkLylFEz3S9ph+1K1C6WrJf3qsrQKAAAAQNfrpZ6iqioXTBExa/tGSXdLqkm6NSIeXraWAQAAAMAaW9KJayNiv6T9y9QWAAAAAL2kD4bkOVbxpFG2j0r6/qo9YX87X9KxtW4E0thevYdt1nvYZr2HbdZb2F5L8yMR0VMzpNn+K7W3e1XHImL3crVnpaxqwYTVY/tgROxc63Ygh+3Ve9hmvYdt1nvYZr2F7YX1asWnFQcAAACAXkXBBAAAAAAdUDCtX3vXugEohe3Ve9hmvYdt1nvYZr2F7YV1iWOYAAAAAKADepgAAAAAoAMKpnXE9lW2H7bdsr1z3n2/ZfuQ7cdsv3Wt2oiXs7272C6HbN+01u3By9m+1fazth+as+w82wdsf7f499y1bCN+yPYltr9u+5Fin/j+YjnbrEvZHrH997b/odhm/71Yfqnt+4r94+dsD611W/FDtmu2v237fxe32V5YlyiY1peHJP17Sd+Yu9D25ZKulvRaSbsl/U/btdVvHuYrtsPHJf28pMslXVNsL3SXT6n92ZnrJklfi4gdkr5W3EZ3mJX0GxFxuaQ3Svr14nPFNute05J+NiJ+QtLrJe22/UZJH5X0sYi4TNJzkq5fuyZiAe+X9Oic22wvrEsUTOtIRDwaEY8tcNceSXdGxHRE/IukQ5KuWN3WoYMrJB2KiMcjYkbSnWpvL3SRiPiGpBPzFu+RdFtx/TZJb1/NNqGziHgqIr5VXD+l9n/oLhbbrGtF20Rxs15cQtLPSvpCsZxt1kVsb5P07yT9WXHbYnthnaJg6g8XS3pyzu3DxTKsPbZN77owIp4qrj8t6cK1bAwWZnu7pJ+UdJ/YZl2tGN71oKRnJR2Q9D1Jz0fEbBFh/9hd/lDSf5XUKm5vEdsL6xQFU4+x/VXbDy1woVcCWCPRnm6UKUe7jO1xSV+U9IGIODn3PrZZ94mIZkS8XtI2tXvff2xtW4RObL9N0rMR8cBatwVYDYNr3QCUExE/V+FhRyRdMuf2tmIZ1h7bpnc9Y/uiiHjK9kVq/yqOLmG7rnaxdEdE/EWxmG3WAyLiedtfl/RvJW22PVj0WrB/7B4/LemXbP+CpBFJmyT9kdheWKfoYeoP+yRdbXvY9qWSdkj6+zVuE9rul7SjmFloSO3JOfatcZuQs0/SdcX16yR9eQ3bgjmKYylukfRoRPzBnLvYZl3K9lbbm4vro5J2qX3s2dclvbOIsc26RET8VkRsi4jtan9v3RsR14rthXWKE9euI7bfIelPJG2V9LykByPircV9/03Sf1R79qgPRMRX1qqdeKniF7o/lFSTdGtE/N7atgjz2f6spDdLOl/SM5J+R9L/kvR5Sa+S9H1JvxIR8yeGwBqw/SZJ/0fSd/TD4yt+W+3jmNhmXcj2j6s9SUBN7R9zPx8RH7L9arUnwzlP0rcl/YeImF67lmI+22+W9JsR8Ta2F9YrCiYAAAAA6IAheQAAAADQAQUTAAAAAHRAwQQAAAAAHVAwAQAAAEAHFEwAAAAA0AEFEwAAAAB0QMEEAAAAAB1QMAEAAABAB/8Pb2ph4rEFSNQAAAAASUVORK5CYII=)
%% 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 Simple
cm_method_params_expl_force = {
'method' : 'cumulant',
'stencil': stencil,
'relaxation_rates': [viscous_rr], # Specify viscous relaxation rate only
'force_model': Simple(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: execute_result
<matplotlib.colorbar.Colorbar at 0x7f4f58d67490>
<matplotlib.colorbar.Colorbar at 0x7f2ea7881df0>
%%%% Output: display_data
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAA0wAAAFoCAYAAABkP3u+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAArbElEQVR4nO3dfZBcV33m8eeZnp4XzUiWLcvGWCYysdisYRNSqBy2wh8kREbJkshs7GDHC65aV8wWcRXUJrvlZAsnS5EqqNqEJBuWKhM7GJfBYAiLaiNiC0yKJZU4lsEJfomDcMxawi96sS3NSDPT0/3bP/o6jMfTo9+989Y9/f1Udan79jNnzvTtvqPfnHPPdUQIAAAAAPBKA2vdAQAAAADoVhRMAAAAANABBRMAAAAAdEDBBAAAAAAdUDABAAAAQAcUTAAAAADQAQUTAAAAgK5ke7ftx20ftH3TAs8P2/5c8fz9trcX2y+z/VBx+3vb78y2+YrvwXWYAAAAAHQb2zVJ/yRpl6RDkh6QdE1EPDon8z5JPx4R/8n21ZLeGRHvsr1B0kxEzNq+QNLfS3q1pDhTm/MxwgQAAACgG10m6WBEPBERM5LukrRnXmaPpNuL+1+Q9DbbjohTETFbbB9Ru1DKtvkyg8vwg6Sde+65sX379tX8lgCwqv7pwSfSWQ+s0N+svAJtrtBkhGi10tnXvem1K9MJAOgCDz744NGI2LrW/Sjj7T8zFseONyt//YP/MP2IpKk5m26JiFvmPL5Q0lNzHh+S9FPzmvmXTDGa9KKkLZKO2v4pSbdJ+hFJ7y6ez7T5MqtaMG3fvl0HDhxYzW8JAKtq18BV6ezA6IZ8w7VaOmovf8VUavp2M//Ls3XqVDq7/8Dd+T4AQI+x/f217kNZR483df892yp/ff2C701FxM5l7NLLRMT9kl5v+19Lut32V6q0s6oFEwAAAID1ItSM/EyBCg5LumjO423FtoUyh2wPSjpL0rGX9TLiMdsTkt6QbPNlOIcJAAAAQGkhqaWofEt4QNIO2xfbHpJ0taS98zJ7JV1X3L9S0n0REcXXDEqS7R+R9GOSnky2+TKMMAEAAADoOsU5RzdKukdSTdJtEfGI7Q9JOhAReyXdKukO2wclHVe7AJKkt0i6yXZDUkvS+yLiqCQt1OZi/aBgAgAAAFBJSys6JU8RsU/Svnnbbp5zf0rSK04gjog7JN2RbXMxFEwAAAAASguFmn1wTVcKJgAAAACVJM9F6mks+gAAAAAAHTDCBAAAAKC0kNTsgxEmCiYA60api8aOjaWz3pC/wGzt0tels83x4XS2NZS/cK1qy3/hWjXzvxAHZvIXrq1NTKezu89/XzobJS6I25qcTGf3t7h4LgDM1Q9T8iiYAAAAAJQWEos+AAAAAEAnK7uoeHdg0QcAAAAA6IARJgAAAAClhYJFHwAAAABgQVFqTaCeRcEEAAAAoLRQf5zDRMEEAAAAoAKrqRW4lEWXYdEHAAAAAOiAESYAAAAApYWkFucwAQAAAMDC+mFKHgUTgK52+fC16Wxtx2vT2eaW8XS2MV7PZzfmD6uzo/lfMs16Phu1dDTNzXy21sj/uXHw9IZ0tv7qjfnsRCOdrR2bSGfLvB/vnb4znQWAXhTqj4KJc5gAAAAAoANGmAAAAABU0or1P8JEwQQAAACgtH6ZkkfBBAAAAKC0kNXsgzN8KJgAAAAAVNIPU/LWf0kIAAAAABUxwgQAAACgNM5hAgAAAICOrGas/wlrFEwAAAAASgtJrT44w4eCCQAAAEAlTMkDgKRdA1els7Ut56Szfv0l6ezpV42ls1Nn19LZxnj+l0FjLJ9tDqejatXz2ZWYHeFWPjvQyL8Gtel8u/XJ/A9Wn8i/YCOb8ztiZGwknd299b3pbPPY8VRuf+vudJsAgOVBwQQAAACgtAjOYQIAAACAjlpMyQMAAACAV2ovK77+R5jW/08IAAAAABUxwgQAAACgAs5hAgAAAIAFcR0mAAAAAFhEM1j0AQAAAABeIWQWfQAAAACAfsYIEwAAAIBKWiz6AKDfXT50TSpX2/HadJuN8zels6fPH05np87JH7Snz8rPuZ4dT0c1OxrpbGs4n416iexAPpvlVv71ciOfHZjOZwdPl8hO5LON8RLZsY3p7Oj4UDpbf3ZzKpf9PErSvTOfTWcBoIp+uQ4TBRMAAACA0kLui0UfzlgS2r7I9tdtP2r7EdvvL7b/ru3Dth8qbr+w8t0FAAAAgNWTGWGalfQbEfEt2xslPWh7f/HcxyLif6xc9wAAAAB0K67DJCkinpb0dHH/pO3HJF240h0DAAAA0L0ipGYfLPpQ6ie0vV3ST0q6v9h0o+1/sH2b7bOXu3MAAAAAupXVWsKtV6QLJtvjkr4o6QMRcULSJyT9qKQ3qj0C9fsdvu4G2wdsHzhy5MjSewwAAABgzYXaI0xVb70i1VPbdbWLpTsj4s8lKSKejYhmRLQkfVLSZQt9bUTcEhE7I2Ln1q1bl6vfAAAAALDizngOk21LulXSYxHxB3O2X1Cc3yRJ75T08Mp0EQAAAEA34jpMbT8t6d2SvmP7oWLbb0u6xvYb1R6Ne1LSe1egfwAAAAC6UMhq9cF1mDKr5H1TWvCsrH3L3x0AAAAAvYIRJgDr0q6Bq9LZ2o7XpnLTF+UXyjz1qno6e3pL/kA8szkd1cymSGeb4818wyP5bG04n60P5rMDA/mfLavVyv8FcXa2ls9Ol8hO5bO1sXy2OVriZxsp0e7wSDq7YSjX7nC6xXKf8/2tu0u0DABtIam1wos32N4t6Y8k1ST9aUR8ZN7zw5I+LelNko5JeldEPGl7l6SPSBqSNCPpv0TEfcXX/JWkCySdLpq5PCKe69QHCiYAAAAAXcd2TdLHJe2SdEjSA7b3RsSjc2LXS3o+Ii6xfbWkj0p6l6Sjkn4xIn5g+w2S7tHLryV7bUQcyPRj/Y+hAQAAAFgBVnMJt4TLJB2MiCciYkbSXZL2zMvskXR7cf8Lkt5m2xHx7Yj4QbH9EUmjxWhUaRRMAAAAAEp7aUpe1Zukc1+6Xmtxu2Het7hQ0lNzHh/Sy0eJXpaJiFlJL0raMi/zy5K+FRHTc7b9me2HbH+wWBW8I6bkAQAAAKgkOVLUydGI2LlcfVmI7derPU3v8jmbr42Iw7Y3qn2t2XerfR7UghhhAgAAANCNDku6aM7jbcW2BTO2ByWdpfbiD7K9TdKXJL0nIr730hdExOHi35OSPqP21L+OKJgAAAAAlBbhpU7JO5MHJO2wfbHtIUlXS9o7L7NX0nXF/Ssl3RcRYXuzpL+QdFNE/PVLYduDts8t7tclvUPSw4t1gil5AAAAACppruCy4hExa/tGtVe4q0m6LSIesf0hSQciYq+kWyXdYfugpONqF1WSdKOkSyTdbPvmYtvlkiYl3VMUSzVJX5X0ycX6QcEEAAAAoLSQ1FraOUxn/h4R+yTtm7ft5jn3pyS94sJzEfFhSR/u0OybyvSBggkAAABABV7REaZuQcEE9KHalnPS2cb5m1K5U6+qp9s8dV7+4Dp9TqSzM5tb6WxsnE1nh8dm0tkNIyWyQ410drSez9acfx2yyvxCPN0o8V6YKZGdGkpnp0fy2anh/K/CVj3/OkStzH8icq/DwEzu8yhJ9eP5zzkAoDMKJgAAAAClta/DtLJT8roBBRMAAACASpp9sOg2BRMAAACA0kJmhAkAAAAAOmn1wQjT+v8JAQAAAKAiRpgAAAAAlBYhNZmSBwAAAAAL4xwmAAAAAFhAe9GH9X+Gz/r/CQEAAACgIkaYAAAAAFTSFFPyAPSIy4evTWf9+kvS2dPnD+dyW/ID1tPnRDo7c04rnfXmmXR24/hUOnv2htPp7ObhfHbTUD47Wmuks3XnX7OsRokpF6eb9XT2xMxoOvvCSD77/FA+O1EfSWdnBobS2TKTONzMZWvTuc+jJA1OvCqdLXP8uHf6znQWwPoW4hwmAAAAAOigP85homACAAAAUEmrD6bkrf+SEAAAAAAqYoQJAAAAQGlcuBYAAAAAFsE5TAAAAACwgPaFaxlhAgAAAIAFsegDAAAAAPQxRpgAAAAAlMaFawEAAABgESz6AGBN7Rq4Kp2t7XhtOnv6VWPp7NQ5uQPhzOZ0k5rZ3EpnvXkmnT1r06l09rzxiXx29GQ6u3Uo3+6mwal0dkNtOp0d8Ww6mzUV+V8Xp5rD6eyJoZF09sjQeDo7MrgxnX1uIP9+fDGdlGZaQ+nsQKOWyk1N5f9jUp/Mf85HJ7els2WOS/tbd6ezAHpQ9MeiD+u/JAQAAACAihhhAgAAAFBaqD9WyaNgAgAAAFBJP0zJo2ACAAAAUBqr5AEAAADAIvqhYGLRBwAAAADogBEmAAAAAKWF+mNZcQomAAAAAJWwSh4AAAAALCT64xwmCiYAAAAApbFKHoA1NzA2ls42t4yns1Nn19LZ6bNyB8KZTZFuMzbOprMbx6fS2fPGJ9LZCze8mM5eMJLPnls/mc6eU8v3d2xgJp2tO//6ZjUi/+tisjWUzh4fzL9vR2uNdLbuVjpbRrOVXyvpRCOfnZnOZWun8/8xqU/kP+dDJY4fAz/IH5cAYD2gYAIAAABQCSNMAAAAALCAflkl74xzAGxfZPvrth+1/Yjt9xfbz7G93/Z3i3/PXvnuAgAAAOgWEa586xWZSdOzkn4jIi6V9GZJv277Ukk3SfpaROyQ9LXiMQAAAACsG2csmCLi6Yj4VnH/pKTHJF0oaY+k24vY7ZKuWKE+AgAAAOhCLbnyrVeUOofJ9nZJPynpfknnR8TTxVPPSDq/w9fcIOkGSXrNa15TuaMAAAAAukf0yXWY0mue2h6X9EVJH4iIE3Ofi4hQeyn2V4iIWyJiZ0Ts3Lp165I6CwAAAKB79MM5TKkRJtt1tYulOyPiz4vNz9q+ICKetn2BpOdWqpMAAAAAug2r5EmSbFvSrZIei4g/mPPUXknXFfevk/Tl5e8eAAAAAKydzAjTT0t6t6Tv2H6o2Pbbkj4i6fO2r5f0fUm/siI9BAAAANCVemlqXVVnLJgi4ptSx2Us3ra83QEwlzdsSGcb4/US2fzBbXY8l2uON9NtDo/NpLNnbzidzp43ejKdvWDkxXR229DxdHbrYL4Pm2uT6eyY869Z3fl9kdWIWjo7GUPp7NhA/ucacSOdLaMR6dN5NTWb/5xNz+TXVZqeyr2+s5P5Nst8zsscP0ZKHJcArG+h/lj0odQqeQAAAAAgSYr2SnnrXf7PagAAAAAwx0pfh8n2btuP2z5o+6YFnh+2/bni+fuLyyDJ9i7bD9r+TvHvz875mjcV2w/a/uNizYaOKJgAAAAAdB3bNUkfl/Tzki6VdI3tS+fFrpf0fERcIuljkj5abD8q6Rcj4t+ovUDdHXO+5hOSfk3SjuK2e7F+UDABAAAAKC204tdhukzSwYh4IiJmJN0lac+8zB5Jtxf3vyDpbbYdEd+OiB8U2x+RNFqMRl0gaVNE/G1xLdlPS7pisU5wDhMAAACACpZ8HaZzbR+Y8/iWiLhlzuMLJT015/EhST81r41/yUTErO0XJW1Re4TpJb8s6VsRMW37wqKduW1euFgnKZgAAAAAVLLERR+ORsTOZerKgmy/Xu1pepdXbYMpeQAAAAC60WFJF815vK3YtmDG9qCksyQdKx5vk/QlSe+JiO/NyW87Q5svQ8EEAAAAoJIVPofpAUk7bF9se0jS1ZL2zsvsVXtRB0m6UtJ9ERG2N0v6C0k3RcRf/7C/8bSkE7bfXKyO9x5JX16sExRMAAAAAEqLWNmCKSJmJd0o6R5Jj0n6fEQ8YvtDtn+piN0qaYvtg5L+s6SXlh6/UdIlkm62/VBxO6947n2S/lTSQUnfk/SVxfrBOUwAAAAAKlniog9nFBH7JO2bt+3mOfenJF21wNd9WNKHO7R5QNIbsn2gYAJW2a6BV3ymO6pd+rp0trEx/3FujOUPbrOjybM5R5rpNjeMzKSzm4dPp7NbhybS2XPrJ/PtDpbJnkhnNw9MpbMbPZvO1lfgd1ejxEm9J6ORztadf9+UMRX1dPZ0M589MTyaz44Mp7PTI0Op3OxoLd1mmc95mePHyNaz09kyx7v9rbvTWQDdY4mLPvQEpuQBAAAAQAeMMAEAAACoJLl4Q0+jYAIAAABQWii92l1Po2ACAAAAUEkfnMJEwQQAAACgguiPKXks+gAAAAAAHTDCBAAAAKCaPpiTR8EEAAAAoJJ+mJJHwQQAAACgEi5cCwAAAAB9jBEmYJUNbNiQzjbHh9PZ2dH8kHgz36xaw7k/HdWGm+k2Nww10tlNQ6fz2cGpdPac2kQ6u7k2mc8OlOjDwGw6O+b837fqJbJZjWjlv3/kfy4p/3o1arV0drI1lM6+OJj/TJZ5P24YGktnTyQ/P9nPoyQ1h/PHhFLHjxLHpTLHOwC9J8SUPAAAAABYWEiiYAIAAACAhfXDOUwUTAAAAACq6YOCiUUfAAAAAKADRpgAAAAAVGAWfQAAAACAjvpgSh4FEwAAAIDyoj+WFeccJgAAAADogBEmAAAAANUwJQ8AAAAAOln/U/IomIDVVqulo62hfLZZzx+wWvV0VFHP/emoPthMtzlab+SztXx2Q206nR0bmMlnnc9u9GyJdvOzojcM5HfaoPLvm6y68/tXrfw+a5R4vU6V2A9l9m+Z902Z92OZ9/lg8vMzk/w8SuU+56WOHyWOSwMljncAehQjTAAAAADQQR8UTCz6AAAAAAAdMMIEAAAAoLyQ1AfLilMwAQAAAKgk+mBKHgUTAAAAgGoomAAAAACggz6YkseiDwAAAADQASNMAAAAACoxU/IAAAAAYAEhzmECAAAAgIW5L85homACVpld4sBSy2ejlm82Spy9GAO5Px0NJHOSVHMrna2XyI54tkS7ZbLNEtl0VHXnd8Sg8ju4VqLdtBJ/QSz3epV5L5RpN79/y71v8v0t8z7Pfn6yn8d2Nh0tdfwoc1wqdbwDgC5FwQQAAACgmj6YknfGvz/Zvs32c7YfnrPtd20ftv1QcfuFle0mAAAAgK4TS7j1iMyA/ack7V5g+8ci4o3Fbd/ydgsAAABA1+uDgumMU/Ii4hu2t69CXwAAAAD0ilBfLPqwlDODb7T9D8WUvbM7hWzfYPuA7QNHjhxZwrcDAAAAgNVVtWD6hKQflfRGSU9L+v1OwYi4JSJ2RsTOrVu3Vvx2AAAAALqNo/qtV1QqmCLi2YhoRkRL0iclXba83QIAAADQ9frgHKZKBZPtC+Y8fKekhztlAQAAAKBXnXHRB9uflfRWSefaPiTpdyS91fYb1a4Nn5T03pXrIgAAAIBu1EtT66rKrJJ3zQKbb12BvgAAAABAVzljwQQAAAAAC+qDZcUpmAAAAACU12OLN1RFwQQAAACgmj4omJZy4VoAAAAAWNcomAAAAABUstIXrrW92/bjtg/avmmB54dtf654/n7b24vtW2x/3faE7T+Z9zV/VbT5UHE7b7E+MCUPAAAAQDUrOCXPdk3SxyXtknRI0gO290bEo3Ni10t6PiIusX21pI9KepekKUkflPSG4jbftRFxINMPRpgAAAAAVBNLuJ3ZZZIORsQTETEj6S5Je+Zl9ki6vbj/BUlvs+2ImIyIb6pdOC0JBRMAAACA0pYyHa+Ykneu7QNzbjfM+xYXSnpqzuNDxbYFMxExK+lFSVsS3f+zYjreB20vujY6U/IAAAAArIWjEbFzDb7vtRFx2PZGSV+U9G5Jn+4UZoQJAAAAQDXh6rczOyzpojmPtxXbFszYHpR0lqRji3Y54nDx70lJn1F76l9HjDABqyyixNmRzXzWzXyzbpXJ5q7g3UrmJKkZ+b/VNEpkpyJ/SGuUytZKZNNRNSK/I+pldvAKnIA7q/z3L/NzlXu9yuyH/P4t977Jvx/LvM+zn5/s57GdTUdLHT/KHJdKHe8A9KaV/Zg/IGmH7YvVLoyulvSr8zJ7JV0n6W8kXSnpvljk4FMUVZsj4qjtuqR3SPrqYp2gYAIAAABQSXZ58CoiYtb2jZLukVSTdFtEPGL7Q5IORMReSbdKusP2QUnH1S6q2n2zn5S0SdKQ7SskXS7p+5LuKYqlmtrF0icX6wcFEwAAAIBqVnggOSL2Sdo3b9vNc+5PSbqqw9du79Dsm8r0gXOYAAAAAKADRpgAAAAAlBcrOyWvW1AwAQAAAKiGggkAAAAAOuiDgolzmAAAAACgA0aYAAAAAFTSD+cwMcIEAAAAAB0wwgQAAACgmj4YYaJgAlZbs5mODszks7VG/og10HA662R2draWbvN0o57PNvPZU83hdHayNZTPRj57MhrpbD1m01m1SrTr/PsmqxGtdHayRPZk5H8NldkPZfZvmfdNmfdjmfd59vOT/TxK0kD+LVPu+FHiuFTmeAegB/XJsuJMyQMAAACADhhhAgAAAFBNH4wwUTABAAAAqIaCCQAAAABeyeqPc5gomAAAAABU0wcFE4s+AAAAAEAHjDABAAAAKK9PlhWnYAIAAABQDQUTAAAAAHTQBwUT5zABAAAAQAeMMAGrrHXqVDpbm5hOZwdPb8i3m29WA9NO5Wana+k2T83U09kTM6P57NBIOnt8cDydHRuYSWfrbqaz0lQ62fBsiT60SvQh+f1L/AXxZOR/tbzQyu+zF5pj6ezxZn7/npjN96HM+7HM+7yZ/PwMJj+PUrnP+eDp/A4uc1xqljjeAehNnMMEAAAAAJ1QMAEAAADAAkIUTAAAAADQST9MyWPRBwAAAADogBEmAAAAANX0wQgTBRMAAACASvphSh4FEwAAAIBqKJgAAAAAYAF9skoeiz4AAAAAQAeMMAEAAAAozcVtvaNgAlbZ/tbd6ezu89+XztZfvTGfncwPLg+ezh0KZ6dq6TZPTQ2lsy+MjKazR4bG09nRWiOdHXE+W0ajVuI180w6W3ezSncW1Yh8XyejxP5tjqWzR2bz7/GjjXz2yEz+ffPCdP79WOZ9ruTnJ/t5lKT6ZH6eTP3kbDqrI8+no2WOdwB6VB9MyaNgAgAAAFBJP6ySd8Y/M9u+zfZzth+es+0c2/ttf7f49+yV7SYAAAAArL7MvJxPSdo9b9tNkr4WETskfa14DAAAAKCfxBJuPeKMBVNEfEPS8Xmb90i6vbh/u6QrlrdbAAAAALpeHxRMVc9hOj8ini7uPyPp/GXqDwAAAIBeEJzDlBIRi9aItm+wfcD2gSNHjiz12wEAAADoFn0wwlS1YHrW9gWSVPz7XKdgRNwSETsjYufWrVsrfjsAAAAAWH1VC6a9kq4r7l8n6cvL0x0AAAAAvcJR/dYrMsuKf1bS30j6V7YP2b5e0kck7bL9XUk/VzwGAAAA0E/6YEreGRd9iIhrOjz1tmXuCwAAAIAe0ksjRVVVXSUPwCqIU6fS2fpEo0S2ns4OTjiVq43V0m1Ojwyls88PjaazI4Mb09m6W+lsGVORf20nW/nXYWxgJp2tezadzWpE/tdFmZ/reHM8nT3ayO/fp6fOSmefO51v9/lT+ffj9GT+dahN5D4/gxPpJlWfyP8vpszxo8xxCQDWAwomAAAAAOX12NS6qiiYAAAAAFRDwQQAAAAAr2RxDhMAAAAAdNYHBVPV6zABAAAAwLpHwQQAAACgEkdUvqXat3fbftz2Qds3LfD8sO3PFc/fb3t7sX2L7a/bnrD9J/O+5k22v1N8zR/bXnRJYAomAAAAAOUt5aK1iXrJdk3SxyX9vKRLJV1j+9J5seslPR8Rl0j6mKSPFtunJH1Q0m8u0PQnJP2apB3Fbfdi/aBgAgAAAFCJo/ot4TJJByPiiYiYkXSXpD3zMnsk3V7c/4Kkt9l2RExGxDfVLpx+2F/7AkmbIuJvIyIkfVrSFYt1goIJAAAAQDUrOMIk6UJJT815fKjYtmAmImYlvShpyxnaPHSGNl+GVfIAAAAArIVzbR+Y8/iWiLhlzXrTAQUT0MVak5PpbO3YRDo7snk4nW2ML3oe5L9ojuZykjQ1nD/0TNRH0tnnBlrpbBmNyA/Gn27W09kXBzeksxtq0+nsiGfT2aypyO+zU838++vEbH7/HpkZT2efO70xn53Itzsxke+vT+Zfs6ETuc/P8Iv59XtHnm+ms2WOH80SxyUA698Sr8N0NCJ2LvL8YUkXzXm8rdi2UOaQ7UFJZ0k6doY2t52hzZdhSh4AAACAalZ2St4DknbYvtj2kKSrJe2dl9kr6bri/pWS7ivOTVq4uxFPSzph+83F6njvkfTlxTrBCBMAAACA8vKLN1RrPmLW9o2S7pFUk3RbRDxi+0OSDkTEXkm3SrrD9kFJx9UuqiRJtp+UtEnSkO0rJF0eEY9Kep+kT0kalfSV4tYRBRMAAACArhQR+yTtm7ft5jn3pyRd1eFrt3fYfkDSG7J9oGACAAAAUM0KjjB1CwomAAAAAKVZKzslr1tQMAEAAACopvP6CusGBRMAAACASvphhIllxQEAAACgA0aYAAAAAJSXv55ST6NgAgAAAFCJW2vdg5VHwQR0sf2tu9PZy4evTWdHxkbS2cbYxlRudqSWbrNVz88GnhkYSmdfTCelZivfh6nZejp7Yng0nd00dDqdHa010tn6Cvz2akT+9TrdLPF6zeRfrxem89nnT+WzExP5z0O8kH8/Dr2Qf82GXsjlRo7n9+3IM5PpbOvJQ+lsmeMSgD7ACBMAAAAALIxFHwAAAACgjzHCBAAAAKC8ENdhAgAAAIBO+mFKHgUTAAAAgGr6oGDiHCYAAAAA6IARJgAAAAClWUzJAwAAAICFRbDoAwAAAAB0wggTAAAAAHRCwQSgV9w7fWc6u3vre9PZ0fGhVK45PJJuM2pl1pvJZ2daub5K0olGvt3pmfyh8sTIcDq7YWgsnR2tN9LZmlvpbFYz8q/X6UY9nT01UyI7ld+/05P5rE/m9+/QC/nXYfi409nRY7l9NvrsdLpNP/VMOntPieMHAPQbCiYAAAAAlTAlDwAAAAAWEpJa679iomACAAAAUM36r5comAAAAABU0w9T8sqceQ0AAAAAfYURJgAAAADVcOFaAAAAAFhYP0zJo2ACAAAAUF6oLxZ94BwmAAAAAOiAESYAAAAApVmSOYcJwHrUPHY8na0/uzmV2zBUK9GDejrpZn4gfKCR78PMdL7d6al8u9MjQ+nsieFmOjs4mM8ODCz/L69Wy+ns7Gz+9WpOl3jflNgPtYl8duhE/mcbeiEd1eixVjq74ZlGKld/9kS6zTKfcwCoLH+o61kUTAAAAAAqYYTpDGw/KemkpKak2YjYuRydAgAAANDl+mTRh+UYYfqZiDi6DO0AAAAAQFdhSh4AAACACqIvLly71GXFQ9K9th+0fcNCAds32D5g+8CRI0eW+O0AAAAAdAtH9VuvWOoI01si4rDt8yTtt/2PEfGNuYGIuEXSLZK0c+fOHnppAAAAACyKEabFRcTh4t/nJH1J0mXL0SkAAAAAXS4kt6rfekXlgsn2mO2NL92XdLmkh5erYwAAAACw1pYyJe98SV+y/VI7n4mIv1yWXgEAAADofn0wJa9ywRQRT0j6iWXsCwAAAIBesv7rJZYVB/rR/tbd6ezlQ9ekcsMlvv/AzKZ0tjadb3lqKj/LuHba6ezsZP5QOTtaS2dbw/nfMjP1fDYGlv+3l1v518uNfHZwukS2xD4bnEhHNfxi/vUaOZ6fdD/67HQ6W3/2RCrX+uf/l26zzOccAKpyH4wwLXVZcQAAAABYtxhhAgAAAFBNH4wwUTABAAAAKC8k9dDy4FVRMAEAAAAozYq+OIeJggkAAABANX1QMLHoAwAAAAB0wAgTAAAAgGoYYQIAAACABby06EPVW4Lt3bYft33Q9k0LPD9s+3PF8/fb3j7nud8qtj9u++1ztj9p+zu2H7J94Ex9YIQJAAAAQCUrueiD7Zqkj0vaJemQpAds742IR+fErpf0fERcYvtqSR+V9C7bl0q6WtLrJb1a0ldtvy4imsXX/UxEHM30g4IJwKLunflsKrdr4Kp0m/Xj56SzgxOvyrc7OZbPTtTS2ca489mxfLY5nM+26umoYgXmDrjEsrEDjXy2Np3P1ifzv5TrE/nsyPPNM4deyj4zmc76qWfS2eax46nc/tbd6TYBYFWs7JS8yyQdjIgnJMn2XZL2SJpbMO2R9LvF/S9I+hPbLrbfFRHTkv7Z9sGivb8p2wmm5AEAAABYC+faPjDndsO85y+U9NScx4eKbQtmImJW0ouStpzha0PSvbYfXOB7vgIjTAAAAAAqiKWOMB2NiJ3L1ZsS3hIRh22fJ2m/7X+MiG90CjPCBAAAAKC8ULtgqno7s8OSLprzeFuxbcGM7UFJZ0k6ttjXRsRL/z4n6UtqT9XriIIJAAAAQDUru0reA5J22L7Y9pDaizjsnZfZK+m64v6Vku6LiCi2X12sonexpB2S/s72mO2NkmR7TNLlkh5erBNMyQMAAADQdSJi1vaNku6RVJN0W0Q8YvtDkg5ExF5Jt0q6o1jU4bjaRZWK3OfVXiBiVtKvR0TT9vmSvtReF0KDkj4TEX+5WD8omAAAAABUspLLiktSROyTtG/etpvn3J+StOBSvRHxe5J+b962JyT9RJk+UDABAAAAqGaFC6ZuQMEEAAAAoLyQ1KJgAgAAAIAFLHlZ8Z7AKnkAAAAA0AEjTACWxf7W3SvS7uXD16azo5Pb0tmhLePpbGO8ns9uzB9WZ0edzjbr+WzU0tE0N/PZWiP/18bB0/ls/eRsPjvRSGdrxybS2daTh9LZe6bvTGcBoGf1wQgTBRMAAACAaiiYAAAAAGABLPoAAAAAAJ2EFK217sSKY9EHAAAAAOiAESYAAAAA1XAOEwAAAAAsgHOYAAAAAGARfTDCxDlMAAAAANABI0wAAAAAqumDESYKJgAAAAAVBAUTAKy1e6fvTGd3DVyVzg78YCydHdmwIZ/denY62xwfTmdbQ7V0VjXns1nN/C/EgZlmOlubmM734cjz6WicOpXONicn09n9rbvTWQBY90JSa/1fh4mCCQAAAEA1fTDCxKIPAAAAANABI0wAAAAAqumDESYKJgAAAAAVBBeuBQAAAIAFhRSx/hd94BwmAAAAAOiAESYAAAAA1TAlDwAAAAA6YNEHAAAAAFhABBeuBQAAAICOGGECgN6xv3X3WndBuwauSmcHNmzIZ2u1dNZ2OpsVZX4hNpv56KlT6Ww37F8AQP+hYAIAAABQSfTBlLwlLStue7ftx20ftH3TcnUKAAAAQLeL9pS8qrceUXmEyXZN0scl7ZJ0SNIDtvdGxKPL1TkAAAAAXSrEsuJncJmkgxHxhCTZvkvSHkkUTAAAAEA/CKbkLeZCSU/NeXyo2PYytm+wfcD2gSNHjizh2wEAAADA6lrSOUwZEXFLROyMiJ1bt25d6W8HAAAAYBWEpGhF5VuvWMqUvMOSLprzeFuxDQAAAMB6F9EXU/KWUjA9IGmH7YvVLpSulvSry9IrAAAAAF2vl0aKqqpcMEXErO0bJd0jqSbptoh4ZNl6BgAAAABrbEkXro2IfZL2LVNfAAAAAPSSPpiS51jFi0bZPiLp+6v2DfvbuZKOrnUnkMb+6j3ss97DPus97LPewv5amh+JiJ5aIc32X6q936s6GhG7l6s/K2VVCyasHtsHImLnWvcDOeyv3sM+6z3ss97DPust7C+sVyu+rDgAAAAA9CoKJgAAAADogIJp/bplrTuAUthfvYd91nvYZ72HfdZb2F9YlziHCQAAAAA6YIQJAAAAADqgYAIAAACADiiY1hHbV9l+xHbL9s55z/2W7YO2H7f99rXqI17J9u5ivxy0fdNa9wevZPs228/ZfnjOtnNs77f93eLfs9eyj/gh2xfZ/rrtR4tj4vuL7eyzLmV7xPbf2f77Yp/992L7xbbvL46Pn7M9tNZ9xQ/Zrtn+tu3/Uzxmf2FdomBaXx6W9O8lfWPuRtuXSrpa0usl7Zb0v2zXVr97mK/YDx+X9POSLpV0TbG/0F0+pfZnZ66bJH0tInZI+lrxGN1hVtJvRMSlkt4s6deLzxX7rHtNS/rZiPgJSW+UtNv2myV9VNLHIuISSc9Lun7tuogFvF/SY3Mes7+wLlEwrSMR8VhEPL7AU3sk3RUR0xHxz5IOSrpsdXuHDi6TdDAinoiIGUl3qb2/0EUi4huSjs/bvEfS7cX92yVdsZp9QmcR8XREfKu4f1Lt/9BdKPZZ14q2ieJhvbiFpJ+V9IViO/usi9jeJunfSfrT4rHF/sI6RcHUHy6U9NScx4eKbVh77JvedX5EPF3cf0bS+WvZGSzM9nZJPynpfrHPuloxveshSc9J2i/pe5JeiIjZIsLxsbv8oaT/KqlVPN4i9hfWKQqmHmP7q7YfXuDGqASwRqJ9fQau0dBlbI9L+qKkD0TEibnPsc+6T0Q0I+KNkrapPfr+Y2vbI3Ri+x2SnouIB9e6L8BqGFzrDqCciPi5Cl92WNJFcx5vK7Zh7bFvetezti+IiKdtX6D2X8XRJWzX1S6W7oyIPy82s896QES8YPvrkv6tpM22B4tRC46P3eOnJf2S7V+QNCJpk6Q/EvsL6xQjTP1hr6SrbQ/bvljSDkl/t8Z9QtsDknYUKwsNqb04x9417hNy9kq6rrh/naQvr2FfMEdxLsWtkh6LiD+Y8xT7rEvZ3mp7c3F/VNIutc89+7qkK4sY+6xLRMRvRcS2iNiu9u+t+yLiWrG/sE65PSsB64Htd0r6n5K2SnpB0kMR8fbiuf8m6T+qvXrUByLiK2vVT7xc8Re6P5RUk3RbRPze2vYI89n+rKS3SjpX0rOSfkfS/5b0eUmvkfR9Sb8SEfMXhsAasP0WSf9X0nf0w/Mrflvt85jYZ13I9o+rvUhATe0/5n4+Ij5k+7VqL4ZzjqRvS/oPETG9dj3FfLbfKuk3I+Id7C+sVxRMAAAAANABU/IAAAAAoAMKJgAAAADogIIJAAAAADqgYAIAAACADiiYAAAAAKADCiYAAAAA6ICCCQAAAAA6+P/CKCVkBN3LaAAAAABJRU5ErkJggg==)
%% 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)
```
......
......@@ -250,8 +250,8 @@ def run(re=6000, eval_interval=0.05, total_time=3.0, domain_size=100, u_0=0.05,
def create_full_parameter_study(gpu=False):
"""Creates a parameter study that can run the Kida vortex flow with entropic, KBC, Smagorinsky and MRT methods."""
opt_cpu = {'target': 'cpu', 'openmp': 4}
opt_gpu = {'target': 'gpu'}
opt_cpu = {'target': ps.Target.CPU, 'openmp': 4}
opt_gpu = {'target': ps.Target.GPU}
mrt_one = [{'method': 'mrt3', 'relaxation_rates': ['viscosity', 1, 1], 'stencil': stencil}
for stencil in ('D3Q19', 'D3Q27')]
......
......@@ -16,7 +16,7 @@ from lbmpy.boundaries import FixedDensity, NoSlip
from lbmpy.geometry import add_pipe_inflow_boundary, add_pipe_walls
from lbmpy.lbstep import LatticeBoltzmannStep
from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
from pystencils import create_data_handling
from pystencils import create_data_handling, Target
from pystencils.slicing import make_slice, slice_from_direction
......@@ -135,4 +135,4 @@ def test_rod_scenario_simple():
if __name__ == '__main__':
# High Re (Entropic method)
rod_simulation(stencil_name=sys.argv[1], re=500, diameter=80, entropic=True, time_to_simulate=17,
parallel=False, optimization_params={'target': 'gpu'})
parallel=False, optimization_params={'target': Target.GPU})
......@@ -124,7 +124,7 @@ def schaefer_turek_2d(cells_per_diameter, u_max=0.3, max_lattice_velocity=0.05,
(u_max, domain_size, dx, dt, omega, re_lattice))
initial_velocity = get_pipe_velocity_field(domain_size, max_lattice_velocity)
scenario = create_channel(domain_size=domain_size, u_max=max_lattice_velocity, kernel_params={'omega_0': omega},
scenario = create_channel(domain_size=domain_size, u_max=max_lattice_velocity, relaxation_rate=omega,
initial_velocity=initial_velocity, **kwargs)
scenario.boundary_handling.set_boundary(NoSlip('obstacle'), mask_callback=geometry_callback)
parameter_info['u_bar_l'] = u_bar_l
......@@ -159,5 +159,5 @@ def test_schaefer_turek():
if __name__ == '__main__':
long_run(entropic=True, method='trt-kbc-n1', compressible=True,
optimization={'target': 'gpu', 'gpuIndexingParams': {'blockSize': (16, 8, 2)}})
long_run(entropic=True, method='trt_kbc_n1', compressible=True,
optimization={'target': Target.GPU, 'gpuIndexingParams': {'blockSize': (16, 8, 2)}})
......@@ -8,6 +8,7 @@ import numpy as np
import sympy as sp
import pytest
from pystencils import Target
from lbmpy.creationfunctions import update_with_default_parameters
from lbmpy.relaxationrates import (
......@@ -166,7 +167,7 @@ def create_full_parameter_study():
'relaxation_rates': ["omega", str(relaxation_rate_from_magic_number(sp.Symbol("omega")))],
'equilibrium_order': eqOrder,
'maxwellian_moments': mbEq,
'optimization': {'target': 'cpu', 'split': True, 'cse_pdfs': True}}
'optimization': {'target': Target.CPU, 'split': True, 'cse_pdfs': True}}
for method in ('srt', 'trt')
for stencil in ('D3Q19', 'D3Q27')
for comp in (True, False)
......@@ -222,7 +223,7 @@ def test_shear_wave():
'stencil': 'D3Q19',
'compressible': True,
"optimization": {"target": "gpu"}
"optimization": {"target": Target.GPU}
}
run(32, nu=1e-2, equilibrium_order=2, method='srt', y_size=1, periodicity_in_kernel=True,
relaxation_rates=[sp.Symbol("omega"), 5, 5], maxwellian_moments=True, **params)
......@@ -6,7 +6,7 @@ Silva, Semiao
python3 scenario_square_channel.py server
python3 scenario_square_channel.py client --host i10staff41 -P '{ "optimization" : { "target" : "gpu"} }'
python3 scenario_square_channel.py client --host i10staff41 -P '{ "optimization" : { "target" : Target.GPU} }'
"""
import numpy as np
......
......@@ -46,7 +46,7 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 3,
"metadata": {},
"outputs": [
{
......@@ -113,7 +113,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.6"
"version": "3.8.2"
}
},
"nbformat": 4,
......
......@@ -144,7 +144,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.4"
"version": "3.8.2"
}
},
"nbformat": 4,
......
%% Cell type:code id: tags:
``` python
import pytest
pytest.importorskip('scipy.ndimage')