{ "cells": [ { "cell_type": "markdown", "id": "782db3d5", "metadata": {}, "source": [ "# 2D polydisperse coarse graining of a Janssen pile" ] }, { "cell_type": "markdown", "id": "3ce2cbf0", "metadata": {}, "source": [ "In this example, we coarse grain a 2D MFiX simulation of a Janssen pile of spherical particles with several particle sizes (no fluid), resolving the coarse-grained fields per particle phase, and plot the vertical pressure.\n", "\n", "The files required for this example can be found in the repository:\n", "\n", "- The DEM simulation data files in `examples/data_examples/jannsen_pile/DEM_data/` (time steps 600 to 614)\n", "- The configuration file `config.ini` in `examples/jannsen_pile/`\n", "\n", "The coarse graining can also be run from the `compute_CG.ipynb` notebook in `examples/jannsen_pile/`." ] }, { "cell_type": "markdown", "id": "98050e46", "metadata": {}, "source": [ "## Import packages" ] }, { "cell_type": "code", "execution_count": 1, "id": "8f9ea4ee", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Hello from pysammos\n", "Loading macroscopic_fields package...\n", "Loading data_read package...\n", "Loading particle_phase package...\n", "Loading spatial_weights package...\n", "Loading neighbour_search package...\n", "Loading grid_generation package...\n", "Loading data_handle package...\n", "Loading data_write package...\n" ] } ], "source": [ "from pysammos.utils.config_loader import load_config\n", "from pysammos.coarse_graining import CoarseGraining\n", "from pysammos.data_write.h5.writer import H5XarrayManager" ] }, { "cell_type": "code", "execution_count": 2, "id": "f611218e", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import pyvista as pv" ] }, { "cell_type": "code", "execution_count": 3, "id": "3cf54bf2", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The history saving thread hit an unexpected error (OperationalError('database is locked')).History will not be written to the database.\n", "Working directory: /exports/csce/datastore/geos/users/s1857688/Coarse_Graining/pysammos/examples/jannsen_pile\n" ] } ], "source": [ "from pathlib import Path\n", "import os\n", "\n", "# The paths inside config.ini (DEM data, output folder) are relative to the example folder,\n", "# so we work from examples/jannsen_pile/ and read its config.ini.\n", "if \"EXAMPLE_DIR\" not in globals(): # resolve once, so the cell can be re-run\n", " EXAMPLE_DIR = (Path.cwd() / \"../../examples/jannsen_pile\").resolve()\n", "os.chdir(EXAMPLE_DIR)\n", "CONFIG_FILE = EXAMPLE_DIR / \"config.ini\"\n", "print(\"Working directory:\", Path.cwd())" ] }, { "cell_type": "markdown", "id": "cc9bcf88", "metadata": {}, "source": [ "## The configuration file\n", "\n", "For a description of every section, see the bedload transport example. The file used here is:" ] }, { "cell_type": "code", "execution_count": 4, "id": "5c2d96e6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[paths]\n", "particles_path = ../data_examples/jannsen_pile/DEM_data/DES_FB1_\n", "contacts_path = ../data_examples/jannsen_pile/DEM_data/FC_\n", "output_path = ./PysammosCG/\n", "\n", "[timesteps]\n", "t0 = 600\n", "tf = 614\n", "dt = 1\n", "t_pad = 4\n", "\n", "[smoothing_function]\n", "type = Lucy\n", "\n", "[flags]\n", "partialignore = False\n", "\n", "[key_mapping]\n", "Global_ID = Particle_ID\n", "Particle_Velocity = Velocity\n", "Particle_Diameter = Diameter\n", "Particle_Density = Density\n", "Particle_Volume = Volume\n", "Particle_Mass = Mass\n", "Particle_Radius = None\n", "Coordination_Number = None\n", "Particle_i_ID = Particle_ID_1\n", "Particle_j_ID = Particle_ID_2\n", "Force_ij = FORCE_CHAIN_FC\n", "Contact_ij = FORCE_CHAIN_CONTACT_POINT\n", "\n", "[DEM_software]\n", "name = mfix\n", " \n", "[grid_info]\n", "grid_dimension = 2\n", "grid_axes = xy\n", "automatic_grid = False\n", "x_min = -0.23\n", "x_max = 0.23\n", "y_min = 0.001\n", "y_max = 0.1\n", "z_min = None\n", "z_max = None\n", "x_transect = None\n", "y_transect = None\n", "z_transect = 0.0\n", "x_axis_periodic = False\n", "y_axis_periodic = False\n", "z_axis_periodic = False\n", "\n", "[fields_to_export]\n", "volume_fraction = True\n", "density_particle = False\n", "density_mixture = True\n", "momentum_density = False\n", "velocity = True\n", "velocity_gradient = False\n", "kinetic_tensor = False\n", "contact_tensor = False\n", "total_stress_tensor = True\n", "pressure = True\n", "granular_temperature = True\n", "granular_temperature_slices = True\n", "fabric_tensor = True\n", "inertial_number = False\n", "coordination_number = True\n", "d43 = False\n", "d32 = False\n", "frictional_coefficient = True\n", "shear_rate_tensor = False\n", "\n", "[output_options]\n", "vkthdf_output = True\n", "h5_output = True\n", "time_average = False\n", "\n" ] } ], "source": [ "print(CONFIG_FILE.read_text())" ] }, { "cell_type": "markdown", "id": "8723a658", "metadata": {}, "source": [ "Compared with the bedload transport example:\n", "\n", "- `partialignore = False`: the particle phases are **not** ignored. Pysammos groups the particles into phases by diameter and density, and computes every field for the bulk and for each phase (*polydisperse* mode). The output files are named `CG_Lucy_Polydisperse*`.\n", "- `x_min = -0.23` ... `y_max = 0.1`, with `z_transect = 0.0`: a 2D grid in the `xy` plane, with no periodic boundaries.\n", "- `coordination_number = True`: it is not in the MFiX output, so Pysammos calculates it from the contacts.\n", "- Time steps 600 to 614 are coarse grained (`t0`, `tf`, `dt`)." ] }, { "cell_type": "markdown", "id": "24e0a521", "metadata": {}, "source": [ "## Initialise the `CoarseGraining` class\n", "\n", "`load_config` returns the arguments of `CoarseGraining` as a dictionary, so they can be passed directly with `**`. You can also give each argument explicitly, as shown in the bedload transport example." ] }, { "cell_type": "code", "execution_count": 5, "id": "26f13b41", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Output path created: ./PysammosCG/\n" ] } ], "source": [ "cfg = load_config(CONFIG_FILE)\n", "CG = CoarseGraining(**cfg)" ] }, { "cell_type": "markdown", "id": "3649403a", "metadata": {}, "source": [ "## Run the coarse graining\n", "\n", "The whole workflow can be run with a single command:\n", "\n", "```python\n", "CG.run()\n", "```\n", "\n", "Here we run it step by step instead, to show what each step does." ] }, { "cell_type": "markdown", "id": "2ed49a2a", "metadata": {}, "source": [ "### 1. Load the size-relevant particle data for the first time step" ] }, { "cell_type": "code", "execution_count": 6, "id": "527baf7f", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Sampling data for the first time step to obtain particle properties and bounds ... \n", "XML-based PolyData detected.\n" ] } ], "source": [ "Bounds_t0, Diameter_t0, Density_t0, Mass_t0, GlobalID_t0 = CG.data_sampling()" ] }, { "cell_type": "markdown", "id": "be8d562b", "metadata": {}, "source": [ "### 2. Calculate the particle size statistics" ] }, { "cell_type": "code", "execution_count": 7, "id": "1bba86fe", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ ">> Particle size statistics: \n", " d43: 0.005067612399400709\n", " dmax: 0.005999645\n", " d50: 0.005061278\n", " d32: 0.005012527178597786\n", " drms: 0.004928972383396669\n" ] } ], "source": [ "CG.get_particle_size_statistics(Diameter_t0, Mass_t0)\n", "print(\">> Particle size statistics: \")\n", "print(\" d43: \", CG.d43)\n", "print(\" dmax: \", CG.dmax)\n", "print(\" d50: \", CG.d50)\n", "print(\" d32: \", CG.d32)\n", "print(\" drms: \", CG.drms)" ] }, { "cell_type": "markdown", "id": "1a09d5a9", "metadata": {}, "source": [ "### 3. Get the phases\n", "\n", "With `partialignore = False`, the particles are clustered into phases by diameter and density (here with up to 8 clusters). The figure shows the clustering." ] }, { "cell_type": "code", "execution_count": 8, "id": "e3c3291f", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe4AAAIhCAYAAACIQVmqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAbRFJREFUeJzt3Xl0E1X/BvBnkibpHrq3QGmL7DsFgQKy04JUFkXhJxYQBJRFVhUVXsBXBVxBVEBkUcrmqyyyWChbAcsmtYCCiEhZW8rWlra0TZP7+wMTCd2SNmk78HzOyTnNzJ2b79xM8nQmk4kkhBAgIiIiWVBUdAFERERkOQY3ERGRjDC4iYiIZITBTUREJCMMbiIiIhlhcBMREckIg5uIiEhGGNxEREQywuAmIiKSkYcyuFesWAFJkvDLL7+YTb9x4wZatmwJV1dXxMbGAgBmzpwJSZKgUCjw999/F+grKysL7u7ukCQJQ4cOLY/y7e7atWuYOnUqGjduDFdXVzg6OqJ27doYP348zp49a2pnHBt72bZtG2bOnGm3/jt16oROnTrZrX9bmTZtGmrUqAEHBwdUqVKlyHbG58N4U6vVCAkJwfjx45GWlmazerKzszFz5kzs3bu3wDzjayspKcmqPm29Le3du7fAWPj4+KBdu3Z4++23ceHCBZs9VlkFBwebvXdcvXoVM2fORGJiYoXVVFaWvraCg4NNz5FCoYBWq0X9+vUxePBg7Nixo9BlJEmy6/uCLdj7vaskDhX2yOXs8uXL6N69O65du4adO3eiTZs2ZvNdXV2xfPly/Pe//zWb/r///Q86nQ4qlao8y7WbI0eOIDIyEkIIjB07FmFhYVCr1Thz5gyio6PRqlUr3L59u1xq2bZtG7744gu7vQC+/PJLu/RrS5s2bcJ7772Ht99+Gz179oRGoylxmZiYGGi1Wty5cwfbtm3D/PnzceTIEcTHx9skHLOzszFr1iwAKPDm3KtXLxw8eBABAQFlfhxbeP/999G5c2fo9XrcvHkThw8fxrJly/Dpp59iyZIlGDRoUEWXiA0bNsDd3d10/+rVq5g1axaCg4PRrFmziiusnLRr1w4fffQRACAzMxNnzpzB2rVrERERgWeeeQZr1qwxe389ePAgqlevXlHlWsTe710leSSC++zZs+jWrRt0Oh3i4uLQuHHjAm0GDBiAb775BrNmzYJC8e+BiKVLl6Jfv3748ccfy7Nku8jIyECfPn3g6OiI+Ph4sxdHp06dMGrUKHz//fcVWKFtZGdnw9nZGQ0aNKjoUkr022+/AQBeffVV+Pr6WrRMixYt4O3tDQDo3r07bt68iZUrVyI+Ph7t2rUrdS1CCOTk5BTbxsfHBz4+PqV+DFurXbu22T/hvXv3xuTJk9GtWzcMHToUTZo0KfT1Xp6aN29eoY9f0apUqWL2HHXr1g1jxozBzJkzMWvWLEybNg1z5841zX9wp+pRYnzvKslDeaj8fomJiWjfvj0cHBxw4MCBIl/Ew4YNw6VLl0yH0AHgzz//xIEDBzBs2LBCl8nIyMCUKVMQEhICtVqNatWqYcKECcjKyjJr98UXX6BDhw7w9fWFi4sLGjdujA8++AA6nc6sXadOndCoUSMcPXoUTzzxBJydnVGzZk3MmTMHBoPB1M5gMODdd99F3bp14eTkhCpVqqBJkyaYP39+sWOxZMkSpKSk4IMPPijyP9r+/fsX20dRh7EePByYnZ1tGhtHR0d4enqiZcuWWLNmDQBg6NCh+OKLL0x9Gm/GQ7BCCHz55Zdo1qwZnJyc4OHhgf79+xf4OMM4Zvv27UPbtm3h7Oxser4ePJyXlJQESZLw0Ucf4ZNPPkFISAhcXV0RFhaGQ4cOFTpederUgUajQYMGDbB69WoMHToUwcHBxY4RcO85+uCDD1CvXj1oNBr4+vpi8ODBuHz5stmYTZs2DQDg5+dX6kOExje6CxcuICcnB5MnT0azZs2g1Wrh6emJsLAwbNq0qcBykiRh7NixWLRoEerXrw+NRoNvvvnGFMyzZs0yPS/G57aoQ+UxMTHo2rUrtFotnJ2dUb9+fcyePbvE2tetW4ewsDC4uLjA1dUVERER+PXXX60eg/t5enpi8eLFyM/Px6effmo27+zZs3j++efh6+sLjUaD+vXrm7ZDI+Nh+DVr1uDtt99G1apV4e7ujm7duuHMmTNmbX/99VdERkaa+qtatSp69epV4Hk2jt/evXvx+OOPAwBefPFF0/jOnDkTK1euhCRJOHjwYIF1euedd6BSqXD16tUi1/uvv/7Ciy++iNq1a8PZ2RnVqlXDU089hZMnT5Z6/YQQ+OCDDxAUFARHR0eEhobip59+KrIGa8ycORMNGzbE559/bvYP44Ovg+vXr2P06NFo0KABXF1d4evriy5dumD//v1m/Rlf3x9++CHmzp2L4OBgODk5oVOnTvjzzz+h0+kwdepUVK1aFVqtFv369UNqamqBukraJsvjvatE4iG0fPlyAUB8+umnQqvVikaNGomrV68W2nbGjBkCgLh+/bp44oknxHPPPWea98Ybb4jg4GBhMBiEi4uLGDJkiGleVlaWaNasmfD29haffPKJ2Llzp5g/f77QarWiS5cuwmAwmNpOnDhRLFy4UMTExIjdu3eLTz/9VHh7e4sXX3zRrJaOHTsKLy8vUbt2bbFo0SIRGxsrRo8eLQCIb775xtRu9uzZQqlUihkzZohdu3aJmJgYMW/ePDFz5sxixyU8PFwolUqRmZlp0Tgax+Z+AMSMGTMKtA0KCjIbn1GjRglnZ2fxySefiD179ogtW7aIOXPmiAULFgghhPjrr79E//79BQBx8OBB0y0nJ0cIIcSIESOESqUSkydPFjExMWL16tWiXr16ws/PT6SkpJiNmaenpwgMDBQLFiwQe/bsEXFxcaZ5HTt2NLU9f/68ACCCg4NFjx49xMaNG8XGjRtF48aNhYeHh0hLSzO1Xbx4sQAgnnnmGbFlyxaxatUqUadOHREUFCSCgoJKHLuRI0cKAGLs2LEiJiZGLFq0SPj4+IjAwEBx/fp1IYQQCQkJYvjw4QKAiImJEQcPHhSXLl0q8fkwLm80ceJEAUDs2LFDpKWliaFDh4qVK1eK3bt3i5iYGDFlyhShUCjMtiEh7j2X1apVE02aNBGrV68Wu3fvFomJiSImJkYAEMOHDzc9L3/99ZcQ4t/X1vnz5039fP3110KSJNGpUyexevVqsXPnTvHll1+K0aNHF6j9fu+9956QJEkMGzZMbNmyRaxfv16EhYUJFxcX8fvvvxc7vnv27BEAxP/+978i2wQEBIjHHnvMdP/3338XWq1WNG7cWHz77bdix44dYvLkyUKhUJi9dox9BwcHi0GDBomtW7eKNWvWiBo1aojatWuL/Px8IYQQmZmZwsvLS7Rs2VJ89913Ii4uTqxbt068/PLL4tSpU6b+7n9tpKenm8Zw2rRppvG9dOmSyM3NFf7+/mLQoEFm66HT6UTVqlXFs88+W+yYxMXFicmTJ4vvv/9exMXFiQ0bNoi+ffsKJycn8ccff1i9fkL8+7wNHz5c/PTTT+Krr74S1apVE/7+/mavraIEBQWJXr16FTl/6tSpAoDYv3+/adqD7zF//PGHeOWVV8TatWvF3r17xZYtW8Tw4cOFQqEQe/bsMbUzvr6DgoLEU089JbZs2SKio6OFn5+fqFOnjoiKihLDhg0TP/30k1i0aJFwdXUVTz31lFk9lmyT5fHeVZKHOrgBCK1WK1JTU4tse/+b4fLly4VGoxE3b94U+fn5IiAgwPSCfjC4Z8+eLRQKhTh69KhZf99//70AILZt21bo4+n1eqHT6cS3334rlEqluHXrlmlex44dBQBx+PBhs2UaNGggIiIiTPcjIyNFs2bNLB4Po3r16gl/f3+L25cluBs1aiT69u1bbP9jxowp0L8QQhw8eFAAEB9//LHZ9EuXLgknJyfx+uuvm6YZx2zXrl0F+ikquBs3bmz25nTkyBEBQKxZs0YIce858vf3F61btzbr78KFC0KlUpUY3KdPnxYAzIJLCCEOHz4sAIi33nrLNK2oMC6MsW1KSorQ6XTi9u3bIjo6Wjg5OYnAwEBx9+7dAsvk5+cLnU4nhg8fLpo3b242z/j6uH8bFEKI69evF/k8Pxjcd+7cEe7u7qJ9+/Zm/6wWVbvRxYsXhYODgxg3bpxZuzt37gh/f3+zf6ALY0lwt27dWjg5OZnuR0REiOrVq4v09HSzdmPHjhWOjo6mcTD2/eSTT5q1++6770xv1kII8csvvwgAYuPGjcXW+uBr4+jRowKAWL58eYG2M2bMEGq1Wly7ds00bd26dQKAxW/qRvn5+SIvL0/Url1bTJw40TTd0vW7ffu2cHR0FP369TNr9/PPPwsANgnuhQsXCgBi3bp1pmlFbXv3r5dOpxNdu3Y1q834+m7atKnQ6/Wm6fPmzRMARO/evc36mTBhggBg2h6s2Sbt/d5Vkof6UHnv3r2Rnp6OCRMmQK/Xl9j+2WefhVqtxqpVq7Bt2zakpKQUeSb5li1b0KhRIzRr1gz5+fmmW0REBCRJMjsj99dff0Xv3r3h5eUFpVIJlUqFwYMHQ6/X488//zTr19/fH61atTKb1qRJE7OzZFu1aoXjx49j9OjR2L59OzIyMiwflHLSqlUr/PTTT5g6dSr27t2Lu3fvWrzsli1bIEkSXnjhBbOx9ff3R9OmTQuc7ezh4YEuXbpY3H+vXr2gVCpN95s0aQIApjE+c+YMUlJS8Nxzz5ktV6NGDYs+Q96zZw8AFNh2WrVqhfr162PXrl0W11oYf39/qFQqeHh44IUXXkBoaChiYmLg6OgI4N4Jle3atYOrqyscHBygUqmwdOlSnD59ukBfXbp0gYeHR6lriY+PR0ZGBkaPHm3ViXHbt29Hfn4+Bg8ebPYcOzo6omPHjoWe0W4tIYTp75ycHOzatQv9+vWDs7Oz2WM++eSTyMnJKfBxSe/evc3uP7id1KpVCx4eHnjjjTewaNEinDp1qsw1v/LKKwDufUxj9Pnnn6Nx48bo0KFDscvm5+fj/fffR4MGDaBWq+Hg4AC1Wo2zZ88W+tyXtH4HDx5ETk5OgRP82rZti6CgIOtXrhD3P0fFWbRoEUJDQ+Ho6Gjapnft2lXoej355JNm5ynVr18fwL3X/f2M0y9evAjANtukvd+7jB7qk9OmT5+OZs2a4Z133oHBYEB0dLTZG/aDXFxcMGDAACxbtgxBQUHo1q1bkRvotWvX8NdffxV5tvmNGzcA3NsonnjiCdStWxfz589HcHAwHB0dceTIEYwZM6ZAoHl5eRXoS6PRmLV788034eLigujoaCxatAhKpRIdOnTA3Llz0bJlyyLXr0aNGjh79iyysrLg4uJSZDtb+Oyzz1C9enWsW7cOc+fOhaOjIyIiIvDhhx+idu3axS577do1CCHg5+dX6PyaNWua3bf2DOcHx9h4JrdxjG/evAkAhT6+n58fzp8/X2z/xuULq6tq1apl/qrSzp07odVqoVKpUL16dbP1Wb9+PZ577jk8++yzeO211+Dv7w8HBwcsXLgQy5YtK9BXWc8Ov379OgBYfRbwtWvXAMD0ee+D7n/jLa2LFy+iatWqAO49J/n5+ViwYAEWLFhQaHvja9aopO1Eq9UiLi4O7733Ht566y3cvn0bAQEBGDFiBKZNm1aqb6L4+flhwIABWLx4MaZOnYrff/8d+/fvx+LFi0tcdtKkSfjiiy/wxhtvoGPHjvDw8IBCocBLL71U6D/Olr4O/P39Cyxb2LTSML4WjM9TYT755BNMnjwZL7/8Mv773//C29sbSqUS06dPLzS4PT09ze6r1epipxs/X7fFNmnv9y6jhzq4gX9PsJk1axYMBgNWrVoFB4eiV3vYsGH4+uuvceLECaxatarIdt7e3nBycir0zdA4HwA2btyIrKwsrF+/3uyfgLJ8h9PBwQGTJk3CpEmTkJaWhp07d+Ktt95CREQELl26VORZiREREdixYwc2b96MgQMHluqxNRoNcnNzC0w3vsiNXFxcMGvWLMyaNQvXrl0z7X0/9dRT+OOPP4p9DG9vb0iShP379xf69agHp9n6u+bGNzTjC/l+KSkpFi+fnJxcINCuXr1q2jZKq2nTpkX2ER0djZCQEKxbt85sXAp7zoCyj53xRLb7T8ayhLH+77//3mZ7b/c7cuQIUlJSMHz4cAD39myUSiWioqIwZsyYQpcJCQmx+nEaN26MtWvXQgiBEydOYMWKFXjnnXfg5OSEqVOnlqr28ePHY+XKldi0aRNiYmJQpUoVi77WFh0djcGDB+P99983m37jxo1irw9QFON2XNg2n5KSYtFJmsURQmDz5s1wcXEpdocjOjoanTp1wsKFC82m37lzp0yP/yBbbJPl9d710Ac3cO/sRYVCgRkzZkAIgdWrVxcZ3mFhYRg2bBjS09PRr1+/IvuMjIzE+++/Dy8vr2Jf8MYn5v4nTAhhdiisLKpUqYL+/fvjypUrmDBhApKSkor8GtTw4cPx4Ycf4vXXX8cTTzyBatWqFWizfv16PP3000U+XnBwME6cOGE2bffu3cjMzCxyGT8/PwwdOhTHjx/HvHnzTF95uP8/fCcnJ1P7yMhIzJkzB1euXClwuLo81K1bF/7+/vjuu+8wadIk0/SLFy8iPj6+2L0DAKZDX9HR0Wb/vR89ehSnT5/G22+/bZ/CAdPFSO5/Q0hJSSn0rPKiPLjnVZy2bdtCq9Vi0aJFGDhwoMVvRBEREXBwcMC5c+fwzDPPWFybJW7duoWXX34ZKpUKEydOBAA4Ozujc+fO+PXXX9GkSRPT3patSJKEpk2b4tNPP8WKFSuQkJBQZNuSxrdFixZo27Yt5s6di99++w0jR4606AiZJEkFgmHr1q24cuUKatWqZcXa3NOmTRs4Ojpi1apVZs9RfHw8Lly4UObgnjVrFk6dOoW33nrL9DFPYQpbrxMnTuDgwYMIDAwsUw33s2abrOj3rkciuAHgP//5DxQKBaZPnw4hBNasWVNkeC9durTE/iZMmIAffvgBHTp0wMSJE9GkSRMYDAZcvHgRO3bswOTJk9G6dWt0794darUa//d//4fXX38dOTk5WLhwYZkucvLUU0+hUaNGaNmyJXx8fHDhwgXMmzcPQUFBxR6G1mq12LRpEyIjI9G8eXOzC7CcPXsW0dHROH78eLHBHRUVhenTp+M///kPOnbsiFOnTuHzzz+HVqs1a9e6dWtERkaiSZMm8PDwwOnTp7Fy5UqEhYWZjggYv5o3d+5c9OzZE0qlEk2aNEG7du0wcuRIvPjii/jll1/QoUMHuLi4IDk52fSVPuNngfagUCgwa9YsjBo1Cv3798ewYcOQlpaGWbNmISAgoMRDZnXr1sXIkSOxYMECKBQK9OzZE0lJSZg+fToCAwNNYWIPkZGRWL9+PUaPHo3+/fvj0qVL+O9//4uAgACzq+IVx83NDUFBQdi0aRO6du0KT09PeHt7F/pG7erqio8//hgvvfQSunXrhhEjRsDPzw9//fUXjh8/js8//7zQxwgODsY777yDt99+G3///Td69OgBDw8PXLt2DUeOHDEdsSnJ2bNncejQIRgMBtMFWJYuXYqMjAx8++23aNiwoant/Pnz0b59ezzxxBN45ZVXEBwcjDt37uCvv/7C5s2bsXv3bovGx2jLli348ssv0bdvX9SsWRNCCKxfvx5paWno3r17kcs99thjcHJywqpVq1C/fn24urqiatWqZv8Qjh8/HgMGDIAkSRg9erRF9URGRmLFihWoV68emjRpgmPHjuHDDz8s9cVMPDw8MGXKFLz77rt46aWX8Oyzz+LSpUuYOXOmVYfK09LSTOcPZGVlmS7Asn//fjz33HMlPs+RkZH473//ixkzZqBjx444c+YM3nnnHYSEhCA/P79U61YYa7bJCn/vsvp0Nhkwnvn64BnfQtw73R+AePrpp0VeXp7FZ/U+eFa5EPe+DjJt2jRRt25doVarTV81mThxotlp/5s3bxZNmzYVjo6Oolq1auK1114TP/30kwBg9nWGjh07ioYNGxZ47CFDhpidyfzxxx+Ltm3bCm9vb6FWq0WNGjXE8OHDRVJSkkXjk5KSIt544w3RsGFD4ezsLDQajahVq5YYNWqUOHnypKldYWeV5+bmitdff10EBgYKJycn0bFjR5GYmFjgzNmpU6eKli1bCg8PD6HRaETNmjXFxIkTxY0bN8z6eumll4SPj4+QJKnA14yWLVsmWrduLVxcXISTk5N47LHHxODBg8Uvv/xS4pgZ5xV2VvmHH35YoC0KOZP1q6++ErVq1RJqtVrUqVNHLFu2TPTp06fA2dmF0ev1Yu7cuaJOnTpCpVIJb29v8cILLxT4uldpziovqe2cOXNEcHCw0Gg0on79+mLJkiVFfkNgzJgxhfaxc+dO0bx5c6HRaAQA03Nb2NfBhBBi27ZtomPHjsLFxUU4OzuLBg0aiLlz5xao/UEbN24UnTt3Fu7u7kKj0YigoCDRv39/sXPnzmLX0XhmtPHm4OAgvLy8RFhYmHjrrbeKfC2cP39eDBs2TFSrVk2oVCrh4+Mj2rZtK959990CfT94xrpx+zGeDf7HH3+I//u//xOPPfaYcHJyElqtVrRq1UqsWLHCbLkHXxtCCLFmzRpRr149oVKpCt32cnNzhUajET169Ch2HO53+/ZtMXz4cOHr6yucnZ1F+/btxf79+wu8DixdPyGEMBgMYvbs2SIwMFCo1WrRpEkTsXnz5gJ9FiUoKMj0HEmSJFxdXUXdunVFVFSU2L59e6HLPDgeubm5YsqUKaJatWrC0dFRhIaGio0bNxZ4Xyzq9V3U+haVE5Zsk/Z+7yqJ9M9AEVEJ0tLSUKdOHfTt2xdfffVVRZdDD7HNmzejd+/e2Lp1K5588smKLocqGQY3USFSUlLw3nvvoXPnzvDy8sKFCxfw6aef4o8//sAvv/xidgiWyFZOnTqFCxcuYPz48XBxcUFCQoJdf+iH5OmR+YybyBoajQZJSUkYPXo0bt26BWdnZ7Rp0waLFi1iaJPdjB49Gj///DNCQ0PxzTffMLSpUNzjJiIikpGH+sppREREDxsGNxERkYwwuImIiGSEJ6dZyGAw4OrVq3Bzc+MJI0REjyghBO7cuYOqVava5Jr6pcHgttDVq1dtenk9IiKSr0uXLpX6qnRlxeC2kJubG4B7T5a7u3up+tDpdNixYwfCw8NL9ctB5Y312pfc6gXkVzPrta9Hsd6MjAwEBgaaMqEiMLgtZDw87u7uXqbgdnZ2hru7u2w2ctZrP3KrF5BfzazXvh7leivyI1OenEZERCQjDG4iIiIZYXATERHJCD/jJiKyESEE8vPzodfrS7W8TqeDg4MDcnJySt1HeXpY61WpVFAqleVYmXUY3ERENpCXl4fk5GRkZ2eXug8hBPz9/XHp0iVZXC/iYa1XkiRUr14drq6u5Vid5RjcRERlZDAYcP78eSiVSlStWhVqtbpUQWYwGJCZmQlXV9cKu7iHNR7GeoUQuH79Oi5fvozatWtXyj1vBjcRURnl5eXBYDAgMDAQzs7OVi0rhDCFvMFgQF5eHhwdHQsEy/3tKovi6q2MLK3Xx8cHSUlJ0Ol0lTK4K/9IExHJhLXhtXv3bvj5+eGrr74qtt1XX30FPz8/7N69uyzlkYUq2z9ID2JwExFVgN27d+PJJ5/EjRs3MGrUqCLDe/HixRg1ahRu3LiBJ598kuFNDG4iovJmDG2dTgchBAAUGt6LFy/Gyy+/DODeoXKdTsfwJgY3EVF5uj+0DQaD2bxXXnkFK1asAGAe2kYGg4HhDWDv3r2QJAlpaWlFtlmxYgWqVKlSbjWVJwY3EVE5KS60jSZOnIgXXnihQGgb2SO8U1JSMG7cONSsWRMajQaBgYF46qmnsGvXLpv0DwCdOnXChAkTbNJX27ZtkZycDK1Wa5P+5IZnlRMRlQMhBAYOHIi8vDzT4fGirFmzptj5xrOjBw4ciGvXrpXpZKqkpCS0a9cOVapUwQcffIAmTZpAp9Nh+/btGDNmDP74449S920tIQT0ej0cHIqPJrVaDX9//3KqqvLhHjcRUTmQJAnvvvtuiaFtKSEE3n333TKfAT169GhIkoQjR46gf//+qFOnDho2bIhJkybh0KFDAID09HSMHDkSvr6+cHd3R5cuXXD8+HFTH7NmzUKzZs2wcuVKBAcHQ6vVYuDAgbhz5w4AYOjQoYiLi8P8+fMhSRIkSUJSUpLpkPf27dvRsmVLaDQa7N+/H7m5uXj11Vfh6+sLR0dHtG/fHkePHjU9XmGHylesWIEaNWrA2dkZ/fr1w82bN83W8/jx4+jatSsCAwNRpUoVtGjRAr/88kuZxq6iMLiJiMrJyJEjsWjRIpv0tXjxYowcObJMfdy6dQsxMTEYM2YMXFxcCsyvUqUKhBDo1asXUlJSsG3bNhw7dgyhoaHo2rUrbt26ZWp77tw5bNy4EVu2bMGWLVsQFxeHOXPmAADmz5+PsLAwjBgxAsnJyUhOTkZgYKBp2ddffx2zZ8/G6dOn0aRJE7z++uv44Ycf8M033yAhIQG1atVCRESE2ePd7/Dhwxg2bBhGjx6NxMREdO7cGe+++65Zm0GDBqFatWrYtWsXjh49iqlTp8rip0gLw0PlRETlaNSoUQBQ5GfYlrBFaAPAX3/9BSEE6tWrV2SbPXv24OTJk0hNTYVGowEAfPTRR9i4cSO+//57DBw4EMC9w/crVqyAm5sbACAqKgq7du3Ce++9B61WC7VaDWdn50IPcb/zzjvo3r07ACArKwsLFy7EihUr0LNnTwDAkiVLEBsbi6VLl+K1114rsPz8+fMRERGBqVOnAgDq1KmD+Ph4xMTEmNpcvHgRkydPRp06deDu7o66deuWZsgqBe5xExGVs1GjRpV6z9tWoQ3AdNi+uMPtx44dQ2ZmJry8vODq6mq6nT9/Hn///bepXXBwsCm0ASAgIACpqakW1dGyZUvT3+fOnYNOp0O7du1M01QqFVq1aoXTp08Xuvzp06cRFhZmNu3B+5MmTcLIkSPRt29fzJ07F+fOnbOotsqIwU1EVAFGjRqF559/3qplBg0aZLPQBoDatWtDkqQiAxG4tycdEBCAxMREs9uZM2cwZcoUU7sHDztLklTkmfMPuv8wfVH/TBR3yVdLzhuYOXMmTp48ifDwcOzevRsNGjTAhg0bLKqvsmFwExFVgMWLF2P16tVWLbNq1aoSL49qDU9PT0REROCLL75AVlZWgflpaWkIDQ1FSkoKHBwcUKtWLbObt7e3xY+lVqst+unPWrVqQa1W48CBA6ZpOp0Ov/zyC+rXr1/oMg0aNDCdSGf04H3g3iH00aNHY/v27Xj66aexfPlyi+uvTBjcRETlrLCLq1iquMujlsaXX34JvV6PVq1a4YcffsDZs2dx+vRpfPbZZwgLC0O3bt0QFhaGvn37Yvv27UhKSkJ8fDymTZtm1VnZwcHBOHz4MJKSknDjxo0i98ZdXFzwyiuv4LXXXkNMTAxOnTqFESNGIDs7G8OHDy90mVdffRUxMTH44IMP8Oeff+Lzzz83+3z77t27GDt2LPbu3YuLFy/i559/xtGjR4v8R6CyY3ATEZWjsoS2kS3DOyQkBAkJCejcuTMmT56MRo0aoXv37ti1axcWLlwISZKwbds2dOjQAcOGDUOdOnUwcOBAJCUlwc/Pz+LHmTJlCpRKJRo0aAAfHx9cvHixyLZz5szBM888g6ioKISGhuKvv/7C9u3b4eHhUWj7Nm3a4Ouvv8aCBQvQrFkz7NixA9OmTTPNVyqVuHnzJoYOHYrHH38cAwcORM+ePTFr1izLB6oyEWSR9PR0AUCkp6eXuo+8vDyxceNGkZeXZ8PK7If12pfc6hVCfjWXV713794Vp06dEnfv3i223eLFiwUAm90WL15s1/UqiV6vF7dv3xZ6vb5C67CUpfUW93zaIgvKinvcRETlQAiBadOm2ewnIyVJwrRp02x2QReSDwY3EVE5kCQJa9euhVqtLvF3u0s621yhUECtVmPt2rWV/rejyfYY3ERE5aRLly7Ytm0bVCpVkeH96aefYuXKlUV+z1uhUEClUmHbtm3o0qWLPculSorBTURUjooL74ULF2Lo0KEACr9IC0ObAAY3EZHNWPp58/3hbTzUXdgV0e4Pb0mSGNrlpLKfN8DgJiIqI+NVw7Kzsy1exhje3t7exV7GdNSoUVi8eDG8vb0Z2uUkLy8PwL2vkVVG/JERIqIyUiqVqFKliuna3M7OzhadNNa2bVtcuHABkiQhJyfH9DvbOTk5ZofRBw8ejKioKFO7yqKoeisrS+o1GAy4fv06nJ2dS/xd8IpSOasiIpIZ469eWfrDGoURQuDu3btwcnKSxdniD2u9CoUCNWrUqLTrxOAmIrIBSZIQEBAAX19f6HS6UvWh0+mwb98+dOjQQRa/Ff2w1mvJV/YqEoObiMiGlEplqT8bVSqVyM/Ph6OjoyyCkPVWjMr7LwUREREVwOAmIiKSEQY3ERGRjDC4iYiIZITBTUREJCMMbiIiIhlhcBMREckIg5uIiEhGGNxEREQywuAmIiKSEQY3ERGRjDC4iYiIZITBTUREJCMMbiIiIhlhcBMREckIg5uIiEhGGNxEREQywuAmIiKSEQY3ERGRjDC4iYiIZITBTUREJCMMbiIiIhlhcBMREckIg5uIiEhGKjS4Z8+ejccffxxubm7w9fVF3759cebMmQLtTp8+jd69e0Or1cLNzQ1t2rTBxYsXTfNzc3Mxbtw4eHt7w8XFBb1798bly5fN+rh9+zaioqKg1Wqh1WoRFRWFtLQ0e68iERGRTVVocMfFxWHMmDE4dOgQYmNjkZ+fj/DwcGRlZZnanDt3Du3bt0e9evWwd+9eHD9+HNOnT4ejo6OpzYQJE7BhwwasXbsWBw4cQGZmJiIjI6HX601tnn/+eSQmJiImJgYxMTFITExEVFRUua4vERFRWTlU5IPHxMSY3V++fDl8fX1x7NgxdOjQAQDw9ttv48knn8QHH3xgalezZk3T3+np6Vi6dClWrlyJbt26AQCio6MRGBiInTt3IiIiAqdPn0ZMTAwOHTqE1q1bAwCWLFmCsLAwnDlzBnXr1rX3qhIREdlEhQb3g9LT0wEAnp6eAACDwYCtW7fi9ddfR0REBH799VeEhITgzTffRN++fQEAx44dg06nQ3h4uKmfqlWrolGjRoiPj0dERAQOHjwIrVZrCm0AaNOmDbRaLeLj4wsN7tzcXOTm5pruZ2RkAAB0Oh10Ol2p1s+4XGmXL2+s177kVi8gv5pZr309ivVWhnWtNMEthMCkSZPQvn17NGrUCACQmpqKzMxMzJkzB++++y7mzp2LmJgYPP3009izZw86duyIlJQUqNVqeHh4mPXn5+eHlJQUAEBKSgp8fX0LPKavr6+pzYNmz56NWbNmFZi+Y8cOODs7l2ldY2Njy7R8eWO99iW3egH51cx67etRqjc7O9uGlZROpQnusWPH4sSJEzhw4IBpmsFgAAD06dMHEydOBAA0a9YM8fHxWLRoETp27Fhkf0IISJJkun//30W1ud+bb76JSZMmme5nZGQgMDAQ4eHhcHd3t27l/qHT6RAbG4vu3btDpVKVqo/yxHrtS271AvKrmfXa16NYr/Hoa0WqFME9btw4/Pjjj9i3bx+qV69umu7t7Q0HBwc0aNDArH39+vVNAe/v74+8vDzcvn3bbK87NTUVbdu2NbW5du1agce9fv06/Pz8Cq1Jo9FAo9EUmK5Sqcq8gdqij/LEeu1LbvUC8quZ9drXo1RvZVjPCj2rXAiBsWPHYv369di9ezdCQkLM5qvVajz++OMFviL2559/IigoCADQokULqFQqs0MfycnJ+O2330zBHRYWhvT0dBw5csTU5vDhw0hPTze1ISIikoMK3eMeM2YMVq9ejU2bNsHNzc30ebNWq4WTkxMA4LXXXsOAAQPQoUMHdO7cGTExMdi8eTP27t1rajt8+HBMnjwZXl5e8PT0xJQpU9C4cWPTWeb169dHjx49MGLECCxevBgAMHLkSERGRvKMciIikpUK3eNeuHAh0tPT0alTJwQEBJhu69atM7Xp168fFi1ahA8++ACNGzfG119/jR9++AHt27c3tfn000/Rt29fPPfcc2jXrh2cnZ2xefNmKJVKU5tVq1ahcePGCA8PR3h4OJo0aYKVK1eW6/oSERGVVYXucQshLGo3bNgwDBs2rMj5jo6OWLBgARYsWFBkG09PT0RHR1tdIxERUWXCa5UTERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJSIUG9+zZs/H444/Dzc0Nvr6+6Nu3L86cOWPWZujQoZAkyezWpk0bszadOnUq0GbgwIFmbW7fvo2oqChotVpotVpERUUhLS3N3qtIRERkUxUa3HFxcRgzZgwOHTqE2NhY5OfnIzw8HFlZWWbtevTogeTkZNNt27ZtBfoaMWKEWZvFixebzX/++eeRmJiImJgYxMTEIDExEVFRUXZdPyIiIltzqMgHj4mJMbu/fPly+Pr64tixY+jQoYNpukajgb+/f7F9OTs7F9nm9OnTiImJwaFDh9C6dWsAwJIlSxAWFoYzZ86gbt26ZVwTIiKi8lGhwf2g9PR0AICnp6fZ9L1798LX1xdVqlRBx44d8d5778HX19eszapVqxAdHQ0/Pz/07NkTM2bMgJubGwDg4MGD0Gq1ptAGgDZt2kCr1SI+Pr7Q4M7NzUVubq7pfkZGBgBAp9NBp9OVav2My5V2+fLGeu1LbvUC8quZ9drXo1hvZVhXSQghKroIABBCoE+fPrh9+zb2799vmr5u3Tq4uroiKCgI58+fx/Tp05Gfn49jx45Bo9EAuLf3HBISAn9/f/z222948803UatWLcTGxgIA3n//faxYsQJ//vmn2WPWqVMHL774It58880C9cycOROzZs0qMH316tVwdna25aoTEZFMZGdn4/nnn0d6ejrc3d0rpIZKs8c9duxYnDhxAgcOHDCbPmDAANPfjRo1QsuWLREUFIStW7fi6aefBnDv8+3729SuXRstW7ZEQkICQkNDAQCSJBV4TCFEodMB4M0338SkSZNM9zMyMhAYGIjw8PBSP1k6nQ6xsbHo3r07VCpVqfooT6zXvuRWLyC/mlmvfT2K9RqPvlakShHc48aNw48//oh9+/ahevXqxbYNCAhAUFAQzp49W2Sb0NBQqFQqnD17FqGhofD398e1a9cKtLt+/Tr8/PwK7UOj0Zj26O+nUqnKvIHaoo/yxHrtS271AvKrmfXa16NUb2VYzwo9q1wIgbFjx2L9+vXYvXs3QkJCSlzm5s2buHTpEgICAops8/vvv0On05nahIWFIT09HUeOHDG1OXz4MNLT09G2bduyrwgREVE5qdA97jFjxmD16tXYtGkT3NzckJKSAgDQarVwcnJCZmYmZs6ciWeeeQYBAQFISkrCW2+9BW9vb/Tr1w8AcO7cOaxatQpPPvkkvL29cerUKUyePBnNmzdHu3btAAD169dHjx49MGLECNPXxEaOHInIyEieUU5ERLJSoXvcCxcuRHp6Ojp16oSAgADTbd26dQAApVKJkydPok+fPqhTpw6GDBmCOnXq4ODBg6YzxtVqNXbt2oWIiAjUrVsXr776KsLDw7Fz504olUrTY61atQqNGzdGeHg4wsPD0aRJE6xcubJC1puIiKi0KnSPu6QT2p2cnLB9+/Zi2wQGBiIuLq7Ex/L09ER0dLRV9REREVU2vFY5ERGRjDC4iYiIZITBTUREJCMMbiIiIhlhcBMREckIg5uIiEhGGNxEREQywuAmIiKSEQY3ERGRjDC4iYiIZITBTUREJCMMbiIiIhlhcBMREckIg5uIiEhGGNxEREQywuAmIiKSEQY3ERGRjDhYu0Bubi6OHDmCpKQkZGdnw8fHB82bN0dISIg96iMiIqL7WBzc8fHxWLBgATZu3Ii8vDxUqVIFTk5OuHXrFnJzc1GzZk2MHDkSL7/8Mtzc3OxZMxER0SPLokPlffr0Qf/+/VGtWjVs374dd+7cwc2bN3H58mVkZ2fj7NmzmDZtGnbt2oU6deogNjbW3nUTERE9kiza4w4PD8f//vc/qNXqQufXrFkTNWvWxJAhQ/D777/j6tWrNi2SiIiI7rEouMeMGWNxhw0bNkTDhg1LXRAREREVrcxnlQshYDAYbFELERERlcDi4M7Pz8e0adPQsWNHzJgxAwDw4YcfwtXVFU5OThgyZAjy8vLsVigRERFZcVb5rFmz8PXXX2PQoEH4/vvvkZqaiq1bt+Krr76CwWDAW2+9hXnz5uH111+3Z71ERESPNIuDe/Xq1fj6668RGRmJV155BXXr1sXq1asxYMAAAICjoyPeeecdBjcREZEdWXyo/OrVq2jatCkAoFatWlCr1ab7ANCyZUtcuHDB9hUSERGRicXBrdVqkZaWZrofGhpqdqGV3NxcSJJk0+KIiIjInMXB3aBBAyQkJJju//zzz6hWrZrp/smTJ1G7dm3bVkdERERmLP6Me9GiRVCpVEXO1+l0/HybiIjIziwO7jp16hQ7//nnny9zMURERFQ8/qwnERGRjFgd3B999JE96iAiIiILWBXcU6dOxZIlS+xVCxEREZXAos+4hRAYOXIkdu7ciX379tm7JiIiIiqCRcHdv39/HDp0CHFxcQgMDLR3TURERFQEi4J7w4YN+Oqrr1CrVi1710NERETFsOgz7gkTJmDy5Mk4evSoveshIiKiYli0x/3JJ5/Ay8sLPXv2xJ49e9C4cWN710VERESFsPgCLG+//TZ8fHzQo0cPXLlyxZ41ERERUREsDm4AGDlyJLy8vOxVCxEREZXA6guwPPPMM/aog4iIiCzAS54SERHJiFWHyo2OHDmCvXv3IjU1FQaDwWzeJ598YpPCiIiIqCCrg/v999/HtGnTULduXfj5+UGSJNO8+/8mIiIi27M6uOfPn49ly5Zh6NChdiiHiIiIimP1Z9wKhQLt2rWzRy1ERERUAquDe+LEifjiiy/sUQsRERGVwOpD5VOmTEGvXr3w2GOPoUGDBlCpVGbz169fb7PiiIiIyJzVwT1u3Djs2bMHnTt3hpeXF09IIyIiKkdWB/e3336LH374Ab169bJHPURERFQMqz/j9vT0xGOPPWaPWoiIiKgEVgf3zJkzMWPGDGRnZ9ujHiIiIiqG1YfKP/vsM5w7dw5+fn4IDg4ucHJaQkKCzYojIiIic1YHd9++fe1QBhEREVnC4uD+888/UadOHcyYMcOe9RAREVExLP6Mu3nz5qhfvz7eeOMNHDx40J41ERERUREsDu6bN2/igw8+wM2bN9GvXz/4+flh+PDh+PHHH5GTk2PPGomIiOgfFge3o6MjnnrqKXz99ddITk7Ghg0b4OPjg6lTp8LLywt9+vTBsmXLkJqaas96iYiIHmlWfx0MuPfznW3btsWcOXNw6tQpJCYmokOHDlixYgUCAwN5LXMiIiI7sfqs8sLUrl0bkydPxuTJk3Hz5k3cunXLFt0SERHRA6wO7h9//LHQ6ZIkwdHREbVr10bt2rXLXBgREREVVKrvcUuSBCGE2XTjNEmS0L59e2zcuBEeHh42K5SIiIhK8Rl3bGwsHn/8ccTGxiI9PR3p6emIjY1Fq1atsGXLFuzbtw83b97ElClT7FEvERHRI83qPe7x48fjq6++Qtu2bU3TunbtCkdHR4wcORK///475s2bh2HDhtm0UCIiIirFHve5c+fg7u5eYLq7uzv+/vtvAPdOVrtx40bZqyMiIiIzVgd3ixYt8Nprr+H69eumadevX8frr7+Oxx9/HABw9uxZVK9e3XZVkl09eL5CWds9Cjhm1uOYVTw+Bw8Hq4N76dKlOH/+PKpXr45atWqhdu3aqF69OpKSkvD1118DADIzMzF9+nSbF0u2t3v3bvj5+eGrr74qtt1XX30FPz8/7N69u5wqq7w4ZtbjmFU8PgcPEVEKBoNB/PTTT2L+/Pli3rx5IiYmRuj1+tJ0JRvp6ekCgEhPTy91H3l5eWLjxo0iLy/PhpWV3q5du4RGoxGSJAkAYvHixWbzjfV+8cUXAoCQJEloNBqxa9euCqq4eOUxviWNmdGiRYtKHLPKtj1YojQ123LMyqPeimSveu31HDyK42uLLCgrq4P74sWLRc47ePBgmYqpzB624Da+kBUKhQBgut3/gs7LyxOvvPKK2XyFQlFpw9ve42vJmAnx75tfSWNWmbYHS1lbs63HzN71VjR71GvP5+BRHF9ZBnfdunXFjRs3Ckw/cOCA0Gq1tqipUnqYgruoF/KDL2jjnvaDt8oa3vYcX0vH7ME3v+LGrLJsD9awpmZ7jJk9660MbF2vvZ+DR3F8K0NwS0JYdxbCiBEjkJCQgL1798LNzQ0AsG/fPjz11FOYOXMmJk6caE13spGRkQGtVov09PRCz6q3xN2kJMSvXoPHu3SBw62byLt+A4YbN6C/mwMHN1dI1avBkHoD+VcuQ5+VBUlSQKqihcrNDUp/fyg0jjCk3Ybk5QWFiwsUjhpArYYhOxtSfj6EkxNw5w4UajVyk5OhO3MWCmdnKD09oHR2hsHFBXtifsLApUuRbzDAUMxT/3SDBlh/6lSR8xWSBAelEmtHjkSnWrWhCAoGMtKhT8+AIecu4KCCUu0Ag9oRDlW0MFy9itzUVEjuWiizMgFHRzg0bAiRmQmFgxIGgwF5h49Cn5UJERAA5c1bUAYEAN5eUGRmQvL0hIOPN9KqVcN1APl6PRwdHaHX63H9+nVoNBo4Ozvj1KlTaNq0Kby8vCCEQHZ2NrKyspCTkwODwYD8/HzodDoA9344Jz8/H/n5+dD/05+7uzvy8/ORkZGB/Px8KJVKnD59Gq+//jry8/OLPWmna9eu2LVrV5HzJUmCg4MD5s+fjzZt2iArKwtnzpxBSEgIlEoltFotcnNzkZWVBWdnZ+j1ety+fRtCCHh5ecHT0xNZWVnIy8sz1Z6ZmQmDwQBJkgAALi4u0Gg00Ol00Ov1MOjzYbh8ErqLp2HQ66FxdYVP3WbIvHEN6ZfPQePiBpWzGyQhoHL3gErtBLWzKyQHDW6eTYRBr4cuJxNKBzWc/ash+3oqLp44Aq23NzwD68AjpB4UCgV0menIlxyQdfUccjPv4Ni5K3hzyQ/I1xuKHbMebUMRE59Q7HamVCqw8J3X0PKxahBKJRxUGgCAUuUApcYd2dcvAgolPAJrI+9uJqA3QOPpi8zkJGTeuIZzZ0+jfvM2cHLXQpeZAbWrFrmZt6HPy4GrXxCcfapBn5MJfVYm8oUEpYMElWsVCL0euXfSkJ93Fw4KB6jdPZFx9S/k3L4Ft2o14ezlCwlOMDhURU6aDhAGODiqIQwGSEoJuRk5yM/TwyNQC4WDAldPpSLvTh48g7QQQoG8zBzo8/RQqpXIyc6DZJDg28ATP8cdRMOaDeFSxQn6XAG1ixJCqUDGlTvIy86Fo6sTHByVuHU1A7o7uXDzd4OkAPKy9HCsooaDQgIUShxOjMeYmS8jX1/8dvtMn/74YdP3RT8HCgUclA5YuWg1Wjd5HDeS0qBxUcEz2AN3Uu/gcNwvaNikEVy1TnDQKJGbrYPSWQUnjQqSSgkJEpROCtz6+zZunL8FRw8n+NfyglKtxN20XBgMBrj7uUGSgLQrGdC4aVClqhv0OffG0cFRCV2OHkqVBEdXR+h1emTfvoubF9Og1+XjbkYO9DoDXLSOcPJygdbfFQqlEgpJID9XD4VKCbWzCjlpuTBIBsQnHEBkv0ioVKoi17k4tsiCsrI6uIUQePbZZ5GamoodO3bg4MGD6N27N959912MHz/eXnVWuLI8WVlbtiD9v+9BXLly7+pydqqxJEIINLuWjFsGA2xxzqgEwFOhQKJfgCk47EEAyHVxQUrdOrjaqBEy/Pygd9QAknTvZkdCCEyZMgWZmZk269PV1RUfffSRXccMej1Ul36F6soJKHIygPuecXtvf0IIvPzjX7iTq7dZn24aJRb1rmXfMbOSEBIMQot8QyDyDfWhN9QEoKzosgDcew5e2z4GmXl3bNanq9oNH0Z8UameA6tJAJRAUMvqePz5pvAKsv7qnpUhuK0+q1ySJKxZswaOjo7o2rUrevfujdmzZ5cqtGfPno3HH38cbm5u8PX1Rd++fXHmzBmzNkOHDoUkSWa3Nm3amLXJzc3FuHHj4O3tDRcXF/Tu3RuXL182a3P79m1ERUVBq9VCq9UiKioKaWlpVtdsrcx165A2eizE5ctABX/FQpIkvObmbpPQBu5FwWtu7nYP7RxXV5zp0hmXQkMhlAoo9Pl2D2wjSZLQp08fm/bZp08f+775CQHltT+gPn8Qipx04IFn3N5boSRJeK6Rj037fK6RT6UJDCHu3fQGHxiEFkrpOjTKn6FSHgNgu39WykKSJPSp19+mffap17/SPAelJgDkAxcOXcb29+OQ+qc8rzdiUXCfOHHC7Hb69GnMmDEDly5dwgsvvIAOHTqY5lkjLi4OY8aMwaFDhxAbG4v8/HyEh4cjKyvLrF2PHj2QnJxsum3bts1s/oQJE7BhwwasXbsWBw4cQGZmJiIjI6HX//siev7555GYmIiYmBjExMQgMTERUVFRVtVrLYNej/SZ7wD6yvFiBoAXXFwxR1vFJn3N0VbBCy6uNumrOJebNUWmjzc0WZmAQUDv6PjvP0Hl8M9Qhw4dMGjQIJv0NWjQIHTo0MEmfRUpNxtC4wqDYxX7Pk4xuj5WBcNb+Nmkr+Et/NH1sSo26ctWDMILQngCcICAIwAJKsVpKKTLJS1abp4I7oLnm7xok74GNXkRTwR3sUlflUXm9UwciU6EPt9Q0aVYzaJD5QqFosAPi9x///4fGNGXIaSuX78OX19fxMXFmd7chg4dirS0NGzcuLHQZdLT0+Hj44OVK1diwIABAICrV68iMDAQ27ZtQ0REBE6fPo0GDRrg0KFDaN26NQDg0KFDCAsLwx9//IG6deuWWFtpDo9krV6DtNffMAsXAfsfqrREdFYmpqanlXr58grtu66u+C2yFyR9PhR6PW5Xr44cd3dIBgOEwz9X7C2nvYB9+/Zh1apVpV6+XEJbCEi3LkJSqqC6mADV9bNFNi2PUdt17jaWHrtW6uUrW2gbX8r5hpoQUEGCDoABgAQJd5FvqIU8fQdUjlf5PfuSdmP1ieWlXv5hDG3j0+Oo1aDHm53gV8/yI0SV4VC5RdcqP3/+vL3rAHAvhAHA09PTbPrevXvh6+uLKlWqoGPHjnjvvffg6+sLADh27Bh0Oh3Cw8NN7atWrYpGjRohPj4eEREROHjwILRarSm0AaBNmzbQarWIj48vNLhzc3ORm5trup+RkQEA0Ol0ppObSnL32LEKPzxeFGPolia8yyu0ASDX1RX5GjWcbmcBkGD4J6wlg7D7Id8HGUO3NOFdLqH9D4Ved++UYAd1uTxecbo+du8zxNKEd2UL7X8pAThAMh0WlwAICDhAIaXjXpBXjs+6AaDDP6FbmvB+KEP7Pvm5ety+mg5PK7YzS9//7cmi4A4KCrJ3HRBCYNKkSWjfvj0aNWpkmt6zZ088++yzCAoKwvnz5zF9+nR06dIFx44dg0ajQUpKCtRqdYGfEPXz80NKSgoAICUlxRT09/P19TW1edDs2bMxa9asAtN37NgBZ2dni9ap6vUbCEDBvezKEuWDXFxxKC8XG+/etXiZfk7OGOTiWm7roNDrIekNEIp7n+pIxn+EKmiHpkOHDjh79iyOHDli8TKtWrUqt9AGcO8IhEIJSVSOQ4BdH/PA6et3EX8xw+Jl2tVwr6ShDdwLZoF7nzQa/wYkGGCAAyrT3rZRh+AuOHvzDxy9ctDiZVpVa/vwhrbxyIk+H8d/S8SZO79bvGh2dradirKcRcF98OBBhIWFWdRhVlYWkpKS0LBhQ6sKGTt2LE6cOIEDBw6YTTce/gaARo0aoWXLlggKCsLWrVvx9NNPF9mf8dC9UWEnVTzY5n5vvvkmJk2aZLqfkZGBwMBAhIeHW3x4RNe0GW7HxUHKz//3MVF5XtbRWZlWhTYAbLibjdZqdbntcbvcvAmntDRke1SBKicXquxs5Dk7w6CsmD2affv2WRXaAHDkyBHUrl273MLboHGHZMiDlJtVcuNysOvcbatCGwB+vpiBej7OlS68Jemf9w1kwQB3SKZ/YSUABugNgSjFOb92ty9pt1WhDQBHrsSjtlfdhzO8pXs3rY8bej7fBU5aR4sXNR59rUgWbWGDBw9G9+7d8d133xX5tZhTp07hrbfeQq1atZCQUPT3Mgszbtw4/Pjjj9izZ0+JP04SEBCAoKAgnD1777M7f39/5OXl4fbt22btUlNT4efnZ2pz7VrBQ3XXr183tXmQRqOBu7u72Q0AVCqVxTfn4CCoWreyaizKS1k+456anoboLNt9Pao4CoMBAadOQxJAjpsrnNIzoMzLAxTl/+ZYls+4V61ahX379tm4okJIEoSLJ6ScO1DeKfrwdHn981iWz7iXHkvBrnNpti3IRhTSjX8+31YD0N8LcuGFfEPNii6tgLJ8xr3qxHLsT3oIr1kuAAe1ErU71YS7t5tV7+ul/f63LVn07nfq1Cn06dMH//nPf+Dh4YGGDRuie/fueOqpp9C+fXt4e3ujRYsWuHDhAmJjYy0+W1sIgbFjx2L9+vXYvXs3QkJCSlzm5s2buHTpEgICAgDc+7UylUqF2NhYU5vk5GT89ttvpt8MDwsLQ3p6utme0uHDh5Genm72u+L24BW9Eg7NmpXbCVSWKOuJaUD5hbcEwPvvvxF09Bc43snE3SpVzPe2ZXJiGlCO4a1UwuDiDYOjG0QhL3E5hLZRZQzve5cPyINScQkK6QYAAb0IRG5+JwAuFVydubKemAY8nOGtclGhUa96aPa0dUeGKwurL8CSkJCA/fv3IykpCXfv3oW3tzeaN2+Ozp07FziprCSjR4/G6tWrsWnTJrMTxLRaLZycnJCZmYmZM2fimWeeQUBAAJKSkvDWW2/h4sWLOH36tOnKba+88gq2bNmCFStWwNPTE1OmTMHNmzdx7NgxKP95k+/ZsyeuXr2KxYsXAwBGjhyJoKAgbN682aJay3om4d2f45Exbx7unPkTTrm5QF7evZu1jBcdcXICVKp7fxsM/37lTAggXwfkFt63LUL7fuV1opoAcMfXFxdDQ5FStzZytdpiQ9v4nX+Doeyf89oitO9Xbieq5efDIfl3OFz9DYrM61CW49kVu86lYemxws8fKQ3LTlRT/HMIVAEY9LDsbBIFFOp7F/MR+ToYT3lUKB3u/aXPh9Ab/jn/zHxbklROkFQByM1tijxdfQihgAH3LrAk/nk5SpIEyUGCQW8A8lEy48fmxrsOCgDi3lUOrfzCzv6k3VhVxtC+35DWL6FttU4Q+f+MkUqCwXBfXcb/Ef95WSpVCihVSigdFNDl6JF/974BkAClWgFhuDdGSrXiXl8GAYVCAbWbChpXDRQOinvv4ZKAg0YFJzc1dLn5uHbmBu5m5BQYE8kBcNQ6wUmrhgIKCCGgVCshAdAbAKVKwh1VGp4cGQGv6tZffAWoHGeVl+rXwWwFhVwbF4BYvny5EEKI7OxsER4eLnx8fIRKpRI1atQQQ4YMKfBDJ3fv3hVjx44Vnp6ewsnJSURGRhZoc/PmTTFo0CDh5uYm3NzcxKBBg8Tt27ctrvVhuFa5wWAQPj4+pl8IKutNkiTh4+MjDAZDhazPg+wxvvYcs4reHkrDkpor03YmtzG2Vb3l9Rw8iuNbGa5VXqFnUYh7P3JS4DZ06FAAgJOTE7Zv347U1FTk5eXhwoULWLFiBQIDA836cXR0xIIFC3Dz5k1kZ2dj8+bNBdp4enoiOjoaGRkZyMjIQHR0NKpUqVJOa1o5SJKEtWvXQq1WQ1HCZ8QlXXBEoVBArVZj7dq18r+aUjE4ZtbjmFU8PgcPt8p3+iPZVZcuXbBt2zaoVKoiX9CLFy9GdHQ0vvjii0LnKxQKqFQqbNu2DV26PIRnnD7AmjFbtGhRofM5ZgVxzOyLz8FDrML29WXmYThUfj/+Hrf1+Hvc/D1ue+PvcdsXD5WTrN3/37jx8NfixYsxcuRIs3YRERGmPW9Jkh7p/74tHbNRo0aZ9mA4Zhyzisbn4OFjdXCX1+VPyf6ML2hvb+9CX8hGI0aMwOLFi+Ht7f3Iv5AtHbNRo0ZxzP7BMat4fA4eLhZdOe1+tWrVQocOHTB8+HD0798fjo6WX3GGKp8uXbrg2rVrJZ50MnLkSIwYMYInp4BjVhocs4rH5+DhYfUe9/Hjx9G8eXNMnjwZ/v7+GDVqlNWXgKTKxdIXKF/I/+KYWY9jVvH4HDwcrA7uRo0a4ZNPPsGVK1ewfPlypKSkoH379mjYsCE++eQTXL9+3R51EhEREcrwdTAHBwf069cP3333HebOnYtz585hypQpqF69OgYPHozk5GRb1klEREQoQ3D/8ssvGD16NAICAvDJJ59gypQpOHfuHHbv3o0rV66gT58+tqyTiIiIUIqT0z755BMsX74cZ86cwZNPPolvv/0WTz75pOkL/iEhIVi8eDHq1atn82KJiIgedVYH98KFCzFs2DC8+OKL8Pf3L7RNjRo1sHTp0jIXR0REROasDu7Y2FjUqFGjwCX0hBC4dOkSatSoAbVajSFDhtisSCIiIrrH6s+4H3vsMdy4caPA9Fu3bln0e9pERERUelYHtyji57szMzN5MRYiIiI7s/hQ+aRJkwDc+2L+f/7zHzg7O5vm6fV6HD58GM2aNbN5gURERPQvi4P7119/BXBvj/vkyZNQq9WmeWq1Gk2bNsWUKVNsXyERERGZWBzce/bsAQC8+OKLmD9/Ptzd3e1WFBERERXO6rPKly9fbo86iIiIyAIWBffTTz+NFStWwN3dHU8//XSxbdevX2+TwoiIiKggi4Jbq9Wafi1Gq9XatSAiIiIqmkXBff/hcR4qJyIiqjhWf4/77t27yM7ONt2/cOEC5s2bhx07dti0MCIiIirI6uDu06cPvv32WwBAWloaWrVqhY8//hh9+vTBwoULbV4gERER/cvq4E5ISMATTzwBAPj+++/h7++PCxcu4Ntvv8Vnn31m8wKJiIjoX1YHd3Z2Ntzc3AAAO3bswNNPPw2FQoE2bdrgwoULNi+QiIiI/mV1cNeqVQsbN27EpUuXsH37doSHhwMAUlNTeVEWIiIiO7M6uP/zn/9gypQpCA4ORuvWrREWFgbg3t538+bNbV4gERER/cvqK6f1798f7du3R3JyMpo2bWqa3rVrV/Tr18+mxREREZE5q4MbAPz9/eHv7282rVWrVjYpiIiIiIpmdXBnZWVhzpw52LVrF1JTU2EwGMzm//333zYrjoiIiMxZHdwvvfQS4uLiEBUVhYCAANOlUImIiMj+rA7un376CVu3bkW7du3sUQ8REREVw+qzyj08PODp6WmPWoiIiKgEVgf3f//7X/znP/8xu145ERERlQ+rD5V//PHHOHfuHPz8/BAcHAyVSmU2PyEhwWbFERERkTmrg7tv3752KIOIiIgsYXVwz5gxwx51EBERkQWs/owbuPdznl9//TXefPNN3Lp1C8C9Q+RXrlyxaXFERERkzuo97hMnTqBbt27QarVISkrCiBEj4OnpiQ0bNph+3pOIiIjsw+o97kmTJmHo0KE4e/YsHB0dTdN79uyJffv22bQ4IiIiMmd1cB89ehSjRo0qML1atWpISUmxSVFERERUOKuD29HRERkZGQWmnzlzBj4+PjYpioiIiApndXD36dMH77zzDnQ6HQBAkiRcvHgRU6dOxTPPPGPzAomIiOhfVgf3Rx99hOvXr8PX1xd3795Fx44dUatWLbi5ueG9996zR41ERET0D6vPKnd3d8eBAwewe/duJCQkwGAwIDQ0FN26dbNHfURERHQfq4PbqEuXLujSpYstayEiIqISWBXcBoMBK1aswPr165GUlARJkhASEoL+/fsjKiqKv81NRERkZxZ/xi2EQO/evfHSSy/hypUraNy4MRo2bIgLFy5g6NCh6Nevnz3rJCIiIlixx71ixQrs27cPu3btQufOnc3m7d69G3379sW3336LwYMH27xIIiIiusfiPe41a9bgrbfeKhDawL3Pu6dOnYpVq1bZtDgiIiIyZ3FwnzhxAj169Chyfs+ePXH8+HGbFEVERESFszi4b926BT8/vyLn+/n54fbt2zYpioiIiApncXDr9Xo4OBT9kbhSqUR+fr5NiiIiIqLCWXxymhACQ4cOhUajKXR+bm6uzYoiIiKiwlkc3EOGDCmxDc8oJyIisi+Lg3v58uX2rIOIiIgsYPWPjBAREVHFYXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZqdDgnj17Nh5//HG4ubnB19cXffv2xZkzZ4psP2rUKEiShHnz5plN79SpEyRJMrsNHDjQrM3t27cRFRUFrVYLrVaLqKgopKWl2WGtiIiI7KdCgzsuLg5jxozBoUOHEBsbi/z8fISHhyMrK6tA240bN+Lw4cOoWrVqoX2NGDECycnJptvixYvN5j///PNITExETEwMYmJikJiYiKioKLusFxERkb04VOSDx8TEmN1fvnw5fH19cezYMXTo0ME0/cqVKxg7diy2b9+OXr16FdqXs7Mz/P39C513+vRpxMTE4NChQ2jdujUAYMmSJQgLC8OZM2dQt25dG60RERGRfVVocD8oPT0dAODp6WmaZjAYEBUVhddeew0NGzYsctlVq1YhOjoafn5+6NmzJ2bMmAE3NzcAwMGDB6HVak2hDQBt2rSBVqtFfHx8ocGdm5uL3Nxc0/2MjAwAgE6ng06nK9X6GZcr7fLljfXal9zqBeRXM+u1r0ex3sqwrpUmuIUQmDRpEtq3b49GjRqZps+dOxcODg549dVXi1x20KBBCAkJgb+/P3777Te8+eabOH78OGJjYwEAKSkp8PX1LbCcr68vUlJSCu1z9uzZmDVrVoHpO3bsgLOzs7WrZ8ZYl1ywXvuSW72A/Gpmvfb1KNWbnZ1tw0pKp9IE99ixY3HixAkcOHDANO3YsWOYP38+EhISIElSkcuOGDHC9HejRo1Qu3ZttGzZEgkJCQgNDQWAQpcXQhTZ75tvvolJkyaZ7mdkZCAwMBDh4eFwd3e3ev2Ae/+pxcbGonv37lCpVKXqozyxXvuSW72A/Gpmvfb1KNZrPPpakSpFcI8bNw4//vgj9u3bh+rVq5um79+/H6mpqahRo4Zpml6vx+TJkzFv3jwkJSUV2l9oaChUKhXOnj2L0NBQ+Pv749q1awXaXb9+HX5+foX2odFooNFoCkxXqVRl3kBt0Ud5Yr32Jbd6AfnVzHrt61GqtzKsZ4UGtxAC48aNw4YNG7B3716EhISYzY+KikK3bt3MpkVERCAqKgovvvhikf3+/vvv0Ol0CAgIAACEhYUhPT0dR44cQatWrQAAhw8fRnp6Otq2bWvjtSIiIrKfCg3uMWPGYPXq1di0aRPc3NxMnzdrtVo4OTnBy8sLXl5eZsuoVCr4+/ubTig7d+4cVq1ahSeffBLe3t44deoUJk+ejObNm6Ndu3YAgPr166NHjx4YMWKE6WtiI0eORGRkJM8oJyIiWanQ73EvXLgQ6enp6NSpEwICAky3devWWdyHWq3Grl27EBERgbp16+LVV19FeHg4du7cCaVSaWq3atUqNG7cGOHh4QgPD0eTJk2wcuVKe6wWERGR3VT4oXJrPfi5dmBgIOLi4kpcztPTE9HR0VY/HhERUWXCa5UTERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iIiIZYXATERHJSIUG9+zZs/H444/Dzc0Nvr6+6Nu3L86cOVNk+1GjRkGSJMybN89sem5uLsaNGwdvb2+4uLigd+/euHz5slmb27dvIyoqClqtFlqtFlFRUUhLS7PDWhEREdlPhQZ3XFwcxowZg0OHDiE2Nhb5+fkIDw9HVlZWgbYbN27E4cOHUbVq1QLzJkyYgA0bNmDt2rU4cOAAMjMzERkZCb1eb2rz/PPPIzExETExMYiJiUFiYiKioqLsun5ERES25lCRDx4TE2N2f/ny5fD19cWxY8fQoUMH0/QrV65g7Nix2L59O3r16mW2THp6OpYuXYqVK1eiW7duAIDo6GgEBgZi586diIiIwOnTpxETE4NDhw6hdevWAIAlS5YgLCwMZ86cQd26de28pkRERLZRocH9oPT0dACAp6enaZrBYEBUVBRee+01NGzYsMAyx44dg06nQ3h4uGla1apV0ahRI8THxyMiIgIHDx6EVqs1hTYAtGnTBlqtFvHx8YUGd25uLnJzc033MzIyAAA6nQ46na5U62dcrrTLlzfWa19yqxeQX82s174exXorw7pWmuAWQmDSpElo3749GjVqZJo+d+5cODg44NVXXy10uZSUFKjVanh4eJhN9/PzQ0pKiqmNr69vgWV9fX1NbR40e/ZszJo1q8D0HTt2wNnZ2eL1KkxsbGyZli9vrNe+5FYvIL+aWa99PUr1Zmdn27CS0qk0wT127FicOHECBw4cME07duwY5s+fj4SEBEiSZFV/QgizZQpb/sE293vzzTcxadIk0/2MjAwEBgYiPDwc7u7uVtVipNPpEBsbi+7du0OlUpWqj/LEeu1LbvUC8quZ9drXo1iv8ehrRaoUwT1u3Dj8+OOP2LdvH6pXr26avn//fqSmpqJGjRqmaXq9HpMnT8a8efOQlJQEf39/5OXl4fbt22Z73ampqWjbti0AwN/fH9euXSvwuNevX4efn1+hNWk0Gmg0mgLTVSpVmTdQW/RRnlivfcmtXkB+NbNe+3qU6q0M61mhZ5ULITB27FisX78eu3fvRkhIiNn8qKgonDhxAomJiaZb1apV8dprr2H79u0AgBYtWkClUpkd+khOTsZvv/1mCu6wsDCkp6fjyJEjpjaHDx9Genq6qQ0REZEcVOge95gxY7B69Wps2rQJbm5ups+btVotnJyc4OXlBS8vL7NlVCoV/P39TSeUabVaDB8+HJMnT4aXlxc8PT0xZcoUNG7c2HSWef369dGjRw+MGDECixcvBgCMHDkSkZGRPKOciIhkpUL3uBcuXIj09HR06tQJAQEBptu6deus6ufTTz9F37598dxzz6Fdu3ZwdnbG5s2boVQqTW1WrVqFxo0bIzw8HOHh4WjSpAlWrlxp61UiIiKyqwrd4xZCWL1MUlJSgWmOjo5YsGABFixYUORynp6eiI6OtvrxiIiIKhNeq5yIiEhGGNxEREQywuAmIiKSEQY3ERGRjDC4iYiIZITBTUREJCMMbiIiIhlhcBMREckIg5uIiEhGGNxEREQywuAmIiKSEQY3ERGRjDC4iYiIZITBTUREJCMMbiIiIhlhcBMREckIg5uIiEhGGNxEREQywuAmIiKSEQY3ERGRjDC4iYiIZITBTUREJCMMbiIiIhlhcBMREckIg5uIiEhGGNxEREQywuAmIiKSEQY3ERGRjDC4iYiIZITBTUREJCMMbiIiIhlhcBMREckIg5uIiEhGGNxEREQywuAmIiKSEQY3ERGRjDC4iYiIZITBTUREJCMMbiIiIhlhcBMREckIg5uIiEhGGNxEREQywuAmIiKSEQY3ERGRjDC4iYiIZITBTUREJCMMbiIiIhlhcBMREckIg5uIiEhGGNxEREQywuAmIiKSEQY3ERGRjDC4iYiIZITBTUREJCMMbiIiIhlxqOgC5EIIAQDIyMgodR86nQ7Z2dnIyMiASqWyVWl2w3rtS271AvKrmfXa16NYrzEDjJlQERjcFrpz5w4AIDAwsIIrISKiinbnzh1otdoKeWxJVOS/DTJiMBhw9epVuLm5QZKkUvWRkZGBwMBAXLp0Ce7u7jau0PZYr33JrV5AfjWzXvt6FOsVQuDOnTuoWrUqFIqK+bSZe9wWUigUqF69uk36cnd3l8VGbsR67Utu9QLyq5n12tejVm9F7Wkb8eQ0IiIiGWFwExERyQiDuxxpNBrMmDEDGo2mokuxCOu1L7nVC8ivZtZrX6y3YvDkNCIiIhnhHjcREZGMMLiJiIhkhMFNREQkIwxuIiIiGWFw/+PLL79ESEgIHB0d0aJFC+zfv7/Y9nFxcWjRogUcHR1Rs2ZNLFq0qECbH374AQ0aNIBGo0GDBg2wYcOGIvubPXs2JEnChAkTzKYLITBz5kxUrVoVTk5O6NSpE37//XezekNDQ/HMM8/A29sbLi4u6N27Ny5fvlzu9ep0Orzxxhto3LgxXFxcULVqVQwePBhXr141q9d49bn7bwMHDqyQ8R06dGiBWtq0aQPAfJvw9fWFVqstcnzLs+YH6zXe+vXrV+5jPHPmzAKP4+/vb9amMm3DJdVb2bZhS8a3qG24Mo4vULm2XwC4cuUKXnjhBXh5ecHZ2RnNmjXDsWPHTPOL2n7vl5ubi3HjxhU7vjYlSKxdu1aoVCqxZMkScerUKTF+/Hjh4uIiLly4UGj7v//+Wzg7O4vx48eLU6dOiSVLlgiVSiW+//57U5v4+HihVCrF+++/L06fPi3ef/994eDgIA4dOlSgvyNHjojg4GDRpEkTMX78eLN5c+bMEW5ubuKHH34QJ0+eFAMGDBBVqlQxq7dx48ZCkiSxatUqkZCQIDp37iyaNm0q8vPzy7XetLQ00a1bN7Fu3Trxxx9/iIMHD4rWrVuLkJAQs3qrVasmHBwcxNGjR0VycrJITk4WaWlpFTK+Q4YMET169DDVkZycLG7evGm2TQwYMEC4uLgIR0dHsXXr1gLjW941319rcnKyWLZsmQBQIWM8Y8YM0bBhQ7N6UlNTzeqtTNtwSfVWtm3YkvEtbBv++uuvK+X4Vrbt99atWyIoKEgMHTpUHD58WJw/f17s3LlT/PXXX8VuvwEBASIjI8PU5uWXXxbVqlUTsbGxhY6vrTG4hRCtWrUSL7/8stm0evXqialTpxba/vXXXxf16tUzmzZq1CjRpk0b0/3nnntO9OjRw6xNRESEGDhwoNm0O3fuiNq1a4vY2FjRsWNHszdpg8Eg/P39xZw5c0zTcnJyhFKpFB06dBBC3HujUalUomrVqqZ6r1y5IhQKhYiJiSnXegtz5MgRAUC88MILpmkdO3YUHh4eFT6+Qtx70+vTp0+BGozbhHF8165da9omHhzf8q75QX369BHu7u5m23B5jfGMGTNE06ZNi6ytsm3DJdVbmIrchi2pt7Bt+P73tMo+vhW5/b7xxhuiffv2RdZW1Par1WrFokWLhBDC7D3CqLD3CFt65A+V5+Xl4dixYwgPDzebHh4ejvj4+EKXOXjwYIH2ERER+OWXX6DT6Ypt82CfY8aMQa9evdCtW7cCj3P+/HmkpKSY9SNJEgwGg+ni9seOHYNOp0NkZKSp76pVq6JRo0am++VVb2Fu3LgBAOjRo4fZ9Lt37+Kjjz5Cw4YNMWXKFNOvr1VEvXv37oWvry/q1KmDESNG4PLly6Ztwji+4eHhpm3iwfGtiJqNrl27hi1btiAzM7NA3+U1xmfPnkXVqlUREhKCgQMH4u+//zbNq4zbcHH1Fqait2FL6r1/Gx42bJjZe1plHt+K3n5//PFHtGzZEs8++yx8fX3RvHlzLFmyxDS/sO1Xo9GgY8eOpn7uf48wKuw9wpYe+eC+ceMG9Ho9/Pz8zKb7+fkhJSWl0GVSUlIKbZ+fn296kRfV5v4+165di4SEBMyePbvIxzEud3+9QghkZ2eb2qjVagQFBZn1ff9jlVe9D8rJycEbb7wBAAgJCTFNHzRoEP7v//4P1atXx/Tp0/HDDz/g6aefNlvv8qq3Z8+eWLVqFXbv3o2PP/4YR48eRdeuXU3bhHF8PTw8zPp/8LEqaoy/+eYbuLq6wmAwmPVdXmPcunVrfPvtt9i+fTuWLFmClJQUtG3bFjdv3jT1YVzOqCK34ZLqfVBFb8OW1PvgNnz48GHo9Xp4eHhU+vGt6O3377//xsKFC1G7dm1s374dL7/8Ml599VV8++23pj6MyxXVz/3vEcU9li3x18H+8eBPdQohiv35zsLaPzi9uD4vXbqE8ePHY8eOHXB0dLSqNkvqffB+edYL3DvJZ+DAgYU+zogRI5Camor4+HgMHDgQtWvXRsuWLZGQkIDQ0NByrXfAgAGmvxs1aoSWLVuiRo0aJfZf2PZR3mMMAMuWLUO/fv2wYsWKch9j4F5oGDVu3BhhYWF47LHH8M0332DSpElF9mNJ37behq2pF6j4bdjSeh/chqtXr47Q0FAcPHgQnTp1KrLvih5foOK3X4PBgJYtW+L9998HADRv3hy///47Fi5ciMGDB1vcT2EsaVNaj/wet7e3N5RKZYH/jFJTUwv8l2Xk7+9faHsHBwd4eXkV28bY57Fjx5CamooWLVrAwcEBDg4OiIuLw2effQYHBwfo9XrT2Zj39+Pt7Q1JkuDk5GR6nLy8PFy8eNGs3vsfq7zqNdLpdHjuuedw/vx57Nixo8TxDQ0NhUqlwtmzZyuk3vsFBAQgKCgICoUCKSkppvG9ffu2Wf8Pbh8VUfP+/ftx5swZjBkzpkLGuDAuLi5o3Lix2eMAlWMbtqReo8qwDVtT7/0aNmwIADh16pTpcSrb+AKVY/sNCAhAgwYNzNrUr18fFy9eNPUBoNh+7n+PKO6xbOmRD261Wo0WLVogNjbWbHpsbCzatm1b6DJhYWEF2u/YsQMtW7aESqUqto2xz65du+LkyZNITEw03Vq2bIlBgwYhMTERSqUSISEh8Pf3L9CPQqGAwWAAALRo0QIqlQpbt2419Z2cnIzffvvNdL+86gX+fcM7e/Ysdu7ciYCAgBLH9/fff4dOp0NAQEC51/ugmzdv4vLlywgODkZsbKxpfGNjY001Pzi+FVXz0qVL0aJFC7Rs2bJCxrgwubm5OH36tOlxKtM2bEm9QOXZhi2t90F37tyBJEm4evUqgMo3vkaVYftt164dzpw5Y9bmzz//RFBQEIDCt9+8vDzExcWZ+rn/PcKosPcIm7LLKW8yY/zqz9KlS8WpU6fEhAkThIuLi0hKShJCCDF16lQRFRVlam/8KsLEiRPFqVOnxNKlSwt8FeHnn38WSqVSzJkzR5w+fVrMmTOnyK/+GBV2BvGcOXOEVqsV69evFydPnhT/93//Z/oqjbHeB7/qERQUJDw8PAp81cPe9ep0OtG7d29RvXp1kZiYaPoax6JFi0z1xsTEiDZt2ggnJyexf/9+sXXrVuHp6Vkh9d65c0dMnjxZxMfHi/Pnz4s9e/aIsLAwUa1aNbF8+XJTzQMGDBCurq7C0dFRbNmyRXTp0kX4+vqanWVcntuEEEKkp6cLZ2dnsXDhQiGE+TZcnmM8efJksXfvXvH333+LQ4cOicjISOHm5mZ67QhRubbhkuqtbNtwSfUWtQ17eHhUyvGtbNvvkSNHhIODg3jvvffE2bNnxapVq4Szs7OIjo4udvst7Otg1atXFzt37hQJCQmiS5cu/DpYefjiiy9EUFCQUKvVIjQ0VMTFxZnmDRkyRHTs2NGs/d69e0Xz5s2FWq0WwcHBpg3wfv/73/9E3bp1hUqlEvXq1RM//PBDsTUU9iZtMBjEjBkzhL+/v9BoNKJDhw7i5MmTZvU2b95c9OvXT3h6egonJydRvXp1s69FlFe958+fFwAKvY0fP95Ur6urq3BzcxNqtVo89thjon79+qJdu3blXm92drYIDw8XPj4+QqVSiRo1aoghQ4aIixcvCiHMtwkfHx/h5uYmnJycRGRkpOjfv3+FbRNCCLF48WLh5ORk9t3W++strzE2fqfV+HWjp59+Wvz+++9mbSrTNlxSvZVtGy6p3uK24co4vkaVZfsVQojNmzeLRo0aCY1GI+rVqye++uors/lFbb/3u3v3rhg7dqxpfCMjI03vI/bAn/UkIiKSkUf+M24iIiI5YXATERHJCIObiIhIRhjcREREMsLgJiIikhEGNxERkYwwuImIiGSEwU1ERCQjDG4iGZAkCRs3bqzoMspNVFSU6RebSmvLli1o3ry56ZroRA8LBjdRBRk6dCgkSYIkSVCpVPDz80P37t2xbNmyAmGTnJxs9pOJFWXmzJlo1qyZXR/jxIkT2Lp1K8aNG1emfiIjIyFJElavXm2jyogqBwY3UQXq0aMHkpOTkZSUhJ9++gmdO3fG+PHjERkZifz8fFM7f39/aDSaCqzUtnQ6XZHzPv/8czz77LNwc3Mr8+O8+OKLWLBgQZn7IapMGNxEFUij0cDf3x/VqlVDaGgo3nrrLWzatAk//fQTVqxYYWr34KHyN954A3Xq1IGzszNq1qyJ6dOnm4Whcc942bJlqFGjBlxdXfHKK69Ar9fjgw8+gL+/P3x9ffHee++Z1ZOeno6RI0fC19cX7u7u6NKlC44fPw4AWLFiBWbNmoXjx4+bjhQYayxuuQfrqVmzJjQaDQr7mQSDwYD//e9/6N27t9n04OBgvPvuuxg8eDBcXV0RFBSETZs24fr16+jTpw9cXV3RuHFj/PLLL2bL9e7dG0eOHMHff/9t1fNCVJkxuIkqmS5duqBp06ZYv359kW3c3NywYsUKnDp1CvPnz8eSJUvw6aefmrU5d+4cfvrpJ8TExGDNmjVYtmwZevXqhcuXLyMuLg5z587FtGnTcOjQIQCAEAK9evVCSkoKtm3bhmPHjiE0NBRdu3bFrVu3MGDAAEyePBkNGzZEcnIykpOTMWDAgBKXM/rrr7/w3Xff4YcffkBiYmKh63XixAmkpaWhZcuWBeZ9+umnaNeuHX799Vf06tULUVFRGDx4MF544QUkJCSgVq1aGDx4sNk/BEFBQfD19cX+/futeQqIKje7/e4YERVryJAhok+fPoXOGzBggKhfv77pPgCxYcOGIvv64IMPRIsWLUz3Z8yYIZydnc1+MzgiIkIEBwcLvV5vmla3bl0xe/ZsIYQQu3btEu7u7iInJ8es78cee0wsXrzY1G/Tpk3N5lu6nEqlEqmpqUWugxBCbNiwQSiVSmEwGMymBwUFmf3+eXJysgAgpk+fbpp28OBBAUAkJyebLdu8eXMxc+bMYh+XSE4cKvj/BiIqhBACkiQVOf/777/HvHnz8NdffyEzMxP5+flwd3c3axMcHGz2ObGfnx+USiUUCoXZtNTUVADAsWPHkJmZCS8vL7N+7t69i3PnzhVZi6XLBQUFwcfHp5i1vreMRqMpdN2bNGliVjcANG7cuMC01NRU+Pv7m6Y7OTkhOzu72MclkhMGN1EldPr0aYSEhBQ679ChQxg4cCBmzZqFiIgIaLVarF27Fh9//LFZO5VKZXbfePb6g9OMZ7AbDAYEBARg7969BR6zSpUqRdZq6XIuLi5F9mHk7e2N7Oxs5OXlQa1Wm827v3ZjsBc27cEz8m/dulXiPwxEcsLgJqpkdu/ejZMnT2LixImFzv/5558RFBSEt99+2zTtwoULZX7c0NBQpKSkwMHBAcHBwYW2UavV0Ov1Vi9nKeNXzU6dOmWTr53l5OTg3LlzaN68eZn7IqoseHIaUQXKzc1FSkoKrly5goSEBLz//vvo06cPIiMjMXjw4EKXqVWrFi5evIi1a9fi3Llz+Oyzz7Bhw4Yy19KtWzeEhYWhb9++2L59O5KSkhAfH49p06aZztYODg7G+fPnkZiYiBs3biA3N9ei5Szl4+OD0NBQHDhwoMzrA9w7OqHRaBAWFmaT/ogqAwY3UQWKiYlBQEAAgoOD0aNHD+zZswefffYZNm3aBKVSWegyffr0wcSJEzF27Fg0a9YM8fHxmD59eplrkSQJ27ZtQ4cOHTBs2DDUqVMHAwcORFJSkunz42eeeQY9evRA586d4ePjgzVr1li0nDVGjhyJVatWlXl9AGDNmjUYNGgQnJ2dbdIfUWUgCVHIlymJiCpITk4O6tati7Vr15ZpT/n69euoV68efvnllyLPFyCSI+5xE1Gl4ujoiG+//RY3btwoUz/nz5/Hl19+ydCmhw73uImIiGSEe9xEREQywuAmIiKSEQY3ERGRjDC4iYiIZITBTUREJCMMbiIiIhlhcBMREckIg5uIiEhGGNxEREQy8v+pNK6WugfGGAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ ">> Phases: \n", " Diameters: [0.00423757 0.00566325 0.00510815 0.00467504]\n", " Densities: [2500. 2500. 2500. 2500.]\n" ] } ], "source": [ "CG.get_particle_phases(Diameter_t0, Density_t0, GlobalID_t0, 8)\n", "print(\">> Phases: \")\n", "print(\" Diameters: \", CG.phases[:,0])\n", "print(\" Densities: \", CG.phases[:,1])" ] }, { "cell_type": "markdown", "id": "aa881edc", "metadata": {}, "source": [ "### 4. Calculate the coarse-graining resolution\n", "\n", "The argument is the reference diameter; a larger value gives a wider kernel (larger `w` and cutoff `c`)." ] }, { "cell_type": "code", "execution_count": 9, "id": "805306e4", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ ">> Coarse Graining resolution: \n", " c: 0.007601418599101063\n", " w: 0.0038007092995505313\n" ] } ], "source": [ "CG.set_resolution(CG.d43)\n", "print(\">> Coarse Graining resolution: \")\n", "print(\" c:\", CG.c)\n", "print(\" w:\", CG.w)" ] }, { "cell_type": "markdown", "id": "166b974d", "metadata": {}, "source": [ "### 5. Generate the coarse-graining grid" ] }, { "cell_type": "code", "execution_count": 10, "id": "08dc944d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Generating Grid with Customised Grid Ranges\n", "Customised grid bounds: x = [-0.23, 0.23], y = [0.001, 0.1], z = [None, None], x_transect = None, y_transect = None, z_transect = 0.0\n", ">> Grid: \n", " Grid Points: (3444, 3) First Point: [-0.23 0.001 0. ]\n", " Nodes: [123 28 1]\n", " Spacing: [0.00377049 0.00366667]\n" ] } ], "source": [ "CG.generate_grid()\n", "print(\">> Grid: \")\n", "print(\" Grid Points: \", CG.GridPoints.shape, \"First Point: \", CG.GridPoints[0])\n", "print(\" Nodes: \", CG.Nodes)\n", "print(\" Spacing: \", CG.Spacing)" ] }, { "cell_type": "markdown", "id": "40015e17", "metadata": {}, "source": [ "### 6. Calculate the coarse-grained fields" ] }, { "cell_type": "code", "execution_count": 11, "id": "ebe31ec1", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \n", "-------------------- Calculating Coarse Grained Fields --------------------\n", " \n", "---> Time step 0: time 0600 ................................................\n", "Loading data ... \n", " Particle data loaded\n", " Repeated pairs in contact data: 0\n", " Contact data loaded and mapped\n", " Coordination number not provided. Calculating it.\n", "... data loaded\n", "Matching particles to grid points ...\n", "... particles assigned to grid nodes\n", "Computing weights ...\n", " Using Lucy kernel with cutoff 0.007601418599101063 and search sampling factor 1000\n", "... weights computed\n", "Computing Coarse Graining fields...\n", " volume fraction done\n", " mixture density done\n", " momentum density done\n", " velocity done\n", " velocity gradient done\n", " kinetic tensor done\n", " contact tensor done\n", " Z done\n", " frabric tensor done\n", " total stress done\n", " pressure done\n", " granular temp done\n", " mu done\n", "... fields computed\n", "Writing results for timestep 600...\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/exports/csce/datastore/geos/users/s1857688/Coarse_Graining/pysammos/pysammos/macroscopic_fields/gridded/derived_fields.py:216: RuntimeWarning: divide by zero encountered in divide\n", " out[f\"frictional_coefficient_{s_name}_P{p_name}\"] = tau / out[f\"pressure_{p_name}\"]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " File successfully updated to ./PysammosCG/CG_Lucy_Polydisperse.h5\n", " Warning: VTKHDFWriter promoting 2D ImageData (nz=1) to thin 3D (nz=2) for VTK compatibility.\n", "File successfully written to ./PysammosCG/CG_Lucy_Polydisperse_0600.vtkhdf\n", "... results written\n", " File successfully updated to ./PysammosCG/CG_GranularTemperature_slices.h5\n", " File successfully written to ./PysammosCG/CG_GranularTemperature_slices.h5\n", ">> time step 0 took 9.809610605239868 to run.\n", " \n", "---> Time step 1: time 0601 ................................................\n", "Loading data ... \n", " Particle data loaded\n", " Repeated pairs in contact data: 0\n", " Contact data loaded and mapped\n", " Coordination number not provided. Calculating it.\n", "... data loaded\n", "Matching particles to grid points ...\n", "... particles assigned to grid nodes\n", "Computing weights ...\n", " Using Lucy kernel with cutoff 0.007601418599101063 and search sampling factor 1000\n", "... weights computed\n", "Computing Coarse Graining fields...\n", " volume fraction done\n", " mixture density done\n", " momentum density done\n", " velocity done\n", " velocity gradient done\n", " kinetic tensor done\n", " contact tensor done\n", " Z done\n", " frabric tensor done\n", " total stress done\n", " pressure done\n", " granular temp done\n", " mu done\n", "... fields computed\n", "Writing results for timestep 601...\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/exports/csce/datastore/geos/users/s1857688/Coarse_Graining/pysammos/pysammos/macroscopic_fields/gridded/derived_fields.py:216: RuntimeWarning: divide by zero encountered in divide\n", " out[f\"frictional_coefficient_{s_name}_P{p_name}\"] = tau / out[f\"pressure_{p_name}\"]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " File successfully updated to ./PysammosCG/CG_Lucy_Polydisperse.h5\n", " Warning: VTKHDFWriter promoting 2D ImageData (nz=1) to thin 3D (nz=2) for VTK compatibility.\n", "File successfully written to ./PysammosCG/CG_Lucy_Polydisperse_0601.vtkhdf\n", "... results written\n", " File successfully updated to ./PysammosCG/CG_GranularTemperature_slices.h5\n", " File successfully written to ./PysammosCG/CG_GranularTemperature_slices.h5\n", ">> time step 1 took 0.6379914283752441 to run.\n", " \n", "---> Time step 2: time 0602 ................................................\n", "Loading data ... \n", " Particle data loaded\n", " Repeated pairs in contact data: 0\n", " Contact data loaded and mapped\n", " Coordination number not provided. Calculating it.\n", "... data loaded\n", "Matching particles to grid points ...\n", "... particles assigned to grid nodes\n", "Computing weights ...\n", " Using Lucy kernel with cutoff 0.007601418599101063 and search sampling factor 1000\n", "... weights computed\n", "Computing Coarse Graining fields...\n", " volume fraction done\n", " mixture density done\n", " momentum density done\n", " velocity done\n", " velocity gradient done\n", " kinetic tensor done\n", " contact tensor done\n", " Z done\n", " frabric tensor done\n", " total stress done\n", " pressure done\n", " granular temp done\n", " mu done\n", "... fields computed\n", "Writing results for timestep 602...\n", " File successfully updated to ./PysammosCG/CG_Lucy_Polydisperse.h5\n", " Warning: VTKHDFWriter promoting 2D ImageData (nz=1) to thin 3D (nz=2) for VTK compatibility.\n", "File successfully written to ./PysammosCG/CG_Lucy_Polydisperse_0602.vtkhdf\n", "... results written\n", " File successfully updated to ./PysammosCG/CG_GranularTemperature_slices.h5\n", " File successfully written to ./PysammosCG/CG_GranularTemperature_slices.h5\n", ">> time step 2 took 0.7271218299865723 to run.\n", " \n", "---> Time step 3: time 0603 ................................................\n", "Loading data ... \n", " Particle data loaded\n", " Repeated pairs in contact data: 0\n", " Contact data loaded and mapped\n", " Coordination number not provided. Calculating it.\n", "... data loaded\n", "Matching particles to grid points ...\n", "... particles assigned to grid nodes\n", "Computing weights ...\n", " Using Lucy kernel with cutoff 0.007601418599101063 and search sampling factor 1000\n", "... weights computed\n", "Computing Coarse Graining fields...\n", " volume fraction done\n", " mixture density done\n", " momentum density done\n", " velocity done\n", " velocity gradient done\n", " kinetic tensor done\n", " contact tensor done\n", " Z done\n", " frabric tensor done\n", " total stress done\n", " pressure done\n", " granular temp done\n", " mu done\n", "... fields computed\n", "Writing results for timestep 603...\n", " File successfully updated to ./PysammosCG/CG_Lucy_Polydisperse.h5\n", " Warning: VTKHDFWriter promoting 2D ImageData (nz=1) to thin 3D (nz=2) for VTK compatibility.\n", "File successfully written to ./PysammosCG/CG_Lucy_Polydisperse_0603.vtkhdf\n", "... results written\n", " File successfully updated to ./PysammosCG/CG_GranularTemperature_slices.h5\n", " File successfully written to ./PysammosCG/CG_GranularTemperature_slices.h5\n", ">> time step 3 took 0.6878252029418945 to run.\n", " \n", "---> Time step 4: time 0604 ................................................\n", "Loading data ... \n", " Particle data loaded\n", " Repeated pairs in contact data: 0\n", " Contact data loaded and mapped\n", " Coordination number not provided. Calculating it.\n", "... data loaded\n", "Matching particles to grid points ...\n", "... particles assigned to grid nodes\n", "Computing weights ...\n", " Using Lucy kernel with cutoff 0.007601418599101063 and search sampling factor 1000\n", "... weights computed\n", "Computing Coarse Graining fields...\n", " volume fraction done\n", " mixture density done\n", " momentum density done\n", " velocity done\n", " velocity gradient done\n", " kinetic tensor done\n", " contact tensor done\n", " Z done\n", " frabric tensor done\n", " total stress done\n", " pressure done\n", " granular temp done\n", " mu done\n", "... fields computed\n", "Writing results for timestep 604...\n", " File successfully updated to ./PysammosCG/CG_Lucy_Polydisperse.h5\n", " Warning: VTKHDFWriter promoting 2D ImageData (nz=1) to thin 3D (nz=2) for VTK compatibility.\n", "File successfully written to ./PysammosCG/CG_Lucy_Polydisperse_0604.vtkhdf\n", "... results written\n", " File successfully updated to ./PysammosCG/CG_GranularTemperature_slices.h5\n", " File successfully written to ./PysammosCG/CG_GranularTemperature_slices.h5\n", ">> time step 4 took 0.7081248760223389 to run.\n", " \n", "---> Time step 5: time 0605 ................................................\n", "Loading data ... \n", " Particle data loaded\n", " Repeated pairs in contact data: 0\n", " Contact data loaded and mapped\n", " Coordination number not provided. Calculating it.\n", "... data loaded\n", "Matching particles to grid points ...\n", "... particles assigned to grid nodes\n", "Computing weights ...\n", " Using Lucy kernel with cutoff 0.007601418599101063 and search sampling factor 1000\n", "... weights computed\n", "Computing Coarse Graining fields...\n", " volume fraction done\n", " mixture density done\n", " momentum density done\n", " velocity done\n", " velocity gradient done\n", " kinetic tensor done\n", " contact tensor done\n", " Z done\n", " frabric tensor done\n", " total stress done\n", " pressure done\n", " granular temp done\n", " mu done\n", "... fields computed\n", "Writing results for timestep 605...\n", " File successfully updated to ./PysammosCG/CG_Lucy_Polydisperse.h5\n", " Warning: VTKHDFWriter promoting 2D ImageData (nz=1) to thin 3D (nz=2) for VTK compatibility.\n", "File successfully written to ./PysammosCG/CG_Lucy_Polydisperse_0605.vtkhdf\n", "... results written\n", " File successfully updated to ./PysammosCG/CG_GranularTemperature_slices.h5\n", " File successfully written to ./PysammosCG/CG_GranularTemperature_slices.h5\n", ">> time step 5 took 0.6912028789520264 to run.\n", " \n", "---> Time step 6: time 0606 ................................................\n", "Loading data ... \n", " Particle data loaded\n", " Repeated pairs in contact data: 0\n", " Contact data loaded and mapped\n", " Coordination number not provided. Calculating it.\n", "... data loaded\n", "Matching particles to grid points ...\n", "... particles assigned to grid nodes\n", "Computing weights ...\n", " Using Lucy kernel with cutoff 0.007601418599101063 and search sampling factor 1000\n", "... weights computed\n", "Computing Coarse Graining fields...\n", " volume fraction done\n", " mixture density done\n", " momentum density done\n", " velocity done\n", " velocity gradient done\n", " kinetic tensor done\n", " contact tensor done\n", " Z done\n", " frabric tensor done\n", " total stress done\n", " pressure done\n", " granular temp done\n", " mu done\n", "... fields computed\n", "Writing results for timestep 606...\n", " File successfully updated to ./PysammosCG/CG_Lucy_Polydisperse.h5\n", " Warning: VTKHDFWriter promoting 2D ImageData (nz=1) to thin 3D (nz=2) for VTK compatibility.\n", "File successfully written to ./PysammosCG/CG_Lucy_Polydisperse_0606.vtkhdf\n", "... results written\n", " File successfully updated to ./PysammosCG/CG_GranularTemperature_slices.h5\n", " File successfully written to ./PysammosCG/CG_GranularTemperature_slices.h5\n", ">> time step 6 took 0.673731803894043 to run.\n", " \n", "---> Time step 7: time 0607 ................................................\n", "Loading data ... \n", " Particle data loaded\n", " Repeated pairs in contact data: 0\n", " Contact data loaded and mapped\n", " Coordination number not provided. Calculating it.\n", "... data loaded\n", "Matching particles to grid points ...\n", "... particles assigned to grid nodes\n", "Computing weights ...\n", " Using Lucy kernel with cutoff 0.007601418599101063 and search sampling factor 1000\n", "... weights computed\n", "Computing Coarse Graining fields...\n", " volume fraction done\n", " mixture density done\n", " momentum density done\n", " velocity done\n", " velocity gradient done\n", " kinetic tensor done\n", " contact tensor done\n", " Z done\n", " frabric tensor done\n", " total stress done\n", " pressure done\n", " granular temp done\n", " mu done\n", "... fields computed\n", "Writing results for timestep 607...\n", " File successfully updated to ./PysammosCG/CG_Lucy_Polydisperse.h5\n", " Warning: VTKHDFWriter promoting 2D ImageData (nz=1) to thin 3D (nz=2) for VTK compatibility.\n", "File successfully written to ./PysammosCG/CG_Lucy_Polydisperse_0607.vtkhdf\n", "... results written\n", " File successfully updated to ./PysammosCG/CG_GranularTemperature_slices.h5\n", " File successfully written to ./PysammosCG/CG_GranularTemperature_slices.h5\n", ">> time step 7 took 0.72312331199646 to run.\n", " \n", "---> Time step 8: time 0608 ................................................\n", "Loading data ... \n", " Particle data loaded\n", " Repeated pairs in contact data: 0\n", " Contact data loaded and mapped\n", " Coordination number not provided. Calculating it.\n", "... data loaded\n", "Matching particles to grid points ...\n", "... particles assigned to grid nodes\n", "Computing weights ...\n", " Using Lucy kernel with cutoff 0.007601418599101063 and search sampling factor 1000\n", "... weights computed\n", "Computing Coarse Graining fields...\n", " volume fraction done\n", " mixture density done\n", " momentum density done\n", " velocity done\n", " velocity gradient done\n", " kinetic tensor done\n", " contact tensor done\n", " Z done\n", " frabric tensor done\n", " total stress done\n", " pressure done\n", " granular temp done\n", " mu done\n", "... fields computed\n", "Writing results for timestep 608...\n", " File successfully updated to ./PysammosCG/CG_Lucy_Polydisperse.h5\n", " Warning: VTKHDFWriter promoting 2D ImageData (nz=1) to thin 3D (nz=2) for VTK compatibility.\n", "File successfully written to ./PysammosCG/CG_Lucy_Polydisperse_0608.vtkhdf\n", "... results written\n", " File successfully updated to ./PysammosCG/CG_GranularTemperature_slices.h5\n", " File successfully written to ./PysammosCG/CG_GranularTemperature_slices.h5\n", ">> time step 8 took 0.7497551441192627 to run.\n", " \n", "---> Time step 9: time 0609 ................................................\n", "Loading data ... \n", " Particle data loaded\n", " Repeated pairs in contact data: 0\n", " Contact data loaded and mapped\n", " Coordination number not provided. Calculating it.\n", "... data loaded\n", "Matching particles to grid points ...\n", "... particles assigned to grid nodes\n", "Computing weights ...\n", " Using Lucy kernel with cutoff 0.007601418599101063 and search sampling factor 1000\n", "... weights computed\n", "Computing Coarse Graining fields...\n", " volume fraction done\n", " mixture density done\n", " momentum density done\n", " velocity done\n", " velocity gradient done\n", " kinetic tensor done\n", " contact tensor done\n", " Z done\n", " frabric tensor done\n", " total stress done\n", " pressure done\n", " granular temp done\n", " mu done\n", "... fields computed\n", "Writing results for timestep 609...\n", " File successfully updated to ./PysammosCG/CG_Lucy_Polydisperse.h5\n", " Warning: VTKHDFWriter promoting 2D ImageData (nz=1) to thin 3D (nz=2) for VTK compatibility.\n", "File successfully written to ./PysammosCG/CG_Lucy_Polydisperse_0609.vtkhdf\n", "... results written\n", " File successfully updated to ./PysammosCG/CG_GranularTemperature_slices.h5\n", " File successfully written to ./PysammosCG/CG_GranularTemperature_slices.h5\n", ">> time step 9 took 0.6694228649139404 to run.\n", " \n", "---> Time step 10: time 0610 ................................................\n", "Loading data ... \n", " Particle data loaded\n", " Repeated pairs in contact data: 0\n", " Contact data loaded and mapped\n", " Coordination number not provided. Calculating it.\n", "... data loaded\n", "Matching particles to grid points ...\n", "... particles assigned to grid nodes\n", "Computing weights ...\n", " Using Lucy kernel with cutoff 0.007601418599101063 and search sampling factor 1000\n", "... weights computed\n", "Computing Coarse Graining fields...\n", " volume fraction done\n", " mixture density done\n", " momentum density done\n", " velocity done\n", " velocity gradient done\n", " kinetic tensor done\n", " contact tensor done\n", " Z done\n", " frabric tensor done\n", " total stress done\n", " pressure done\n", " granular temp done\n", " mu done\n", "... fields computed\n", "Writing results for timestep 610...\n", " File successfully updated to ./PysammosCG/CG_Lucy_Polydisperse.h5\n", " Warning: VTKHDFWriter promoting 2D ImageData (nz=1) to thin 3D (nz=2) for VTK compatibility.\n", "File successfully written to ./PysammosCG/CG_Lucy_Polydisperse_0610.vtkhdf\n", "... results written\n", " File successfully updated to ./PysammosCG/CG_GranularTemperature_slices.h5\n", " File successfully written to ./PysammosCG/CG_GranularTemperature_slices.h5\n", ">> time step 10 took 0.7532615661621094 to run.\n", " \n", "---> Time step 11: time 0611 ................................................\n", "Loading data ... \n", " Particle data loaded\n", " Repeated pairs in contact data: 0\n", " Contact data loaded and mapped\n", " Coordination number not provided. Calculating it.\n", "... data loaded\n", "Matching particles to grid points ...\n", "... particles assigned to grid nodes\n", "Computing weights ...\n", " Using Lucy kernel with cutoff 0.007601418599101063 and search sampling factor 1000\n", "... weights computed\n", "Computing Coarse Graining fields...\n", " volume fraction done\n", " mixture density done\n", " momentum density done\n", " velocity done\n", " velocity gradient done\n", " kinetic tensor done\n", " contact tensor done\n", " Z done\n", " frabric tensor done\n", " total stress done\n", " pressure done\n", " granular temp done\n", " mu done\n", "... fields computed\n", "Writing results for timestep 611...\n", " File successfully updated to ./PysammosCG/CG_Lucy_Polydisperse.h5\n", " Warning: VTKHDFWriter promoting 2D ImageData (nz=1) to thin 3D (nz=2) for VTK compatibility.\n", "File successfully written to ./PysammosCG/CG_Lucy_Polydisperse_0611.vtkhdf\n", "... results written\n", " File successfully updated to ./PysammosCG/CG_GranularTemperature_slices.h5\n", " File successfully written to ./PysammosCG/CG_GranularTemperature_slices.h5\n", ">> time step 11 took 0.6785123348236084 to run.\n", " \n", "---> Time step 12: time 0612 ................................................\n", "Loading data ... \n", " Particle data loaded\n", " Repeated pairs in contact data: 0\n", " Contact data loaded and mapped\n", " Coordination number not provided. Calculating it.\n", "... data loaded\n", "Matching particles to grid points ...\n", "... particles assigned to grid nodes\n", "Computing weights ...\n", " Using Lucy kernel with cutoff 0.007601418599101063 and search sampling factor 1000\n", "... weights computed\n", "Computing Coarse Graining fields...\n", " volume fraction done\n", " mixture density done\n", " momentum density done\n", " velocity done\n", " velocity gradient done\n", " kinetic tensor done\n", " contact tensor done\n", " Z done\n", " frabric tensor done\n", " total stress done\n", " pressure done\n", " granular temp done\n", " mu done\n", "... fields computed\n", "Writing results for timestep 612...\n", " File successfully updated to ./PysammosCG/CG_Lucy_Polydisperse.h5\n", " Warning: VTKHDFWriter promoting 2D ImageData (nz=1) to thin 3D (nz=2) for VTK compatibility.\n", "File successfully written to ./PysammosCG/CG_Lucy_Polydisperse_0612.vtkhdf\n", "... results written\n", " File successfully updated to ./PysammosCG/CG_GranularTemperature_slices.h5\n", " File successfully written to ./PysammosCG/CG_GranularTemperature_slices.h5\n", ">> time step 12 took 0.798590898513794 to run.\n", " \n", "---> Time step 13: time 0613 ................................................\n", "Loading data ... \n", " Particle data loaded\n", " Repeated pairs in contact data: 0\n", " Contact data loaded and mapped\n", " Coordination number not provided. Calculating it.\n", "... data loaded\n", "Matching particles to grid points ...\n", "... particles assigned to grid nodes\n", "Computing weights ...\n", " Using Lucy kernel with cutoff 0.007601418599101063 and search sampling factor 1000\n", "... weights computed\n", "Computing Coarse Graining fields...\n", " volume fraction done\n", " mixture density done\n", " momentum density done\n", " velocity done\n", " velocity gradient done\n", " kinetic tensor done\n", " contact tensor done\n", " Z done\n", " frabric tensor done\n", " total stress done\n", " pressure done\n", " granular temp done\n", " mu done\n", "... fields computed\n", "Writing results for timestep 613...\n", " File successfully updated to ./PysammosCG/CG_Lucy_Polydisperse.h5\n", " Warning: VTKHDFWriter promoting 2D ImageData (nz=1) to thin 3D (nz=2) for VTK compatibility.\n", "File successfully written to ./PysammosCG/CG_Lucy_Polydisperse_0613.vtkhdf\n", "... results written\n", " File successfully updated to ./PysammosCG/CG_GranularTemperature_slices.h5\n", " File successfully written to ./PysammosCG/CG_GranularTemperature_slices.h5\n", ">> time step 13 took 0.7418897151947021 to run.\n", " \n", "---> Time step 14: time 0614 ................................................\n", "Loading data ... \n", " Particle data loaded\n", " Repeated pairs in contact data: 0\n", " Contact data loaded and mapped\n", " Coordination number not provided. Calculating it.\n", "... data loaded\n", "Matching particles to grid points ...\n", "... particles assigned to grid nodes\n", "Computing weights ...\n", " Using Lucy kernel with cutoff 0.007601418599101063 and search sampling factor 1000\n", "... weights computed\n", "Computing Coarse Graining fields...\n", " volume fraction done\n", " mixture density done\n", " momentum density done\n", " velocity done\n", " velocity gradient done\n", " kinetic tensor done\n", " contact tensor done\n", " Z done\n", " frabric tensor done\n", " total stress done\n", " pressure done\n", " granular temp done\n", " mu done\n", "... fields computed\n", "Writing results for timestep 614...\n", " File successfully updated to ./PysammosCG/CG_Lucy_Polydisperse.h5\n", " Warning: VTKHDFWriter promoting 2D ImageData (nz=1) to thin 3D (nz=2) for VTK compatibility.\n", "File successfully written to ./PysammosCG/CG_Lucy_Polydisperse_0614.vtkhdf\n", "... results written\n", " File successfully updated to ./PysammosCG/CG_GranularTemperature_slices.h5\n", " File successfully written to ./PysammosCG/CG_GranularTemperature_slices.h5\n", ">> time step 14 took 0.7545323371887207 to run.\n", " \n" ] } ], "source": [ "CG.fields_in_time()" ] }, { "cell_type": "markdown", "id": "ec928561", "metadata": {}, "source": [ "## Plotting the results\n", "\n", "The coarse-grained fields are saved in `examples/jannsen_pile/PysammosCG/`." ] }, { "cell_type": "markdown", "id": "1d494717", "metadata": {}, "source": [ "### 1. Read the coarse-grained data" ] }, { "cell_type": "code", "execution_count": 12, "id": "2df026c8", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Processing key: coordination_number\n", "Processing key: density_mixture\n", "Processing key: fabric_tensor\n", "Processing key: frictional_coefficient_01_Pxy\n", "Processing key: frictional_coefficient_01_Pxyz\n", "Processing key: frictional_coefficient_01_Py\n", "Processing key: frictional_coefficient_Dxy_Pxy\n", "Processing key: frictional_coefficient_Dxy_Pxyz\n", "Processing key: frictional_coefficient_Dxy_Py\n", "Processing key: frictional_coefficient_Dxyz_Pxy\n", "Processing key: frictional_coefficient_Dxyz_Pxyz\n", "Processing key: frictional_coefficient_Dxyz_Py\n", "Processing key: granular_temperature_x\n", "Processing key: granular_temperature_xyz\n", "Processing key: granular_temperature_y\n", "Processing key: granular_temperature_z\n", "Processing key: phases\n", "Processing key: positions\n", "Processing key: pressure_x\n", "Processing key: pressure_xy\n", "Processing key: pressure_xyz\n", "Processing key: pressure_y\n", "Processing key: pressure_z\n", "Processing key: time\n", "Processing key: total_stress_tensor_xy_dev\n", "Processing key: total_stress_tensor_xy_dev_mag\n", "Processing key: total_stress_tensor_xyz\n", "Processing key: total_stress_tensor_xyz_dev\n", "Processing key: total_stress_tensor_xyz_dev_mag\n", "Processing key: total_stress_tensor_xyz_mag\n", "Processing key: velocity\n", "Processing key: volume_fraction\n" ] } ], "source": [ "manager = H5XarrayManager(\"./PysammosCG/CG_Lucy_Polydisperse.h5\")\n", "janssen_CG = manager.h5_to_xarray()" ] }, { "cell_type": "markdown", "id": "00990725", "metadata": {}, "source": [ "In polydisperse mode the fields have a `phase` dimension: index 0 is the bulk, followed by one entry per phase, in the order of `CG.phases`." ] }, { "cell_type": "code", "execution_count": 13, "id": "1c0d7e54", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'phase' (phase: 5)> Size: 700B\n",
       "array(['Bulk', 'Phase_[4.2375745e-03 2.5000000e+03]',\n",
       "       'Phase_[5.6632482e-03 2.5000000e+03]',\n",
       "       'Phase_[5.108148e-03 2.500000e+03]',\n",
       "       'Phase_[4.6750368e-03 2.5000000e+03]'], dtype='<U35')\n",
       "Coordinates:\n",
       "  * phase    (phase) <U35 700B 'Bulk' ... 'Phase_[4.6750368e-03 2.5000000e+03]'
" ], "text/plain": [ " Size: 700B\n", "array(['Bulk', 'Phase_[4.2375745e-03 2.5000000e+03]',\n", " 'Phase_[5.6632482e-03 2.5000000e+03]',\n", " 'Phase_[5.108148e-03 2.500000e+03]',\n", " 'Phase_[4.6750368e-03 2.5000000e+03]'], dtype='" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# coarse-grained data: time average of pressure_y, shape (points, phases)\n", "x = np.asarray(janssen_CG.positions[:, 0])\n", "y = np.asarray(janssen_CG.positions[:, 1])\n", "p_y = np.asarray(janssen_CG.pressure_y.values).mean(axis=0)\n", "p_y_bulk = p_y[:, 0]\n", "\n", "# DEM particles of the last time step\n", "last_step = int(janssen_CG.time[-1])\n", "janssen_DEM = pv.read(f\"../data_examples/jannsen_pile/DEM_data/DES_FB1_{last_step:04d}.vtp\")\n", "\n", "fig, ax = plt.subplots(1, 2, figsize=(13, 5), gridspec_kw={\"width_ratios\": [3, 1]})\n", "\n", "# (a) DEM particles and bulk pressure on the grid\n", "particles = ax[0].scatter(janssen_DEM.points[:, 0], janssen_DEM.points[:, 1], s=5,\n", " c=janssen_DEM[\"Diameter\"], cmap=\"bone_r\", zorder=3, label=\"DEM particles\")\n", "cax = ax[0].inset_axes([0.02, 0.93, 0.4, 0.035])\n", "cbar = fig.colorbar(particles, cax=cax, orientation=\"horizontal\")\n", "cbar.set_label(\"Particle diameter [m]\", fontsize=10)\n", "cbar.ax.tick_params(labelsize=9)\n", "\n", "log_p = np.log10(np.where(p_y_bulk > 0, p_y_bulk, np.nan)) # empty nodes have p = 0\n", "sc = ax[0].scatter(x, y, s=65, c=log_p, cmap=\"viridis\", marker=\"s\", zorder=1,\n", " vmin=-0.5, label=\"Pysammos CG\")\n", "cbar_sc = fig.colorbar(sc, ax=ax[0], pad=0.02, extend=\"min\")\n", "cbar_sc.set_label(r\"log$_{10}$(P$_{y,\\mathrm{bulk}}$) [Pa]\")\n", "\n", "x_center, tolerance = 0.0, 0.004\n", "ax[0].axvline(x=x_center, color=\"blue\", ls=\"--\", lw=2, label=\"Profile location\")\n", "\n", "# (b) vertical profiles along x = x_center: bulk and phases\n", "bulk_profile, y_profile = vertical_profile(x, y, p_y_bulk, x_center, tolerance)\n", "ax[1].plot(bulk_profile, y_profile, \"--\", color=\"blue\", lw=2, label=\"Bulk\")\n", "phase_colors = [\"#FFE686\", \"#BD5340\", \"orange\", \"#B3E9FF\", \"#6A8D73\", \"#8E7DBE\", \"#C9A227\", \"#555555\"]\n", "for p, d_phase in enumerate(CG.phases[:, 0], start=1):\n", " profile, _ = vertical_profile(x, y, p_y[:, p], x_center, tolerance)\n", " ax[1].plot(profile, y_profile, \"-\", lw=2, color=phase_colors[(p - 1) % len(phase_colors)],\n", " label=f\"d = {d_phase * 1e3:.1f} mm\")\n", "\n", "# hydrostatic reference: bulk density and fill height set by hand for this pile\n", "rho_bulk, fill_height, g = 800.0, 0.11, 9.81\n", "h = np.linspace(0, fill_height)\n", "ax[1].plot(rho_bulk * g * (fill_height - h), h, \"k--\", label=r\"$\\rho g h$\")\n", "\n", "for a, label in zip(ax, [\"(a)\", \"(b)\"]):\n", " a.set_title(label, loc=\"left\", fontweight=\"bold\")\n", " a.set_ylim(0, 0.15)\n", " a.grid()\n", "ax[0].legend(fontsize=10, loc=\"upper right\")\n", "ax[1].legend(fontsize=9, loc=\"upper right\")\n", "ax[0].set_xlabel(\"Distance [m]\")\n", "ax[0].set_ylabel(\"Height [m]\")\n", "ax[1].set_xlabel(r\"P$_y$ [Pa]\")\n", "ax[1].set_xlim(-1, 900)\n", "plt.subplots_adjust(wspace=0.15)\n", "plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Pysammos_Env", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.11" } }, "nbformat": 4, "nbformat_minor": 5 }