phasefield-gravity-wave.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import time\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "import numpy as np\n",
    "import math\n",
    "\n",
    "from pystencils.session import *\n",
    "from lbmpy.session import *\n",
    "\n",
    "from pystencils.simp import sympy_cse\n",
    "from pystencils.boundaries import BoundaryHandling\n",
    "\n",
    "from lbmpy.phasefield_allen_cahn.parameter_calculation import AllenCahnParameters\n",
    "from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle\n",
    "from lbmpy.phasefield_allen_cahn.kernel_equations import *\n",
    "\n",
    "from lbmpy.advanced_streaming import LBMPeriodicityHandling\n",
    "from lbmpy.boundaries import NoSlip, LatticeBoltzmannBoundaryHandling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If `pycuda` is installed the simulation automatically runs on GPU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No pycuda installed\n"
     ]
    }
   ],
   "source": [
    "try:\n",
    "    import pycuda\n",
    "except ImportError:\n",
    "    pycuda = None\n",
    "    gpu = False\n",
    "    target = ps.Target.CPU\n",
    "    print('No pycuda installed')\n",
    "\n",
    "if pycuda:\n",
    "    gpu = True\n",
    "    target = ps.Target.GPU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gravity wave simulated with a phase-field model for immiscible fluids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Geometry Setup\n",
    "\n",
    "First of all the stencils for the phase-field LB step as well as the stencil for the hydrodynamic LB step are defined. According to the stencils the simulation runs either in 2D or 3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "stencil_phase = LBStencil(Stencil.D2Q9)\n",
    "stencil_hydro = LBStencil(Stencil.D2Q9)\n",
    "assert(stencil_hydro.D == stencil_phase.D)\n",
    "\n",
    "dimensions = stencil_phase.D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "timesteps = 238800\n",
      "Re = 10\n",
      "Pe = 0.41876046901171776\n",
      "Cn = 0.1\n",
      "domain_width = 50\n",
      "fluid_depth = 25.0\n",
      "amplitude = 0.5\n",
      "relaxation_rate_heavy = 1.99\n",
      "mobility = 0.02\n",
      "interface_width = 5\n",
      "density_heavy = 1.0\n",
      "density_light = 0.001\n",
      "density_ratio = 1000\n",
      "kinematic_viscosity_heavy = 0.0008375209380234356\n",
      "kinematic_viscosity_light = 0.0008375209380234356\n",
      "kinematic_viscosity_ratio = 1\n",
      "dynamic_viscosity_heavy = 0.0008375209380234356\n",
      "dynamic_viscosity_light = 8.375209380234357e-07\n",
      "dynamic_viscosity_ratio = 1000.0\n",
      "wavelength = 50\n",
      "wavenumber = 0.12566370614359174\n",
      "wave_frequency = 0.00033500837520937423\n",
      "surface_tension = 0\n",
      "gravitational_acceleration = -8.964447065459421e-07\n",
      "data_extract_frequency = 298\n",
      "vtk_output_frequency = 2985\n"
     ]
    }
   ],
   "source": [
    "# user defined input\n",
    "\n",
    "Re = 10 # Reynolds number\n",
    "domain_width = 50\n",
    "fluid_depth = 0.5 * domain_width\n",
    "amplitude = 0.01 * domain_width\n",
    "relaxation_rate_heavy = 1.99\n",
    "mobility = 0.02 # phase field mobility\n",
    "interface_width = 5 # phase field interface width\n",
    "density_heavy = 1.0 # density of heavy phase\n",
    "density_ratio = 1000\n",
    "density_light = density_heavy / density_ratio # density of light phase\n",
    "kinematic_viscosity_ratio = 1\n",
    "\n",
    "kinematic_viscosity_heavy = 1 / 3 * (1 / relaxation_rate_heavy - 0.5)\n",
    "kinematic_viscosity_light = kinematic_viscosity_heavy / kinematic_viscosity_ratio\n",
    "wavelength = domain_width\n",
    "wavenumber = 2.0 * np.pi / domain_width\n",
    "wave_frequency = Re * kinematic_viscosity_heavy / domain_width / amplitude # angular wave frequency\n",
    "surface_tension = 0\n",
    "gravitational_acceleration = - wave_frequency**2 / wavenumber / np.tanh(wavenumber * fluid_depth)\n",
    "Pe = domain_width * amplitude * wave_frequency / mobility # Peclet number\n",
    "Cn = interface_width / domain_width # Cahn number\n",
    "dynamic_viscosity_heavy = kinematic_viscosity_heavy * density_heavy\n",
    "relaxation_time_heavy = 3.0 * kinematic_viscosity_heavy\n",
    "kinematic_viscosity_light = kinematic_viscosity_heavy / kinematic_viscosity_ratio\n",
    "dynamic_viscosity_light = kinematic_viscosity_light * density_light\n",
    "relaxation_time_light = 3.0 * kinematic_viscosity_light\n",
    "\n",
    "timesteps = int(80 / wave_frequency)\n",
    "\n",
    "data_extract_frequency = int(0.1 / wave_frequency)\n",
    "vtk_output_frequency = int(1 / wave_frequency)\n",
    "vtk_output_path = \"vtk_out/gravity-wave\"\n",
    "vtk_base_directory = vtk_output_path.split(\"/\")[0] # create directory for vtk-output if it does not yet exist\n",
    "if not os.path.exists(vtk_base_directory):\n",
    "    os.mkdir(os.getcwd() + \"/\" + vtk_base_directory)\n",
    "\n",
    "domain_size = (domain_width, domain_width)\n",
    "filename = \"pf-re-\" + str(Re) + \"-resolution-\" + str(domain_width) + \".txt\"\n",
    "\n",
    "print(\"timesteps =\", timesteps)\n",
    "print(\"Re =\", Re)\n",
    "print(\"Pe =\", Pe)\n",
    "print(\"Cn =\", Cn)\n",
    "print(\"domain_width =\", domain_width)\n",
    "print(\"fluid_depth =\", fluid_depth)\n",
    "print(\"amplitude =\", amplitude)\n",
    "print(\"relaxation_rate_heavy =\", relaxation_rate_heavy)\n",
    "print(\"mobility =\", mobility)\n",
    "print(\"interface_width =\", interface_width)\n",
    "print(\"density_heavy =\", density_heavy)\n",
    "print(\"density_light =\", density_light)\n",
    "print(\"density_ratio =\", density_ratio)\n",
    "print(\"kinematic_viscosity_heavy =\", kinematic_viscosity_heavy)\n",
    "print(\"kinematic_viscosity_light =\", kinematic_viscosity_light)\n",
    "print(\"kinematic_viscosity_ratio =\", kinematic_viscosity_ratio)\n",
    "print(\"dynamic_viscosity_heavy =\", dynamic_viscosity_heavy)\n",
    "print(\"dynamic_viscosity_light =\", dynamic_viscosity_light)\n",
    "print(\"dynamic_viscosity_ratio =\", dynamic_viscosity_heavy/dynamic_viscosity_light)\n",
    "print(\"wavelength =\", wavelength)\n",
    "print(\"wavenumber =\", wavenumber)\n",
    "print(\"wave_frequency =\", wave_frequency)\n",
    "print(\"surface_tension =\", surface_tension)\n",
    "print(\"gravitational_acceleration =\", gravitational_acceleration)\n",
    "print(\"data_extract_frequency =\", data_extract_frequency)\n",
    "print(\"vtk_output_frequency =\", vtk_output_frequency)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "parameters = AllenCahnParameters(density_heavy=density_heavy, density_light=density_light,\n",
    "                                 dynamic_viscosity_heavy=dynamic_viscosity_heavy,\n",
    "                                 dynamic_viscosity_light=dynamic_viscosity_light,\n",
    "                                 gravitational_acceleration=gravitational_acceleration,\n",
    "                                 surface_tension=surface_tension, mobility=mobility,\n",
    "                                 interface_thickness=interface_width)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <table style=\"border:none; width: 100%\">\n",
       "            <tr style=\"border:none\">\n",
       "                <th style=\"border:none\" >Name</th>\n",
       "                <th style=\"border:none\" >SymPy Symbol </th>\n",
       "                <th style=\"border:none\" >Value</th>\n",
       "            </tr>\n",
       "            <tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">Density heavy phase</td>\n",
       "                            <td style=\"border:none\">$\\rho_{H}$</td>\n",
       "                            <td style=\"border:none\">$1.0$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">Density light phase</td>\n",
       "                            <td style=\"border:none\">$\\rho_{L}$</td>\n",
       "                            <td style=\"border:none\">$0.001$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">Relaxation time heavy phase</td>\n",
       "                            <td style=\"border:none\">$\\tau_{H}$</td>\n",
       "                            <td style=\"border:none\">$0.00251256281407031$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">Relaxation time light phase</td>\n",
       "                            <td style=\"border:none\">$\\tau_{L}$</td>\n",
       "                            <td style=\"border:none\">$0.00251256281407031$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">Relaxation rate Allen Cahn LB</td>\n",
       "                            <td style=\"border:none\">$\\omega_{\\phi}$</td>\n",
       "                            <td style=\"border:none\">$1.78571428571429$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">Gravitational acceleration</td>\n",
       "                            <td style=\"border:none\">$F_{g}$</td>\n",
       "                            <td style=\"border:none\">$-8.96444706545942 \\cdot 10^{-7}$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">Interface thickness</td>\n",
       "                            <td style=\"border:none\">$W$</td>\n",
       "                            <td style=\"border:none\">$5$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">Mobility</td>\n",
       "                            <td style=\"border:none\">$M_{m}$</td>\n",
       "                            <td style=\"border:none\">$0.02$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">Surface tension</td>\n",
       "                            <td style=\"border:none\">$\\sigma$</td>\n",
       "                            <td style=\"border:none\">$0$</td>\n",
       "                         </tr>\n",
       "\n",
       "        </table>\n",
       "        "
      ],
      "text/plain": [
       "<lbmpy.phasefield_allen_cahn.parameter_calculation.AllenCahnParameters at 0x12d1bd550>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fields\n",
    "\n",
    "As a next step all fields which are needed get defined. To do so we create a `datahandling` object. More details about it can be found in the third tutorial of the [pystencils framework]( http://pycodegen.pages.walberla.net/pystencils/). Basically it holds all fields and manages the kernel runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a datahandling object\n",
    "dh = ps.create_data_handling((domain_size), periodicity=(True, False), parallel=False, default_target=target)\n",
    "\n",
    "# fields \n",
    "g = dh.add_array(\"g\", values_per_cell=len(stencil_hydro))\n",
    "dh.fill(\"g\", 0.0, ghost_layers=True)\n",
    "h = dh.add_array(\"h\",values_per_cell=len(stencil_phase))\n",
    "dh.fill(\"h\", 0.0, ghost_layers=True)\n",
    "\n",
    "g_tmp = dh.add_array(\"g_tmp\", values_per_cell=len(stencil_hydro))\n",
    "dh.fill(\"g_tmp\", 0.0, ghost_layers=True)\n",
    "h_tmp = dh.add_array(\"h_tmp\",values_per_cell=len(stencil_phase))\n",
    "dh.fill(\"h_tmp\", 0.0, ghost_layers=True)\n",
    "\n",
    "u = dh.add_array(\"u\", values_per_cell=dh.dim)\n",
    "dh.fill(\"u\", 0.0, ghost_layers=True)\n",
    "\n",
    "rho = dh.add_array(\"rho\", values_per_cell=1)\n",
    "dh.fill(\"rho\", 1.0, ghost_layers=True)\n",
    "\n",
    "C = dh.add_array(\"C\")\n",
    "dh.fill(\"C\", 0.0, ghost_layers=True)\n",
    "C_tmp = dh.add_array(\"C_tmp\")\n",
    "dh.fill(\"C_tmp\", 0.0, ghost_layers=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set the frequency of the file outputs\n",
    "vtk_writer = dh.create_vtk_writer(vtk_output_path, [\"C\", \"u\", \"rho\"], ghost_layers=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho_L = parameters.symbolic_density_light\n",
    "rho_H = parameters.symbolic_density_heavy\n",
    "density = rho_L + C.center * (rho_H - rho_L)\n",
    "\n",
    "body_force = [0, 0, 0]\n",
    "body_force[1] = parameters.symbolic_gravitational_acceleration * density"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definition of the lattice Boltzmann methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_c = parameters.symbolic_omega_phi\n",
    "\n",
    "config_phase = LBMConfig(stencil=stencil_phase, method=Method.MRT, compressible=True,\n",
    "                         delta_equilibrium=False,\n",
    "                         force=sp.symbols(\"F_:2\"), velocity_input=u,\n",
    "                         weighted=True, relaxation_rates=[0, w_c, w_c, 1, 1, 1, 1, 1, 1],\n",
    "                         output={'density': C_tmp}, kernel_type='stream_pull_collide')\n",
    "\n",
    "method_phase = create_lb_method(lbm_config=config_phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "omega = parameters.omega(C)\n",
    "config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.MRT, compressible=False,\n",
    "                         weighted=True, relaxation_rates=[omega, 1, 1, 1],\n",
    "                         force=sp.symbols(\"F_:2\"),\n",
    "                         output={'velocity': u}, kernel_type='collide_stream_push')\n",
    "\n",
    "method_hydro = create_lb_method(lbm_config=config_hydro)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <table style=\"border:none; width: 100%\">\n",
       "            <tr>\n",
       "                <th colspan=\"3\" style=\"text-align: left\">\n",
       "                    Moment-Based Method\n",
       "                </th>\n",
       "                <td>Stencil: D2Q9</td>\n",
       "                <td>Zero-Centered Storage: &#10003;</td>\n",
       "                <td>Force Model: Guo</td>\n",
       "            </tr>\n",
       "        </table>\n",
       "        \n",
       "        <table style=\"border:none; width: 100%\">\n",
       "            <tr>\n",
       "                <th colspan=\"3\" style=\"text-align: left\">\n",
       "                    Continuous Hydrodynamic Maxwellian Equilibrium\n",
       "                </th>\n",
       "                <td rowspan=\"2\" style=\"width: 50%; text-align: center\">\n",
       "                    $f (\\rho, \\left( u_{0}, \\  u_{1}\\right), \\left( v_{0}, \\  v_{1}\\right)) \n",
       "                        = \\frac{3 \\delta_{\\rho} e^{- \\frac{3 v_{0}^{2}}{2} - \\frac{3 v_{1}^{2}}{2}}}{2 \\pi} - \\frac{3 e^{- \\frac{3 v_{0}^{2}}{2} - \\frac{3 v_{1}^{2}}{2}}}{2 \\pi} + \\frac{3 e^{- \\frac{3 \\left(- u_{0} + v_{0}\\right)^{2}}{2} - \\frac{3 \\left(- u_{1} + v_{1}\\right)^{2}}{2}}}{2 \\pi}$\n",
       "                </td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                <td>Compressible: &#10007;</td>\n",
       "                <td>Deviation Only: &#10003;</td>\n",
       "                <td>Order: 2</td>\n",
       "            </tr>\n",
       "        </table>\n",
       "        \n",
       "        <table style=\"border:none; width: 100%\">\n",
       "            <tr> <th colspan=\"3\" style=\"text-align: left\"> Relaxation Info </th> </tr>\n",
       "            <tr>\n",
       "                <th>Moment</th>\n",
       "                <th>Eq. Value </th>\n",
       "                <th>Relaxation Rate</th>\n",
       "            </tr>\n",
       "        <tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">$1$</td>\n",
       "                            <td style=\"border:none\">$\\delta_{\\rho}$</td>\n",
       "                            <td style=\"border:none\">$0.0$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">$x$</td>\n",
       "                            <td style=\"border:none\">$u_{0}$</td>\n",
       "                            <td style=\"border:none\">$0.0$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">$y$</td>\n",
       "                            <td style=\"border:none\">$u_{1}$</td>\n",
       "                            <td style=\"border:none\">$0.0$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">$x^{2} - y^{2}$</td>\n",
       "                            <td style=\"border:none\">$u_{0}^{2} - u_{1}^{2}$</td>\n",
       "                            <td style=\"border:none\">$\\frac{2}{2 {C}_{(0,0)} \\left(\\tau_{H} - \\tau_{L}\\right) + 2 \\tau_{L} + 1}$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">$x y$</td>\n",
       "                            <td style=\"border:none\">$u_{0} u_{1}$</td>\n",
       "                            <td style=\"border:none\">$\\frac{2}{2 {C}_{(0,0)} \\left(\\tau_{H} - \\tau_{L}\\right) + 2 \\tau_{L} + 1}$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">$3 x^{2} + 3 y^{2} - 2$</td>\n",
       "                            <td style=\"border:none\">$3 u_{0}^{2} + 3 u_{1}^{2}$</td>\n",
       "                            <td style=\"border:none\">$1$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">$3 x^{2} y - y$</td>\n",
       "                            <td style=\"border:none\">$0$</td>\n",
       "                            <td style=\"border:none\">$1$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">$3 x y^{2} - x$</td>\n",
       "                            <td style=\"border:none\">$0$</td>\n",
       "                            <td style=\"border:none\">$1$</td>\n",
       "                         </tr>\n",
       "<tr style=\"border:none\">\n",
       "                            <td style=\"border:none\">$9 x^{2} y^{2} - 3 x^{2} - 3 y^{2} + 1$</td>\n",
       "                            <td style=\"border:none\">$0$</td>\n",
       "                            <td style=\"border:none\">$1$</td>\n",
       "                         </tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x12d4a1e50>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "method_hydro"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[['spatialInner0'], ['spatialInner1']]\n",
      "[['spatialInner0'], ['spatialInner1']]\n"
     ]
    }
   ],
   "source": [
    "h_updates = initializer_kernel_phase_field_lb(method_phase, C, u, h, parameters)\n",
    "g_updates = initializer_kernel_hydro_lb(method_hydro, 1, u, g)\n",
    "\n",
    "h_init = ps.create_kernel(h_updates, target=dh.default_target, cpu_openmp=True).compile()\n",
    "g_init = ps.create_kernel(g_updates, target=dh.default_target, cpu_openmp=True).compile()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize the phase field\n",
    "def initialize_phasefield():\n",
    "    Nx = domain_size[0]\n",
    "    Ny = domain_size[1]\n",
    "    \n",
    "    for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):\n",
    "        # get x as cell center coordinate, i.e., including shift with 0.5\n",
    "        x = np.zeros_like(block.midpoint_arrays[0])\n",
    "        x[:, :] = block.midpoint_arrays[0]\n",
    "        \n",
    "        # get y as cell center coordinate, i.e., including shift with 0.5\n",
    "        y = np.zeros_like(block.midpoint_arrays[1])\n",
    "        y[:, :] = block.midpoint_arrays[1] \n",
    "        \n",
    "        tmp = fluid_depth + amplitude * np.cos(x / (Nx - 1) * np.pi * 2 + np.pi)\n",
    "        \n",
    "        # initialize diffuse interface with tanh profile\n",
    "        init_values = 0.5 - 0.5 * np.tanh((y - tmp) / (interface_width / 2))\n",
    "        \n",
    "        block[\"C\"][:, :] = init_values\n",
    "        block[\"C_tmp\"][:, :] = init_values\n",
    "        \n",
    "    if gpu:\n",
    "        dh.all_to_gpu()            \n",
    "    \n",
    "    dh.run_kernel(h_init, **parameters.symbolic_to_numeric_map)\n",
    "    dh.run_kernel(g_init)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6gAAAFlCAYAAADmqMVrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABMiElEQVR4nO3deXxU1cHG8edkMkkmEBh2SNgRosgWDJtoXdu01ipS9w1Q3FutbdO3dH27WKypWm3dEcRda2Ns3VIXtK5AIEAQCPuWAAnLhEAmyWTmvH8k8rIkECDkzvL7fj58TM7MTR7kZuY+ufeeY6y1AgAAAADAaXFOBwAAAAAAQKKgAgAAAADCBAUVAAAAABAWKKgAAAAAgLBAQQUAAAAAhAUKKgAAAAAgLMQ7HaAxnTt3tn379nU6BgAAAACghS1YsGC7tbZLY4+FZUHt27evCgoKnI4BAAAAAGhhxpgNTT3GJb4AAAAAgLBAQQUAAAAAhAUKKgAAAAAgLFBQAQAAAABhgYIKAAAAAAgLFFQAAAAAQFigoAIAAAAAwgIFFQAAAAAQFiioAAAAAICwEO90AAAAolVeYYly8otV6vMr1etRdla6JmSktdr2AABEGgoqAAAnQF5hiablFskfCEqSSnx+TcstkqRmlczj3R4AgEhEQQUA4AS4L3/FvnL5NX8gqD+8uUztk91H3P4Pby5rdPuc/GIKKgAgalFQAQBoRHMur7XWqnxPjdZvr9K67Xu0dvterd++V+u271Wpr7rRr7tjb62mzJp/zLlKfH7d9vwC9e3cRv06t1H/zm3Ut3MbdWqTIGPMMf09AAAIFxRUAAAO0tjltf/zzyWav2GHOrVJ0rrte7Vu+x6t316lPTV1+7Zzu4z6dGqjvp3aqNRXfcBjX+vSNlFPXn/aETPc/OwCle+pOWQ8KT5Oxdsq9d6ybaoL2X3jKUnx+8pqv4Y/G3bs1aMfrVF1ILTv78FlwgCAcGastUd+VivLzMy0BQUFTscAAMSoMX96X9t2H1oOJSnOSGkdPOrXuW19IeyUrH5d6j9O9Xrkiqs/i3lwyZUkj9ul6ROHHtM9qAdvXxcMafMuv9bt2Kt15Xu1fkf9mdu15XtVWuHX4d7e07weffbzc5v5fwMAgJZljFlgrc1s7DHOoAIAYl4gGFLB+l2aU1ymOSvKmiynRtLyP3xbifGuI37Nr0vosV5ee6Tt411x6ttwxvSc9AO3rQ4EtXFnlb714H8b/dolPr/ue3eFzjm5qzJ6eRXvYtU5AEB44AwqACAmlVVW6+Pics0pLtMnK7ersqZObpfR6H4dtbSkQhX+Qy/PjbQzj+Pv/VAlPv8h4wmuOAWtVTBk1d7j1jcGddG5J3fRWYO6qmObBAeSAgBiCWdQAQAx5+DJgX76zUHq17WtPlxRpo+Ky7Rkc4UkqWtKoi4Y2kPnnNxV40/qpJQkd5OX12ZnpTf17cJSdlZ6k5cJn3NyV326ars+XFGmj1eW6d+LS2WMNKKXV+ekd9W5J3fVqant9MaiUiZZAgC0Gs6gAgCiTmMF82vGSBkNJeychhIWzbPfNufvEQpZFZVU7CvvixvKe0qiS1WBkIL7TcZ0NPfRAgDQmMOdQaWgAgCizqg/vt/oDLgdkt364CdncxnrEZRX1ujjleX6dV6R/A0zAO8v0i51BgCEl8MVVGZFAABEhdq6kN5YVKJLHv2s0XIqSb6qAOW0GbqkJOrS03ruW57mYCU+v94p2qK6YOOPAwBwrLgHFQAQ0bbvqdFLczfquS83qKyyRv06t1F7T3yjkxylej0OJIxcqV5Po5MsueKMbnthoVLbJ+m6cX111ehe8iZT/AEAx4+CCgCISEtLKjTrs/X69+JS1QZD+sagLvrz9/vqrEFd9K/FpVExyZHTmppk6Z4JQ9QmKV7PfLZef353hR76YKUuyUjTpNP76uTu7RxMDACIdBRUAEDEqAuGlP/VNj3z+TrNX79LyQkuXTGqlyad3lcndW2773nHuwYp6h3p/2PWqd21Yutuzf58vXIXluileZt0+oBOmnx6X513Sje54g6dfAoAgMNhkiQAQNg5eObZ284eoN3VAT33xQZtqahWr44eTRrXV5dl9lJ7j9vpuJC0a2+tXp6/Sc99sV6lFdXq2aH+3+jyUb00Z0UZvywAAOzDLL4AgIhxuCVixp/USZNP76dzT+7K2bkwVRcM6T/LtumZz9Zr3vqdcscZhSSWqgEA7MMsvgCAiJGTX9xoOe2akqgXpo7VNwdz6Wg4i3fF6YKhPfTqreP05g/PkDs+7oByKkn+QFA5+cUOJQQAhDMKKgAgbJRX1jQ6a+zXjyGyDElrL3/tob9skKTSJv6dAQCxjYIKAHBcdSCoR+as1tk5c5p8DkvERKbD/bvN+GStautYSxUA8P8oqAAAx1hr9a/FpTrv/o+Vk1+scQM66xcXnCyP23XA81giJnJlZ6Uf8u+ZGB+nQd3a6o9vLde3HvxY7y7dqnCcEwMA0PpYZgYA4IgFG3bpj28tU+FGnwb3aKecS4fp9JM6S5K6piQx62uUONxSNR8Vl+met5br1ucXaEy/jvr1hYM1JK29w4kBAE5iFl8AQKvatLNKf353hd5cskVdUhKVnZWu74/sycRHMaouGNJL8zfpwfdWaldVrSZm9FR2Vrq6t09yOhoA4ARhmRkAgOMqqwN69KM1evrTdYoz0s1n9tctZw1Qm0Qu5oG0uzqgR+as1qxP18sVZ3TzN/rrlrP6KzmB/QMAog0FFQDQavIKSw64nPMn3xwof11ID/xnpXbsrdXEjDT9NCudSY/QqE07q3Tvuyv01pIt6tYuUdlZJ2tiRpr+tbiUy74BIEpQUAEArSKvsETTcosOWMfUSLKSRvXtoF9fOFjDenqdiocIUrB+p/7w5jIt3lyhnt4kle2pPWDGX4/bpekTh1JSASACHa6gMosvAKDF5OQXH1BOpfpy2jHZrVdvGUc5RbNl9u2o128fr4euHKEtFdWHLEfjDwSVk1/sUDoAwIlCQQUAtJhSn7/R8V1VARnDJEg4OnFxRhePSFOoiYu9mtrfAACRi4IKAGgRZburlRDf+NsK95vieDS1/zDTLwBEnyMWVGNML2PMHGPMcmPMV8aYuxrG/9cYU2KMWdTw54LDfA2XMabQGPNmS4YHAISHN5eU6lt//a/qgiG5XQeeKfW4XcrOSncoGaJBdla6PG7XIeOV1QG9v2ybA4kAACdKc86g1kn6ibX2FEljJd1hjBnc8NiD1toRDX/ePszXuEvS8uPMCgAIM76qWt35UqF+8GKh+nRMVv7dZynn0uFK83pkJKV5PUxkg+M2ISNN0ycOPWC/ys5KV88OyZr6bIH+57UlqqwOOB0TANACjnoWX2PMG5L+Lmm8pD3W2r8c4fk9Jc2WdI+kH1trLzzS92AWXwAIfx+vLNfPXlusHXtqded5A3X72QMU7+LOEbSemrqgHnp/lR7/eI1SvR795bLhGtu/k9OxAABH0GKz+Bpj+krKkDS3YegHxpglxpiZxpgOTWz2V0k/kxRq4nEAQASpqq3Tr/KKNGnmPLVLcuv128frzvMGUk7R6hLjXfrZt0/WP24dJ1ec0VVPfak/vrlM1QfNJA0AiBzNPpowxrSV9E9JP7LW7pb0mKQBkkZI2iLp/ka2uVBSmbV2QTO+/s3GmAJjTEF5eXlzYwEAWtGCDTv1nYc+0QtzN+qmM/vp3z88Q0N7tnc6FmLcaX066u07z9Q1Y3prxqfr9L2/faqizRVOxwIAHINmXeJrjHFLelNSvrX2gUYe7yvpTWvtkIPGp0u6TvX3sSZJaicp11p77eG+H5f4AkB4qakL6q/vr9ITH69Rj/Ye3X85l1IiPO1/6fkPzx2o288ZIDdn9wEgrBzuEt8jFlRTv3DdbEk7rbU/2m+8h7V2S8PHd0saY6298jBf52xJP+UeVACILMu37NbdryzSiq2VuiKzl3514SlKSXI7HQtoUkVVQL/511K9sahUw3u21/2Xj9BJXds6HQsA0OBwBTW+GduPV/1Z0CJjzKKGsV9IusoYM0KSlbRe0i0N3yxV0gxrbZPLzgAAwldeYYly8otV6vMrJSlee2rq1LFNomZcn6nzB3dzOh5wRO2T3Xroygx9c3A3/Spvqb778Cf67tDu+nLdTm3xVSu1YRZgZpcGgPBz1LP4tgbOoAKAM/IKSzQtt0j+/SaZiTPS7y8eomvH9nEwGXBsynZXa9LMeVq+tfKAcY/bxRJIAOCQFpvFFwAQ3XLyiw8op5IUstJjH61xKBFwfLq2S9LuRtZI9QeCyskvdiARAOBwKKgAAElSKGRV4vM3+lhpE+NAJCj1VTcxzn4NAOGGggoA0O7qgG56tulbK1K9nlZMA7SspvbfJHecqmrrWjkNAOBwKKgAEONWbavUhL9/po9Xluv7I9PkcR/41uBxu5Sdle5QOuD4ZWely+N2HTAWH2fkD4Q08dHPtWHHXoeSAQAORkEFgBiW/9VWTXjkM+2uDujFm8bq/stHaPrEYUrzemQkpXk9TCSDiDchI03TJw49YL/+y2XD9ewNo7WloloXNfyCBgDgPGbxBYAYFApZPfj+Sv3tw9Ua3rO9Hr/uNPVoz2W8iD0bd1Tp5ucKtHJbpbKzTtatZ/VX/RLwAIAThVl8AQD7VPgDunH2fP3tw9W67LSeeuWWcZRTxKzenZKVe/vpumBoD/353RX6wYuF2lvDfakA4JR4pwMAAFrPqm2Vuvm5Bdq0s0p/mDBE147pzdkixLzkhHj97aoMDevZXve+s0Kry/boietOU9/ObZyOBgAxhzOoABAj3l26RRMe+UyV1XV66eaxum5sH8op0MAYo5u/MUCzbxitbZXVuujvn+qj4jKnYwFAzKGgAkCUC4as/pJfrFufX6iB3VL05g/P0Ki+HZ2OBYSlMwd20b9/cIbSOiRryjPz9cic1QrH+ToAIFpRUAEgin19v+nf56zWlaN66ZVbxqp7+ySnYwFhrVfHZOXedrq+NyxVOfnFuv2FhdyXCgCthHtQASCK5BWWKCe/WKU+v7qkJCoUCqmiuk73XDJEV4/mflOguTwJLj105QgN69lef3p7udaU79Flp/XSM5+vV6nPr1SvR9lZ6SzBBAAtjIIKAFEir7BE03KL5A8EJUlllTWSpDvPO0nXjOnjZDQgIhljNPXM/jqlRzvdNHu+7nl7+b7HSnx+TcstkiRKKgC0IC7xBYAokZNfvK+c7u+fC0ocSANEj/EndVaKx33IuD8QVE5+sQOJACB6UVABIEqU+vxHNQ6g+cp21zQ6zs8XALQsCioARIHqQFCJ8Y2/pKd6Pa2cBog+Tf0c9WDSMQBoURRUAIhwvqpaXTtjrqrrQnK7DpwEyeN2KTsr3aFkQPTIzkqXx+06ZLxtYrz2MMMvALQYCioARLBSn1+XPf6Flmyu0CNXj1TOpcOV5vXISErzejR94lAmcAFawISMNE2fOPSAn6+rRvfSmu17deWTX6i8svFLgAEAR8eE4+LTmZmZtqCgwOkYABDWirdWatLMedpbU6cnr8/UuAGdnI4ExJw5xWW6/fmF6pKSqNk3jFa/zm2cjgQAYc8Ys8Bam9nYY5xBBYAINHftDl36+Oeysnr11nGUU8Ah56R31Us3j9Wemjp9/7HPtXiTz+lIABDRKKgAEGHeKdqi62bOU9eURP3zttN1So92TkcCYtqIXl69dus4JSe4dOWTX+qj4jKnIwFAxKKgAkAEefaL9br9xYUaktpOr916unp2SHY6EgBJ/bu0Ve7tp6t/lzaaOrtAry3Y7HQkAIhIFFQAiADWWuXkr9Bv3vhK553cTS9MHasObRKcjgVgP11TkvTyzWM1pn9H/fQfi/XoR6sVjnN9AEA4o6ACQJgLBEP62WtL9MicNbpqdC89fu1IeRIOXe4CgPNSktyaNXm0LhqeqvveLdbv/r1MwRAlFQCaK97pAACAplXV1umOFxZqTnG57jpvoH50/kAZY468IQDHJMTH6a9XjFDXlETN+HSdyitrdP/lw5XUyDqqAIADUVABIEzt2FOjG2YXqGizT/dcMkTXjOnjdCQAzRQXZ/SrCwerW7sk3fP2cm3fU6Mnr89Ue4/b6WgAENYoqAAQRvIKS5STX6xSn19xcUayVo9fe5q+dWp3p6MBOAY3faO/urZL1E//sVjffvBjWRlt212tVK9H2VnpmpCR5nREAAgrFFQACBN5hSWallskfyAoSQqGrBLi41RVG3Q4GYDjcfGINK3YUqnHPl6zb6zE59e03CJJoqQCwH6YJAkAwkROfvG+cvq12rqQcvKLHUoEoKX8a3HpIWP+QJCfbwA4CAUVAMJEqc9/VOMAIgc/3wDQPBRUAAgDX6zZ0eRjqV5PKyYBcCI09XPcuW1iKycBgPBGQQUAh326arumPDNPXVMSlRR/4Muyx+1Sdla6Q8kAtJTsrHR5Dlpmxkja7a/Vgg07nQkFAGGIggoADvp4ZblunD1ffTu10Vt3nal7vz9MaV6PjKQ0r0fTJw5lAhUgCkzISNP0iUMP+Pn+7UWDldohWdc/PU/z1lFSAUCSjLXW6QyHyMzMtAUFBU7HAIATas6KMt3y3AKd1LWtnp86Rh3bJDgdCUAr27a7Wlc/9aVKfdWaOXmUxg3o5HQkADjhjDELrLWZjT3GGVQAcMB7y7bp5ucKlN49RS/eRDkFYlW3dkl6+eZx6tnBoynPzNOnq7Y7HQkAHEVBBYBW9u7Srbrt+QUanNpez08dI28y5RSIZV1SEvXyzWPVt1Mb3Th7vj5eWe50JABwDAUVAFrRW0u26I4XF2pYz/Z67sbRau9xOx0JQBjo1DZRL940VgO6tNVNswv04YptTkcCAEdQUAGglbyxqER3vlyokb29evbGMWqXRDkF8P86tknQizeNUXr3FN3y3AK9t4ySCiD2UFABoBXkLtysu19ZpFF9O+iZKaPVNjHe6UgAwpA3OUHPTx2jwantddvzC/Tu0i1ORwKAVkVBBYAT7NWCTfrJPxZr3IBOmjV5tNpQTgEcRnuPW8/dOFrDerbXHS8W6s0lpU5HAoBWQ0EFgBPopXkb9bPXluiMkzrr6Umj5ElwOR0JQARol+TWszeO0cjeXt35UqHeWFTidCQAaBUUVAA4QZ77coOm5RbpnPQueur6TCW5KacAmq9tYryemTJao/t11N2vLFLuws1ORwKAE47rzACgheQVlignv1ilPr/aeeJV4a/T+ad00yPXZCgxnnIK4Oi1SYzXrMmjNfXZ+frJPxZr3rqd+mTVdpX6/Er1epSdla4JGWlOxwSAFsMZVABoAXmFJZqWW6QSn19WUoW/TnFG+vap3SinAI6LJ8GlpyeN0qBuKXp5/qZ9rzMlPr+m5RYpr5DLfwFEDwoqALSAnPxi+QPBA8ZCVnrw/VUOJQIQTZLcLlX6A4eM+wNB5eQXO5AIAE4MCioAtIBSn/+oxgHgaG2pqG50nNcZANGEggoALcCb7G50PNXraeUkAKJVU68nvM4AiCYUVAA4Tu8UbdGuqoDizIHjHrdL2VnpzoQCEHWys9LlaWQ28O8O6+5AGgA4MSioAHAc5qwo050vFyqzTwfdO3Go0rweGUlpXo+mTxzK7JoAWsyEjDRN3+91pkf7JKV5kzT78w36Ys0Op+MBQIsw1lqnMxwiMzPTFhQUOB0DAA7r8zXbNWXWfA3s1lYv3jRW7ZIav8wXAE6UnXtrdcUTX6jE59fzU8doZO8OTkcCgCMyxiyw1mY29hhnUAHgGCzYsEtTZxeoT6dkPXvDGMopAEd0bJOg56eOUZeURE2eOU9flVY4HQkAjssRC6oxppcxZo4xZrkx5itjzF0N4/9rjCkxxixq+HNBc7cFgEi2tKRCk2fNU9eURD1/4xh1bJPgdCQAMaxbuyS9MHWM2ibG6/qn52l12R6nIwHAMWvOGdQ6ST+x1p4iaaykO4wxgxsee9BaO6Lhz9tHuS0ARJzVZZW6fuY8pSTG6/mpY9S1XZLTkQBAPTsk6/mpY2SM0TUzvtTGHVVORwKAY3LEgmqt3WKtXdjwcaWk5ZKaNevH8WwLAOFmw469uvqpuYozRi/cNFY9OyQ7HQkA9unfpa2enzpaNXUhXfP0l9pSwfqoACLPUd2DaozpKylD0tyGoR8YY5YYY2YaYw57V34j2x78+M3GmAJjTEF5efnRxAKAE67U59fVT81VbTCkF6aOUb/ObZyOBACHOLl7Oz17w2jt2hvQNTPmavueGqcjAcBRaXZBNca0lfRPST+y1u6W9JikAZJGSNoi6f6j2PYQ1tonrbWZ1trMLl26NP9vAAAnWHllja6dMVe7/QE9d8MYpXdPcToSADRpWE+vZk4epVKfX9fOmCtfVa3TkQCg2ZpVUI0xbtUXzBestbmSZK3dZq0NWmtDkp6SNLq52wJApPBV1eq6p+dqS0W1Zk4ZpaE92zsdCQCOaHS/jnryukytLd+rSbPma09NndORAKBZmjOLr5H0tKTl1toH9hvvsd/TLpG0tLnbAkAkqKwOaNKs+VpbvldPXZ+pUX07Oh0JAJrtG4O66O9XZ2hpSYVueGa+/LVBpyMBwBE15wzqeEnXSTr3oCVl7jPGFBljlkg6R9LdkmSMSTXGvH2EbQEgrPlrg7pxdoGWllTokWtG6oyBnZ2OBABH7VundtcDlw/X/PU7devzC1RTR0kFEN6MtdbpDIfIzMy0BQUFTscAEGPyCkuUk1+sUp9fCfFxqqkL6eGrMnTR8FSnowHAcXll/kb9zz+LNCytnbbvrdUWX7VSvR5lZ6VrQgYLLABoXcaYBdbazMYei2/tMAAQjvIKSzQtt0j+QP3ZhZq6kNwuo1Ao/H6JBwBH64pRvfXlmh16fVHpvrESn1/TcoskiZIKIGwc1TIzABCtcvKL95XTrwWCVjn5xQ4lAoCWNW/9rkPG/IEgr3MAwgoFFQBUv87p0YwDQKThdQ5AJKCgAoCklKTG73hI9XpaOQkAnBhNvZ7xOgcgnFBQAcS8l+dt1O7qOrmMOWDc43YpOyvdoVQA0LKys9LlcbsOGT/3lK4OpAGAxlFQAcS0/K+26hevF+msQV1036VDleb1yEhK83o0feJQJg4BEDUmZKRp+sT/f51LbZ+k/p3b6OV5G/XJqnKn4wGAJJaZARDD5q7doetmztMpPdrpxalj1CaRic0BxJYKf0BXPPGFNu6s0ss3j9Wwnl6nIwGIAYdbZoYzqABi0vItuzX12QL17ODRrMmjKKcAYlJ7j1uzbxitjm0SNHnWfK0t3+N0JAAxjoIKIOZs2lmlSTPnqU1CvJ5tODADgFjVrV2Snr1htCTp+pnztG13tcOJAMQyCiqAmLJjT42unzlP1YGgZt8wWj07JDsdCQAc179LWz0zZZR27a3VpJnzVOEPOB0JQIyioAKIGXtq6jTlmfkq9fk1c/IopXdPcToSAISNYT29evy607SmfI9uml2g6kDQ6UgAYhAFFUBMqK0L6bbnF+ir0t165OqRyuzb0elIABB2zhzYRfdfPkLzN+zUnS8Vqi4YcjoSgBhDQQUQ9UIhq5/+Y7E+WbVd0ycO1fmDuzkdCQDC1kXDU/XbCwfrP8u26ddvLFU4rvgAIHoxbSWAqGat1e/fXKZ/LS7V/3z7ZF2e2cvpSAAQ9iaP76fte2r19zmr1blton7yrXSnIwGIERRUAFHt0Y/W6JnP1+uG8f1061n9nY4DABHjJ98apO17avS3D1erU5sETR7fz+lIAGIABRVA1Hpl/kbl5Bfr4hGp+tV3T5ExxulIABAxjDH644Qh2rG3Vr97c5k6tU3U94anOh0LQJTjHlQAUem9Zds0LbdI3xjURTmXDldcHOUUAI5WvCtOf7sqQ6P6dNSPX12kT1dtdzoSgChHQQUQdeav36kfvLhQQ3t69dg1I5UQz0sdAByrJLdLT03K1IAubXXLcwUq2lzhdCQAUcyE48xsmZmZtqCgwOkYACJIXmGJcvKLVerzS5I6pyQo/0dnqWObBIeTAUB02La7WhMf/Vy+qlq1TYpX2e4apXo9ys5K14SMNKfjAYggxpgF1trMxh7jtAKAiJdXWKJpuUUq8fllJVlJu/11+u/KcqejAUDU6NYuSZNO76O9tUFt210jK6nE59e03CLlFZY4HQ9AlKCgAoh4OfnF8geCB4zV1IWUk1/sUCIAiE6zP99wyJg/EOT1FkCLoaACiHhfX9bb3HEAwLHh9RbAiUZBBRDRQiGrRHfjL2WpXk8rpwGA6NbU62qqN6mVkwCIVhRUABHtnreXqzoQktt14DIyHrdL2VnpDqUCgOiUnZUuj9t1yPiQ1PYOpAEQjSioACLWjE/W6ulP12nK+L667/vDlOb1yEhK83o0feJQZpUEgBY2ISNN0ycO3e/1Nklj+nVQ/rJtenneRqfjAYgC8U4HAIBj8daSLbrn7eX6zpDu+tV3B8sVZ3TJyJ5OxwKAqDchI+2AXwAGgiFNnV2gX+YtVbd2STrn5K4OpgMQ6TiDCiDizFu3U3e/ukin9e6gB68YIVecOfJGAIATwu2K06PXjNQpPVJ0x4sLVbS5wulIACIYBRVARFldVqmbni1Qzw4ePXV9ppIauRcKANC62iTGa+bkUerYJkFTnpmvTTurnI4EIEJRUAFEjLLd1Zo0c77crjjNnjJaHdokOB0JANCga0qSnpkyWoFgSJNmzdOuvbVORwIQgSioACLCnpo6TXlmvnZV1WrW5FHq1THZ6UgAgIOc1LWtZkzK1OZdfk19tkDVgaDTkQBEGAoqgLAXCIZ0xwsLtWJrpR65eqSG9mQ5AwAIV6P6dtRfrxihhRt36UcvL1IwZJ2OBCCCUFABhDVrrX75epE+XlmueyYMYXZIAIgAFwztoV99d7De/Wqr/vDmMllLSQXQPCwzAyCsPfTBKr1asFl3njdQV47u7XQcAEAz3XhGP5X6/Hr603Xq2cGjqWf2dzoSgAhAQQUQtl6dv0l/fX+VLj2tp+4+f6DTcQAAR+mXF5yirRXV+uNby9W9fZIuHJbqdCQAYY6CCiAsfVRcpmmvF+nMgZ01feJQGcNapwAQaeLijO6/fLjKKqv141cWq0vbRI3p38npWADCGPegAgg7S0sqdPsLC5XeLUWPXXua3C5eqgAgUiW5XXrq+kz16ujRTc8WaNW2SqcjAQhjJhxvWs/MzLQFBQVOxwDQSvIKS5STX6xSn19d2yVqb02d2nsSlHv76erWLsnpeACAFrBpZ5UmPva5AnUhJbld2ra7Wqlej7Kz0jUhI83peABakTFmgbU2s7HHOC0BwFF5hSWallukEp9fVtK23TXaUxPUdWN7U04BIIr06pis68f1kc8f0Nbd1bKSSnx+TcstUl5hidPxAIQJCioAR+XkF8vfyELuz3250YE0AIAT6eV5mw4Z8weCyskvdiANgHBEQQXgqFKf/6jGAQCRi9d8AEdCQQXgqFSv56jGAQCRi9d8AEdCQQXgqHEDDl1uwON2KTsr3YE0AIATKTsrXR6365Dxs0/u7EAaAOGIggrAMR8s36bchZs1uEeKUr1JMpLSvB5NnziUGR0BIApNyEjT9IlDleb1yEhKbZ+kgV3b6NX5m/X5mu1OxwMQBlhmBoAjijZX6PInvtBJXdvqlVvGKjkh3ulIAAAHVPgDuuzxz7Wlolq5t52ugd1SnI4E4ARjmRkAYWXzrirdMHu+OrZJ0NOTMymnABDD2nvcmjl5lJLcLk2eNV9lldVORwLgIAoqgFZVURXQ5FnzVR0I6pkpo9Q1hbVOASDW9eyQrFmTR2lXVa1ueGa+9tbUOR0JgEMoqABaTU1dULc8X6ANO/bqietO4zIuAMA+Q9La65GrR2pZ6W798KVC1QVDTkcC4AAKKoBWYa3Vz/9ZpC/X7lTOpcN1+gBmbAQAHOick7vq9xcP0YcryvS///5K4ThXCoATixu/ALSKB95bqdcLS/TTbw1ihl4AQJOuHdtHm3ZV6YmP16pXh2TdctYApyMBaEUUVAAn3KvzN+lvH67WlaN66Y5zTnI6DgAgzP1P1skq2eXX9HdWKK2DRxcOS3U6EoBWcsRLfI0xvYwxc4wxy40xXxlj7moY/19jTIkxZlHDnwua2P7bxphiY8xqY8zPW/ovACC8/Xdluaa9XqRvDOqiP0wYImOM05EAAGEuLs7oL5cN16i+HfTjVxerYP1OpyMBaCXNuQe1TtJPrLWnSBor6Q5jzOCGxx601o5o+PP2wRsaY1ySHpH0HUmDJV2137YAotyy0t26/YWFGti1rR65OkNuF7e9AwCaJ8nt0pPXZaqn16OpzxZobfkepyMBaAVHPFq01m6x1i5s+LhS0nJJzb2BbLSk1dbatdbaWkkvS7r4WMMCiBxbKvy64Zn5apsYr1lTRiklye10JABAhOnQJkHPTBktlzGaPGu+duypcToSgBPsqE5nGGP6SsqQNLdh6AfGmCXGmJnGmA6NbJImadN+n29W88stgAhVWR3QlFnztaemTrOmjFKP9h6nIwEAIlTvTsmaMSlTZZXVunF2gfy1QacjATiBml1QjTFtJf1T0o+stbslPSZpgKQRkrZIur+xzRoZa3S+cGPMzcaYAmNMQXl5eXNjAQgzgWBIt7+wUKvL9uixa0fqlB7tnI4EAIhwGb076KErM7R4s08/eqVQwRDLzwDRqlmz+Bpj3Kovpy9Ya3MlyVq7bb/Hn5L0ZiObbpbUa7/Pe0oqbex7WGuflPSkJGVmZvKqA0SQvMIS5eQXq9TnlyfBparaoO67dJjOHNjF6WgAgCiRdWp3/fq7g/X7N5fphlnztLp8r0p9fqV6PcrOSmcJMyBKHLGgmvopN5+WtNxa+8B+4z2stVsaPr1E0tJGNp8vaaAxpp+kEklXSrr6uFMDCBt5hSWallskf6D+kquq2qDi44wSmBAJANDCbjijnz4qLtPHq7bvGyvx+TUtt0iSKKlAFGjOEeR4SddJOvegJWXuM8YUGWOWSDpH0t2SZIxJNca8LUnW2jpJP5CUr/rJlV611n51Iv4iAJyRk1+8r5x+rS5klZNf7FAiAEA0W93IbL7+QJD3HSBKHPEMqrX2UzV+L+khy8o0PL9U0gX7ff52U88FEPlKff6jGgcA4Hhs8VU3Os77DhAduAYPwHHpkpLY6Hiql5l7AQAtr6n3F953gOhAQQVwzLbtrlZt3aHT/XvcLmVnpTuQCAAQ7bKz0uVxuw4YM5JuO3uAM4EAtCgKKoBjsqemTjc8M1+BoNVPvjVIaV6PjKQ0r0fTJw5logoAwAkxISNN0ycO3fe+07ltguKM9HphiaoDrJEKRDpjbfit6JKZmWkLCgqcjgGgCXXBkG6cXaBPV2/X05MydXZ6V6cjAQBi2DtFW3T7iwv1nSHd9ferRiourrHpUwCEC2PMAmttZmOPcQYVwFGx1upXeUv18cpy3TNhCOUUAOC47wztoV9ecIreLtqq6e8sdzoOgONwxFl8AWB/j360Ri/P36QfnnuSrhzd2+k4AABIkm48o5827/LrqU/WKc3r0eTx/ZyOBOAYUFABNNvrhZuVk1+siRlp+vE3BzkdBwCAfYwx+vWFg1Xq8+t3by5Tqtejb53a3elYAI4Sl/gCaJbP12zXz15bonH9O+ne7w+TMdzfAwAIL644o4euzNDwnl7d+XKhCjfucjoSgKNEQQVwRCu3VeqW5xaoX+c2evy605QQz0sHACA8eRJcenpSprq1S9LU2QXasGOv05EAHAWOMgEc1rbd1Zo8c548bpdmTRmt9h6305EAADisTm0TNWvyKIWs1eRZ87Vzb63TkQA0EwUVQJO+Xuu0wh/QzMmjlOb1OB0JAIBm6d+lrWZMylSJz6+bni1gjVQgQlBQATSqLhjSHS8s1IqtlXrkmpEaktbe6UgAAByV0/p01ENXjNDCjbt09yuLFApZpyMBOAIKKoBDsNYpACBafL1G6jtLt+pPb7NGKhDuWGYGwCFY6xQAEE2+XiN1xqfrlNbBoymskQqELQoqgAOw1ikAINrsv0bq7xvWSM1ijVQgLFFQgRiXV1iinPxilfr86tQ2QTv31rLWKQAg6ny9RupVT32pO18q1G1nD9A/Cjar1OdXqtej7Kx0TchIczomEPO4BxWIYXmFJZqWW6QSn19W0vY9tbJW+t6IHqx1CgCIOl+vkdo20aW/vr9q3/tfic+vablFyisscToiEPM4AgViWE5+sfwHTbtvJT3y4RpnAgEAcIJ1apsoV9yhh8D+QFA5+cUOJAKwPwoqEMNKff6jGgcAIBqUV9Y0Os77H+A8CioQw3q0T2p0PNXraeUkAAC0nqbe53j/A5xHQQViVChk1b3doQXV43YpOyvdgUQAALSO7Kx0edyuA8ZccUY//Raz1wNOo6ACMWr6O8u1cJNPFw7roTSvR0ZSmtej6ROHMoshACCqTchI0/SJQ/e9/7VNjFcwZLVhZ5XT0YCYxzIzQAya8claPfXJOk0+va9++73BLCcDAIg5EzLS9v1C1lqr7NeW6K/vr1L3dkm6cnRvh9MBsYuCCsSYfy0u1R/fWq4LhnbXry+knAIAYIzR9IlDVV5Zo1/mLVWXlESdd0o3p2MBMYlLfIEY8vma7frpq4s1ul9HPXD5CLniKKcAAEiS2xWnR68ZqVNT2+mOFxdq4cZdTkcCYhIFFYgRy0p365ZnF6hv52Q9dV2mkg6aHAIAgFjXJjFeMyePUrd2SbrxmflaW77H6UhAzKGgAjFg864qTZ41T22T4jX7htFqn+x2OhIAAGGpc9tEzZ4yWnHG6PqZ81RWWe10JCCmUFCBKOerqtWkmfNUHQjqmSmj1aM9a7wBAHA4fTu30czJo7RjT62mzJqvyuqA05GAmEFBBaJYdSCoG2cXaNMuv566PlPp3VOcjgQAQEQY3surR68dqRVbK3Xb8wtVWxdyOhIQEyioQJQKhqzufKlQCzfu0kNXjNCY/p2cjgQAQEQ5J72r7p04VJ+u3q6fvbZYoZB1OhIQ9VhmBohC1lr95o2l+s+ybfrdRafqO0N7OB0JAICIdFlmL5VV1ignv1jd2idp2ndOcToSENUoqEAUemTOar0wd6NuPWuAJp3e1+k4AABEtNvPHqCtFdV64uO16t4uSVPG93M6EhC1KKhAlHm1YJP+8p+VmpiRpv/5drrTcQAAiHjGGP3vRaeqrLJav39zmbqmJOm7w7g6CTgRuAcViCJzVpRpWm6RzhzYWX++dJiMMU5HAgAgKrjijB66MkOn9e6gu19ZpC/X7nA6EhCVOIMKRLi8whLl5Ber1OeXJKV5k/TYtafJ7eL3TwAAtKQkt0szJmXq0se/0KSZc9Xek6Dyyhqlej3KzkrXhIw0pyMCEY8jWCCC5RWWaFpukUp8fllJVtL2PbV6f9k2p6MBABCVvMkJumZMb9XUWZVV1shKKvH5NS23SHmFJU7HAyIeBRWIYDn5xfIHggeMVdeFlJNf7FAiAACi34xP1h0y5g8Eef8FWgAFFYhgX1/W29xxAABw/Hj/BU4cCioQoapq6xTvanwSpFSvp5XTAAAQO5p6n+3RPqmVkwDRh4IKRKDaupBueW6BAkGrhINKqsftUnYWy8sAAHCiZGely+N2HTLeuW2C6oIhBxIB0YOCCkSYYMjq7lcX6ZNV23Xf94fpvkuHK83rkZGU5vVo+sShzCIIAMAJNCEjTdMnDj3g/feSEalaUrJbv3i9SNZapyMCEYtlZoAIYq3Vr/KW6q0lW/TLC07R5aN6SRKFFACAVjYhI+2Q999endro4Q9WyZucoGnfOZn1yIFjQEEFIkhOfrFemrdRt589QDd9o7/TcQAAwH7uPn+gKqpq9eR/18qb7NbtZ5/kdCQg4lBQgQjxxMdr9OhHa3T1mN7cYwoAQBgyxui33ztVPn9A971brPYet64Z08fpWEBEoaACEeCV+Rs1/Z0VunBYD/3h4iFcMgQAQJiKizP6y2XDVVldp1/lLVW7JLe+NzzV6VhAxGCSJCDMvVO0RdNyi3TWoC564PIRcsVRTgEACGduV5wevWakRvXpqB+/ukgfFZc5HQmIGBRUIIx9umq77np5kTJ6d9Bj145UQjw/sgAARIIkt0szJmdqYNcU3fr8Ai3YsNPpSEBE4GgXCFOFG3fp5ucK1L9LG82cNErJCVyRDwBAJGmX5NazN45Wj/YeTZk1X8u37HY6EhD2KKhAGCreWqnJs+arS0qinr1htNonu52OBAAAjkHntol67sbRSk6I13VPz9P67XudjgSENQoqEGY27azSdU/PVWJ8nJ6/cYy6tktyOhIAADgOPTsk6/mpoxUMhXTt03O1bXe105GAsHXEgmqM6WWMmWOMWW6M+coYc9dBj//UGGONMZ2b2P7uhu2WGmNeMsZwtA00oayyWtc+PVc1dSE9d+MY9eqY7HQkAADQAk7qmqLZN4zWrr21uu7pufJV1TodCQhLzTmDWifpJ9baUySNlXSHMWawVF9eJX1T0sbGNjTGpEm6U1KmtXaIJJekK1siOBBtKqoCuv7peSqvrNGsKaOU3j3F6UgAAKAFDevp1VOTMrV+R5Umz5qvvTV1TkcCws4RZ12x1m6RtKXh40pjzHJJaZKWSXpQ0s8kvXGE7+ExxgQkJUsqPd7QQLTIKyxRTn6xSn1+uV1xqguFNPuG0RrZu4PT0QAAwAlw+oDO+ttVGbrt+QW65JHPtKemTlsqqpXq9Sg7K10TMtKcjgg46qjuQTXG9JWUIWmuMeYiSSXW2sVNPd9aWyLpL6o/w7pFUoW19j/HHheIHnmFJZqWW6QSn19WUm0wpPg4ox17uOQHAIBolnVqd10xqpdWlu1RaUW1rKQSn1/TcouUV1jidDzAUc0uqMaYtpL+KelHqr/s95eSfnOEbTpIulhSP0mpktoYY65t4rk3G2MKjDEF5eXlzY0FRKyc/GL5A8EDxmqDVjn5xQ4lAgAAreW/K7cfMuYPBDkOQMxrVkE1xrhVX05fsNbmShqg+tK52BizXlJPSQuNMd0P2vR8SeusteXW2oCkXEmnN/Y9rLVPWmszrbWZXbp0Oba/DRBBSn3+oxoHAADRg+MAoHHNmcXXSHpa0nJr7QOSZK0tstZ2tdb2tdb2lbRZ0khr7daDNt8oaawxJrnh65wnaXmL/g2ACBQMWSW5XY0+lur1tHIaAADQ2pp6v0/1suAFYltzzqCOl3SdpHONMYsa/lzQ1JONManGmLclyVo7V9JrkhZKKmr4fk8ef2wgcoVCVj97bYn8gaDi48wBj3ncLmVnpTuUDAAAtJbsrHR5Gvlldb/ObWStdSAREB5MOP4AZGZm2oKCAqdjAC0uFLKallukVwo26e7zB6lPp+R9s/gyex8AALFl/9n8U71JGtStreYUb9cN4/vp1xeeovoLEIHoY4xZYK3NbOyxIy4zA6BlWGv1qzeW6pWCTbrz3JN01/kDJYlCCgBAjJqQkXbAcYC1Vr9/c5lmfrZO8S6jad85mZKKmENBBVqBtVa//ddXenHuRt129gDd/c1BTkcCAABhxhij31w4WMGQ1ZP/Xav4OKPsrHRKKmIKBRU4wb7+beizX2zQzd/or5/xRgMAAJpgjNHvLjpVdSGrRz9ao3hXnH7ML7YRQyiowAlkrdX0d1Zo1mfrNWV8Xy7VAQAAR2SM0R8vHqJg0OrhD1YpPs7ozvMGOh0LaBUUVOAEsdbqvvxiPfnftbp+XB/95sLBlFMAANAscXFG0ycOVV3I6oH3VsoVZ3THOSc5HQs44SiowAny4Hsr9dhHa3T1mN763UWnUk4BAMBRiYszuu/SYQqGQsrJL5bbZXTzNwY4HQs4oSiowAnw0Pur9PCHq3VFZi/98eIhlFMAAHBMXHFGf7lsuIJW+tPbK+SKi9ONZ/RzOhZwwlBQgRb2yJzVevD9lfr+yJ6aPnGo4uIopwAA4NjFu+L04OXDFQyF9Ic3lyk+zmjS6X2djgWcEBRUoAU98fEa5eQXa8KIVN136TDKKQAAaBHxrjg9dGWG6oIL9dt/fSVXnNG1Y/s4HQtocXFOBwCixYxP1mr6Oyv0veGp+stlw+WinAIAgBbkdsXp71eP1PmndNWv8pbq5XkbnY4EtDhjrXU6wyEyMzNtQUGB0zGAI8orLFFOfrFKfH5J0vCe7fTP28Yr3sXvfgAAwIlRUxfUrc8t0JzicnmT3aqoCijV61F2VromZKQ5HQ84ImPMAmttZmOPcRQNHKO8whJNyy3aV04lqXjbHr25ZIuDqQAAQLRLjHfpgqE9FGckX1VAVlKJz69puUXKKyxxOh5wXCiowDHKyS+WPxA8YKw6UD8NPAAAwIn01/dXKXTQhZD+QJDjEEQ8CipwDKy1B5w53V9pE+MAAAAtpanjDY5DEOkoqMBRstbqnreWN/l4qtfTimkAAEAsaup4o21SvMJxjhmguSiowFEIhqym5RZpxqfrdObAzvK4D/wR8rhdys5KdygdAACIFdlZ6fK4XQeMuYxRZXWd/vjWckoqIhbroALNFAiGdPcri/Tmki364bkn6cffHKQ3FpUqJ79YpT4/s+cBAIBW8/Xxxv7HIT/95iAtKa3Q05+u057qOv1p4lCWvUPEYZkZoBmqA0Hd/sJCfbiiTNO+c7JuOWuA05EAAAAOYa3Vg++v0sMfrNJ3h/XQg5ePUEI8F00ivBxumRnOoAJHUFkd0NTZBZq3fqfuuWSIrhnTx+lIAAAAjTLG6MffHKR2SfH641vLtbemTo9dc5o8Ca4jbwyEAX6dAhzGrr21unbGXBVs2KW/XjGCcgoAACLC1DP7696JQ/XxynJNmjVPldUBpyMBzUJBBZpQtrtaVzz5hZZvrdQT156mi0dwbykAAIgcV47urYevzNDCDbt0zYy52rm31ulIwBFRUIFGbNpZpUsf/0Ilu/x6ZsoonT+4m9ORAAAAjtr3hqfqyetPU/HWSl3xxBfatrva6UjAYVFQgYOsLqvUZY9/oQp/QM9PHaPTB3R2OhIAAMAxO/fkbnpmymiV+vy69PHPtXFHldORgCZRUIH9LC2p0OVPfKm6kNUrt4xVRu8OTkcCAAA4buMGdNILN41VZXWdLnvic63aVul0JKBRFFSgwfz1O3XVk1/K43bpH7eO08nd2zkdCQAAoMWM6OXVKzePU8hKlz/xhYo2VzgdCTgE66AiZuUVluxb3LpjmwRV+GvVu2MbPT91jFK9HqfjAQAAnBDrt+/VNTPmarc/oMnj+yp3YYlKfX6lej3KzkrXhAwmhsSJdbh1UDmDipiUV1iiablFKvH5ZSXt2FurYEiackZfyikAAIhqfTu30Wu3jZMnIU5/+3D1vuOhEp9f03KLlFdY4nRExDAKKmJSTn6x/IHgAWNW0uMfrXUmEAAAQCvq0d6juLhDq4A/EFROfrEDiYB6FFTEpFKf/6jGAQAAos22isaXnOF4CE6ioCLm+GuDSnQ3vutzeS8AAIgVTR339Gif1MpJgP9HQUVMKdtdrSuf/ELVgZDcLnPAYx63S9lZ6Q4lAwAAaF3ZWenyuF2HjCe5Xdq1t9aBRAAFFTFk+ZbdmvDIZ1pVtkdPXZ+pnEuHK83rkZGU5vVo+sShzFoHAABixoSMNE2fOPSA46HrxvXRZp9flzz6mdaU73E6ImIQy8wgJsxZUaYfvLhQKUluzZiUqSFp7Z2OBAAAEJYWbNilm58tUF3I6vFrT9O4AZ2cjoQowzIziGnPfLZON86er76d2yjvjvGUUwAAgMM4rU8H5d0xXl1TEnXd03P1asEmpyMhhlBQEbXqgiH95o2l+t9/L9N5p3TTP24dp+7c9A8AAHBEvTom65+3n65xAzrpZ68t0b3vrFAoFH5XXiL6UFARlSqrA7pxdoGe/WKDbv5Gfz1+7WlKToh3OhYAAEDEaJfk1szJo3TNmN56/OM1uv2FhfLXBo+8IXAcKKiIOpt3VenSx77QZ6u3a/rEofrFBafIFWeOvCEAAAAO4HbF6Y8ThujXFw5W/rKtuuLJL1S2u/H1U4GWQEFFVCncuEsTHvlcpRV+zb5htK4a3dvpSAAAABHNGKMbz+inp67L1OqyPZrwyGdaVrrb6ViIUhRURI23lmzRlU9+qeQEl16//XSNP6mz05EAAACixvmD6+f0CFnpssc/14crtjkdCVGIZWYQsfIKS5STX6xSn18pSfHaXV2nzD4d9MR1p6lT20Sn4wEAAESlbburNXV2gb4qrdCEEamau26nSn3VSvV6lJ2VzrryOCKWmUHUySss0bTcIpX4/LKSdlfXyWWMrhzVi3IKAABwAnVrl6RXbhmrwT3aKbewVCW+allJJT6/puUWKa+wxOmIiGAUVESknPxi+QMHziIXtFYPvr/KoUQAAACxIzkhXjurag8Z9weCyskvdiARogUFFRGpxOdvdLy0iXEAAAC0rC2+xmfz5XgMx4OCiogSDFk9+N7KJh9P9XpaMQ0AAEDsauq4KyUpXqFQ+M1zg8hAQUXE2La7WtfM+FIPfbBKmX28SnIfuPt63C5lZ6U7lA4AACC2ZGely+N2HTAWZ+rnBrlh9nzt2FPjUDJEMgoqIsJ/V5brgoc+0eJNFfrLZcP12m3jde/EYUrzemQkpXk9mj5xKLPGAQAAtJIJGWmaPnHoAcdj9182XH+YMESfr9mhCx7+RHPX7nA6JiIMy8wgrNUFQ3rgvZV69KM1Su+Wor9fnaGB3VKcjgUAAIDD+Kq0Qj94sVAbduzV3ecP0u3nnCRXnHE6FsIEy8wgIpX6/LryyS/16EdrdOWoXsq7YzzlFAAAIAKcmtpe//7hGfre8FTd/95KTZo5T+WVXPKLI6OgIix9uGKbLnj4Ey3fslsPXTlC935/mDwJriNvCAAAgLDQNjFef71ihP78/aEq2LBT33noE322ervTsRDmKKgIK7V1Id3z1jLd8EyBUtt79O8fnqGLR3BfKQAAQCQyxuiKUb31xh1nyJvs1rVPz9UD761UkFl+0QQKKsLGpp1VuvyJL/TUJ+t03dg+yr39dPXv0tbpWAAAADhO6d1T9K8fjNelI3vq4Q9W6eqnvtTWisbXUUVsO2JBNcb0MsbMMcYsN8Z8ZYy566DHf2qMscaYzk1s7zXGvGaMWdHwNca1VHhEj3eXbtV3H/5Ea8r26NFrRuoPE4Yoyc0lvQAAANEiOSFeOZcN1wOXD1dRSYUuePgTfVRc5nQshJn4ZjynTtJPrLULjTEpkhYYY96z1i4zxvSS9E1JGw+z/UOS3rXWXmqMSZCUfPyxEenyCkuUk1+sUp9fyYku7a0JaljP9vr7VSPVuxO7CAAAQLSaOLKnhvX06gcvLtTkWfN13sldtXzrbm3xVSvV61F2VjpLB8awI55BtdZusdYubPi4UtJySV/vMQ9K+pmkRi8iN8a0k/QNSU83bF9rrfUdf2xEsrzCEk3LLVKJzy8raW9NUK44o+vH9qGcAgAAxICTurZV3h3jNW5AR32wokylvmpZSSU+v6blFimvsMTpiHDIUd2DaozpKylD0lxjzEWSSqy1iw+zSX9J5ZJmGWMKjTEzjDFtmvjaNxtjCowxBeXl5UcTCxHmvvwV8geCB4wFQ1YPvr/KoUQAAABobUlulzbu8B8y7g8ElZNf7EAihINmF1RjTFtJ/5T0I9Vf9vtLSb85wmbxkkZKesxamyFpr6SfN/ZEa+2T1tpMa21mly5dmhsLEWbhxl0q9TV+Q3yp79AXKAAAAESvpo7/OC6MXc0qqMYYt+rL6QvW2lxJAyT1k7TYGLNeUk9JC40x3Q/adLOkzdbauQ2fv6b6wooYs7s6oF/nLdX3H/tccabx56R6Pa0bCgAAAI5q6vjPSvrjm8u0t6audQPBcc2Zxdeo/h7S5dbaByTJWltkre1qre1rre2r+iI60lq7df9tGz7fZIxJbxg6T9KylvwLILxZa/V20Radf//HemHuBk0+va/+dMlQeQ6aodfjdik7K72JrwIAAIBolJ2VfshxYZI7TuMGdNSMT9fpWw/+Vx8s3+ZQOjihObP4jpd0naQiY8yihrFfWGvfbuzJxphUSTOstRc0DP1Q0gsNM/iulTTl+CIjUmzeVaXfvPGVPlxRplNT22nGpEwN6+mVVH/Pwdez+DJbGwAAQGz6+vivsePCgvU79YvXi3Tj7AJdMLS7fvu9U9WtXZLDiXGiGWsbnYDXUZmZmbagoMDpGDhGdcGQZn22Xg+8t1LGSD/+5iBNPr2v4l1HNScXAAAAYlxtXUhPfbJWD3+wSgmuOP3s2+m6ekwfuZq6ZwwRwRizwFqb2ehjFFS0pMWbfJqWW6RlW3br/FO66ncXD1Ea95YCAADgOGzYsVe/yluqT1Zt14heXv3pkqEanNrO6Vg4RhRUnHCV1QHd/5+Vmv3FenVNSdTvLjpVWad2V/0tzAAAAMDxsdbqX4tL9ft/L5PPH9DUM/rprvMHKjmhOXctIpwcrqDyr4mjlldYst99Akn61qnd9U7RVm2rrNb1Y/vop1npSklyOx0TAAAAUcQYo4tHpOmsQV107zsr9MR/1+qtoi36w4QhqqgKML9JlOAMKo5KXmGJpuUWyR8IHjCe2j5Jj157mkb08joTDAAAADFl3rr6SZRWl+2Ry0jB/WqNx+3S9IlDKalh6nBnUJm1BkclJ7/4kHIqSTKinAIAAKDVjO7XUW/feaZSkuIPKKeS5A8ElZNf7EwwHBcKKprNV1WrEp+/0ce2+KpbOQ0AAABiXUJ8nPZU1zX6WGkTx60IbxRUHNGuvbX6S36xzvjznCafk8pMvQAAAHDA4Y5DH/5glXZXB1oxDY4XBRVN2rW3Vjn5K3TGnz/UIx+t1lnpXfQ/306Xx+064Hket0vZWekOpQQAAEAsy8469Pg0MT5Op6a20wPvrdQZ936oh95fpQo/RTUSMIsvDrFzb61mfLJWsz9fr6pAUN8d2kN3njdQg7qlSJJ6tPcwSxoAAADCwtfHoY0dny4tqdDDH6zSg++v1IxP1+rGM/ppyvh+au9hxYlwxSy+2Gfn3lo91VBM/Y0UUwAAACASfV1U/7Nsm1KS4nXD+H664QyKqlMON4svBTUGHbiOqUe3nd1fm3dV69kv6ovphcNSdee5J2kgxRQAAABR5KvS+qKa/9U2pSTGa8r4vrrxjP6aU1zGFYKtiIKKfZpax1SSLhqeqjvPO0kndaWYAgAAIHotK92thz9YpXe/2qpEl1HQSnWh/+9FrKN6YrEOKvbJyV/RaDntmpKoh6/KoJwCAAAg6g1ObafHrztN79x1pkycOaCcSqyj6iQmSYoRWyr8emNRqUqaWK+0vLKmlRMBAAAAzjqlRzvVBEKNPlbi86usslpdU5JaOVVso6BGscrqgN5dulWvF5boi7U7ZK3kdhkFgode1s06pgAAAIhFqV6PSnz+Rh8b+6cPdMbALrokI1VZp3ZXcgL16UTj/3CUCQRD+mRVuV4vLNV7y7aqOhBSn07Juuu8gbokI02FG32H3IPKOqYAAACIVdlZ6Y0eH//o/IGqrK7T64UluvuVxUpOWKpvn9pdEzLSNP6kznLFGQdTRy8KaoQ5eAbe7Kx0XTwiVUs2V+j1whL9e3GpduytVYdkty47rZcuGZmmjF5eGVP/A9SnUxtJja8TBQAAAMSaw62jKkk//uYgFWzYpdcLN+vNJVuUW1iirimJunhEqiZkpGlwj3YyxjR6nM4x9tFjFt8I0tgMvPFxRh3auFVeWauE+Didf0pXXZLRU2cN6qKEeObAAgAAAFpKdSCoOSvKlFtYoo+KyxQIWqV3S9HAbm31/rJtqq77//tZmQm4aSwzEyXGTf9AWyoOneQowRWn3198qr4ztAeLDQMAAACtYNfeWr1ZtEWvL9yshRt9jT4nzevRZz8/t3WDRQAKagQKhazWbt+jhRt9Kty4S4UbfVqxtbLR5xpJ6+79busGBAAAACBJ6vfzt9RUq7pmTG+N7N1BGb296te5zb5b72LZ4Qoq96C2sqauTfdV1apwk0+FDYV00SafKqvrJEntkuKV0buDSnz+fWP7YwZeAAAAwDlNzQScGB+nfy0q1QtzN0qSvMlujejlVUavDhrZx6vhvbxql1R/BST3sNbjDOpRON6dprF7SF3GqGPb+ntIJSnOSOnd2ymjt1cZvbzK6N1B/Tu3UVycaXR7rm0HAAAAnHW44/SLhqdqdfmefVdFFm70aWVZpayVjJFO6tJWHdu4tXCj74DlII/lOD9SSi6X+LaApna6P1x8qs4c1EW7qmq1a29AFf5a7aoKaFdVrXxVAe3aWyufP1B/hnSjT3WhQ/9/J8XH6YfnDVRGb6+G9fSqbWLTJ7YjZacDAAAAYsnRHKdXVge0eFNFfWnd5NNHxWVqpCYowRWns9K7qEOyWx2SE+RNTlCHZLe8yQnyNox9/fnbRVsi5mQWBbUFjL/3wyYX8G1KYnxcw47kljfZrS/X7mz0edxDCgAAAMSuw93DenL3lPqTYVUB1e43S/DBjNTo1wjHiZq4B7UFlB6mnP5xwpB9v71ov+83GQnyJLgOeF5TJZd7SAEAAIDY1dQ9rGlej9790TckSdZa+QPB+qs199aqwh/YV1x9e2t1/3srG/3ah+sx4YiC2kyH22muHdunWV8jOyu90dPu2VnpLZYTAAAAQGRpTk8wxig5IV7JCfFKa+QE18vzN0XFybA4pwNEiuysdHncB54RPdpyOSEjTdMnDlWa1yOj+nIbjteEAwAAAGg9LdETWqKvhAPuQT0KTFAEAAAAIFxFSl9hkiQAAAAAQFg4XEHlEl8AAAAAQFigoAIAAAAAwgIFFQAAAAAQFiioAAAAAICwQEEFAAAAAIQFCioAAAAAICxQUAEAAAAAYYGCCgAAAAAICxRUAAAAAEBYoKACAAAAAMKCsdY6neEQxphySRucznEYnSVtdzoEcBD2S4Qb9kmEI/ZLhCP2S4SbE71P9rHWdmnsgbAsqOHOGFNgrc10OgewP/ZLhBv2SYQj9kuEI/ZLhBsn90ku8QUAAAAAhAUKKgAAAAAgLFBQj82TTgcAGsF+iXDDPolwxH6JcMR+iXDj2D7JPagAAAAAgLDAGVQAAAAAQFigoB4lY8y3jTHFxpjVxpifO50HsckYM9MYU2aMWbrfWEdjzHvGmFUN/+3gZEbEFmNML2PMHGPMcmPMV8aYuxrG2S/hCGNMkjFmnjFmccM++buGcfZJOM4Y4zLGFBpj3mz4nP0SjjLGrDfGFBljFhljChrGHNkvKahHwRjjkvSIpO9IGizpKmPMYGdTIUY9I+nbB439XNIH1tqBkj5o+BxoLXWSfmKtPUXSWEl3NLw+sl/CKTWSzrXWDpc0QtK3jTFjxT6J8HCXpOX7fc5+iXBwjrV2xH7LyziyX1JQj85oSauttWuttbWSXpZ0scOZEIOstf+VtPOg4YslzW74eLakCa2ZCbHNWrvFWruw4eNK1R94pYn9Eg6x9fY0fOpu+GPFPgmHGWN6SvqupBn7DbNfIhw5sl9SUI9OmqRN+32+uWEMCAfdrLVbpPqyIKmrw3kQo4wxfSVlSJor9ks4qOEyykWSyiS9Z61ln0Q4+Kukn0kK7TfGfgmnWUn/McYsMMbc3DDmyH4Z3xrfJIqYRsaYBhkAGhhj2kr6p6QfWWt3G9PYyybQOqy1QUkjjDFeSa8bY4Y4HAkxzhhzoaQya+0CY8zZDscB9jfeWltqjOkq6T1jzAqngnAG9ehsltRrv897Sip1KAtwsG3GmB6S1PDfMofzIMYYY9yqL6cvWGtzG4bZL+E4a61P0keqv3effRJOGi/pImPMetXfKnauMeZ5sV/CYdba0ob/lkl6XfW3NjqyX1JQj858SQONMf2MMQmSrpT0L4czAV/7l6RJDR9PkvSGg1kQY0z9qdKnJS231j6w30Psl3CEMaZLw5lTGWM8ks6XtELsk3CQtXaatbantbav6o8jP7TWXiv2SzjIGNPGGJPy9ceSviVpqRzaL421XKF6NIwxF6j+3gGXpJnW2nucTYRYZIx5SdLZkjpL2ibpt5LyJL0qqbekjZIus9YePJEScEIYY86Q9ImkIv3/fVW/UP19qOyXaHXGmGGqn9TDpfpfyL9qrf29MaaT2CcRBhou8f2ptfZC9ks4yRjTX/VnTaX6W0BftNbe49R+SUEFAAAAAIQFLvEFAAAAAIQFCioAAAAAICxQUAEAAAAAYYGCCgAAAAAICxRUAAAAAEBYoKACAAAAAMICBRUAAAAAEBYoqAAAAACAsPB/AE6XVB6PidEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot initial profile\n",
    "x = np.arange(0, domain_size[0], 1)   # start,stop,step\n",
    "y = fluid_depth + amplitude * np.cos(x / (domain_size[0] - 1) * np.pi * 2 + np.pi)\n",
    "\n",
    "plt.plot(x,y, marker='o')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "initialize_phasefield()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA54AAAFlCAYAAACDRTcUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXmElEQVR4nO3db6ykZ3kf4N99zpp/oSS4XjsudrNEcmkoLZCuXFqkqo1x6iQIW5VAIBGtWiSrEk1JlSo15EM/VXLVKgWp6QcLaLcKDTgEZAupKe4mKIqUUNZAAtRJsBzHOGy9JxQ3VJUwu3v3wxmatX3e2TNn59nZmXNdkjUzz7x/7jnP+pzzm3fO/VR3BwAAAEbZWnUBAAAAbDbBEwAAgKEETwAAAIYSPAEAABhK8AQAAGAowRMAAIChjlzJk1133XV97NixK3lK4BD6g4cfW3UJAFedv/TXf3DVJQCHwMMPP/wn3X30ueNXNHgeO3Ysp0+fvpKnBA6h27feuuoSAK46D53+5VWXABwCVfVHe437qC0AAABDCZ4AAAAMJXgCAAAwlOAJAADAUPtqLlRVjyf5VpLzSc519/GqujbJx5IcS/J4krd19zfHlAkAAMC6WuSK59/t7td19/HZ43uSnOruW5Kcmj0GAACAZ7mcj9remeTk7P7JJHdddjUAAABsnP0Gz07y6ap6uKruno3d0N1nkmR2e/1eO1bV3VV1uqpO7+zsXH7FAAAArJV9/Y1nkjd299er6vokD1XV7+33BN19X5L7kuT48eN9gBoBAABYY/u64tndX5/dnk3yySS3Jnmqqm5Mktnt2VFFAgAAsL4uGTyr6nuq6s99936SH03y5SQPJjkx2+xEkgdGFQkAAMD62s9HbW9I8smq+u72/7m7f7WqPpfk/qp6V5Inkrx1XJkAAACsq0sGz+5+LMlr9xj/RpLbRhQFAADA5ric5VQAAADgkgRPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoQRPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoQRPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoQRPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoQRPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKH2HTyraruqvlBVn5o9vraqHqqqr85uXz6uTAAAANbVIlc835PkkYse35PkVHffkuTU7DEAAAA8y76CZ1XdlOQnknzwouE7k5yc3T+Z5K6lVgYAAMBG2O8Vz/cn+dkkFy4au6G7zyTJ7Pb65ZYGAADAJrhk8KyqNyc5290PH+QEVXV3VZ2uqtM7OzsHOQQAAABrbD9XPN+Y5C1V9XiSjyb5kar6xSRPVdWNSTK7PbvXzt19X3cf7+7jR48eXVLZAAAArItLBs/ufm9339Tdx5K8Pcmvdfc7kzyY5MRssxNJHhhWJQAAAGvrctbxvDfJ7VX11SS3zx4DAADAsxxZZOPu/kySz8zufyPJbcsvCQAAgE1yOVc8AQAA4JIETwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoQRPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoQRPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoQRPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIYSPAEAABhK8AQAAGAowRMAAIChjqy6AICNUrXqCtivWrP3XvvCqitgv7pXXQHAVWfNfuoCAACwbgRPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIaynAqweba2V13BRqitq3BpmHVbAmWplvjv+ipdmqUvbMgyJFfh/zoAq3aYf4IDAABwBQieAAAADCV4AgAAMJTgCQAAwFCCJwAAAENdsqttVb0oyW8keeFs+49397+oqmuTfCzJsSSPJ3lbd39zXKkA+1PbV2lX2xV2ia1a4rm3lvie5dVa14Kmvr7dK+7SemGJ3WuX+Vom6jrIv4aVfo03pQsvwBWwn5/S307yI9392iSvS3JHVb0hyT1JTnX3LUlOzR4DAADAs1wyePau/zN7eM3sv05yZ5KTs/GTSe4aUSAAAADrbV+fS6qq7ar6YpKzSR7q7s8muaG7zyTJ7Pb6iX3vrqrTVXV6Z2dnSWUDAACwLvYVPLv7fHe/LslNSW6tqtfs9wTdfV93H+/u40ePHj1gmQAAAKyrhToxdPfTST6T5I4kT1XVjUkyuz277OIAAABYf/vpans0yXe6++mqenGSNyX5V0keTHIiyb2z2wdGFgqwX1svftGqS1jMQbrd1oIdXKfOcZCuslPnnvM6Fu6qO7X9QTrXHug1Lqfb7oGOcpAurYvuM6/b7YLHmuwqO6/jay/YbXdOTbVoZ9lFz53oXguwBJcMnkluTHKyqraze4X0/u7+VFX9VpL7q+pdSZ5I8taBdQIAALCmLhk8u/t3k7x+j/FvJLltRFEAAABsjtWttg0AAMChIHgCAAAwlOAJAADAUIInAAAAQ+2nqy3AWqnvfdn4kxxoCZQlLSky57meqmtqGZJ551jwWD3nWL292HIuk8eaOs7cY03vMmnB+Z2qtw6yNMoBlu6oqV2mzn9+zvIkU/sseKy5r31qOZep136AY00us7Lo65vnCs0vwCZwxRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoXS1BTbOd2768+NPMtX0dF5n1wX3mewEm+nutT3xduLUsSa74CaTb01eOLL4sRava7HjzN1nak7mvPTpuZreZ2ETzU2nO9ROH2qqg2xNNI+dGp+7z1T32slzTBe8dW7iuQMca1l1TR0nyWT32oW7CSdz5xFgk7niCQAAwFCCJwAAAEMJngAAAAwleAIAADCU4AkAAMBQutoCG+fpV71kacdatLvp5PbJ5Ft9U91YL8zpajvZcXZ74hwLjs/dZ7JD7eLHurC9d4vPyWPN62o7NSdT+2xNtxedO4+DTXZKvTBdVE08N32sOec/v/f41vm9v5BT20+NJ/M65x7gWJP1Tuwwce6tOV1tJ7sAT72OeZ1rF+1mDLAhXPEEAABgKMETAACAoQRPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEspwJsnG+8donrEiy6bMqcJTqml1OZWl9h+lCZWIZkary29173oaaOk2R7Yp+tifGp7ZPkyNQ+E2tITG6/NX2OqeemvoxH5hxrSi1xzYtecM2Wcxem3yuequr8xD5T40lybmLZlPMT9Z6f2n5iPEkuLLhPn5/+WvXUeab2mRqfM7VTy9VMLkszb+mbqfNYTgXYcK54AgAAMJTgCQAAwFCCJwAAAEMJngAAAAwleAIAADCUrrbAxvmLrzmz8D5bC3Yrndp+u6Y7pU7tM9Vd9UidnzzW1D4v2Np7n2smxo/MqfeFW+cmzj1xjjn1Tj03VdfU9ltzWn9eM1Hv9gHahW7N+bqMdqEXf0/4/ETv3u9c2PvH/IU5LZO/09sTx5oYn9p+YjxJzk0c69sT9Z6b8zWZquuZifGp7sDn5ta79z4Xpjr9zql3ap9lbZ8s/v0M4EpwxRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoXS1BTbO3/8LXxx+jqmup9tZvKvt1D7zOuROdX2dOtbk9vO68E4c6wWTHWenjzV1nqmOswf5+i7avXaVnWsPYpndbs/Ped956jyTx5rY/sKcczwz0UF2ap95XWKnuudOvcbJ7eecY+pYk11tD/D1Bdh0vvsBAAAwlOAJAADAUIInAAAAQwmeAAAADCV4AgAAMJSutsDG+RsveXRpx1q0U+pBjjXV7fZAx5oY317qORa3vXfzz0kHOscB9tl0e/cfnm/RXr/nD/C/yNQ5pjrnzj//3vtcmOzou/g5prvXLn6sKcs8FsDVyBVPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIYSPAEAABjKcirAxrl5+9sL77Ndq1vKYJnvAG5PLsmw+OvbquVVNl3X8mx5L/V5Liy8OMrizi9xyaELvbxjnV/iax//VZx2/gBfk1V+PwOY4qc0AAAAQwmeAAAADCV4AgAAMJTgCQAAwFCCJwAAAENdsqttVd2c5D8l+f7sNna7r7s/UFXXJvlYkmNJHk/ytu7+5rhSAfbne7desOoSluJq7Uype+w62V51AXua7La7wn/yB+kee0Vcnd8GABa2n98eziX5me7+oSRvSPLuqnp1knuSnOruW5Kcmj0GAACAZ7lk8OzuM939+dn9byV5JMkrktyZ5ORss5NJ7hpUIwAAAGtsoc9LVdWxJK9P8tkkN3T3mWQ3nCa5fmKfu6vqdFWd3tnZucxyAQAAWDf7Dp5V9dIkv5Lkp7v7T/e7X3ff193Hu/v40aNHD1IjAAAAa2xfwbOqrslu6PxId39iNvxUVd04e/7GJGfHlAgAAMA6209X20ryoSSPdPfPX/TUg0lOJLl3dvvAkAoBFnRNra6T55YWlFed7dr8Lrzne6JL7FVq6yrstntkif/rXshV2iEXYIUuGTyTvDHJTyb5UlV9cTb2vuwGzvur6l1Jnkjy1iEVAgAAsNYuGTy7+zczvYrUbcstBwAAgE2z+Z8/AgAAYKUETwAAAIYSPAEAABhK8AQAAGCo/XS1BVgrq1xOBVbhMCwZs058BwJ4Pj+pAAAAGErwBAAAYCjBEwAAgKEETwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoQRPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoQRPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoQRPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoQRPAAAAhhI8AQAAGOqSwbOqPlxVZ6vqyxeNXVtVD1XVV2e3Lx9bJgAAAOtqP1c8/2OSO54zdk+SU919S5JTs8cAAADwPJcMnt39G0n+13OG70xycnb/ZJK7llsWAAAAm+Kgf+N5Q3efSZLZ7fVTG1bV3VV1uqpO7+zsHPB0AAAArKvhzYW6+77uPt7dx48ePTr6dAAAAFxlDho8n6qqG5Nkdnt2eSUBAACwSQ4aPB9McmJ2/0SSB5ZTDgAAAJtmP8up/FKS30ryqqp6sqreleTeJLdX1VeT3D57DAAAAM9z5FIbdPc7Jp66bcm1AAAAsIGGNxcCAADgcBM8AQAAGErwBAAAYCjBEwAAgKEETwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoQRPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoQRPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoQRPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoQRPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIa6rOBZVXdU1e9X1aNVdc+yigIAAGBzHDh4VtV2kl9I8mNJXp3kHVX16mUVBgAAwGa4nCuetyZ5tLsf6+5nknw0yZ3LKQsAAIBNcTnB8xVJvnbR4ydnYwAAAPD/XU7wrD3G+nkbVd1dVaer6vTOzs5lnA4AAIB1dDnB88kkN1/0+KYkX3/uRt19X3cf7+7jR48evYzTAQAAsI4uJ3h+LsktVfXKqnpBkrcneXA5ZQEAALApqvt5n47d/85VP57k/Um2k3y4u//lJbbfSfJHBz7h8l2X5E9WXQRXlDk/fMz54WTeDx9zfjiZ98PHnF/9fqC7n/dR18sKnuuuqk539/FV18GVY84PH3N+OJn3w8ecH07m/fAx5+vrcj5qCwAAAJckeAIAADDUYQ+e9626AK44c374mPPDybwfPub8cDLvh485X1OH+m88AQAAGO+wX/EEAABgsEMbPKvqn1VVV9V1F429t6oerarfr6q/t8r6WK6q+tdV9XtV9btV9cmq+r6LnjPvG6qq7pjN66NVdc+q62H5qurmqvr1qnqkqr5SVe+ZjV9bVQ9V1Vdnty9fda0sV1VtV9UXqupTs8fmfMNV1fdV1cdnP88fqaq/ad43W1X909n39i9X1S9V1YvM+fo6lMGzqm5OcnuSJy4ae3WStyf5K0nuSPLvq2p7NRUywENJXtPdfy3JHyR5b2LeN9lsHn8hyY8leXWSd8zmm81yLsnPdPcPJXlDknfP5vmeJKe6+5Ykp2aP2SzvSfLIRY/N+eb7QJJf7e6/nOS12Z1/876hquoVSf5JkuPd/Zok29n9nc2cr6lDGTyT/NskP5vk4j9wvTPJR7v72939h0keTXLrKopj+br70919bvbwt5PcNLtv3jfXrUke7e7HuvuZJB/N7nyzQbr7THd/fnb/W9n9RfQV2Z3rk7PNTia5ayUFMkRV3ZTkJ5J88KJhc77BquplSf52kg8lSXc/091Px7xvuiNJXlxVR5K8JMnXY87X1qELnlX1liR/3N2/85ynXpHkaxc9fnI2xub5h0n+y+y+ed9c5vaQqapjSV6f5LNJbujuM8luOE1y/QpLY/nen903kC9cNGbON9sPJtlJ8h9mH7H+YFV9T8z7xuruP07yb7L7CcUzSf53d3865nxtHVl1ASNU1X9L8v17PPVzSd6X5Ef32m2PMS1/18i8ee/uB2bb/Fx2P5r3ke/utsf25n0zmNtDpKpemuRXkvx0d/9p1V7TzyaoqjcnOdvdD1fV31lxOVw5R5L8cJKf6u7PVtUH4iOWG232t5t3JnllkqeT/HJVvXOlRXFZNjJ4dveb9hqvqr+a3X+8vzP7peSmJJ+vqluzezXk5os2vym7l/NZE1Pz/l1VdSLJm5Pc1n+2jpB531zm9pCoqmuyGzo/0t2fmA0/VVU3dveZqroxydnVVciSvTHJW6rqx5O8KMnLquoXY8433ZNJnuzuz84efzy7wdO8b643JfnD7t5Jkqr6RJK/FXO+tg7VR227+0vdfX13H+vuY9n9JvbD3f0/kzyY5O1V9cKqemWSW5L89xWWyxJV1R1J/nmSt3T3/73oKfO+uT6X5JaqemVVvSC7DQkeXHFNLFntvov4oSSPdPfPX/TUg0lOzO6fSPLAla6NMbr7vd190+zn+NuT/Fp3vzPmfKPNflf7WlW9ajZ0W5L/EfO+yZ5I8oaqesnse/1t2f07fnO+pjbyiudBdPdXqur+7H4TO5fk3d19fsVlsTz/LskLkzw0u9r92939j8z75uruc1X1j5P81+x2wvtwd39lxWWxfG9M8pNJvlRVX5yNvS/JvUnur6p3ZfeXl7eupjyuIHO++X4qyUdmbyY+luQfZPciinnfQLOPVH88yeez+zvaF5Lcl+SlMedrqf7sE4cAAACwfIfqo7YAAABceYInAAAAQwmeAAAADCV4AgAAMJTgCQAAwFCCJwAAAEMJngAAAAwleAIAADDU/wPRyj9r/6QIQQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scalar_field(dh.cpu_arrays[\"C\"])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6IAAAFlCAYAAADxilWiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxyklEQVR4nO3de3zV9Z3n8ffnnFwhQCCEQBIQEBJAQIGordom2o54G6UUt9rdbjs7W+uM7XZ3prSwtp1exuos3dl2H7V1aKc7086MTK2I15ZabdB6B1HuAUQUEu4QIHBCknO++0cOGMIJOSEn53sur2cfPJLzPb9v8raP7wPyzu/7+/3MOScAAAAAAJIl4DsAAAAAACC7UEQBAAAAAElFEQUAAAAAJBVFFAAAAACQVBRRAAAAAEBSUUQBAAAAAEmV4+sbjxw50o0fP97Xt4/LiRMnNHjwYN8xkIVYe/CFtQdfWHvwhbUHnzJ9/a1Zs+agc6401nveiuj48eO1evVqX98+LvX19aqrq/MdA1mItQdfWHvwhbUHX1h78CnT15+ZvdfTe2zNBQAAAAAkFUUUAAAAAJBUFFEAAAAAQFJRRAEAAAAASUURBQAAAAAkFUUUAAAAAJBUFFEAAAAAQFJRRAEAAAAASUURBQAAAAAkVU5vB5jZzyXdImm/c256jPdN0g8l3STppKTPOefeTHTQZFqxtlFLVjaosTmkilef18K51Zo3q6LP85uaQyovLkyr+emcPZPmZ+PaS8R8AAAApIdei6ikf5L0I0m/6OH9GyVNjv65UtJPoh/T0oq1jVq8fL1C7WFJUmNzSIuXr5ekuH4gTuf56Zyd+cz3XYJ9zwcAAEgnvRZR59wLZjb+PIfcJukXzjkn6VUzKzazMc65PYkKmUxLVjac+UH4tFB7WN99apOGDcrtdf53n9qUtvPTOTvzzz+/2PP84YPzZD3MMzv//L99epNKh+R/MN9Of7Az8/+47YCWvviu2joikjpL7FcfXacdB1tUW1V6ZqJZ53Qzi37s/Dr1Dfv1oz9s16ku87/26DrtPRrSx6eNVsCkgJkC1vk1AgE7M2Ymrdy4V/c9vVmt7R/MX7R8nSIRp0/MrpBZT//1nfpbwk9/DYosAABIF9bZH3s5qLOIPtXD1tynJD3gnPtj9PVzkr7mnFsd49i7JN0lSWVlZXOWLVvWv/QD4HO/PeE7AoAME7TOPwGTgoHTn1vneEA6GHKKxPirOC8gzRkdVH7AlBeU8oLRj2ded47tONKh3+8KK9qjz8z93PQ8XVXe+y8RJOnlpnY9urVdh1ojKikI6JNVuXHPBRKhpaVFRUVFvmMgC7H24FOmr79rr712jXOuJtZ78WzN7U2sX/XHbLfOuaWSlkpSTU2Nq6urS8C3T6yKV59XY3PonPHSonwt/c9zep1/1y/W6EDLqbScn87ZmX/++f8Qx/wvDND8kUX5+ofPzI45p+vvwe7+lzU62NIWY36eHvx053zXbZ6Ljnz6p6/1mOsX/+UKOUnORY92nfOc6/w6TtLnf3HO783O+L93zpJzThHnFIlIEdc5N+I6y2PEOX19xYYe5/+36yapI+IUjji1h53CkYjaI07hsFNHxKkjEtHjbzXFnNsWkXa35inUFlFre1ih9g6FYzXWHub+bH2bnt+bpxGD8zSyqPPjiMH5XT7P08iifL327iH9cvNmhdqdJNOhVqdfbg5r2tRpnFVF0tTX1ysVfy5A5mPtwadsXn+JKKK7JY3t8rpSUuyfqtLAwrnVZ22Rk6TC3KDuvXmqZo0b3uv8e2+emrbz0zk7888/f7bH+V+/earmXDSi1/lfv3laD/On6cqJJeedW1FcGPMXSBXFhfroma25Fzb/1kvLe53/k/p3epz/V9dX9zp/9c4jPc5/8avXnTXWHo4o1B5Wa3tYrW2dn9/wgxdi/vYv4qTJo4p0qKVNDXuP6/CJNh052d5rHumDbdF11aUqHpQX1xwAAIB4JaKIPiHpi2a2TJ03KTqarteHSh9cj3XmzqV9vNaq6/wLuVbL5/x0zp5p87Nt7fV3fk+/QFo4t/cSmG7zc4MB5QYDGlrwwbbZ8vMU6Z/8p7PPZneEIzpysl2HT7TpUMspHTrRpi89vDZmroMtbbrsO89qzLACTR0zVFPHDIl+HKrxJYMVDHywIYZrVAEAQF/0eo2omT0sqU7SSEn7JP2NpFxJcs49FH18y48k3aDOx7f8WazrQ7urqalxq1f3ephX2XyqHH6x9vrO911rfT92KVaRvX/+jLi+xtUPxL4koWRwnu766ERt3nNMm/cc1zsHWtQR3RpckBtQ9eihmjp6iNrCET21bs+Zm0X19fsDEn/vwR/WHnzK9PVnZhd+jahz7s5e3neS7rnAbACQEPNmVfSr9KTz/P6eje7pjOw3bjn7GtFTHWFt29dypphu2XtMKzfujbndN9Qe1pKVDRRRAAAQUyK25gIAPEtUkT3ftvD8nKCmVwzT9IphZ8acc5q4+JmY16g2Noe0eudhzbloeK+PsAEAANmFIgoAOFNk+7pFyMx6vEbVJC146BVNHlWkT10+Vp+cXanhg7nxEQAAkAK+AwAA0tvCudUqzA2eNVaYG9QD82fo7z45Q4Pzc/S3T2/Wld97Tl96eK1e3n5QkTgfQwMAADITZ0QBAP3S2zWqn7p8nDbvOaZ/f2OXlr+5W0++3aSLSgbpU5eP1YI5lRo1pIC77gIAkGUoogCAfuvtGtWpY4bqW7deokU3TtFvNuzRw6/t0v/6bYP+/ndbNXXMEDXsbVFbuPOuu43NIS1evv7M1wUAAJmHIgoASJqC3KA+MatSn5hVqe37W/Tvb7yvf/zju+q+U5e77gIAkNm4RhQA4MWkUUW69+Zp6ulx1k0xboAEAAAyA0UUAOBVeXFhzPHRwwqSnAQAACQLRRQA4FWsu+5K0qmOsBr2HveQCAAADDSKKADAq3mzKnT//BmqKC6USaooLtR///hk5QQCmv/jl/S7jXt9RwQAAAnGzYoAAN7FuuvuHZeP0xd+uVp3/XKNvnJ9le65dpLMzFNCAACQSJwRBQCkpNHDCvTvX/iw5l1Wru//bqu+9PBahdrCvmMBAIAE4IwoACBlFeQG9X8+dZmmjBmqv/vtFu08dEJLP1PT4w2OAABAeuCMKAAgpZmZ7q69WP/42Rq9d/Ckbv3RS1rz3mHfsQAAQD9QRAEAaeG6KWV67J6rVJQf1B1LX9Wv3tjlOxIAALhAFFEAQNqYNGqIHr/nGn1oYom++ug6fefJTeoIR3zHAgAAfcQ1ogCAtDJsUK7+3+cu1/ee2aKfv/Sutu0/rhunj9aDf3hHTc0hlRcXauHc6nPuwgsAAFIHRRQAkHZyggF980+nacroIVq0fJ3+uO2gXPS9xuaQFi9fL0mUUQAAUhRbcwEAaes/XD5WIwbnnSmhp4Xaw1qyssFLJgAA0DuKKAAgrR1qaYs53tQcSnISAAAQL4ooACCt9fRMUZ41CgBA6qKIAgDS2sK51SrMDZ41lp8T0MK51Z4SAQCA3nCzIgBAWjt9Q6IlKxvU1BySmTRiUK5unjnGczIAANATiigAIO3Nm1VxppA+u2mfPv+L1frpizv0l3WTPCcDAACxsDUXAJBR/mRamW6cPlo//P027Tx4wnccAAAQA0UUAJBxvnXrJcoLBnTvivVyrvvDXQAAgG8UUQBAxikbWqCv3ThFL20/pOVvNvqOAwAAuqGIAgAy0qevGKc5Fw3X3z69SYdaTvmOAwAAuqCIAgAyUiBgun/+DLWc6tB9T2/2HQcAAHRBEQUAZKyqsiH6i9qLtXxto17cdsB3HAAAEEURBQBktL+8dpImjhysex/boFBb2HccAAAgiigAIMMV5Ab1vfkz9P7hk/rhc9t8xwEAAKKIAgCywIcmluhTNWP10xd3aFPTMd9xAADIehRRAEBWWHzTFA0flKvFy9cpHOHZogAA+EQRBQBkheJBefrmn16it3cf1S9e2ek7DgAAWY0iCgDIGn86c4zqqku1ZGWDGptDvuMAAJC1KKIAgKxhZvrubdPlnPTNFRvkHFt0AQDwgSIKAMgqY0cM0l9fX6XntuzXbzbs9R0HAICsRBEFAGSdz101XtMrhupvntioo6F233EAAMg6FFEAQNbJCQb0wPyZOtRySn/32y2+4wAAkHUoogCArDS9Ypj+/JoJ+rfX3tcbOw/7jgMAQFbJ8R0AAABf/sefVOmZ9Xt1z7++qZyAac/RVpUXF2rh3GrNm1XhOx4AABmLM6IAgKw1KC9HN88crf3HT6npaKucpMbmkBYvX68Vaxt9xwMAIGNRRAEAWe3pdefeOTfUHtaSlQ0e0gAAkB0oogCArNbUHOrTOAAA6D+KKAAgq5UXF/ZpHAAA9B9FFACQ1RbOrVZhbvCsscLcoBbOrfaUCACAzMddcwEAWe303XGXrNyixuZW5QUDun/+DO6aCwDAAOKMKAAg682bVaGXFn1MX7m+Sm3hiGaPG+47EgAAGY0iCgBA1PzZlTKTfv3mbt9RAADIaHEVUTO7wcwazGy7mS2K8f4wM3vSzN42s41m9meJjwoAwMAqLy7UNZNG6tE1uxWJON9xAADIWL0WUTMLSnpQ0o2Spkm608ymdTvsHkmbnHOXSqqT9L/NLC/BWQEAGHC314xVY3NIr+w45DsKAAAZK54zoldI2u6c2+Gca5O0TNJt3Y5xkoaYmUkqknRYUkdCkwIAkATXTyvTkIIc/XoN23MBABgo8RTRCkm7urzeHR3r6keSpkpqkrRe0pedc5GEJAQAIIkKcoO69dJy/WbDHh1rbfcdBwCAjBTP41ssxlj3C2fmSnpL0nWSLpb0rJm96Jw7dtYXMrtL0l2SVFZWpvr6+r7mTaqWlpaUz4jMxNqDL6y9ThdbWK3tEf39I/WqG5vrO05WYO3BF9YefMrm9RdPEd0taWyX15XqPPPZ1Z9JesA55yRtN7N3JU2R9HrXg5xzSyUtlaSamhpXV1d3gbGTo76+XqmeEZmJtQdfWHudap3Twzte0PqWXH2r7irfcbICaw++sPbgUzavv3i25r4habKZTYjegOgOSU90O+Z9SR+TJDMrk1QtaUcigwIAkCxmpgVzKrXmvSN650CL7zgAAGScXouoc65D0hclrZS0WdKvnHMbzexuM7s7eth3JV1lZuslPSfpa865gwMVGgCAgfaJ2RUKBoybFgEAMADi2Zor59wzkp7pNvZQl8+bJF2f2GgAAPgzakiB6qpKtfzN3frK9dUKBmLdMgEAAFyIeLbmAgCQlRbMqdS+Y6f04rYDvqMAAJBRKKIAAPTgY1PLNHxQrh5hey4AAAlFEQUAoAd5OQHddlmFnt24T80n23zHAQAgY1BEAQA4j9trKtUWjuiJt7s/uQwAAFwoiigAAOdxSfkwTR0zlLvnAgCQQBRRAAB6cfucSq3bfVQNe4/7jgIAQEagiAIA0It5syqUGzQ9snqX7ygAAGQEiigAAL0YMThPH5tSphVvNao9HPEdBwCAtEcRBQAgDgvmVOpgS5vqG3imKAAA/UURBQAgDnXVpRpZlM/2XAAAEoAiCgBAHHKCAc2fXaHnt+zXwZZTvuMAAJDWKKIAAMRpwZxKdUScHn+LZ4oCANAfFFEAAOJUVTZEl1YO0yOrd8k55zsOAABpiyIKAEAfLKgZqy17j2tj0zHfUQAASFsUUQAA+uDWmeXKywlw0yIAAPqBIgoAQB8MG5Sr66eV6fG3m3SqI+w7DgAAaYkiCgBAH91eM1bNJ9v13Ob9vqMAAJCWKKIAAPTRNZNGasywArbnAgBwgSiiAAD0UTBgmj+7Qqu2HtC+Y62+4wAAkHYoogAAXIBPzq5UxEmPrW30HQUAgLRDEQUA4AJMLC1SzUXDeaYoAAAXgCIKAMAFur2mUu8cOKG1u5p9RwEAIK1QRAEAuEA3zRijgtyAfr1mt+8oAACkFYooAAAXaEhBrm6aPkZPvt2k1naeKQoAQLxyfAcAACCdLaip1PK1jfrw/c+p+WS7yosLtXButebNqvAdDQCAlEURBQCgH/Y1t8okHTnZLklqbA5p8fL1kkQZBQCgB2zNBQCgH77/7FZ1v2duqD2sJSsbvOQBACAdUEQBAOiHpuZQn8YBAABFFACAfikvLuzTOAAAoIgCANAvC+dWqzA3eNZYYW5QC+dWe0oEAEDq42ZFAAD0w+kbEv3t05t0sKVNJYPz9I1bpnGjIgAAzoMzogAA9NO8WRV6adF1GpQX1E0zxlBCAQDoBUUUAIAEyM8J6qqLS1S/db+c634fXQAA0BVFFACABKmtKtWuwyHtPHTSdxQAAFIaRRQAgASprRolSapv2O85CQAAqY0iCgBAgowrGaSJIwdr1dYDvqMAAJDSKKIAACTQR6tK9eqOQ2ptD/uOAgBAyqKIAgCQQHXVpWptj+i1dw/7jgIAQMqiiAIAkEAfmlii/JyAVjWwPRcAgJ5QRAEASKCC3KCunFiiVVu5YREAAD2hiAIAkGB1VaV658AJ7TrMY1wAAIiFIgoAQILVVpdKEnfPBQCgBxRRAAASbOLIwaocXkgRBQCgBxRRAAASzMxUV12ql7cfVFtHxHccAABSDkUUAIABUFs1Sifawlr9Ho9xAQCgO4ooAAAD4MMXlyg3aGzPBQAgBoooAAADoCg/RzUXjeB5ogAAxEARBQBggNRVl2rL3uPae7TVdxQAAFIKRRQAgAFy+jEuL7A9FwCAs1BEAQAYINVlQzR6aIHqt+73HQUAgJRCEQUAYICYmWqrSvXitoPqCPMYFwAATouriJrZDWbWYGbbzWxRD8fUmdlbZrbRzFYlNiYAAOmptrpUx1s79NauZt9RAABIGb0WUTMLSnpQ0o2Spkm608ymdTumWNKPJd3qnLtE0u2JjwoAQPq5etJIBQOmeu6eCwDAGfGcEb1C0nbn3A7nXJukZZJu63bMpyUtd869L0nOOS6GAQBA0rDCXM0eV8zzRAEA6MKcc+c/wGyBpBucc/81+vozkq50zn2xyzE/kJQr6RJJQyT90Dn3ixhf6y5Jd0lSWVnZnGXLliXoP2NgtLS0qKioyHcMZCHWHnxh7Q2MJ95p0/Jt7fq/1w7S0HzzHSclsfbgC2sPPmX6+rv22mvXOOdqYr2XE8f8WP9idm+vOZLmSPqYpEJJr5jZq865rWdNcm6ppKWSVFNT4+rq6uL49v7U19cr1TMiM7H24Atrb2CUTDqq5dv+qI7SyaqbXek7Tkpi7cEX1h58yub1F8/W3N2SxnZ5XSmpKcYxv3XOnXDOHZT0gqRLExMRAID0dkn5UI0symN7LgAAUfEU0TckTTazCWaWJ+kOSU90O+ZxSR8xsxwzGyTpSkmbExsVAID0FAiYPjq5VC9sPaBw5PyXxAAAkA16LaLOuQ5JX5S0Up3l8lfOuY1mdreZ3R09ZrOk30paJ+l1ST9zzm0YuNgAAKSX2upSHTnZrg2NR31HAQDAu3iuEZVz7hlJz3Qbe6jb6yWSliQuGgAAmeMjk0tlJtU3HNClY4t9xwEAwKt4tuYCAIB+GjE4TzMri7VqK084AwCAIgoAQJLUVpXqrV3Naj7Z5jsKAABeUUQBAEiSuupSRZz04raDvqMAAOAVRRQAgCS5tLJYwwpzeYwLACDrUUQBAEiSYMD0kckjtWrrATnHY1wAANmLIgoAQBLVVY/SgeOntGnPMd9RAADwhiIKAEASfXTySEliey4AIKtRRAEASKJRQws0bcxQrWqgiAIAshdFFACAJKurLtWa947oeGu77ygAAHhBEQUAIMlqq0rVEXF6afsh31EAAPCCIgoAQJLNvmi4huTncJ0oACBrUUQBAEiy3GBAV08aqVUN+3mMCwAgK1FEAQDwoLa6VE1HW7V9f4vvKAAAJB1FFAAAD2qrSiXxGBcAQHaiiAIA4EF5caGqyopUz2NcAABZiCIKAIAntVWlev3dwzrZ1uE7CgAASUURBQDAk9qqUWoLR/TqDh7jAgDILhRRAAA8uXzCcBXmBtmeCwDIOhRRAAA8yc8J6qqLS7hhEQAg61BEAQDwqLa6VO8dOqmdB0/4jgIAQNJQRAEA8KiuapQkqb5hv+ckAAAkD0UUAACPxpUM0siiPN3/my2asOhpXf3A81qxttF3LAAABlSO7wAAAGSzFWsbdeRku8IRJ0lqbA5p8fL1kqR5syp8RgMAYMBwRhQAAI+WrGw4U0JPC7WHtWRlg6dEAAAMPIooAAAeNTWH+jQOAEAmoIgCAOBReXFhn8YBAMgEFFEAADxaOLdahbnBs8YKc4NaOLfaUyIAAAYeNysCAMCj0zck+t4zm7X/+CkVF+bqW7dewo2KAAAZjTOiAAB4Nm9WhV77nx9TRXGhrpgwghIKAMh4FFEAAFKAmam2ulQvbT+oto6I7zgAAAwoiigAACmirqpUJ9rCWvPeEd9RAAAYUBRRAABSxFWTRionYFq19YDvKAAADCiKKAAAKaIoP0c144ervmG/7ygAAAwoiigAACmkrnqUtuw9rn3HWn1HAQBgwFBEAQBIIbVVpZLE9lwAQEajiAIAkEKmjB6isqH5WtVAEQUAZC6KKAAAKcTMVFtVqhe3HVBHmMe4AAAyE0UUAIAUU1s1SsdaO/T27mbfUQAAGBAUUQAAUsw1k0YqYGJ7LgAgY1FEAQBIMcMG5Wr2uOGq54ZFAIAMRREFACAF1VaVat3uozrYcsp3FAAAEo4iCgBACqqt7nyMyx+3HfScBACAxKOIAgCQgqaXD1PJ4DzVN+z3HQUAgISjiAIAkIICAdNHq0r1wraDikSc7zgAACQURRQAgBRVW1WqwyfatKHpqO8oAAAkFEUUAIAU9ZHJI2Um1fMYFwBAhqGIAgCQokqK8jWzYphW8RgXAECGoYgCAJDCaqtKtfb9Izp6st13FAAAEoYiCgBACqutHqWIk17czllRAEDmoIgCAJDCLq0cpmGFuVrFdaIAgAwSVxE1sxvMrMHMtpvZovMcd7mZhc1sQeIiAgCQvXKCAV0zeaRWbT0g53iMCwAgM/RaRM0sKOlBSTdKmibpTjOb1sNxfydpZaJDAgCQzeqqSrX/+Clt3nPcdxQAABIinjOiV0ja7pzb4Zxrk7RM0m0xjvuSpEcl7U9gPgAAsl5tVakkcfdcAEDGyInjmApJu7q83i3pyq4HmFmFpE9Iuk7S5T19ITO7S9JdklRWVqb6+vo+xk2ulpaWlM+IzMTagy+svdQ1dkhAj7++VVPP+ic5c7D24AtrDz5l8/qLp4hajLHuF6n8QNLXnHNhs1iHRyc5t1TSUkmqqalxdXV18aX0pL6+XqmeEZmJtQdfWHup65bWLfrpCzs050NXa0hBru84Ccfagy+sPfiUzesvnq25uyWN7fK6UlJTt2NqJC0zs52SFkj6sZnNS0RAAADQuT23I+L08juHfEcBAKDf4imib0iabGYTzCxP0h2Snuh6gHNugnNuvHNuvKRfS/pL59yKRIcFACBbzR43XEX5OVwnCgDICL1uzXXOdZjZF9V5N9ygpJ875zaa2d3R9x8a4IwAAGS9vJyArrq4RKsaOh/jcr5LYQAASHXxXCMq59wzkp7pNhazgDrnPtf/WAAAoLu66lH63aZ9eudAiyaNGuI7DgAAFyyerbkAACAFfLRqpCSpvoHtuQCA9EYRBQAgTVQOH6RJo4q4ThQAkPYoogAApJG6qlK9tuOwTrZ1+I4CAMAFo4gCAJBGaqtL1RaO6LUdh31HAQDgglFEAQBII5ePH6HC3KDqG/b7jgIAwAWjiAIAkEYKcoP68MUlXCcKAEhrFFEAANJMbVWpdh46qZ0HT/iOAgDABaGIAgCQZmqrSiWJs6IAgLRFEQUAIM2MHzlY40sGUUQBAGmLIgoAQBqqrSrVK+8cUmt72HcUAAD6jCIKAEAaqq0uVag9rDd28hgXAED6oYgCAJCGPjSxRHk5Aa1qYHsuACD9UEQBAEhDg/JydOWEEVwnCgBISxRRAADSVG1Vqbbtb1Fjc8h3FAAA+oQiCgBAmqqrjj7Ghe25AIA0QxEFACBNXVxapIriQq3aut93FAAA+oQiCgBAmjIzfbSqVC9tP6S2jojvOAAAxI0iCgBAGqurLlXLqQ69+f4R31EAAIgbRRQAgDR25MQpSdIdS1/V1Q88rxVrGz0nAgCgdxRRAADS1Iq1jfr2k5vPvG5sDmnx8vWUUQBAyqOIAgCQppasbFCoPXzWWKg9rCUrGzwlAgAgPhRRAADSVFMPzw/taRwAgFRBEQUAIE2VFxf2aRwAgFRBEQUAIE0tnFutwtzgWWP5OQEtnFvtKREAAPHJ8R0AAABcmHmzKiR1Xiva1BySk3T5+OFnxgEASFUUUQAA0ti8WRVniueXHl6rF7Ye0KmOsPJzgr3MBADAH7bmAgCQIW6fU6mjoXb9ftN+31EAADgviigAABni6kkjNWZYgR5Zs8t3FAAAzosiCgBAhggGTPNnV+iFrQe092ir7zgAAPSIIgoAQAZZMGesIk56bG2j7ygAAPSIIgoAQAaZMHKwLh8/XI+s2SXnnO84AADERBEFACDDLJhTqR0HTujN95t9RwEAICaKKAAAGebmmeUqzA3q19y0CACQoiiiAABkmKL8HN04Y7SeenuPQm1h33EAADgHRRQAgAx0+5yxOn6qQys37vUdBQCAc1BEAQDIQFdOGKHK4YU8UxQAkJIoogAAZKBAwLRgTqVefueQdh856TsOAABnoYgCAJChPjm7Us5Jy9/kmaIAgNRCEQUAIEONHTFIH55Yol+v2a1IhGeKAgBSB0UUAIAMdntNpd4/fFKv7zzsOwoAAGdQRAEAyGA3Th+jovwcPbJ6t+8oAACcQREFACCDFeYFdcvMMfrNhj06carDdxwAACRRRAEAyHgL5lTqZFtYT6/f4zsKAACSKKIAAGS8ORcN18SRg/VrtucCAFIERRQAgAxnZvrknEq9vvOwdh484TsOAAAUUQAAssEnZ1cqYNKjb3JWFADgH0UUAIAsMHpYga6ZXKpH1+xWmGeKAgA8o4gCAJAlbp9TqaajrXr5nYO+owAAshxFFACALPEn08o0tIBnigIA/KOIAgCQJQpyg7rtsgqt3LhXR0PtvuMAALJYXEXUzG4wswYz225mi2K8/x/NbF30z8tmdmniowIAgP5aMKdSpzoiempdk+8oAIAs1msRNbOgpAcl3ShpmqQ7zWxat8PelVTrnJsp6buSliY6KAAA6L+ZlcNUVVbE9lwAgFfxnBG9QtJ259wO51ybpGWSbut6gHPuZefckejLVyVVJjYmAABIBDPT7XPG6q1dzdq+/7jvOACALBVPEa2QtKvL693RsZ78uaTf9CcUAAAYOPNmVSgYMD2yhrOiAAA/cuI4xmKMxXwAmZldq84iek0P798l6S5JKisrU319fXwpPWlpaUn5jMhMrD34wtrLHjNKAlr26ru6In+vgoFY/9QnF2sPvrD24FM2r794iuhuSWO7vK6UdM4dDsxspqSfSbrROXco1hdyzi1V9PrRmpoaV1dX19e8SVVfX69Uz4jMxNqDL6y97NE6cq/u/pc1svJpqptS5jsOaw/esPbgUzavv3i25r4habKZTTCzPEl3SHqi6wFmNk7Sckmfcc5tTXxMAACQSNdNGaURg/O4aREAwItei6hzrkPSFyWtlLRZ0q+ccxvN7G4zuzt62DcllUj6sZm9ZWarBywxAADot7ycgOZdVqHfb96nIyfafMcBAGSZuJ4j6px7xjlX5Zy72Dl3X3TsIefcQ9HP/6tzbrhz7rLon5qBDA0AAPpvwZxKtYedHn+r0XcUAECWiauIAgCAzDOtfKguKR/K3XMBAElHEQUAIItVlxVpY9MxTVj0tK5+4HmtWMvZUQDAwKOIAgCQpVasbdQzG/ZK6nwuW2NzSIuXr6eMAgAGHEUUAIAstWRlg1rbI2eNhdrDWrKywVMiAEC2oIgCAJClmppDfRoHACBRKKIAAGSp8uLCPo0DAJAoFFEAALLUwrnVKswNnjUWMOkr11d5SgQAyBYUUQAAstS8WRW6f/4MVRQXyiQVF+Yq4qT2sPMdDQCQ4XJ8BwAAAP7Mm1WhebMqJEmRiNMdP31V9z2zWddOGaXSIfme0wEAMhVnRAEAgCQpEDB97xMzFGoL67tPbfIdBwCQwSiiAADgjEmjinTPtZP0xNtN+kPDft9xAAAZiiIKAADOcnfdRE0aVaSvP7ZBJ9s6fMcBAGQgiigAADhLfk5Q98+focbmkP7Ps1t9xwEAZCCKKAAAOMfl40fo01eO0z/+8V1taDzqOw4AIMNQRAEAQExfu2GKSorytWj5OnWEI77jAAAyCEUUAADENKwwV9++9RJtaDymf3p5p+84AIAMQhEFAAA9unH6aH186ij9799t1a7DJ33HAQBkCIooAADokZnpO7dNV8Ckr6/YIOec70gAgAxAEQUAAOdVXlyor8yt1qqtB/Tkuj2+4wAAMgBFFAAA9Oo/f3i8Lq0cpu88uVHNJ9t8xwEApDmKKAAA6FUwYLp//kwdOdmu+5/Z4jsOACDNUUQBAEBcppUP1ec/MlH/vnqXXnnnkO84AIA0RhEFAABx+/LHJmvciEG697H1am0P+44DAEhTFFEAABC3wryg7vvEdO04eEI//sN233EAAGmKIgoAAPrkI5NLNX9WhX6y6h1t3XfcdxwAQBqiiAIAgD679+apKsrP0eLl6xWJ8GxRAEDf5PgOAAAA0k9JUb6+fvM0/fUjb2vWd5/VsVC7yosLtXButebNqvAdDwCQ4iiiAADgggSs88/RULskqbE5pMXL10sSZRQAcF5szQUAABfk+7/bqu67ckPtYS1Z2eAnEAAgbVBEAQDABWlqDvVpHACA0yiiAADggpQXF8YcHzU0P8lJAADphiIKAAAuyMK51SrMDZ4z3tLaodffPewhEQAgXVBEAQDABZk3q0L3z5+hiuJCmaSK4kItvmmKyoYW6NM/fVUPv/6+74gAgBTFXXMBAMAFmzer4pw75N5RM05fWrZWi5ev15Y9x/T1W6YpN8jvvgEAH+BfBQAAkFDDBuXq55+t0ec/MkH//Mp7+uzPX9eRE22+YwEAUghFFAAAJFxOMKB7b56m799+qVbvPKLbHnxJW/cd9x0LAJAiKKIAAGDALJhTqWVf+JBC7WF94sGX9PtN+3xHAgCkAIooAAAYULPHDdcTX7xaE0uL9PlfrtaDf9gu55zvWAAAjyiiAABgwI0ZVqhH7v6wbr20XEtWNui/LXtLobaw71gAAE+4ay4AAEiKgtygfvCpyzRl9FD9r5Vb9O7BFn1ydqV+9uK7amwOqeLV57VwbvU5d+EFAGQeiigAAEgaM9Nf1F2sqrIi/eW/rNG3n9x05r3G5pAWL18vSZRRAMhwbM0FAABJ97GpZRo2KO+c8VB7WEtWNnhIBABIJoooAADw4sDxUzHHm5pDSU4CAEg2iigAAPCivLgw5riT9B/+4RUtf3O3Wtu5oREAZCKKKAAA8GLh3GoV5gbPGivICeiWmWO0/1ir/upXb+uK+36vv3l8gzbvOeYpJQBgIHCzIgAA4MXpGxItWdnQedfc4sIzd811zumVHYe07PVdevj1XfrnV97TZWOLdecVY3XLzHINzudHGABIZ/wtDgAAvJk3q0LzZlWovr5edXV1Z8bNTFddPFJXXTxSR060afnaRi17/X197dH1+s6Tm3TrZRW684qxemd/i77/u61qag6pvEuRBQCkNoooAABIacMH5+nPr5mg/3L1eL35/hE9/PouPbZ2tx5+/X2ZOq8plXj8CwCkE64RBQAAacHMNOeiEfr+7Zfq9Xs/rmGFuWdK6Gmh9rC+sWKDfrN+j3YePKFIpPsRAIBUwBlRAACQdoYW5OpYqD3me8dPdegv/vVNSdLgvKCqRw/RlDFDNXXMUE0bM0TVo4eqKHqN6Yq1jVqysoGtvQCQZBRRAACQlsqLC9UY45mj5cMK9NBn5mjznmPavOe4Nu05piffbtK/vfb+mWMuKhmkoQU52rznuDqiZ007t/aukxT/1t7+FlmKMIBsFVcRNbMbJP1QUlDSz5xzD3R736Lv3yTppKTPOefeTHBWAACAMxbOrdbi5esV6vKs0cLcoL56wxTNrCzWzMriM+POOTU2h7Rlz/HOgrr3mFZu3Kdwt627ofaI/upXb+nBP2xXSVGeSgbna8TgvOjneSopir4enKfX3j2s+57epFB7RFLfr1FdsbbxrPwXco2r7yLsc36ivndjc0gVrz6fVv/tzM+c+Rey/jLlF1i9FlEzC0p6UNKfSNot6Q0ze8I5t6nLYTdKmhz9c6Wkn0Q/AgAADIiuj3/p7QcyM1Pl8EGqHD5IH59WJkmasOjpmF834qSJpYN1+ESbNu89pkMtbTrawzbg7kLtYX3t0XV68u0mFeQFVZATVGFeQIW5QRVE/5z+fMnKLWeV6NPz73t6syaNKlJuMKBgwJQTMOUETTmBzte5QVMwYPrN+r365hMb1OqpCPucn87Zmc98n2s/lZhz57+I38w+LOlbzrm50deLJck5d3+XY/5BUr1z7uHo6wZJdc65PT193ZqaGrd69er+/xcMoO63kgeShbUHX1h78MXH2rv6gedjbu2tKC7US4uuO2usPRzRkRNtOnSiTYdPtOlgyyl9edlbPX7tS8qHKtQe1qn2iELtYYXawmrtCKuXH7sSpiA3oICZAmYySWZSIGDRsc5ifqjllGLdyykYMFUUF8pM0bmdX0PR14qO7Tx44sy25q5yAqZJo4rUuWEuenyX908Pb913XO3hc+fnBk1TRg+V2Tlvnflam/Yc63HutPJhsSd2sanpaI/zL4lj/sbzzJ9e0fv8DY3MZ/6Fze9pbqy/t1KBma1xztXEei+erbkVknZ1eb1b557tjHVMhaSziqiZ3SXpLkkqKytTfX19HN/en5aWlpTPiMzE2oMvrD344mPt3TwurH86JrVFPhjLC3SO95ZlmKSSAtOh1nN/ICwpMC2cefpMZyD6J1fOObVHpLaw1B5x+vYrrWo+de78IbnS56bnK+KksJPCEffB506KRDo/Lmto6zHftZVBOefknBSR5JzkTv/Pdb6uPx57bjjiVJF/6swdiTvnfsBFx2KVUKlzvDByMuZ7p4u4k2L+MK3ouLW1xA7nep8bCfXwH9btuJ7GO/o5v/0E85k/cPN7mtvYHEq7f7/jKaKxfh/V/f+BeI6Rc26ppKVS5xnRVP+tO2cG4AtrD76w9uCLj7VXJ2laP661+saws7fISZ3XqH7jthmqi+NrWFns+d+dPyOuDC+e54zuT77Q+5mR850RfvjL/Zv/2F/1b/4Tf33++eeb++RX+ve9n+rv/IXMZ/7AzT/f3HT79zue54juljS2y+tKSU0XcAwAAEBKmTerQi8tuk7vPnCzXlp0XZ+usZo3q0L3z5/RuY1VnT8I3h9niUzE/IVzq1WYGzxrrDA3qIVzqzN+fjpnZz7zfa79VBLPGdE3JE02swmSGiXdIenT3Y55QtIXzWyZOrftHj3f9aEAAACZYN6sin7dIKQ/8/tys6ZMm5/I793YHFJFGv23Mz+z5vd1/fX3e6eSXm9WJElmdpOkH6jz8S0/d87dZ2Z3S5Jz7qHo41t+JOkGdT6+5c+cc+e9ExE3KwJ6xtqDL6w9+MLagy+sPfiU6euvvzcrknPuGUnPdBt7qMvnTtI9/QkJAAAAAMgO8VwjCgAAAABAwlBEAQAAAABJRREFAAAAACQVRRQAAAAAkFQUUQAAAABAUlFEAQAAAABJRREFAAAAACQVRRQAAAAAkFQUUQAAAABAUplzzs83Njsg6T0v3zx+IyUd9B0CWYm1B19Ye/CFtQdfWHvwKdPX30XOudJYb3grounAzFY752p850D2Ye3BF9YefGHtwRfWHnzK5vXH1lwAAAAAQFJRRAEAAAAASUURPb+lvgMga7H24AtrD76w9uALaw8+Ze364xpRAAAAAEBScUYUAAAAAJBUFNEYzOwGM2sws+1mtsh3HmQ2M/u5me03sw1dxkaY2bNmti36cbjPjMg8ZjbWzP5gZpvNbKOZfTk6ztrDgDOzAjN73czejq6/b0fHWX8YcGYWNLO1ZvZU9DXrDklhZjvNbL2ZvWVmq6NjWbv+KKLdmFlQ0oOSbpQ0TdKdZjbNbypkuH+SdEO3sUWSnnPOTZb0XPQ1kEgdkv7aOTdV0ock3RP9u461h2Q4Jek659ylki6TdIOZfUisPyTHlyVt7vKadYdkutY5d1mXR7Zk7fqjiJ7rCknbnXM7nHNtkpZJus1zJmQw59wLkg53G75N0j9HP/9nSfOSmQmZzzm3xzn3ZvTz4+r8oaxCrD0kgevUEn2ZG/3jxPrDADOzSkk3S/pZl2HWHXzK2vVHET1XhaRdXV7vjo4ByVTmnNsjdRYGSaM850EGM7PxkmZJek2sPSRJdHvkW5L2S3rWOcf6QzL8QNJXJUW6jLHukCxO0u/MbI2Z3RUdy9r1l+M7QAqyGGPcWhhARjKzIkmPSvrvzrljZrH+CgQSzzkXlnSZmRVLeszMpnuOhAxnZrdI2u+cW2NmdZ7jIDtd7ZxrMrNRkp41sy2+A/nEGdFz7ZY0tsvrSklNnrIge+0zszGSFP2433MeZCAzy1VnCf1X59zy6DBrD0nlnGuWVK/Oa+VZfxhIV0u61cx2qvPSq+vM7F/EukOSOOeaoh/3S3pMnZcEZu36o4ie6w1Jk81sgpnlSbpD0hOeMyH7PCHps9HPPyvpcY9ZkIGs89TnP0ra7Jz7+y5vsfYw4MysNHomVGZWKOnjkraI9YcB5Jxb7JyrdM6NV+fPd8875/6TWHdIAjMbbGZDTn8u6XpJG5TF68+cY9dpd2Z2kzqvIQhK+rlz7j6/iZDJzOxhSXWSRkraJ+lvJK2Q9CtJ4yS9L+l251z3GxoBF8zMrpH0oqT1+uBaqf+pzutEWXsYUGY2U5035Qiq85fiv3LOfcfMSsT6QxJEt+Z+xTl3C+sOyWBmE9V5FlTqvDzy35xz92Xz+qOIAgAAAACSiq25AAAAAICkoogCAAAAAJKKIgoAAAAASCqKKAAAAAAgqSiiAAAAAICkoogCAAAAAJKKIgoAAAAASCqKKAAAAAAgqf4/0IJ9l7oxZvUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot initial interface profile\n",
    "x = np.arange(0, domain_size[1]+2, 1)   # start,stop,step\n",
    "y = dh.cpu_arrays[\"C\"][int(domain_size[0] * 0.5), :]\n",
    "\n",
    "plt.plot(x,y, marker='o')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "force_h = interface_tracking_force(C, stencil_phase, parameters)\n",
    "hydro_force = hydrodynamic_force(g, C, method_hydro, parameters, body_force)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definition of the LB update rules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[['spatialInner0'], ['spatialInner1']]\n"
     ]
    }
   ],
   "source": [
    "lbm_optimisation = LBMOptimisation(symbolic_field=h, symbolic_temporary_field=h_tmp)\n",
    "allen_cahn_update_rule = create_lb_update_rule(lbm_config=config_phase,\n",
    "                                               lbm_optimisation=lbm_optimisation)\n",
    "\n",
    "allen_cahn_update_rule = add_interface_tracking_force(allen_cahn_update_rule, force_h)\n",
    "\n",
    "ast_kernel = ps.create_kernel(allen_cahn_update_rule, target=dh.default_target, cpu_openmp=True)\n",
    "kernel_allen_cahn_lb = ast_kernel.compile()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[['spatialInner0'], ['spatialInner1']]\n"
     ]
    }
   ],
   "source": [
    "force_Assignments = hydrodynamic_force_assignments(g, u, C, method_hydro, parameters, body_force)\n",
    "\n",
    "lbm_optimisation = LBMOptimisation(symbolic_field=g, symbolic_temporary_field=g_tmp)\n",
    "hydro_lb_update_rule = create_lb_update_rule(lbm_config=config_hydro,\n",
    "                                             lbm_optimisation=lbm_optimisation)\n",
    "\n",
    "hydro_lb_update_rule = add_hydrodynamic_force(hydro_lb_update_rule, force_Assignments, C, g, parameters)\n",
    "\n",
    "ast_kernel = ps.create_kernel(hydro_lb_update_rule, target=dh.default_target, cpu_openmp=True)\n",
    "kernel_hydro_lb = ast_kernel.compile()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boundary Conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# periodic Boundarys for g, h and C\n",
    "periodic_BC_C = dh.synchronization_function(C.name, target=dh.default_target, optimization = {\"openmp\": True})\n",
    "\n",
    "periodic_BC_g = LBMPeriodicityHandling(stencil=stencil_hydro, data_handling=dh, pdf_field_name=g.name,\n",
    "                                       streaming_pattern='push')\n",
    "periodic_BC_h = LBMPeriodicityHandling(stencil=stencil_phase, data_handling=dh, pdf_field_name=h.name,\n",
    "                                       streaming_pattern='pull')\n",
    "\n",
    "# No slip boundary for the phasefield lbm\n",
    "bh_allen_cahn = LatticeBoltzmannBoundaryHandling(method_phase, dh, 'h',\n",
    "                                                 target=dh.default_target, name='boundary_handling_h',\n",
    "                                                 streaming_pattern='pull')\n",
    "\n",
    "# No slip boundary for the velocityfield lbm\n",
    "bh_hydro = LatticeBoltzmannBoundaryHandling(method_hydro, dh, 'g' ,\n",
    "                                            target=dh.default_target, name='boundary_handling_g',\n",
    "                                            streaming_pattern='push')\n",
    "\n",
    "contact_angle = BoundaryHandling(dh, C.name, stencil_hydro, target=dh.default_target)\n",
    "contact = ContactAngle(90, interface_width)\n",
    "\n",
    "wall = NoSlip()\n",
    "if dimensions == 2:\n",
    "    bh_allen_cahn.set_boundary(wall, make_slice[:, 0])\n",
    "    bh_allen_cahn.set_boundary(wall, make_slice[:, -1])\n",
    "\n",
    "    bh_hydro.set_boundary(wall, make_slice[:, 0])\n",
    "    bh_hydro.set_boundary(wall, make_slice[:, -1])\n",
    "    \n",
    "    contact_angle.set_boundary(contact, make_slice[:, 0])\n",
    "    contact_angle.set_boundary(contact, make_slice[:, -1])\n",
    "else:\n",
    "    bh_allen_cahn.set_boundary(wall, make_slice[:, 0, :])\n",
    "    bh_allen_cahn.set_boundary(wall, make_slice[:, -1, :])\n",
    "\n",
    "    bh_hydro.set_boundary(wall, make_slice[:, 0, :])\n",
    "    bh_hydro.set_boundary(wall, make_slice[:, -1, :])\n",
    "    \n",
    "    contact_angle.set_boundary(contact, make_slice[:, 0, :])\n",
    "    contact_angle.set_boundary(contact, make_slice[:, -1, :])\n",
    "\n",
    "\n",
    "bh_allen_cahn.prepare()\n",
    "bh_hydro.prepare()\n",
    "contact_angle.prepare()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Full timestep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# definition of the timestep for the immiscible fluids model\n",
    "def timeloop():\n",
    "    # Solve the interface tracking LB step with boundary conditions\n",
    "    periodic_BC_h()\n",
    "    bh_allen_cahn()    \n",
    "    dh.run_kernel(kernel_allen_cahn_lb, **parameters.symbolic_to_numeric_map)\n",
    "    dh.swap(\"C\", \"C_tmp\")\n",
    "    \n",
    "    # apply the three phase-phase contact angle\n",
    "    contact_angle()\n",
    "    # periodic BC of the phase-field\n",
    "    periodic_BC_C()\n",
    "    \n",
    "    # solve the hydro LB step with boundary conditions\n",
    "    dh.run_kernel(kernel_hydro_lb, **parameters.symbolic_to_numeric_map)\n",
    "    periodic_BC_g()\n",
    "    bh_hydro()\n",
    "    \n",
    "    # compute density (only for vtk output)\n",
    "    # must be done BEFORE swapping fields to avoid having outdated values\n",
    "    compute_density()\n",
    "    \n",
    "    # field swaps\n",
    "    dh.swap(\"h\", \"h_tmp\")\n",
    "    dh.swap(\"g\", \"g_tmp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initialize_hydrostatic_pressure():   \n",
    "    for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):       \n",
    "        \n",
    "        # get y as cell center coordinate, i.e., including shift with 0.5\n",
    "        y = np.zeros_like(block.midpoint_arrays[1])\n",
    "        y[:, :] = block.midpoint_arrays[1]\n",
    "        \n",
    "        # compute hydrostatic density\n",
    "        rho_hydrostatic = 1 + 3 * gravitational_acceleration * (y - fluid_depth)\n",
    "        \n",
    "        # subtract 1 because PDFs in incompressible LBM are centered around zero in lbmpy\n",
    "        rho_hydrostatic -= 1\n",
    "\n",
    "        # set equilibrium PDFs with velocity=0 and rho; \n",
    "        for i in range(0, stencil_hydro.Q, 1):\n",
    "            block[\"g\"][:,:,i] = method_hydro.weights[i] * rho_hydrostatic[:,:]\n",
    "            block[\"g_tmp\"][:,:,i] = method_hydro.weights[i] * rho_hydrostatic[:,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_density():   \n",
    "    for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):\n",
    "        # PDFs in incompressible LBM are centered around zero in lbmpy\n",
    "        # => add 1 to whole sum, i.e., initialize with 1\n",
    "        block[\"rho\"].fill(1);\n",
    "        \n",
    "        # compute density\n",
    "        for i in range(block[\"g\"].shape[-1]):\n",
    "            block[\"rho\"][:,:] += block[\"g\"][:,:,i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA54AAAFlCAYAAACDRTcUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUa0lEQVR4nO3dYYxl5Xkf8P+zu7MQG6wUsZAVoK4jbdO4bm2nK+QWqUqDSWmCDF9AjuQItUiokps6VaIUnA+VKlWiSpXaUtMPK9vtVnFjiGMLFKmp6bZWVCmhHhwnrkMSI0Iw9paduLhx1cY2zNMPc9xcYOnM7sy7d+6Z309C55733uE+0oOG+5/n3PdUdwcAAABGObTsAgAAAJg3wRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoY5czje79tpr+8SJE5fzLYED6A+efGbZJQDsO3/hr37vsksADoAnn3zyj7v72KvXL2vwPHHiRNbX1y/nWwIH0G2H7l52CQD7zuPrv7zsEoADoKr+6ELrLrUFAABgKMETAACAoQRPAAAAhrqs3/EEuBzqiiuWXQIAAAt2FDyr6tkk30jycpKXuvtUVV2T5OEkJ5I8m+Se7n5xTJkAAACsqou51PZvdvfbu/vUdP5AkrPdfTLJ2ekcAAAAXmE33/G8M8mZ6fGZJHftuhoAAABmZ6fBs5N8uqqerKr7p7Xru/tckkzH6y70g1V1f1WtV9X6xsbG7isGAABgpex0c6FbuvurVXVdkser6vd2+gbdfTrJ6SQ5depUX0KNAAAArLAdTTy7+6vT8XySTyW5OckLVXU8Sabj+VFFAgAAsLq2nXhW1RuTHOrub0yPfzjJP0nyWJJ7kzw0HR8dWSjATh1yOxUAgH1lJ5faXp/kU1X1ndf/u+7+tar6bJJHquq+JM8luXtcmQAAAKyqbYNndz+T5G0XWP9akltHFAUAAMB87OZ2KgAAALAtwRMAAIChBE8AAACG2ul9PAFWh11tAQD2FRNPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIayqy0wO3XF0WWXAADAAhNPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIayqy0wP0fXll0BAAALTDwBAAAYSvAEAABgKMETAACAoQRPAAAAhhI8AQAAGErwBAAAYCi3UwFmp69wOxUAgP3ExBMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoexqC8xOH/WrDQBgPzHxBAAAYCjBEwAAgKEETwAAAIYSPAEAABhK8AQAAGAoWz8Cs2NXWwCA/cXEEwAAgKEETwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACGcs8BYHY2jx5edgkAACww8QQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKLvaArOzedTf1AAA9hOfzgAAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYase72lbV4STrSb7S3XdU1TVJHk5yIsmzSe7p7hdHFAlwMV62qy0AwL5yMZ/O3p/kqYXzB5Kc7e6TSc5O5wAAAPAKOwqeVXVjkh9N8uGF5TuTnJken0ly155WBgAAwCzsdOL5wSQ/k2RzYe367j6XJNPxur0tDQAAgDnYNnhW1R1Jznf3k5fyBlV1f1WtV9X6xsbGpfwrAAAAWGE7mXjekuTdVfVsko8n+aGq+sUkL1TV8SSZjucv9MPdfbq7T3X3qWPHju1R2QAAAKyKbXe17e4HkzyYJFX1g0l+urvfW1U/l+TeJA9Nx0fHlQmwc5trdrUFANhPdvPp7KEkt1XVl5LcNp0DAADAK+z4Pp5J0t2fSfKZ6fHXkty69yUBAAAwJ65HAwAAYCjBEwAAgKEETwAAAIYSPAEAABjqojYXAlgFm0dr2SUAALDAxBMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoexqC8zOy2t2tQUA2E9MPAEAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYyq62wOxsri27AgAAFpl4AgAAMJTgCQAAwFCCJwAAAEMJngAAAAwleAIAADCU4AkAAMBQbqcCzM7La7XsEgAAWGDiCQAAwFCCJwAAAEMJngAAAAwleAIAADCU4AkAAMBQdrUFZmdzbdkVAACwyMQTAACAoQRPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKHsagvMzuZaLbsEAAAWmHgCAAAwlOAJAADAUIInAAAAQwmeAAAADCV4AgAAMJRdbYHZ2VxbdgUAACwy8QQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAodxOBZgdt1MBANhfTDwBAAAYSvAEAABgKMETAACAoQRPAAAAhhI8AQAAGGrbXW2r6sokv57kiun1n+juf1xV1yR5OMmJJM8muae7XxxXKsDObNqvGwBgX9nJxPObSX6ou9+W5O1Jbq+qdyZ5IMnZ7j6Z5Ox0DgAAAK+wbfDsLf9rOl2b/ukkdyY5M62fSXLXiAIBAABYbTv6jmdVHa6qzyc5n+Tx7n4iyfXdfS5JpuN1r/Oz91fVelWtb2xs7FHZAAAArIodBc/ufrm7357kxiQ3V9Vbd/oG3X26u09196ljx45dYpkAAACsqova1ba7v57kM0luT/JCVR1Pkul4fq+LAwAAYPXtZFfbY0m+3d1fr6rvSvKuJP8syWNJ7k3y0HR8dGShADu1ubbsCgAAWLSTmw4cT3Kmqg5na0L6SHf/alX9RpJHquq+JM8luXtgnQAAAKyobYNnd/9OkndcYP1rSW4dURQAAADzcVHf8QQAAICLJXgCAAAwlOAJAADAUIInAAAAQ+1kV1uAlbK51ssuAQCABSaeAAAADCV4AgAAMJTgCQAAwFCCJwAAAEMJngAAAAxlV1tgdnpt2RUAALDIxBMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAoexqC8zO5pFedgkAACww8QQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAodxOBZidzbVlVwAAwCITTwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACGsqstMDu9trnsEgAAWGDiCQAAwFCCJwAAAEMJngAAAAwleAIAADCU4AkAAMBQdrUF5udIL7sCAAAWmHgCAAAwlOAJAADAUIInAAAAQwmeAAAADCV4AgAAMJRdbYH5WdtcdgUAACww8QQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKMETAACAodxOBZidQ26nAgCwr5h4AgAAMJTgCQAAwFCCJwAAAEMJngAAAAwleAIAADDUtrvaVtVNSf5tku9JspnkdHd/qKquSfJwkhNJnk1yT3e/OK5UgJ05fOTlZZcAAMCCnUw8X0ryU939/UnemeR9VfWWJA8kOdvdJ5Ocnc4BAADgFbYNnt19rrs/Nz3+RpKnktyQ5M4kZ6aXnUly16AaAQAAWGEX9R3PqjqR5B1JnkhyfXefS7bCaZLrXudn7q+q9apa39jY2GW5AAAArJodB8+quirJryT5ye7+k53+XHef7u5T3X3q2LFjl1IjAAAAK2xHwbOq1rIVOj/W3Z+cll+oquPT88eTnB9TIgAAAKtsJ7vaVpKPJHmqu39+4anHktyb5KHp+OiQCgEu0tqaXW0BAPaTbYNnkluS/HiSL1TV56e1D2QrcD5SVfcleS7J3UMqBAAAYKVtGzy7+78kqdd5+ta9LQcAAIC5uahdbQEAAOBiCZ4AAAAMJXgCAAAwlOAJAADAUDvZ1RZgpRw94nYqAAD7iYknAAAAQwmeAAAADCV4AgAAMJTgCQAAwFCCJwAAAEPZ1RaYnSvWXlp2CQAALDDxBAAAYCjBEwAAgKEETwAAAIYSPAEAABhK8AQAAGAou9oCs3PlEbvaAgDsJyaeAAAADCV4AgAAMJTgCQAAwFCCJwAAAEMJngAAAAxlV1tgdq488u1llwAAwAITTwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYyu1UgNl5w5FvLbsEAAAWmHgCAAAwlOAJAADAUIInAAAAQwmeAAAADCV4AgAAMJRdbYHZecORby+7BAAAFph4AgAAMJTgCQAAwFCCJwAAAEMJngAAAAwleAIAADCUXW2B2XnjkW8uuwQAABaYeAIAADCU4AkAAMBQgicAAABDCZ4AAAAMJXgCAAAwlOAJAADAUG6nAszOVYfdTgUAYD8x8QQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgqG13ta2qjya5I8n57n7rtHZNkoeTnEjybJJ7uvvFcWUC7NxVR+xqCwCwn+xk4vlvktz+qrUHkpzt7pNJzk7nAAAA8BrbBs/u/vUk/+NVy3cmOTM9PpPkrr0tCwAAgLm41O94Xt/d55JkOl73ei+sqvurar2q1jc2Ni7x7QAAAFhVwzcX6u7T3X2qu08dO3Zs9NsBAACwz1xq8Hyhqo4nyXQ8v3clAQAAMCfb7mr7Oh5Lcm+Sh6bjo3tWEcAuXXX4T5ddAgAAC7adeFbVLyX5jSTfV1XPV9V92Qqct1XVl5LcNp0DAADAa2w78ezuH3udp27d41oAAACYoeGbCwEAAHCwCZ4AAAAMJXgCAAAw1KXuaguwb119yK62AAD7iYknAAAAQwmeAAAADCV4AgAAMJTgCQAAwFCCJwAAAEMJngAAAAzldirA7Fx92O1UAAD2ExNPAAAAhhI8AQAAGErwBAAAYCjBEwAAgKEETwAAAIayqy0wO1cf+j/LLgEAgAUmngAAAAwleAIAADCU4AkAAMBQgicAAABDCZ4AAAAMZVdbYHbedOhPl10CAAALTDwBAAAYSvAEAABgKMETAACAoQRPAAAAhhI8AQAAGErwBAAAYCi3UwFm52q3UwEA2FdMPAEAABhK8AQAAGAowRMAAIChBE8AAACGEjwBAAAYyq62wOxcfejbyy4BAIAFJp4AAAAMJXgCAAAwlOAJAADAUIInAAAAQwmeAAAADGVXW2B2rq5edgkAACww8QQAAGAowRMAAIChBE8AAACGEjwBAAAYSvAEAABgKLvaArNz9SG/2gAA9hMTTwAAAIYSPAEAABhK8AQAAGAowRMAAIChBE8AAACG2lXwrKrbq+r3q+rpqnpgr4oCAABgPi75ngNVdTjJLyS5LcnzST5bVY919+/uVXEAl+KqQ1cuuwQAABbsZuJ5c5Knu/uZ7v5Wko8nuXNvygIAAGAudhM8b0jy5YXz56c1AAAA+H92EzzrAmv9mhdV3V9V61W1vrGxsYu3AwAAYBXtJng+n+SmhfMbk3z11S/q7tPdfaq7Tx07dmwXbwcAAMAq2k3w/GySk1X15qo6muQ9SR7bm7IAAACYi+p+zdWxO//hqh9J8sEkh5N8tLv/6Tav30jyR5f8hnvv2iR/vOwiuKz0/ODR84NJ3w8ePT+Y9P3g0fP9789392sudd1V8Fx1VbXe3aeWXQeXj54fPHp+MOn7waPnB5O+Hzx6vrp2c6ktAAAAbEvwBAAAYKiDHjxPL7sALjs9P3j0/GDS94NHzw8mfT949HxFHejveAIAADDeQZ94AgAAMNiBDZ5V9dNV1VV17cLag1X1dFX9flX9rWXWx96qqp+rqt+rqt+pqk9V1XcvPKfvM1VVt099fbqqHlh2Pey9qrqpqv5zVT1VVV+sqvdP69dU1eNV9aXp+OeWXSt7q6oOV9VvVdWvTud6PnNV9d1V9Ynp/+dPVdVf0/d5q6p/OP1u/29V9UtVdaWer64DGTyr6qYktyV5bmHtLUnek+QvJbk9yb+qqsPLqZABHk/y1u7+K0n+IMmDib7P2dTHX0jyt5O8JcmPTf1mXl5K8lPd/f1J3pnkfVOfH0hytrtPJjk7nTMv70/y1MK5ns/fh5L8Wnf/xSRvy1b/9X2mquqGJP8gyanufmuSw9n6zKbnK+pABs8k/yLJzyRZ/ILrnUk+3t3f7O4/TPJ0kpuXURx7r7s/3d0vTae/meTG6bG+z9fNSZ7u7me6+1tJPp6tfjMj3X2uuz83Pf5Gtj6I3pCtXp+ZXnYmyV1LKZAhqurGJD+a5MMLy3o+Y1X1piR/I8lHkqS7v9XdX4++z92RJN9VVUeSvCHJV6PnK+vABc+qeneSr3T3b7/qqRuSfHnh/Plpjfn5u0n+/fRY3+dLbw+YqjqR5B1JnkhyfXefS7bCaZLrllgae++D2foD8ubCmp7P2/cm2Ujyr6dLrD9cVW+Mvs9Wd38lyT/P1hWK55L8z+7+dPR8ZR1ZdgEjVNV/TPI9F3jqZ5N8IMkPX+jHLrBmy98V8v/re3c/Or3mZ7N1ad7HvvNjF3i9vs+D3h4gVXVVkl9J8pPd/SdVF2o/c1BVdyQ5391PVtUPLrkcLp8jSX4gyU909xNV9aG4xHLWpu9u3pnkzUm+nuSXq+q9Sy2KXZll8Ozud11ovar+crb+4/3t6UPJjUk+V1U3Z2sactPCy2/M1jifFfF6ff+Oqro3yR1Jbu0/u4+Qvs+X3h4QVbWWrdD5se7+5LT8QlUd7+5zVXU8yfnlVcgeuyXJu6vqR5JcmeRNVfWL0fO5ez7J8939xHT+iWwFT32fr3cl+cPu3kiSqvpkkr8ePV9ZB+pS2+7+Qndf190nuvtEtn6J/UB3//ckjyV5T1VdUVVvTnIyyX9dYrnsoaq6Pck/SvLu7v7fC0/p+3x9NsnJqnpzVR3N1oYEjy25JvZYbf0V8SNJnurun1946rEk906P703y6OWujTG6+8HuvnH6//h7kvyn7n5v9HzWps9qX66q75uWbk3yu9H3OXsuyTur6g3T7/pbs/U9fj1fUbOceF6K7v5iVT2SrV9iLyV5X3e/vOSy2Dv/MskVSR6fpt2/2d1/T9/nq7tfqqq/n+Q/ZGsnvI929xeXXBZ775YkP57kC1X1+WntA0keSvJIVd2XrQ8vdy+nPC4jPZ+/n0jysemPic8k+TvZGqLo+wxNl1R/IsnnsvUZ7beSnE5yVfR8JdWfXXEIAAAAe+9AXWoLAADA5Sd4AgAAMJTgCQAAwFCCJwAAAEMJngAAAAwleAIAADCU4AkAAMBQgicAAABD/V+ajALsc8W5GgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "if (gravitational_acceleration != 0):\n",
    "    initialize_hydrostatic_pressure()\n",
    "    compute_density()\n",
    "    plt.scalar_field(dh.cpu_arrays[\"rho\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "================================= start of the simulation ===================================\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8c9dbae7a6154aabae39e76d60d326b2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/238800 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "time needed for the calculation: 99.6110 seconds\n",
      "MLUPS: 5.99\n"
     ]
    }
   ],
   "source": [
    "print(\"================================= start of the simulation ===================================\")\n",
    "\n",
    "start = time.time()\n",
    "\n",
    "pbar = tqdm(total=timesteps)\n",
    "\n",
    "timestep = []\n",
    "surface_position = []\n",
    "symmetry_norm = []\n",
    "mass = []\n",
    "\n",
    "for i in range(0, timesteps):\n",
    "    \n",
    "    sum_c_2 = 0.0\n",
    "    sum_delta_c_2 = 0.0\n",
    "    \n",
    "    # write vtk output\n",
    "    if(i % vtk_output_frequency == 0):\n",
    "        if gpu:\n",
    "            dh.to_cpu(\"C\")\n",
    "            dh.to_cpu(\"u\")\n",
    "            dh.to_cpu(\"rho\")\n",
    "        vtk_writer(i)\n",
    "    \n",
    "    # extract data (to be written to file)\n",
    "    if(i % data_extract_frequency == 0):\n",
    "        pbar.update(data_extract_frequency)\n",
    "            \n",
    "        timestep.append(i)\n",
    "        \n",
    "        ny = domain_size[1]\n",
    "\n",
    "        # get index containing phase field value < 0.5\n",
    "        i1 = np.argmax(dh.cpu_arrays[\"C\"][int(domain_size[0] * 0.5), :] < 0.5)\n",
    "        i0 = i1 - 1 # index containing phase field value >= 0.5    \n",
    "\n",
    "        f0 = dh.cpu_arrays[\"C\"][int(domain_size[0] * 0.5), i0] # phase field value >= 0.5\n",
    "        f1 = dh.cpu_arrays[\"C\"][int(domain_size[0] * 0.5), i1] # phase field value < 0.5\n",
    "\n",
    "        # coordinate of cell center is index+0.5-1 (0.5 to get to cell center; -1 to remove ghost layer from index)\n",
    "        y0 = i0 + 0.5 - 1\n",
    "        y1 = i1 + 0.5 - 1\n",
    "\n",
    "        #interpolate\n",
    "        surface_position.append( y0 + (y1 - y0) / (f1 - f0) * (0.5 - f0) )\n",
    "        \n",
    "        # evaluate symmetry in x-direction\n",
    "        for y in range(0, domain_size[1] - 1):\n",
    "            for x in range(0, domain_size[0] - 1):\n",
    "                if (x >= domain_size[0] * 0.5):\n",
    "                    continue\n",
    "                    \n",
    "                x_mirrored = domain_size[0] - 1 - x;\n",
    "                sum_c_2 += dh.cpu_arrays[\"C\"][x, y]**2\n",
    "                sum_delta_c_2 += (dh.cpu_arrays[\"C\"][x, y] - dh.cpu_arrays[\"C\"][x_mirrored, y])**2\n",
    "                \n",
    "        symmetry_norm.append( (sum_delta_c_2 / sum_c_2)**2 )\n",
    "        \n",
    "        mass.append(np.sum(dh.cpu_arrays[\"C\"]))\n",
    "        \n",
    "    timeloop()\n",
    "\n",
    "pbar.close()\n",
    "end = time.time()\n",
    "sim_time = end - start\n",
    "print(\"\\n\")\n",
    "print(\"time needed for the calculation: %4.4f\" % sim_time , \"seconds\")\n",
    "\n",
    "nrOfCells = np.prod(domain_size)\n",
    "mlups = nrOfCells * timesteps / sim_time * 1e-6\n",
    "\n",
    "print(\"MLUPS: %4.2f\" % mlups)"
   ]