{ "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", "< [6.0 Adjoint-Based Local Sensitivity Analysis](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.00-Adjoint-Based-Local-Sensitivity-Analysis.html) | [Contents](toc.html) | [6.2 A Simple Example of Adjoint Sensitivity Analysis](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.02-Contributed-Example.html)

\"Open

\"Download\"" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "5F6q8D_2xysC", "nbpages": { "level": 1, "link": "[6.1 Nonlinear Diffusion-Reaction Equation](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1-Nonlinear-Diffusion-Reaction-Equation)", "section": "6.1 Nonlinear Diffusion-Reaction Equation" } }, "source": [ "# 6.1 Nonlinear Diffusion-Reaction Equation\n", "\n", "Created by Shanghui Huang (shuang7@nd.edu)\n", "\n", "This example was adapted from:\n", "\n", "McClarren, Ryan G (2018). Uncertainty Quantification and Predictive Computational Science: A Foundation for Physical Scientists and Engineers, Chapter 6: Adjoints-Based Local Sensitivity Analysis, Springer, https://doi.org/10.1007/978-3-319-99525-0_6\n", "\n", "And special acknowledgement to *Han Gao* for a fruitful discussion" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "nVkXoGG5MVjd", "nbpages": { "level": 2, "link": "[6.1.1 An Overview](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.1-An-Overview)", "section": "6.1.1 An Overview" } }, "source": [ "## 6.1.1 An Overview\n", "This chapter introduces an adjoint-based method to compute the local sensitivity of a QoI against a parameter. In this method, adjoints gives information of perturbed quantity with only one solve of the forward system and one solve of the adjoint equations.\n", "\n", "Take the steady ADR equation for example: \n", "\\begin{equation}\n", "v\\frac{\\partial u}{\\partial x} - \\omega \\frac{\\partial ^2 u}{\\partial x^2} + \\kappa u = q\n", "\\end{equation}\n", "the operator $L$ is defined as:\n", "\\begin{equation}\n", "L = v\\frac{\\partial}{\\partial x} - \\omega \\frac{\\partial ^2}{\\partial x^2} + \\kappa\n", "\\end{equation}\n", "Then the adjoints satisfies:\n", "\\begin{equation}\n", "(Lu,u^{\\ast})=(u,L^{\\ast}u^{\\ast})\n", "\\end{equation}\n", "and\n", "\\begin{equation}\n", "L^{\\ast}u^{\\ast} = w\n", "\\end{equation}\n", "After a series of mathematical derivations, we can calculate the local sensitivities:\n", "\\begin{equation}\n", "\\frac{\\partial Q}{\\partial \\theta} = (\\frac{\\partial q}{\\partial \\theta},u^{\\ast}) + (u,\\frac{\\partial w}{\\partial \\theta}) - (\\frac{\\partial L}{\\partial \\theta}u,u^{\\ast})\n", "\\end{equation} " ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "r0IDoJRGXT_-", "nbpages": { "level": 2, "link": "[6.1.2 Example of nonlinear diffusion-reaction equation](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.2-Example-of-nonlinear-diffusion-reaction-equation)", "section": "6.1.2 Example of nonlinear diffusion-reaction equation" } }, "source": [ "## 6.1.2 Example of nonlinear diffusion-reaction equation" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "8Eyf7kcNbGU9", "nbpages": { "level": 2, "link": "[6.1.2 Example of nonlinear diffusion-reaction equation](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.2-Example-of-nonlinear-diffusion-reaction-equation)", "section": "6.1.2 Example of nonlinear diffusion-reaction equation" } }, "source": [ "Solving nonlinear diffusion-reaction equation:\n", "\\begin{equation}\n", "F(u,\\dot{u}) = \\rho\\dot{u}-\\omega\\frac{\\partial^2u^4}{\\partial x^2} + \\kappa u^4 - S\n", "\\end{equation}\n", "in which \n", "\\begin{equation}\n", "\\kappa (x) =\n", "\\begin{cases}\n", "\\kappa_h, & 1 \\le x \\lt 1.5 \\\\\n", "\\kappa_l, & otherwise\n", "\\end{cases}\n", "\\end{equation}\n", "\\begin{equation}\n", "\\S(x) = \n", "\\begin{cases}\n", "q, & 0.5 \\le x \\le 1.5 \\\\\n", "0, & otherwise\n", "\\end{cases}\n", "\\end{equation}\n", "with nominal values \n", "\\begin{equation}\n", "\\rho = 1, \\omega = 0.1, \\kappa _l = 0.1 \\kappa _h =2, q = 1\n", "\\end{equation}" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "OFnYIf6MX_bw", "nbpages": { "level": 2, "link": "[6.1.3 Codes](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3-Codes)", "section": "6.1.3 Codes" } }, "source": [ "## 6.1.3 Codes" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 228 }, "colab_type": "code", "executionInfo": { "elapsed": 5038, "status": "ok", "timestamp": 1593523486948, "user": { "displayName": "黄尚晖", "photoUrl": "", "userId": "16496223364462243774" }, "user_tz": 240 }, "id": "L1kWRZFJ9M3n", "nbpages": { "level": 2, "link": "[6.1.3 Codes](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3-Codes)", "section": "6.1.3 Codes" }, "outputId": "2e20d19a-a99a-40f3-85e2-b6e7069e5221" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Collecting pyDOE\n", " Downloading https://files.pythonhosted.org/packages/bc/ac/91fe4c039e2744466621343d3b8af4a485193ed0aab53af5b1db03be0989/pyDOE-0.3.8.zip\n", "Requirement already satisfied, skipping upgrade: numpy in /usr/local/lib/python3.6/dist-packages (from pyDOE) (1.18.5)\n", "Requirement already satisfied, skipping upgrade: scipy in /usr/local/lib/python3.6/dist-packages (from pyDOE) (1.4.1)\n", "Building wheels for collected packages: pyDOE\n", " Building wheel for pyDOE (setup.py) ... \u001b[?25l\u001b[?25hdone\n", " Created wheel for pyDOE: filename=pyDOE-0.3.8-cp36-none-any.whl size=18178 sha256=eae884fee2b8e02ed0591f4026fbd6d99f8c3937a2ef7b62e330bfd5bc6672e9\n", " Stored in directory: /root/.cache/pip/wheels/7c/c8/58/a6493bd415e8ba5735082b5e0c096d7c1f2933077a8ce34544\n", "Successfully built pyDOE\n", "Installing collected packages: pyDOE\n", "Successfully installed pyDOE-0.3.8\n" ] } ], "source": [ "!pip install --upgrade pyDOE\n", "!pip install tabulate" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "colab": {}, "colab_type": "code", "executionInfo": { "elapsed": 498, "status": "ok", "timestamp": 1593523488826, "user": { "displayName": "黄尚晖", "photoUrl": "", "userId": "16496223364462243774" }, "user_tz": 240 }, "id": "OKtmTiFh6I9k", "nbpages": { "level": 2, "link": "[6.1.3 Codes](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3-Codes)", "section": "6.1.3 Codes" } }, "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", "from scipy.stats import gamma\n", "import math\n", "from scipy.stats.distributions import norm\n", "from scipy import stats\n", "from pyDOE import *\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "colab": {}, "colab_type": "code", "executionInfo": { "elapsed": 342, "status": "ok", "timestamp": 1593523492495, "user": { "displayName": "黄尚晖", "photoUrl": "", "userId": "16496223364462243774" }, "user_tz": 240 }, "id": "6H4HdAhJ6I9p", "nbpages": { "level": 2, "link": "[6.1.3 Codes](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3-Codes)", "section": "6.1.3 Codes" } }, "outputs": [], "source": [ "def ADRSource(Lx, Nx, Source, omega, v, kappa):\n", " #Solves the diffusion equation with Generalized Source\n", " A = sparse.dia_matrix((Nx,Nx),dtype=\"complex\")\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)\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", " #solve A x = Source\n", " Solution = linalg.spsolve(A,Source)\n", " Q = integrate.trapz(Solution*kappa,dx=dx)\n", " return Solution, Q" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "h1lNWlMYZDQn", "nbpages": { "level": 3, "link": "[6.1.3.1A function calculating solution $u$](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.1A-function-calculating-solution-$u$)", "section": "6.1.3.1A function calculating solution $u$" } }, "source": [ "### 6.1.3.1A function calculating solution $u$" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "colab": {}, "colab_type": "code", "executionInfo": { "elapsed": 528, "status": "ok", "timestamp": 1593523515395, "user": { "displayName": "黄尚晖", "photoUrl": "", "userId": "16496223364462243774" }, "user_tz": 240 }, "id": "07NSOk9t6I9t", "nbpages": { "level": 3, "link": "[6.1.3.1A function calculating solution $u$](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.1A-function-calculating-solution-$u$)", "section": "6.1.3.1A function calculating solution $u$" } }, "outputs": [], "source": [ "def run_calc(Lx, Nx, rho, Source, omega, v, \n", " kappa, dt, steps,\n", " tmin,tmax,xmin,xmax, IC):\n", " dx = Lx/Nx\n", " assert IC.size==Nx\n", " irho = 1/rho\n", " dt2 = dt/2\n", " idx2 = 1/dx**2\n", " solution = np.zeros((Nx+2,steps+1))\n", " solution[1:Nx+1,0] = 1.0*IC #copy IC\n", " xs = np.linspace(-dx/2,Lx+dx/2,Nx+2)\n", " Q = 0\n", " for step in range(steps):\n", " #predictor\n", " pred = solution[:,step]*1.0\n", " \n", " pred[1:Nx+1] = (pred[1:Nx+1]-irho*kappa*dt2*pred[1:Nx+1]**4 +\n", " irho*dt2*Source +\n", " irho*dt2*omega*idx2*(pred[2:Nx+2]**4 - 2*pred[1:Nx+1]**4 + pred[0:Nx]**4))\n", " \n", " #corrector\n", " solution[1:Nx+1,step+1] = (solution[1:Nx+1,step]-irho*kappa*dt*pred[1:Nx+1]**4 + \n", " irho*dt*Source +\n", " irho*dt*omega*idx2 * (pred[2:Nx+2]**4 - 2*pred[1:Nx+1]**4 + pred[0:Nx]**4))\n", " if ((step+1)*dt >= tmin) and ((step+1)*dt <= tmax):\n", " Q += np.sum(dt*solution[1:Nx+1,step][(xs[1:Nx+1]>=xmin)*(xs[1:Nx+1]<=xmax)]*\n", " kappa[(xs[1:Nx+1]>=xmin)*(xs[1:Nx+1]<=xmax)]*dx)\n", " return solution,Q,xs" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "mazDyC_6Yo14", "nbpages": { "level": 3, "link": "[6.1.3.2 Setting up the nonlinear diffusion-reaction equation](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.2-Setting-up-the-nonlinear-diffusion-reaction-equation)", "section": "6.1.3.2 Setting up the nonlinear diffusion-reaction equation" } }, "source": [ "### 6.1.3.2 Setting up the nonlinear diffusion-reaction equation" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "colab": {}, "colab_type": "code", "executionInfo": { "elapsed": 3314, "status": "ok", "timestamp": 1593523526457, "user": { "displayName": "黄尚晖", "photoUrl": "", "userId": "16496223364462243774" }, "user_tz": 240 }, "id": "kfi1jUS36I9y", "nbpages": { "level": 3, "link": "[6.1.3.2 Setting up the nonlinear diffusion-reaction equation](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.2-Setting-up-the-nonlinear-diffusion-reaction-equation)", "section": "6.1.3.2 Setting up the nonlinear diffusion-reaction equation" } }, "outputs": [], "source": [ "#set up simple advection problem\n", "Lx = 2\n", "Nx = 200\n", "dx = Lx/Nx\n", "dt = dx**2\n", "#print('Timestep = ',dt)\n", "tfinal = 2\n", "numsteps = math.ceil(tfinal/dt)\n", "xs = np.linspace(dx/2,Lx-dx/2,Nx)\n", "Source_func = lambda x, q: (x>0.25*Lx)*(x<0.75*Lx)*q\n", "kappa_func = lambda x, kappal, kappah: kappah + (kappal-kappah)*(x>0.5*Lx)*(x<0.75*Lx)\n", "#print('Number of steps = ', numsteps)\n", "IC = np.zeros(Nx)\n", "IC[Nx//4:3*Nx//4] = 1*0\n", "source = IC*0\n", "source[Nx//4:3*Nx//4] = 10\n", "\n", "omega_nom = .1\n", "v_nom = 0.0\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", "rho = 1.0\n", "\n", "sol,Q,xs2 = run_calc(Lx, Nx, rho, Source_func(xs,1), omega_nom, v_nom, \n", " kappa_func(xs,kappal_nom, kappah_nom), dt, numsteps, 1.8,tfinal,1.5,1.9,IC)\n", "#print(Q)\n", "xs3 = xs2\n", "#plt.plot(xs2,sol[:,numsteps])" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "I3y7UScJOOQV", "nbpages": { "level": 3, "link": "[6.1.3.3 A function calculating adjoint solution $u^{\\ast}$](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.3-A-function-calculating-adjoint-solution-$u^{\\ast}$)", "section": "6.1.3.3 A function calculating adjoint solution $u^{\\ast}$" } }, "source": [ "### 6.1.3.3 A function calculating adjoint solution $u^{\\ast}$" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "colab": {}, "colab_type": "code", "executionInfo": { "elapsed": 317, "status": "ok", "timestamp": 1593523531990, "user": { "displayName": "黄尚晖", "photoUrl": "", "userId": "16496223364462243774" }, "user_tz": 240 }, "id": "sF1qggY56I93", "nbpages": { "level": 3, "link": "[6.1.3.3 A function calculating adjoint solution $u^{\\ast}$](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.3-A-function-calculating-adjoint-solution-$u^{\\ast}$)", "section": "6.1.3.3 A function calculating adjoint solution $u^{\\ast}$" } }, "outputs": [], "source": [ "def run_adjoint_calc(Lx, Nx, rho, Source, omega, v, kappa, dt, steps,tmin,tmax,xmin,xmax, IC, forward):\n", " dx = Lx/Nx\n", " assert IC.size==Nx\n", " dt2 = dt/2\n", " irho = 1/rho\n", " solution = np.zeros((Nx+2,steps+1))\n", " solution[1:Nx+1,steps] = 1.0*IC #copy IC\n", " xs = np.linspace(-dx/2,Lx+dx/2,Nx+2)\n", " Q = 0\n", " ks = kappa*(xs[1:Nx+1]>=xmin)*(xs[1:Nx+1]<=xmax)\n", " for step in range(steps-1,-1,-1):\n", " #predictor\n", " pred = solution[:,step+1]*1.0\n", " pred[1:Nx+1] = (pred[1:Nx+1]*(1- irho*4*kappa*dt2*forward[1:Nx+1,step]**3) + \n", " irho*dt2*ks*((step+1)*dt >= tmin)*((step+1)*dt <= tmax) +\n", " irho*4*dt2*omega/(dx**2)*forward[1:Nx+1,step]**3 * (pred[2:Nx+2]-2*pred[1:Nx+1] + pred[0:Nx])) \n", " \n", " \n", " #corrector\n", " solution[1:Nx+1,step] = (solution[1:Nx+1,step+1]-irho*4*kappa*dt*pred[1:Nx+1]*forward[1:Nx+1,step]**3 + \n", " irho*dt*ks*((step+1)*dt >= tmin)*((step+1)*dt <= tmax)+\n", " irho*4*dt*omega/(dx**2)*forward[1:Nx+1,step]**3 * (pred[2:Nx+2]-2*pred[1:Nx+1] + pred[0:Nx])) \n", " return solution,xs" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "colab": {}, "colab_type": "code", "executionInfo": { "elapsed": 2359, "status": "ok", "timestamp": 1593523536988, "user": { "displayName": "黄尚晖", "photoUrl": "", "userId": "16496223364462243774" }, "user_tz": 240 }, "id": "BD4mRrzr6I97", "nbpages": { "level": 3, "link": "[6.1.3.3 A function calculating adjoint solution $u^{\\ast}$](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.3-A-function-calculating-adjoint-solution-$u^{\\ast}$)", "section": "6.1.3.3 A function calculating adjoint solution $u^{\\ast}$" } }, "outputs": [], "source": [ "sol_adj,xs2 = run_adjoint_calc(Lx, Nx,rho, Source_func(xs,1), omega_nom, v_nom, \n", " kappa_func(xs,kappal_nom, kappah_nom), dt, numsteps, 1.8,tfinal,1.5,1.9,IC*0, sol)\n", "\n", "output = np.vstack([xs2,sol[:,numsteps], sol_adj[:,round(1.8/dt)]]).transpose()\n", "np.savetxt('nonlindr_foradj.csv', output, delimiter=',')" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "_pLnbneDZTUX", "nbpages": { "level": 3, "link": "[6.1.3.4 Results of foward and adjoint solutions](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.4-Results-of-foward-and-adjoint-solutions)", "section": "6.1.3.4 Results of foward and adjoint solutions" } }, "source": [ "### 6.1.3.4 Results of foward and adjoint solutions" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 283 }, "colab_type": "code", "executionInfo": { "elapsed": 432, "status": "ok", "timestamp": 1593523575994, "user": { "displayName": "黄尚晖", "photoUrl": "", "userId": "16496223364462243774" }, "user_tz": 240 }, "id": "vxYcVslCV0NR", "nbpages": { "level": 3, "link": "[6.1.3.4 Results of foward and adjoint solutions](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.4-Results-of-foward-and-adjoint-solutions)", "section": "6.1.3.4 Results of foward and adjoint solutions" }, "outputId": "f0daaca1-a26d-4875-c8c4-a0f6c696370a" }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 16, "metadata": { "tags": [] }, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3xUVfr48c+ZSSa9J5QAKUAooRdBEBGUKgqIfdm1L6tr13V1Rf26fnVX3bWt7ftDV9cG2FkUUEFAQUTpLZSEnhBIgQTSk5nz++NOYAgJJDCTOzN53q9XXpm5986dh5vhyclzzj1Haa0RQgjh+yxmByCEEMI9JKELIYSfkIQuhBB+QhK6EEL4CUnoQgjhJwLMeuP4+HidkpJi1tsLIYRPWrNmTYHWOqG+faYl9JSUFFavXm3W2wshhE9SSu1taJ+UXIQQwk9IQhdCCD8hCV0IIfzEGWvoSql3gMuAPK11z3r2K+AV4FKgDLhJa73W3YEKIVqe6upqsrOzqaioMDuUZhccHEz79u0JDAxs9Gsa0yn6H+A14P0G9o8H0pxfg4E3nd+FEOKcZGdnExERQUpKCkbbsWXQWlNYWEh2djapqamNft0ZSy5a6x+Bw6c5ZBLwvjasBKKVUm0bHYEQQjSgoqKCuLi4FpXMAZRSxMXFNfkvE3fU0NsB+12eZzu3nUIpNU0ptVoptTo/P98Nby2E8HctLZnXOpt/d7OOQ9dazwBmAAwcOFDm7RUtSnFZNfsOl5FfUkFRWTVHyqqpqLZjd2jsDk2ARREWFEBYkJWoEBtto4JpExVMfHgQVkvLTGqiadyR0HOADi7P2zu3CdEiORyarQePsmF/MZsPFLPlwFF255dwtKLmrM4XFGChc6twurSOIL1tJANSYuiZGIUtQAapiZO5I6HPBe5SSs3G6Awt1lrnuuG8QviMvKMVfJdxiOWZBazcXUhRWTUAEcEB9EyMYlLfdiTHhdIhNpRWEUHEhNqIDg0kONBKgEVhtSiq7ZqyqhpKKmsoKqsmt7iCg8Xl7C0sY0deCSt3FfLlOqOtZAuwMCAphpHdEri4Wys6JYS32NKEOKExwxZnASOAeKVUNvA/QCCA1vr/gPkYQxazMIYt3uypYIXwJodLq5izLod5m3JZu+8IWkO76BBGd2/N0M5xDEiKpUNsSKMTrS1AYQuwER1qo30M9GwXdcoxeccqWLv3CKv3HGF5VgF/m7+Nv83fRmp8GJP6JnJFv3Ykx4W5+58qGlBeXs64ceNYvHgxVqu13mOqqqoYNWoUixcvJiDAs1XuM55da339GfZr4E63RSSEF9Na8/POQmb+uo/vthyiyu6ge9tI7h/VhfE929C5lWdbyq0ighnXsy3jehoDyXKKylmyLY95G3N55ftMXl6UyXkpMdw4NIWxPdoQaJWyjCe98847TJkypcFkDmCz2bjkkkv4+OOPmTp1qkfjMW1yLiF8SY3dwYLNB/l/P+5kc85RokMDmXp+Etee14FubSJNi6tddAi/PT+Z356fzIGicv67/gCzft3HXTPX0TYqmBuGpPC7IcmEB/n+f/W/frWFjANH3XrO9MRI/ufyHmc8bsiQIcycOZPU1FRycnKYOHEia9as4aOPPmLmzJnHjxs5ciSPPvooo0eP5rHHHqO4uJhXX32VyZMn85e//EUSuhBmW7I9j6e+ymB3QSkd48N4dkovJvdrR3Bgw60yMyRGh3DHiE5MG96RJdvyeHfFbp77ZhtvLdvFtOEduWFIMqE2+S/fVA6Hg71791I73ffGjRvp3bs3VVVV7Nq1C9dpwP/617/yxBNPkJeXx7p165g7dy4APXv2ZNWqVR6PVX66QjRg/+Eynvo6g4UZh+iYEMb//bY/Y9LbYPHyIYRWi2JUemtGpbdm3b4jvLwok2cXbOPtZbt5eFxXruzf3uv/DfVpTEvaE3bu3ElqaurxUtrGjRvp1asXBQUFREdHn3Ts8OHD0Vrz4osvsnTp0uOlGKvVis1m49ixY0RERHgsVimwCVFHRbWdf32fyagXf2B5ZgEPj+vGN/cOZ1zPtj6XCPslxfDeLYP4/I4hdIgN4aHPNjLlzRVs2F9kdmg+Y9OmTfTq1ev489WrV9OrVy9CQkJOuZNz06ZN5ObmYrPZTknclZWVBAcHezRWSehCuFiyLY+xL//Iiwt3MKp7a75/8CLuGNHJ58d8D0iO5fPbh/LC1X3IKSpn8hs/8bf5W6motpsdmtc7fPjw8Zb41q1bmTdvHr179yYmJga73X48qefm5jJ16lT++9//Eh4ezjfffHP8HIWFhcTHxzdpoq2zISUXITi1vPLhrYMZlhZvdlhuZbEorhzQnjE9WvO3+duY8eMuFm/L459X96Fvh+gzn6CFGjt2LK+99hr79++na9euxMXF0bp1awDGjBnD8uXLGTp0KFOmTOGFF16ge/fuPP744zz88MOMGzcOgCVLljBhwgSPx6qMUYfNb+DAgVqWoBNmq6i2M+PHXby+JAurRXHPJWncckGqz7fIG+OHHfk8/NlG8ksq+cv4btw6LNXrbk7aunUr3bt3NzuMBq1du5aXXnqJDz744LTHTZkyhWeffZYuXbo06fz1/fuVUmu01gPrO15a6KJF0lqzeFseT32dwd7CMib0asv0Cd1JjA4xO7Rmc1GXBL69fzgPfbqBp+dtZc3eIzx/VW8igj1bFvAn/fv3Z+TIkdjt9tPeWDR58uQmJ/OzIS100eJk5R3jqa+38uOOfDomhPHUxJ5+V15pCq01by3bxXPfbCc5NpR3bjqPlHjvuNvU21voniYtdCEaUFxWzSvfZ/L+z3sIsVl5/LJ0bhiS3OLvplRKMW14J/q0j+b2D9cw5c0VvH3jQPonxZgdmmiilv1JFi2C3aH56Je9jHxhKe+u2M0153Vg6Z9GcOuw1BafzF0N7hjHF3+8gIjgAH7z1kq+3XLQ7JBEE8mnWfgtrTXfbz3Epa8sY/qXm0lrFc7Xdw/jb1f0Ii48yOzwvFJqfBhf3DGUbm0iuePDNXy5LtvskEQTSMlF+KU1ew/z7IJtrNpzhNT4MN6Y2p/xPdt43SgObxQXHsSs35/Pre+t4sFPNgBwRb/2JkclGkMSuvArWXnHeP6b7XyXcYj48CCentyTa8/rIKWVJgqxWfn3jedx2/ureOCTDWgNU/pLUvd2ktCFX9hbWMpri7P4fG02obYA/jSmC7cMS5XJqM5BiM3K2zcYSf3BTzcQagtgXM82ZofldTIzM4mJiSE+3vyRUvJpFz5td0Epry7O5L/rDxBgUdw0NJU7R3aSGrmb1Cb1699ayb2z1zHz9+czIFlGvwDU1NTw1VdfUVxcTJcuXVi2bBlXXHGFqTHJ36HCJ2XllXD/x+u55IWlzN+Uy01DU1j255E8cXm6JHM3M8ovA2kbFcxt761id0Gp2SF5hdrVh55++mkefPBB2rQ58ddLeXk5F110EXZ7w3PlVFVVMXz4cGpqzm6t2fpIQhc+JSvvGPfMWsfol35gweZcbh2Wyo9/Hsnjl6XTKtKzM9m1ZHHhQfzn5kEopbjp3V85Ulpldkimq6mp4euvv6Z379507dr1pMm4mrqSkbtIQhc+YcehY9w1cy2jX/qRhRmHmHZhR5Y/fDHTJ6TTKkISeXNIiQ/j7RsHkltUwf2frMfhMOcuczMMGTKE3bt3A5CTk8OAAQMICAjgjTfeYOLEidx11108+uijx4//6KOPmDRp0vHnI0eOZOHChQA89thj3H333QBMnjyZjz76yG1xSg1deLVtB4/yr+8zmb/pIGE2K7df1InbhqVKWcUk/ZNieOLydB6bs5nXlmRxzyVpzffmCx6Bg5vce842vWD8s6c9pKEViwCCgoK46aabTjrezJWMJKELr5RxwEjk32w5SHhQAHeO7MRtwzoSE2YzO7QWb+rgJNbsPcJLi3bQLymaC9MSzA7JoxpasaghZq5kJAldeJXNOcW88n0mCzMOEREUwD0Xd+aWYalEh0oi9xZKKZ65oidbDhRz7+z1LLj3Qlo3R//FGVrSnlLfikXTpk1r8PjTrWQUFxfn0ZWMpIYuvMKm7GJue28Vl726nJW7Crn3kjSWP3wxD4zpKsncC4XaAnjztwMoq6rhoc82Ytasrc2hoRWLGmLmSkbSQhemWr+/iH99n8nibXlEhQTywOgu3HRBCpEyJ7fX65QQzvRLu/P4f7fw4S/7+N35yWaH5BGnW7GoIWatZCTzoQtTbMwu4qWFO1iyPZ/o0EBuG5bKjUNTZHEFH6O15sZ3V7Fq92Hm33shqW6eR91X50N310pGTZ0PXUouolltzinmtvdWM/G1n1i7r4iHxnZl+cMXc9fFaZLMfZBSin9c1RtbgIX7P15Pjd1hdkhewXUlo4Z4YiUjKbmIZrHt4FFeXmiMWokMDuBBZ2lFkrjvax0ZzNOTe3L3rHW8uXQndzfnUEYvdsstt5x2v81m44YbbnDre0pCFx6VeegYLy/KZN6mXCKCArj3kjRuGZZKVIgkcn9yeZ9EFmYc4pXvMxnRtRW92keZHVKLJAldeMShoxW8+N0OPl2zn5BAK3eN7MxtF8rwQ3/2v5N68uvuw9z/yXq+vnsYwYEN3/YuPEMSunCrksoaZvywk7eW7abG4eCmoancdXFnYuWGIL8XFRrIP67uze/+/SvPf7OdJy5Pd8t5tdYtcmGSsxmwIglduEWN3cHsVft5edEOCkqquKx3W/48thtJcaFmhyaa0YVpCdw4JJl3ftrNqO6tGNr53OYIDw4OprCwkLi4uBaV1LXWFBYWNvmGo0YldKXUOOAVwAq8rbV+ts7+JOA9INp5zCNa6/lNikT4rBVZBTwxdwtZeSUMSonlrRu60U9WjG+xHhnfnWWZBfzp0w0suG/4OfWXtG/fnuzsbPLz890YoW8IDg6mffumrRJ1xnHoSikrsAMYDWQDq4DrtdYZLsfMANZprd9USqUD87XWKac7r4xD932Hjlbw9LytfLXhAMlxoUy/tDuj01u3qJaUqN/6/UVc+eYKJvVJ5MVr+5odjl853Tj0xrTQBwFZWutdzpPNBiYBGS7HaCDS+TgKOHD24QpvV2N38J8Ve3h5USZVdgf3jUrj9os6SSeYOK5vh2juHNmZf32fyej01ozv1dbskFqExiT0dsB+l+fZwOA6xzwJfKeUuhsIA0bVdyKl1DRgGkBSUlJTYxVeYPvBYzz46Xo25xxlZNcEnpzYg+Q4994dKPzD3Rd3Zun2PB79chMDkmNkAZJm4K47Ra8H/qO1bg9cCnyglDrl3FrrGVrrgVrrgQkJ/j3lpr+psTt4fUkWl7+6nNyiCt6c2p93bjpPkrloUKDVwovX9KWsys7Dn/v3BF7eojEJPQfo4PK8vXObq1uBTwC01j8DwYD5S2ALt8jKO8aVb67gH99uZ1R6K767fzjje7WVWrk4o86twnlkfDeWbM9n1q/7z/wCcU4ak9BXAWlKqVSllA24Dphb55h9wCUASqnuGAm95XVL+xmtNR/9spdL/7WcvYfLePX6frwxdYCsFiSa5MYhKVzQOY6n52Wwt1AWmPakMyZ0rXUNcBfwLbAV+ERrvUUp9ZRSaqLzsAeB3yulNgCzgJu0/H3l045VVHP3rHVM/3Izg1Nj+e7+4VzeJ9HssIQPslgU/7iqD1aL4oFPNmBvQWuRNjeZPlecYsehY0x7fzX7Dpfx4Jiu3HFRJywWKa+IczNnXQ73fbyeh8Z25c6Rnc0Ox2ed67BF0YIsyjjEfR+vJzjQyuxpQxiUGmt2SMJPTOprTOD18qIdXNQlgZ7tZAIvd5P50AVg1MvfWJrF7z9YTWp8GF/dfYEkc+FWSimentyTmFAb93+8norqhucKF2dHErrA7tA8Nmczz3+znct6J/LJH4bQNirE7LCEH4oJs/HcVb3JzCvhn99uNzscvyMJvYWrrLFz96y1fPTLPu4Y0Yl/XdeXEJvc8Sk8Z2TXVkwdnMS/f9rNzzsLzQ7Hr0hCb8FKK2u4+d1VzN90kMcmdOfhcd1kbLloFtMndCc5NpQ/fbqBoxXVZofjNySht1DlVXZufW8Vv+w+zIvX9OG2CzuaHZJoQUJtAbx4bV9yi8v569yMM79ANIok9BaootrOtA9WH0/mU/o3bYpOIdyhf1IMfxzRmc/XZvPN5oNmh+MXJKG3MFU1Du74cA3LMgt4/sreTOrbzuyQRAt2zyVp9EiM5NEvN5F3rMLscHyeJPQWRGvN9C83sWR7Ps9c0ZOrB3Y484uE8CBbgIWXr+1LSWUNf/l8k0zgdY4kobcgbyzdyadrsrnnkjSmDk42OxwhAEhrHcGfx3bl+215fLxKJvA6F5LQW4h5G3P5x7fbmdgnkftHpZkdjhAnueWCVIZ0jON/v84gp6jc7HB8liT0FmBTdjEPfLKeAckxPH9VbxmaKLyOxaJ4/qreODQ8+oWUXs6WJHQ/V1xezR9nriE2zMaM3w2QZeKE1+oQG8qfx3Xlhx35fLG27pILojEkofsxrTUPfbqB3KIKXvtNf5nHXHi9G4ekMDA5hqe+zpBRL2dBErof+/fy3XyXcYhHxndjQHKM2eEIcUYWi+K5q3pTXm3niTlbzA7H50hC91Nr9h7m2QXbGJPemluHpZodjhCN1ikhnPtGpfHNloPM35Rrdjg+RRK6HzpcWsVdM9eRGB3CP67uI52gwudMu7AjvdpF8cR/N3OktMrscHyGJHQ/43Bo7vt4PYUlVbwxtT9RIYFmhyREkwVYLTx3ZW+Kyqp56muZ66WxJKH7mTeWZvHjjnyeuDxdVoQRPi09MZI/jujEl+tyWLItz+xwfIIkdD+yYmcBLy7cwcQ+iUwdnGR2OEKcszsv7kyX1uE8+uUmSiprzA7H60lC9xN5Ryu4Z9Z6UuPD+PuUXlI3F34hKMDK36f0Jre4gpcX7jA7HK8nCd0P2B2ae2evp6SymjemDiAsSNb+Fv5jQHIM1w9K4t0Ve8g4cNTscLyaJHQ/8MaSLH7eVchTk3rStU2E2eEI4XYPj+tKdEggj83ZhMMh0wI0RBK6j1u95zAvf5/JxD6JXD1AFqoQ/ik61Majl3Zn7b4iPl4tMzI2RBK6Dysuq+be2etpFx3CM1f0lLq58GtT+rdjcGoszy7YRkFJpdnheCVJ6D5Ka80jX2zk0NEK/nV9PyKCZby58G9KKZ65oidlVTX8bf5Ws8PxSpLQfdSsX/ezYPNB/jS2K307RJsdjhDNonOrCKYN78gXa3NYuavQ7HC8jiR0H7Tj0DH++tUWLkyLZ9qFHc0OR4hmddfINNpFh/Dk3C3U2B1mh+NVJKH7mIpqO3fPXEdEcAAvXNMHi0Xq5qJlCbFZmT6hO9sOHmO2LFl3EknoPuaZeVvZfugYL1zTl1YRwWaHI4Qpxvdsw+DUWF74bjvFZdVmh+M1JKH7kCXb8/hg5V5+f2EqF3VJMDscIUyjlOJ/Lu9BcXk1Ly2SO0hrNSqhK6XGKaW2K6WylFKPNHDMNUqpDKXUFqXUTPeGKY6UVvHnzzbStXUED47panY4QpguPTGS6wcl8cHKvWQeOmZ2OF7hjAldKWUFXgfGA+nA9Uqp9DrHpAF/AS7QWvcA7vNArC2W1prH5mymqKyKF6/tI+uCCuH04JiuhNqs/H3BNrND8QqNaaEPArK01ru01lXAbGBSnWN+D7yutT4CoLWWuS7daO6GA8zblMt9o7rQI1GmxBWiVmyYjT+O6MzibXkyjJHGJfR2gGtXcrZzm6suQBel1E9KqZVKqXH1nUgpNU0ptVoptTo/P//sIm5hcovLeXzOZvonRfOH4TJEUYi6br4ghbZRwfx9wTa0btnzvLirUzQASANGANcDbymlTrnbRWs9Q2s9UGs9MCFBOvXOxOHQPPTpRqrtmhev6UuAVfqwhagrONDK/aO7sGF/EfM3HTQ7HFM1JkPkAB1cnrd3bnOVDczVWldrrXcDOzASvDgHH/26j+VZBUyf0J2U+DCzwxHCa13Zvz1dW0fw/LfbqKppuTcbNSahrwLSlFKpSikbcB0wt84xczBa5yil4jFKMLvcGGeLk1tcznMLtnFB5zhZfUiIM7BaFA+P78rewjJmr9pndjimOWNC11rXAHcB3wJbgU+01luUUk8ppSY6D/sWKFRKZQBLgIe01tJDcZa01jw+ZzM1Dgd/v6K3zKIoRCOM7NqK81JieH1JFhXVdrPDMUWjirJa6/la6y5a605a62ec257QWs91PtZa6we01ula615a69meDNrfzduUy6KteTw4uitJcaFmhyOET1BKcf+oLhw6WsnsX1tmK1162bxMUVkVT87dQq92Udx8QYrZ4QjhU4Z0imNQaixvLN3ZIlvpktC9zDPztnKkrJpnr+wlo1qEaCKlFPeNSiPvWCWzWmArXTKGF1meWcCna7L5w/COcgOREGdpaKd4BqfG8mYLbKVLQvcSFdV2ps/ZRGp8GPdcIiM+hTgX943q0iJb6ZLQvcSMH3ext7CMpyf3lLlahDhHQzrFMSgllreX7W5Ri2BIQvcC2UfKeH1JFhN6teWCzvFmhyOEX5g2vCM5ReXM25RrdijNRhK6F3hm3lYsSvHohO5mhyKE37i4Wys6JoTx1rJdLWaOF0noJluWmc+CzQe5c2Qn2kWHmB2OEH7DYlH8/sKObM45ys8tZCZGSegmqrY7eHLuFpLjQrlNFnsWwu2u6NeO+HAbM35sGTORSEI30cer9rMzv5Tpl3aXjlAhPCA40MqNQ1JYuj2f7Qf9f1UjSegmKams4eVFOxiUEsvo9NZmhyOE3/rt+cmEBFp5Z/lus0PxOEnoJnnrx10UlFTxl0u7yeRbQnhQTJiNiX0SmbvhAMcqqs0Ox6MkoZsg71gFby3bxaW92tAvKcbscITwe78ZnER5tZ3/rj9gdigeJQndBK8syqSqxsFDY7uZHYoQLULv9lGkt41k5i/7/HoIoyT0ZravsIzZq/bzm8FJpMoqREI0C6UU1w9OIiP3KBuzi80Ox2MkoTezN5ZmYbUo7hzZ2exQhGhRJvdNJCTQysxf/Hd+F0nozSj7SBmfrcnm+vM60Doy2OxwhGhRIoID/b5zVBJ6M3pz6U4sSnH7iE5mhyJEi1TbOTrHTztHJaE3kwNF5Xyyej9XD2xP2yi5xV8IM/RuH0W3NhF8vibb7FA8QhJ6M/l/P+xEa7hDWudCmEYpxRX92rF+fxF7CkrNDsftJKE3g4KSSmat2s+U/u1oHyOLPgthpol9E1EKvxyTLgm9Gcz6ZR9VNQ6mDZfWuRBmaxsVwuDUWOasz/G7MemS0D2sqsbBByv3MrxLAp1bhZsdjhACmNy3HbsLStly4KjZobiVJHQPW7A5l7xjldw8NMXsUIQQTmN6tMFqUXyz+aDZobiVJHQPe/enPaTGh3FRlwSzQxFCOMWG2RicGsuCzf61PJ0kdA9at+8I6/cXceOQZCwWmVFRCG8yrmcbduaXkpXnP/OkS0L3oP+s2EN4UABXDmhvdihCiDrG9mgDwIJN/lN2kYTuIQUllczbmMtVA9oTERxodjhCiDpaRwbTLymaRdvyzA7FbSShe8iXa3OocWimDk4yOxQhRANGdm3FxuwiCkoqzQ7FLSShe4DWmk9W76dfUjRprSPMDkcI0YCRXVuhNfy4I9/sUNxCEroHbMguJjOvhKsHdDA7FCHEafRIjCQ+PIil21tQQldKjVNKbVdKZSmlHjnNcVcqpbRSaqD7QvQ9n6/JJjjQwmV92podihDiNCwWxUVdEvhhRz52h+/fNXrGhK6UsgKvA+OBdOB6pVR6PcdFAPcCv7g7SF9SY3cwb1Muo7q3JlI6Q4XweiO6JlBcXs3G7CKzQzlnAY04ZhCQpbXeBaCUmg1MAjLqHPe/wHPAQ26N0Mf8tLOQw6VVXN4n0exQhBCNMKRTHAArdx2mX9tg2PEt7P8FjuZAVRlUl4Gj5sQLel0Ng35vUrSn15iSSztgv8vzbOe245RS/YEOWut5pzuRUmqaUmq1Ump1fr5/1Kzq+mrDASKCAuTOUCF8RHx4EF0SQolY/xa83As+vRFWvwt526CsAFAQGGJ85W+HzZ+bHXKDGtNCPy2llAV4EbjpTMdqrWcAMwAGDhzo+wWrOipr7Hy7+SBjerQhONBqdjhCiMaoPMYr6p90L1qOI3UElgvvh+QLwFpPyfSDK6DSe+8sbUwLPQdwHa7R3rmtVgTQE1iqlNoDnA/MbYkdoz/vLORYZQ0TercxOxQhRGNUl8PM6+h67GeerL6BDSPehY4j6k/mAJZAsHvveqSNSeirgDSlVKpSygZcB8yt3am1LtZax2utU7TWKcBKYKLWerVHIvZii7YeIiTQytBO8WaHIoQ4E4cDPr8N9v5EyaWv8x/7OFbuPnL611gDT66ne5kzJnStdQ1wF/AtsBX4RGu9RSn1lFJqoqcD9BVaaxZl5DG8S7yUW4TwBT+/Ctu+hrHPEHne9XRuFc6qPYdP/xqL1atb6I2qoWut5wPz62x7ooFjR5x7WL5ny4GjHDxawYPdu5gdihDiTPavgkV/hfRJcP4fAejbIZql2/PQWqNUA7OjWgLB4b0JXe4UdZOFGYewKLi4WyuzQxFCnE5VGXz5B4hsBxNfBWfy7tM+ioKSKnKKyht+rTUQ7D5cchGNsywznz4dookLDzI7FCHE6Sx+Gg7vhEmvQXDU8c19OkQDsGF/ccOvtQT4dg1dnNmximo2ZBdzgXSGCuHdctbAyjdg4C3Q8aKTdnVrE4nNajn9HaNWKbn4vV92Hcbu0FzQWRK6EF7LXgNf3QvhrWHUk6fstgVY6J4Yyfr9p0nofjBsUZzBTzsLCA600D852uxQhBAN+eX/4OAmGP/cSaUWV33bR7Epp7jhibqk5OL/VmQVcl5KLEEBMlxRCK9UtB+W/A3SxhojWxrQu300ZVV2duWX1H+ANUBa6P6soKSS7YeOyc1EQngrrWH+Q4CGCf88PqqlPl3bGAvS7DjUQEKXYYv+bd0+o952XkqMyZEIIeq1dS7sWAAj/gLRp18SslNCOEpBZl4D87VYA0E7jLtMvZAk9HO0fv8RAiyKnu3qry2cqQMAABTpSURBVMkJIUxUXgTz/wxteh+/geh0QmxWOsSEkpnXUAvdeS+ml9bRz3m2xZZu3b4iurWNkNv9hfBGi56E0jz4zWyj/t0Iaa3CyWqo5FI7aZejGrC5JUR3khb6ObA7NBuzi+nbQUa3COF19v4Ma941WuaJ/Rr9ss6tw9lVUEKNvZ6yisWZ0L20Y1QS+jnYmV9CSWUNfTtI/VwIr1JTaYw5j0oyaudNkNYqgmq7Zu/hslN3ennJRRL6OVjv7BCVFroQXmb5y1CwHSa8AEHhTXppWivj+Mz6yi61ZRtpofuf9dlFRAQH0DE+zOxQhBC18nfAsn9Czyuhy5gmv7yTM6Fn1TfSxeJaQ/c+ktDPQeahY3RrE4HF0vC4ViFEM7LXwJzbITAUxj17VqcIDwqgXXRI/SNdjneKSsnF72TlldC5VdP+nBNCeNBPLxkTcE14AcLPfirrpNhQ9p+uhu6lU+hKQj9LhSWVHCmrplOCJHQhvELuBlj6LPSYAr2uOqdTtY0OJre44tQdVim5+KUs559j0kIXwgvUVMKXt0NovNE6P0eJUSEcOlpx6tBFi3SK+qWsfEnoQniNRX+FvAxjBaLQ2HM+XWJ0CA4NeccqT95hkRq6X8rKKyEk0EpiVIjZoQjRsm2bDytfh0HTzmpUS33aRgcDkFtcZzk6Gbbon7LySujUKkxGuAhhpqJ9MOcOaNsHxjztttPWNtRyiurU0WXYon/amVdCZ+kQFcI89mr47BZw2OGqdyHAfev5Hm+h110w2suHLcrkXGehtLKGA8UVUj8XwkwLn4DsVUYyj+vk1lNHBgcSHhRw6kiX43O5eGdClxb6Wcg+YvzWTo6TO0SFMMW6D43FngffDj2neOQtEqODOXBKC712LhcpufiNghKj5zshwn1/4gkhGmnfL/D1/dBxBIx5xmNv0zYqpJ4WunSK+p3ahB4fLgldiGZVnA0fT4Wo9kappZFznJ+NxOjgU0e5yLBF/5PvHJuaIAldiOZTfgQ+vAqqK+D62W4Zb346baNCKCipoqLafmKjDFv0PwUlVdisFiJDpE9ZiGZRVQYzr4PDO+G6jyChq8ffsm2UMdLloGvZRYYt+p+Ckkriwm2o06weLoRwE3sNfHYz7P8FprwFHS9qlret7SMrLK06sVGGLfqf2oQuhPAwh924cWjHN8YcLT0mN9tbh9qM9HhSycXLhy1KQj8LhSVV0iEqhKfZa+DLP8Dmz+Dix+G825r17UOcC7+XV9VTQ/flkotSapxSartSKksp9Ug9+x9QSmUopTYqpb5XSiW7P1TvUVBSKQldCE+yV8PntxrJfNSTMPxPzR5CiM1Ij2UntdB9vFNUKWUFXgfGA+nA9Uqp9DqHrQMGaq17A58Bz7s7UG+htZYWuhCeVFUGH/8OMuYY87MMu9+UMIKdLfSKqnpKLj7cQh8EZGmtd2mtq4DZwCTXA7TWS7TWtct7rATauzdM73G0vIYqu4N4qaEL4X4l+fDeZUbN/NJ/wtC7TQultoZeftKwRd+vobcD9rs8zwYGn+b4W4EF9e1QSk0DpgEkJSU1MkTvki93iQrhGYU74cMr4VguXPshdL/M1HCO19BdE7pSoKwtY5SLUuq3wECg3nFFWusZwAyAgQMHane+d3ORu0SF8IDtC+CLP4DFCjd+DR3OMzsiggKcNXTXkgsYrXQvLbk0JqHnAB1cnrd3bjuJUmoUMB24SGtdWXe/vygsMcakyrBFIdzAYYelf4cf/2HMaX7N+xCTYnZUAFgsiuBAy8nDFsGoo/twyWUVkKaUSsVI5NcBv3E9QCnVD/h/wDitdZ7bo/Qi0kIXwk2K9htjzPcsg36/hUtfgMBgs6M6Sagt4ORhi2D8FeGrLXStdY1S6i7gW8AKvKO13qKUegpYrbWeC/wDCAc+dd49uU9rPdGDcZumoKQSi4KYUGmhC3FWtIYNs2DBw6AdMPE16P87s6OqV0ig9eQaOhglFy8dttioGrrWej4wv862J1wej3JzXF6roKSS2LAgrLL0nBBNd3i3kcgzv4WkoXDFm15TYqlPcKClnha6b9fQhYsjpdXEhgWaHYYQvqW6HJa/DMtfMlq4Y/9mLE5hsZod2WmF2gLqaaEH+HQNXbgoraohLEgumxCNUlMF6z6AZS/A0RzoeRWM+V+ITDQ7skYJCbQ20EKXhO4XyqrshNnksglxWlVlsHG20SIv2gcdBsOUGZAyzOzImiTYZqW4vE55xceHLQoXpZU1xIVJh6gQ9TqyF1b/G9a8BxVFkNgfJrwEnS8xbsrxMSGBFg4V+9ewReFCSi5C1FFaaMy7sulT2PczKAt0v9yokScN8clEXqveGrovD1sUJyurtBNq8+6OHCE8Sms4uBF2fAeZ30H2KkBDQjdjmtve10J0hzOexhcE++OwRXFCaVUN4dJCFy1N5THYucRI4JkLoeQgoKBdfxjxCHSbAK17+nRrvD7SKerH7A5NRbXj+CxsQvgtraEg05nAv4W9PxtlhqAo6HwxpI2FzqMgPMHsSD0qxGahvNqO1vrEkpPWAGmh+4PSKuO3cliQlFyEH6ouhz3LnUn8Oziyx9jeKh2G3AlpY4zRKtaWkzZCbQHYHZpqu8YW4EzolkBjFI8Xajk/GTcoqzT+9JIWuvALWkNeBuxcDFnfw94VYK+EwFBIvQiG3mMkcT+ph5+NYJcpdG3O2Rdl2KKfkBa68HmlhbBriZHEdy425h4HSOhurNnZ+WJIHuZ1k2SZxXVd0agQ5x3iFrlT1C/UttDlxiLhM2qqjFEoO783WuG5GwANITHQcaQxPrzjSIhqZ3akXql2RFt53XVFpYXu+2pb6KHSQhferHDniRb47h+hqsRYZafDIBg53WiFt+3r9fOoeINglxb6cTJs0T+UVjpLLtJCF96kJM9I3LuWwu4fjFvtwZjFsPe10OliSB0OwZFmRumTQuptocuwRb9Q6vwtLTV0YaqKo7D3J9j1g5HA8zKM7cFRkHKh0ZnZ6WKI62RunH4gpN4WeoAkdH9Q5myhyygX0ayqKyD71xMJPGctaDsEBBu31ve+xhiV0raPlFHcrP4aupRc/MLxFrokdOFJDjvkrj+RwPethJoKow7ebgBc+ICRwDsMggBZCtGTXIctHifDFv3D8Ra6lFyEO2kNBTtOJPA9y6Ci2NjXKh0G3AwdR0DyUKmDN7PjNfQqlxKLDFv0DyVVNdgCLARaLWaHInyZ1pC/3Ujce5Yb9fDSfGNfdBJ0n2gk8NThEN7KzEhbvHpr6DJs0T+UVdoJk5kWRVNpDfnbjOS9Zxns+QnKCox9ke2g0yXGwg8pwyA21dxYxUlO1NAdJzbKsEX/UFpVIx2i4swcjpMT+N4VLgm8PaSNhuQLjAQek+J3MxT6kyDn7f6ndIqijb4OL+uEluzUBGWVdhmyKE7lcED+VmcCd5ZQygqNfVEdjARe2wKPTpYE7kOUUs4pdF1q5rWTkzlqJKH7MlmtSAB1ErizhFJ+2NgXlWRMLVubwGOSzY1VnLMQm7WeFjpG2cXLRhlJdmqC0soaGbLYEtlr4OAGo3Sy92fYtwLKjxj7opOg63gjeSdfIAncDxkt9Do1dPDKjlHJTk1QVmUnPty7fiMLD6guh5w1zgS+Avb/CtWlxr7YjtB1grMFfoGR0IVfC7FZqag7ORd45dBFSehNICUXP3XsIBxYZyTuvSvgwFqwVwEKWveAflONOzKTh0JEG7OjFc0sJNB6fGI+4ERClxa6b5MFov1AaYExheyBtXBgvXEb/bEDxj5LACT2g/PvgKShkDTYmGZWtGjtokPYfujYiQ1Wlxq6l5GE3gSyQLQPKS8y1sTMy4C8rZC3xfheewMPQGwno2yS2N9Y7LhNb7CFmhez8ErpiZF8m3GQkkrn///aTlEvnKBLslMj1dgdskC0N9Ha6Jg8vBsO73J+7TzxuHbYIBhLqiV0M0aftE43bqdP7Ach0ebFL3xGj8RItIZtuUcZmBJ78rBFLyPZqZHKqmXq3GajtTGXydEDcDQHirNPPD6aA8U5xvPajspake2NOy27X250XsZ1NpJ3dDJYZLoGcXbSE435czJqE7pFSi4+TxaIPgc1VUaCrigySiHlR4w7J0vznV/1PLZXnXwOZYHwNhCZaLSy00Ybt83HdjS+YpIhMMScf5/wa20ig4kJDSTjwFFjg693iiqlxgGvAFbgba31s3X2BwHvAwOAQuBarfUe94ZqrsLSSqAFtdAdDqgug6pSYwmzqlKXrxKXbS77KoqNhF1RdPLj6rKG3ycgGMJaQVgchLeG1j0hLB7CEoyEHdnOWO8yvPWJzighmpFSivTESDJynQn9eKeoD5ZclFJW4HVgNJANrFJKzdVaZ7gcditwRGvdWSl1HfAccK0nAjZDVY2Dx+dsJtRmpX+Sh0c9OBzGb357lfEnnb2qgcf1HVNljKGuqYSa2u8VxgIJNRUnb6922V/jsr+63PlVeuZYaykL2MKNFXOCo43adGzHE4+Do419tY9Doo2EHZYAtjC5FV54vfS2kbz3816q7Q4CfbyFPgjI0lrvAlBKzQYmAa4JfRLwpPPxZ8BrSimltdZujBWAVV+8QqvNb6FwPbXx2DUtqOPbnMe5hKLqHAOn2wdaa17VmtiwQEL+Yz3pXJxyfpd9jdmmnZP81CZk7XIDgztYAoxWcEAQBIQ4vwdDYLDxPTgSAlob2wNDThwXFG4kW1uYkaxPehx+8vaAIEnKwq/1SIyiqsbBmJd+pJ8jgxeBg/+5kQp1djcaFg64jwETbnNvkDQuobcD9rs8zwYGN3SM1rpGKVUMxAEFrgcppaYB0wCSks7uDruAiHgKQjvWnvH4du18fPJvEOc2peo5vs4x9ZzLVUJkCCEJ4S6ncDmm7vnVyb9azrjNEmj8GWcNBKvN5Xudx5aA+rfXvvZ44g4+8WWVmr8Q52pE1wSu7N/euGPU0Yef88cR5Cg/6/PZwmPdGN0Jzfq/XWs9A5gBMHDgwLNqvfcbPRVGT3VrXEIIcTrRoTZeuKaPy5ZhpsVyOo0Zy5UDdHB53t65rd5jlFIBQBRG56gQQohm0piEvgpIU0qlKqVswHXA3DrHzAVudD6+Cljsifq5EEKIhp2x5OKsid8FfIsxbPEdrfUWpdRTwGqt9Vzg38AHSqks4DBG0hdCCNGMGlVD11rPB+bX2faEy+MK4Gr3hiaEEKIp5H5oIYTwE5LQhRDCT0hCF0IIPyEJXQgh/IQya3ShUiof2HsOp4inzp2oXkLiahpvjMsbYwKJq6n8Na5krXVCfTtMS+jnSim1Wms90Ow46pK4msYb4/LGmEDiaqqWGJeUXIQQwk9IQhdCCD/hywl9htkBNEDiahpvjMsbYwKJq6laXFw+W0MXQghxMl9uoQshhHAhCV0IIfyEVyZ0pdQ4pdR2pVSWUuqRevYHKaU+du7/RSmV4rLvL87t25VSY5sxpgeUUhlKqY1Kqe+VUsku++xKqfXOr7pTD3s6rpuUUvku73+by74blVKZzq8b677Ww3G95BLTDqVUkcs+T16vd5RSeUqpzQ3sV0qpfznj3qiU6u+yzyPXqxExTXXGskkptUIp1cdl3x7n9vVKqdXuiqmRcY1QShW7/KyecNl32p+/h+N6yCWmzc7PU6xznyevVwel1BJnHtiilLq3nmM8+/nSWnvVF8YUvTuBjoAN2ACk1znmj8D/OR9fB3zsfJzuPD4ISHWex9pMMY0EQp2P76iNyfm8xMRrdRPwWj2vjQV2Ob/HOB/HNFdcdY6/G2NaZo9eL+e5hwP9gc0N7L8UWICxTuD5wC/NcL3OFNPQ2vcCxtfG5Hy+B4g36VqNAL4+15+/u+Oqc+zlGOszNMf1agv0dz6OAHbU8//Ro58vb2yhH1+UWmtdBdQuSu1qEvCe8/FnwCVKKeXcPltrXam13g1kOc/n8Zi01ku01mXOpysxVnbytMZcq4aMBRZqrQ9rrY8AC4FxJsV1PTDLTe99WlrrHzHm7G/IJOB9bVgJRCul2uLB63WmmLTWK5zvCc332WrMtWrIuXwu3R1Xc362crXWa52PjwFbMdZbduXRz5c3JvT6FqWue1FOWpQaqF2UujGv9VRMrm7F+C1cK1gptVoptVIpNdkN8TQ1riudf959ppSqXU7QU9eqSed2lqZSgcUumz11vRqjodg9eb2aou5nSwPfKaXWKGMR9uY2RCm1QSm1QCnVw7nNK66VUioUIyl+7rK5Wa6XMsrA/YBf6uzy6OdLloR3M6XUb4GBwEUum5O11jlKqY7AYqXUJq31zmYK6Stglta6Uin1B4y/bC5upvdujOuAz7TWdpdtZl4vr6WUGomR0F1XKB7mvFatgIVKqW3OFmxzWIvxsypRSl0KzAHSmum9G+Ny4CettWtr3uPXSykVjvFL5D6t9VF3nvtMvLGFfi6LUjfmtZ6KCaXUKGA6MFFrXVm7XWud4/y+C1iK8ZvbHc4Yl9a60CWWt4EBjX2tJ+NycR11/iT24PVqjIZi9+T1OiOlVG+Mn98krfXxBdhdrlUe8CXuKTE2itb6qNa6xPl4PhColIrH5Gvl4nSfLY9cL6VUIEYy/0hr/UU9h3j28+WJzoFz7FgIwOgQSOVEh0qPOsfcycmdop84H/fg5E7RXbinU7QxMfXD6AhKq7M9BghyPo4HMnFTB1Ej42rr8vgKYKU+0Qmz2xlfjPNxbHPF5TyuG0YnlWqO6+XyHik03NE3gZM7rX719PVqRExJGP1BQ+tsDwMiXB6vAMY147VqU/uzw0iM+5zXrVE/f0/F5dwfhVFnD2uu6+X8t78PvHyaYzz6+XLbBXbzD+tSjB7incB057anMFq+AMHAp84P+a9AR5fXTne+bjswvhljWgQcAtY7v+Y6tw8FNjk/1JuAW5v5Wv0d2OJ8/yVAN5fX3uK8hlnAzc0Zl/P5k8CzdV7n6es1C8gFqjHqlLcCtwO3O/cr4HVn3JuAgZ6+Xo2I6W3giMtna7Vze0fnddrg/BlPb+ZrdZfLZ2slLr9w6vv5N1dczmNuwhgg4fo6T1+vYRg1+o0uP6tLm/PzJbf+CyGEn/DGGroQQoizIAldCCH8hCR0IYTwE5LQhRDCT0hCF0IIPyEJXQgh/IQkdCGE8BP/H7foUGwBS0lxAAAAAElFTkSuQmCC\n", "text/plain": [ "

" ] }, "metadata": { "needs_background": "light", "tags": [] }, "output_type": "display_data" } ], "source": [ "plt.plot(xs3,sol[:,numsteps],label=r'$u(x)$')\n", "plt.plot(xs2,sol_adj[:,0],label=r'$u^{\\ast}(x)$')\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "igOTigTMeFBq", "nbpages": { "level": 3, "link": "[6.1.3.5 Compare sensitivities to parameters resulted from finite differences and the adjoint methodology: $\\rho$, $\\omega$,$\\kappa_l$,$\\kappa_h$,$q$ respectively](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.5-Compare-sensitivities-to-parameters-resulted-from-finite-differences-and-the-adjoint-methodology:-$\\rho$,-$\\omega$,$\\kappa_l$,$\\kappa_h$,$q$-respectively)", "section": "6.1.3.5 Compare sensitivities to parameters resulted from finite differences and the adjoint methodology: $\\rho$, $\\omega$,$\\kappa_l$,$\\kappa_h$,$q$ respectively" } }, "source": [ "### 6.1.3.5 Compare sensitivities to parameters resulted from finite differences and the adjoint methodology: $\\rho$, $\\omega$,$\\kappa_l$,$\\kappa_h$,$q$ respectively" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "executionInfo": { "elapsed": 3333, "status": "ok", "timestamp": 1593494147396, "user": { "displayName": "黄尚晖", "photoUrl": "", "userId": "16496223364462243774" }, "user_tz": 240 }, "id": "ggbDJFrh6I9-", "nbpages": { "level": 3, "link": "[6.1.3.5 Compare sensitivities to parameters resulted from finite differences and the adjoint methodology: $\\rho$, $\\omega$,$\\kappa_l$,$\\kappa_h$,$q$ respectively](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.5-Compare-sensitivities-to-parameters-resulted-from-finite-differences-and-the-adjoint-methodology:-$\\rho$,-$\\omega$,$\\kappa_l$,$\\kappa_h$,$q$-respectively)", "section": "6.1.3.5 Compare sensitivities to parameters resulted from finite differences and the adjoint methodology: $\\rho$, $\\omega$,$\\kappa_l$,$\\kappa_h$,$q$ respectively" }, "outputId": "400717a8-bffc-462f-9ffd-d9dd73760fc4" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-0.09947973321489112 -0.09948429352239387\n" ] } ], "source": [ "sens_adj = np.zeros(5)\n", "sens_for = np.zeros(5)\n", "#rho sens\n", "delta = 1.0e-6\n", "sol_pert,Q_v,xs2 = run_calc(Lx, Nx, rho+delta,Source_func(xs,1), omega_nom, v_nom, \n", " kappa_func(xs,kappal_nom, kappah_nom), dt, numsteps,1.8,tfinal,1.5,1.9,IC)\n", "sens_for[0] = (Q_v-Q)/delta\n", "\n", "dudx = np.zeros((Nx+2,numsteps+1))\n", "dudx[:,0:(numsteps)] = (sol[:,1:(numsteps+1)] - sol[:,0:(numsteps)])/(dt)\n", "#dudx[Nx-1] = (0-sol[Nx-2])/(2*dx)\n", "#dudx[0,:] = (sol[0,:]-0)/(dx)\n", "sens_adj[0] = np.sum(-dudx*sol_adj*dx*dt)\n", "print(sens_for[0],sens_adj[0])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "executionInfo": { "elapsed": 4014, "status": "ok", "timestamp": 1593494539024, "user": { "displayName": "黄尚晖", "photoUrl": "", "userId": "16496223364462243774" }, "user_tz": 240 }, "id": "vK0MHKop6I-B", "nbpages": { "level": 3, "link": "[6.1.3.5 Compare sensitivities to parameters resulted from finite differences and the adjoint methodology: $\\rho$, $\\omega$,$\\kappa_l$,$\\kappa_h$,$q$ respectively](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.5-Compare-sensitivities-to-parameters-resulted-from-finite-differences-and-the-adjoint-methodology:-$\\rho$,-$\\omega$,$\\kappa_l$,$\\kappa_h$,$q$-respectively)", "section": "6.1.3.5 Compare sensitivities to parameters resulted from finite differences and the adjoint methodology: $\\rho$, $\\omega$,$\\kappa_l$,$\\kappa_h$,$q$ respectively" }, "outputId": "483944ee-a8ea-46d8-bbf6-70f962609667" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.28897547187101136 0.288993704223791\n" ] } ], "source": [ "#omega sens\n", "delta = 1.0e-6\n", "sol_pert,Q_v,xs2 = run_calc(Lx, Nx,rho, Source_func(xs,1), omega_nom+delta, v_nom, \n", " kappa_func(xs,kappal_nom, kappah_nom), dt, numsteps, 1.8,tfinal,1.5,1.9,IC)\n", "sens_for[1] = (Q_v-Q)/delta\n", "\n", "dudx2 = np.zeros((Nx+2,numsteps+1))\n", "dx = Lx/Nx\n", "dudx2[1:(Nx+1),:] = (sol[0:Nx,:]**4 - 2*sol[1:(Nx+1),:]**4 + sol[2:(Nx+2),:]**4)/(dx**2)\n", "sens_adj[1] = np.sum(dudx2*sol_adj*dx*dt)\n", "print(sens_for[1],sens_adj[1])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "executionInfo": { "elapsed": 3386, "status": "ok", "timestamp": 1593495888556, "user": { "displayName": "黄尚晖", "photoUrl": "", "userId": "16496223364462243774" }, "user_tz": 240 }, "id": "QPbNZ1RD6I-D", "nbpages": { "level": 3, "link": "[6.1.3.5 Compare sensitivities to parameters resulted from finite differences and the adjoint methodology: $\\rho$, $\\omega$,$\\kappa_l$,$\\kappa_h$,$q$ respectively](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.5-Compare-sensitivities-to-parameters-resulted-from-finite-differences-and-the-adjoint-methodology:-$\\rho$,-$\\omega$,$\\kappa_l$,$\\kappa_h$,$q$-respectively)", "section": "6.1.3.5 Compare sensitivities to parameters resulted from finite differences and the adjoint methodology: $\\rho$, $\\omega$,$\\kappa_l$,$\\kappa_h$,$q$ respectively" }, "outputId": "c5ba1397-1180-438b-9776-f1945666236d" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-0.03022436832755826 -0.030226185976361303\n" ] } ], "source": [ "low_area = lambda x: (x>0.5*Lx)*(x<0.75*Lx)\n", "low_area_abs = lambda x: (x>1.5)*(x<1.9)*(x>0.5*Lx)*(x<0.75*Lx)\n", "delta = 1e-6\n", "sens_adj[2] = np.sum(dx*(sol[low_area_abs(xs2),:])*dt)\n", "sens_adj[2] += np.sum(-sol[low_area(xs2),:]**4*sol_adj[low_area(xs2),:]*dt*dx)\n", "sol_pert,Q_v,xs2 = run_calc(Lx, Nx,rho, Source_func(xs,1), omega_nom, v_nom, \n", " kappa_func(xs,kappal_nom+delta, kappah_nom), dt, numsteps, 1.8,tfinal,1.5,1.9,IC)\n", "sens_for[2] = (Q_v-Q)/delta\n", "print(sens_for[2],sens_adj[2])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "executionInfo": { "elapsed": 3403, "status": "ok", "timestamp": 1593495893175, "user": { "displayName": "黄尚晖", "photoUrl": "", "userId": "16496223364462243774" }, "user_tz": 240 }, "id": "PSCchuzT6I-G", "nbpages": { "level": 3, "link": "[6.1.3.5 Compare sensitivities to parameters resulted from finite differences and the adjoint methodology: $\\rho$, $\\omega$,$\\kappa_l$,$\\kappa_h$,$q$ respectively](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.5-Compare-sensitivities-to-parameters-resulted-from-finite-differences-and-the-adjoint-methodology:-$\\rho$,-$\\omega$,$\\kappa_l$,$\\kappa_h$,$q$-respectively)", "section": "6.1.3.5 Compare sensitivities to parameters resulted from finite differences and the adjoint methodology: $\\rho$, $\\omega$,$\\kappa_l$,$\\kappa_h$,$q$ respectively" }, "outputId": "dd8fcc29-1b80-4206-fce0-0f349ea27f96" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.03215599961303717 0.03215767861360508\n" ] } ], "source": [ "delta = 1e-6\n", "high_area_abs = lambda x: (x>1.5)*(x<1.9)*(x>1.75*Lx)+(x>1.5)*(x<1.9)*(x<1.5*Lx)\n", "high_area = lambda x: (x<0.5*Lx)+(x>0.75*Lx)\n", "ts = np.linspace(0,tfinal,numsteps+1)\n", "intime = lambda t: (t>=1.8)*(t<=tfinal)*1+0\n", "sens_adj[3] = np.sum(dx*(sol[high_area_abs(xs2),:])*intime(ts)*dt)\n", "sens_adj[3] += np.sum(-sol[high_area(xs2),:]**4*sol_adj[high_area(xs2),:]*dt*dx)\n", "sol_pert,Q_v,xs2 = run_calc(Lx, Nx, rho, Source_func(xs,1), omega_nom, v_nom, \n", " kappa_func(xs,kappal_nom, kappah_nom+delta), dt, numsteps, 1.8,tfinal,1.5,1.9,IC)\n", "sens_for[3] = (Q_v-Q)/delta\n", "print(sens_for[3],sens_adj[3])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "colab_type": "code", "executionInfo": { "elapsed": 3307, "status": "ok", "timestamp": 1593495897636, "user": { "displayName": "黄尚晖", "photoUrl": "", "userId": "16496223364462243774" }, "user_tz": 240 }, "id": "qVvHnV1q6I-I", "nbpages": { "level": 3, "link": "[6.1.3.5 Compare sensitivities to parameters resulted from finite differences and the adjoint methodology: $\\rho$, $\\omega$,$\\kappa_l$,$\\kappa_h$,$q$ respectively](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.5-Compare-sensitivities-to-parameters-resulted-from-finite-differences-and-the-adjoint-methodology:-$\\rho$,-$\\omega$,$\\kappa_l$,$\\kappa_h$,$q$-respectively)", "section": "6.1.3.5 Compare sensitivities to parameters resulted from finite differences and the adjoint methodology: $\\rho$, $\\omega$,$\\kappa_l$,$\\kappa_h$,$q$ respectively" }, "outputId": "f9a2d5e9-4bbd-4a84-b589-142c7f8dcb3e" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.09638222440988553 0.09638749590621687\n" ] } ], "source": [ "sens_adj[4] = np.sum(dx*sol_adj[(xs2>0.25*Lx)*(xs2<0.75*Lx),:]*dt)\n", "sol_pert,Q_v,xs2 = run_calc(Lx, Nx, rho,Source_func(xs,1+delta), omega_nom, v_nom, \n", " kappa_func(xs,kappal_nom, kappah_nom), dt, numsteps, 1.8,tfinal,1.5,1.9,IC)\n", "sens_for[4] = (Q_v-Q)/delta\n", "print(sens_for[4],sens_adj[4])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 156 }, "colab_type": "code", "executionInfo": { "elapsed": 556, "status": "ok", "timestamp": 1593495917592, "user": { "displayName": "黄尚晖", "photoUrl": "", "userId": "16496223364462243774" }, "user_tz": 240 }, "id": "OTgiGJ_Y6I-L", "nbpages": { "level": 3, "link": "[6.1.3.5 Compare sensitivities to parameters resulted from finite differences and the adjoint methodology: $\\rho$, $\\omega$,$\\kappa_l$,$\\kappa_h$,$q$ respectively](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.5-Compare-sensitivities-to-parameters-resulted-from-finite-differences-and-the-adjoint-methodology:-$\\rho$,-$\\omega$,$\\kappa_l$,$\\kappa_h$,$q$-respectively)", "section": "6.1.3.5 Compare sensitivities to parameters resulted from finite differences and the adjoint methodology: $\\rho$, $\\omega$,$\\kappa_l$,$\\kappa_h$,$q$ respectively" }, "outputId": "b6640462-a0dd-4c1c-8d97-c60614b6f152" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "---------- -------------------- --------------------- -----------------------\n", "Parameters Finite Difference Adjoint Estimate Abs. Difference (10^-5)\n", "ρ -0.09947973321489112 -0.09948429352239387 4.584157350819606\n", "ω 0.28897547187101136 0.288993704223791 6.309308074348965\n", "κₗ -0.03022436832755826 -0.030226185976361303 6.013852079040496\n", "κₕ 0.03215599961303717 0.03215767861360508 5.221422403645427\n", "q 0.09638222440988553 0.09638749590621687 5.469365708888537\n", "---------- -------------------- --------------------- -----------------------\n" ] } ], "source": [ "output = np.vstack([sens_for,sens_adj,np.abs(sens_for-sens_adj)/np.abs(sens_for)/1e-5]).transpose()\n", "output1 = [[\"Parameters\" ,'Finite Difference','Adjoint Estimate','Abs. Difference (10^-5)'],['\\N{GREEK SMALL LETTER RHO}',output[0,0],output[0,1],output[0,2]],['\\N{GREEK SMALL LETTER omega}',output[1,0],output[1,1],output[1,2]],['\\N{GREEK SMALL LETTER kappa}\\u2097',output[2,0],output[2,1],output[2,2]],['\\N{GREEK SMALL LETTER kappa}\\u2095',output[3,0],output[3,1],output[3,2]],['q',output[4,0],output[4,1],output[4,2]]]\n", "import tabulate\n", "print(tabulate.tabulate(output1, tablefmt=\"Latex\", floatfmt=\".5f\"))\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "collapsed": true, "id": "p_WUJZ0v6I-P", "nbpages": { "level": 3, "link": "[6.1.3.6 Comparison for Reliability](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.6-Comparison-for-Reliability)", "section": "6.1.3.6 Comparison for Reliability" } }, "source": [ "### 6.1.3.6 Comparison for Reliability" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 1000 }, "colab_type": "code", "executionInfo": { "elapsed": 3861, "status": "ok", "timestamp": 1593523668753, "user": { "displayName": "黄尚晖", "photoUrl": "", "userId": "16496223364462243774" }, "user_tz": 240 }, "id": "OpwOPqfq6I-P", "nbpages": { "level": 3, "link": "[6.1.3.6 Comparison for Reliability](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.6-Comparison-for-Reliability)", "section": "6.1.3.6 Comparison for Reliability" }, "outputId": "5323499e-108c-456a-e6d6-98518be80bc2" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0. 0. 0. 0. 0.]\n", " [0. 0. 0. 0. 0.]\n", " [0. 0. 0. 0. 0.]\n", " [0. 0. 0. 0. 0.]\n", " [0. 0. 0. 0. 0.]]\n", "[[ 1. 0.1 -0.05 0. 0. ]\n", " [ 0.1 1. -0.4 0.3 0.5 ]\n", " [-0.05 -0.4 1. 0.2 0. ]\n", " [ 0. 0.3 0.2 1. -0.1 ]\n", " [ 0. 0.5 0. -0.1 1. ]]\n", "0.9922977210635758 0.010566538975926807 0.10279367186712812\n", "0.05021961567550527 2.4905472687980364e-05 0.004990538316452481\n", "0.10022954616433 9.927306454631912e-05 0.009963586931738946\n", "2.015168892538771 0.04332608207069289 0.20814918224843665\n", "1.0005110913540498 0.010210495960569357 0.10104699877071736\n", "[[ 1.05771161e-02 5.54030481e-05 -3.96931120e-05 3.54289587e-04\n", " -6.92092191e-05]\n", " [ 5.54030481e-05 2.49304031e-05 -2.20896335e-05 2.87449455e-04\n", " 2.44165203e-04]\n", " [-3.96931120e-05 -2.20896335e-05 9.93724370e-05 3.48562887e-04\n", " 7.25861375e-06]\n", " [ 3.54289587e-04 2.87449455e-04 3.48562887e-04 4.33694515e-02\n", " -2.45241544e-03]\n", " [-6.92092191e-05 2.44165203e-04 7.25861375e-06 -2.45241544e-03\n", " 1.02207167e-02]] 0.00024568041409483445\n", "(array([[0., 0., 0., ..., 0., 0., 0.],\n", " [0., 0., 0., ..., 0., 0., 0.],\n", " [0., 0., 0., ..., 0., 0., 0.],\n", " ...,\n", " [0., 0., 0., ..., 0., 0., 0.],\n", " [0., 0., 0., ..., 0., 0., 0.],\n", " [0., 0., 0., ..., 0., 0., 0.]]), 0.06981831556455251, array([-0.005, 0.005, 0.015, 0.025, 0.035, 0.045, 0.055, 0.065,\n", " 0.075, 0.085, 0.095, 0.105, 0.115, 0.125, 0.135, 0.145,\n", " 0.155, 0.165, 0.175, 0.185, 0.195, 0.205, 0.215, 0.225,\n", " 0.235, 0.245, 0.255, 0.265, 0.275, 0.285, 0.295, 0.305,\n", " 0.315, 0.325, 0.335, 0.345, 0.355, 0.365, 0.375, 0.385,\n", " 0.395, 0.405, 0.415, 0.425, 0.435, 0.445, 0.455, 0.465,\n", " 0.475, 0.485, 0.495, 0.505, 0.515, 0.525, 0.535, 0.545,\n", " 0.555, 0.565, 0.575, 0.585, 0.595, 0.605, 0.615, 0.625,\n", " 0.635, 0.645, 0.655, 0.665, 0.675, 0.685, 0.695, 0.705,\n", " 0.715, 0.725, 0.735, 0.745, 0.755, 0.765, 0.775, 0.785,\n", " 0.795, 0.805, 0.815, 0.825, 0.835, 0.845, 0.855, 0.865,\n", " 0.875, 0.885, 0.895, 0.905, 0.915, 0.925, 0.935, 0.945,\n", " 0.955, 0.965, 0.975, 0.985, 0.995, 1.005, 1.015, 1.025,\n", " 1.035, 1.045, 1.055, 1.065, 1.075, 1.085, 1.095, 1.105,\n", " 1.115, 1.125, 1.135, 1.145, 1.155, 1.165, 1.175, 1.185,\n", " 1.195, 1.205, 1.215, 1.225, 1.235, 1.245, 1.255, 1.265,\n", " 1.275, 1.285, 1.295, 1.305, 1.315, 1.325, 1.335, 1.345,\n", " 1.355, 1.365, 1.375, 1.385, 1.395, 1.405, 1.415, 1.425,\n", " 1.435, 1.445, 1.455, 1.465, 1.475, 1.485, 1.495, 1.505,\n", " 1.515, 1.525, 1.535, 1.545, 1.555, 1.565, 1.575, 1.585,\n", " 1.595, 1.605, 1.615, 1.625, 1.635, 1.645, 1.655, 1.665,\n", " 1.675, 1.685, 1.695, 1.705, 1.715, 1.725, 1.735, 1.745,\n", " 1.755, 1.765, 1.775, 1.785, 1.795, 1.805, 1.815, 1.825,\n", " 1.835, 1.845, 1.855, 1.865, 1.875, 1.885, 1.895, 1.905,\n", " 1.915, 1.925, 1.935, 1.945, 1.955, 1.965, 1.975, 1.985,\n", " 1.995, 2.005]))\n" ] } ], "source": [ "means = [rho, omega_nom, kappal_nom, kappah_nom, q_nom]\n", "varmat = np.zeros((5,5))\n", "#fill in diagonal\n", "corrmat = np.ones((5,5))\n", "corrmat[0,:] = (1,.1,-0.05,0,0)\n", "corrmat[1,:] = (.1,1,-.4,.3,.5)\n", "corrmat[2,:] = (-0.05,-.4,1,.2,0)\n", "corrmat[3,:] = (0,.3,0.2,1,-.1)\n", "corrmat[4,:] = (0,.5,0,-.1,1)\n", "print(corrmat-corrmat.transpose())\n", "print(corrmat)\n", "\n", "\n", "samps = 10**3\n", "test = norm.cdf(np.random.multivariate_normal(np.zeros(5), corrmat, samps))\n", "\n", "def gen_samps(samps,test):\n", " #rho will have v_nom = 1 v_var = .1 which makes alpha = 10 beta = 10 or theta = 1/10\n", " vsamps = gamma.ppf(test[:,0], a = 100, scale = 1/100)\n", " print(np.mean(vsamps), np.var(vsamps), np.std(vsamps))\n", " #plt.hist(vsamps)\n", "\n", " #omega will have omega_nom = .1, var = .1 which makes alpha = 100 beta = 5, theta = 1/5\n", " omegasamps = gamma.ppf(test[:,1], a = 100, scale = 1/2000)\n", " print(np.mean(omegasamps), np.var(omegasamps), np.std(omegasamps))\n", " #plt.hist(omegasamps)\n", "\n", " #kappa_l will have kappa_l = 0.1 var = (0.01)^2 this makes alpha = 100 and theta = 1/1000\n", " kappalsamps = gamma.ppf(test[:,2], a = 100, scale = 1/1000)\n", " print(np.mean(kappalsamps), np.var(kappalsamps), np.std(kappalsamps))\n", " #plt.hist(kappalsamps)\n", "\n", " #kappa_h will have kappa_h = 2 var = .04 this makes alpha = 100 and theta = 1/50\n", " kappahsamps = gamma.ppf(test[:,3], a = 100, scale = 1/50)\n", " print(np.mean(kappahsamps), np.var(kappahsamps), np.std(kappahsamps))\n", " #plt.hist(kappahsamps)\n", "\n", " #q will have q = 1 var = 0.01 this makes alpha = 100 and theta = 1/100\n", " qsamps = gamma.ppf(test[:,4], a = 100, scale = 1/100)\n", " print(np.mean(qsamps), np.var(qsamps),np.std(qsamps))\n", " #plt.hist(qsamps)\n", " \n", " return vsamps,omegasamps,kappalsamps,kappahsamps,qsamps\n", "vsamps,omegasamps,kappalsamps,kappahsamps,qsamps = gen_samps(samps,test)\n", "\n", "var_list = [vsamps,omegasamps,kappalsamps,kappahsamps,qsamps]\n", "cormat_emp = np.zeros((5,5))\n", "tmp = np.vstack((var_list[0],var_list[1],var_list[2],var_list[3],var_list[4]))\n", "cormat_emp = np.cov(tmp)\n", "sens = np.array([-0.0994797332704,0.288975471816,-0.0302243683414,0.0321559995853,0.096382224285])\n", "print(cormat_emp, np.dot(sens,np.dot(cormat_emp,sens)))\n", "\n", "print(run_calc(Lx, Nx, np.mean(vsamps), Source_func(xs,np.mean(qsamps)), np.mean(omegasamps), v_nom, \n", " kappa_func(xs,np.mean(kappalsamps), np.mean(kappahsamps)), dt, numsteps, 1.8,tfinal,1.5,1.9,IC))" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "colab": {}, "colab_type": "code", "executionInfo": { "elapsed": 2896055, "status": "ok", "timestamp": 1593526573639, "user": { "displayName": "黄尚晖", "photoUrl": "", "userId": "16496223364462243774" }, "user_tz": 240 }, "id": "jgdcjUIg6I-S", "nbpages": { "level": 3, "link": "[6.1.3.6 Comparison for Reliability](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.6-Comparison-for-Reliability)", "section": "6.1.3.6 Comparison for Reliability" } }, "outputs": [], "source": [ "outputs=np.zeros(samps)\n", "for i in range(samps):\n", " sol,Q,xs2 = run_calc(Lx, Nx, vsamps[i], Source_func(xs,qsamps[i]), omegasamps[i], v_nom, \n", " kappa_func(xs,kappalsamps[i], kappahsamps[i]), dt, numsteps, 1.8,tfinal,1.5,1.9,IC)\n", " outputs[i] = Q\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", "id": "i1irw8Sk6I-U", "nbpages": { "level": 3, "link": "[6.1.3.6 Comparison for Reliability](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.6-Comparison-for-Reliability)", "section": "6.1.3.6 Comparison for Reliability" } }, "outputs": [], "source": [ "print(np.var(outputs),outputs,vsamps)\n", "\n", "np.savetxt(\"nonlinear_ref_adr.csv\", outputs, delimiter=\",\")" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 352 }, "colab_type": "code", "executionInfo": { "elapsed": 541, "status": "ok", "timestamp": 1593527144314, "user": { "displayName": "黄尚晖", "photoUrl": "", "userId": "16496223364462243774" }, "user_tz": 240 }, "id": "yNBINPK66I-W", "nbpages": { "level": 3, "link": "[6.1.3.6 Comparison for Reliability](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.6-Comparison-for-Reliability)", "section": "6.1.3.6 Comparison for Reliability" }, "outputId": "8b75efbf-4620-46fb-c47a-300b24294701" }, "outputs": [ { "data": { "text/plain": [ "(array([ 12., 41., 121., 210., 245., 199., 108., 49., 13., 2.]),\n", " array([0.03171055, 0.04013366, 0.04855677, 0.05697988, 0.06540299,\n", " 0.0738261 , 0.08224921, 0.09067232, 0.09909542, 0.10751853,\n", " 0.11594164]),\n", " )" ] }, "execution_count": 20, "metadata": { "tags": [] }, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD4CAYAAAAKA1qZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAOjElEQVR4nO3df8xeZX3H8fdn1B9T3Cz2WYOl7FFT/8BMK3lkJLqNjUwRkhUTw9BMG2dSk0EiiW4r+odsCUtnVBazjaVGYllUYKixCUSGnYlzmT8KwUpBpEKRdpVWMP4iUVu/++M51dtSnufpfe4fT3u9X8md+9zXOee+vucK5dPr3OecpqqQJLXpN6ZdgCRpegwBSWqYISBJDTMEJKlhhoAkNWzFtAsAWLVqVc3Ozk67DEk6qdx1113fq6qZPt+xLEJgdnaWnTt3TrsMSTqpJHmk73d4OkiSGrZoCCRZm+QLSe5LsjvJO7v2a5LsT3JP97p4YJ+rk+xJ8kCS143zACRJw1vK6aDDwLuq6u4kzwPuSnJnt+66qvrA4MZJzgEuB14GvBD4fJKXVtWRURYuSepv0ZlAVR2oqru75R8B9wNrFthlA3BTVf20qh4G9gDnjaJYSdJondBvAklmgVcCX+markyyK8kNSVZ2bWuARwd228dxQiPJpiQ7k+w8dOjQCRcuSepvySGQ5HTgU8BVVfVD4HrgJcB64ADwwRPpuKq2VtVcVc3NzPS6wkmSNKQlhUCSZzAfAB+vqk8DVNVjVXWkqn4BfIRfnfLZD6wd2P2srk2StMws5eqgAB8F7q+qDw20nzmw2RuAe7vl7cDlSZ6V5EXAOuCroytZkjQqS7k66NXAW4BvJLmna3sP8KYk64EC9gLvAKiq3UluAe5j/sqiK7wySJKWp0VDoKq+BOQ4q25fYJ9rgWt71CX9mtnNt02t771bLpla39K4ecewJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGrZi2gXo5DK7+bZplyBphJwJSFLDDAFJapghIEkNMwQkqWGGgCQ1bNEQSLI2yReS3Jdkd5J3du1nJLkzyYPd+8quPUk+nGRPkl1Jzh33QUiShrOUmcBh4F1VdQ5wPnBFknOAzcCOqloH7Og+A7weWNe9NgHXj7xqSdJILBoCVXWgqu7uln8E3A+sATYA27rNtgGXdssbgBtr3peB5yc5c+SVS5J6O6HfBJLMAq8EvgKsrqoD3arvAqu75TXAowO77evajv2uTUl2Jtl56NChEyxbkjQKS75jOMnpwKeAq6rqh0l+ua6qKkmdSMdVtRXYCjA3N3dC+0qTNK27pPduuWQq/aotS5oJJHkG8wHw8ar6dNf82NHTPN37wa59P7B2YPezujZJ0jKzlKuDAnwUuL+qPjSwajuwsVveCHx2oP2t3VVC5wM/GDhtJElaRpZyOujVwFuAbyS5p2t7D7AFuCXJ24FHgMu6dbcDFwN7gCeBt420YknSyCwaAlX1JSBPs/rC42xfwBU965IkTYB3DEtSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMWDYEkNyQ5mOTegbZrkuxPck/3unhg3dVJ9iR5IMnrxlW4JKm/pcwEPgZcdJz266pqffe6HSDJOcDlwMu6ff41yWmjKlaSNFqLhkBVfRF4YonftwG4qap+WlUPA3uA83rUJ0kaoz6/CVyZZFd3umhl17YGeHRgm31d21Mk2ZRkZ5Kdhw4d6lGGJGlYw4bA9cBLgPXAAeCDJ/oFVbW1quaqam5mZmbIMiRJfQwVAlX1WFUdqapfAB/hV6d89gNrBzY9q2uTJC1DQ4VAkjMHPr4BOHrl0Hbg8iTPSvIiYB3w1X4lSpLGZcViGyT5JHABsCrJPuB9wAVJ1gMF7AXeAVBVu5PcAtwHHAauqKoj4yldktTXoiFQVW86TvNHF9j+WuDaPkVJkibDO4YlqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhi36bwxr+ZndfNu0S5B0inAmIEkNMwQkqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhPjZCWqam9XiQvVsumUq/mg5nApLUMENAkhq2aAgkuSHJwST3DrSdkeTOJA927yu79iT5cJI9SXYlOXecxUuS+lnKTOBjwEXHtG0GdlTVOmBH9xng9cC67rUJuH40ZUqSxmHREKiqLwJPHNO8AdjWLW8DLh1ov7HmfRl4fpIzR1WsJGm0hv1NYHVVHeiWvwus7pbXAI8ObLeva3uKJJuS7Eyy89ChQ0OWIUnqo/cPw1VVQA2x39aqmququZmZmb5lSJKGMGwIPHb0NE/3frBr3w+sHdjurK5NkrQMDRsC24GN3fJG4LMD7W/trhI6H/jBwGkjSdIys+gdw0k+CVwArEqyD3gfsAW4JcnbgUeAy7rNbwcuBvYATwJvG0PNkqQRWTQEqupNT7PqwuNsW8AVfYuSJE2GdwxLUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDVvTZOcle4EfAEeBwVc0lOQO4GZgF9gKXVdX3+5UpSRqHUcwE/riq1lfVXPd5M7CjqtYBO7rPkqRlaByngzYA27rlbcClY+hDkjQCfUOggP9McleSTV3b6qo60C1/F1jdsw9J0pj0+k0AeE1V7U/yO8CdSb45uLKqKkkdb8cuNDYBnH322T3LkDQqs5tvm1rfe7dcMrW+W9VrJlBV+7v3g8BngPOAx5KcCdC9H3yafbdW1VxVzc3MzPQpQ5I0pKFDIMlzkzzv6DLwWuBeYDuwsdtsI/DZvkVKksajz+mg1cBnkhz9nk9U1eeSfA24JcnbgUeAy/qXKUkah6FDoKoeAl5xnPbHgQv7FCVJmgzvGJakhhkCktSwvpeINm2al9JJ0ig4E5CkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSpYSumXYAkHTW7+bap9Lt3yyVT6Xc5OOlDYFr/0UjSqcDTQZLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhYwuBJBcleSDJniSbx9WPJGl4Y7lPIMlpwL8AfwrsA76WZHtV3TeO/iSpj2nebzTtG9XGNRM4D9hTVQ9V1c+Am4ANY+pLkjSkcd0xvAZ4dODzPuD3BzdIsgnY1H38cZIHxlTLuKwCvjftIpY5x2hhjs/Cmhif/OPQu64Cfrdv/1N7bERVbQW2Tqv/vpLsrKq5adexnDlGC3N8Fub4LKwbn9m+3zOu00H7gbUDn8/q2iRJy8i4QuBrwLokL0ryTOByYPuY+pIkDWksp4Oq6nCSK4E7gNOAG6pq9zj6mqKT9lTWBDlGC3N8Fub4LGwk45OqGsX3SJJOQt4xLEkNMwQkqWGGwHEs9siLJM9KcnO3/itJZo9Zf3aSHyd596RqnqQ+45Pk5Un+N8nuJN9I8uxJ1j4pw45Rkmck2daNzf1Jrp507ZOwhPH5wyR3Jzmc5I3HrNuY5MHutXFyVU/OsOOTZP3An69dSf580c6qytfAi/kfsr8NvBh4JvB14Jxjtvkr4N+65cuBm49ZfyvwH8C7p308y2l8mL8QYRfwiu7zC4DTpn1My2yM3gzc1C0/B9gLzE77mKYwPrPAy4EbgTcOtJ8BPNS9r+yWV077mJbR+LwUWNctvxA4ADx/of6cCTzVUh55sQHY1i3fClyYJABJLgUeBk61q6GO6jM+rwV2VdXXAarq8ao6MqG6J6nPGBXw3CQrgN8Efgb8cDJlT8yi41NVe6tqF/CLY/Z9HXBnVT1RVd8H7gQumkTREzT0+FTVt6rqwW75/4CDwMxCnRkCT3W8R16sebptquow8APgBUlOB/4W+LsJ1DktQ48P839LqSR3dFPZv5lAvdPQZ4xuBX7C/N/gvgN8oKqeGHfBE7aU8RnHvieLkRxjkvOYn0l8e6HtpvbYiFPUNcB1VfXjbmKgX7cCeA3wKuBJYEeSu6pqx3TLWlbOA44wP5VfCfx3ks9X1UPTLUsnkyRnAv8ObKyqY2dTv8aZwFMt5ZEXv9ymm7b/NvA48w/Je3+SvcBVwHu6m+ZOJX3GZx/wxar6XlU9CdwOnDv2iievzxi9GfhcVf28qg4C/wOcas/P6fNYmRYeSdPrGJP8FnAb8N6q+vJi2xsCT7WUR15sB45elfBG4L9q3h9U1WzNP9Tpn4B/qKp/nlThEzL0+DB/B/nvJXlO9z++PwJOxX9jos8YfQf4E4AkzwXOB745kaonp89jZe4AXptkZZKVzP/OdMeY6pyWocen2/4zwI1VdeuSepv2L+HL8QVcDHyL+XNp7+3a/h74s2752cxf/bMH+Crw4uN8xzWcglcH9R0f4C+Y/9H8XuD90z6W5TZGwOld+27mA/Kvp30sUxqfVzE/c/wJ8zOk3QP7/mU3bnuAt037WJbT+HR/vn4O3DPwWr9QXz42QpIa5ukgSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIa9v+95DfyN7jJNwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light", "tags": [] }, "output_type": "display_data" } ], "source": [ "plt.hist(outputs)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": {}, "colab_type": "code", "id": "fVIu49VO6I-Z", "nbpages": { "level": 3, "link": "[6.1.3.6 Comparison for Reliability](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.6-Comparison-for-Reliability)", "section": "6.1.3.6 Comparison for Reliability" } }, "outputs": [], "source": [ "#lhd will have the values in 0 to 1\n", "samps =10**3\n", "lhd = lhs(5, samples=samps)\n", "#now i need to turn these into samples from N(0,Corrmat)\n", "#do cholesky fact\n", "chol = np.linalg.cholesky(corrmat)\n", "\n", "lhs_unif = np.zeros((samps,5))\n", "for i in range(samps):\n", " lhs_unif[i,:] = np.dot(chol,norm.ppf(lhd[i,:]))\n", "test_lhs = norm.cdf(lhs_unif)\n", "\n", "vsamps,omegasamps,kappalsamps,kappahsamps,qsamps = gen_samps(samps,test_lhs)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 303 }, "colab_type": "code", "executionInfo": { "elapsed": 2830651, "status": "ok", "timestamp": 1593488077187, "user": { "displayName": "黄尚晖", "photoUrl": "", "userId": "16496223364462243774" }, "user_tz": 240 }, "id": "vj98tkLP6I-c", "nbpages": { "level": 3, "link": "[6.1.3.6 Comparison for Reliability](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.01-Contributed-Example.html#6.1.3.6-Comparison-for-Reliability)", "section": "6.1.3.6 Comparison for Reliability" }, "outputId": "8eab34ee-852f-4d4d-ff38-563c5b0238b1" }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAOVklEQVR4nO3df6zdd13H8edrKz9kGOlsaeY2uMWUhE6hwmUjEU3NIhtroCOSuaHQIEk1jkSMKAX/YJqQVKKSkOhIDWNdorA5xS1sDkc1gj/4cYtbWZmTMjq3WtbLxq8xBDve/nG/dad3t7s/zjn3fPfx+UhOzvd8vt/v+b521vu6n/v9nnNvqgpJUltOm3QASdLoWe6S1CDLXZIaZLlLUoMsd0lq0JpJBwBYt25dTU1NTTqGJD2l7N+//2tVtX6hdb0o96mpKWZmZiYdQ5KeUpLcd6p1npaRpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QG9eITqtJipnbdMrFjH969bWLHllbKmbskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDVq03JOcm+QfknwxycEkv9GNn5nk9iRf6u7XduNJ8v4kh5IcSPLScf9HSJJOtpSZ+3Hgt6pqM/AK4Mokm4FdwL6q2gTs6x4DvBrY1N12AlePPLUk6UktWu5VdbSqPt8tfxu4Gzgb2A7s7TbbC1zaLW8Hrqs5nwaek+SskSeXJJ3Sss65J5kCfgr4DLChqo52q74KbOiWzwbuH9jtgW5MkrRKllzuSZ4N/BXwtqr61uC6qiqglnPgJDuTzCSZmZ2dXc6ukqRFrFnKRkmexlyx/3lV/XU3/GCSs6rqaHfa5Vg3fgQ4d2D3c7qxk1TVHmAPwPT09LK+MUiraWrXLRM57uHd2yZyXLVhKe+WCfBB4O6q+uOBVTcDO7rlHcBNA+Nv6t418wrgmwOnbyRJq2ApM/efBt4IfCHJHd3Yu4DdwA1J3gLcB1zWrbsVuAQ4BDwKvHmkiSVJi1q03Kvqn4CcYvWFC2xfwJVD5pIkDcFPqEpSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGrSkv6EqnTCpvycqaXmcuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNWrTck1yT5FiSuwbGrkpyJMkd3e2SgXXvTHIoyT1JLhpXcEnSqS1l5n4tcPEC4++rqi3d7VaAJJuBy4Hzun3+NMnpoworSVqaRcu9qj4JPLzE59sOfKSqvldVXwEOAecPkU+StALDnHN/a5ID3Wmbtd3Y2cD9A9s80I09QZKdSWaSzMzOzg4RQ5I030rL/Wrgx4EtwFHgj5b7BFW1p6qmq2p6/fr1K4whSVrIisq9qh6sqseq6gfAn/H4qZcjwLkDm57TjUmSVtGKyj3JWQMPXweceCfNzcDlSZ6RZCOwCfjscBElScu1ZrENknwY2AqsS/IA8G5ga5ItQAGHgV8FqKqDSW4AvggcB66sqsfGE12SdCqLlntVXbHA8AefZPv3AO8ZJpQkaTh+QlWSGmS5S1KDLHdJapDlLkkNWvSCqqTJmNp1y0SOe3j3tokcV6PlzF2SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAYtWu5JrklyLMldA2NnJrk9yZe6+7XdeJK8P8mhJAeSvHSc4SVJC1vKzP1a4OJ5Y7uAfVW1CdjXPQZ4NbCpu+0Erh5NTEnScixa7lX1SeDhecPbgb3d8l7g0oHx62rOp4HnJDlrVGElSUuz0nPuG6rqaLf8VWBDt3w2cP/Adg90Y0+QZGeSmSQzs7OzK4whSVrI0BdUq6qAWsF+e6pquqqm169fP2wMSdKAlZb7gydOt3T3x7rxI8C5A9ud041JklbRSsv9ZmBHt7wDuGlg/E3du2ZeAXxz4PSNJGmVrFlsgyQfBrYC65I8ALwb2A3ckOQtwH3AZd3mtwKXAIeAR4E3jyGzJGkRi5Z7VV1xilUXLrBtAVcOG0qSNBw/oSpJDbLcJalBlrskNchyl6QGLXpBVf0zteuWSUeQ1HPO3CWpQZa7JDXIcpekBlnuktQgy12SGuS7ZSSdZJLvxjq8e9vEjt0aZ+6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1KA1w+yc5DDwbeAx4HhVTSc5E7gemAIOA5dV1deHiylJWo5RzNx/rqq2VNV093gXsK+qNgH7useSpFU0jtMy24G93fJe4NIxHEOS9CSGLfcC/i7J/iQ7u7ENVXW0W/4qsGGhHZPsTDKTZGZ2dnbIGJKkQUOdcwdeWVVHkjwXuD3Jvw+urKpKUgvtWFV7gD0A09PTC24jSVqZoWbuVXWkuz8GfBQ4H3gwyVkA3f2xYUNKkpZnxeWe5IwkP3xiGXgVcBdwM7Cj22wHcNOwISVJyzPMaZkNwEeTnHiev6iq25J8DrghyVuA+4DLho8pSVqOFZd7Vd0LvGSB8YeAC4cJJUkajp9QlaQGWe6S1KBh3wr5/9rUrlsmHUGSFuTMXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgP8QkqTcm9cHAw7u3TeS44+TMXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDXrK/4HsSf1BXUnqM2fuktSgp/zMXZKGNckzAId3bxvL8zpzl6QGja3ck1yc5J4kh5LsGtdxJElPNJZyT3I68CfAq4HNwBVJNo/jWJKkJxrXzP184FBV3VtV3wc+Amwf07EkSfOM64Lq2cD9A48fAC4Y3CDJTmBn9/CRJPes4DjrgK+tKOH49TVbX3OB2VbKbCvTi2z5gwWHl5rt+adaMbF3y1TVHmDPMM+RZKaqpkcUaaT6mq2vucBsK2W2lWk927hOyxwBzh14fE43JklaBeMq988Bm5JsTPJ04HLg5jEdS5I0z1hOy1TV8SRvBT4OnA5cU1UHx3CooU7rjFlfs/U1F5htpcy2Mk1nS1WNIogkqUf8hKokNchyl6QG9bLcF/vVBUmekeT6bv1nkkx14+cnuaO73ZnkdX3JNrD+eUkeSfL2vmRLMpXkuwOv3Qf6kq1b9+Ik/5rkYJIvJHlmH7Il+aWB1+yOJD9IsqUn2Z6WZG/3et2d5J09yfX0JB/qct2ZZOsocy0x288m+XyS40leP2/djiRf6m47epbttiTfSPKxJR2sqnp1Y+4C7JeBFwBPB+4ENs/b5teBD3TLlwPXd8vPAtZ0y2cBx048nnS2gfU3An8JvL1Hr9sUcFdP/5+uAQ4AL+ke/yhweh+yzdvmJ4Ev9+h1ewPwkW75WcBhYKoHua4EPtQtPxfYD5y2yq/ZFPBi4Drg9QPjZwL3dvdru+W1fcjWrbsQeA3wsaUcr48z96X86oLtwN5u+UbgwiSpqker6ng3/kxg1FeLV5wNIMmlwFeAcbxzaKhsYzZMtlcBB6rqToCqeqiqHutJtkFXdPuO0jDZCjgjyRrgh4DvA9/qQa7NwN8DVNUx4BvAKD9ItGi2qjpcVQeAH8zb9yLg9qp6uKq+DtwOXNyTbFTVPuDbSz1YH8t9oV9dcPaptunK/JvMzehIckGSg8AXgF8bKPuJZkvybOAdwO+NMM9IsnXrNib5tyT/mORnepTthUAl+Xj34+rv9CjboF8EPtyjbDcC3wGOAv8J/GFVPdyDXHcCr02yJslG4GWc/IHH1cg2jn378Pwnae6PdVTVZ4DzkrwI2Jvkb6vqvyedC7gKeF9VPbI6k+VlOQo8r6oeSvIy4G+SnFdVo5rpDWMN8Erg5cCjwL4k+7tZTC8kuQB4tKrumnSWAecDjwE/xtwphk8l+URV3TvZWFwDvAiYAe4D/oW5nBqxPs7cl/KrC/5vm+7Hzh8BHhrcoKruBh4BfqIn2S4A3pvkMPA24F2Z+6DXxLNV1feq6iGAqtrP3HnBF/YhG3Ozm09W1deq6lHgVuClPcl2wuWMftY+bLY3ALdV1f90pz/+mdGd/hjm39rxqvrNqtpSVduB5wD/MaJcS802jn378PwnG9XFghFedFjD3IWMjTx+0eG8edtcyckXa27oljfy+AXV5wP/BazrQ7Z521zF6C+oDvO6rae7SMncxZ4jwJk9ybYW+DzdxXLgE8C2PmTrHp/WvV4v6NnXwjt4/MLlGcAXgRf3INezgDO65Z9n7hv3qr5mA9teyxMvqH6l+ze3tlte1a+DU2UbGN/KEi+ojvQf4whfhEuY+27+ZeB3u7HfB17bLT+TuXecHAI+e+ILC3gjcxcr7+gK4dK+ZJv3HFcx4nIf8nX7hXmv22v6kq1b98tdvruA9/Ys21bg0z38Wnh2N36QuWL/7Z7kmgLuAe5m7hv18yfwmr2cuZ8Iv8PcTzkHB/b9lS7zIeDNPcv2KWAW+G63zUVPdix//YAkNaiP59wlSUOy3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KD/hdoeSNcojH+mgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light", "tags": [] }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "0.06851216544161239 0.09131389762764693 0.046238313965240176 0.013419289882262094 -0.2772953252193724 0.07380452982344808\n" ] } ], "source": [ "QLHS = np.zeros(samps)\n", "outputs_LHS=np.zeros(samps)\n", "for i in range(samps):\n", " sol,QLHS[i],xs2 = run_calc(Lx, Nx, vsamps[i], Source_func(xs,qsamps[i]), omegasamps[i], v_nom, \n", " kappa_func(xs,kappalsamps[i], kappahsamps[i]), dt, numsteps, 1.8,tfinal,1.5,1.9,IC)\n", " outputs_LHS[i] = QLHS[i]\n", "\n", "\n", "plt.hist(QLHS)\n", "plt.show()\n", "print(np.mean(QLHS),stats.scoreatpercentile(QLHS,95),stats.scoreatpercentile(QLHS,5),np.std(QLHS), \n", " stats.kurtosis(QLHS), stats.skew(QLHS))\n", "#print(np.mean(Qs),stats.scoreatpercentile(Qs,95) ,stats.scoreatpercentile(Qs,5),np.std(Qs), stats.kurtosis(Qs), stats.skew(Qs))\n", "\n", "np.savetxt(\"nonlinear_ref_adr_LHS.csv\", outputs_LHS, delimiter=\",\")\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "< [6.0 Adjoint-Based Local Sensitivity Analysis](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.00-Adjoint-Based-Local-Sensitivity-Analysis.html) | [Contents](toc.html) | [6.2 A Simple Example of Adjoint Sensitivity Analysis](https://ndcbe.github.io/cbe67701-uncertainty-quantification/06.02-Contributed-Example.html)

\"Open

\"Download\"" ] } ], "metadata": { "colab": { "collapsed_sections": [ "mazDyC_6Yo14", "I3y7UScJOOQV", "_pLnbneDZTUX", "igOTigTMeFBq" ], "name": "Nonlinear ADR.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 }