diff --git a/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/ssd-cube.ipynb b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/ssd-cube.ipynb
index b77bc7d1d80..624ed8ad900 100644
--- a/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/ssd-cube.ipynb
+++ b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/ssd-cube.ipynb
@@ -5,12 +5,13 @@
"id": "34c87f77-b604-4200-b102-da8290ef81b3",
"metadata": {},
"source": [
+ "+++\n",
"title = \"SteadyStateDiffusion Cube Test\"\n",
"date = \"2021-11-09\"\n",
"author = \"Lars Bilke\"\n",
"web_subsection = \"elliptic\"\n",
"draft = true\n",
- "\n"
+ "+++\n"
]
},
{
diff --git a/Tests/Data/HydroMechanics/SeabedResponse/Stationary_waves.ipynb b/Tests/Data/HydroMechanics/SeabedResponse/Stationary_waves.ipynb
index a7d6399627e..ded6fcfdbc6 100644
--- a/Tests/Data/HydroMechanics/SeabedResponse/Stationary_waves.ipynb
+++ b/Tests/Data/HydroMechanics/SeabedResponse/Stationary_waves.ipynb
@@ -5,11 +5,12 @@
"id": "09d7a481-d07d-465f-8f3b-a098b650295d",
"metadata": {},
"source": [
+ "+++\n",
"title = \"Seabed response to water waves\"\n",
"date = \"2023-01-11\"\n",
"author = \"Linda Günther\"\n",
"web_subsection = \"hydro-mechanics\"\n",
- ""
+ "+++\n"
]
},
{
@@ -139,7 +140,7 @@
"\\end{align}\n",
"$$\n",
"\n",
- "where $p$ is the pore pressure, $\\tilde{p}$ is the amplitude of the applied load, $\\sigma'_{yy}$ is the effective vertical stress and $\\sigma'_{xy}$ is the effective shear stress. The boundary condition of the pore pressure describes the space- and time-dependent water wave. Compared to Arnold Verruijt's solution, in this example there is a phase shift of $-\\frac{\\pi}{2}$ in the time-dependent part. The phase shift is necessary to obtain a water wave that starts oszillating from the equilibrium state (sine instead of cosine) and thus to be able to set an initial condition of $p=0$ Pa on the whole domain in the numerical solution. \n",
+ "where $p$ is the pore pressure, $\\tilde{p}$ is the amplitude of the applied load, $\\sigma'_{yy}$ is the effective vertical stress and $\\sigma'_{xy}$ is the effective shear stress. The boundary condition of the pore pressure describes the space- and time-dependent water wave. Compared to Arnold Verruijt's solution, in this example there is a phase shift of $-\\frac{\\pi}{2}$ in the time-dependent part. The phase shift is necessary to obtain a water wave that starts oszillating from the equilibrium state (sine instead of cosine) and thus to be able to set an initial condition of $p=0$ Pa on the whole domain in the numerical solution.\n",
"\n",
"With these boundary conditions, the following four constants are determined:\n",
"\n",
@@ -180,13 +181,15 @@
"import numpy as np\n",
"\n",
"import matplotlib.pyplot as plt\n",
- "plt.rc ('font', size = 8)\n",
- "plt.rc ('axes', titlesize = 10)\n",
- "plt.rc ('axes', labelsize = 10)\n",
+ "\n",
+ "plt.rc(\"font\", size=8)\n",
+ "plt.rc(\"axes\", titlesize=10)\n",
+ "plt.rc(\"axes\", labelsize=10)\n",
"\n",
"import gmsh\n",
"\n",
"import pyvista as pv\n",
+ "\n",
"pv.set_plot_theme(\"document\")\n",
"pv.set_jupyter_backend(\"static\")"
]
@@ -200,37 +203,90 @@
},
"outputs": [],
"source": [
- "def compute_pressure_and_stresses(t,x,z):\n",
- " \n",
- " n=0.4\n",
- " G=100e3 # [Pa]\n",
- " K=2/3*G # [Pa] (with ny=0)\n",
- " ny=0 # E = 3K(1-2ny) = 2G(1+ny)\n",
- " Cf=0 # in the book: Cf = 0.001/K\n",
- " Cs=0\n",
- " Cm=1/K\n",
- " my=1.3e-3 # [Pa*s]\n",
- " kappa=1e-11 # [m²] (medium sand, kf=10e-4 m/s) \n",
- " gamma_w = 9.81e3 # [Pa/m]\n",
- " lam=2*np.pi*0.1*10/100\n",
- " omega=2*np.pi*0.1\n",
- " \n",
- " k = kappa*gamma_w/my # Gl. (1.33)\n",
- " alpha = 1-Cs/Cm # Gl. (4.15)\n",
- " S = n*Cf + (alpha-n)*Cs # Gl. (1.28)\n",
- " theta = S*G/alpha**2 # Gl. (4.13)\n",
- " m = 1/(1-2*ny) # = K+1/3*G/G # Gl. (4.5)\n",
- " cv = k*G*(1+m) / (alpha**2*(1+theta+m*theta)*gamma_w) # Gl. (4.12)\n",
- " xi_2 = complex(lam**2, (omega/cv)) # Gl. (4.19)\n",
- " \n",
- " B1 = (1+m)*(xi_2-lam**2)-2*lam*(np.sqrt(xi_2)-lam)\n",
- " B2 = 2*m*theta*lam*np.sqrt(xi_2)+theta*((1+m)*(xi_2-lam**2)-2*lam*(np.sqrt(xi_2)-lam))\n",
- " B3 = 2*m*theta*lam\n",
- " D = 2*lam*(2*lam*(np.sqrt(xi_2)-lam)-(1+m)*(1+m*theta)*(xi_2-lam**2))\n",
- " p_rel = np.real((-2*lam*B1*np.exp(-lam*z) - (1+m)*(xi_2-lam**2)*B3*np.exp(-np.sqrt(xi_2)*z))/D * np.exp((omega*t-np.pi*0.5)*1j)*np.cos(lam*x))\n",
- " sig_xx_rel = np.real(((-2*(m-1)*lam*theta + 2*lam*(1+m*theta)*lam*z)*B1*np.exp(-lam*z) - 2*lam*B2*np.exp(-lam*z) + ((m-1)*(xi_2-lam**2) - 2*lam**2)*B3*np.exp(-np.sqrt(xi_2)*z))/D * np.exp((omega*t-np.pi*0.5)*1j)*np.cos(lam*x))\n",
- " sig_zz_rel = np.real(((-2*(m+1)*lam*theta - 2*lam*(1+m*theta)*lam*z)*B1*np.exp(-lam*z) + 2*lam*B2*np.exp(-lam*z) + ((m-1)*(xi_2-lam**2) + 2*xi_2)*B3*np.exp(-np.sqrt(xi_2)*z))/D * np.exp((omega*t-np.pi*0.5)*1j)*np.cos(lam*x))\n",
- " sig_xz_rel = np.real(((-2*lam*(1+m*theta)*lam*z-2*lam*theta)*B1*np.exp(-lam*z) + 2*lam*B2*np.exp(-lam*z) + 2*np.sqrt(xi_2)*lam*B3*np.exp(-np.sqrt(xi_2)*z))/D * np.exp((omega*t-np.pi*0.5)*1j)*np.sin(lam*x))\n",
+ "def compute_pressure_and_stresses(t, x, z):\n",
+ " n = 0.4\n",
+ " G = 100e3 # [Pa]\n",
+ " K = 2 / 3 * G # [Pa] (with ny=0)\n",
+ " ny = 0 # E = 3K(1-2ny) = 2G(1+ny)\n",
+ " Cf = 0 # in the book: Cf = 0.001/K\n",
+ " Cs = 0\n",
+ " Cm = 1 / K\n",
+ " my = 1.3e-3 # [Pa*s]\n",
+ " kappa = 1e-11 # [m²] (medium sand, kf=10e-4 m/s)\n",
+ " gamma_w = 9.81e3 # [Pa/m]\n",
+ " lam = 2 * np.pi * 0.1 * 10 / 100\n",
+ " omega = 2 * np.pi * 0.1\n",
+ "\n",
+ " k = kappa * gamma_w / my # Gl. (1.33)\n",
+ " alpha = 1 - Cs / Cm # Gl. (4.15)\n",
+ " S = n * Cf + (alpha - n) * Cs # Gl. (1.28)\n",
+ " theta = S * G / alpha**2 # Gl. (4.13)\n",
+ " m = 1 / (1 - 2 * ny) # = K+1/3*G/G # Gl. (4.5)\n",
+ " cv = (\n",
+ " k * G * (1 + m) / (alpha**2 * (1 + theta + m * theta) * gamma_w)\n",
+ " ) # Gl. (4.12)\n",
+ " xi_2 = complex(lam**2, (omega / cv)) # Gl. (4.19)\n",
+ "\n",
+ " B1 = (1 + m) * (xi_2 - lam**2) - 2 * lam * (np.sqrt(xi_2) - lam)\n",
+ " B2 = 2 * m * theta * lam * np.sqrt(xi_2) + theta * (\n",
+ " (1 + m) * (xi_2 - lam**2) - 2 * lam * (np.sqrt(xi_2) - lam)\n",
+ " )\n",
+ " B3 = 2 * m * theta * lam\n",
+ " D = (\n",
+ " 2\n",
+ " * lam\n",
+ " * (\n",
+ " 2 * lam * (np.sqrt(xi_2) - lam)\n",
+ " - (1 + m) * (1 + m * theta) * (xi_2 - lam**2)\n",
+ " )\n",
+ " )\n",
+ " p_rel = np.real(\n",
+ " (\n",
+ " -2 * lam * B1 * np.exp(-lam * z)\n",
+ " - (1 + m) * (xi_2 - lam**2) * B3 * np.exp(-np.sqrt(xi_2) * z)\n",
+ " )\n",
+ " / D\n",
+ " * np.exp((omega * t - np.pi * 0.5) * 1j)\n",
+ " * np.cos(lam * x)\n",
+ " )\n",
+ " sig_xx_rel = np.real(\n",
+ " (\n",
+ " (-2 * (m - 1) * lam * theta + 2 * lam * (1 + m * theta) * lam * z)\n",
+ " * B1\n",
+ " * np.exp(-lam * z)\n",
+ " - 2 * lam * B2 * np.exp(-lam * z)\n",
+ " + ((m - 1) * (xi_2 - lam**2) - 2 * lam**2)\n",
+ " * B3\n",
+ " * np.exp(-np.sqrt(xi_2) * z)\n",
+ " )\n",
+ " / D\n",
+ " * np.exp((omega * t - np.pi * 0.5) * 1j)\n",
+ " * np.cos(lam * x)\n",
+ " )\n",
+ " sig_zz_rel = np.real(\n",
+ " (\n",
+ " (-2 * (m + 1) * lam * theta - 2 * lam * (1 + m * theta) * lam * z)\n",
+ " * B1\n",
+ " * np.exp(-lam * z)\n",
+ " + 2 * lam * B2 * np.exp(-lam * z)\n",
+ " + ((m - 1) * (xi_2 - lam**2) + 2 * xi_2) * B3 * np.exp(-np.sqrt(xi_2) * z)\n",
+ " )\n",
+ " / D\n",
+ " * np.exp((omega * t - np.pi * 0.5) * 1j)\n",
+ " * np.cos(lam * x)\n",
+ " )\n",
+ " sig_xz_rel = np.real(\n",
+ " (\n",
+ " (-2 * lam * (1 + m * theta) * lam * z - 2 * lam * theta)\n",
+ " * B1\n",
+ " * np.exp(-lam * z)\n",
+ " + 2 * lam * B2 * np.exp(-lam * z)\n",
+ " + 2 * np.sqrt(xi_2) * lam * B3 * np.exp(-np.sqrt(xi_2) * z)\n",
+ " )\n",
+ " / D\n",
+ " * np.exp((omega * t - np.pi * 0.5) * 1j)\n",
+ " * np.sin(lam * x)\n",
+ " )\n",
" return p_rel, sig_xx_rel, sig_zz_rel, sig_xz_rel"
]
},
@@ -241,7 +297,7 @@
"source": [
"By evaluating these equations at different times $t$ and depths $y$, we gain a better understanding of the pressure and stress distribution in the seabed. The below plot illustrates the pore pressure and the amplitude of the effective stresses as a function of depth directly underneath an anti-node of the standing water wave (the place where the amplitude is at maximum, i.e. for $x = k \\cdot \\frac{L}{2}$, where $k=0, 1, 2,$ ...).\n",
"\n",
- "Along the top edge of the seabed, the pore pressure is always as large as the applied load and the effective stresses are zero. This means, that all the change in pressure is absorbed by the fluid while the soil particles remain in their initial stress state (in this case zero, since body forces are being disregarded). The increased pore pressure at the top edge cannot propagate freely downwards into the seabed because seepage is limited by the hydraulic conductivity of the soil. Consequently, the pore pressure decreases with depth as the soil matrix gradually takes up the remaining share of the total stress in the seabed. "
+ "Along the top edge of the seabed, the pore pressure is always as large as the applied load and the effective stresses are zero. This means, that all the change in pressure is absorbed by the fluid while the soil particles remain in their initial stress state (in this case zero, since body forces are being disregarded). The increased pore pressure at the top edge cannot propagate freely downwards into the seabed because seepage is limited by the hydraulic conductivity of the soil. Consequently, the pore pressure decreases with depth as the soil matrix gradually takes up the remaining share of the total stress in the seabed."
]
},
{
@@ -249,6 +305,7 @@
"execution_count": 3,
"id": "4f48e224-330e-4f3f-89aa-734981c8fdd4",
"metadata": {
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [
@@ -264,29 +321,56 @@
}
],
"source": [
- "y = np.linspace(0,100,1000)\n",
- "y_rel = y/100\n",
- "colors = {0:\"orangered\", 2:\"gold\", 4:\"blueviolet\", 6:\"forestgreen\", 8:\"darkorange\", 10:\"royalblue\"}\n",
- "\n",
- "fig, ax = plt.subplots(ncols=2, figsize=(15,7))\n",
- "for idx in (0,1):\n",
+ "y = np.linspace(0, 100, 1000)\n",
+ "y_rel = y / 100\n",
+ "colors = {\n",
+ " 0: \"orangered\",\n",
+ " 2: \"gold\",\n",
+ " 4: \"blueviolet\",\n",
+ " 6: \"forestgreen\",\n",
+ " 8: \"darkorange\",\n",
+ " 10: \"royalblue\",\n",
+ "}\n",
+ "\n",
+ "fig, ax = plt.subplots(ncols=2, figsize=(15, 7))\n",
+ "for idx in (0, 1):\n",
" ax[idx].grid(True)\n",
" ax[idx].set_ylabel(\"$y$ / $L$\")\n",
- " ax[idx].set_xlim(-1.1,1.1)\n",
+ " ax[idx].set_xlim(-1.1, 1.1)\n",
"\n",
- "for t in [0,2,4,6,8,10]:\n",
- " ax[0].plot(compute_pressure_and_stresses(t,0,y)[0], -y_rel, color=colors[t], label= \"t = %.1f s\" %t)\n",
+ "for t in [0, 2, 4, 6, 8, 10]:\n",
+ " ax[0].plot(\n",
+ " compute_pressure_and_stresses(t, 0, y)[0],\n",
+ " -y_rel,\n",
+ " color=colors[t],\n",
+ " label=\"t = %.1f s\" % t,\n",
+ " )\n",
"\n",
- "t=2.5\n",
+ "t = 2.5\n",
"ax[0].set_xlabel(\"$p$ / $\\\\tilde{p}$\")\n",
"ax[0].legend()\n",
- "ax[1].plot(compute_pressure_and_stresses(t,0,y)[1], -y_rel, color = colors[6], label = r\"$\\sigma'_{xx}/(\\alpha\\tilde{p})$\")\n",
- "#ax[1].plot(compute_pressure_and_stresses(t,0,y)[1]+compute_pressure_and_stresses(t,0,y)[0], -y_rel, linestyle = \"--\", color = colors[3], label = \"$\\\\sigma_{xx}$/$\\\\alpha\\\\tilde{p}$\") # Total horizontal stress\n",
- "ax[1].plot(compute_pressure_and_stresses(t,0,y)[2], -y_rel, color = colors[2], label = r\"$\\sigma'_{yy}/(\\alpha\\tilde{p})$\")\n",
- "#ax[1].plot(compute_pressure_and_stresses(t,0,y)[2]+compute_pressure_and_stresses(t,0,y)[0], -y_rel, linestyle = \"--\", color = colors[1], label = \"$\\\\sigma_{yy}$/$\\\\alpha\\\\tilde{p}$\") # Total vertical stress\n",
- "ax[1].plot(compute_pressure_and_stresses(t,0,y)[3], -y_rel, color = colors[4], label = r\"$\\sigma'_{xy}/(\\alpha\\tilde{p})$\")\n",
+ "ax[1].plot(\n",
+ " compute_pressure_and_stresses(t, 0, y)[1],\n",
+ " -y_rel,\n",
+ " color=colors[6],\n",
+ " label=r\"$\\sigma'_{xx}/(\\alpha\\tilde{p})$\",\n",
+ ")\n",
+ "# ax[1].plot(compute_pressure_and_stresses(t,0,y)[1]+compute_pressure_and_stresses(t,0,y)[0], -y_rel, linestyle = \"--\", color = colors[3], label = \"$\\\\sigma_{xx}$/$\\\\alpha\\\\tilde{p}$\") # Total horizontal stress\n",
+ "ax[1].plot(\n",
+ " compute_pressure_and_stresses(t, 0, y)[2],\n",
+ " -y_rel,\n",
+ " color=colors[2],\n",
+ " label=r\"$\\sigma'_{yy}/(\\alpha\\tilde{p})$\",\n",
+ ")\n",
+ "# ax[1].plot(compute_pressure_and_stresses(t,0,y)[2]+compute_pressure_and_stresses(t,0,y)[0], -y_rel, linestyle = \"--\", color = colors[1], label = \"$\\\\sigma_{yy}$/$\\\\alpha\\\\tilde{p}$\") # Total vertical stress\n",
+ "ax[1].plot(\n",
+ " compute_pressure_and_stresses(t, 0, y)[3],\n",
+ " -y_rel,\n",
+ " color=colors[4],\n",
+ " label=r\"$\\sigma'_{xy}/(\\alpha\\tilde{p})$\",\n",
+ ")\n",
"ax[1].set_xlabel(r\"$\\sigma'/(\\alpha\\tilde{p})$\")\n",
- "ax[1].legend();"
+ "ax[1].legend()"
]
},
{
@@ -302,6 +386,7 @@
"execution_count": 4,
"id": "9cc304d2-bd0a-443f-98dc-9c8488989d61",
"metadata": {
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [
@@ -317,28 +402,43 @@
}
],
"source": [
- "t=np.linspace(0,20,200)\n",
- "colors = {1:\"gold\", 2:\"blueviolet\", 3:\"forestgreen\", 4:\"royalblue\"}\n",
+ "t = np.linspace(0, 20, 200)\n",
+ "colors = {1: \"gold\", 2: \"blueviolet\", 3: \"forestgreen\", 4: \"royalblue\"}\n",
"\n",
- "fig, ax = plt.subplots(ncols=2, figsize=(15,7))\n",
- "for idx in (0,1):\n",
+ "fig, ax = plt.subplots(ncols=2, figsize=(15, 7))\n",
+ "for idx in (0, 1):\n",
" ax[idx].grid(True)\n",
" ax[idx].set_xlabel(\"$t$ / s\")\n",
"\n",
- "for y in (np.linspace(0,100,6)):\n",
- " ax[0].plot( t, compute_pressure_and_stresses(t,0,y)[0], color = colors[4])\n",
+ "for y in np.linspace(0, 100, 6):\n",
+ " ax[0].plot(t, compute_pressure_and_stresses(t, 0, y)[0], color=colors[4])\n",
" ax[0].set_ylabel(\"$p/\\\\tilde{p}$\")\n",
- " ax[1].plot(t, compute_pressure_and_stresses(t,0,y)[1], color = colors[3], label = r\"$\\sigma'_{xx}/(\\alpha\\tilde{p})$\")\n",
- " ax[1].plot(t, compute_pressure_and_stresses(t,0,y)[2], color = colors[1], label = r\"$\\sigma'_{yy}/(\\alpha\\tilde{p})$\")\n",
- " ax[1].plot(t, compute_pressure_and_stresses(t,0,y)[3], color = colors[2], label = r\"$\\sigma'_{xy}/(\\alpha\\tilde{p})$\")\n",
+ " ax[1].plot(\n",
+ " t,\n",
+ " compute_pressure_and_stresses(t, 0, y)[1],\n",
+ " color=colors[3],\n",
+ " label=r\"$\\sigma'_{xx}/(\\alpha\\tilde{p})$\",\n",
+ " )\n",
+ " ax[1].plot(\n",
+ " t,\n",
+ " compute_pressure_and_stresses(t, 0, y)[2],\n",
+ " color=colors[1],\n",
+ " label=r\"$\\sigma'_{yy}/(\\alpha\\tilde{p})$\",\n",
+ " )\n",
+ " ax[1].plot(\n",
+ " t,\n",
+ " compute_pressure_and_stresses(t, 0, y)[3],\n",
+ " color=colors[2],\n",
+ " label=r\"$\\sigma'_{xy}/(\\alpha\\tilde{p})$\",\n",
+ " )\n",
" if y == 0:\n",
" ax[1].legend(loc=\"upper right\")\n",
"\n",
"ax[1].set_ylabel(\"$\\sigma$'/$\\\\alpha\\\\tilde{p}$\")\n",
"\n",
- " \n",
+ "\n",
"ax[0].set_title(\"Pore pressure over time\")\n",
- "ax[1].set_title(\"Effective stresses over time\");"
+ "ax[1].set_title(\"Effective stresses over time\")"
]
},
{
@@ -354,6 +454,7 @@
"execution_count": 5,
"id": "739678d6-b1a4-4d0d-94e8-1774a2eb5b17",
"metadata": {
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [
@@ -369,28 +470,28 @@
}
],
"source": [
- "x, y = np.meshgrid(np.linspace(0,200,1000),np.linspace(0,100,1000))\n",
+ "x, y = np.meshgrid(np.linspace(0, 200, 1000), np.linspace(0, 100, 1000))\n",
"t = 2.5\n",
"\n",
- "fig, ax = plt.subplots(ncols=2, nrows=2, figsize=(15,7))\n",
- "l1=ax[0][0].contourf(x,-y, compute_pressure_and_stresses(t,x,y)[0], 15)\n",
- "l2=ax[0][1].contourf(x,-y, compute_pressure_and_stresses(t,x,y)[1], 15)\n",
- "l3=ax[1][1].contourf(x,-y, compute_pressure_and_stresses(t,x,y)[2], 15)\n",
- "l4=ax[1][0].contourf(x,-y, compute_pressure_and_stresses(t,x,y)[3], 15)\n",
- "fig.colorbar(l1,ax=ax[0][0])\n",
- "fig.colorbar(l2,ax=ax[0][1])\n",
- "fig.colorbar(l3,ax=ax[1][1])\n",
- "fig.colorbar(l4,ax=ax[1][0])\n",
- "for i in (0,1):\n",
- " for j in (0,1):\n",
- " ax[i][j].set_aspect('equal')\n",
- " ax[i][j].set_xlabel('$x$ / m')\n",
- " ax[i][j].set_ylabel('$y$ / m')\n",
+ "fig, ax = plt.subplots(ncols=2, nrows=2, figsize=(15, 7))\n",
+ "l1 = ax[0][0].contourf(x, -y, compute_pressure_and_stresses(t, x, y)[0], 15)\n",
+ "l2 = ax[0][1].contourf(x, -y, compute_pressure_and_stresses(t, x, y)[1], 15)\n",
+ "l3 = ax[1][1].contourf(x, -y, compute_pressure_and_stresses(t, x, y)[2], 15)\n",
+ "l4 = ax[1][0].contourf(x, -y, compute_pressure_and_stresses(t, x, y)[3], 15)\n",
+ "fig.colorbar(l1, ax=ax[0][0])\n",
+ "fig.colorbar(l2, ax=ax[0][1])\n",
+ "fig.colorbar(l3, ax=ax[1][1])\n",
+ "fig.colorbar(l4, ax=ax[1][0])\n",
+ "for i in (0, 1):\n",
+ " for j in (0, 1):\n",
+ " ax[i][j].set_aspect(\"equal\")\n",
+ " ax[i][j].set_xlabel(\"$x$ / m\")\n",
+ " ax[i][j].set_ylabel(\"$y$ / m\")\n",
"ax[0][0].set_title(\"$p/\\\\tilde{p}$\")\n",
"ax[0][1].set_title(\"$\\\\sigma'_{xx}/\\\\alpha\\\\tilde{p}$\")\n",
"ax[1][1].set_title(\"$\\\\sigma'_{yy}/\\\\alpha\\\\tilde{p}$\")\n",
"ax[1][0].set_title(\"$\\\\sigma'_{xy}/\\\\alpha\\\\tilde{p}$\")\n",
- "fig.tight_layout();"
+ "fig.tight_layout()"
]
},
{
@@ -436,7 +537,7 @@
"import os\n",
"\n",
"# out_dir will contain all data we produce\n",
- "out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', '_out')\n",
+ "out_dir = os.environ.get(\"OGS_TESTRUNNER_OUT_DIR\", \"_out\")\n",
"os.makedirs(out_dir, exist_ok=True)"
]
},
@@ -449,10 +550,10 @@
},
"outputs": [],
"source": [
- "def generate_mesh_axb(a,b,Nx,Ny,P):\n",
+ "def generate_mesh_axb(a, b, Nx, Ny, P):\n",
" output_file = f\"{out_dir}/square_{a}x{b}.msh\"\n",
- " \n",
- " lc=0.5\n",
+ "\n",
+ " lc = 0.5\n",
"\n",
" # Before using any functions in the Python API, Gmsh must be initialized:\n",
" gmsh.initialize()\n",
@@ -464,13 +565,12 @@
" dim2 = 2\n",
"\n",
" # Outer points (ccw)\n",
- " gmsh.model.geo.addPoint(0, -b, 0, lc, 1)\n",
- " gmsh.model.geo.addPoint(a, -b, 0, lc, 2)\n",
- " gmsh.model.geo.addPoint(a, -b/2, 0, lc, 3)\n",
- " gmsh.model.geo.addPoint(a, 0, 0, lc, 4)\n",
- " gmsh.model.geo.addPoint(0, 0, 0, lc, 5)\n",
- " gmsh.model.geo.addPoint(0, -b/2, 0, lc, 6)\n",
- "\n",
+ " gmsh.model.geo.addPoint(0, -b, 0, lc, 1)\n",
+ " gmsh.model.geo.addPoint(a, -b, 0, lc, 2)\n",
+ " gmsh.model.geo.addPoint(a, -b / 2, 0, lc, 3)\n",
+ " gmsh.model.geo.addPoint(a, 0, 0, lc, 4)\n",
+ " gmsh.model.geo.addPoint(0, 0, 0, lc, 5)\n",
+ " gmsh.model.geo.addPoint(0, -b / 2, 0, lc, 6)\n",
"\n",
" # Outer lines (ccw)\n",
" gmsh.model.geo.addLine(1, 2, 1)\n",
@@ -481,11 +581,10 @@
" gmsh.model.geo.addLine(6, 1, 6)\n",
" gmsh.model.geo.addLine(6, 3, 7)\n",
"\n",
- "\n",
- " # The third elementary entity is the surface. In order to define a surface \n",
+ " # The third elementary entity is the surface. In order to define a surface\n",
" # from the curves defined above, a curve loop has first to be defined (ccw).\n",
- " gmsh.model.geo.addCurveLoop([ 1, 2, -7, 6], 1)\n",
- " gmsh.model.geo.addCurveLoop([ 7, 3, 4, 5], 2)\n",
+ " gmsh.model.geo.addCurveLoop([1, 2, -7, 6], 1)\n",
+ " gmsh.model.geo.addCurveLoop([7, 3, 4, 5], 2)\n",
"\n",
" # Add plane surfaces defined by one or more curve loops.\n",
" gmsh.model.geo.addPlaneSurface([1], 1)\n",
@@ -494,13 +593,13 @@
" gmsh.model.geo.synchronize()\n",
"\n",
" # Prepare structured grid\n",
- " gmsh.model.geo.mesh.setTransfiniteCurve( 1, Nx)\n",
- " gmsh.model.geo.mesh.setTransfiniteCurve( 2, int(Ny*0.3))\n",
- " gmsh.model.geo.mesh.setTransfiniteCurve( 3, Ny, \"Progression\", -P)\n",
- " gmsh.model.geo.mesh.setTransfiniteCurve( 4, Nx)\n",
- " gmsh.model.geo.mesh.setTransfiniteCurve( 5, Ny, \"Progression\", P)\n",
- " gmsh.model.geo.mesh.setTransfiniteCurve( 6, int(Ny*0.3))\n",
- " gmsh.model.geo.mesh.setTransfiniteCurve( 7, Nx)\n",
+ " gmsh.model.geo.mesh.setTransfiniteCurve(1, Nx)\n",
+ " gmsh.model.geo.mesh.setTransfiniteCurve(2, int(Ny * 0.3))\n",
+ " gmsh.model.geo.mesh.setTransfiniteCurve(3, Ny, \"Progression\", -P)\n",
+ " gmsh.model.geo.mesh.setTransfiniteCurve(4, Nx)\n",
+ " gmsh.model.geo.mesh.setTransfiniteCurve(5, Ny, \"Progression\", P)\n",
+ " gmsh.model.geo.mesh.setTransfiniteCurve(6, int(Ny * 0.3))\n",
+ " gmsh.model.geo.mesh.setTransfiniteCurve(7, Nx)\n",
"\n",
" gmsh.model.geo.mesh.setTransfiniteSurface(1, \"Alternate\")\n",
" gmsh.model.geo.mesh.setTransfiniteSurface(2, \"Alternate\")\n",
@@ -530,10 +629,10 @@
"\n",
" gmsh.model.mesh.generate(dim2)\n",
" # gmsh.option.setNumber('Mesh.SecondOrderIncomplete', 1) # serendipity elements\n",
- " gmsh.model.mesh.setOrder(2) # higher order elements (quadratic)\n",
+ " gmsh.model.mesh.setOrder(2) # higher order elements (quadratic)\n",
" gmsh.write(output_file)\n",
"\n",
- " gmsh.finalize()\n"
+ " gmsh.finalize()"
]
},
{
@@ -541,6 +640,7 @@
"execution_count": 8,
"id": "137f610d-adbc-4e3f-9ea0-8ee958bbb817",
"metadata": {
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [
@@ -579,7 +679,7 @@
}
],
"source": [
- "generate_mesh_axb(200,100,25,45,1.07)"
+ "generate_mesh_axb(200, 100, 25, 45, 1.07)"
]
},
{
@@ -587,6 +687,7 @@
"execution_count": 9,
"id": "f0aaf7e1-9c00-4c16-a0ef-605cbf877377",
"metadata": {
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [],
@@ -599,6 +700,7 @@
"execution_count": 10,
"id": "a362da7c-dcc4-47f9-9e43-c5871f1006ed",
"metadata": {
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [
@@ -662,6 +764,7 @@
"execution_count": 11,
"id": "9c745eea-75c9-49db-8b40-833aca73b436",
"metadata": {
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [
@@ -681,7 +784,7 @@
"pv.set_jupyter_backend(\"static\")\n",
"\n",
"mesh = pv.read(f\"{out_dir}/square_200x100_domain.vtu\")\n",
- "plotter = pv.Plotter(window_size = [1000, 800])\n",
+ "plotter = pv.Plotter(window_size=[1000, 800])\n",
"plotter.add_mesh(mesh, show_edges=True, show_scalar_bar=False, color=None, scalars=None)\n",
"\n",
"plotter.show_bounds(ticks=\"outside\", xlabel=\"x / m\", ylabel=\"y / m\")\n",
@@ -702,7 +805,7 @@
"\n",
"In this example, the boundary conditions are defined as follows:\n",
"\n",
- "**Top:** \n",
+ "**Top:**\n",
"$$\n",
"\\begin{align}\n",
"p(y=0)&=\\tilde{p}\\cdot \\sin(\\omega \\cdot t) \\cdot \\cos(\\frac{2 \\pi}{L} \\cdot x) \\\\\n",
@@ -768,74 +871,87 @@
"outputs": [],
"source": [
"## Helper Functions\n",
- "def read_timestep_mesh(a,time):\n",
+ "def read_timestep_mesh(a, time):\n",
" reader = pv.PVDReader(f\"{out_dir}/square_{a}x100.pvd\")\n",
- " reader.set_active_time_point(int(time*4)) # time [s], delta t = 0.25 s\n",
+ " reader.set_active_time_point(int(time * 4)) # time [s], delta t = 0.25 s\n",
" mesh = reader.read()[0]\n",
" return mesh\n",
"\n",
+ "\n",
"def slice_along_line(mesh, start_point, end_point):\n",
- " line = pv.Line(start_point, end_point, resolution = 2)\n",
+ " line = pv.Line(start_point, end_point, resolution=2)\n",
" return mesh.slice_along_line(line)\n",
"\n",
+ "\n",
"def get_pressure_sorted(mesh):\n",
" pressure = mesh.point_data[\"pressure_interpolated\"]\n",
- " depth = mesh.points[:,1]\n",
+ " depth = mesh.points[:, 1]\n",
" indices_sorted = np.argsort(depth)\n",
" pressure_sorted = pressure[indices_sorted]\n",
" return pressure_sorted\n",
"\n",
+ "\n",
"def get_stresses_sorted(mesh):\n",
" sigma = mesh.point_data[\"sigma\"]\n",
- " depth = mesh.points[:,1]\n",
+ " depth = mesh.points[:, 1]\n",
" indices_sorted = np.argsort(depth)\n",
- " sigma_xx = - sigma[indices_sorted, 0] # switching sign convention\n",
- " sigma_yy = - sigma[indices_sorted, 1]\n",
- " #sigma_zz = - sigma[indices_sorted, 2]\n",
- " sigma_xy = + sigma[indices_sorted, 3]\n",
- " return sigma_xx, sigma_yy, sigma_xy #,sigma_zz\n",
+ " sigma_xx = -sigma[indices_sorted, 0] # switching sign convention\n",
+ " sigma_yy = -sigma[indices_sorted, 1]\n",
+ " # sigma_zz = - sigma[indices_sorted, 2]\n",
+ " sigma_xy = +sigma[indices_sorted, 3]\n",
+ " return sigma_xx, sigma_yy, sigma_xy # ,sigma_zz\n",
+ "\n",
"\n",
"def get_depth_sorted(mesh):\n",
- " depth = mesh.points[:,1]\n",
+ " depth = mesh.points[:, 1]\n",
" indices_sorted = np.argsort(depth)\n",
" return depth[indices_sorted]\n",
"\n",
+ "\n",
"def compute_abs_and_rel_pressure_error(pressures, depth, t, x):\n",
" num_points = pressures.shape[0]\n",
" f_abs = np.zeros(num_points)\n",
" f_rel = np.zeros(num_points)\n",
- " \n",
- " for pt_idx in range(num_points): \n",
- " y = -depth[pt_idx]\n",
- " pressure_ana = compute_pressure_and_stresses(t,x,y)[0] # returns pressure normalised to the pressure amplitude\n",
- " pressure_num = pressures[pt_idx]/0.1e5 # absolute pressure divided by pressure amplitude\n",
+ "\n",
+ " for pt_idx in range(num_points):\n",
+ " y = -depth[pt_idx]\n",
+ " pressure_ana = compute_pressure_and_stresses(t, x, y)[\n",
+ " 0\n",
+ " ] # returns pressure normalised to the pressure amplitude\n",
+ " pressure_num = (\n",
+ " pressures[pt_idx] / 0.1e5\n",
+ " ) # absolute pressure divided by pressure amplitude\n",
" f_abs[pt_idx] = pressure_num - pressure_ana\n",
- " \n",
+ "\n",
" if pressure_ana == 0:\n",
" f_rel[pt_idx] = f_abs[pt_idx] / 1e-2\n",
" else:\n",
" f_rel[pt_idx] = f_abs[pt_idx] / pressure_ana\n",
- " \n",
+ "\n",
" return f_abs, f_rel\n",
"\n",
+ "\n",
"def compute_abs_and_rel_stress_error(sigmas, depth, t, x):\n",
" num_points = depth.shape[0]\n",
" f_abs = np.zeros((3, num_points))\n",
" f_rel = np.zeros((3, num_points))\n",
- " \n",
- " for stress_idx in (0,1,2):\n",
- " \n",
+ "\n",
+ " for stress_idx in (0, 1, 2):\n",
" for pt_idx in range(num_points):\n",
" y = -depth[pt_idx]\n",
- " sigma_ana = compute_pressure_and_stresses(t,x,y)[stress_idx+1] # returns stresses normalised to the pressure amplitude\n",
- " sigma_num = sigma[stress_idx][pt_idx]/0.1e5 # absolute stresses divided by pressure amplitude\n",
+ " sigma_ana = compute_pressure_and_stresses(t, x, y)[\n",
+ " stress_idx + 1\n",
+ " ] # returns stresses normalised to the pressure amplitude\n",
+ " sigma_num = (\n",
+ " sigma[stress_idx][pt_idx] / 0.1e5\n",
+ " ) # absolute stresses divided by pressure amplitude\n",
" f_abs[stress_idx][pt_idx] = sigma_num - sigma_ana\n",
"\n",
" if sigma_ana == 0:\n",
- " f_rel[stress_idx][pt_idx] = f_abs[stress_idx][pt_idx]/ 1e-2\n",
+ " f_rel[stress_idx][pt_idx] = f_abs[stress_idx][pt_idx] / 1e-2\n",
" else:\n",
" f_rel[stress_idx][pt_idx] = f_abs[stress_idx][pt_idx] / sigma_ana\n",
- " \n",
+ "\n",
" return f_abs, f_rel"
]
},
@@ -844,6 +960,7 @@
"execution_count": 14,
"id": "848ce73f-688f-4b19-94ef-a4718ad61570",
"metadata": {
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [
@@ -857,9 +974,10 @@
}
],
"source": [
- "model = ogs.OGS(INPUT_FILE=\"seabed_response_200x100.prj\", PROJECT_FILE=\"seabed_response_200x100.prj\")\n",
- "model.run_model(logfile=f\"{out_dir}/out.txt\",\n",
- " args=f\"-o {out_dir} -m {out_dir}\")"
+ "model = ogs.OGS(\n",
+ " INPUT_FILE=\"seabed_response_200x100.prj\", PROJECT_FILE=\"seabed_response_200x100.prj\"\n",
+ ")\n",
+ "model.run_model(logfile=f\"{out_dir}/out.txt\", args=f\"-o {out_dir} -m {out_dir}\")"
]
},
{
@@ -875,6 +993,7 @@
"execution_count": 15,
"id": "ed7c61b9-18e9-48d5-87d7-2fb3f932adf9",
"metadata": {
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [
@@ -892,14 +1011,21 @@
"source": [
"time = 2.5 # [s]\n",
"reader = pv.get_reader(f\"{out_dir}/square_200x100.pvd\")\n",
- "reader.set_active_time_point(int(time*4))\n",
+ "reader.set_active_time_point(int(time * 4))\n",
"mesh = reader.read()[0]\n",
"\n",
"plotter = pv.Plotter()\n",
"\n",
- "sargs = dict(title=\"p / Pa\" , height=0.25, position_x=0.2, position_y=0.02)\n",
- "plotter.add_mesh(mesh, scalars = \"pressure_interpolated\", show_edges=False, show_scalar_bar=True, label=\"p\", scalar_bar_args=sargs)\n",
- "plotter.show_bounds(ticks=\"outside\", xlabel = \"x / m\", ylabel = \"y / m\")\n",
+ "sargs = dict(title=\"p / Pa\", height=0.25, position_x=0.2, position_y=0.02)\n",
+ "plotter.add_mesh(\n",
+ " mesh,\n",
+ " scalars=\"pressure_interpolated\",\n",
+ " show_edges=False,\n",
+ " show_scalar_bar=True,\n",
+ " label=\"p\",\n",
+ " scalar_bar_args=sargs,\n",
+ ")\n",
+ "plotter.show_bounds(ticks=\"outside\", xlabel=\"x / m\", ylabel=\"y / m\")\n",
"plotter.add_axes()\n",
"plotter.view_xy()\n",
"plotter.show()"
@@ -910,9 +1036,9 @@
"id": "d4631296-5f3a-42fa-96d1-5ce0556bf206",
"metadata": {},
"source": [
- "For a more detailed comparison between the analytical and the numerical solution, both solutions are evaluated along the vertical line directly underneath an anti-node of the standing wave. As before, the pore pressure and the amplitude of the effective stresses are illustrated as a function of depth. The results of the numerical solution are marked as dots in the same color as the analytical solution. Additionally, the absolute errors $\\Delta p = p_{numerical}-p_{analtical}$ and $\\Delta \\sigma_{i}' = \\sigma_{i, numerical}'-\\sigma_{i, analytical}'$ are illustrated on the right. \n",
+ "For a more detailed comparison between the analytical and the numerical solution, both solutions are evaluated along the vertical line directly underneath an anti-node of the standing wave. As before, the pore pressure and the amplitude of the effective stresses are illustrated as a function of depth. The results of the numerical solution are marked as dots in the same color as the analytical solution. Additionally, the absolute errors $\\Delta p = p_{numerical}-p_{analtical}$ and $\\Delta \\sigma_{i}' = \\sigma_{i, numerical}'-\\sigma_{i, analytical}'$ are illustrated on the right.\n",
"\n",
- "The plot shows that the absolute errors are very small at about $2 \\%$ of the wave's amplitude. They can mostly be ascribed to the space- and time-discretization. Close to the top boundary of the domain, larger errors occur. These errors could originate in the definition of both a pressure and displacement (Neumann-) boundary condition along the top edge. "
+ "The plot shows that the absolute errors are very small at about $2 \\%$ of the wave's amplitude. They can mostly be ascribed to the space- and time-discretization. Close to the top boundary of the domain, larger errors occur. These errors could originate in the definition of both a pressure and displacement (Neumann-) boundary condition along the top edge."
]
},
{
@@ -920,6 +1046,7 @@
"execution_count": 16,
"id": "2e06842d-1f73-4182-9ffa-02498a9f9341",
"metadata": {
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [
@@ -936,28 +1063,55 @@
],
"source": [
"x = 0\n",
- "y = np.linspace(0,100,1000)\n",
- "y_rel = y/100\n",
- "colors = {0:\"orangered\", 2:\"gold\", 4:\"blueviolet\", 6:\"forestgreen\", 8:\"darkorange\", 10:\"royalblue\"}\n",
- "\n",
- "fig, ax = plt.subplots(ncols=2, nrows=2, figsize=(15,15))\n",
+ "y = np.linspace(0, 100, 1000)\n",
+ "y_rel = y / 100\n",
+ "colors = {\n",
+ " 0: \"orangered\",\n",
+ " 2: \"gold\",\n",
+ " 4: \"blueviolet\",\n",
+ " 6: \"forestgreen\",\n",
+ " 8: \"darkorange\",\n",
+ " 10: \"royalblue\",\n",
+ "}\n",
+ "\n",
+ "fig, ax = plt.subplots(ncols=2, nrows=2, figsize=(15, 15))\n",
"\n",
"## Plotting analytical solution\n",
- "for t in [2,4,6,8,10]:\n",
- " ax[0][0].plot(compute_pressure_and_stresses(t,x,y)[0], -y_rel, color=colors[t], label= \"analytical, t = %.1f s\" %t)\n",
- "\n",
- "ax[1][0].plot(compute_pressure_and_stresses(2.5,x,y)[1], -y_rel, color = colors[6], label = \"analytical, $\\\\sigma'_{xx}/\\\\alpha\\\\tilde{p}$\")\n",
- "ax[1][0].plot(compute_pressure_and_stresses(2.5,x,y)[2], -y_rel, color = colors[2], label = \"analytical, $\\\\sigma'_{yy}/\\\\alpha\\\\tilde{p}$\")\n",
- "ax[1][0].plot(compute_pressure_and_stresses(2.5,x,y)[3], -y_rel, color = colors[4], label = \"analytical, $\\\\sigma'_{xy}/\\\\alpha\\\\tilde{p}$\")\n",
+ "for t in [2, 4, 6, 8, 10]:\n",
+ " ax[0][0].plot(\n",
+ " compute_pressure_and_stresses(t, x, y)[0],\n",
+ " -y_rel,\n",
+ " color=colors[t],\n",
+ " label=\"analytical, t = %.1f s\" % t,\n",
+ " )\n",
+ "\n",
+ "ax[1][0].plot(\n",
+ " compute_pressure_and_stresses(2.5, x, y)[1],\n",
+ " -y_rel,\n",
+ " color=colors[6],\n",
+ " label=\"analytical, $\\\\sigma'_{xx}/\\\\alpha\\\\tilde{p}$\",\n",
+ ")\n",
+ "ax[1][0].plot(\n",
+ " compute_pressure_and_stresses(2.5, x, y)[2],\n",
+ " -y_rel,\n",
+ " color=colors[2],\n",
+ " label=\"analytical, $\\\\sigma'_{yy}/\\\\alpha\\\\tilde{p}$\",\n",
+ ")\n",
+ "ax[1][0].plot(\n",
+ " compute_pressure_and_stresses(2.5, x, y)[3],\n",
+ " -y_rel,\n",
+ " color=colors[4],\n",
+ " label=\"analytical, $\\\\sigma'_{xy}/\\\\alpha\\\\tilde{p}$\",\n",
+ ")\n",
"\n",
"## Plotting numerical solution\n",
- "p1 = (x+1e-6, 0, 0)\n",
- "p2 = (x+1e-6, -100, 0)\n",
+ "p1 = (x + 1e-6, 0, 0)\n",
+ "p2 = (x + 1e-6, -100, 0)\n",
"\n",
"for t_num in (2, 2.5, 4, 6, 8, 10):\n",
" mesh = read_timestep_mesh(200, t_num)\n",
- " \n",
- " line_mesh = slice_along_line(mesh, p1, p2) \n",
+ "\n",
+ " line_mesh = slice_along_line(mesh, p1, p2)\n",
" pressure = get_pressure_sorted(line_mesh)\n",
" sigma = get_stresses_sorted(line_mesh)\n",
" depth = get_depth_sorted(line_mesh)\n",
@@ -965,31 +1119,76 @@
" f_abs_sigma = compute_abs_and_rel_stress_error(sigma, depth, t_num, x)[0]\n",
"\n",
" if t_num != 2.5:\n",
- " ax[0][0].plot(pressure/0.1e5, depth/100, \"o\", markevery=10, color=colors[t_num], label= \"numerical, t = %.1f s\" %t_num) \n",
+ " ax[0][0].plot(\n",
+ " pressure / 0.1e5,\n",
+ " depth / 100,\n",
+ " \"o\",\n",
+ " markevery=10,\n",
+ " color=colors[t_num],\n",
+ " label=\"numerical, t = %.1f s\" % t_num,\n",
+ " )\n",
" ax[0][0].set_xlabel(\"$p$ / $\\\\tilde{p}$\")\n",
"\n",
- " ax[0][1].plot(f_abs_pressure, depth/100, color=colors[t_num], label = \"t = %.1f s\" %t_num)\n",
+ " ax[0][1].plot(\n",
+ " f_abs_pressure, depth / 100, color=colors[t_num], label=\"t = %.1f s\" % t_num\n",
+ " )\n",
" ax[0][1].set_xlabel(\"$\\\\Delta p /\\\\tilde{p}$\")\n",
- " \n",
+ "\n",
" if t_num == 2.5:\n",
- " ax[1][0].plot(sigma[0]/0.1e5, depth/100, \"o\", markevery=10, color = colors[6], label = \"numerical, $\\\\sigma'_{xx}/\\\\alpha\\\\tilde{p}$\")\n",
- " ax[1][0].plot(sigma[1]/0.1e5, depth/100, \"o\", markevery=10, color = colors[2], label = \"numerical, $\\\\sigma'_{yy}/\\\\alpha\\\\tilde{p}$\")\n",
- " ax[1][0].plot(sigma[2]/0.1e5, depth/100, \"o\", markevery=10, color = colors[4], label = \"numerical, $\\\\sigma'_{xy}/\\\\alpha\\\\tilde{p}$\")\n",
+ " ax[1][0].plot(\n",
+ " sigma[0] / 0.1e5,\n",
+ " depth / 100,\n",
+ " \"o\",\n",
+ " markevery=10,\n",
+ " color=colors[6],\n",
+ " label=\"numerical, $\\\\sigma'_{xx}/\\\\alpha\\\\tilde{p}$\",\n",
+ " )\n",
+ " ax[1][0].plot(\n",
+ " sigma[1] / 0.1e5,\n",
+ " depth / 100,\n",
+ " \"o\",\n",
+ " markevery=10,\n",
+ " color=colors[2],\n",
+ " label=\"numerical, $\\\\sigma'_{yy}/\\\\alpha\\\\tilde{p}$\",\n",
+ " )\n",
+ " ax[1][0].plot(\n",
+ " sigma[2] / 0.1e5,\n",
+ " depth / 100,\n",
+ " \"o\",\n",
+ " markevery=10,\n",
+ " color=colors[4],\n",
+ " label=\"numerical, $\\\\sigma'_{xy}/\\\\alpha\\\\tilde{p}$\",\n",
+ " )\n",
" ax[1][0].set_xlabel(\"$\\\\sigma$'/$\\\\alpha\\\\tilde{p}$\")\n",
- " \n",
- " ax[1][1].plot(f_abs_sigma[0], depth/100, color = colors[6], label = \"$\\\\Delta\\\\sigma'_{xx}/\\\\alpha\\\\tilde{p}$\")\n",
- " ax[1][1].plot(f_abs_sigma[1], depth/100, color = colors[2], label = \"$\\\\Delta\\\\sigma'_{yy}/\\\\alpha\\\\tilde{p}$\")\n",
- " ax[1][1].plot(f_abs_sigma[2], depth/100, color = colors[4], label = \"$\\\\Delta\\\\sigma'_{xy}/\\\\alpha\\\\tilde{p}$\")\n",
+ "\n",
+ " ax[1][1].plot(\n",
+ " f_abs_sigma[0],\n",
+ " depth / 100,\n",
+ " color=colors[6],\n",
+ " label=\"$\\\\Delta\\\\sigma'_{xx}/\\\\alpha\\\\tilde{p}$\",\n",
+ " )\n",
+ " ax[1][1].plot(\n",
+ " f_abs_sigma[1],\n",
+ " depth / 100,\n",
+ " color=colors[2],\n",
+ " label=\"$\\\\Delta\\\\sigma'_{yy}/\\\\alpha\\\\tilde{p}$\",\n",
+ " )\n",
+ " ax[1][1].plot(\n",
+ " f_abs_sigma[2],\n",
+ " depth / 100,\n",
+ " color=colors[4],\n",
+ " label=\"$\\\\Delta\\\\sigma'_{xy}/\\\\alpha\\\\tilde{p}$\",\n",
+ " )\n",
" ax[1][1].set_xlabel(\"$\\\\Delta\\\\sigma$'/$\\\\alpha\\\\tilde{p}$\")\n",
- " \n",
- " #ax[1][0].plot(sigma[3]/0.1e5, depth/100, \"o\", markevery=10, color = colors[4], label = \"numerical, $\\\\sigma'_{zz}/\\\\alpha\\\\tilde{p}$\")\n",
+ "\n",
+ " # ax[1][0].plot(sigma[3]/0.1e5, depth/100, \"o\", markevery=10, color = colors[4], label = \"numerical, $\\\\sigma'_{zz}/\\\\alpha\\\\tilde{p}$\")\n",
"\n",
"## layout settings\n",
- "ax[0][0].set_title('Comparison numerical and analytical solution')\n",
- "ax[0][1].set_title('Absolute error')\n",
+ "ax[0][0].set_title(\"Comparison numerical and analytical solution\")\n",
+ "ax[0][1].set_title(\"Absolute error\")\n",
"\n",
- "for idx_1 in (0,1):\n",
- " for idx_2 in (0,1):\n",
+ "for idx_1 in (0, 1):\n",
+ " for idx_2 in (0, 1):\n",
" ax[idx_1][idx_2].grid(True)\n",
" ax[idx_1][idx_2].set_ylabel(\"$y$ / $L$\")\n",
" ax[idx_1][0].set_xlim(-1.1, 1.1)\n",
diff --git a/Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole_convergence_analysis.ipynb b/Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole_convergence_analysis.ipynb
index d5ff704aebc..1a801d39c0e 100644
--- a/Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole_convergence_analysis.ipynb
+++ b/Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole_convergence_analysis.ipynb
@@ -1,15 +1,16 @@
{
"cells": [
{
- "cell_type": "raw",
+ "cell_type": "markdown",
"id": "bb0907b4-4e26-4c4e-ab1f-22b5330cb1d2",
"metadata": {},
"source": [
+ "+++\n",
"title = \"Linear elasticity: disc with hole convergence study\"\n",
"date = \"2022-09-15\"\n",
"author = \"Linda Günther, Sophia Einspänner, Robert Habel, Christoph Lehmann and Thomas Nagel\"\n",
"web_subsection = \"small-deformations\"\n",
- ""
+ "+++"
]
},
{
@@ -76,6 +77,7 @@
"jupyter": {
"source_hidden": true
},
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [],
@@ -107,7 +109,8 @@
"metadata": {
"jupyter": {
"source_hidden": true
- }
+ },
+ "lines_to_next_cell": 2
},
"outputs": [],
"source": [
@@ -307,6 +310,7 @@
"jupyter": {
"source_hidden": true
},
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [],
@@ -322,6 +326,7 @@
"jupyter": {
"source_hidden": true
},
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [
@@ -562,6 +567,7 @@
"jupyter": {
"source_hidden": true
},
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [],
@@ -597,6 +603,7 @@
"jupyter": {
"source_hidden": true
},
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [],
@@ -615,6 +622,7 @@
"jupyter": {
"source_hidden": true
},
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [
@@ -666,6 +674,7 @@
"jupyter": {
"source_hidden": true
},
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [],
@@ -920,6 +929,7 @@
"jupyter": {
"source_hidden": true
},
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [],
@@ -947,6 +957,7 @@
"jupyter": {
"source_hidden": true
},
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [],
@@ -975,6 +986,7 @@
"jupyter": {
"source_hidden": true
},
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [],
@@ -1042,6 +1054,7 @@
"jupyter": {
"source_hidden": true
},
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [
@@ -1642,6 +1655,7 @@
"jupyter": {
"source_hidden": true
},
+ "lines_to_next_cell": 2,
"tags": []
},
"outputs": [
diff --git a/Tests/Data/Mechanics/Linear/SimpleMechanics.ipynb b/Tests/Data/Mechanics/Linear/SimpleMechanics.ipynb
index bce8e212961..c69a1a6da29 100644
--- a/Tests/Data/Mechanics/Linear/SimpleMechanics.ipynb
+++ b/Tests/Data/Mechanics/Linear/SimpleMechanics.ipynb
@@ -5,11 +5,12 @@
"id": "96f29a77",
"metadata": {},
"source": [
+ "+++\n",
"title = \"SimpleMechanics\"\n",
"date = \"2021-09-10\"\n",
"author = \"Lars Bilke, Jörg Buchwald\"\n",
"web_subsection = \"small-deformations\"\n",
- ""
+ "+++\n"
]
},
{
@@ -24,7 +25,9 @@
"cell_type": "code",
"execution_count": 19,
"id": "420713a5-74d6-47ad-815c-4ca9e7e914bf",
- "metadata": {},
+ "metadata": {
+ "lines_to_next_cell": 2
+ },
"outputs": [],
"source": [
"import os\n",
@@ -38,7 +41,9 @@
"cell_type": "code",
"execution_count": 27,
"id": "8da3a8e8-be97-4092-88a9-1fb7792fa644",
- "metadata": {},
+ "metadata": {
+ "lines_to_next_cell": 2
+ },
"outputs": [
{
"name": "stdout",
@@ -179,7 +184,9 @@
"cell_type": "code",
"execution_count": 26,
"id": "1d730e79",
- "metadata": {},
+ "metadata": {
+ "lines_to_next_cell": 2
+ },
"outputs": [
{
"data": {
diff --git a/Tests/Data/Mechanics/PLLC/PLLC.ipynb b/Tests/Data/Mechanics/PLLC/PLLC.ipynb
index ee751646246..5c5de5df16c 100644
--- a/Tests/Data/Mechanics/PLLC/PLLC.ipynb
+++ b/Tests/Data/Mechanics/PLLC/PLLC.ipynb
@@ -5,11 +5,12 @@
"id": "bb0907b4-4e26-4c4e-ab1f-22b5330cb1d2",
"metadata": {},
"source": [
+ "+++\n",
"title = \"Power Law Linear Creep\"\n",
"date = \"2023-01-02\"\n",
"author = \"Florian Zill\"\n",
"web_subsection = \"small-deformations\"\n",
- ""
+ "+++\n"
]
},
{
@@ -27,7 +28,9 @@
"cell_type": "code",
"execution_count": null,
"id": "7962f42f-fd53-4fc1-b966-a8ba924aca6c",
- "metadata": {},
+ "metadata": {
+ "lines_to_next_cell": 2
+ },
"outputs": [],
"source": [
"import contextlib\n",
@@ -39,16 +42,18 @@
"from ogs6py import ogs\n",
"\n",
"prj_name = \"uniax_compression\"\n",
- "data_dir = os.environ.get('OGS_DATA_DIR', str(os.getcwd()).split(\"/Data/\")[0] + \"/Data/\")\n",
+ "data_dir = os.environ.get(\n",
+ " \"OGS_DATA_DIR\", str(os.getcwd()).split(\"/Data/\")[0] + \"/Data/\"\n",
+ ")\n",
"input_file = f\"{data_dir}/Mechanics/PLLC/{prj_name}.prj\"\n",
- "out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', f'{data_dir}/Mechanics/PLLC/_out')\n",
+ "out_dir = os.environ.get(\"OGS_TESTRUNNER_OUT_DIR\", f\"{data_dir}/Mechanics/PLLC/_out\")\n",
"\n",
"if not os.path.exists(out_dir):\n",
" os.makedirs(out_dir)\n",
"os.chdir(out_dir)\n",
"\n",
"prj_file = f\"{out_dir}/{prj_name}_out.prj\"\n",
- "ogs_model = ogs.OGS(INPUT_FILE=input_file, PROJECT_FILE=prj_file)\n"
+ "ogs_model = ogs.OGS(INPUT_FILE=input_file, PROJECT_FILE=prj_file)"
]
},
{
@@ -66,18 +71,123 @@
"cell_type": "code",
"execution_count": null,
"id": "5a20a14e",
- "metadata": {},
+ "metadata": {
+ "lines_to_next_cell": 2
+ },
"outputs": [],
"source": [
"# Unfortunately the source for the WIPP data has gone missing - will be added if it's found again\n",
"ExData = {\n",
- " \"WIPP CS 25\": (25, \"^\", [[9.87970002, 2.013560846E-05], [11.84642707, 3.178356756E-05], [7.87388785, 1.66059726E-06]]),\n",
- " \"WIPP CS 60\": (60, \"^\", [[3.98589289, 5.7824853E-06], [5.94266985, 2.075776623E-05], [7.87388785, 1.953209818E-05], [9.96978837, 5.841438703E-05], [11.84642707, 0.00011762092257], [13.94911482, 0.00026749321794], [17.9857158, 0.00111804208073], [1.9814251, 8.7645834E-07], [3.91418422, 4.01350889E-06], [5.88897108, 3.34371363E-06], [7.87388785, 1.129440706E-05], [9.87970002, 2.99068674E-05], [11.84642707, 7.681792203E-05], [13.82306874, 0.00011067584933], [15.83934389, 0.00052247037957]]),\n",
- " \"DeVries 1988 25\": (25, \"s\", [[4.99, 2.10816E-06], [4.99, 2.4192E-06], [5, 1.8144E-06], [9.99, 2.2032E-05], [14.96, 9.2448E-05], [14.98, 0.000216]]),\n",
- " \"DeVries 1988 100\": (100, \"s\", [[4.95, 9.6768E-05], [6.77, 0.000292896], [7.46, 0.000324], [8.55, 0.000664416], [8.92, 0.00091584], [8.98, 0.0009936], [9.91, 0.00124416], [10.1, 0.00139968], [10.22, 0.00093312], [10.27, 0.00132192], [12.1, 0.00216], [12.3, 0.00409536], [12.35, 0.00320544], [12.37, 0.00292032], [12.39, 0.00253152], [12.4, 0.0026784], [12.46, 0.0025056], [12.49, 0.00347328], [13.57, 0.00273024], [13.78, 0.00242784], [14.7, 0.00482112], [16.87, 0.0095904], [17.2, 0.0123552], [19.96, 0.030672]]),\n",
- " \"DeVries 1988 200\": (200, \"s\", [[3.47, 0.00117504], [4.71, 0.0032832], [6.67, 0.0104544], [6.78, 0.0132192], [9.86, 0.214272]]),\n",
- " \"Berest 2015 14.3\": (14.3, \"P\", [[0.09909639, 8.944207E-08], [0.19575886, 1.4118213E-07], [0.29452325, 1.4118213E-07], [0.49411031, 9.799173E-08]]),\n",
- " \"Berest 2017 7.8\": (7.8, \"P\", [[0.19575886,2.2285256E-07], [0.19575886,9.505469E-08], [0.19754389,2.5947583E-07], [0.19754389,2.647936E-08], [0.39379426,4.9162047E-07], [0.39738509,6.801413E-08], [0.59247161,4.0957628E-07], [0.59247161,5.7241269E-07], [0.59787408,1.0735864E-07], [1.0591736,1.11804208E-06]])}"
+ " \"WIPP CS 25\": (\n",
+ " 25,\n",
+ " \"^\",\n",
+ " [\n",
+ " [9.87970002, 2.013560846e-05],\n",
+ " [11.84642707, 3.178356756e-05],\n",
+ " [7.87388785, 1.66059726e-06],\n",
+ " ],\n",
+ " ),\n",
+ " \"WIPP CS 60\": (\n",
+ " 60,\n",
+ " \"^\",\n",
+ " [\n",
+ " [3.98589289, 5.7824853e-06],\n",
+ " [5.94266985, 2.075776623e-05],\n",
+ " [7.87388785, 1.953209818e-05],\n",
+ " [9.96978837, 5.841438703e-05],\n",
+ " [11.84642707, 0.00011762092257],\n",
+ " [13.94911482, 0.00026749321794],\n",
+ " [17.9857158, 0.00111804208073],\n",
+ " [1.9814251, 8.7645834e-07],\n",
+ " [3.91418422, 4.01350889e-06],\n",
+ " [5.88897108, 3.34371363e-06],\n",
+ " [7.87388785, 1.129440706e-05],\n",
+ " [9.87970002, 2.99068674e-05],\n",
+ " [11.84642707, 7.681792203e-05],\n",
+ " [13.82306874, 0.00011067584933],\n",
+ " [15.83934389, 0.00052247037957],\n",
+ " ],\n",
+ " ),\n",
+ " \"DeVries 1988 25\": (\n",
+ " 25,\n",
+ " \"s\",\n",
+ " [\n",
+ " [4.99, 2.10816e-06],\n",
+ " [4.99, 2.4192e-06],\n",
+ " [5, 1.8144e-06],\n",
+ " [9.99, 2.2032e-05],\n",
+ " [14.96, 9.2448e-05],\n",
+ " [14.98, 0.000216],\n",
+ " ],\n",
+ " ),\n",
+ " \"DeVries 1988 100\": (\n",
+ " 100,\n",
+ " \"s\",\n",
+ " [\n",
+ " [4.95, 9.6768e-05],\n",
+ " [6.77, 0.000292896],\n",
+ " [7.46, 0.000324],\n",
+ " [8.55, 0.000664416],\n",
+ " [8.92, 0.00091584],\n",
+ " [8.98, 0.0009936],\n",
+ " [9.91, 0.00124416],\n",
+ " [10.1, 0.00139968],\n",
+ " [10.22, 0.00093312],\n",
+ " [10.27, 0.00132192],\n",
+ " [12.1, 0.00216],\n",
+ " [12.3, 0.00409536],\n",
+ " [12.35, 0.00320544],\n",
+ " [12.37, 0.00292032],\n",
+ " [12.39, 0.00253152],\n",
+ " [12.4, 0.0026784],\n",
+ " [12.46, 0.0025056],\n",
+ " [12.49, 0.00347328],\n",
+ " [13.57, 0.00273024],\n",
+ " [13.78, 0.00242784],\n",
+ " [14.7, 0.00482112],\n",
+ " [16.87, 0.0095904],\n",
+ " [17.2, 0.0123552],\n",
+ " [19.96, 0.030672],\n",
+ " ],\n",
+ " ),\n",
+ " \"DeVries 1988 200\": (\n",
+ " 200,\n",
+ " \"s\",\n",
+ " [\n",
+ " [3.47, 0.00117504],\n",
+ " [4.71, 0.0032832],\n",
+ " [6.67, 0.0104544],\n",
+ " [6.78, 0.0132192],\n",
+ " [9.86, 0.214272],\n",
+ " ],\n",
+ " ),\n",
+ " \"Berest 2015 14.3\": (\n",
+ " 14.3,\n",
+ " \"P\",\n",
+ " [\n",
+ " [0.09909639, 8.944207e-08],\n",
+ " [0.19575886, 1.4118213e-07],\n",
+ " [0.29452325, 1.4118213e-07],\n",
+ " [0.49411031, 9.799173e-08],\n",
+ " ],\n",
+ " ),\n",
+ " \"Berest 2017 7.8\": (\n",
+ " 7.8,\n",
+ " \"P\",\n",
+ " [\n",
+ " [0.19575886, 2.2285256e-07],\n",
+ " [0.19575886, 9.505469e-08],\n",
+ " [0.19754389, 2.5947583e-07],\n",
+ " [0.19754389, 2.647936e-08],\n",
+ " [0.39379426, 4.9162047e-07],\n",
+ " [0.39738509, 6.801413e-08],\n",
+ " [0.59247161, 4.0957628e-07],\n",
+ " [0.59247161, 5.7241269e-07],\n",
+ " [0.59787408, 1.0735864e-07],\n",
+ " [1.0591736, 1.11804208e-06],\n",
+ " ],\n",
+ " ),\n",
+ "}"
]
},
{
@@ -95,18 +205,27 @@
"cell_type": "code",
"execution_count": null,
"id": "8066a6d3",
- "metadata": {},
+ "metadata": {
+ "lines_to_next_cell": 2
+ },
"outputs": [],
"source": [
- "A1 = 0.18 # d^-1\n",
- "Q1 = 54e3 # kJ / mol\n",
- "A2 = 6.5e-5 # m^3 K d^−1\n",
- "Q2 = 24.5e3 # kJ / mol\n",
- "dGrain = 5e-2 # m\n",
- "sref = 1. # MPa\n",
- "BGRa = lambda sig, T: A1 * np.exp(-Q1/(8.3145*(273.15+T))) * np.power(sig/sref,5.)\n",
- "PLLC = lambda sig, T: A1 * np.exp(-Q1/(8.3145*(273.15+T))) * np.power(sig/sref,5.) + \\\n",
- " A2 * np.exp(-Q2/(8.3145*(273.15+T))) * sig/sref / np.power(dGrain, 3) / (273.15+T)"
+ "A1 = 0.18 # d^-1\n",
+ "Q1 = 54e3 # kJ / mol\n",
+ "A2 = 6.5e-5 # m^3 K d^−1\n",
+ "Q2 = 24.5e3 # kJ / mol\n",
+ "dGrain = 5e-2 # m\n",
+ "sref = 1.0 # MPa\n",
+ "BGRa = (\n",
+ " lambda sig, T: A1\n",
+ " * np.exp(-Q1 / (8.3145 * (273.15 + T)))\n",
+ " * np.power(sig / sref, 5.0)\n",
+ ")\n",
+ "PLLC = lambda sig, T: A1 * np.exp(-Q1 / (8.3145 * (273.15 + T))) * np.power(\n",
+ " sig / sref, 5.0\n",
+ ") + A2 * np.exp(-Q2 / (8.3145 * (273.15 + T))) * sig / sref / np.power(dGrain, 3) / (\n",
+ " 273.15 + T\n",
+ ")"
]
},
{
@@ -124,29 +243,36 @@
"cell_type": "code",
"execution_count": null,
"id": "7e2e294c-e803-4f02-b5ab-9bdfef94b00f",
- "metadata": {},
+ "metadata": {
+ "lines_to_next_cell": 2
+ },
"outputs": [],
"source": [
"lo_stresses = np.array([0.2e6, 0.6e6])\n",
"hi_stresses = np.array([2e6, 10e6])\n",
- "Exps = {7.8: ('blue', lo_stresses), 14.3: ('orange', lo_stresses),\n",
- " 25: ('lime', hi_stresses), 60: ('red', hi_stresses),\n",
- " 100: ('gray', hi_stresses), 200: ('mediumpurple', hi_stresses)}\n",
+ "Exps = {\n",
+ " 7.8: (\"blue\", lo_stresses),\n",
+ " 14.3: (\"orange\", lo_stresses),\n",
+ " 25: (\"lime\", hi_stresses),\n",
+ " 60: (\"red\", hi_stresses),\n",
+ " 100: (\"gray\", hi_stresses),\n",
+ " 200: (\"mediumpurple\", hi_stresses),\n",
+ "}\n",
"\n",
"fig, ax = plt.subplots(1, 1, figsize=(8, 6))\n",
- "ax.set_xlabel('$\\\\sigma_\\\\mathrm{ax}$ / MPa')\n",
- "ax.set_ylabel('$\\\\dot{\\\\epsilon}_{zz}$ / d$^{-1}$')\n",
+ "ax.set_xlabel(\"$\\\\sigma_\\\\mathrm{ax}$ / MPa\")\n",
+ "ax.set_ylabel(\"$\\\\dot{\\\\epsilon}_{zz}$ / d$^{-1}$\")\n",
"ax.set_xlim(0.15, 30)\n",
"ax.set_ylim(1e-15, 1e1)\n",
- "ax.grid(visible=True, which='both')\n",
- "points = {'pt0': (1., 1., 1.)}\n",
+ "ax.grid(visible=True, which=\"both\")\n",
+ "points = {\"pt0\": (1.0, 1.0, 1.0)}\n",
"\n",
"sigs = np.logspace(-1, 2, 100)\n",
- "for temp, (col, stresses) in Exps.items(): \n",
+ "for temp, (col, stresses) in Exps.items():\n",
" # plot analytical curves\n",
" if temp >= 25:\n",
- " ax.plot(sigs, BGRa(sigs, temp), color=col, ls='--')\n",
- " ax.plot(sigs, PLLC(sigs, temp), color=col, ls='-')\n",
+ " ax.plot(sigs, BGRa(sigs, temp), color=col, ls=\"--\")\n",
+ " ax.plot(sigs, PLLC(sigs, temp), color=col, ls=\"-\")\n",
"\n",
" # simulation in ogs and plot results\n",
" eps_dot = []\n",
@@ -156,14 +282,15 @@
" ogs_model.write_input()\n",
" # hide output\n",
" with contextlib.redirect_stdout(None):\n",
- " ogs_model.run_model(logfile=f\"{out_dir}/out.txt\", \n",
- " args=\"-m \" + f\"{data_dir}/Mechanics/PLLC/\")\n",
+ " ogs_model.run_model(\n",
+ " logfile=f\"{out_dir}/out.txt\", args=\"-m \" + f\"{data_dir}/Mechanics/PLLC/\"\n",
+ " )\n",
" pvdfile = vtuIO.PVDIO(f\"{prj_name}.pvd\", dim=3)\n",
" eps_zz = pvdfile.read_time_series(\"epsilon\", points)[\"pt0\"][:, 2]\n",
" eps_zz_dot = np.abs(np.diff(eps_zz)) / np.diff(pvdfile.timesteps)\n",
" # omit the first timestep\n",
" eps_dot += [np.mean(eps_zz_dot[1:])]\n",
- " ax.loglog(1e-6*stresses, eps_dot, 'o', c=col, markeredgecolor=\"k\")\n",
+ " ax.loglog(1e-6 * stresses, eps_dot, \"o\", c=col, markeredgecolor=\"k\")\n",
"\n",
"# plot experimental data points\n",
"for Ex, (temp, m, Data) in ExData.items():\n",
@@ -171,20 +298,23 @@
" ax.loglog(stresses, eps_dot, m, c=Exps[temp][0])\n",
"\n",
"# create legend\n",
- "patches = [mpl.patches.Patch(color=col, label=str(temp) + '°C')\n",
- " for temp, (col, _) in Exps.items() if temp >= 25][::-1]\n",
- "addLeg = lambda **args : patches.append(mpl.lines.Line2D([], [], **args))\n",
- "addLeg(c='k', label='PLLC')\n",
- "addLeg(c='k', ls='--', label='BGRa')\n",
- "addLeg(c='w', ls='None', marker='o', mec=\"k\", label='OGS')\n",
- "addLeg(c='k', ls='None', marker='s', label='DeVries (1988)')\n",
- "addLeg(c='k', ls='None', marker='^', label='WIPP CS')\n",
- "addLeg(c='b', ls='None', marker='P', label='Bérest (2017) 7.8°C')\n",
- "addLeg(c='orange', ls='None', marker='P', label='Bérest (2015) 14.3°C')\n",
- "ax.legend(handles=patches, loc='best')\n",
+ "patches = [\n",
+ " mpl.patches.Patch(color=col, label=str(temp) + \"°C\")\n",
+ " for temp, (col, _) in Exps.items()\n",
+ " if temp >= 25\n",
+ "][::-1]\n",
+ "addLeg = lambda **args: patches.append(mpl.lines.Line2D([], [], **args))\n",
+ "addLeg(c=\"k\", label=\"PLLC\")\n",
+ "addLeg(c=\"k\", ls=\"--\", label=\"BGRa\")\n",
+ "addLeg(c=\"w\", ls=\"None\", marker=\"o\", mec=\"k\", label=\"OGS\")\n",
+ "addLeg(c=\"k\", ls=\"None\", marker=\"s\", label=\"DeVries (1988)\")\n",
+ "addLeg(c=\"k\", ls=\"None\", marker=\"^\", label=\"WIPP CS\")\n",
+ "addLeg(c=\"b\", ls=\"None\", marker=\"P\", label=\"Bérest (2017) 7.8°C\")\n",
+ "addLeg(c=\"orange\", ls=\"None\", marker=\"P\", label=\"Bérest (2015) 14.3°C\")\n",
+ "ax.legend(handles=patches, loc=\"best\")\n",
"\n",
"fig.tight_layout()\n",
- "plt.show()\n"
+ "plt.show()"
]
},
{
@@ -197,7 +327,7 @@
"\n",
"Zill, Florian, Wenqing Wang, and Thomas Nagel. Influence of THM Process Coupling and Constitutive Models on the Simulated Evolution of Deep Salt Formations during Glaciation. The Mechanical Behavior of Salt X. CRC Press, 2022. https://doi.org/10.1201/9781003295808-33.\n",
"\n",
- "Li, Shiyuan, and Janos Urai. Numerical Studies of the Deformation of Salt Bodies with Embedded Carbonate Stringers. Online, print. Publikationsserver der RWTH Aachen University, 2012. http://publications.rwth-aachen.de/record/211523/files/4415.pdf "
+ "Li, Shiyuan, and Janos Urai. Numerical Studies of the Deformation of Salt Bodies with Embedded Carbonate Stringers. Online, print. Publikationsserver der RWTH Aachen University, 2012. http://publications.rwth-aachen.de/record/211523/files/4415.pdf"
]
}
],
diff --git a/Tests/Data/Notebooks/SimplePETSc.ipynb b/Tests/Data/Notebooks/SimplePETSc.ipynb
index f3990b10e54..34c1fc3519a 100644
--- a/Tests/Data/Notebooks/SimplePETSc.ipynb
+++ b/Tests/Data/Notebooks/SimplePETSc.ipynb
@@ -5,11 +5,12 @@
"id": "bb0907b4-4e26-4c4e-ab1f-22b5330cb1d2",
"metadata": {},
"source": [
+ "+++\n",
"title = \"SimplePETSc\"\n",
"date = \"2021-11-09\"\n",
"author = \"Lars Bilke\"\n",
"web_subsection = \"elliptic\"\n",
- ""
+ "+++\n"
]
},
{
@@ -24,15 +25,17 @@
"cell_type": "code",
"execution_count": 8,
"id": "7962f42f-fd53-4fc1-b966-a8ba924aca6c",
- "metadata": {},
+ "metadata": {
+ "lines_to_next_cell": 2
+ },
"outputs": [],
"source": [
"import os\n",
"\n",
"prj_name = \"square_1e1_neumann\"\n",
- "data_dir = os.environ.get('OGS_DATA_DIR', '../../../Data')\n",
+ "data_dir = os.environ.get(\"OGS_DATA_DIR\", \"../../../Data\")\n",
"prj_file = f\"{data_dir}/EllipticPETSc/{prj_name}.prj\"\n",
- "out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', '_out')\n",
+ "out_dir = os.environ.get(\"OGS_TESTRUNNER_OUT_DIR\", \"_out\")\n",
"\n",
"if not os.path.exists(out_dir):\n",
" os.makedirs(out_dir)\n",
@@ -42,6 +45,7 @@
"! mpirun -np 2 ogs {prj_file} > out.txt\n",
"\n",
"from datetime import datetime\n",
+ "\n",
"print(datetime.now())"
]
},
@@ -79,10 +83,11 @@
"\n",
"pvdfile = vtuIO.PVDIO(f\"{prj_name}.pvd\", dim=2)\n",
"time = pvdfile.timesteps\n",
- "points={'pt0': (0.3,0.5,0.0), 'pt1': (0.24,0.21,0.0)}\n",
+ "points = {\"pt0\": (0.3, 0.5, 0.0), \"pt1\": (0.24, 0.21, 0.0)}\n",
"pressure_linear = pvdfile.read_time_series(\"pressure\", points)\n",
"\n",
"import matplotlib.pyplot as plt\n",
+ "\n",
"plt.plot(time, pressure_linear[\"pt0\"], \"b-\", label=\"pt0 linear interpolated\")\n",
"plt.plot(time, pressure_linear[\"pt1\"], \"r-\", label=\"pt1 linear interpolated\")\n",
"plt.legend()\n",
diff --git a/Tests/Data/Notebooks/testrunner.py b/Tests/Data/Notebooks/testrunner.py
index 27710943c0b..3cc73eb667f 100644
--- a/Tests/Data/Notebooks/testrunner.py
+++ b/Tests/Data/Notebooks/testrunner.py
@@ -11,61 +11,49 @@
import toml
from pathlib import Path
import jupytext
+import subprocess
def save_to_website(exec_notebook_file, web_path):
- from nb2hugo.writer import HugoWriter
-
- output_path = "docs/benchmarks"
+ output_path_arg = ""
notebook = nbformat.read(exec_notebook_file, as_version=4)
first_cell = notebook.cells[0]
- if is_jupytext:
- if "Tests/Data" not in exec_notebook_file:
- output_path = str(Path(exec_notebook_file).parent.parent)
- else:
- lines = first_cell.source.splitlines()
- toml_begin = lines.index("+++")
- toml_end = max(loc for loc, val in enumerate(lines) if val == "+++")
- toml_lines = lines[toml_begin + 1 : toml_end]
- parsed_frontmatter = toml.loads("\n".join(toml_lines))
- output_path = (
- Path(build_dir)
- / Path("web/content")
- / Path(output_path)
- / Path(parsed_frontmatter["web_subsection"])
- )
- elif first_cell.cell_type == "raw":
+ if "Tests/Data" in exec_notebook_file:
lines = first_cell.source.splitlines()
- last_line = lines[-1]
- if "" not in last_line:
- print(
- f"Warning: {exec_notebook_file} does not contain '' as the "
- "last line in the RAW cell!"
- )
- parsed_frontmatter = toml.loads("\n".join(lines[:-1]))
- if "web_subsection" not in parsed_frontmatter:
- print(
- f"Error: {exec_notebook_file} frontmatter does not contain "
- "'web_subsection'!"
- )
- output_path = os.path.join(output_path, parsed_frontmatter["web_subsection"])
- output_path = Path(build_dir) / (Path("web/content") / Path(output_path))
- else:
- print(
- f"Warning: {exec_notebook_file} does not contain a RAW cell as its first "
- "cell!"
+ toml_begin = lines.index("+++")
+ toml_end = max(loc for loc, val in enumerate(lines) if val == "+++")
+ toml_lines = lines[toml_begin + 1 : toml_end]
+ parsed_frontmatter = toml.loads("\n".join(toml_lines))
+ output_path = (
+ Path(build_dir)
+ / Path("web/content/docs/benchmarks")
+ / Path(parsed_frontmatter["web_subsection"])
)
- output_path = os.path.join(output_path, "notebooks")
- writer = HugoWriter()
- writer.convert(
- exec_notebook_file,
- web_path,
- output_path,
- os.path.join(
- os.path.dirname(os.path.abspath(__file__)),
- "nbconvert_templates/collapsed.md.j2",
- ),
+ output_path_arg = (
+ f"--output-dir={Path(output_path) / Path(exec_notebook_file).stem}"
+ )
+
+ template = os.path.join(
+ os.path.dirname(os.path.abspath(__file__)),
+ "nbconvert_templates/collapsed.md.j2",
)
+ subprocess.run(
+ [
+ "jupyter",
+ "nbconvert",
+ "--to",
+ "markdown",
+ f"--template-file={template}",
+ "--output=index",
+ output_path_arg,
+ exec_notebook_file,
+ ],
+ check=True,
+ )
+
+ if not "Tests/Data" in exec_notebook_file:
+ return
+
for subfolder in ["figures", "images"]:
figures_path = os.path.abspath(
os.path.join(os.path.dirname(notebook_file_path), subfolder)
@@ -140,7 +128,7 @@ def save_to_website(exec_notebook_file, web_path):
nb = nbformat.read(f, as_version=4)
ep = ExecutePreprocessor(kernel_name="python3")
- # 1. Run the notebook
+ # Run the notebook
print(f"[Start] {notebook_filename}")
start = timer()
try:
@@ -162,15 +150,7 @@ def save_to_website(exec_notebook_file, web_path):
pass
end = timer()
- # 2. Instantiate the exporter. We use the `classic` template for now; we'll get
- # into more details later about how to customize the exporter further.
- html_exporter = HTMLExporter()
- html_exporter.template_name = "classic"
-
- # 3. Process the notebook we loaded earlier
- (body, resources) = html_exporter.from_notebook_node(nb)
-
- # 4. Write new notebook
+ # Write new notebook
with open(convert_notebook_file, "w", encoding="utf-8") as f:
repo = "https://gitlab.opengeosys.org/ogs/ogs"
branch = "master"
@@ -178,19 +158,36 @@ def save_to_website(exec_notebook_file, web_path):
repo = os.environ["CI_MERGE_REQUEST_SOURCE_PROJECT_URL"]
branch = os.environ["CI_MERGE_REQUEST_SOURCE_BRANCH_NAME"]
+ # Check frontmatter has its own cell
+ first_cell = nb["cells"][0]
+ if (
+ first_cell.cell_type == "markdown"
+ and first_cell.source.startswith("+++")
+ and not first_cell.source.endswith("+++")
+ ):
+ print(
+ f"Error: {notebook_filename} notebook metadata is not a separate cell (in markdown: separate by two newlines)!"
+ )
+ success = False
+
+ # Check second cell is markdown
+ second_cell = nb["cells"][1]
+ if second_cell.cell_type != "markdown":
+ print(
+ f"Error: {notebook_filename} first cell after the frontmatter needs to be a markdown cell! Move the first Python cell below."
+ )
+ success = False
+
# Modify metadata
- meta_cell = nb["cells"][0]
- if is_jupytext:
- if meta_cell.source.startswith("---"):
- print(
- f"Error: {notebook_filename} frontmatter is not in TOML format! Use +++ delimitiers!"
- )
- success = False
- meta_cell.source = meta_cell.source.replace(
- "+++\n", "+++\nnotebook = true\n", 1
+ first_cell = nb["cells"][0]
+ if first_cell.source.startswith("---"):
+ print(
+ f"Error: {notebook_filename} frontmatter is not in TOML format! Use +++ delimitiers!"
)
- else:
- meta_cell.source = f"notebook = true\n{meta_cell.source}"
+ success = False
+ first_cell.source = first_cell.source.replace(
+ "+++\n", "+++\nnotebook = true\n", 1
+ )
# Insert Jupyter header with notebook source and binderhub link
binder_link = f"https://mybinder.org/v2/gh/bilke/binder-ogs-requirements/master?urlpath=git-pull%3Frepo={repo}%26urlpath=lab/tree/ogs/{notebook_file_path_relative}%26branch={branch}"
@@ -216,41 +213,10 @@ def save_to_website(exec_notebook_file, web_path):
src="https://img.shields.io/static/v1?label=&message=Launch notebook&color=5c5c5c&logo=data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAABwAAAAcCAYAAAByDd+UAAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAABmJLR0QA/wD/AP+gvaeTAAAACXBIWXMAAC4jAAAuIwF4pT92AAAAB3RJTUUH4gsEADkvyr8GjAAABQZJREFUSMeVlnlsVFUUh7/7ZukwpQxdoK2yGGgqYFKMQkyDUVBZJECQEERZVLQEa4iKiggiFjfqbkADhVSgEVkETVSiJBATsEIRja1RoCwuU5gC7Qww03Zm3rzrH/dOfJSZUm4y6Xt9957vnnN/55wruI7RVjMNQAA3AiX6bxw4BTQAQQDvnF1pbYjrAAEUAmXADGAQ0AOQwCWgHqgGdgCRdNBrAm2wW4A1wN2ACZwG/gbcQBFwg/Z2I/AS0JoKanQzmoXAamA0cBx4EhgDTAYmAvcArwNhYD6wHHDbNts9D20LlgMrgWPAXKAO/j8rPc8A5uiNAUwH9tjnddfDAn1mFkJWyoRR58hsv8KIfraAz/QvC3golf2UwEBZBYGyCoJfj/LFz/ceDxRJ09Hccbz/6dDu0ozg7lICZRVXrNFQEyWaDmAkkNslMAnSE59x9IrsMVt8awBP4rI3P9acs83hC3+BkFMAd2eoHn8BrdpG77RA2+IiYDPwHnAbEAOkMGQMcAKTdNheBXqmgDoBhw6xda2Q9tGHPhE4hRTlrrxQGRB29IqE3IUtTyDFu9rQC8AiwAiUVdgFNhTIA85oT68G2nb5ODABJf25niL/emfexX1AA0IWeIr8xWbY+yKwBJVzC4FSm71MlFIdwH505UnnYT5KWRawCvgp0eYBCKEqSBwpFuVMqp2a5Q1WO6TcakiZ55DWwyVVKxDC8gLPA1OAJh32q8qcHTgEKEbl2ncAua99lPy2FdgskH2FlFXNI8IVewcO8P+WUyjr8vqPfmvt+plhmVltIJeilLoK+CWVopy250LAgyrELcl/9nB/ixkbF3GKyOJ/rJs8hxNDZx1KDFvsz+9jJvINAQz1EKvxR7OddzrroyXGiRV5zvp1WPlSzN7bJVCmEtKDF38khguQeR5iBRYGFoaZaUUv9YsEc+KGYfq9vssN1qDsP2MDHRZiYBRXpoEMwa1XAe3Gm4A2YDDQ1z7JTbyvG3O1hXEvcNI0xFPzTh5ZueB4HeXH6hoGR1onC2SlhQgD5RnEl7kwXTOqfu4SeBT4Q5/jVIBtL29KfnsUGAecsISY++W+mpohwQujXJYlPAnzh2HBc7Uxw1iGSpU2VAu7C6Az1A68gEr4ZI6NXT78Pkxh9JEwU4JlGsYbO3a+c7g50/esFGIqcBb4fEzgNBlWwgI2AVsAH13V0oL1K5LvNcBOYACwsfb7qiX3n2mcmGXGirPjHf8uPHqw/Xy/IeuAV/TG3gaOAGyfPwJUbm4HosAdpKilzk7vIVT1iAPTTWG8Of5MY/vIFn8Pt2UVZkfbqi0hvFrFlcBaQNo2DKoxt6CqjQ84nzKktkV+YIE+hz1OaUVyou0iKx41BAR02KYB7wMdnWBJm4aOgOz8MWUDTpa6/NazGdUlo8c2ZuVukdBWfOnCtHlffXAwdPsEK2o47Ju0i2MysAt1xxkLtOpwpwzpFd4+sOHXKHDAIa16YNTJrJzS3x9ZVdvoy+WbecNTLfUCs7Xd/aQr3umGy0rgshIhQ8pNhpSmIeVzTZm9pnjNuLDLXT97gKdRKXUWXUvt3qUNqX1oYz2Bj1H3mXPABh22JlRnuBl4DHWPAVgKfAjIzkDntYB6hIHFKPXO0gbLUQp0oO49Xv1eCXySCtYtDzt56kU159moQulDqfEccAD4FDgEJFLBrgtog4I6r36oG0IC1d0DqNZEOhjAfzgw6LulUF3CAAAAJXRFWHRkYXRlOmNyZWF0ZQAyMDE4LTExLTA0VDAwOjU3OjQ3LTA0OjAwLtN9UwAAACV0RVh0ZGF0ZTptb2RpZnkAMjAxOC0xMS0wNFQwMDo1Nzo0Ny0wNDowMF+Oxe8AAAAASUVORK5CYII=" />
"""
text += f"""
\n\n"""
+ second_cell.source = text + second_cell.source
- for cell in nb["cells"]:
- # Check frontmatter has its own cell
- if (
- cell.cell_type == "markdown"
- and cell.source.startswith("+++")
- and not cell.source.endswith("+++")
- ):
- print(
- f"Error: {notebook_filename} notebook metadata is not a separate cell (in markdown: separate by two newlines)!"
- )
- success = False
- # Get first regular markdown cell
- if cell.cell_type == "markdown" and not cell.source.startswith("+++"):
- first_markdown_cell = cell
- break
-
- first_markdown_cell.source = text + first_markdown_cell.source
nbformat.write(nb, f)
- # 5. Symlink images or figures subfolder
- for subfolder in ["figures", "images"]:
- figures_path = os.path.abspath(
- os.path.join(os.path.dirname(notebook_file_path), subfolder)
- )
- symlink_figures_path = os.path.join(notebook_output_path, subfolder)
- if os.path.exists(figures_path) and not os.path.exists(
- symlink_figures_path
- ):
- print(
- f"{subfolder} folder detected, symlink {figures_path} to "
- f"{symlink_figures_path}"
- )
- os.symlink(figures_path, symlink_figures_path)
-
status_string = ""
if notebook_success:
status_string += "[Passed] "
diff --git a/Tests/Data/Notebooks/thermo-osmosis.run-skip.ipynb b/Tests/Data/Notebooks/thermo-osmosis.run-skip.ipynb
index 43ce56d5cfa..2f57bd7e3ae 100644
--- a/Tests/Data/Notebooks/thermo-osmosis.run-skip.ipynb
+++ b/Tests/Data/Notebooks/thermo-osmosis.run-skip.ipynb
@@ -1,284 +1,317 @@
{
- "cells": [
- {
- "cell_type": "raw",
- "metadata": {},
- "source": [
- "author = \"Jörg Buchwald\"\n",
- "date = \"2022-05-27T12:39:58+01:00\"\n",
- "title = \"Thermo-Osmosis in a one-dimensional column\"\n",
- "weight = 70\n",
- "web_subsection = \"thermo-hydro-mechanics\"\n",
- ""
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "tags": []
- },
- "source": [
- "## Problem description\n",
- "\n",
- "The problem describes a one-dimensional column at $T$=300 K in sudden contact with a temperature reservoir at one side at $T_1$ = 350 K.\n",
- "\n",
- "Thermo-osmotic and filtration effects are described by contributions to the hydraulic flux $J^w$\n",
- "\\begin{equation}\n",
- "J^w=-\\rho_w \\frac{\\mathbf{k}}{\\mu}\\left(\\nabla p-\\rho_w \\mathbf{g} \\right)-\\rho_w \\mathbf{k}_{pT} \\nabla T,\n",
- "\\end{equation}\n",
- "and the conductive heat flux $I$\n",
- "\\begin{equation}\n",
- "I=- \\mathbf{\\lambda}_s (1-\\phi)+\\mathbf{\\lambda}_w \\phi)- \\mathbf{k}_{Tp} \\nabla p,\n",
- "\\end{equation}\n",
- "\n",
- "where $\\mathbf{k}_{pT}$ is the phenomenological coefficient of thermo-osmosis and $\\mathbf{k}_{Tp}$ the phenomenological coefficient of thermo-filtration. \n",
- "It can be shown that $\\mathbf{k}_{Tp}=T*\\mathbf{k}_{pT}$ (Zhou et al. 1998)."
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "## Get benchmark results\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 1,
- "metadata": {},
- "outputs": [
- {
- "name": "stderr",
- "output_type": "stream",
- "text": [
- "/home/buchwalj/.local/lib/python3.10/site-packages/vtuIO.py:147: PerformanceWarning: DataFrame is highly fragmented. This is usually the result of calling `frame.insert` many times, which has poor performance. Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`\n",
- " df[\"r_\"+str(i)] = (df[x]-val[x]) * (df[x]-val[x]) + (df[y]-val[y]) * (df[y]-val[y])\n"
- ]
- }
- ],
- "source": [
- "import os\n",
- "import vtuIO\n",
- "import numpy as np\n",
- "\n",
- "filename = \"expected_Column_ts_68_t_7200000.000000.vtu\"\n",
- "data_dir = os.environ.get('OGS_DATA_DIR', '../../Data')\n",
- "file = {}\n",
- "file[\"THM\"] = f\"{data_dir}/ThermoHydroMechanics/Linear/ThermoOsmosis/{filename}\"\n",
- "file[\"TR\"] = f\"{data_dir}/ThermoRichardsFlow/ThermoOsmosis/{filename}\"\n",
- "file[\"TRM\"] = f\"{data_dir}/ThermoRichardsMechanics/ThermoOsmosis/{filename}\"\n",
- "x=np.array([i*0.1 for i in range(200)])\n",
- "r = np.array([[i,0.5,0.0] for i in x])\n",
- "resp = {}\n",
- "respvars = [\"temperature\", \"pressure\"]\n",
- "for model in file:\n",
- " resp[model] = {}\n",
- " f = vtuIO.VTUIO(file[model], dim=2)\n",
- " for var in respvars:\n",
- " if \"M\" in model:\n",
- " resp[model][var] = f.get_set_data(f\"{var}_interpolated\",pointsetarray=r)\n",
- " else:\n",
- " resp[model][var] = f.get_set_data(f\"{var}\",pointsetarray=r)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Read-in the analytical solution\n",
- "\n",
- "An analytical solution was provided by Zhou et al. 1998 and can be obtained via [github](https://github.com/joergbuchwald/thermo-osmosis_analytical_solution).\n",
- "For this example we used $\\mathbf{k}_{pT}=2.7e-10\\, m^2/(s K)$ and a fully saturated material. More details on model parameters can be found in the corresponding project files.\n",
- "The Thermo-Richards (TR) model uses a correction to account for mechanical effects in the mass-balance equation. See Buchwald et al. 2021 for further details."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 3,
- "metadata": {},
- "outputs": [],
- "source": [
- "import zhou_solution_thermo_osmosis\n",
- "aTO = zhou_solution_thermo_osmosis.ANASOL(0,50,100)\n",
- "aNoTO = zhou_solution_thermo_osmosis.ANASOL(0,50,100)\n",
- "aNoTO.Sw = 0\n",
- "t=7.2e6"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Plot temperature and pressure along the column"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 9,
- "metadata": {},
- "outputs": [
- {
- "data": {
- "image/png": "iVBORw0KGgoAAAANSUhEUgAAAw8AAAJ/CAYAAAA+ie+jAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAACcZUlEQVR4nOzdd3gVZfr/8feTkB4g9E4ogRBIQ0Cw0hT9KiCiC6tr3UVQ97cWlC1WRF0Liq4F6wrqrqtRgdWsZa1rARSUmARilCBI7wFCevL8/jjFk5CEEJLMSfJ5Xde55pyZZ2bumZRz7vM0Y61FRERERETkaAKcDkBERERERJoGJQ8iIiIiIlIrSh5ERERERKRWlDyIiIiIiEitKHkQEREREZFaUfIgIiIiIiK1ouRBRERERERqRcmDiIiIiIjUipIHEZE6Msb0McZY92OM0/HIsfH52V3hdCwiIk2FkgcREREREakVJQ8iIiIiIlIrSh5ERERERKRWlDyIiBwjY8wYY4wFfvJZ/YlPG3prjFlcxX5TjTGpxpidxphiY8wuY8x/jTFXGGMCqyj/aaVjbnSvn2KM+dIYc8gYs8MYs9QYk+izX39jzD/d2/KNMd8YYy6u5lps5biNMRHGmLuMMRnuc+QZY742xsw0xtT4vmGM6WaMmW+MWeve77AxJtsY86Qxpn9197LS4wpjTBtjzF+NMd8bYwqq6p9gjEk0xswzxqw0xuw3xpQYY/YYYz4xxlxljGlVxfm8/VR8Vi+qdP5P3WXvqRxbFcfbU6nM3Pq6Pvf+ocaYG4wxy93XWGSM2WyMeU39bETECcbaI/4XiohIDYwx7YBhQBfgH+7VNwPf+RTbZq1d5y4fDrwCnAfkAI/gSjx6A9cCCcDnwGRrba7PeYYB7YDLgEuBTcATwDnAc4AFpgNTgEPA6UAh8ALwErADOA24AdeXRb+z1r5Q6VrOcD99GEgE/g3EALuAF4E9QJL7+toBqcAF1triKu7L2UAKEAEsBv4DlANnADPd8V5hrf1XFfcS4AOfeznDfU/eA4KBP7vjuNJau9h9b1a7j/kPd7lDQDRwhfuYnwLnWGsLfM4XCpxa6Xzzgf/6XMp+a+037mSnLzABmANgrTWVrvl0d3z/wPX7cJe1du7xXp97337AO0As8AmwCNjnLncj0BF4DLjB6s1cRBqLtVYPPfTQQ486PIA+uD68WmBMDeX+5S7zI9Cu0rYw4Gv39iXV7D/XvT0f+BAI9NlmcH0otcDbwFKgS6X973Rv3+67b6Uyn/pcyyu4v1zy2T4YOOzePr+K/ZOAAvf2P1Sx/QL3tmJgWDUxeM5/ENeHaN9tQ/gl+QAY7n79p2qOtdS9/aEafi7W95g1lLvCU7aGMhvdZebW4ny1ub7WwA/uda9X8fPoB+S6t1/n9N+CHnro0XIearYkItKA3E1Lfu1+eZu1dr/vduv6VvwW98vzjTEn1nC4MOCv1toyn/0t8Ib75UQgw1q7s9J+Ke5lV1wf8mtSBtzkPq5vnOuAZ90vbzDGdK+039+AUCALV+1IBdbaN4EvgSDg7qPEsMFau6jSuizgKvcxALYBdwELqznG39zLmVU1X3JYba7vJmAAUEoVNQvW2g38cp/vNMaENWC8IiJeSh5ERBrW79zLYlw1A1X5DChxP/91NWXA9cH+yyrWb/J5/r8qtm/0ed63huMDpFlrt1ez7d/uZSvgV56V7uY9o90vl1T+oOvjI/fybGNMVA0xvFd5hbW23Fr7vLX2R/frbdbaudbaQz5xBLn7CITiqmUB1zf4A2s4lxOOen388nvzlbV2azXH8dzP9sCZ9RyjiEiV/O3bGBGR5uYU9zIHKHd/sK3KXlw1A8Oq2Q6w21pbVMX6Qp/nWypvtNYWGONtqh9Rc7gVOoFXluXzfITP81N8nq+r4Rp3uZcGGIqrHX9VNtYUoIcxpj1wPTAJGISrZqYqHWpzvEa0saaNxpjeQE/3y8wa7ucen+fDgLeOPzQRkZopeRARaVjd3Ms4XH0CjqZLDduqShyOtczRapwP17DNt8mVb5zdfJ7/8yjHr2r/ygpr2AaAMWYwrm/eu+LqG3Cze5nnE9MS9/MjRrJy2NGuz/d+znI/jqam+ykiUm+UPIiINCxPE55MXG3aj6Y2CYJTbDXPfd0AfFWLY+UcZywv4UocNuHqgJ3nu9EY0+c4j+8vXgKeqkW53Q0diIgIKHkQEWlo23GNjBNgrV3pdDC1EF7DtvY+z3f5PN/m+7yhr9M9hKmnedc/KicODcCbKBljTDV9Ouqzw7Lv/cxrIr83ItJCqMO0iEjD8nRw7m+MCa6ukDGmpzHmamPMuEaKqzr9atg22Of5Kp/nX1ZT5gjGNcHd1caYNnUJzq2rz/PqOhMfrW/HsfBtyhVZeaO7T0LH+jqZtXYzsNn98mj380T3/RxUX+cXEamJkgcRkbrznSjN267eGHOaMeZpY8wI4Hn36hBcHXurcyOu5il96jvIY5RsjOlazbYp7mUprrkHAO+woZ7Oz1ONT+9sX8aYDsCrwF9wTehWVzt8nsdUUyaxmvW+PCNc+f7sBrt/dmf7lPMdzaqq0arGUP/vp57fm1Nq+HmAa6jahbh+JiIiDU7Jg4hI3e3kl07Qvh/wJuLq5NraWvsZv3QifsAY06nyQYwxpwC/x9Xh9x+VtzeyYuChyiuNMUNwzYoM8Ki1dlulItfjuheJuBKhyvsH4pr5OgS4o4bhXI/Knayku19eaYzx7WDsmdH79locaqN76fuzG43rZ+c7j0UartGwAKZWOlcQrnk6atMZ/lgsAL7HNS/GQvf9q8AYczOu5lsvWmvX1/P5RUSqpD4PIiJ1ZK0tM8a8AVwK/NEYU4hr1JurcX1b/YW76FW4PjRfCHxnjHkC14ffNsBpwG9xfZs+2Vrrrc1wjyjUnV+aEoUaY85wn/tDY0w7XB8efSd+O8UYEwN8g+sD7amVwh7sPsYG94fwyt4EOhpjPgIWA/vcx78ZV3+IVODWKu5FhjHmPFw1Eg8bY07FNdpRrjv+mbhmUX7QWvuizzWG1hAjwBfW2qpGJ5qBa7SldsAaY8wCIBvoDfyBis2WhrkniltbaQ6LV3ElGbOMMRuAYOA2XDNAe+fksNaWGGPmAo8Dd7jnqPjMfY5r3PesNxAN9HPHvt9a+01dr89am2eM+T9c9/t84GtjzPO4fq+64UpizgE+Bv5fFfdHRKRBmOP48kdEpMVzt91/AJgMdMZVG/E5cIu19qdKZSfhmvxrJK428gW4vl3+N/CY74Rn7vKLgcurOq+11rhnr65uroSxuL5Zr27ehrustXN9zvUprm/dX3THeD1wCa5mQQHAOuDvwHPW2vJqjokxpjOuEZfOxZU0hOC6J8uBJ901Mb7l+9QQI0Bfa+3Gas7VH1cicyau2gPP/Xwd+A+wttIuV1prF/vsH4xrlupfAz1wJUqrcf3s0ivtizHmElz3ZQiuJk/fAQ9Za98yxmzElTx4/M9aO+Z4rs99zhBcidI0IAHXpHe5uGpDXsbVYbzan4eISH1T8iAiIhWSB2vtFc5GIyIi/kp9HkREREREpFaUPIiIiIiISK2ow7SISAvm03G3nXvZzb2uwFr7ZTW7iYhIC6U+DyIiLZgxpro3gU3W2j6NGYuIiPg/JQ8iIiIiIlIrarZUT84++2z73nvvOR2GiIiIiDR/xqkTq8N0PdmzZ4/TIYiIiIiINCglDyIiIiIiUitKHkREREREpFaUPIiIiIiISK0oeRARERERkVpR8iAiIiIiIrWi5EFERERERGpF8zyIiLQABw8eZNeuXZSUlDgdioiIVCMoKIjOnTvTpk0bp0OplpIHEZFm7uDBg+zcuZMePXoQFhaGMY7NLSQiItWw1lJQUMDWrVsB/DaBULMlEZFmbteuXfTo0YPw8HAlDiIifsoYQ3h4OD169GDXrl1Oh1MtJQ8iIs1cSUkJYWFhTochIiK1EBYW5tdNTJU8iIi0AKpxEBFpGvz9/7WSBxERERERqRUlDyIiIiIiUitKHkREREREpFaUPIiISJPRp08fBg8eTHJyMsnJyXTt2hVjzBHrevbsSXJyMsYY2rdvT3JyMnv37vUe5/zzz6d3794YY4iJieHqq68GYP78+d79goOD2bx5c7WxvPvuuxWOn52d3eDXLyLiNL9LHowxfYwx1xlj3jPGbDfGlBhjDhlj0owxdxtj2lWxzxXGGHuUR2QN55xkjPnEGJPrPtdKY8zlDXulIiJSF++88w5paWmkpaV5P/RXXjdjxgzS0tIAmDx5MmlpaXTo0MF7jKVLlzJv3jwAnn/+eZ5++mkA5syZ492vtLSUhx9+uNo4HnzwwQrHj42Nre9LFRHxO36XPAArgUeBtcBkIAaYAGQAtwHfGGM6VrFfAZBdw6O8qpMZY24H3gL2AWOAE4E0YLEx5rn6uSQREakPo0ePPuqwszExMcTExBz3uSZPnsxzzz3Hnj17jti2fPly2rU74rssEZFmzx+TB4CF1tqbrLWrrLWbrLUrrLWXAv8D+gJXV7HP19baQTU88ivvYIwZDcwD1gDTrLVp1tosa+3VwNvADGPMZQ14nSIiTdYr73/f6Od88cUX6dKlS41lLrnkEi655JLjPtdf/vIX8vPzeeyxx47Ydt999/GnP/3puM8hItLU+GPyMAt4oJpt37iXnevpXHe6l49Za8sqbVvgXt5RT+cSEWlW/vXf5t3Gf+TIkYwePZonnniCQ4cOeddnZGSQn5/PyJEjHYxORMQZrZwOoDJr7b+rWm9cM2ac6H758fGexxjTGRjtfvlRFUW+BIqA/saYYdbab6ooIyIifu6tt94iOTn5iPX79u076r5//vOf+b//+z+eeeYZbr75ZgDuv/9+/vjHP9Z3mCIiTYI/1jxUYIwJMcYkAS8DJwH3WmuXVVG0tTHmTmPMN8aYXcaYLe5O15cYY6q6zmG4rv+wtfaI4TSstSXABvfLEUeLc1uuaxSP3PQMtiypKjwREXGCp0Nz5Yenw3RNzj77bJKTk1mwYAFFRUVs2LCBH374gbPOOqsRIhcR8T9+nTwYY1YAhbg6MA8CTrHW3lZN8RNw1Uzciqvj8wwgEFfS8bYxJrhS+f7u5c4aQtjuXvY7WqzlZQXkpmeQPX8BkTH9j1ZcRESaiD/96U9s376dF198kfnz53trIEREWiK/Th6AacAQYCqQByx3D9daOe51wE3W2nOtte9Za9dZa98DzgZWA+cA91bap417eURHah8F7mXbowXaPgi+e/ppYufMJiox4WjFRUSkifjVr35F//79uffee/n888+58MILnQ5JRMQxfp08WGs3uxOBpcA4YAWu4VofrFTua2vtgir2LwP+6n75e2NM6DGGYDyHqnKjMTONMauNMasBNsdGMTNrISmZqcd4GhGRpueiCS1jXoPAwEBuvvlmfv75Z/7whz8QGBjodEgiIo7x6+TBl7W2HLjL/fIPxpioWu76rXsZBgz1WX/QvQyvYV9PsnGwqo3W2mettcOttcMBurQN4tm4a5kWP7GWoYmINF0XnzXI6RAazZVXXskHH3zAFVdc4XQoIiKOajLJg1uGexlMxUSgJr59Gnxn9MlxL2saMLybe7mhhjJerTsXkvXQI+SmZxy9sIiIHJeRI0d6Z4Y+55xzuPfeX1qnPv/8894RljyjLe3du9e7/fzzz+eOO1wjcc+YMcM7U7XvfsnJySxatAiAkJAQzjjjDEJCQgBYtGjREcfftm1bg12riIi/MNZW2SLHEcaYWGCktfalara3AQ64X55trX3fGBMGjAc+sdYermKfaGCj++Up1trl7vWdcXWIDgB6Vx5xyRgThKvGIRQYYa1dXVPssTHd7SsLf0MwY2mXV0zPqVNqdc0iIg0tKyuLuLg4p8MQEZFaqsX/bVPTxobkbzUPJwGLauibMNjn+Xr3sguu2aCrG07VU0NRhGvUJgCstbuAz90vx1ex3ym4EoefjpY4AHgGczpUvl+Jg4iIiIg0S/6WPIArphnVbLvdvVxhrc2ptO3SyoXdozL92f3yOWtt5ZGVPH0orjPGVO4Bd6N7efSBwIFWQWEAlBZtrE1xEREREZEmx9+Sh1L38mFjzHxjzChjTF9jzJnGmP/iGnJ1C3C5zz5l7uVvjTEvGGNOMcb0NsachqtGYiTwKXDEdKDW2k9wJRBDgRRjTJIxJs4Y8xQwGVhsrV1cm8BDQsMoLzeEh+ZycP/RZy0VEREREWlq/Cp5sNb+AxgLvABMAP4L/Ai8jmtehtuBeGvtjz77bAZigbtxTST3Nq4Ozv/GNZLSTOAMa20BVbDWzgWmAB2Az4BVuCac+6219sraxm4CAsgv7Igx8NP3abW+ZhERERGRpqKV0wFUZq39FFdNwbHs8wNwh/tRl3P+G1eycVxCI/tB+W4O7MnGNS2FiIiIiEjz4Vc1D01d6O4iAFqZrdy+1DW8X256BluWLHMwKhERERGR+qHkoR71iounqDCIkOAiDuWtIzc9g+z5C4iM6e90aCIiIiIix03JQz1qn5xEcXFnAEYeDCV7/gJi58wmKjHB4chERERERI6fkod6lJKZylf5rhlMO7c1fNm7lJlZC0nJTHU4MhERERGR4+d3HaabsmnxE9mS15EdB96kdYdCRm0N4bfTryUqXjUPIiIiItL0qeahHuWmZ7D1yWfIz29DYKAlbPpEsucvIDc9w+nQRERERESOm5KHepS3PofYObNpFdbP9bp8H7FzZpO3vvJk2CIiIiIiTY+aLdWjnlOnANAlsIzcrWmYsp+JSvytOkyLiNSTPn36EB4eTnBwMAA7duxg586dxMXFVVjXqlUrOnbsyHfffUe7du3o3bs3AIcOHaJVq1ZcccUVzJkzh1at9DYoInIs9F+zAUQPHMKeTYGEhR5m9/ZtdOrW3emQRESajXfeeYc+ffoAMHfuXO66664j1nmWxhgmT57M4sWLK+w/adIkCgoKmDdvXuMGLyLSxKnZUgNoFRREYUlXAH7+YY3D0YiINB+jR48mLCysxjIxMTHExMRUu/2cc84hPj6el156qb7DExFp9lTz0EAio2KgaCt5uT86HYqISINIyUxlWvzERj3niy++eNQyl1xyyVHLlJaWsm/fvvoISUSkRVHNQwPpPXAoACFBOykrKXE4GhGR+vfG2v84HcIxs9by4osvsm7dOk477TSnwxERaXKUPDSQohWrKCgIJ6hVKY+//gLgGsp1y5JlzgYmItLCvPXWWyQnJxMXF0dYWBgzZ87k3HPP5dlnn3U6NBGRJkfJQwOJjOnP4W0GgNZlO8hNzyB7/gIiY/o7HJmISN2lZKYy7bVrmPbaNQDe5ymZqQ5HVr3JkyeTlpZGVlYWr7/+OoMGDeKee+6hR48eTocmItLkqM9DA4lKTKDT1hwsX9IntIzs+QuInTNbw7aKSJM2LX6it5/DtNeuIWX6Uw5HdGwmTZpESkoK559/PtnZ2d7hXUVEpHZU89BAUjJTWbD/Q8rKDW3bFLCybyAzsxb69bdzIiItwR//+Ec2btzIyy+/7HQoIiJNjpKHBjItfiJPxM8ib28YxsCgwE48G3dto49MIiLSUC4ccq7TIdRJQkICY8eO5cEHH6S8vNzpcEREmhQlDw3E08chOLwPACFDupI9fwG56RnOBiYiUk+a8pch119/PT/88ANLlixxOhQRkSZFyUMDyVufQ+yc2UQnnQpAcOhuBtx0A3nrcxyOTESkeRg5ciRPP/004Jr47d577/Vue/7550lOTgZ+GW1p+fLl3u2TJk2iX79+zJo1i+TkZHJzcxszdBGRJstYa52OoVkYPny4Xb169RHry8vLWfnO7YQEF9NlwEx69h3gQHQi0pJlZWURFxfndBgiIlJLtfi/bRorlspU89DAAgICKC13DQe4NSfd4WhEREREROpOyUMjaNspFoDCvA0ORyIiIiIiUndKHhpB30FDsRbCQnZTVJDvdDgiIiIiInWi5KERtGnXnvzCtgQGWnKyvnM6HBERERGROlHy0Ai2LFkGpZ0ByPz+S8A1lOuWJcucC0pERERE5BgpeWgEkTH9Kfn6BwA6heV654CIjOnvcGQiIiIiIrWn5KERRCUmkPTrSygtDaRNRBHfPfU0sXNmE5WY4HRoIiIiIiK1puShEaRkpnLtD8+wd18IAD8PimJm1kJSMlMdjkxEREREpPZaOR1ASzAtfiITyqP56j+vQWfo2jaIZ+OuJSpeNQ8iIiIi0nSo5qERePo4DDh1AgCtOxWS9dAj5KZnOByZiIiIiEjtKXloBHnrc4idM5t+p51OXn4oQUFltL70V+Stz3E6NBERERGRWlOzpUbQc+oU7/OAoD7A9xws30f81N84FZKIiIiIyDFTzUMj69A1DoDSwo3OBiIi0gT16dOHwYMHk5ycTHJyMl27dsUYc8S6nj17kpycjDGG9u3be7f179+f2NhY7rvvPkpLS73H/ec//0lycjLBwcEYY/jqq6+qjWHdunUEBAQQGRlJcnIyn3zySWNcuoiIX1Dy0Mj6xSVTXm4ID83lYG6u0+GIiDQ577zzDmlpaaSlpXH11VdXuW7GjBmkpaUBMHnyZO+2nJwcHnnkEW677TbmzZvnPeZvfvMb0tLS6N69O8YY7rvvvmrP/8ADDwAwfPhw0tLSGDt2bMNdrIiIn1Hy0MhCw8PJL+yIMbAx61unwxERaVJGjx5NWFhYjWViYmKIiYmpdvs555xDfHw8L730UpXbJ0+ezFtvvcXatWuP2LZp0yaysrLo3bv3sQUuItJMKHlwQEhkPwBy92Q7HImIyLHZsmTZESPF5aZnsGXJskY5/4svvkiXLl1qLHPJJZdwySWX1FimtLSUffv2VbntpptuolWrVt4aBl8PPfQQN910U+0DFhFpZpQ8OCB0dxEArcxWbl+6CGjcN18RkbqKjOlP9vwF3gTCMxR1ZEx/hyOrHWstL774IuvWreO0006rskyvXr24+OKL+de//sXGjRu963fv3s3//vc/LrzwwkaKVkTE/yh5cEDvuHiKCoMICS7iUN66JvfmKyItV1RiArFzZpM9fwGb/vkvsucvIHbObKIS/XfSy7feeovk5GTi4uIICwtj5syZnHvuuTz77LPV7vOnP/2JsrIyHnroIe+6Rx99lN///vcEBgY2RtgiIn5JyYMD2icnUVzcGYCRB0OaxJuviIhHVGICXc+ewJaUN+h69gS//9/l6TCdlZXF66+/zqBBg7jnnnvo0aNHtfvExcUxefJk/v73v7Nr1y4OHjzIkiVLuOKKKxovcBERP6TkwQEpmal8le9qa9u5TQBf9i5lZtZCUjJTHY5MROToctMz2PHef+k57UJ2vPffI/pA+LNJkyaRmJjI+eefT3FxcY1l//KXv1BYWMijjz7KU089xRVXXEFISEgjRSoi4p+UPDhgWvxELh54DtZC646FjNoawrNx1zItfqLToYmI1MjTzDJ2zmyif3ORtwlTU0og/vjHP7Jx40ZefvnlGsuNHDmS0aNHs3DhQl544QWuueaaRopQRMR/KXlwQG56BtuefIb8/NYEBlrCpk9qcm++ItIy5a3PqdDM0tMHIm99jsOR1V5CQgJjx47lwQcfpLy8vMayf/7znzlw4ABTp06lTZs2jRShiIj/UvLgAM+bb2BYXwAO271N7s1XRFqmnlOnHNHHISoxgZ5TpzgTUB1df/31/PDDDyxZsqTGcmeffTYfffQRf/7znxspMhER/6bkwQGeN98uPeNdK0o3N8k3XxERJ40cOZKnn34acE38du+993q3Pf/88yQnJwO/jLa0fPly7/ZJkybRr18/Zs2aRXJyMq+88grJycls27btiGONGzeOtm3bAvD+++97y61evZrk5GRWr17dCFcrIuIfjLXW6RiaheHDh9tjfQMpLSnhmw9up1WrMqITb6Rj1+4NFJ2ItGRZWVnExcU5HYaIiNRSLf5vm8aKpTLVPDioVVAQhSVdAdj0Q5qzwYiIiIiIHIWSB4dFRMUAkLf/R4cjERERERGpmZIHh/UeMBSAkKAdlJWUOByNiIiIiEj1lDw4rHP3HhQUhhPUqpRNP2Y5HY6IiIiISLWUPDhsy5JllBZ3AmDlmg8B1zwQW5YsczAqEREREZEjKXlwWGRMf4rTNgPQMSTXO3trZEx/hyMTEREREalIyYPDohITSJzyK8rLDe1bF5Dx2BMVZm8VEREREfEXSh4clpKZynU5f2d/bhjGwPrBUczMWkhKZqrToYmIiIiIVNDK6QBaumnxE5lQHs3KD16H9tCjdRjPxl1CVLxqHkRERETEv6jmwWGePg59hp4GQGTXUr6fv4Dc9AyHIxMRERERqUjJg8Py1ucQO2c2A0ePo7AoiJCQItr/9lLy1uc4HZqIiIiISAVKHhzWc+oUohITCAgMpMz2AGB/eS49p05xNjAREQFg27ZtJCcnExkZyZgxY5wOp1ZmzJhB7969McawcePGBj/f008/zeDBgzHGsHjx4gY/n4g4R8mDH2nTMRaAwjzVOoiI+Ivu3buTlpbG8OHD67T/smXLePTRR49Yv2bNGtq3b8/XX399nBEe6fnnn2fevHn1ftzqXH311bzzzjvHdQwn7pOIHDslD36k76ChAIQF76aoIN/haEREpD5U96E4IiKC6OhoIiIiGj8oP6T7JNI0aLQlP9K2fQcOF7QlIuwAG77PIG7oSKdDEhGRBjJw4EDWrFnjdBh+T/dJxL+o5sHPBIZEA7B3+1qHIxER8U+ffPIJkyZN4oQTTiApKYmRI0ce0WTmnHPOoWvXrhhjWLVqFWeddRZ9+vRh1KhRrF1b8f/rm2++yfjx4xk+fDiJiYmMGTOGlStX1hjDP//5T3r06IExhvj4eJYuXQrAk08+Sb9+/Wjfvj3z5s3jrLPO4q233vL2m0hOTub+++/n/fffJzk5GWMMc+fOrXDsdevWMWnSJKKjo0lKSuLEE0/k/vvvJy8vD4ANGzbw29/+luTkZIYOHUpycjIPP/wwZWVldbqfW7duZfr06SQmJjJ06FBGjRrFAw88UKHM5s2bueiii4iOjqZ///6cfPLJfPTRRzUe9/PPPz/iGnNzc0lOTiY4OJgrrrjCW7Yu9ykzM5OJEyfSp08f+vbty4QJE/j222+92337YTz55JPMnDmTpKQk+vTpwxNPPFGneyUigLVWj3p4DBs2zB6vzW8utenv/8eufv9m+8mS26211u7/Lt1ufnPpcR9bRFqudevWOR1CvZo1a5a95ZZbbHl5ubXW2i+//NKGhYXZVatWVSh35513WsBed911tqyszJaUlNjTTz/dnnTSSRXKnXXWWfaZZ57xvn7jjTdsRESE/fnnnyuUGz16tB09erT39eeff24Bu3Tp0grlbrnlFvvEE094X19++eU2Ojq6ymsB7J133ul9vX79ehsVFWVnz57tvb4lS5ZYY4xds2aNtdbaf/3rX3bs2LG2oKDAWmvt9u3b7YABA+zDDz9c4diLFi2ygP3pp5+qPLfH+PHj7VVXXeU9X2pqqnV9PHDZs2eP7dWrl50+fbotKSmx1lr7wgsv2MDAQPv+++97y/30008WsIsWLarxGq21Njo62l5++eUV1h3Lffrxxx9tmzZt7E033eSN+84777QRERE2PT39iJgSEhLshg0brLXWPvPMM9YYY7Oysmq8LyJOqsX/bcc+86rZkh+JjOnP5oceIeDXnWkdUcDGL75g5zN/J3bObKdDE5FmZt5NbzsdAgB3PDzpmPe55ZZb6NSpE8YYAE4++WQSExP5+9//XmWn5iuvvJKAgAACAgKYNGkSc+bMoaioiJCQEAAef/xx+vXr5y1/wQUX8Pvf/55XXnmFP/3pT9XGccopp9CvXz9efvllpkyZAri+kEtJSTlqzUV15s6dS1lZGXfffbf3+s4//3xOPfVUAgJcjQXOOussxo0bR2hoKABdu3Zl6tSpPPfcc8yefezvFytXrmTcuHHe85177rnccsst3u2PPPIIW7ZsYf78+bRq5frYcOWVV/L4449z8803k56eXqdrPR6eWgjf+3Trrbfy1FNPceutt/LWW29VKD9u3Dj69u0LwNSpU5k1axafffYZgwYNatS4RZoDJQ9+JCoxgbibb2T1mjdp1+0w2Z+9x8g5s4lK1GzTIiIeERER3HbbbXz66aeUlJQQEBDA+vXradu2bZXlBw4c6H3evn17AHbt2kWvXr0ACA0N5ZprrmHVqlWUl5djjGHfvn1s2LChxjiMMVxyySXcf//97N+/n3bt2vHpp58SHx9Phw4d6nRtH3zwAUOGDCE8PLzC+s8++8z7vHXr1jz11FO8+uqrHDhwgFatWrFjxw72799fp3Oefvrp3HXXXWzevJlLL72UUaNGce+993q3f/jhh3Tt2tV7vzxOPPFEnnnmGXbu3EmXLl3qdO66+vDDDxkyZAhhYWHedUFBQQwdOpQPP/wQa603qYCqfwd27tzZeAGLNCNKHvxISmYqb2T9h8kH2tGuG5T2CmRm1kIuDDiXafETnQ5PRJqRunzj7w/Ky8uZNGkSBw4c4P3336dnz54AjBkzhqKioir38f0g7vn23tM/4PDhw4wdO5aePXvy8ccf065dOwD69OlT7fF8XXrppcybN4/XXnuNq6++mhdffJHLLruszte3Z88ehg0bVmOZ2267jb/97W989NFHnHzyyYDrm/i77rqrTud84403eOihh3j++ed5+umn6d27N7fffjszZszwxuS5L748H8L37NnT6MlDdfepffv2FBQUkJ+fX2F0ppp+B0Tk2KjDtB+ZFj+RZ+Oupdf3BwDo0L6IhQNnKXEQEXFbv349K1as4He/+503cTgeX375JTk5OVx33XVVfkA+mpiYGEaNGsXLL79Mfn4+n376Keeee26d4+nYseNRaxBeeuklzjzzTG/icLzCw8O544472LRpEx999BHR0dFcddVVfPjhh96Y9u3bd8R+nnUdO3as8fgBAQFYayusO3z48HHFXFNMYWFhR9TciEj9UfLgR3LTM8iev4Ckq6/m4OFQgoLKSHv1n+SmZzgdmoiIX/DUBvg2SQHYsWNHvR2vvLyc3bt31/oYl112GcuXL+fBBx/k3HPPJTg4uML2oKAg74fnw4cPH9Ee39eZZ57J2rVrKSgoqLB+2rRpfPrpp96Y6+v6AS666CLAdQ/GjRvHsmXLALx9Gc444wx27tzJzz//XGG/VatWkZCQcNRah86dO1dIiPbu3cvevXuPKHcs9+mMM85g7dq15Of/MidSaWkpaWlpnHHGGUfcHxGpP0oe/Eje+hxi3X0cdhe42u4GnRhD3nrNOC0iAjBo0CD69evHokWLvB9IX3/9dbKzs+t0vJNPPpmoqCgWLlxIYWEhAA8//HCFD6VHM336dIKDg7nnnnuqbLLUt29f9uzZQ1FREcuXL+eGG26o9lhz584lICCAuXPnej9I/+Mf/2DNmjWMHOma++fcc8/lgw8+ICPD9cXSDz/8wGuvvVbreCt79dVXWbJkiff1F198QWBgIKeffjoAN954Iz179mTOnDmUlpYCrtqP7777joceeuioxx89ejQffPCBd6jZRx55hMjIyCPKHct9uvPOOzHGcPvtt3vv01//+lcOHTpUob+GiDQAJ4d6ak6P+hiq1VfWmq/t6vdvtp8tm1evxxWRlqe5DdWamZlpx44da7t06WJHjx5tb7jhBjts2DAbERFhk5KSbFFRkb344ottly5dLGCTkpJsenq6ve+++2yvXr0sYOPi4uzrr79urXUNuTpixAjbvXt3O2bMGHvXXXfZHj162Hbt2tnx48fbrVu32qSkJBsREeE9x86dOyvENGXKFDtw4MAq4925c6cdM2aMHTBggB0yZIj997//bd977z2blJRkAdulSxc7adIkb/m1a9fac8891/bu3dsmJSXZ8847z+bk5Hi379u3z1522WW2S5cudtSoUXbatGn2sssu817rl19+aX/3u99VuNbnnnuu2vv5wAMP2BEjRtjExESbmJhoTzzxRLts2bIKZTZt2mSnT59ue/XqZfv162dPOukk+8EHH3i3P/XUUzYuLs4CtlevXnbWrFnebT///LMdP3687dGjhx0zZox99913bXR0tG3Xrp31fe881vuUkZFhzznnHNu7d28bHR1tzzjjDLt69Wrv9ldffbVCTPfcc49du3ZtheNdeuml1d4XESf581CtxlZqhyh1M3z4cLt69ep6O15xYSFpn9xJYGA5/YbOoV2nzvV2bBFpWbKysoiLi3M6jGZt3rx5BAYGcuuttzodiog0A7X4v+1Y2zw1W/JTwaGhFBS7Eoafvv/2KKVFRMRJS5cu5ZJLLnE6DBGRBqfkwY+FtxkAwKF9PzgciYiIVDZ69GiKior4/PPP6datG9HR0U6HJCLS4JQ8+LGwHa6h7EJabeeOpYsA14hMW5YsczAqEREB1+hEgwYN4qabbuKRRx5xOhwRkUahSeL8WM8h8WSsX0t4ZBEle7O9Q7nGzpntdGgiIi2eZ+hUEZGWRDUPfiwqMYEy2w2AofnB3sQhKjHB4chEREREpCVS8uDHUjJT+fTwNgA6t7d82buUmVkLSclMdTgyEREREWmJ1GzJj02Ln8iuom5s3JVC66hChu+O5Ldx1xIVr5oHEREREWl8qnnwY7npGfz06BPk57cDIGTSGLLnLyA3PcPhyERERESkJVLy4Mfy1ucQO2c24e0GAVBkdxM7ZzZ563McjkxEREREWiI1W/JjPadOAaD35jZszfqS4MCttB4yWB2mRURERMQRqnloAjr36EVhYRhBQaX8/GOW0+GIiIiISAvVbJMHY8wSY4w1xnzqdCzHKyAggPKA3gBs35TucDQiIi3Ltm3bSE5OJjIykjFjxjgdTq3MmDGD3r17Y4xh48aNDX6+p59+msGDB2OMYfHixQ1+vuburbfeIjk5GWMMc+fOdTockQr8LnkwxvQxxlxnjHnPGLPdGFNijDlkjEkzxtxtjGlXi2NMAc4/hnNOMsZ8YozJdZ9rpTHm8uO5jvrWodtgAEoLNzobiIhIC9O9e3fS0tIYPnx4nfZftmwZjz766BHr16xZQ/v27fn666+PM8IjPf/888ybN6/ej1udq6++mnfeeee4juHEffIHixcvPiLhmjx5MmlpaY7EI3WXmppKp06d2Lx5s9OhNCi/Sx6AlcCjwFpgMhADTAAygNuAb4wxHavb2RjTBngC+Lk2JzPG3A68BewDxgAnAmnAYmPMc3W8hnrXLy6Z8nJDeOh+Du7f53Q4IiJSS9V9KI6IiCA6OpqIiIjGD8oPtdT7VFXyIE1TmzZtiI6OJiQkxOlQGpQ/Jg8AC621N1lrV1lrN1lrV1hrLwX+B/QFrq5h3/uBMuDBo53EGDMamAesAaZZa9OstVnW2quBt4EZxpjLjvtq6sGe9/7L4cPtMAb+9c4rgGso1y1LljkbmIiI1MnAgQNZs2YNQ4YMcToUv6b7JE3F6aefzurVq+ncubPToTQof0weZgEPVLPtG/eyyp+KMeZkXInFtcDhWpzrTvfyMWttWaVtC9zLO2pxnAYXGdOfog35ALQ3e8lNzyB7/gIiY/o7HJmISOP65JNPmDRpEieccAJJSUmMHDnyiCYz55xzDl27dsUYw6pVqzjrrLPo06cPo0aNYu3atRXKvvnmm4wfP57hw4eTmJjImDFjWLlyZY0x/POf/6RHjx4YY4iPj2fp0qUAPPnkk/Tr14/27dszb948zjrrLN566y1vv4nk5GTuv/9+3n///WrbtK9bt45JkyYRHR1NUlISJ554Ivfffz95eXkAbNiwgd/+9rckJyczdOhQkpOTefjhhykrq/w2Vjtbt25l+vTpJCYmMnToUEaNGsUDD1R8G968eTMXXXQR0dHR9O/fn5NPPpmPPvqoxuN+/vnnR1xjbm4uycnJBAcHc8UVV3jL1uU+ZWZmMnHiRPr06UPfvn2ZMGEC3377rXe7bz+MJ598kpkzZ5KUlESfPn144oknjnpfkpOTad++PX369OHdd99l3Lhx9OzZkzPPPJMtW7ZUKFteXs4DDzzAwIEDiY2NJSYmhrlz51JaWlrt8cvKykhOTmb16tWsXr3ae90vvfTSEeX+9Kc/MWzYMHr27Mmtt956xLF+/vlnpk2bRnR0NDExMYwdO7ZCMy/fe7Fw4UJmzZrFsGHDCAwM5IYbbqjQR+bDDz9k0qRJ9O3blxNPPJGMjAy2b9/OtGnT6N+/PyNHjmTdunVHxPD3v/+d+Ph4YmNj6dOnDzfccAOHD9fmo5irj8eIESMYMGAA0dHRXHnllezatatCmcWLFzN06FCGDh1KYmIil156qbdpV+Xru+qqq4iPj6dfv368/vrrlJaWMnv2bBITE+nfvz9Lliw5IoYvvviC0aNH069fP6Kjo7ngggvIyak4NP5//vMfRo4cyQknnEBiYiJTp07l008/BWDRokVV9vupzd9Xk2OtbRIPwACfAxaYUsX2YFxNnVLcr69wl/20muN1xlVDYYFeVWwPAgrd24cdLb5hw4bZhpb9yUd29fs32y/f+rNdcckVdv936Q1+ThFp+tatW+d0CPVq1qxZ9pZbbrHl5eXWWmu//PJLGxYWZletWlWh3J133mkBe91119mysjJbUlJiTz/9dHvSSSdVKHfWWWfZZ555xvv6jTfesBEREfbnn3+uUG706NF29OjR3teff/65BezSpUsrlLvlllvsE0884X19+eWX2+jo6CqvBbB33nmn9/X69ettVFSUnT17tvf6lixZYo0xds2aNdZaa//1r3/ZsWPH2oKCAmuttdu3b7cDBgywDz/8cIVjL1q0yAL2p59+qvLcHuPHj7dXXXWV93ypqanW9fHAZc+ePbZXr152+vTptqSkxFpr7QsvvGADAwPt+++/7y33008/WcAuWrSoxmu01tro6Gh7+eWXV1h3LPfpxx9/tG3atLE33XSTN+4777zTRkRE2PT0X94bPTElJCTYDRs2WGutfeaZZ6wxxmZlZdV4XzwxtWnTxt5+++3WWmsPHTpkBw4caC+66KIK5a655hrbtWtXm52d7T1v79697aWXXnrUc1T+vap83dHR0Xb16tXWWmvff/99C1S4756fz3nnnWeLi4uttdbOnz/fhoeHV7hGz72IjY21GRkZ1lprH374YXv99ddba3/5fZk+fbotLCy0JSUl9rTTTrMJCQn2jjvu8K475ZRT7CmnnFIhzgceeMBGRETY5cuXW2ut3bt3rx06dKgdM2aMLSsrq/H6X331VRsYGGiXLVtmrbW2oKDAnnXWWTYuLs7m5eVZa6397LPPbEhIiM3JybHWWpuXl2dHjx5d4XfCc32JiYl206ZN1lpr//KXv9igoCB7++23e9f98Y9/tBEREXb//v3efT///HMbHBxsH3vsMWuttWVlZfbKK6+0Xbp0sdu2bbPWuv42g4OD7RdffGGttba4uNhedNFFFX6Pq/obONrfV3Vq8X/bsc/kfj/PgzEmBBgEzAFOAu611i6rouifgR7A+FoeehiumpfD1tojerZYa0uMMRuAOGAEv9R6OCIlM5U3dvyH68LaExZWzKpBbViQtZALA85lWvxEJ0MTkSbom//OcToEAIZNmH/M+9xyyy106tQJYwwAJ598MomJifz973+vslPzlVdeSUBAAAEBAUyaNIk5c+ZQVFTkbZf8+OOP069fP2/5Cy64gN///ve88sor/OlPf6o2jlNOOYV+/frx8ssvM2XKFMD1hVxKSspRay6qM3fuXMrKyrj77ru913f++edz6qmnEhDgaixw1llnMW7cOEJDQwHo2rUrU6dO5bnnnmP27NnHfM6VK1cybtw47/nOPfdcbrnlFu/2Rx55hC1btjB//nxatXJ9bLjyyit5/PHHufnmm0lPb/xRAD21EL736dZbb+Wpp57i1ltv5a233qpQfty4cfTt2xeAqVOnMmvWLD777DMGDRp01HMdOnSIG264AYDIyEjOPPPMCt9c//jjjzz99NPcfffdDBw4EIA+ffpw0003cf3113PDDTdwwgkn1Plak5OTGTZsGAATJkwgMjKSTz/9lAkTJgCun8/mzZv5+OOPCQoKAuDGG2/k0Ucf5f777z+iP8X48eOJj48HYNasWRw6dKjC9osvvtj7t3Heeedx8803c88993jXTZkyhT/+8Y8UFxcTHBzMgQMHuOuuu7j44os56aSTAGjfvj1z587lvPPOY9myZUydOrXKa7PWMmfOHMaOHct5550HQGhoKA8++CBJSUk888wzzJ49m6+++oqQkBC6d+8OuPrB3H333VXWbIwfP57evV0jVF5wwQXcd9995OXledf96le/4sEHH2TVqlWceeaZAPz5z3+md+/e/OEPfwBco1w++OCD/POf/+S+++7jscceY82aNRQXF3t/j4KCgrj11lv57rvvavz5He3vqynyx2ZLXsaYFbi+/U/DlUCcYq29rYpyg4BbgD9Za3fU8vCe9j47ayiz3b3sV0OZRjEtfiLPDr6W/B2ufwx9wtvybNy1ShxEpMWJiIjgtttuY9iwYSQmJpKcnExmZiYbNmyosrznAx24PtQAFZpEhIaGcs011zB06FCSkpJITk5m37591R7PwxjDJZdcQmpqKvv37wfg008/JT4+ng4dOtTp2j744AOGDBlCeHh4hfWfffYZiYmJALRu3ZrXXnuNU045hfj4eJKTk1m8ePFR463O6aefzl133cU111zD8uXLKS8v59577/Vu//DDD+natSu9evWqsJ+nScvOnTW9jTaMDz/8kCFDhhAWFuZdFxQUxNChQ/nwww89LQi8qvodqG3cHTt29O7j2d93348++ghrLSNGjKiw34knngi4fqbHwzd2gHbt2lU4/4cffkjnzp2JiYnxrgsMDGTQoEHeJjW+4uLivM8jIiLo2rVrhe2+x/Fct++6Dh06YK31xrBixQry8/PrdP3Z2dls3rz5iH0TExMJDQ317nvqqaeSl5fHiSeeyLPPPsvu3bs57bTTOPvss484Zm3iB9ixw/VxMT8/nxUrVhwRQ8eOHenbt683hhEjRhAWFsYpp5zCggUL2Lx5M0OGDOHiiy+u9vrg6H9fTZG/1zxMA1oDscD1wHJjzF+BO6215QDGlco9i6tm4NljOHYb9zK/hjIF7mXbYwm6IXj6OHT49URgJWHRwWTPX0DsnNmacVpEjlldvvH3B+Xl5UyaNIkDBw7w/vvv07NnTwDGjBlDUVFRlfv4fhD3fHvv6R9w+PBhxo4dS8+ePfn4449p1841GnifPn2qPZ6vSy+9lHnz5vHaa69x9dVX8+KLL3LZZXUfZ2PPnj3eb5mrc9ttt/G3v/2Njz76iJNPPhlwfRN/11131emcb7zxBg899BDPP/88Tz/9NL179+b2229nxowZ3pg898WX54PZnj176NKlS53OXVfV3af27dtTUFBAfn5+hdGZavodOJrKiVxAQADl5eUVYgGOuEe+9+d4VHV+39j37NnDwYMHSU5OrlAuNzf3iCQKXLUntT2f59vyqtZ5Yqjt9c+YMYPVq1d7tz///PMUFhZWua9nnWffUaNG8b///Y8HHniA3//+91x77bWcd955PPbYY/To0eO44t+3bx/l5eXV/o7/+OOPAERHR/PVV19x3333ceutt3LTTTcxbtw4Hn/8cQYPHnzEvh5H+/tqivy65sFau9lau85auxQYB6zANVyr70hKM4FRwExb1V/J8TGeUKrcaMxMY8xqY8zq3bt31/OpK8pbn0PsnNnEjZ5AebkhIuIAvf7f1eStzzn6ziIizcT69etZsWIFv/vd77yJw/H48ssvycnJ4brrrqvyw8PRxMTEMGrUKF5++WXy8/P59NNPOffcc+scT8eOHb21GNV56aWXOPPMM72Jw/EKDw/njjvuYNOmTXz00UdER0dz1VVX8eGHH3pj2rfvyCHCPes6dqx29HTA9WG38ttzbTvSVqemmMLCwo74wN2QPNdfOZ7a3p/6OH+3bt1IS0ur8Ni4cSObNm1q0HN7zg9Hv/7nn3++QnzDhw+vdl+A/fv3V7h3p556Km+//TZbtmxh3rx5vPvuu0yfPv2442/fvj0BAQHV/j75xpCQkMArr7zCjh07ePLJJ0lLS+Pss8+ukExWdrS/r6bIr5MHX+6aBs/XKn8wxkQZY7rhGpnpAWvt2ur3rtJB97Km/zChlcpWjulZa+1wa+3wTp06HePpj03PqVOISkwgPLI1+YUdMAb22Hx6Tp3SoOcVEfEnntoAz7eHHp4mCPVxvPLyco7lC6HLLruM5cuX8+CDD3LuuecSHBxcYXtQUJD3w/Phw4ePaI/v68wzz2Tt2rUUFBRUWD9t2jRvE5SioqJ6u36Aiy66CHDdg3HjxrFs2TIAb1+GM844g507d/LzzxWnT1q1ahUJCQlHrXXo3LlzhYRo79697N2794hyx3KfzjjjDNauXUt+/i+NB0pLS0lLS+OMM8444v40pPHjx3tH9fLlee1pV18d3+vevXv3MX+o9Iz+VPnD7/vvv8/tt99+TMeqi5NOOonw8PA6XX9sbCy9evU6Yt+MjAwKCwu9+77yyiu8/fbbAHTp0oVbbrmFGTNm1Et/m/DwcE466aQjYti7dy8//fSTN4aPPvqI559/HoC2bdty7bXXcuutt7J582Zyc3OrPf7R/r6aoiaTPLhluJfBwFBck8e1BW4yxuT5PoCn3WVP81nvm2B4vrKv6b9eN/eybg1JG0hIpKu7Ru7u7x2ORESkcQ0aNIh+/fqxaNEi7wfS119/nezs7Dod7+STTyYqKoqFCxd6m1A8/PDDFT6UHs306dMJDg7mnnvuqbLJUt++fdmzZw9FRUUsX77c2/m2KnPnziUgIIC5c+d6P1D+4x//YM2aNYwcORJwdbj84IMPyMhwvSX+8MMPvPbaa7WOt7JXX321QgfgL774gsDAQE4//XTA1fm2Z8+ezJkzxzv06EsvvcR3333HQw89dNTjjx49mg8++MA71OwjjzxSZdOZY7lPd955J8YYbr/9du99+utf/8qhQ4cavT35gAEDuPrqq3nyySf54YcfANfQqQsWLODSSy89amfpvn37snXrVqy1LFu2jL/+9a/HdP4bb7yR7t27c8MNN1BcXAy4hta9/vrrSUpKqttFHYO2bdty55138uqrr7JixQrAVWswd+5cxowZ4x1MoCrGGObPn88nn3zi/VBdWFjIn/70JwYNGsSsWbMA1+/4Aw88wMGDru9yCwoK+Pbbbxk3bly9XMP999/Pzz//zGOPPQa4vkD485//TLt27fjLX/4CuO7pAw884O3rUVpayldffUViYmKFPjGVHe3vq0lycqinyg9cfRsuq2F7G1xNiCxwFq7+EDHVPP7oLveVz7pon2PVZqjWAvf24UeLvTGGavX4OecH15Ctb99y1CHQRESa21CtmZmZduzYsbZLly529OjR9oYbbrDDhg2zERERNikpyRYVFdmLL77YdunSxQI2KSnJpqen2/vuu8/26tXLAjYuLs6+/vrr1lrXMI0jRoyw3bt3t2PGjLF33XWX7dGjh23Xrp0dP3683bp1q01KSrIRERHec+zcubNCTFOmTLEDBw6sMt6dO3faMWPG2AEDBtghQ4bYf//73/a9996zSUlJFrBdunSxkyZN8pZfu3atPffcc23v3r1tUlKSPe+887xDVFpr7b59++xll11mu3TpYkeNGmWnTZtmL7vsMu+1fvnll/Z3v/tdhWt97rnnqr2fDzzwgB0xYoRNTEy0iYmJ9sQTT/QOm+mxadMmO336dNurVy/br18/e9JJJ9kPPvjAu/2pp56ycXFxFrC9evWys2bN8m77+eef7fjx422PHj3smDFj7Lvvvmujo6Ntu3btrO9757Hep4yMDHvOOefY3r172+joaHvGGWd4hzS11jUEqG9M99xzj127dm2F49U0lOqYMWNsu3btbFBQkE1KSrJ79uyxf/jDHyr8Xn355ZfWWtfQnvfdd5+NiYmxAwcOtP369bN33nmnd2jbmmRnZ9thw4bZQYMG2eTkZLtixQr72WefVYjz2muvtfv377dJSUk2KCjI+7vpsXnzZnvRRRfZnj172uTkZDtq1Cj72muvVXsvkpKSbGlpqXf7zTffXOH35b333rP33HPPUde9+uqr3mM899xzdsiQIXbgwIE2OjraXnfddfbQoUNHvX5rrV22bJkdNmyYjYmJsb169bKXX355hb+xtLQ0+5vf/MbGxcXZpKQkGxcXZ2fOnGn37t1b5fXdc8899r333jvquptvvtl7js8++8yefvrptm/fvrZ37972/PPPt+vXr/du37Bhg501a5YdPHiwN4Zf//rX3iFgX3jhhSr/Bmrz91UVfx6q1dh67yZQd8aYK4C/AxHW2sIqto/C1e8BIMZaW22Df/exFgH/s9aOqabMp8Bo4Epr7eJK28YAnwA/WWuPOtrS8OHDrW9HoIZUXl7OynfuICS4iE79ZtA7JrZRzisiTVNWVlaFEVak/s2bN4/AwMAqJ/ASETlWtfi/3Xht8yrxx2ZLAUB1XdA9jfdW1JQ4HANPH4rrjDGBlbbd6F7Oq4fz1KuAgABKrauj4Lafmm6bORGR5mLp0qVccsklTochItLg/C158Mzj/rAxZr4xZpQxpq8x5kxjzH+Bc4AtwOXVHcAY08kY05VfhlcNNsZ0dT/CfMtaaz/BlUAMBVKMMUnGmDhjzFPAZGBx5RoJf7BlyTLCjKv3/8EDWYBrKNctS5Y5GJWISMsyevRoioqK+Pzzz+nWrRvR0dFOhyQi0uD8ap4Ha+0/jDFbgOm4OkPPwjUaUh7wPa6ah8ettQdqOMwqwPc/+En8MtnblcDiSueca4xZg6um4TMgEFgL/NZau+h4r6khRMb0p+hvjxPyq7a0izzE9q9X8fPjC4mdc+wzi4qISN0YYxg0aBCdOnXi5ZdfdjocEZFG4Vd9HpqyxuzzAK6ahm/Xvk7bDgXs/bIVJ15wsSaLE5Eqqc+DiEjToj4PUq9SMlOZmbWQXftdiV9Bn2BmZi0kJTPV4chEREREpDnzq2ZLUjvT4icyoTyabxf9HWKC6NShhKd7X037+IYfz1lEREREWi7VPDRBuekZZM9fQPLlV1JQGERIaAlpLy0mNz3j6DuLSIukJqoiIk2Dv/+/VvLQBOWtzyF2zmzaJyex63AUAEGnxpO3vj5GrxWR5iYoKIiCggKnwxARkVooKCggKCjI6TCqpeShCeo5dYq3c3RszMkA2MCd9Jw6xcGoRMRfde7cma1bt5Kfn+/332iJiLRU1lry8/PZunUrnTt3djqcaqnPQxPXb8hQMj97i/DQ/Rzcv4827do7HZKI+Jk2bdoAsG3bNkpKShyORkREqhMUFESXLl28/7f9kZKHJi4sPIL8ok5Ehu1iQ9a3JJ98htMhiYgfatOmjV+/GYmISNOgZkvNQFjrGAAO7vne4UhEREREpDlT8tAM9B5wAgDBgVspU5MEEREREWkgSh6agdJVaygoCCMoqJTHX38BcA3numXJMmcDExEREZFmRclDMxAZ05/DW10/yrblO7zzQETG9Hc4MhERERFpTpQ8NANRiQl07utqutQntJTs+QuInTPbO5yriIiIiEh9UPLQDKRkpjI/931KSwNo06aQr2OCmZm1kJTMVKdDExEREZFmREO1NgPT4icyoTya1WvepF23wwwM6sSlcdOJilfNg4iIiIjUH9U8NAOePg4RUQMACBnYluz5C8hNz3A4MhERERFpTpQ8NAN563OInTObASPHAxAWvpd+N15H3vochyMTERERkeZEzZaagZ5Tp3ifZ33VmvCwQ+wNKGWQz3oRERERkeOlmodmJiC4LwB7tqnJkoiIiIjULyUPzUzX3u5O0qWbnA1ERERERJodJQ/NTJ+BQygpaUVYaD47NiuBEBEREZH6o+Shmdn+9n8oLOgEwHufLQNcozFtWbLMuaBEREREpFlQ8tDMRMb0pyhrDwAdg/Z7h3GNjOnvcGQiIiIi0tRptKVmJioxgcH5h9l54E06t8ln7aN/Y8ic2UQlasI4ERERETk+qnloZlIyU7lp04scOBBGYIAle3B7ZmYtJCUz1enQRERERKSJU81DMzMtfiITyqP56qM3IAp6tg7n2bhLiYpXzYOIiIiIHB/VPDQznj4OvZNOASCieynfz19AbrrmfRARERGR46PkoZnJW59D7JzZxI4eT2FREKEhRbT77SXkrc9xOjQRERERaeLUbKmZ6Tl1ivd5Gb2ADewv3Ufs1F85FpOIiIiINA+qeWjG2nUZDEBxvmodREREROT4KXloxvrFnUB5uSEsZC95Bw84HY6IiIiINHFKHpqxiNatyS/qQEAAbFj7jdPhiIiIiEgTp+ShmQuNiAEgd3eWw5GIiIiISFOn5KEZ27JkGe1btQOgVcAWysvKyE3PYMuSZc4GJiIiIiJNkpKHZiwypj97//4yhYWhhASXkvXJh2TPX0BkTH+nQxMRERGRJkjJQzMWlZjAoDmzydsaCMCWdcuJnTObqETNNi0iIiIix07JQzOWkpnKzKyFbD50GIDg7uXMzFpISmaqw5GJiIiISFOkSeKasWnxE5lQHs26dx6lNL4jbdoU8mDEb+gTf6rToYmIiIhIE6Sah2YsNz2D7PkLGHzTDWw/GAZA9mfvkZue4XBkIiIiItIUKXloxvLW53j7OOwp6wBAyMDW5K3XjNMiIiIicuyUPDRjPadO8XaOvvCsXwMQHr6XjmdPcDIsEREREWmilDy0EFEdO3G4IIrAQEvOum+dDkdEREREmiAlDy1IUJhrfod9O9Y6HImIiIiINEVKHlqQHv2GAhDIz5SXlTkcjYiIiIg0NUoeWpLv1lFYFEJIcDELUp4HXCMybVmyzNm4RERERKRJUPLQgrQZEMPhra6pPaJKd3iHco2M6e9wZCIiIiLSFCh5aEGiEhPoFJ0MQL/QUrLnL/AO5SoiIiIicjRKHlqQlMxUHsr9gLKyANq2LeSr/kHMzFpISmaq06GJiIiISBPQyukApPFMi5/IhPJoVq95k3bdDhMb3JnL4qYTFa+aBxERERE5OtU8tCCePg4RUQMACBnYhuz5C8hNz3A4MhERERFpCpQ8tCB563OInTObQSefCUB4xF763vD/yFuf43BkIiIiItIUqNlSC9Jz6hTv88MFbYkIO8Aeihjss15EREREpDqqeWihWoW6hmfduz3T4UhEREREpKlQ8tBC9eibDECg/Zny8nJngxERERGRJkHJQwvVKyaWouIQQkKK2JLzo9PhiIiIiEgToOShhdq27C2KCzsD8PnX7wCu0Zi2LFnmYFQiIiIi4s+UPLRQkTH9Kf5uKwCdQ/d7h3GNjOnvcGQiIiIi4q+UPLRQUYkJJE2ZRllZAO1bF5Dx5FPEzplNVKImjBMRERGRqil5aKFSMlP5w/rn2bsvFICNg9sxM2shKZmpDkcmIiIiIv5K8zy0UNPiJzKhPJqvUl+DTtC9XSuejbuWqHjVPIiIiIhI1VTz0EJ5+jgMHnsu1kLrTgWsffRv5KZnOB2aiIiIiPgpJQ8tVN76HGLnzKbXSSexPy+SwEBL6AX/R976HKdDExERERE/pWZLLVTPqVO8z1u3GQz2aw6X7aTn1GudC0pERERE/JpqHoTeA4YBEBS4hbKSEoejERERERF/peRB6Nq7DwWF4QQHlfDTD2udDkdERERE/JSSByEgIABa9QVgx6Y0Z4MREREREb+l5EHYsmQZbVt1BqC8OIfy8nJy0zPYsmSZs4GJiIiIiF9R8iBExvQn7x9vUFzSiojwQnI++5Ts+QuIjOnvdGgiIiIi4keUPAhRiQnE3Xwjh7aHAJDzzafEzplNVKImjBMRERGRXyh5EFIyU5mZtZCtBwoAaNULZmYtJCUz1eHIRERERMSfaJ4HYVr8RCaUR7P2/b9RNrg9UW0LeKjNpfSOP9np0ERERETEj6jmQchNzyB7/gKG3HA9Ow+GYwxkffIOuekZTocmIiIiIn5EyYOQtz7H28dhV0kHAEJi25C3PsfhyERERETEnyh5EHpOneLtHD317F8DEB6+l45nT3AyLBERERHxM0oepIJ2HTtxuCCKwMBy1q/9xulwRERERMSPKHmQIwSFDwBg345MhyMREREREX+i5EGO0HvAMACCzGbKysocjkZERERE/IWSBzlC2TfpFBSEEhxcwuOvPQe4RmTasmSZs4GJiIiIiKOUPMgR2gyI4fDWQACiynd6h3KNjOnvcGQiIiIi4iQlD3KEqMQEusWMAKBfRAnfz1/gHcpVRERERFouJQ9yhJTMVO7d9x+Ki1sREVHEtwMjmJm1kJTMVKdDExEREREHtXI6APE/0+InMqE8mq+/epMO0aX0i2zPxXG/IipeNQ8iIiIiLZlqHuQInj4OHXokARDepxXZ8xeQm57hcGQiIiIi4iS/Sx6MMX2MMdcZY94zxmw3xpQYYw4ZY9KMMXcbY9rVxz6V9p9kjPnEGJPr3m+lMebyhrtK/5a3PofYObNJGPd/lJYGEBFxiM6zfkfe+hynQxMRERERB/ld8gCsBB4F1gKTgRhgApAB3AZ8Y4zpWA/7AGCMuR14C9gHjAFOBNKAxcaY5+rtqpqQnlOnEJWYQHBoKIUl3QHYXbyHnlOnOBuYiIiIiDjKX/s8LLTW3uTzehOwwhjTCxgNXA3cc7z7GGNGA/OANcA0a61nRrSrjTHdgRnGmM+ttS/V14U1NW07DaEsbwuFh7KdDkVEREREHOaPNQ+zgAeq2faNe9m5HvYBuNO9fMwncfBY4F7eUc1xW4QBCSMpLzeEh+7hwP59TocjIiIiIg7yu+TBWvtva+3myuuNMQZXkyKAj493H2NMZ1w1EgAfVRHKl0AR0N8YM+yYLqIZiWjdmvyizhgDOZlfOx2OiIiIiDjI75KHyowxIcaYJOBl4CTgXmvtsnrYZxiu6z9cVeJhrS0BNrhfjjiui2jCtixZRjDdANi+3VWJk5uewZYlyxyMSkRERESc4NfJgzFmBVCIqwPzIOAUa+1t9bRPf/dyZw2H2+5e9qt91M1LZEx/St/5AmuhU+sD7Fi9muz5C4iM6X/0nUVERESkWfHr5AGYBgwBpgJ5wHL30Ks1xV3bfdq4l/k1HKvAvWx7zJE3E1GJCST8/loO7Q8jMNCSnvomsXNmE5WoCeNEREREWhq/Th6stZutteustUuBccAKXEOvPlif+9TAeA5b5UZjZhpjVhtjVu/evbsOh/d/KZmpzMxayI69rltQ3DeYmVkLSclMdTgyEREREWlsfp08+LLWlgN3uV/+wRgTdZz7HHQvw2s4RGilspWP/6y1dri1dninTp2OFk6TNC1+Is/GXUu3rFwAOnYsYuGAmUyLn+hsYCIiIiLS6JpM8uCW4V4GA0OPcx/PdMldati3m3u5oYYyzVpuegbZ8xeQPHMmuXmhBAWVsebNV8lNzzj6ziIiIiLSrPhV8mCMiTXGXFZDEd/+CcF13cftG6AciHBPJFc5liCgr/vl6hoDb8by1ud4+zjsLGgHQNDQHuStzznKniIiIiLS3PhV8oBrWNVFxpjQarYP9nm+/jj2wVq7C/jc/XJ8FfudgqvZ0k/W2habPPScOsXbOfq0Ua6mSiGhO+l23iQnwxIRERERB/hb8gCumGZUs+1293KFtdb3q++67AO/9Ie4zhgTWGnbje7lvKPE22L06BtDYWEYwcHFbMxe63Q4IiIiItLI/C15KHUvHzbGzDfGjDLG9DXGnGmM+S9wDrAFuPw49wHAWvsJrgRiKJBijEkyxsQZY54CJgOLrbWLG+RKm6CAgABsK9eUF9s3futwNCIiIiLS2PwqebDW/gMYC7wATAD+C/wIvI5rXobbgXhr7Y/Hs0+lc84FpgAdgM+AVcAJwG+ttVfW9zU2ZVuWLKNtq84AlJfkUF5ertmmRURERFoQY22VUxjIMRo+fLhdvbp5d43ITc/g+4cWYKd1Izi4hMhWZ5O76B+aNE5ERESkcZmjF2kYflXzIP4tKjGBQTfP5tA216BVGzM+U+IgIiIi0oIoeZBa88w2/XNuHgAhvco127SIiIhIC9LK6QCk6ZgWP5EJ5dFk/ecRiod0pnVkIfeGTGNA/FinQxMRERGRRqCaB6k1z2zTcTffyOYDYQDkfPOJZpsWERERaSGUPEit+c42vZtOAIT3DdRs0yIiIiIthJIHqTXf2aavuvB3lJYGEhGRR9CJJzgcmYiIiIg0BiUPUifBwSEUlfYE4Kfvv3Y4GhERERFpDEoepM7ad0sCoPjw9w5HIiIiIiKNQcmD1NmAhBGupkthB9i9favT4YiIiIhIA1PyIHW2+533yM93dZx+9+M3ANeITFuWLHMwKhERERFpKEoepM4iY/pTkrkHgC4h+7xDuUbG9Hc4MhERERFpCEoepM6iEhNImnghZWUBdGyTT/rCp71DuYqIiIhI86PkQeosJTOVP+Q8z+49oQBsHtKOmVkLSclMdTgyEREREWkIrZwOQJquafETmVAezddL/wVdoFvHQJ6Nu5aoeNU8iIiIiDRHqnmQOvP0cUg453zKygxt2ueT8eRT5KZnOB2aiIiIiDQAJQ9SZ3nrc4idM5tuI0aw+2AUxkDg/51E3vocp0MTERERkQbQYMmDMaaDMWZdQx1fnNdz6hRv5+gePYcDUGK20XPqFAejEhEREZGGctTkwRhTZozpXIdjtwJi67CfNEEDEkdRXm6ICN3Dgf37nA5HRERERBpAbWoeDBBYh2PfXod9pImKaN2G/KIuGAPrM1Y6HY6IiIiINIDaNlt6/FgOaox5HLj22MORpmrLkmWE0A2AXTu+BTTbtIiIiEhzU9vkYZwx5tHaFHQnDr93v7yuLkFJ0xMZ05+Stz+jvBw6tjnA1pVfabZpERERkWamtsnDr4ErjDE31lSoUuLw/6y1Tx5PcNJ0RCUmkHDd7zm4N5yAAMj8+G3NNi0iIiLSzNQmebgS+B9wIfBXY8wFVRWqlDj83lq7sH5ClKYgJTOVmVkL2ba7xLWib6BmmxYRERFpZoy1tvaFjbkEeAY401q73Gd95cThqXqNsgkYPny4Xb16tdNhOCo3PYOMx54g9MK2GGPp0uYCep00yumwRERERJob49SJj2meB2vtP4C7gbeMMQPgiMTh2paYOIjPbNPX/T92HAjHGFj7v1TNNi0iIiLSjNRmnoc7jDHhntfW2vuB14D3jDHP40ocLK4ah6d99oswxtzRADGLH/LMNh2VmMCO0g4AhA0M12zTIiIiIs3IUZstGWPKgG7W2l0+6wywBDgPKMfVOfrpSvt1AbZZa+syR0STo2ZLvzh86BDrvrwbYyz9T5hDu051mWNQRERERKrhWLOlVrUoY4CbjDGHK63PBkrcy85V1DJE1kN80gRFtG5NflF3IsO28mPGck4cN8XpkERERESkHtQmeQC4uZr1BogHhlSzrfa9saVZieqcSOmhrRQeXAdMcTocEREREakHtU0eXgcKjvHY4biGd5UWqPWW/eyJDCAifD9/TXmeW6bNIDc9g7z1OfScOsXp8ERERESkDmqbPFzn2+ehNowxXVHy0GK1jx1IzjcZtO9xmKiiLd7RmGLnzHY6NBERERGpo9oM1fo/oLgOxy4CPqvDftIMRCUm0K6za3bp2NByb+KgGadFREREmq6jJg/W2rHW2txjPbC1dr+1dmydopImLyUzlUcOfkxpaQBRbQtYFROqGadFREREmrjaNlsSOSbT4icyoTyaVavepH2vw8S07sQlcb8iKl41DyIiIiJN1THNMC1SW54+Du27JgEQ3jeQ7PkLNOO0iIiISBOm5EEahGfG6YQzz6GkJJCIiDw6zLhcM06LiIiINGFqtiQNwnc41uKy3gQF/cSuoh2cPPXXzgUlIiIiIsdFNQ/S4Dr2HApAWUG2w5GIiIiIyPFQ8iANLmLDdkpKWhEelseDrz4HuPpEbFmyzNnAREREROSYHDV5MMacboxRkiF11nbgAA5uDQGgQ+k2b2fqyJj+DkcmIiIiIseiNknBp8BOY8wiY8z5xpjwBo5JmpmoxAS6RA8DYEBECd9rwjgRERGRJqk2ycMtwA/ApcAbwB5jzNvGmKuMMV0bNDppFlIyU/nr/ncoKmpFREQRaXFtNWGciIiISBNkrLW1K2hMZ2Cy+zEeCAPKgdXAMuBta+3ahgnT/w0fPtyuXr3a6TD8Vm56Bl998SYdYw6zZ30EI0+9QDUPIiIiInVjnDpxrfsyWGt3WWuft9ZOBjoCU4GXgL7AX4F0Y8x6Y8xDxpjR6ichHp4+Dr1iTwKgde9ish56RBPGiYiIiDQxdfqAb60tsNYus9b+FugKnAYsAEqB2cDHwC5jzGJjzFRjTES9RSxNjmfCuEFjzyQvP5SQ4BJCL5qiCeNEREREmpjjrh2wLl9aa+dYawcBcbj6SWQDl/BLP4nfHe+5pGnqOXUKUYkJBAQEEBQ2GIBDJVsqTCQnIiIiIv6v3psWWWuzrbUPWGtPAboBVwHvA+3q+1zS9PQbfDIAIa1+pqgg3+FoRERERORYNGi/BGvtbmvt3621U6y1DzXkuaRp6NormsMFUbRqVUb2d185HY6IiIiIHAN1apZGtWXJMgLLewGwedMKQLNNi4iIiDQVSh6kUUXG9Kf8vRVYC53b5rLtq68027SIiIhIE6HkQRpVVGICCb+/lgO7wwkIsGR8/LZmmxYRERFpIpQ8SKNKyUxlZtZCtu4pBiCgX4BmmxYRERFpIo6aPBhjPjPGXNoYwUjzNy1+Is/GXUv/jD2UlQXQvl0BD/b4DdPiJzodmoiIiIgcRW1qHk4FxjZ0INIyeGabHnL9dWw5EA7A91/9V7NNi4iIiDQBarYkjcoz23RUYgI7yjoBEN4/ULNNi4iIiDQBSh6kUXlmmwa46lczKClpRUR4Hgwd4nBkIiIiInI0Sh7EMUHBwZTYvgBs+n6lw9GIiIiIyNEoeRBHhR0IcT0p/YHbly4CNGmciIiIiL+qbfJwtjHmBWPMH4wxpxpjWjdoVNJi9BuSTGFBMKEhhZTlZXs7VGvSOBERERH/06qW5boAVwCXu19bY8wGYI3vw1q7q94jlGatfXISpTkrgBxOLAkhe/4CTRonIiIi4qdqmzysBL4EhgLJQAcgxv240FPIGLODiglFmrV2Qz3GK81MSmYqX+bncFlr6NqliOXRloezFnJhwLma+0FERETEz9Q2eci21v7R88IY0wtXIuH76AV0cz/+z13UHsM5pAWaFj+RCeXRpH3/Oq2jChgc1I0r4y4mKl41DyIiIiL+pk4dpq21m621b1lr77LWTrHWRgMdgTOBPwKvAj/UY5zSTHn6OAQHxwAQEh9F9vwFmjRORERExA/VW62AtXYf8JH7AYAxJqy+ji/Nk2fSuIBePcj+KpPIyD10vuYq8tbnqN+DiIiIiJ9p0CZF1tqChjy+NH09p07xPs8v6k5k2Fa2529n5NQLnAtKRERERKpUm2ZL6UBJQwci0qH7MACKD6+lvLzc4WhEREREpLKjJg/W2mRr7azGCEZattikkRSXBBEedojNOdlOhyMiIiIilWiGafEbO1PfoaigBwBfr0oFNNu0iIiIiD9R8iB+IzKmP2XL1wPQo80+9ny7RrNNi4iIiPgRJQ/iN6ISExh62RUcPhRKSHApa5a+ptmmRURERPyIkgfxGymZqVyd/TRbd7heFw8MYWbWQlIyU50NTEREREQAzf4sfsQz23TGk09hYyLo3KmAR6KuoEf8SKdDExERERFU8yB+xDPbdMLvr2H7gQgCAiwZ/0vVbNMiIiIifkLJg/gNz2zTUYkJbCvpCED4wGDy1uc4HJmIiIiIgJIH8SM9p07xdo6+4oIZlJYGEhFxkIATEh2OTERERERAyYP4qd3vvEdBfncAPlu+FNCcDyIiIiJOU/Igfikypj+lX/8MQM/We9m7Jk1zPoiIiIg4TMmD+KWoxAROuPhS8g+HEBZSwrdv/EtzPoiIiIg4TMmD+CXPnA/bdhgAigeFas4HEREREYdpngfxS545H9IXPoXt55rz4eG2l9MrfpTToYmIiIi0WKp5EL/kmfMh8dpr2JbrmvNh3Rf/0ZwPIiIiIg5S8iB+yXfOhy2lnQEIHxjAwR/XOxyZiIiISMvld8mDMaaPMeY6Y8x7xpjtxpgSY8whY0yaMeZuY0y7KvaJcW/7yhhzwBhTbIzZaox50xgzrhbnnGSM+cQYk+s+10pjzOUNc4VSG75zPsycdhXFxcGEh+VTGqfRlkRERESc4nfJA7ASeBRYC0wGYoAJQAZwG/CNMaajp7AxZhKQDdwALAXGAPHAX4BRwEfGmHuqO5kx5nbgLWCfe98TgTRgsTHmuXq8LqmjVkFBlAfGArAl50uHoxERERFpufwxeQBYaK29yVq7ylq7yVq7wlp7KfA/oC9wtU/ZDriuY6a19n5r7Rpr7Q/W2peAs4FS4FZjzOjKJ3GvmwesAaZZa9OstVnW2quBt4EZxpjLGvZSpTbaHAgCILTVRu564wVAk8aJiIiINDZ/TB5mAQ9Us+0b97JzpfWHgJTKha21GcBX7pcXVnG8O93Lx6y1ZZW2LXAv76gxWmkUvQbHc2BPOIGB5bTJ3+TtUK1J40REREQaj98lD9baf1trN1deb4wxuJoUAXzss+kVoEcVH/49triX7SsdrzPgqY34qIr9vgSKgP7GmGG1DF8aSFRiAhGRgwFIDLNkz1+gSeNEREREGpnfJQ+VGWNCjDFJwMvAScC91tplnu3W2mJr7aEaDtHNvcystH4Yrus/XFWyYq0tATa4X46oY/hST1IyU3ns8GeUlAQS1baAb2IjNWmciIiISCPz60nijDErcHV6BleTpVOstV/VsEvl/dsBI4FC4IVKmz3tXXbWcIjtQBzQr7bnlIbhmTTuq+Vv0rHfYfq378DFcRcQFa+aBxEREZHG4u81D9OAIcBUIA9Y7h6StbZx3wiEALdYaysnCW3cy/wa9i9wL9tWtdEYM9MYs9oYs3r37t21DEnqwtPHoffAUwBoG13Euocf1aRxIiIiIo3Ir5MHa+1ma+06a+1SYBywAtdwrQ8ebV9jzEhcw7W+gWvo17ownlCqie9Za+1wa+3wTp061fEUUhueSePixp1Jbl4EQUGlBP3qLPLW5zgdmoiIiEiL4dfJgy9rbTlwl/vlH4wxUdWVNcYMAlKBD4HfWGur+vB/0L0Mr+G0oZXKikN8J41r22E4APmlG+k5dYqDUYmIiIi0LE0meXDztFEJBoZWVcAYE4sraVgBTLHWFldzLM9X1l1qOJ+ns/WGGspII4s74XTKygKIDN/Nrm1H9HUXERERkQbiV8mDMSb2KJOy+fZPCK5i/yG4JpJbCVxgrS2q4VjfAOVAhDGmVxXHCsI1IR3A6qPFLo1n/wcfk3+4KwAffvI6oAnjRERERBqDXyUPuIZiXWSMCa1m+2Cf5+t9NxhjkoFPcc3ZMN091Kpn25nGmBd9y1trdwGfu1+Or+Jcp+BqtvSTtVbJgx+JjOlP6aqtAPRus4e9a9I0YZyIiIhII/C35AFcMc2oZtvt7uUKa623p6wx5kRcE8e9DVxaxYRxPfhlQjhfnj4U1xljAittu9G9nFfbwKVxRCUmcMJFl5KfF0JYSAnfLH1NE8aJiIiINAJ/m+eh1L182BgTDbyJax6GGGAOcCauGaMv9+zgThw+AFoDScDXrsmoK+hQ1cmstZ8YY+4C7gRSjDHzgGLgOmAysNhau7herkzqTUpmKm9k/4cpuzsQOxBsXDAzsxZyYcC5TIuf6HR4IiIiIs2WqXogIucYY8YA04GTcfU5CMc1x8P3uEZQetxae8Cn/FxcH/6PZpO1tk815zwPV03DUCAQWAs8ba1dVNu4hw8fblevVuumxpKbnkHmY48TfEE7AgPLiQqdSP/Tq6pcEhEREWl2jvimvNFO7G/JQ1Ol5KHxeCaMi50zm/czlxDTMY89P0Yw8rQL1HRJREREWgLHkgd/7PMgUiPPhHFRiQlsNa5Rl9r2KWJ/9g8ORyYiIiLSvCl5kCbHd8K42RfN4nBBW4KCSsnt1trhyERERESaNyUP0qRtWbKMYPoBsHfnCkBzPoiIiIg0FCUP0qRFxvSnfOnHlJS0ol3rPLI+/q/mfBARERFpIEoepEmLSkxg8OzrObApBIDN65drzgcRERGRBqLkQZq0lMxUZmYtZP2e3QBE9Srg+sxnSclMdTgyERERkeZHyYM0adPiJ/Js3LUM/zGfvXvDCQwsZ2b4CE0WJyIiItIAlDxIk+Y758PqknIAysp+YF/adw5HJiIiItL8KHmQJs13zofCyH4UFQcT0bqInLVpTocmIiIi0uwoeZAmzXfOh3lTr4SgIQDkRx5wMiwRERGRZknJgzQr7XIDKS+H8JAt3PP684DmfRARERGpL0oepFnpFjeYAzsjCAiwdC7c6u0ToXkfRERERI6fkgdpVqISE+jYZRgAQ9oWkfXwI5r3QURERKSeKHmQZiUlM5W/5r7DobwQQkNLyE7uzMyshZr3QURERKQetHI6AJH6NC1+IhPKo/nqP6/Rehj06R7Ms3HXEhWvmgcRERGR46WaB2lWPH0ckiecT3FJIK3bFfDNS4vITc9wOjQRERGRJk/JgzQrnnkfugwbxqbcTgAEjOpB3vochyMTERERafqUPEiz4jvvwxnjf0N5OURG7iDs1FEORyYiIiLS9Cl5kGarcPlKDud1JSAA/vvhPwHN+SAiIiJyPJQ8SLMVGdOf0q+2AdCnzR52fvON5nwQEREROQ5KHqTZikpM4ITfXMah3DCCg0tJe3+p5nwQEREROQ5KHqTZSslM5ersp9m0tRiAkEGGmes054OIiIhIXWmeB2m2PHM+ZL3zCIX9u9I6spBbys8hOX6i06GJiIiINEmqeZBmyzPnQ9xNN7L2QAgAe3Z9ozkfREREROpIyYM0W545H6ISE9gZ2pOyMkNUl8NsXpfpdGgiIiIiTZKSB2m2fOd8uDKwA/mHe2AMfF++BdCwrSIiIiLHSsmDtAiRMf0p//InAHpH7WH716s0bKuIiIjIMVLyIC1CVGICJ1z5Ow7sCadVq3K+++xtDdsqIiIicoyUPEiLkJKZysyshfy0NR+AyIGlXLP2KQ3bKiIiInIMNFSrtAieYVu/T13Aof49aR1ZyPWcxigN2yoiIiJSa6p5kBbBM2zroDmzWZ1nACguWce+tO8cjkxERESk6VDyIC2C77CtByP7UFQUQkTrQn7I/Mbp0ERERESaDCUP0iL4Dts6b+pvKTnQHYDi1ru4fdliQEO3ioiIiByNkgdpkWIHnUhJSSARYXspPvi9t1mThm4VERERqZ6SB2mRug4fTnFhXwBGB4SQPX+Bhm4VEREROQolD9IipWSm8lLhOsrKDJ075fPNwHBmZi3U0K0iIiIiNdBQrdIieYZu/XrFm3Toe5iYLp24OO4CouJV8yAiIiJSHdU8SIvk6ePQb9AYrIV2vfPJePIpctMznA5NRERExG8peZAWyTN0a8zoMWzf346AAEv5hETy1uc4HZqIiIiI31LyIC2S79Ct8cmTAQgN20j7CeOdDEtERETEryl5kBYvaN168vLa0apVGW/8+++A5nwQERERqYqSB2nxImP6U/zNAQD6tdnF7m/XaM4HERERkSooeZAWLyoxgWG/upi8g6GEhZTw7ftvas4HERERkSooeZAWLyUzlauzn+anTcUARAy2XLP2Kc35ICIiIlKJ5nmQFs8z58P3qQs41LcHrSOLuI5TOCl+otOhiYiIiPgV1TxIi+eZ82HQnNmszDMAlNl17FmT5mxgIiIiIn5GyYO0eJ45H6ISE+i0sZz8gnDCwov56Iv/ABp5SURERMRDyYO0eL5zPlww5dfkf1cCQNfeBexdk6aRl0RERETc1OdBxEdUYgLDy8pY99MbREYWsTo1hZEaeUlEREQEUM2DSAUpmalc88Mz/PRzGQCh8QHMWrdQIy+JiIiIoJoHkQo8Iy9l/ecRDvfpRmREETe1GsuJGnlJRERERDUPIr48Iy/F3Xwjqw8GAlBYnM5ejbwkIiIiouRBxJfvyEu5kX0oLAolonUR369d7XRoIiIiIo5T8iDiw3fkpXlTf0vZgR4AlLfexu1LFwEaulVERERaLiUPIjWIGzySwoIgwsMOEXoox9usSUO3ioiISEuk5EGkBp1OGIqxcQCcEmH4fv4Cb7MmERERkZZGyYNIDVIyU3miYCUFBcG0bl1I9gldmJmloVtFRESkZdJQrSI18Azd+tXbrxI2Avr3CWRh9Ew6xg91OjQRERGRRqeaB5EaePo4DD/nVxw8HEJYRDHf/HcJuekZTocmIiIi0uiUPIjUwDN0a4ehyWzfFAZA6yHlLHkrBdDISyIiItKyKHkQqYHv0K1nj5nIodwwQkKKKY4O1MhLIiIi0uIoeRCppfbJSUS1GQFAcvtCMh97XCMviYiISIui5EGkllIyU/lr7jvs3RdOUFAZ20d01chLIiIi0qIYa63TMTQLw4cPt6tXr3Y6DGlguekZfPPSIqImBFJaGkC3dufR+6STnQ5LREREWhbj1IlV8yBSS54+DsMuu5It+yJo1aqcrLT/auQlERERaTGUPIjUkmfkpajEBLLpgbXQvm8+W9etdTo0ERERkUah5EGklnxHXro4uBN5h7oTEGD5kU2Ahm0VERGR5k/Jg0gdRMb0h49zKC839Gy3hx8++UjDtoqIiEizp+RBpA6iEhNIvnoW+34KxxjYvOUzDdsqIiIizZ6SB5E6SMlMZWbWQrJ3bKWkJJCoLvk8mvaqhm0VERGRZk3Jg0gdTIufyLNx1zJyQyk/bQwCYHybVkyNPcvhyEREREQajpIHkTrwDNsaO2c2b7ffR0FBGOHh+Xz9n1edDk1ERESkwSh5EKkD32Fb+4eNIOBQNACm1Trmvf4CoNGXREREpPlR8iBSB77Dtt495QoGxQ/nwJ5wgoNL6VO02VszodGXREREpDlR8iBSD9onJ9Gt66lYC3EdD5P23LMafUlERESaHSUPIvUgJTOVO3cvYev2CAICLPkj2zMza6FGXxIREZFmRcmDSD3wjL7UfsUuSksD6NI5n1ujzmVa/ESnQxMRERGpN0oeROqBp49Dwu+vYcPPIa51B79iwcL53u3qPC0iIiJNnZIHkXrgO/rSzvI2FBYEE9m2kA6l+9V5WkRERJoNY611OoZmYfjw4Xb16tVOhyF+4uu3UwgMWUVhURAs282Q669T52kRERGpL8apE6vmQaSepWSm8tDhT9ifG0ZoSAk7T+quztMiIiLSLDS75MEYc7YxZqsxRlUq4ohp8RN5dvC1BKw8iLXQu1c+8zpNVedpERERafL8LnkwxvQxxlxnjHnPGLPdGFNijDlkjEkzxtxtjGlXzX4RxpingHeA7sd4zknGmE+MMbnuc600xlxeH9cjLY+nj8MJV/yOrD2RGAPbd33OvrTvnA5NRERE5Lj4XfIArAQeBdYCk4EYYAKQAdwGfGOM6ei7gzEmBkgDxgLTj+VkxpjbgbeAfcAY4ET3sRYbY56r81VIi+XbebpkYxlFRcG0aV/AB1+4mi1p5CURERFpqvwxeQBYaK29yVq7ylq7yVq7wlp7KfA/oC9wdaXyg4F3gaHAqtqexBgzGpgHrAGmWWvTrLVZ1tqrgbeBGcaYy+rjgqTl6Dl1irdz9NRJv+LQt+UA9O57iM0rV2rkJREREWmy/DF5mAU8UM22b9zLzpXWp1prr7PWFhzjue50Lx+z1pZV2rbAvbzjGI8p4hWVmMCI835N7u5wgoPKyMp4z1srISIiItLU+F3yYK39t7V2c+X1xhiDq0kRwMeV9ik/1vMYYzoDo90vP6qiyJdAEdDfGDPsWI8vAq6Rl67Ofprvc3ZTVm7oEH2YB777p0ZeEhERkSbJ75KHyowxIcaYJOBl4CTgXmvtsno49DBc13+4qmTFWlsCbHC/HFEP55MWaFr8RJ6Nu5bhPxawcWMoAOe0DmVq7FkORyYiIiJy7Pw6eTDGrAAKcXVgHgScYq29rZ4O72l0vrOGMtvdy371dE5pYTwjL8XOmc2/2+6lID+E8PDDrEz9J7cvW+wtow7UIiIi0hT4dfIATAOGAFOBPGC5e7jW+oi7jXuZX0MZTx+KtlVtNMbMNMasNsas3r17dz2EJM2N78hL/cJGEBEyHICgkCx25X3nTS7UgVpERESaglZOB1ATn+ZE64wx/wY+wzVcaxhwcyOE4Jn6u8oJ56y1zwLPAgwfPlyT0skRek6d4n1+95QrAPj0jfW0brOT82jjrZVQB2oRERFpCvy95sHL3Sn6LvfLPxhjoo7zkAfdy/AayoRWKityXFIyU/lX4VZKSwPo0jmfdcmdmZm1UB2oRUREpEnw65qHKmS4l8G45nT45DiOleNedqmhTDf3ckMNZURqbVr8RCaUR/PVOyl0PKGQ2IGBnNX5CrrHj3Q6NBEREZGj8quaB2NM7FEmZfPtnxB8nKf7BigHIowxvaqIJQjXhHQAq4/zXCLALx2oR5x9Ifv2hxESUkJm+rs88tRD3u3qPC0iIiL+yq+SB1xDsS4yxoRWs32wz/P1x3Mia+0u4HP3y/FVFDkFV7Oln6y1Sh6kXng6ULdPTuLHfW0oK3PN/RBxeKc6T4uIiIjf87fkAVwxzahm2+3u5QprbU41ZY6Fpw/FdcaYwErbbnQv59XDeUQAVwdqT+foa6+5maKCOADiBgSz9tG/qfO0iIiI+DV/Sx5K3cuHjTHzjTGjjDF9jTFnGmP+C5wDbAEur7yjMaaTMaYr0MlnXVf3o1Pl8gDW2k9wJRBDgRRjTJIxJs4Y8xQwGVhsrV1cr1co4paSmcoTxas4cDCU8LBidp3aQ52nRURExK8Za/1rhFFjzBhgOnAyrj4H4bjmePgeSAUet9YeqGK/jUB0NYfdZK3tU8M5z8NV0zAUCATWAk9baxfVNu7hw4fb1avVukmOTW56Bt/+/e+0/r9gAgIsgeWnk3z2JKfDEhEREf9mjl6kgU7sb8lDU6XkQY6V7+zTr323hOFd8igsCKJft/NZuG0td0+5gtz0DPLW51SYL0JERERaPMeSB39rtiTSYvjOPr01sjd5h9sQGlbCuuyPyC76Sh2oRURExO+o5qGeqOZBjte2TRvYsu5pAgMt330HQ9L3qAO1iIiIVEU1DyIt3ReH1rF6t2vC80GDWrGyXyt1oBYRERG/0tRmmBZptqbFT2RfaS/S1i2hbcd8hvTvwaWx59M+Psnp0EREREQA1TyI+I3c9Ax+fPhR+kWfSUlJIO26HmbVu69r9mkRERHxG0oeRPyEpwN1n1NOJXtrRwDaJZcQsP8ndZ4WERERv6AO0/VEHaalvn3yxiO0abONAwdCCUvdStxNN6rztIiIiIA6TIuIr5TMVBYVbyC/IJi2bQvZdmovdZ4WERERxyl5EPFD0+In8lj8TIo/L6C83NAnOp8/tzmLjPV7APV/EBEREWcoeRDxQ54+DsMuu5LMLaEAFJV+hf15tfo/iIiIiGOUPIj4Id/Zpzd06MuhQx0IDillTKf2fP/QAk0eJyIiIo5Q8iDih3pOneJNDhIGdOLlwu0UFgXRrkMBm09R/wcRERFxhkZbqicabUkaUm56BqteeYl2Y8EYMGWncML/TXE6LBEREXGGRlsSkap5+jiMuPgyVu+KBKC4bCU/ffEZty9b7C2jDtQiIiLS0JQ8iPg53/4PWyN7c+hQJ4KCy/h5x8fkFKxSB2oRERFpNGq2VE/UbEkaS96BXL7738OEhhayeUs43T75WR2oRUREWhY1WxKR2nln8xcsPVBEWZmhV898fhzeXR2oRUREpFEoeRBpYqbFT+S2pEvZ/00QALGxJdzRfrImkBMREZEGp+RBpInx9HEYOfnXrN8WQUCA5WDB55pATkRERBqckgeRJsa3A/XaNtEcOtSOkNASxnVtx7oFj6r/g4iIiDQYJQ8iTYzvBHLxAzvxUuEO8guCaRtVyL6xPZm5Tv0fREREpGEoeRBpwqbFT2RBwlWUfHyY0tIAunXL5/8FD1P/BxEREWkQSh5EmjBPH4cTfvc7lm80WAsRkdl03vGD+j+IiIhIvVPyINKE+fZ/2Nd9EEX58QCMjC7jm5cWqf+DiIiI1CslDyJNmG//h4SYjvytaCWbt4QTGFhO2JgQbk1frP4PIiIiUm+UPIg0E9PiJ/Ls4Gvp8ulmdu8JJyS4lItah3NWj5O5fdliQH0gRERE5PgoeRBpJjx9HOJuvpF/spe8vEjCwgpJ+/gJNuSvUh8IEREROW5KHkSaCd/+D33Ch5N0+rUUFoYQ2foQl9iOZD30iPpAiIiIyHEx1lqnY2gWhg8fblevXu10GCJeKZmpfJH2KRdFBREcXMqOreG8GLaLDq168NSFt5GbnkHe+hx6Tp3idKgiIiJybIxTJ1bNg0gzNS1+IvMSLyf/kyJKSwPo2iOfa8oHsrd0q5owiYiISJ0oeRBppjwJwrDLryR1TyDl5YY2HbcxfU9HsucvUBMmEREROWZKHkSaKd8+EHvCw/lgVwjWQp+BheSM7MHMrIUaxlVERESOifo81BP1eRB/l5uewVdvv0rHEcUAFOUn8i5l3D3lCvV/EBERaVrU50FEGo6nCdPISb/mm42hAISEp9Nj53r1fxAREZFaU/Ig0gL4NmHa3Lk/hYfjARjep4iv3npV/R9ERESkVpQ8iLQAPadO8SYHCTEd+VvxSn5cH4ox0GFEMa+ueI2UzFTNRC0iIiI1UvIg0sJMi5/Is3HX0nfFFn740ZVAjOhbRK+cA2QXfaVmTCIiIlItJQ8iLYwnOYidM5ulHfdQcHgIAMFhaZy3q72GcRUREZFqKXkQaWF8+z/Ehoxk24COrNkUBsCg2GK2nNiTmesWcs0b9wBqwiQiIiK/0FCt9URDtUpTlpuewVepr9FhWBHGwIE9PXg6IJtnB1+rmggRERH/o6FaRcQZ3mFcJ07nvztDKS83tO24lUsPdybroUeUOIiIiIiXkgeRFq7CMK7hIfxnVyvKygLo3iufAxN68f8yn1ETJhEREQGUPIi0eL7DuD514W3ckHwRhz4qpbi4FZ0653NVSCeKi3ZrFCYRERFR8iAiv/AkCMMuu4JX9pdQUBBCZOs8Lg+PJO2559SESUREpIVT8iAiXr5NmErDolicd4iDB0KJiCgm5Kxwnv/iZU0mJyIi0oIpeRARr8pNmB5KmEHwu9vZtTuc4OAyRveD3hsOajI5ERGRFkrJg4hUyZMcDLnhehYH7iLvYD8CAixBoWu4dH8n1t4/n9g5s3l4wzfe8qqFEBERad6UPIhIlXybMA0MHcnOQb1Y/VMI5eWG7v0KKJrYh/n/e1q1ECIiIi2IJomrJ5okTlqC3PQMVr+8mLDTQggNLaGoMIjMdQUkrc1VZ2oREZHGo0niRMS/eWoXhl96BS8cOsye3HBCQksYmhzEpmHdmbluoeaDEBERaeaUPIhIrfg2Y+oeMZThA85n9w/hBARY+g8u5mr6kFe8Q02YREREmjElDyJSK74jMd3UbxjrF/yNUadfyPs7QiktDaRt+z3MCIni22eepucF5xOVmKAhXUVERJoZJQ8icsx8ayE2h4fwj32l5OWFENm6iMiJEaz5fiVL//2sOlOLiIg0M+owXU/UYVpastz0DNY++jd2n9qDXj3zATiwO5wdq7bSd2cRcX/5Iw9v+Ia7p1xBbnoGeetz6Dl1irNBi4iINF3qMC0iTZPvfBCvROxi055oiopa0bZTPv3O7Mj+mC7M+/RR1UKIiIg0A0oeROS4+DZhig0ZydSL/x8dC4exf0c4QUFldDnN8H8RvZjyWT5Z9z2oieVERESaMCUPInJcfDtSe5ol7X/rP5yQOJXPdoRTWhZA+16HiT6rB/v7dmLeJ6qFEBERaaqUPIhIvfLURLRPTmJfm750jDiHg/vCCA0tocuYQM6J6slp31uy7rlPozKJiIg0MUoeRKReVR7Sdc9zi0mMncqXO8IpKQ2kXbd8ThzRli2n9GPtv15m2VvPqSZCRESkiVDyICINxlML0WFoMrvb9KVb1CT2bY2gVaty+sUVETp9CK3/l8X5H+1XfwgREZEmQMmDiDSYyrUQO556nhHDLiB1RzAFhSFERB6m7YRA2g8dSF4rNCqTiIiIn1PyICKNwndUptI2Axh+xi0c3NMTa6Fj/wLaXNSdCaXdmfLpIdVCiIiI+CklDyLSKCqPylS4Pofwd9OICv0/tux3NWXqFF9I9Nm92Rfbjbs1KpOIiIjfUfIgIo7w1EQMGDOO1SF9aGVHcyg3lNCwErqebJncvTtnbgzXqEwiIiJ+RMmDiDiicn+Iwn+8SeLAC/loRygFhcFEti3khKEB5J0fw5ovPmTpv58lcs3HbF32FtnzF/D6zvWAEgkREZHGpORBRBznOzfEoTb9GdhjKrvTQykpCaR9hwK6TQihbf4Wugd15KfFL9HzgvNZ0T5HTZpEREQamZIHEXFc5VqInx59glETprHwwGG27u5GcXEgbTsUMPjEQIouHsyalR8z8rs8NWkSERFpZK2cDkBExJfvqEzRG4YxecoVbHhjCT9kfUubuDKiogqIGhdEZG4fNkaVcPDVf7AqYCeRaz5nK+3Z8uZSlp88hBtxJRJ563PoOXWK05clIiLSLKjmQUT8SuVRmXLTM9j977cZddZ0Fh48RN7uvhQVtqJ1VCEDk8to/ZsYQnO20DE/iI1q0iQiItKglDyIiF/zrYkYdqAPoe+tonvJiWR+H8Shw+GEhJTQaUg+CeM7sH/yENI+ekczVouIiDQQY611OoZmYfjw4Xb16tVOhyHSrG1ZsozImP7ePg53Tb6Mb179BwdKcojqmo8xrnIH94aR/2Mh/4s6wMYuAfx1TyK7v1hO3F/+yMMbvvHWaKhJk4iINFHGqROr5kFEmozKTZoOZq6l/N1PGD70QhbtLufAnh6UlATSpkMBXUdZLoiJ4pKCLvy08UfvMTTcq4iISN0peRCRJsu3SVPMoW5EvPsdXQpHkpEVxJ7ccFq1KqfH/2/vzsOkqs48jn/fqupulqZBdhAQEAUiqMQFcTdRo6gILphVCZPJZiaJcXSGGU0yMYkzcdTETNSYiTGazSQiqIlGk9HEREFBDEtcosgmIMi+9Vb1zh/3dlMUVc2t6mq6uuv3eZ5+qrrOuaduHd576LfOPfcO3kXvc+I0XDGM53//S/rXd2teG1G/frESCRERkTwoeRCRDit9JuLyAaMYfd2XGHHpNN487EgmjrmMHQ/vYvXrXajdU0H3bvX0H1fL0Wf1oO4jY1n093n0r+u6zyLrN753F6/c/C2qRx2uS7+KiIhkoeRBRDqFzHtFvHbLbRz7mc/ws37vMmbYdDY97Wxe3Y3Gxhg9e9Yy6JgGjn5fTXMicfqrsPHPzzWfRKrTm0RERPan5EFEOp3005lGV00kHo9Ts+IdRtQO4o5tO6nfOJZNq7ruk0hMmlRN3eWHsfzcw5n3wL303pbU6U0iIiIZdLWlItHVlkRKV/pVmm6/6785+bllDLl0Gj99eTanHXE6O/asoGZQPRUVyeZtGhribH63itSKPWzYtYnRr29j+IwruabucV29SURE2lu7XW1Jd5gWkU4v/Q/6yweMovq6D9Dr6PFUrHud2BPPctSl0/jJCw9x/imXs/b1F+k2KEn3HnUMGLQbBsEgerPzpEEsXr+A83fXsHbZS1SG7VUv+r/mO1svH9iL88Ob0s2e+yAzb/y6kgoREelUNPNQJJp5EOl4Mu8bce3I43j15m/R99RTuLXX63yII6lPrqVH3zoSidQ+227f2YW69TE279jDwFc2csTFU7l99SNMf74eB341qZLPxo/XDIWIiLQF3edBRORgy7xvxM433mTMrOsZdfWnOXL7YKqemM/YmuN5Zv562HQMG1+uZNO73UgmY9RU19Jv1G5GT3B6frgvq+PzObv3ELYeN5ydQ/oydEMy5wLsZ5//I1sXL2Hr4iXce9MNgNZQiIhIx6CZhyLRzINI55JrncRPFv2amRdfw8s/u5/k4Gp8QIxePeuprGzcZ3t32Lmriq27E1Sua2RT7W6GvrKRI6Zfzq2r5uacoZg990EuufgKQKc+iYhITlrzICJSSnKtk9jJZioSCXq8/S59R4zm3xJ/5daG83j94UfYMLofvXtUUdUnRfeaWnpU19Gjug76wwDiMHEgK/bM45w+Q9h6CtiWekbv2sPaVxc0r6FY2nMd77n5WzhQPzjF23MeYc1DD/PcyUdx+ew5WCyGp1I8uWyBkgwRETnoSm7mwcyGA1OAycAxQF+gFngTeBS4zd235Nj2IuBLwAQgDiwD7nL3Hx/gPQvaLp1mHkTKR65ZiUf/9CDHrovhwEOTujDzqMtZNe8ZrHcF8d5Ojx77r51osqe2grodFWyvha6bktiOBtb5Tka8vplRH/4Q19Q9zu1V57PivvsZPuPKfWYv/jo4xUWnXbE3yRgwqjnJ+H5sq9ZaiIh0PlrzkGYe8G2CP+CnAKOAc4ElwA3AQjPrm7mRmd0IPAJsBs4ETgReBu4zsx/kerNCtxOR8pXtztaHTp1C5SEjGDPresbOup73bBvA4H6D6LV0NYenBnEPGxlUdzJbfrmFVxfCxqVd2LS6G9u3dyGZjNG1SwO9+u1m2NDd9Du2jr6npRh/ejeqPzGElfHn+aT142+bXqL+8tNZunQ+Y/f0YnfvGpIVFSQafZ87ZVssxor77sdisf1udrdm9pxgNmP2HO696YbmtRdzPvsZrcMQEZEDKsWZh/XAr939c1nKngHOAG5096+nvX4G8AywCDjB3ZNpZY8AFwFXufv9Ge0VtF02mnkQkUxRZih+NamSfz/tn1j8/e/T0LuajYOr6BuLE6+GRHWSbl3ricdbHqcbG+PU1iZo3BNjV0OMqh1O196DqH1tBettN4eu3sXwyRcwq/Epbu+SffbimdHGWa/5AddhPHP3nZz56c8CwelS5x51fNZTqTLr6bQqEZGi0pqHNJ8CXspRtpAgeeif8fpXwsc70hOA0G0EScCXgcwkoNDtREQOKNe6icolLzBmRvBH9ri5DxKPx+m6bSdDxx/DE8kXmf6XvYuprxk6hb/Pfpg1h/diQKI7lSMG0Fi7Gbqn6FKVpKpLI4lEkurqJFRDr+Z3fJPu/aAPVUAVW3mef071YEXdPJJXjuPVXQs5q+cQtp0KXptiaEWKraOTxOqTjN7VyKq3lxGv7kZtfR1La/auw3hltDEwfL50UiUfCGc5hs+4cp/1Gpn13vjeXc3JyL033dDqZKSQejqFS0Sk9Upu5iEXMzPgT8CpwDR3nxO+3h9YR3AK1jB3X52xXQWwA6gCjnf3ha3ZLhfNPIhIodJnKNL/sH78gXsZuX5r81WePjrhsv1mDVLAo5O688H+Z7B+wTy2DOxGr3gllQNr8Ibt0NWpqkhSUdmYc73FgaRS0NiYINkQoz4Zg7oYqQbYDVTVGpU1h1C3YTNbKlIcsi2FNaZY2z3Foe+msMZG3uyV4si1jcTrGzniEzO5fd4PuWR+LbjvN+Nx7bCLI82MFFKvqQ+bb+j3sZlAy8lIWyc0nbFeKe5TqdcrxX0q9XqluE9tXS9jBlczD7mYWRUwBrgOmAR8oylxCB1HkADsykwAANy9wcyWA2OBEwhmL1qznYhIUaV/Az7zxuYzMjlt0hnNScVONuOpFMNnXImnUozbNogxs4L/YFY+cC+1Lz7JuDDJOC0tybim7nFu53xW3HM/i0d05cjaGmpOP5l3ly5mXb8Y/evj0CVGbbXRHYhVgFc5lfEU8coUFYkk8XgquBRtJXRN2+9Dmp/tgr7QL61s78K0BL0BwutJrdv2az44tieM7UkqBe9LxUgdH8NTxjQ3lidfwK8ax9/rF3D+gEPZdaHhSTjaYMdQ8JRxetxZuu4lKqafwtI1C5lUdShbJwIpGNHF2VoNpJzjzFny6gt0v+wcli6bz7Dug1my4E/UnHcqG7a/wQu/uB8ctgyNsa52O1teeJHeE09kY+86Ftx3L47zzsg4G+t28e6fnqX/Gaez9pBdvPSD7+MOK4fEWXT33bg7q96TYHtjPet/+ziDJk9mec+tLL7ju7g7fz8Met/+HXDn9fdWcGYqxeqfP8hhH5zOKzXvsOxbt+I4r44y+v/Xf+PA306s4Gxg5f0PcNjHPsqymvWM/c9bgrIjYEABsz8Hq14p7lOp1yvFfSr1eqW4T21db+viJbx2y22Mvu5L+/9HchCV9MyDmT0PnBT+uhC42t3nZ9T5HPBdYLm7H56jnT8A7wNucffrW7NdLpp5EJH2knmn7E+lejV/i/Wrd97Y5x4V2WYvsn1j/+Uzv8grN38LA54bGeO4dZWkqip4bWgFI7cl8KoEG3rH6VcXp8uQgdRt2sDOrlCdAotDQxV0MYjFHSqcRMyxRIpEPEUs5sRiKazdvjfrvNL/S9/7fN+O9n2eZJRl/ZPAWi5zx7Odfm3pb9aCUq9Xwhyw8DO4ZTwv5XqE9awpfvZul/V5S/WK0UaEepnMjCh/Qxez3o5lCSZOnt50wQ7NPOQwHegBjAa+ADxnZt8EvuLuTfPvNeHj7hba2RM+9kx7rdDtmpnZJ4FPAgwbNqyFZkRE2k76zMVNU2fsU3b57Dn73KMi1+zFO3ffyZhZwfT4uLkPAsH/TH1PPYU1yRc5eXktvqeWLQOMmlXBf6y/m1TJtaOD04KOiJiMAPx5JJy0HIjFeGFknAlvxyAeZ8mQGEdtquSQkyeyacEClveDw7bEIG6s7R1j0E7D4zE21Rh99sToOmQwtevXsaWr06vWIAY7ukKPBsMMdneB7g1GomcPGnfuoK4SKpMpUgkjkXIaK4yYQSLlJBNGMh48j6cgFQ8+f9ydVAxScSOecmI4Hgv+6oyx94+l4PWg08w9SIwMrOnPagNr+sup6fXm//rTX88UlKUnWi0lXdnrtfQHSQf/61mkjKzskeCTr9zJZbELmD7uwnbbj5JOHtJOJ/qbmc0lWPNwA8HM+T/n0VRLyWPB27n7PcA9EMw85Nm2iEibaymxmJlWNvXOu/a+fvR41syew5hZ19Pr6PGMu2lNziSj0GQktrweT6XYE3eq9gRJxvJDKpl2THA1qDEzruQ3q+Yy/uUgGVlXYYwJk5En05KWcU1Jy4tBvUVpScvv0tY8jG46hSu8V8Yrh1Vy1IZ4OLMCp6xKMOTSaax56GH+MqyRScvZv+zh/cuy1svSxqTlQd9GqRe1vabnJ69KcOglU3l79hz+MrSRSSuC/7rmjYCJbwX15o0wJq2JM/jiKayd+yjzhjZywgrDgBeGwwkrg3+gF4fBiWsTDL7gAtb+5re8eGiS41YF/3oLh8F7VwftvTQMjlsbZ+B5H2DdE0/x0uBGjg3LFg2F965NMPC8c1ifpayk6v3uKRYNbuSYsOzloezzfMLaBAM/cHDq5dPGgHPP5p0nf8+iwckW6sVLqt6ioTBhTXAcLBoSft5z38/6J//AosGN+5ZFrVeMNiLWO3ZN+DmGwIR1CQae837WP/UHFg1q3KesLesdv8q555ov0GvceNpTKd7nIatwpuE/wl//ycx6hc+3h4/dWti8S0bd1mwnItLppd/LYuaNX6fX0ePpdfR4pt55V/PzmTd+nSGXTOXQqVMYcsnUFuvtfONNxsy6nlFXfzpMMoL7YYzd0L35+bhtg7IkI8WrN2nz4ax56GGGz7iSxrg1n7WyZkAlQy6d1nyvjDUDKiOV1VVEa6Ot6w29dBqrfvwAQy+dxtsDK4m5Y+7UVQQzJ+bOmgEJhk2bypr7f8KwaRezun+CRDJJLJmkNpEi0Zgk3pBkdf84w6dMYe0DP2X4lItY1c+oqG8gUV/PnngDlXX1JOrqWdkXRlx4Ie/85OeMvHAyK/s6VbV1VNTWsbKvM/LCyWzIUrY7Vl9a9S6YzIo+TtWeOir21LHL6pufr+jjjLzg4NTLt42NP/1FWC9F1Z5aKvbUssvqmp+v6JMquXq7rY7K3bVU7K5lZZ8Uh19wPht/+gsOv+B8VvZJNZdFrVeMNvKpV5Veb3JYb3JQr+og1Tvqmi/w2i23sXXxkoP4v8H+SnrmIYum3qokuBv00wR3ngYY0MJ2g8LH5WmvFbqdiIjkKdei8MwZj3QtzYwUUi/9FK7My+XmmkHJZ3ZF9dRn6jP1WVvW63X0eEZf9yV2vvFm85c77aGkFkyb2WhgYq6bsplZDbAt/PU8d/9dxEuubieYRTjB3ReErxe0XS5aMC0iIiIiB0m7LZgutdOWJgE/MrMuOcrfk/b8DQB33wA8G772/izbnEKQALyVngAUup2IiIiISLkqteQBgn36RI6yG8PH5939zbTXm9ZCfN7M4hnbXBM+fi1Le4VuJyIiIiJSdkoteWgMH281s1vM7CQzG2Fm55jZk8BkYA1wVfpG7v40QSIwAfilmR1jZmPN7C5gCnCfu9+X+WaFbiciIiIiUo5Kas0DgJmdCVwBnAyMILga0k7gVeAx4Lvuvi3HthcTzBhMAOLAMuBud//RAd6zoO3Sac2DiIiIiBwk7bbmoeSSh45KyYOIiIiIHCRaMC0iIiIiIqVNyYOIiIiIiESi5EFERERERCJR8iAiIiIiIpEoeRARERERkUiUPIiIiIiISCRKHkREREREJBIlDyIiIiIiEomSBxERERERiUTJg4iIiIiIRKLkQUREREREIlHyICIiIiIikSh5EBERERGRSJQ8iIiIiIhIJEoeREREREQkEnP39t6HTsHMdgCvtfd+dCJ9gXfbeyc6EfVn8agvi0v9WVzqz+JRXxaX+rO4urj7uPZ440R7vGkn9Zq7H9/eO9FZmNkC9WfxqD+LR31ZXOrP4lJ/Fo/6srjUn8VlZgva67112pKIiIiIiESi5EFERERERCJR8lA897T3DnQy6s/iUn8Wj/qyuNSfxaX+LB71ZXGpP4ur3fpTC6ZFRERERCQSzTyIiIiIiEgkSh5ERCQvZnaemb1tZpq6LgL1Z3GpP4tHfSnZKHnIwswuMrOnzWyrme0ws3lmdlUr2jvVzB4zs3fNbLeZ/dXMrjGzeDH3u9SY2Sgzu8nM5pvZNjOrDwehh8zsfQW0d6aZ+QF+2uWaxweDmc2I8PmrC2i37OLTzIZH6Mumny9GbLPTx6eZdTezu4DfAoPz2K6oY2rYZoeP23z7s9hjathmp4nbAvqzTcbUsO0OHZ/59GVbjKdhu50iNltz3Jbq2Kn7PGQwsxuBrwGzgTOBOuALwH1mdqq7/2Oe7V0F3Av8GbgI2AhcCdwKnGtmF7l7Y/E+QWkws4uAOcBu4BvA74BdwEnAzcAlZvYNd78hz6YbgTdbKK/Lf287lD3AqhbKU/k0Vq7xmWY50JCjrA/BTY1ezaO9ThufZjYKeBxIAlcAv4y4XVHH1LDNDh+3+fZnG46p0AnittD4pMhjargvHTo+W9GXxR5PoYPHZmuO25IeO91dP+EPcAbgwEtAPKPskbDsyjzaOwKoB94GqjPK7gjb+3J7f+426ssZ4ef7UJay8QQDjANn5NHmmcCK9v5s7dynzxSxvXKOz+Hh5xveQp2ngNcJLywRoc1OHZ/AlDAuuqb1nx9gm6KOqeF2nSJu8+3PthhTw207RdwWGJ9FHVPDNjt8fBYQm0UfT8NtOnxsFnrclvrYqdOW9vWV8PEOd09mlN0WPn45j/ZmARXAD9x9Z0bZ7eHjdWbWLb/d7DB2kOUbC3dfAswPf73soO6RpCvn+KwDFpLjWyszGwOcDdzp4cgqPObun3f3PXlsU+wxFTpP3BbSnxpTcyukP9tCZ4jPfPtS42nLCjluS3rsVPIQMrP+BJkewB+yVPkLwYFxuJkdF6G9ODAtV3vu/hbwFlANnF/IPpe4nwGHZgn6JmvCx94HaX8kTbnHp7uvc/fj3X1djiqfI5ha/tFB3K2S5u75nhJX1DE1bLPTxG2+/YnG1BYV0J9F11niM9++1HjaoryP244wdip52Os4gv7Y5e6rMwvdvYHgfD6AEyK0dyTQK3ye6xy/ptejtNehuHu9u+9oocqg8HFpnk1XmNkXw0VD681srZk9Y2ZXm1lVgbvbkfQws6+Y2UIz22Bma8zsCTP7qJnlczyXdXy2xMx6EJwD+hN335bn5uUen+mKPaZCGcdtG46pUN5xW6wxFco4PnNp5XgKHTw2CzxuS37sVPKw1+Hh4zst1GnKqkfm0V7S3TcWob1Ow8wOASYCtQQLd/IxGJgO/BfwfuAjwHrgf4A/h213Zu8FTgT+neB80E8AceAB4FEzq4zYjuIzt6uAHgQxla9yj890xR5T09tU3KZp5ZgK5R23xRpTQfGZTWvGU+jEsdnCcVvyY6eutrRXTfi4u4U6Tef/9cyjvZbOGcynvc7kGqAK+JK7t3RwZFoD/AfwjTDzBlgGPB0O8NMIbtd+eTF3toT8DbjW3W9Lf83MngLmAZMJruZwXYS2FJ+5XU2wiDLfb3DLPT4zFXtMTW9TcbuvQsdUKO+4LeaYCorPbAodT6Hzx2au47bkx07NPOTHwsdiLfgpdnslz8wmEiza+TXw7Xy2dfc33P2raYNIupvCx8vMbHirdrJEufsLGf/JNb2eBL4Z/nq1mXUp0luWY3yeA4yhgG/Jyj0+C9QWMVZWcduaMRXKO27bYUyFMorP1oyn0Lljs7XHLe08dip52Gt7+NjSKvOmAWR7C3Uy2+tapPY6vPCKC48Bvwc+UuSrLiwmuB40wKQitttRvBQ+dgUmRKiv+MzucwTfds0pcrvlGJ/FHlPT6yluafMxFcozbpvkO6aC4jNTW42n0IFjM8JxW/Jjp5KHvZpuQjKghTpNC1uWt1Ans724mfUrQnsdmpmNJjhQngemunt9MdsPvynaFP7aYc+BbIX0Kc8on1/xmcHMDgMuBO5u4coYBSnT+Cz2mJreZtnHbVuPqVC2cdsk3zEVFJ/N2nI8hY4bmxGP25IfO5U87LWQ4E6S3c1saGahmVUAI8JfF0Ro73Wg6coCY3LUaXo9SnsdlpkdBfyR4BzSS929oDtCmtmFZtY3R1mc4A6WAFsLab+UmVnX8PN3z1ElfZDZGqFJxef+Pktww54fFLJxOcdnDsUeU0FxCxRvTA3bKsu4bYMxFRSf6Vo1nkLni808jtuSHzuVPITcfQPwbPjr+7NUOYVgSuctdz9gx4ZZ8Zxc7ZnZCIJ//F0Et4HvlMzsWOAZgusKX5F+7qKZnWNmP86juUcJvsnIZjx7LwAwL/89LXkDCD5/rkuoNU2r1wEvH6gxxee+wnOa/wH4ZTgWFKKc43M/xR5TwzbLPm6LPKZC+cZtUcdUUHw2KdJ4Cp0oNvM5bjvE2Ol53Nq6s/8AZ5H7duBzw7IZGa9/nODGGrOytHckuW8F/p2wva+29+duw/48EdhMcAmyWJbyGWTcev4A/enA73O816/C8kfb+3O3UV8ODz/fD7OUxQgGTwe+m0d/lnV8ZnzemeHnPfEA9RSfvk88+gHq5T2mRujnThe3efRn3mNqhP7sdHEbpT8LHVMj9Genis+osZmxTaTxNEJfdorYLOS4LfWxs907tdR+gK+GHfgQcAwwFrgrfO1HWeovDct25GhvJpAkmKo6CRhFcOmxFPAkUNHen7mN+vFEgimyFMEU3IIsP29lOWBy9ifB4igPD5yzwkHtRODH4etLgH7t/dnbqD+HNg3g4QB0CjAMOA34Tfj600BXxWdB/bsQeCFCvbKOT6AfMJDg29qmeBwY/mT9bPmOqQfq57C8U8RtPv1Z6JhaTnGbZ38WNKaWS3wWcqynbRtpPC2H2GzlcftVSnTsbPeOLcUf4GKC6aVtwE5gPvDxHHWvBXYAt7XQXtNgtJngOrqLw+0S7f1Z27APm4L+QD8rovYnMAT413BA3xAOLFsJFh5dm22Q70w/BN8afA14LoylxvDxaeAfyfh2QvEZuV9PDmPxygh1yzo+gRVRj+WM7SKPqQfq57Q6HT5u8+nPQsfUcorbfOOzkDG1XOKzFcd65PG0HGKzNcdtuH1Jjp0WNiQiIiIiItIiLZgWEREREZFIlDyIiIiIiEgkSh5ERERERCQSJQ8iIiIiIhKJkgcREREREYlEyYOIiIiIiESi5EFERERERCJR8iAiIiIiIpEoeRARERERkUiUPIiIyEFlZv9qZm5mZ7b3voiISH7M3dt7H0REpIyY2XPAaGCAuze29/6IiEh0mnkQEZGDxswGABOB3yhxEBHpeJQ8iIjIwXQRwf89j7T3joiISP6UPIiIyH7M7JvhuoSnspSZmf00LP+tmVXk0fQUoA544gDvf2rYftPPfWY23Mxmm9lWM9toZg+EMxmY2Wgze8zMtpvZ5rB+z7w+tIiIHJCSBxERyeZmYANwtpmdnVH2XeDDwLPApe7eEKVBM+sGnA38n7vvPED1+cAg4Ivh7/2A/wXuAE4AvgN8FHjCzA4FvgF8jeCUqF8AVwH3R9kvERGJTgumRUQkKzP7LPA9YIG7nxC+9jXgRmAh8D53355He1OBh4HPuPvdEbeZAfwo/HWCu7+cVvYn4DTgJeB8d98Qvh4DVgBDgWHuvjrqPoqISMs08yAiIrncA7wKHG9ml5nZFwgSh1eA8/JJHEJTAKew9Q7L0hOH0MLwcV1T4gDg7imChALg2ALeS0REcki09w6IiEhpcvdGM/sXYC5wF9CH4Bv9c9z93XzaCmcDLiSYxVhbwO68leW17S2UbQsfexXwXiIikoNmHkREJCd3fwRYBvQFNgJnu/vbBTR1MsG6hbkF7srmbLsXoSxe4PuJiEgWSh5ERCQnM/s8cFT4axf2ftufr4vDx0Iv0drSAj0t3hMROUiUPIiISFZmdhXwbeBt4FGgBvhKgc1NAd5y9yXF2TsREWkPSh5ERGQ/ZjYN+CHBKUHnAFcDtcCnzOzIPNsaAxxJ4acsiYhIiVDyICIi+wjv6/BzYDfBVZVeCS93+j8EF9r4zzybbO0pSyIiUiKUPIiISDMzOwmYE/56sbsvSCu+meAqRtPM7NQ8mp0CbCG4qVzU/Yib2UCg6S7RXc1soJl1TSurDsuqw7LKsHwg0DUs6xmWaeG0iEgR6CZxIiICgJmNB/4I9AAuc/f9TjMys1nAN4H57n5ShDb7A+uAn7n7x/LYl+FkvwTrx4FncpSdBQxn703l0o1w9xVR319ERLJT8iAiIm3GzP4B+F9gurv/qr33R0REWkenLYmISFuaAtQDT7T3joiISOvpDtMiItKW/gI85u472ntHRESk9XTakoiIiIiIRKLTlkREREREJBIlDyIiIiIiEomSBxERERERiUTJg4iIiIiIRKLkQUREREREIlHyICIiIiIikSh5EBERERGRSJQ8iIiIiIhIJP8PlhlaT45BE4AAAAAASUVORK5CYII=",
- "text/plain": [
- "