{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "*This notebook contains material from [cbe67701-uncertainty-quantification](https://ndcbe.github.io/cbe67701-uncertainty-quantification);\n", "content is available [on Github](https://github.com/ndcbe/cbe67701-uncertainty-quantification.git).*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "< [4.1 Local sensitivity analysis and difference approximation](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.01-Local-sensitivity-analysis.html) | [Contents](toc.html) | [5.0 Regression Approximations to Estimate Sensitivities](https://ndcbe.github.io/cbe67701-uncertainty-quantification/05.00-Regression-Approximations-to-Estimate-Sensitivities.html)

\"Open

\"Download\"" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "RUPnOS_kKsBa", "nbpages": { "level": 1, "link": "[4.2 Advection-Diffusion-Reaction (ADR) Example](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.02-Simple-ADR-Example.html#4.2-Advection-Diffusion-Reaction-(ADR)-Example)", "section": "4.2 Advection-Diffusion-Reaction (ADR) Example" } }, "source": [ "# 4.2 Advection-Diffusion-Reaction (ADR) Example\n", "Noah Wamble (nwamble@nd.edu) \n", "June 22nd, 2020" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 1, "link": "[4.2 Advection-Diffusion-Reaction (ADR) Example](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.02-Simple-ADR-Example.html#4.2-Advection-Diffusion-Reaction-(ADR)-Example)", "section": "4.2 Advection-Diffusion-Reaction (ADR) Example" } }, "source": [ "This example was adapted from code from:\n", "\n", "McClarren, Ryan G (2018). *Uncertainty Quantification and Predictive Computational Science: A Foundation for Physical Scientists and Engineers*, *Chapter 4: Local Sensitivity Analysis Based on Derivative Approximations*, Springer, https://doi.org/10.1007/978-3-319-99525-0_4" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "pm6095qCKz8I", "nbpages": { "level": 3, "link": "[4.2.1 **What do we want to know?** ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.02-Simple-ADR-Example.html#4.2.1-**What-do-we-want-to-know?**)", "section": "4.2.1 **What do we want to know?** " } }, "source": [ "### 4.2.1 **What do we want to know?** \n", "\n", "* How would the QoIs vary for small, expected perturbations in the inputs if given knowledge of the QoI at a particular value? To accomplish this, we are interested in the derivative of the QoI at a nominal input value (95)." ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "AUwChLScLAYK", "nbpages": { "level": 3, "link": "[4.2.2 **Why do we want to know that?**](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.02-Simple-ADR-Example.html#4.2.2-**Why-do-we-want-to-know-that?**)", "section": "4.2.2 **Why do we want to know that?**" } }, "source": [ "### 4.2.2 **Why do we want to know that?**\n", "\n", "\n", "* This information will allow us **to see which parameters have a larger effect** on the QoI. \n", "\n", "* Knowing the sensitivities allows the researcher **to narrow the focus** to a smaller number of parameters and then apply more rigorous techniques to only these parameters. \n", "\n", "* If the variabilities in these distributions are small, it may be possible to **quantify the uncertainty** in the QoI to these parameters using only local information (95).\n" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "GAePL8DTL-Vb", "nbpages": { "level": 3, "link": "[4.2.3 **4.3.1 Simple ADR Example** ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.02-Simple-ADR-Example.html#4.2.3-**4.3.1-Simple-ADR-Example**)", "section": "4.2.3 **4.3.1 Simple ADR Example** " } }, "source": [ "### 4.2.3 **4.3.1 Simple ADR Example** \n", "We will study the steady advection-diffusion-reaction equation in one-spatial dimension with a spatially constant, but uncertain, diffusion coefficient, a linear reaction term, and a prescribed uncertain source.\n", "\n", "$$\n", "\\nu \\frac{du}{dx} - \\omega \\frac{d^2u}{dx^2} + \\kappa(x)u = S(x)\n", "$$\n", "With initial condition: \n", "$$\n", " u(0)=u(10)=0\n", "$$\n", "Where:\n", "$$\n", "\\kappa(x) = \\begin{cases}\n", " \\kappa_l & x \\: \\epsilon \\: (5,7.5)\\\\\n", " \\kappa_h & \\text{otherwise}\\\\\n", " \\end{cases} \n", "$$\n", "$$\n", "$$\n", "$$ \n", "S(x)= qx(10-x)\n", "$$\n", "\n", "With $\\mu_v =10$, $\\mu_\\omega= 20$, Var$(\\mu_v)=0.0723493$, Var$(\\omega)= 0.3195214$\n", "\n", "\n", "\n", " $\\mu_{\\kappa_h} = 2$, Var$( \\kappa_h ) = 0.002778142$ , and $\\mu_{k_l} = 0.1$, Var$( \\kappa_l )= 8.511570 \\times 10^{-6}$, $\\mu_q = 1,$ Var$(q)= 7.062353 \\times 10^{-4}$ \n", "\n", "We also prescribe that the parameters ordered as $(\\nu,\\omega,\\kappa_l,\\kappa_u,q)$ have the following correlation matrix: \n", "\n", "$$\n", "R= \n", "\\begin{pmatrix}\n", "1.00 & 0.10 & -0.05 & 0.00 & 0.00\\\\\n", "0.10 & 1.00 & -0.40 & 0.30 & 0.50 \\\\\n", "0.00 & 0.30 & 0.20 & 1.00 & -0.10 \\\\\n", "0.00 & 0.50 & 0.00 & -0.10 & 1.00 \\\\\n", "\\end{pmatrix}\n", "$$\n", "The QoI is total reaction rate: \n", "$$\n", "Q= \\int^{10}_{0} \\kappa(x_ u) u(x)dx\n", "$$\n" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "dwRcQkuIQozE", "nbpages": { "level": 3, "link": "[4.2.3 **4.3.1 Simple ADR Example** ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.02-Simple-ADR-Example.html#4.2.3-**4.3.1-Simple-ADR-Example**)", "section": "4.2.3 **4.3.1 Simple ADR Example** " } }, "source": [ "**Import Libraries**\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "colab": {}, "colab_type": "code", "id": "uGpgzLB_PCOh", "nbpages": { "level": 3, "link": "[4.2.3 **4.3.1 Simple ADR Example** ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.02-Simple-ADR-Example.html#4.2.3-**4.3.1-Simple-ADR-Example**)", "section": "4.2.3 **4.3.1 Simple ADR Example** " } }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import scipy.sparse as sparse\n", "import scipy.sparse.linalg as linalg\n", "import scipy.integrate as integrate\n", "import math" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "ms4kt8rRWl9p", "nbpages": { "level": 3, "link": "[4.2.3 **4.3.1 Simple ADR Example** ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.02-Simple-ADR-Example.html#4.2.3-**4.3.1-Simple-ADR-Example**)", "section": "4.2.3 **4.3.1 Simple ADR Example** " } }, "source": [ "**Numerical Method to Solve ADR Equation** " ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "colab": {}, "colab_type": "code", "id": "TLf2mdKkKy_z", "nbpages": { "level": 3, "link": "[4.2.3 **4.3.1 Simple ADR Example** ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.02-Simple-ADR-Example.html#4.2.3-**4.3.1-Simple-ADR-Example**)", "section": "4.2.3 **4.3.1 Simple ADR Example** " } }, "outputs": [], "source": [ "def ADRSource(Lx, Nx, Source, omega, v, kappa):\n", " \"\"\"\n", " This function solves the diffusion equation with a generalized source. \n", " \n", " Inputs \n", " Lx, Nx: Used to make the size of the dx term \n", " Source: Source function\n", " omega: diffusion coefficient \n", " v: velocity field \n", " kappa: reaction coefficient equation \n", "\n", "\n", " Outputs:\n", " Solution: u(x) values \n", " Q: Total reaction rate \n", " \n", " \"\"\"\n", " \n", " #Initialize A matrix that will be solved \n", " A = sparse.dia_matrix((Nx,Nx),dtype=\"complex\")\n", " #Initialize size of dx step \n", " dx = Lx/Nx\n", " i2dx2 = 1.0/(dx*dx)\n", " #fill diagonal of A\n", " A.setdiag(2*i2dx2*omega + np.sign(v)*v/dx + kappa)\n", " #fill off diagonals of A\n", " A.setdiag(-i2dx2*omega[1:Nx] +\n", " 0.5*(1-np.sign(v[1:Nx]))*v[1:Nx]/dx,1)\n", " A.setdiag(-i2dx2*omega[0:(Nx-1)] -\n", " 0.5*(np.sign(v[0:(Nx-1)])+1)*v[0:(Nx-1)]/dx,-1)\n", " \n", " #solving equation A x = Source\n", " Solution = linalg.spsolve(A,Source)\n", " #Use trapezoidal integrater to find integral\n", " Q = integrate.trapz(Solution*kappa,dx=dx)\n", "\n", " return Solution, Q" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "kCjKFR0FYHQx", "nbpages": { "level": 3, "link": "[4.2.3 **4.3.1 Simple ADR Example** ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.02-Simple-ADR-Example.html#4.2.3-**4.3.1-Simple-ADR-Example**)", "section": "4.2.3 **4.3.1 Simple ADR Example** " } }, "source": [ "**Define Parameters**" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 385 }, "colab_type": "code", "id": "SKObIHN3KRhC", "nbpages": { "level": 3, "link": "[4.2.3 **4.3.1 Simple ADR Example** ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.02-Simple-ADR-Example.html#4.2.3-**4.3.1-Simple-ADR-Example**)", "section": "4.2.3 **4.3.1 Simple ADR Example** " }, "outputId": "c075e88b-3977-4adf-b684-a568bc0002dd" }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/anaconda3/lib/python3.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.py:133: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format\n", " SparseEfficiencyWarning)\n", "/anaconda3/lib/python3.7/site-packages/numpy/core/numeric.py:538: ComplexWarning: Casting complex values to real discards the imaginary part\n", " return array(a, dtype, copy=False, order=order)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaoAAAEWCAYAAAA3h9P4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd8VvX5//HXlcFI2BDCDCFM2SOAoKLiqIp1L5ClKNraaq39Vm1/tbbVVjusraMOtojbqlXcG0VIGLL3nknYSci+fn+cQxtjAgnkzufcua/n48GDO/c67zOv+3zOOZ8jqooxxhgTVFGuAxhjjDHHYoXKGGNMoFmhMsYYE2hWqIwxxgSaFSpjjDGBZoXKGGNMoJ1UoRKR+0VkVnWFOcEMKiKdT/CzZ4jImurOVMlh9xCR9Eq+93URuSDUmcoZ7lkisr2mh1seEZkgInNd5wgXruadiCSKyBciclhE/lbTww8Vl9sKc5xCJSLZpf6ViMiRUn9fX1Mhq0vZoqaqX6pqN0dx/gD8tZLvfQh4MIRZQk5EPhORmxwNe7qIPOBi2BFoEpAFNFLVu1yHqUhVf2Q73lY4E5R155iFSlUbHP0HbAV+WOq552smYu0jIq2Bs4E3KvN+VV0ANBKR1JAGM+bkdQBWaoB7EhCRGNcZTlS4Za+2vKpaqX/AZuDcMs/dD7wMzAQOAyuA1FKvtwFeAzKBTcDtx/j+i4CV/vfsAH5R6rWbgfXAPuAtoE2p1xTo7D/+DLip1GsTgLn+4y/89+YA2cC1wFnA9lLvP8X/jgP+uFxS6rXpwBPAO37G+UCnCsblO99bdvoB44CPSr3WyR+3AaWmWxZwVqn3PAv89hjT72JgiZ/9a6CP//w9wKtl3vsP4J/+4xuAVf44bQRuqWg8Sk/rUtPkAf9xU+Btf17v9x+38197ECgG8vxp/7j/fHfgQ3/c1wDXlPru5v68PgQswNsDnXuM8X8F2A0c9Od1T//5SUAhUOAP+z8VfF6BHwPr/GnxB3++zPMzvAzUOd70LjXNN/jfsxK4vOwyibc3vR9vvbiwgkw1Mu+ONz7l5BoGpPnTOg0YVuo7S0/rc8v57GdUsI6WynmrPx/2461zUmZbcHScV/LddabcbQ3edupVYJY/L3/iZyz0c357AtNzM/ALYKk/HV4C6lUwvSYAXwGP+e9dDZxT6vXjDhe4G2/5fo5jrGulpvED/nzMBv6Dtz49749/GpBc6v3lrodUsO5UcVrfBAwG0v2/9wCPVLRsVbjMVfqNFReqPLwiEw38CfjGfy0KWAjcB9QBUvyZ8IMKvn8XcEapjd7RBXAE3kZ7AFDXn9lflLcCUrmVoPTK+t+FD4jFK4a/8vOO8BecbqVWwn3+RI/xZ/qLFYzLf7+3vOkH/AV4oszrR1fAOOB94K9lXv858HoFwxsAZABD/Pkw3h9eXbxfuLl4TTH4r+8CTvX/Hom3QRbgTP+9A8obj3Km33T+V6iaA1f6+RviFY43jrGBige24a2kMf44ZPG/AvMiXnGIB3rh/Xg5VqG60R9uXeBRYEl5OY/xecUrjI2AnkA+8DHectsYb6M4/njT23/9aryVOQrvB1EO0LrUMlnoz+9o4EfATkptjEtlqql5d8zxKZOpGd7Gcaw/30b5fzevzLQuZzmYwPfX0beBJkAS3sbwglLTdQcwyB/nzv40Oua2Bm87VQhc5r+3vv/crDLZqjI9N+P9gGrjT5NVwK0VjPMEoAi4E287cy1ewWpWyeEWAQ/jLdv1qdy6tt7/zqPL7lrgXH+ezQSmVXI9/M78PMFpPQ8Y67/eAH/5rcq/6jjrb66qzlHVYrxq39d/fhCQoKq/V9UCVd2It1dwXQXfUwj0EJFGqrpfVRf5z18PTFXVRaqaD9wLDBWR5GrIXtqpeBPxIT/vJ3grzKhS73ldVReoahFeoep3gsNqglcE/0tVn8X7FTkfaA38usxnDvufK8/NwNOqOl9Vi1V1Bt6G9lRV3QIswltwwCvAuar6jT/cd1R1g3o+Bz4AzqjqCKnqXlV9TVVzVfUw3l7Umcf4yMXAZlWdpqpF/vx+DbhKRKLxVsT7VDVHVZcDM44z/KmqethfRu4H+opI4yqOxsOqekhVVwDLgQ9UdaOqHgTeBfr776twevtZXlHVnapaoqov4c3XwaWGs0VVn/XXmRl48zuxnHGqkXl3vPEpYySwTlWf8+fbC3h7CD88geFW5CFVPaCqW4FP+d96dhPwZ1VN88d5vT+NKrOtmaeqb/jz5Eh5Az2B6flPfz7vw9trOdb2IAN4VFUL/WViDd60rMxwS/BaU/JV9Ugl17Vp/nceXXY3qOpH/rbrFf63LFe4HlYwHicyrQuBziLSQlWzjy6/VVEdhWp3qce5QD2/XbID0EZEDhz9h7e38r0V0ncl3p7ZFhH5XESG+s+3AbYcfZOqZgN7gbbVkL20NsA2VS0p9dyWMsMpO64NTnBY+/F+CZX1LN7ew2P+Bre0hnjNMuXpANxVZlq3xxsngNn8r+CO9v8GQEQuFJFvRGSf/7mLgBZVHSERiRORp0Vki4gcwmt+a+IXnYoyDymT+XqgFZCA9+tuW6n3bynnO44OO1pEHhKRDf6wN/svVXU89pR6fKScv4/O72NObxEZJyJLSr3Wq0yW/y5HqprrP6xoWQr5vDve+JTxnfXRV3Y9OVkVrWft8ZpUy6rMtmZbOZ/7jhOYnlXZHuxQ9XYpfFv43/JyvOFmqmpeqZyVWdeqsixXtB6W50Sm9USgK7BaRNJE5OIKvrtCoTwwtw3YpKpdKvNmVU0DLhWRWLw25JfxFsydeBMHABGJx9v13VHO1+Tg7Q4fVdHELs9OoL2IRJUqVkl4u8xV9Z0c/gKUUOr1pXjNK5R6TwO8JqspwP0i8pr/S+2oU4BvKxjeNuBBVa3ozMBXgL+JSDvgcmCoP8y6eL+exgFvqmqhiLyB1wRRnly+P32PngJ9F9ANGKKqu0WkH7C41HeVXkmPZv5cVc8rOxB/ehXhzf/V/tNJFWQCbwN+KV7Txma85o79xxj2yapweotIB7wfHOfg/bIsFpElVDxNj6cm5t3xlp/SvrM++pKA9yrxWTi5dXQbXnNWec8fb1tTdhn4zt8nMD2rqq2ISKlilQS8Vcnhls1+vHWtKipcDysYdpWntaquA0aJSBRwBfCqiDRX1ZzKhgzlBb8LgEMicreI1Pd/9fYSkUFl3ygidUTkehFprKqFeAfdiv2XZwM3iEg/f6b+EZivqpvLGeYS4Ar/F0dnvEpe2h68NtXyzMdbiX4pIrEichZec8aLVRprz1q8PcuRfuH9f3jty0d9CAwQkXqlnvsHsFBVb8I7YeOpMt95Jt4ufHmeBW4VkSHiifeH3RBAVTPx2q2n4S1kq/zP1fFzZQJFInIhcP4xxmsJMNqflxfw3eaGhni/1A6ISDPgt2U+W3bavw10FZGx/vSOFZFBInKK3yT2Ol7BjhORHpQp7GU0xGuq2ou3EfzjcYZ9so41vePxVtRMABG5AW+P6oTU0Lw75vJTxhy8+TZaRGJE5FqgB978rIzjraPHMhn4hYgM9HN29n8YVHpbU8oeINnfeELVp2dVtQRu95fzq/F+eM45weEeb12rigrXQ//1sutOlae1iIwRkQR/B+Boq1BxRe8vT8gKlb+x+SFeu+0mvAN0k/F+7ZZnLLDZ35W9FRjjf8/HwG/wfnXswvtFVdFxrr/jnaGyB6/tv+wp9PcDM/xd1mvK5C0ALgEu9LM+CYxT1dVUkd8u/GO88d2BVwC3l3p9D/AJ3l4AInIpcIE/3uCdODFA/GvV/IUgR73T1MsbXjrecYbH8fYk1uMdwC1tNt4ex+xSnzsM3I6397ofb8/krWOM2h148/Ro80Dp0+sfxTtwmgV8w/d/Yf8D7/jTfhH5pz/s8/Hm5U68ZpSjB4zB26tu4D8/HW9DXZGZeE0pO/AOHJdtA5+Cd/zzgP9r9aQca3qr6krgb3gHkPcAvfHO+DoZIZ13lVx+jr53L95xjbvwfhj8ErhYVbMqOS7HW0crpKqv4B2PmY13zPYNvBMSqrqtAW9PFWCviCw6gelZVfOBLn62B4Gr/GNNJzLc461rlVaJ9fA7684JTusLgBUiko23HbiudFNmZch3m01NTfH3EmYAg/U4M0FEXgOmqOqcGglnjKk2IjIB70zH011nCVdhdfFYbeL/8j5W00Tp914Z4jjGGBNY1imtMcaYQLOmP2OMMYFme1TGGGMCLWyPUbVo0UKTk5NdxzDGmLCxcOHCLFVNOP47gyVsC1VycjLp6ZW6nZMxxhhARCrs4SXIrOnPGGNMoFmhMsYYE2hWqIwxxgSaFSpjjDGBZoXKGGNMoFmhMsYYE2hWqIwxxgRa2F5HZYw5OXuz89mQmUPm4Xz25eRzKK+Io12q1YuNpkWDurRoUJeUhHhaN66HSHXdQ9CYqrFCZUwEKCgqYcm2A8zfuJcFm/exYuch9uUUVPrzDevG0KNNI4Z2as6wTi0YkNSEmGhrkDE1wwqVMbVUXmExX6zN5N3lu/lo5R4O5xcB0L1VQ847JZEuiQ3o3LIBrRrXo1l8HRrXjyVKBFU4UlBMVk4+GYfyWZ+Zzbo9h1my7QD//Hgdj360jmbxdRjZuzWX9W/DgKSmtrdlQipse09PTU1V60LJmO/bnJXD7AVbeSV9G/tzC2kSF8v5PRI555REBic3o2l8nRP+7oO5hczbmMXbS3fx0ao95BWW0KttI246PYWRfVoTa3tZgSYiC1U11XWOqrJCZUwtsXjrfh7/ZD0fr84gOko4v0ci1w1OYlin5iEpIDn5RbyxZAdT525iQ2YOSc3iuOv8rvywTxuiomwPK4isUNUwK1TGeBZt3c+jH63ji7WZNImLZfzQZEYPSSKxUb0aGX5JifLpmgz++sFaVu06RM82jfjDZb0YkNS0RoZvKs8KVQ2zQmUi3fb9uTz07mreXrqL5vF1uHl4CmNO7UCDum4OPZeUKG99u5OH31vN7kN5jB+azC9+0M1ZHvN94VqoArEEiUg34KVST6UA96nqo44iGRNYeYXFPPnpep7+YiMAt5/ThVvPTCGujtvVOSpKuKx/W87tkchf3lvNjHmb+WjVHh4b1Z/+tndlTkLg9qhEJBrYAQxR1QrvnWJ7VCYSpW3ex92vLmVjVg6X9G3D3Rd2p22T+q5jlSt98z7ueHEJew7lcfcF3Zl4ekc7duWY7VFVn3OADccqUsZEmpz8Iv783mpmfrOFtk3q89zEwZzRJdg3ak1Nbsac28/gl699y4NzVrFk2wH+enVf6teJdh3NhJkgFqrrgBfKe0FEJgGTAJKSkmoykzHOLN9xkNtfWMymvTmMH5rM//2gG/FhctyncVwsT40ZyLNfbuRP765m+/5cnh2XSssaOtHD1A6BavoTkTrATqCnqu451nut6c/UdqrKtK8289C7q2kaH8uj1/ZnaKfmrmOdsA9W7OaOF5fQNC6WWTcNISWhgetIESdcm/6CdnXehcCi4xUpY2q7Q3mF3DxzIb9/eyXDu7bg3TuGh3WRAji/ZyteuXUo+UUlXPP0N6zefch1JBMmglaoRlFBs58xkWJDZjaXPfEVn63J4L6Le/DsuFSanURvEkHSq21jXrplKNFRcN0z37Bs+0HXkUwYCEyhEpE44DzgdddZjHHlk9V7uOzxrziYW8jzNw3hxtM71rp+9Dq3bMArtwyjQd0Yxk6dz9o9h11HMgEXmEKlqrmq2lxV7SeWiTiqytOfb2DijHQ6tIjjrZ+ezpCU8G7qO5ak5nHMvulU6sZEMWbyfLbuzXUdyQRYYAqVMZGquES5/60V/Ond1VzUuzWv3DIssNdGVaek5nHMmjiEwuISRk/+ht0H81xHMgFlhcoYh/IKi/nx8wuZMW8LN5/Rkceu6x9R1xl1SWzIjBsHsz+ngAnTFnA4r9B1JBNAVqiMceRAbgFjJs/ng5V7+M3FPfj1yB4R2XNDn3ZNeHLMQNZlZHPb7MUUFpe4jmQCxgqVMQ7szc5n1LPzWbr9II+PGsDE0zu6juTUmV0TeOCyXnyxNpP73lxOkK7vNO6Fx+XtxtQiGYfyGD15Ptv35zJ5fCrDuwa7K6SaMmpwEtv25fLkZxtIahbPj87q5DqSCQgrVMbUoJ0HjjD62W/IOJzP9BsGc2otPrPvRPzi/G5s23+Eh99bTUpCPD/o2cp1JBMA1vRnTA3Zti+Xa56ex97sAp6baEWqPFFRwl+u6kPfdo35+UtLWJ9h11gZK1TG1IhdB48w6tlvOJxXxPM3D2Fgh2auIwVWvdhonho7kPp1opk0cyGH7EzAiGeFypgQyzycz/XPzudgbiHPTRxMn3ZNXEcKvNaN6/PE6AFs3ZfLnS8uoaTETq6IZFaojAmh/TkFjJ0yn10H85h2wyArUlUwJKU5v7m4Bx+vzuDRj9a6jmMcskJlTIgcyitk/LQFbMzKYfL4VFKTrbmvqsYN7cCVA9rx2Kfr+XJdpus4xhErVMaEQF5hMTfNSGflzkM8NWYAp3Vu4TpSWBIR/nBZTzonNODOl5aQcdi6WYpEVqiMqWYlJcpdL3/Lgk37eOTafozonug6UliLqxPDE9cPIDu/iDtfWkKxHa+KOFaojKlmD85ZxTvLdvHri07hkr5tXMepFbomNuR3l/Tkq/V7+ddn613HMTXMCpUx1WjylxuZMncTE4Ylc9MZkd0tUnW7JrU9l/RtwyMfriVt8z7XcUwNskJlTDV5e+lOHnhnFRf2asVvLu5R62546JqI8ODlvWjXNI6fv7yE7Pwi15FMDbFCZUw1mL9xLz9/6VtSOzTl79f2IzoCe0GvCQ3rxfLINX3Zvv8ID76zynUcU0OsUBlzktbtOczNM9Np36w+k8enUi82cu4n5UJqcjMmDU/hhQVb+XR1hus4pgYEplCJSBMReVVEVovIKhEZ6jqTMcez51AeE6alUTc2muk3DKZJXB3XkSLCz8/rSrfEhvzytaXszylwHceEWGAKFfAP4D1V7Q70BWy/3gTa4bxCJkxL40BuAdMmDKJ9szjXkSJG3ZhoHrm2LwdyC/h/by53HceEWCAKlYg0AoYDUwBUtUBVD7hNZUzFCopK+PHzi1i75zBPjhlIr7aNXUeKOD3bNOZn53blnaW7eG/5LtdxTAgFolABKUAmME1EFovIZBGJL/smEZkkIukikp6Zad2pGDdUlXteX8qX67J46IrenGk3PnTmluEp9GjdiPveXMHBI9bLem0VlEIVAwwA/qWq/YEc4J6yb1LVZ1Q1VVVTExJs42DceOTDtby+aAc/P68rV6e2dx0nosVER/HwlX3Iys7n4fdWu45jQiQohWo7sF1V5/t/v4pXuIwJlOfnb+GxT9YzanB7fjqis+s4BujdrjETT+/I7Plbmb9xr+s4JgQCUahUdTewTUS6+U+dA6x0GMmY7/l41R5+88Zyzu6WwB8u7WUX9AbIned1pV3T+tz772XkFRa7jmOqWSAKle+nwPMishToB/zRcR5j/mvJtgP8ZPZierVtzOOjBxATHaRVx8TVieGPl/dmY2YOT3xqfQHWNoFZ21R1iX/8qY+qXqaq+11nMgZgc1YOE6enkdCwLlPGDyK+bozrSKYcw7smcHn/tjz9+UY2ZeW4jmOqUWAKlTFBtDc7nwnTFlCiyvQbBpHQsK7rSOYY7r2oO3Viovjdf1agarcDqS2sUBlTgSMFxUyckc6ug3lMHj+IlIQGriOZ42jZsB4/O7cLn63J5KNV1r1SbWGFyphyFBWX8NMXFrF0+wEeG9WfgR2auo5kKmn8sGS6Jjbg92+vsBMragkrVMaUoarc/58VfLQqg/sv6cn5PVu5jmSqIDY6ivsv6cm2fUd46vMNruOYamCFypgynvxsA7O+2cqtZ3Zi3NBk13HMCRjWqQUX92nNvz7bwLZ9ua7jmJNkhcqYUt5YvIO/vL+Gy/q14Zc/6Hb8D5jA+vXIU4gS4SHrsSLsWaEyxjdvw17+79VvOTWlGQ9f1Ycou/lhWGvduD6ThqfwztJdLNxiV7uEMytUxgDrMw5zy3PpdGgez9NjUqkbYzc/rA0mDU8hoWFd/jhnlZ2uHsasUJmIl3k4nwnT0qgTE820CYNoHBfrOpKpJvF1Y7jrvK4s3LKf95bvdh3HnCArVCai5RYUMXFGGnuzC5g6IdVuflgLXZ3anm6JDXnovdUUFJW4jmNOgBUqE7GKS5TbX1jM8h0HeWxUf/q0a+I6kgmB6Cjh3ou6s2VvLrO+2eI6jjkBVqhMRFJVfl/qWqlzeyS6jmRC6MyuCZzRpQX//GSd3WAxDFmhMhFpytxNzJi3hZvP6GjXSkUAEeHuC7pzILeQKV9udB3HVJEVKhNx3l22iwfnrOLCXq2498JTXMcxNaRX28Zc1LsVU+ZuYl9Oges4pgqsUJmIsnDLfn720hL6t2/C36/tZ9dKRZifn9eVI4XF1rVSmLFCZSLG5qwcbp6ZTqvG9Xh2XCr1Yu1aqUjTuWVDLuvflhlfb2bPoTzXcUwlWaEyEWF/TgE3TE9DVZl+w2CaN7D7SkWqn53TleIS5fFP7E7A4cIKlan18gqLuXlmOjsOHOHZcal0bBHvOpJxKKl5HNcMas+LaVutw9owEZhCJSKbRWSZiCwRkXTXeUztUFKi3PXKt6Rv2c/fr+lHanIz15FMAPx0RGdEhMc+Wec6iqmEwBQq39mq2k9VU10HMbXDw++v5p2lu7j3wu6M7NPadRwTEK0b12f04CReX7TD9qrCQNAKlTHVZtY3W3j6841cPySJScNTXMcxAXPLmSlEidgZgGEgSIVKgQ9EZKGITCrvDSIySUTSRSQ9MzOzhuOZcPLp6gzue3M5I7q35HeX9ETETkM339W6cX2uSm3HK+nb2X3QzgAMsiAVqtNUdQBwIXCbiAwv+wZVfUZVU1U1NSEhoeYTmrCwfMdBbpu9iB5tGvHYqP7ERAdpMTdB8qMzO1GsyjNfWG8VQRaYNVhVd/r/ZwD/Bga7TWTC0a6DR7hxehpN6scydfwg4uvGuI5kAqx9szgu69eW2Qu2kJWd7zqOqUAgCpWIxItIw6OPgfOB5W5TmXCTk1/ExOnp5BYUM/WGQbRsVM91JBMGfnx2J/KLSpgyd5PrKKYCgShUQCIwV0S+BRYA76jqe44zmTBSUqL87KUlrN59iMdG9ad7q0auI5kw0SmhASN7t2bm15s5kGt9AAZRIAqVqm5U1b7+v56q+qDrTCa8PPz+aj5cuYffXNyDs7u3dB3HhJmfjOhMTkExM+fZ/aqCKBCFypiT8XLaNp7+fCNjTk1iwrBk13FMGOreqhFndUtgxtebySssdh3HlGGFyoS1eRv28qt/L+OMLi24/4d2Gro5cZOGp7A3p4B/L97hOoopwwqVCVubsnK4ddZCklvE8/joAXYaujkpQ1Oa06ttI579ciMlJeo6jinF1mwTlg7kFjBxehpRAlPHD6Jx/VjXkUyYExFuPiOFjZk5fLw6w3UcU4oVKhN2CotL+NGsRWzff4RnxqWS1DzOdSRTS1zUuzVtm9TnWbsAOFCsUJmwoqr85o3lzNu4lz9d0ZtB1hu6qUax0VHceHpHFmzex+Kt+13HMT4rVCasTJm7iRfTtnHb2Z24cmA713FMLXTtoPY0rBfD5C/tAuCgsEJlwsZHK/fw4JxVXNirFXed1811HFNLNagbw/VDOvDu8l12C5CAsEJlwsLKnYe4/cXF9GrTmEeu6UdUlJ2GbkJn3NAOiAizvrELgIPACpUJvIxDedw0I43G9WOZPD6V+nWiXUcytVybJvU5v0ciL6Zt40iBXQDsmhUqE2hHCoq5eWY6B44UMnl8KonW0aypIeOHJXPwSCFvLrELgF2zQmUCq6REueuVJSzdcZBHr+1HzzaNXUcyEWRIx2Z0b9WQ6V9vRtUuAHbJCpUJrL9/tJY5y3Zz74XdOb9nK9dxTIQREcYPS2b17sMs2LTPdZyIZoXKBNK/F2/nsU/Wc21qe24+I8V1HBOhLuvXlsb1Y5kxb7PrKBHNCpUJnPTN+7j71WWcmtKMP1zWyzqaNc7UrxPNdYPa8/6KPew8cMR1nIhlhcoEyta9uUx6biFtm9bnqTEDqRNji6hxa8ypHVBVnp9vp6q7YlsBExiH8gqZOCON4hJlyvhUmsTVcR3JGNo3i2NE95a8lLadwuIS13EikhUqEwhFxSX8ZPZiNmXl8K8xA0hJaOA6kjH/NXpIElnZ+Xy0co/rKBEpUIVKRKJFZLGIvO06i6lZf3h7JV+szeSBy3oxrFML13GM+Y4zu7akTeN6zF6w1XWUiBSoQgXcAaxyHcLUrBlfb2bGvC1MGp7CdYOTXMcx5nuio4RrByXx5bostu61/v9qWmAKlYi0A0YCk11nMTXnszUZ/O4/Kzj3lETuvqC76zjGVOiaQe2IEngxzfaqalpgChXwKPBLoMKjlSIySUTSRSQ9MzOz5pKZkFiz+zA/mb2Ybq0a8Y/r+hFtHc2aAGvduD4juifycrqdVFHTAlGoRORiIENVFx7rfar6jKqmqmpqQkJCDaUzoZCVnc/EGWnUrxPNlPGpxNeNcR3JmOMaPaS9nVThQMgKlYg0FZGeIpIiIscbzmnAJSKyGXgRGCEis0KVzbiVV1jMLc8tJPNwPpPHpdKmSX3XkYypFDupwo1qLVQi0lhEfiUiy4BvgKeBl4EtIvKKiJxd3udU9V5VbaeqycB1wCeqOqY6s5lgUFXueW0pC7fs55Fr+tG3fRPXkYypNDupwo3q3qN6FdgGnKGq3VT1dL+prj3wEHCpiEys5mGaMPL4J+t5Y8lOfnF+V0b2ae06jjFVdvSkilcXbnMdJWJU64EBVT3vGK8tBI55DMp/32fAZ9WXygTF20t38rcP13JF/7bcdnZn13GMOSGtG9fn9C4JvLZoBz87t6vdbboGhOQYVdm9Jv9C3t+GYlgmPCzZdoC7Xv6W1A5N+dOVva2jWRPWrhrYjh0HjjBv417XUSJCqE6mOEdE5ohIaxHphXe8qmGIhmUCbseBI9w0I52Wjery9NiB1I2xW8mb8HZ+j0Qa1otdXCwnAAAWi0lEQVTh1YXbXUeJCCEpVKo6GpgBLAPmAD9T1V+EYlgm2LLzi5g4PY38wmKmjh9E8wZ1XUcy5qTVi43mh33b8O7yXRzOK3Qdp9YLVdNfF7zukF4DNgNjRSQuFMMywVVcotzxwmLWZWTz+PUD6JJoO9Wm9rhqYDvyCkuYs2yX6yi1Xqia/v4D/EZVbwHOBNYBaSEalgmoP81ZxcerM7j/hz04s6tdoG1ql/7tm9ApId6a/2pAqArVYFX9GEA9fwMuC9GwTADNnr+VyXM3MWFYMmOHJruOY0y1ExGuGtietM372ZyV4zpOrVbdF/yeDqCqh8q+pqrrRKSRf3KFqcW+Wp/FfW8u58yuCfy/kae4jmNMyFzevy1RAq8tsr2qUKruPaorReRrEblPREaKyGARGS4iN4rIc8DbgPWXU4ttyMzmR7MWkpIQz2Oj+xMTHYjuJI0JiVaN63FGlwReW7idkhJ1HafWqtatiKreiXerjl3A1cDvgTuBzsBTqjpcVe1YVS21P6eAidPTiI2OYsr4QTSqF+s6kjEhd8WAtuw8mMeCzftcR6m1qv3nrqruB9oAm4B5wBIgDzinuodlgqOgqIRbZy1k58E8nhk3kPbN7CRPExnO65FIXJ1o3lyy03WUWitU7TLZpf4VARcCySEalnFMVfn1v5cxf9M+/nJVHwZ2aOY6kjE1Jq5ODOf3SGTOsl0UFNl9qkIhVBf8/q3UvweBs4C2oRiWce/pLzbyysLt3H5OFy7tZ7PZRJ5L+7Xl4JFCPl9rN3QNhZo60h0HpNTQsEwNem/5bh5+bzUX92nNned2cR3HGCdO79KCZvF1eHPJDtdRaqWQ3FbVvx/V0VNgooEEvBMrTC2yfMdB7nxpCX3aNeGvV/e1jmZNxIqNjmJk79a8snAb2flFNLA7VlerUE3Ni0s9LgL2qGpRiIZlHNh9MI+JM9JoGhfLs+MGUi/WOpo1ke3Sfm147pstfLBiN1cMaOc6Tq0SqmNUW0r922FFqnbJLSjipplpZOcVMWXCIFo2rOc6kjHODUhqStsm9XnDzv6rdnY1pqmSkhLl5y99y8qdh3hsdH9Oad3IdSRjAiEqSri0Xxu+Wp9F5uF813FqFStUpkr++sEa3luxm1+P7MGI7omu4xgTKJf2a0txiVqP6tUsEIVKROqJyAIR+VZEVojI71xnMt/36sLtPPnZBkYPSeLG05JdxzEmcLq1aki3xIa89a01/1WnQBQqIB8Yoap9gX7ABSJyquNMppQFm/Zx7+tLOa1zc353SU87w8+YCozs05qFW/az+2Ce6yi1RiAKlX8rkGz/z1j/n/XwGBCbs3K45bl02jeL48nRA4m1jmaNqdBFvVsD8O5ya/6rLoHZ4ohItIgsATKAD1V1fjnvmSQi6SKSnplpV4DXhIO5hdw4Iw0Fpo4fROM462jWmGPp3LIB3RIb2nGqahSYQqWqxaraD2gHDC7vvlWq+oyqpqpqakKC3TE21AqLS7ht9iK27cvlqTEDSW4R7zqSMWHhot6tSd+ynz2HrPmvOgSmUB2lqgeAz4ALHEeJaKrKb99awdz1Wfzx8t6cmtLcdSRjwsbIPq1QhXdtr6paBKJQiUiCiDTxH9cHzgVWu00V2aZ9tZnZ87fyo7M6cXVqe9dxjAkrnVs2pGtiA+Ys2+06Sq0QiEIFtAY+FZGlQBreMaq3HWeKWJ+s3sMD76zkBz0T+b/zu7mOY0xYuqh3a9K27CPDmv9OWiAKlaouVdX+qtpHVXupqnVg68iqXYf46ezF9GjTiL9f24+oKDsN3ZgTMbJ3a6/5b7ntVZ2sQBQqEwwZh/KYOD2NhvVimTxuEHF1rAdoY05Ul8SGdGnZgHfsONVJs0JlADhSUMzNM9PZn1vI5PGptGpsHc0ac7Iu6t2atM3W/HeyrFAZr6PZl5ewdMdB/jmqP73aNnYdyZhaYWQfr/nv/ZV7XEcJa1aoDH/9YA3vLt/Nry86hfN6WEezxlSXLi0b0LFFPB9aoTopVqgi3Mvp23jysw2MGpzExNM7uo5jTK0iIpzfI5F5G7I4lFfoOk7YskIVweZt2Muv/72M0zu34PeXWkezxoTCeT0SKSxWPltj3b6dKCtUEWpjZja3zlpIh+bxPHH9AOto1pgQ6Z/UlBYN6vDBCjtN/UTZ1ikCHcgtYOKMdKKjxOtotr51NGtMqERHCeeekshnazLJLyp2HScsWaGKMAVFJdzy3EJ27D/CM2MHktQ8znUkY2q983smkp1fxDcb97mOEpasUEUQVeVX/17G/E37+PNVfUhNbuY6kjERYVinFsTVibbmvxNkhSqC/OvzDby6cDu3n9OFy/q3dR3HmIhRLzaas7ol8OHKPZSU2D1hq8oKVYSYs2wXf35vDZf0bcOd53ZxHceYiHN+j1ZkHM7n2+0HXEcJO1aoIsC32w5w50tLGJDUhD9f1cdOQzfGgbO7tSQmSvjALv6tMitUtdyOA0e4aWY6CQ3r8sy4VOrFRruOZExEahwXy6kpzXnfjlNVmRWqWuxQXiE3Tksjr6CYaRMG0aJBXdeRjIlo557Sko2ZOWzOynEdJaxYoaqlCotL+PGsRWzIzOapsQPpktjQdSRjIt6I7l5fmp+sznCcJLxYoaqFVJV7X1/G3PVZ/OmK3pzWuYXrSMYYIKl5HJ0S4q1QVZEVqlrosU/W//c09KtT27uOY4wpZUT3lszftJfs/CLXUcJGIAqViLQXkU9FZJWIrBCRO1xnClevL9rOIx+u5YoBbe00dGMCaER3r5PaueuyXEcJG4EoVEARcJeqngKcCtwmIj0cZwo7X6/P4u7XljKsU3MeusJOQzcmiFKTm9KwXgyfWvNfpQWiUKnqLlVd5D8+DKwCrOuEKli35zC3zFpIcvN4/jVmIHViAjFrjTFlxEZHMbxLAp+uybBeKiopcFszEUkG+gPzy3ltkoiki0h6Zqbd2+WojMN5TJiWRr3YaKbdYL2hGxN0Z3dvScbhfFbsPOQ6SlgIVKESkQbAa8DPVPV7c1BVn1HVVFVNTUhIqPmAAZSTX8TE6ensyylg6vhBtGtqvaEbE3RndUtAxE5Tr6zAFCoRicUrUs+r6uuu84SD4hLl9hcWs2LnQR4f3Z/e7Rq7jmSMqYQWDerSt10TPlljhaoyAlGoxDvqPwVYpaqPuM4TDlSV3/1nBR+vzuB3l/TknFMSXUcyxlTBiO4tWbr9AFnZ+a6jBF4gChVwGjAWGCEiS/x/F7kOFWTPfrmRmfO2MGl4CmOHJruOY4ypohHdW6IKn62x4+3HE+M6AICqzgXsXOpKemPxDv44ZzUj+7Tmngu6u45jjDkBPds0omXDuny6OoOrBrZzHSfQgrJHZSrpy3WZ/OKVbzk1pRmPXNOXqCir78aEIxFheNcE5q7PothOUz8mK1RhZPmOg9z63EI6t2zAM+NSqRtjt+wwJpwN75rAwSOFLLWbKR6TFaowsXVvLhOmLaBJXB1m3DiYRvXsWiljwt0ZnVsgAl+ste6UjsUKVRjIys5n3NT5FJUoM24cTGKjeq4jGWOqQdP4OvRp25gv1tkJFcdihSrgvAt609h9KI8p4wfRuWUD15GMMdXojC4JLNl2gINHCl1HCSwrVAFWWFzCbbMXsWzHQR4fNYCBHZq6jmSMqWbDuyZQXKLM22DNfxWxQhVQqso9ry3jszWZ/PHy3pzbwy7oNaY26p/UhAZ1Y/jcjlNVyApVQP3l/TW8tmg7d57blesGJ7mOY4wJkdjoKIZ1as4XazNRtdPUy2OFKoAmf7mRJz/bwKjBSdx+TmfXcYwxITa8awI7DhxhU1aO6yiBZIUqYF5O38YD76ziot6teOCyXnbzQ2MiwJldvbtBfLHWzv4rjxWqAHlv+W7ueW0pZ3Rpwd+v7Ue09TphTERo3yyO5OZxfGG3py+XFaqA+Gp9Fre/sJi+7Zvw1JiB1uuEMRFmeNcE5m3YS35RsesogWOFKgAWb93PzTPT6dginmkTBhFfNxB9BRtjatDwLgkcKSxm4Zb9rqMEjhUqx9buOcwN09No0aAuz00cTJO4Oq4jGWMcGJLSjOgo4ev1e11HCRwrVA5t25fL2CnzqRMdxayJQ2hpXSMZE7Ea1oulb7vGfGUX/n6PFSpHMg7nMWbKfPIKS3hu4hCSmse5jmSMcWxYpxYs3X6Qw3nWnVJpVqgc2Judz/XPzifjUD7TbhhEt1YNXUcyxgTAsM7NKS5RFmza5zpKoFihqmEHcgsYO2UBW/flMmV8KgOSrP8+Y4xnQFJT6sZE8ZUdp/qOQBQqEZkqIhkistx1llA6lFfI+KkLWJ+RzdNjBzKscwvXkYwxAVIvNprU5KZ8bcepviMQhQqYDlzgOkQo5eQXccO0NFbsPMST1w/grG4tXUcyxgTQsE4tWL37MFnZ+a6jBEYgCpWqfgHU2kbZIwXFTJyRxpJtB3hsVH/rCd0YU6HT/JaWeRus+e+oQBSqyhKRSSKSLiLpmZnh0SdWXmExk55LZ/6mfTxyTV8u7N3adSRjTID1atOIhnVjrPmvlLAqVKr6jKqmqmpqQkKC6zjHVVBUwm3PL+LLdVk8fGUfLu3X1nUkY0zAxURHMSSluZ1QUUpYFapwkl9UzI+fX8jHqzN44LJeXJPa3nUkY0yYOK1zc7buy2XbvlzXUQLBClUI5BUWc+tzC/loVQZ/uLQnY07t4DqSMSaMDOtkx6lKC0ShEpEXgHlANxHZLiITXWc6UXmFxdzy3EI+XZPJg5f3YuzQZNeRjDFhpmtiA1o0qGvdKfkC0U23qo5ynaE65BUWc/PMdOauz+LhK3tz7SC7hbwxpupEhGGdmvP1hr2oasTfQDUQe1S1wZGCYm6cnsbc9Vn8+co+VqSMMSfl1JTmZB7Ot9vTY4WqWuTkF3HD9AV8s3Evj1zTl6vtxAljzEka3LEZAPOt3z8rVCfr4JFCxk1dwIJN+/j7tf24vH8715GMMbVAp4R4WjSoy/yNdkJFII5RhavMw/mMm7qA9RmHeWL0ALuY1xhTbUSEIR2bMX/Tvog/TmV7VCdo+/5crnl6HpuzcpgyfpAVKWNMtRuS0oxdB/PYtu+I6yhOWaE6Aeszsrn6qXlkZecz66bBDO8a/F4yjDHhZ0jH5gDM3xTZzX9WqKpo+Y6DXPP0PAqLS3hp0lAGdmjmOpIxppbq0rIBTeJiI/6ECjtGVQXzNuxl0sx0GtWPZdZNQ+jYIt51JGNMLRYVJQxObmZ7VK4DhIu3vt3J+KkLSGxcj1duHWpFyhhTI4akNGfbviPsPBC5x6msUB2HqvLMFxu4/YXF9Etqwmu3DqNNk/quYxljIsSQ/15PFbl7VVaojqG4RPndf1byxzmrGdmnNTNvHEzjuFjXsYwxEeSU1o1oWC+GBRF8nMqOUVUgr7CYO15czPsr9nDzGR2598JTiIqK3OsYjDFuREcJg5KbMX9j5BYq26MqR+bhfEY/+w0frNzDfRf34Ncje1iRMsY4M6RjMzZm5ZBxKM91FCesUJWxYudBLn18Lit3HeJf1w/gxtM7uo5kjIlwQ1KOXk8VmXtVVqhKeW/5Lq761zwUePXWYVzQy3qbMMa416tNI+LrREfscSo7RoV3Zt/jn6znbx+upV/7JjwzdiAtG9VzHcsYYwCIiY5iQIempG2OzEIV8XtUuQVF/PSFxfztw7Vc3r8tL0461YqUMSZwUjs0Y82ewxw8Uug6So2L6D2qDZnZ/GjWQtZlZHP3Bd259cyUiO6h2BgTXIOSm6IKi7bs5+zuLV3HqVGB2aMSkQtEZI2IrBeRe0I9vDnLdnHJY3PJyi5g5o2D+dFZnaxIGWMCq19SE2KiJCKb/wKxRyUi0cATwHnAdiBNRN5S1ZXVPazC4hIeenc1U+Zuon9SE54YPcB6mjDGBF5cnRh6tm1M+ub9rqPUuEAUKmAwsF5VNwKIyIvApUC1Fqrs/CImTF1A+pb9TBiWzK8uOoU6MYHZqTTGmGMa1KEpz32zhYKikojadgWlULUFtpX6ezswpOybRGQSMAkgKSmpygOJrxNNcot4xg7twKX92p5gVGOMcWPS8BR+fHbniCpSEJxCVd7BIf3eE6rPAM8ApKamfu/14w5EhL9e3bfq6YwxJgAi9YzkoJTl7UD7Un+3A3Y6ymKMMSZAglKo0oAuItJRROoA1wFvOc5kjDEmAALR9KeqRSLyE+B9IBqYqqorHMcyxhgTAIEoVACqOgeY4zqHMcaYYAlK058xxhhTLitUxhhjAs0KlTHGmECzQmWMMSbQRLXK180GgohkAltO8OMtgKxqjBMObJxrv0gbX7BxrqoOqppQnWFqQtgWqpMhIumqmuo6R02yca79Im18wcY5UljTnzHGmECzQmWMMSbQIrVQPeM6gAM2zrVfpI0v2DhHhIg8RmWMMSZ8ROoelTHGmDBhhcoYY0ygRVShEpELRGSNiKwXkXtc5wk1EWkvIp+KyCoRWSEid7jOVFNEJFpEFovI266z1AQRaSIir4rIan9+D3WdKdRE5E5/uV4uIi+ISK27q6CITBWRDBFZXuq5ZiLyoYis8/9v6jJjTYiYQiUi0cATwIVAD2CUiPRwmyrkioC7VPUU4FTgtggY56PuAFa5DlGD/gG8p6rdgb7U8nEXkbbA7UCqqvbCuz3QdW5ThcR04IIyz90DfKyqXYCP/b9rtYgpVMBgYL2qblTVAuBF4FLHmUJKVXep6iL/8WG8jVdbt6lCT0TaASOBya6z1AQRaQQMB6YAqGqBqh5wm6pGxAD1RSQGiKMW3hVcVb8A9pV5+lJghv94BnBZjYZyIJIKVVtgW6m/txMBG+2jRCQZ6A/Md5ukRjwK/BIocR2khqQAmcA0v7lzsojEuw4VSqq6A/grsBXYBRxU1Q/cpqoxiaq6C7wfo0BLx3lCLpIKlZTzXEScmy8iDYDXgJ+p6iHXeUJJRC4GMlR1oessNSgGGAD8S1X7AznU8uYg/7jMpUBHoA0QLyJj3KYyoRJJhWo70L7U3+2ohU0FZYlILF6Rel5VX3edpwacBlwiIpvxmndHiMgst5FCbjuwXVWP7i2/ile4arNzgU2qmqmqhcDrwDDHmWrKHhFpDeD/n+E4T8hFUqFKA7qISEcRqYN34PUtx5lCSkQE77jFKlV9xHWemqCq96pqO1VNxpvHn6hqrf6lraq7gW0i0s1/6hxgpcNINWErcKqIxPnL+TnU8hNISnkLGO8/Hg+86TBLjYhxHaCmqGqRiPwEeB/vDKGpqrrCcaxQOw0YCywTkSX+c79S1TkOM5nQ+CnwvP8jbCNwg+M8IaWq80XkVWAR3tmti6mFXQuJyAvAWUALEdkO/BZ4CHhZRCbiFeyr3SWsGdaFkjHGmECLpKY/Y4wxYcgKlTHGmECzQmWMMSbQrFAZY4wJNCtUxhhjAs0KlTHGmECzQmWMMSbQrFAZUw1EZJCILBWReiIS798nqZfrXMbUBnbBrzHVREQeAOoB9fH63vuT40jG1ApWqIypJn73RWlAHjBMVYsdRzKmVrCmP2OqTzOgAdAQb8/KGFMNbI/KmGoiIm/h3VqkI9BaVX/iOJIxtULE9J5uTCiJyDigSFVni0g08LWIjFDVT1xnMybc2R6VMcaYQLNjVMYYYwLNCpUxxphAs0JljDEm0KxQGWOMCTQrVMYYYwLNCpUxxphAs0JljDEm0P4/QtQrk7tXv6cAAAAASUVORK5CYII=\n", "text/plain": [ "

" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#Define desired spatial map and step\n", "Lx = 10\n", "Nx = 2000\n", "dx = Lx/Nx\n", "\n", "#Define functions that will be used by numerical solver \n", "#source function \n", "Source_func = lambda x, q: q*x*(10-x)\n", "#reaction constants \n", "kappa_func = lambda x, kappal, kappah: kappah + (kappal-kappah)*(x>5)*(x<7.5)\n", "#Initialize these to the correct size \n", "v_func = lambda x,v: v*np.ones(x.size)\n", "omega_func = lambda x,omega: omega*np.ones(x.size)\n", "\n", "#Define array for plotting purposes \n", "xs = np.linspace(dx/2,Lx-dx/2,Nx)\n", "\n", "# evaluate source and kappa\n", "source = Source_func(xs, 1)\n", "kappa = kappa_func(xs, 0.1, 2)\n", "\n", "#nominal values and given variances \n", "omega_nom = 20\n", "omega_var = 0.3195214\n", "v_nom = 10\n", "v_var = 0.0723493\n", "kappal_nom = 0.1\n", "kappal_var = 8.511570e-6\n", "kappah_nom = 2\n", "kappah_var = 0.002778142\n", "q_nom = 1\n", "q_var = 7.062353e-4\n", "vs = v_func(xs, v_nom)\n", "# print(vs)\n", "\n", "#Call Numerical Function above to solve equation \n", "sol,Q = ADRSource(Lx, Nx, source, omega_func(xs, omega_nom), vs, kappa)\n", "# print(Q)\n", "\n", "#Plot solution \n", "plt.plot(xs,sol)\n", "plt.xlabel('x')\n", "plt.ylabel('u(x)')\n", "plt.title('The solution u(x) evaluated at mean value of uncertain parameters ')\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "T6ugw9X1ZF4z", "nbpages": { "level": 3, "link": "[4.2.3 **4.3.1 Simple ADR Example** ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.02-Simple-ADR-Example.html#4.2.3-**4.3.1-Simple-ADR-Example**)", "section": "4.2.3 **4.3.1 Simple ADR Example** " } }, "source": [ "**Sensitivity Analysis** " ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "2M96vCvlfyQ4", "nbpages": { "level": 3, "link": "[4.2.3 **4.3.1 Simple ADR Example** ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.02-Simple-ADR-Example.html#4.2.3-**4.3.1-Simple-ADR-Example**)", "section": "4.2.3 **4.3.1 Simple ADR Example** " } }, "source": [ "Sensitivity using Forward Difference Derivative Approximation: \n", "\n", "$$\n", "\\frac{\\partial Q}{\\partial x_i} \\bigg|_{\\overline{x}} \\approx \\frac{Q(\\overline{x} + \\delta_i\\hat{e_i}) - Q(\\overline{x})}{\\delta_i}\n", "$$\n", "Scaled Sensitivity Coefficient: \n", "\n", "$$\n", "\\overline{x_i}\\frac{\\delta Q}{\\delta x_i}\\bigg|_\\overline{x}\n", "$$\n", "\n", "Sensitivity Index:\n", "$$\n", "\\sigma_i\\frac{\\delta Q}{\\delta x_i}\\bigg|_\\overline{x}\n", "$$" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 530 }, "colab_type": "code", "id": "Gkfvtj6rZEiI", "nbpages": { "level": 3, "link": "[4.2.3 **4.3.1 Simple ADR Example** ](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.02-Simple-ADR-Example.html#4.2.3-**4.3.1-Simple-ADR-Example**)", "section": "4.2.3 **4.3.1 Simple ADR Example** " }, "outputId": "2249e1f4-ad01-40ce-8fa0-3023e8a82a30" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "######################### Sensitivity Analysis #################\n", "v sensitivity: (-1.7406387549101512+0j)\n", "v scaled sensitivity: (-17.406387549101513+0j)\n", "v Sensitivity index: (-0.46819396950072134+0j)\n", "\n", "\n", "omega sensitivity: (-0.970203846506479+0j)\n", "omega scaled sensitivity: (-19.40407693012958+0j)\n", "omega sensitivity index: (-0.5484195995587706+0j)\n", "\n", "\n", "Kappa_l sensitivity: (12.867971008745371+0j)\n", "Kappa_l scaled sensitivity: (1.2867971008745371+0j)\n", "Kappa_l sensitivity index: (0.0375417844104318+0j)\n", "\n", "\n", "Kappa_h sensitivity: (17.761271102756382+0j)\n", "Kappa_h scaled sensitivity: (35.522542205512764+0j)\n", "Kappa_h sensitivity index: (0.9361625491891298+0j)\n", "\n", "\n", "q sensitivity: (52.39025306025269+0j)\n", "q scaled sensitivity: (52.39025306025269+0j)\n", "q sensitivity index: (1.3922755832423062+0j)\n", "\n", "\n" ] } ], "source": [ "\n", "#v sensitivity\n", "#finite difference step amount \n", "delta = 1e-6\n", "vs_pert = v_func(xs, v_nom*(1+delta))\n", "#Calculate second point for finite difference\n", "sol,Q_new = ADRSource(Lx, Nx, source, omega_func(xs, omega_nom), vs_pert, kappa)\n", "#finite difference Sensitivity Calculation \n", "sens_v = (Q_new - Q)/(v_nom*delta)\n", "#Scaled Sensitivity \n", "scaled_v = sens_v*v_nom\n", "#Sensitivity Index \n", "index_v = math.sqrt(v_var)*sens_v\n", "#Print Commands \n", "print('######################### Sensitivity Analysis #################')\n", "print(\"v sensitivity: \", sens_v)\n", "print(\"v scaled sensitivity: \", scaled_v)\n", "print(\"v Sensitivity index: \", index_v)\n", "print('\\n')\n", "\n", "#omega sensitivity\n", "#Calculate second point for finite difference \n", "sol,Q_new = ADRSource(Lx, Nx, source, omega_func(xs, omega_nom*(1+delta)), vs, kappa)\n", "#Sensitivity calculation using finite difference \n", "sens_omega = (Q_new - Q)/(omega_nom*delta)\n", "#Scaled sensitivity Calc\n", "scaled_v = sens_omega*omega_nom\n", "#Sensitivity Index Calc\n", "index_v = math.sqrt(omega_var)*sens_omega\n", "#Print Commands \n", "print(\"omega sensitivity: \", sens_omega)\n", "print(\"omega scaled sensitivity: \", scaled_v)\n", "print(\"omega sensitivity index: \", index_v)\n", "print('\\n')\n", "\n", "#kappal sensitivity\n", "#Calculate second point value for finite difference \n", "sol,Q_new = ADRSource(Lx, Nx, source, omega_func(xs, omega_nom), vs, kappa_func(xs, kappal_nom*(1+delta), 2))\n", "#Sensitivity Calculation using finite difference \n", "sens_kappal = (Q_new - Q)/(kappal_nom*delta)\n", "#Scaled Sensitivity Calc \n", "scaled_v = sens_kappal*kappal_nom\n", "#Sensitivity Index Calc\n", "index_v = math.sqrt(kappal_var)*sens_kappal\n", "#Print Commands \n", "print(\"Kappa_l sensitivity: \", sens_kappal)\n", "print(\"Kappa_l scaled sensitivity: \", scaled_v)\n", "print(\"Kappa_l sensitivity index: \", index_v)\n", "print('\\n')\n", "\n", "#kappah sensitivity\n", "#Calculate second point value for finite difference \n", "sol,Q_new = ADRSource(Lx, Nx, source, omega_func(xs, omega_nom), vs, kappa_func(xs, kappal_nom, kappah_nom*(1+delta)))\n", "#Sensitivity Calculation using finite difference \n", "sens_kappah = (Q_new - Q)/(kappah_nom*delta)\n", "#Scaled Sensitivity Calc\n", "scaled_v = sens_kappah*kappah_nom\n", "#Sensitivity Index Calc\n", "index_v = math.sqrt(kappah_var)*sens_kappah\n", "#Print Commands \n", "print(\"Kappa_h sensitivity: \", sens_kappah)\n", "print(\"Kappa_h scaled sensitivity: \", scaled_v)\n", "print(\"Kappa_h sensitivity index: \", index_v)\n", "print('\\n')\n", "\n", "#q sensitivity\n", "#Calculate second point value for finite difference \n", "sol,Q_new = ADRSource(Lx, Nx, Source_func(xs, q_nom*(1+delta)), omega_func(xs, omega_nom), vs,kappa_func(xs, kappal_nom, kappah_nom))\n", "sens_q = (Q_new - Q)/(q_nom*delta)\n", "#Scaled Sensitivity Calc\n", "scaled_v = sens_q*q_nom\n", "#Sensitivity Index Calc\n", "index_v = math.sqrt(q_var)*sens_q\n", "#print commands\n", "print(\"q sensitivity: \", sens_q)\n", "print(\"q scaled sensitivity: \", scaled_v)\n", "print(\"q sensitivity index: \", index_v)\n", "print('\\n')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "< [4.1 Local sensitivity analysis and difference approximation](https://ndcbe.github.io/cbe67701-uncertainty-quantification/04.01-Local-sensitivity-analysis.html) | [Contents](toc.html) | [5.0 Regression Approximations to Estimate Sensitivities](https://ndcbe.github.io/cbe67701-uncertainty-quantification/05.00-Regression-Approximations-to-Estimate-Sensitivities.html)

\"Open

\"Download\"" ] } ], "metadata": { "colab": { "collapsed_sections": [], "name": "04.02.Simple ADR Example.ipynb", "provenance": [] }, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 1 }