simulators/teasers/PyReSim_teaser1.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "PyReSim-teaser1.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true,
      "authorship_tag": "ABX9TyPMGIOHo7tU6VdUztPvymSL",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/yohanesnuwara/pyresim/blob/master/simulators/teasers/PyReSim_teaser1.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "a4I79FO4rRFb",
        "colab_type": "text"
      },
      "source": [
        "# Python Reservoir Simulation (PyReSim) Teaser #1"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pt353HJYrX7b",
        "colab_type": "text"
      },
      "source": [
        "Read the simulation case description [here](https://github.com/yohanesnuwara/pyresim/blob/master/docs/teaser.md#teaser-1)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "J-Sf8NNT3kbo",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import pandas as pd"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HrKrhulN4urc",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 138
        },
        "outputId": "e5bbefdc-0d24-4d82-c4b3-cec40d420928"
      },
      "source": [
        "!git clone https://github.com/yohanesnuwara/pyresim"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Cloning into 'pyresim'...\n",
            "remote: Enumerating objects: 211, done.\u001b[K\n",
            "remote: Counting objects: 100% (211/211), done.\u001b[K\n",
            "remote: Compressing objects: 100% (207/207), done.\u001b[K\n",
            "remote: Total 980 (delta 133), reused 0 (delta 0), pack-reused 769\n",
            "Receiving objects: 100% (980/980), 348.26 KiB | 587.00 KiB/s, done.\n",
            "Resolving deltas: 100% (579/579), done.\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lmU_Khbj4zAc",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import sys\n",
        "sys.path.append('/content/pyresim/functions')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ovm13Nqb46Lh",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from input_output import read_input\n",
        "from boundary import boundary2d_location, boundary_flow2d_constant_pressuregrad, boundary_floweq2d\n",
        "from gridding import fill_active_blocks, source2d\n",
        "from solver import lhs_coeffs2d_welltype, rhs_constant2d_welltype, fill2d_lhs_mat, fill2d_rhs_mat\n",
        "from transmissibility import transmissibility2d, transmissibility2d_boundary\n",
        "from wellblock import prodrate1d, fraction_wellblock_geometric_factor\n",
        "from wellsimulation import solution_well2d"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "h1CgJw9E-VPJ",
        "colab_type": "text"
      },
      "source": [
        "## Simulator"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AKdI5GFu5zKv",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 374
        },
        "outputId": "79d4c7f1-71dd-4502-c4f2-eceb01ea8b5a"
      },
      "source": [
        "import warnings\n",
        "\n",
        "warnings.filterwarnings(\"ignore\", category=FutureWarning)\n",
        "\n",
        "\"\"\"\"\"\"\"\"\"\"\"\n",
        "PUT THE INPUT TXT DATA\n",
        "\"\"\"\"\"\"\"\"\"\"\"\n",
        "\n",
        "filepath = '/content/pyresim/input/teasers/teaser1_input.txt'\n",
        "\n",
        "# call the read_input function to read input data\n",
        "reservoir_input, wells, west_boundary, east_boundary, south_boundary, north_boundary = read_input(filepath)\n",
        "\n",
        "\"\"\"\"\"\"\"\"\"\"\"\n",
        "INPUT PROCESSING\n",
        "\"\"\"\"\"\"\"\"\"\"\"\n",
        "\n",
        "# number of blocks in x and y\n",
        "xi = reservoir_input['xi']\n",
        "yi = reservoir_input['yi']\n",
        "\n",
        "# define parameters, in this case blocks are homogeneous and same in size\n",
        "dx = reservoir_input['dx'] \n",
        "dy = reservoir_input['dy']\n",
        "dz = reservoir_input['dz']\n",
        "kx = reservoir_input['kx']\n",
        "ky = reservoir_input['ky']\n",
        "B = reservoir_input['B']\n",
        "mu = reservoir_input['mu']\n",
        "\n",
        "# well information (dataframe)\n",
        "well_name = wells['well_name']\n",
        "well_rw = wells['well_rw'] / 12  # wellbore radius, inch to ft\n",
        "well_skin = wells['well_skin']\n",
        "well_loc = wells['well_loc']\n",
        "well_condition = wells['well_condition']\n",
        "well_value = wells['well_value']\n",
        "well_config = wells['well_config']\n",
        "\n",
        "\"\"\"\"\"\"\"\"\"\"\"\n",
        "GRIDDING\n",
        "\"\"\"\"\"\"\"\"\"\"\"\n",
        "\n",
        "\" Create grid points \"\n",
        "# meshgrid the original points\n",
        "x_ = np.arange(1, xi + 1)\n",
        "y_ = np.arange(1, yi + 1)\n",
        "\n",
        "x, y = np.meshgrid(x_, y_, indexing='ij')\n",
        "\n",
        "\" Handle irregularities \"\n",
        "# NOTHING. The grid is regular\n",
        "\n",
        "\" Classify the location of boundary with codes (1, 12, 13, etc) \"\n",
        "bound_loc = boundary2d_location(x, y, xi, yi)\n",
        "\n",
        "\"\"\"\"\"\"\"\"\"\"\"\n",
        "WELL INFORMATION PROCESSING\n",
        "\"\"\"\"\"\"\"\"\"\"\"\n",
        "\n",
        "# source block (production or injection well)\n",
        "qsc = source2d(well_value, np.array(well_loc), xi, yi)  # call function HERE\n",
        "\n",
        "# Calculate wellblock geometric factor\n",
        "kh, r_eq, Gw = [], [], []\n",
        "\n",
        "for i in range(len(well_config)):\n",
        "    kh_, r_eq_, Gw_ = fraction_wellblock_geometric_factor(dx, dy, kx, ky, well_skin[i],\n",
        "                                                          well_rw[i], dz, well_config[i])\n",
        "    kh.append(kh_);\n",
        "    r_eq.append(r_eq_);\n",
        "    Gw.append(Gw_)\n",
        "\n",
        "# create WELLBLOCK PROPERTY \n",
        "# inject (fill) the operating wells to 'qsc' grid that we created before\n",
        "# the values won't be filled in the 'qsc' grid that has value = 0\n",
        "qsc = qsc.astype('float')\n",
        "qsc[qsc == 0] = np.nan\n",
        "\n",
        "well_value = qsc\n",
        "well_condition = fill_active_blocks(well_condition, qsc)\n",
        "well_rw = fill_active_blocks(well_rw, qsc)\n",
        "well_Gw = fill_active_blocks(np.array(Gw), qsc)\n",
        "well_mu = fill_active_blocks(np.array(mu), qsc)\n",
        "well_B = fill_active_blocks(np.array(B), qsc)\n",
        "\n",
        "\"\"\"\"\"\"\"\"\"\"\"\n",
        "SIMULATION\n",
        "\"\"\"\"\"\"\"\"\"\"\"\n",
        "\n",
        "def lookup(i, j, bound_dict):\n",
        "    \"\"\"\n",
        "    Determine the type and value of boundary at a grid block\n",
        "    given the boundary dictionary\n",
        "    \"\"\"\n",
        "    if bound_dict['loc'] == 'all':\n",
        "        bound_type = bound_dict['type']\n",
        "        bound_value = bound_dict['value']\n",
        "    else:\n",
        "        x, y = i + 1, j + 1\n",
        "\n",
        "        id1 = (np.where(np.all(bound_dict['loc'] == np.array([x, y]), axis=-1)))[0][0]\n",
        "        id2 = (np.where(np.all((bound_dict['loc'][id1]) == np.array([x, y]), axis=-1)))[0][0]\n",
        "\n",
        "        bound_type = (bound_dict['type'])[id1]\n",
        "        bound_value = (bound_dict['value'])[id1]\n",
        "    return bound_type, bound_value\n",
        "\n",
        "# Create block index\n",
        "block_index = np.arange(1, (xi * yi) + 1)\n",
        "block_index = (np.reshape(block_index, (-1, xi))).T\n",
        "\n",
        "# Create empty matrix for LHS and RHS\n",
        "lhs_mat = np.array([[0.] * (xi * yi)] * (xi * yi))\n",
        "rhs_mat = np.array([[0.] * 1] * (xi * yi))\n",
        "\n",
        "for i in range(xi):\n",
        "    for j in range(yi):\n",
        "\n",
        "        if bound_loc[i, j] == 13:\n",
        "            \"\"\"\"\"\"\"\"\"\"\"\n",
        "            MODIFY THIS W.R.T. THE GRID BLOCK LOCATION \n",
        "            \"\"\"\"\"\"\"\"\"\"\"\n",
        "\n",
        "            # get the TYPE and VALUE of the boundary from the BOUNDARY DICTIONARY\n",
        "            b1_type, b1_value = lookup(i, j, west_boundary)\n",
        "            b3_type, b3_value = lookup(i, j, south_boundary)\n",
        "\n",
        "            # create boundary dictionary\n",
        "            boundary = {'loc': np.array([1, 3]),\n",
        "                        'type': np.array([b1_type, b3_type]),\n",
        "                        'value': np.array([b1_value, b3_value])}\n",
        "\n",
        "            # create non-boundary dictionary\n",
        "            non_boundary = {'loc': np.array([2, 4])}\n",
        "\n",
        "            # create well dictionary\n",
        "            well = {'condition': well_condition[i, j],\n",
        "                    'value': well_value[i, j],\n",
        "                    'rw': well_rw[i, j],\n",
        "                    'Gw': well_Gw[i, j]}\n",
        "\n",
        "        if bound_loc[i, j] == 14:\n",
        "            \"\"\"\"\"\"\"\"\"\"\"\n",
        "            MODIFY THIS W.R.T. THE GRID BLOCK LOCATION \n",
        "            \"\"\"\"\"\"\"\"\"\"\"\n",
        "\n",
        "            # get the TYPE and VALUE of the boundary from the BOUNDARY DICTIONARY\n",
        "            b1_type, b1_value = lookup(i, j, west_boundary)\n",
        "            b4_type, b4_value = lookup(i, j, north_boundary)\n",
        "\n",
        "            # create boundary dictionary\n",
        "            boundary = {'loc': np.array([1, 4]),\n",
        "                        'type': np.array([b1_type, b4_type]),\n",
        "                        'value': np.array([b1_value, b4_value])}\n",
        "\n",
        "            # create non-boundary dictionary\n",
        "            non_boundary = {'loc': np.array([2, 3])}\n",
        "\n",
        "            # create well dictionary\n",
        "            well = {'condition': well_condition[i, j],\n",
        "                    'value': well_value[i, j],\n",
        "                    'rw': well_rw[i, j],\n",
        "                    'Gw': well_Gw[i, j]}\n",
        "\n",
        "        if bound_loc[i, j] == 23:\n",
        "            \"\"\"\"\"\"\"\"\"\"\"\n",
        "            MODIFY THIS W.R.T. THE GRID BLOCK LOCATION \n",
        "            \"\"\"\"\"\"\"\"\"\"\"\n",
        "\n",
        "            # get the TYPE and VALUE of the boundary from the BOUNDARY DICTIONARY\n",
        "            b2_type, b2_value = lookup(i, j, east_boundary)\n",
        "            b3_type, b3_value = lookup(i, j, south_boundary)\n",
        "\n",
        "            # create boundary dictionary\n",
        "            boundary = {'loc': np.array([2, 3]),\n",
        "                        'type': np.array([b2_type, b3_type]),\n",
        "                        'value': np.array([b2_value, b3_value])}\n",
        "\n",
        "            # create non-boundary dictionary\n",
        "            non_boundary = {'loc': np.array([1, 4])}\n",
        "\n",
        "            # create well dictionary\n",
        "            well = {'condition': well_condition[i, j],\n",
        "                    'value': well_value[i, j],\n",
        "                    'rw': well_rw[i, j],\n",
        "                    'Gw': well_Gw[i, j]}\n",
        "\n",
        "        if bound_loc[i, j] == 24:\n",
        "            \"\"\"\"\"\"\"\"\"\"\"\n",
        "            MODIFY THIS W.R.T. THE GRID BLOCK LOCATION \n",
        "            \"\"\"\"\"\"\"\"\"\"\"\n",
        "\n",
        "            # get the TYPE and VALUE of the boundary from the BOUNDARY DICTIONARY\n",
        "            b2_type, b2_value = lookup(i, j, east_boundary)\n",
        "            b4_type, b4_value = lookup(i, j, north_boundary)\n",
        "\n",
        "            # create boundary dictionary\n",
        "            boundary = {'loc': np.array([2, 4]),\n",
        "                        'type': np.array([b2_type, b4_type]),\n",
        "                        'value': np.array([b2_value, b4_value])}\n",
        "\n",
        "            # create non-boundary dictionary\n",
        "            non_boundary = {'loc': np.array([1, 3])}\n",
        "\n",
        "            # create well dictionary\n",
        "            well = {'condition': well_condition[i, j],\n",
        "                    'value': well_value[i, j],\n",
        "                    'rw': well_rw[i, j],\n",
        "                    'Gw': well_Gw[i, j]}\n",
        "\n",
        "        if bound_loc[i, j] == 1:\n",
        "            \"\"\"\"\"\"\"\"\"\"\"\n",
        "            MODIFY THIS W.R.T. THE GRID BLOCK LOCATION \n",
        "            \"\"\"\"\"\"\"\"\"\"\"\n",
        "\n",
        "            # get the TYPE and VALUE of the boundary from the BOUNDARY DICTIONARY\n",
        "            b1_type, b1_value = lookup(i, j, west_boundary)\n",
        "\n",
        "            # create boundary dictionary\n",
        "            boundary = {'loc': np.array([1]),\n",
        "                        'type': np.array([b1_type]),\n",
        "                        'value': np.array([b1_value])}\n",
        "\n",
        "            # create non-boundary dictionary\n",
        "            non_boundary = {'loc': np.array([2, 3, 4])}\n",
        "\n",
        "            # create well dictionary\n",
        "            well = {'condition': well_condition[i, j],\n",
        "                    'value': well_value[i, j],\n",
        "                    'rw': well_rw[i, j],\n",
        "                    'Gw': well_Gw[i, j]}\n",
        "\n",
        "        if bound_loc[i, j] == 2:\n",
        "            \"\"\"\"\"\"\"\"\"\"\"\n",
        "            MODIFY THIS W.R.T. THE GRID BLOCK LOCATION \n",
        "            \"\"\"\"\"\"\"\"\"\"\"\n",
        "\n",
        "            # get the TYPE and VALUE of the boundary from the BOUNDARY DICTIONARY\n",
        "            b2_type, b2_value = lookup(i, j, east_boundary)\n",
        "\n",
        "            # create boundary dictionary\n",
        "            boundary = {'loc': np.array([2]),\n",
        "                        'type': np.array([b2_type]),\n",
        "                        'value': np.array([b2_value])}\n",
        "\n",
        "            # create non-boundary dictionary\n",
        "            non_boundary = {'loc': np.array([1, 3, 4])}\n",
        "\n",
        "            # create well dictionary\n",
        "            well = {'condition': well_condition[i, j],\n",
        "                    'value': well_value[i, j],\n",
        "                    'rw': well_rw[i, j],\n",
        "                    'Gw': well_Gw[i, j]}\n",
        "\n",
        "        if bound_loc[i, j] == 3:\n",
        "            \"\"\"\"\"\"\"\"\"\"\"\n",
        "            MODIFY THIS W.R.T. THE GRID BLOCK LOCATION \n",
        "            \"\"\"\"\"\"\"\"\"\"\"\n",
        "\n",
        "            # get the TYPE and VALUE of the boundary from the BOUNDARY DICTIONARY\n",
        "            b3_type, b3_value = lookup(i, j, south_boundary)\n",
        "\n",
        "            # create boundary dictionary\n",
        "            boundary = {'loc': np.array([3]),\n",
        "                        'type': np.array([b3_type]),\n",
        "                        'value': np.array([b3_value])}\n",
        "\n",
        "            # create non-boundary dictionary\n",
        "            non_boundary = {'loc': np.array([1, 2, 4])}\n",
        "\n",
        "            # create well dictionary\n",
        "            well = {'condition': well_condition[i, j],\n",
        "                    'value': well_value[i, j],\n",
        "                    'rw': well_rw[i, j],\n",
        "                    'Gw': well_Gw[i, j]}\n",
        "\n",
        "        if bound_loc[i, j] == 4:\n",
        "            \"\"\"\"\"\"\"\"\"\"\"\n",
        "            MODIFY THIS W.R.T. THE GRID BLOCK LOCATION \n",
        "            \"\"\"\"\"\"\"\"\"\"\"\n",
        "\n",
        "            # get the TYPE and VALUE of the boundary from the BOUNDARY DICTIONARY\n",
        "            b4_type, b4_value = lookup(i, j, north_boundary)\n",
        "\n",
        "            # create boundary dictionary\n",
        "            boundary = {'loc': np.array([4]),\n",
        "                        'type': np.array([b4_type]),\n",
        "                        'value': np.array([b4_value])}\n",
        "\n",
        "            # create non-boundary dictionary\n",
        "            non_boundary = {'loc': np.array([1, 2, 3])}\n",
        "\n",
        "            # create well dictionary\n",
        "            well = {'condition': well_condition[i, j],\n",
        "                    'value': well_value[i, j],\n",
        "                    'rw': well_rw[i, j],\n",
        "                    'Gw': well_Gw[i, j]}\n",
        "\n",
        "        if bound_loc[i, j] == 0.0:\n",
        "            \"\"\"\"\"\"\"\"\"\"\"\n",
        "            MODIFY THIS W.R.T. THE GRID BLOCK LOCATION \n",
        "            \"\"\"\"\"\"\"\"\"\"\"\n",
        "            # create boundary dictionary\n",
        "            # no boundary\n",
        "            boundary = {'loc': None,\n",
        "                        'type': None,\n",
        "                        'value': None}\n",
        "\n",
        "            # create non-boundary dictionary\n",
        "            non_boundary = {'loc': np.array([1, 2, 3, 4])}\n",
        "\n",
        "            # create well dictionary\n",
        "            well = {'condition': well_condition[i, j],\n",
        "                    'value': well_value[i, j],\n",
        "                    'rw': well_rw[i, j],\n",
        "                    'Gw': well_Gw[i, j]}            \n",
        "\n",
        "        \"\"\"\"\"\"\"\"\"\"\"\n",
        "        DO NOT MODIFY THIS\n",
        "        \"\"\"\"\"\"\"\"\"\"\"\n",
        "\n",
        "        # calculate INTER-BLOCK transmissibilities\n",
        "        T = transmissibility2d(dx, dy, dz, kx, ky, mu, B)\n",
        "        T = np.array(T)\n",
        "\n",
        "        # calculate BOUNDARY transmissibilities\n",
        "\n",
        "        if bound_loc[i, j] == 0.0:\n",
        "            # Interior blocks. All side are non-boundary, so flow are all inter-block\n",
        "            T = T\n",
        "            boundary['T'] = None\n",
        "\n",
        "        else:\n",
        "            # Boundary blocks.\n",
        "            bound_T = []\n",
        "            for k in range(len(boundary['loc'])):\n",
        "                a = boundary['loc'][k]\n",
        "                _ = transmissibility2d_boundary(boundary['loc'][k], boundary['type'][k],\n",
        "                                                dx, dy, dz, kx, ky, mu, B)\n",
        "                T[a - 1] = _\n",
        "                bound_T.append(_)\n",
        "\n",
        "            boundary['T'] = bound_T\n",
        "\n",
        "        # calculate LHS coefficients\n",
        "        px_min, px_plus, py_min, py_plus, p = lhs_coeffs2d_welltype(bound_loc[i, j],\n",
        "                                                                    well, T, mu, B, \n",
        "                                                                    solver='incompressible') \n",
        "\n",
        "        # calculate potential term \n",
        "        # Equals 0. the grid block is not elevated\n",
        "        potential_term = 0\n",
        "\n",
        "        # calculate RHS constants\n",
        "        rhs = rhs_constant2d_welltype(boundary, well, 0, dx, dy, \n",
        "                                      dz, kx, ky, mu, B, solver='incompressible')\n",
        "\n",
        "        # fill in LHS matrix\n",
        "        lhs_mat = fill2d_lhs_mat(bound_loc[i, j], block_index[i, j], xi, lhs_mat,\n",
        "                                 px_min, px_plus, py_min, py_plus, p)\n",
        "\n",
        "        # fill in RHS matrix\n",
        "        rhs_mat = fill2d_rhs_mat(block_index[i, j], rhs_mat, rhs)\n",
        "\n",
        "\"\"\"\"\"\"\"\"\"\"\"\n",
        "PRESSURE SOLVER\n",
        "\"\"\"\"\"\"\"\"\"\"\"\n",
        "p_sol = np.linalg.solve(lhs_mat, rhs_mat)\n",
        "\n",
        "p_sol = p_sol.T.reshape(-1)\n",
        "\n",
        "\"\"\"\"\"\"\"\"\"\"\"\n",
        "DISPLAY PRESSURE\n",
        "\"\"\"\"\"\"\"\"\"\"\"\n",
        "\n",
        "# display well\n",
        "xsc = np.array(well_loc)[:,0]\n",
        "ysc = np.array(well_loc)[:,1]\n",
        "\n",
        "p_sol = p_sol.reshape(-1, xi)\n",
        "plt.figure(figsize=(10,6))\n",
        "plt.imshow(p_sol, origin='lower')\n",
        "plt.colorbar()\n",
        "# plt.scatter(xsc, ysc, color='black')\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbQAAAFlCAYAAACDVh3MAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3df4xlZZ3n8c+3blU3CCq/XJbtZgd2ZGNwsoOTXmTi/sFi1FaJMIljUNftTMiym2AWM7qOuH8wOsNGk1mRyTomvcKKE1dk0Vk6xl3CIsQ1WcFGEQVibH+FJgiDgAJKN1313T/uaS3a+/123afOPffc57xfSaWrzrnPeZ576lY/95zzOd9r7i4AABbd0rwHAABAG5jQAABVYEIDAFSBCQ0AUAUmNABAFZjQAABVWO6ys9Hxx/nyySdOXrlmExdbdlfB2uTFbbcJ1wXbStuk/QQrs+fT2fbiRvFzLbglJGtTMLb0uZY0mvcY2r7LJtjf5d2U/M6LO6vCc3pWB/3A5P8AZ+AN//I4/9kTq8Xt77nvwK3uvrPFIbWm0wlt+eQT9Q//4xUT1y39cvLB4ui5+Pe8/KvJ60YH4jEsBetGB+O/qmh7WZul54M2zydtgu0tPR/PnEvZ9oJ2tpptL1h3KG5jq5P/OGw1+Z/qUPAHtZb0U9BG0XPNJs7g+UiSh9srGUPSZi2YaLJxR/uhoE3eTzZ5J88pbFIyCU7fT9EbrA7c5bd32t/jT6zqrlu3F7dfOe0Hp7Q4nFZxyhEAUIVOj9AAAPPmWi05wl0ATGgAMCAuaa3SC5dMaAAwMGtZom2B9WdCC67meXKVby0YvWUBnijlGKQsx+uixFkWTIreAZWEmZKdkMYzg8Wr8Rh8KUibLsdjWDo0eZ1nQZJR0E8SJPHl0eQVScglCj1kwRh5/GcRtsuCKQVhjSiUEQVw0u0VhDjSV2myvTBMkoV94o1lo5i6n65MXfC9u4Bj9fozoQEAZs7lWu1p4nOzmNAAYGC4hgYAWHguaZUJDQBQg1qP0LixGgBQhd4coflS8I4hSN1Jko+C8kBBgm68brrlkrQWJfKW43c5a1FtyiQhthalJgvfTLlNfr+SltIK2mRpyiggGiUZJUmHJq/zJOVoUYKtoE1YwkoqKhWV/V7j9GGS/AuTte32E9ZyLCx9FScWC/f3tErKaGUKbkCeNrNoz3d7XOESoRAAQB3mf3PDbDChAcCAuJxQCACgAp6epV9ohEIAAFVgQgOAARkXJy7/2igzG5nZt8zsS83PZ5rZXWa2z8w+b2ZbmuVbm5/3NevPWLeNK5vl3zOzNxytTyY0ABgU0+omvqZwhaQH1/38UUnXuPvLJT0p6dJm+aWSnmyWX9M8TmZ2tqRLJL1S0k5Jf2NmSR6962topiT+HcS4s/0XpdyTNlGx46wIchj1z+rEBm2iOP+4UVCQNqt7m7zAoqe0lr2PCW6fsCS2vxTt0yxOHxVBzmLXq0FH6SdWF7RJVoW3AZTE3Avi9EW3BxRE/dN+MiVjiLRdaLivUfVnui1O7Gr/7oYjmdl2SW+WdLWkPzUzk3SBpHc0D7lB0p9L+qSki5rvJelmSf+lefxFkm509wOSfmRm+ySdK+n/Rf1yhAYAmMYpZrZ33ddlEx7zcUnv12/eHp4s6Sl3P9T8vF/Stub7bZIekqRm/c+bx/96+YQ2E5FyBICBmfLU4ZEed/cd0Uozu1DSY+5+j5mdv5mOpsWEBgADMi5OPNPTnK+R9BYze5OkYyS9RNK1kk4ws+XmKGy7pIebxz8s6XRJ+81sWdJLJf1s3fLD1reZiFOOADAwa27FX0fj7le6+3Z3P0PjUMdX3P2dku6Q9NbmYbsk3dJ8v6f5Wc36r/j4ovQeSZc0KcgzJZ0l6e6sb47QAGBAOjhCi/yZpBvN7C8lfUvSdc3y6yT9bRP6eELjSVDufr+Z3STpAUmHJF3unkXx+jShRfs3Sde5RVVxW+w/WZemKcPnU9BPUqA5S25FSctkl4apyexQfi0oaGxRwWkpjFnlycgokTd9AWtby2KtBUnLkiRhVq6hzaLBJQWNM0mbonRky2OY67ZKRP+PVcDd75R0Z/P9DzVOKR75mOck/XHQ/mqNk5Ib0p8JDQAwcy7TaqVXm5jQAGBgNnItbBExoQHAgMzxGtrMMaEBwKCYVrPSSAuszmcFABgcjtAAYEDG1fbrPJZhQtuINhPCBVH/rE2W+I3i/tnZBguebFYEObx7IrljxIJGWRHkMMqeRf2jJHsSL7f0Voj2igZHt1WkYyiJ4GfFltssJqzktoK2q+ES298UrqEBABaeO9fQAADoNY7QAGBgsksIi4wJDQAGZHwfWp0n55jQAGBQ6r2G1p8JLQoaJRG/MBCXhJbabGNZeixYV9Sm7YBY+loO0ofJDvLo9EWW3oqKE6f7JxhbUgQ5CrClScaS1GSaPiwYQ7B/0tdCQWIx2t/Z2DorJpy8Fkqkz2meOk451hzbr/NZAQAGpz9HaACATqxSnBgAsOj4+BgAQDXWCIUAABZdzbH9Op8VAGBwen+ElseUgzYl0fiskO6hyR0tHYrbROuWskh4GPWfPhKeSgsaT98oivSnye+oOHFQUHm8vSjKPv3tAVmbtm8DiOLnWdS/q9sDwtsn4hat3zpQpCDS39PQflhEfGb9yQiFAADqUOt9aExoADAg7qq2UkidzwoAMDgcoQHAoBjV9gEAi89V7ynHbic0V1xsOFweb66oAHCUPnw+bjM6GCw/EA9uFGzPClKOpYqKsUa/huT1H6W00ucTFolOo5FTN9EoGFtJUi9rlyRRywpiR2nB6dtkBZ+j55rtn9aTlkPXcXFiqd770DhCA4ABcZnWKo3t1zlNAwAGhyM0ABgYTjkCABaei+LEAIAqmFaJ7QMAFh1HaJLMbCRpr6SH3f1CMztT0o2STpZ0j6R3uXsQcN+ENF49fVw8KkI8Ohh3tPzc5HUrz8YdLT8bdJQV3w3W+ShskhY2zdqFbaIIccEbuuxvJoyYp5WTg20V3B5QWnw3TJ9nxY4LblEIiyentwcU3PoS/R6yNsmtJ3Fx4ulfQK0H8eYQj9+Qng5rEU0zTV8h6cF1P39U0jXu/nJJT0q6tM2BAQBmY7U57Vjy1WcbmtDMbLukN0v6VPOzSbpA0s3NQ26QdPEsBggAaI+7ac2Xir/6bKOnHD8u6f2SXtz8fLKkp9z9cN2N/ZK2tTw2AMAM1Fr66qjPyswulPSYu99T0oGZXWZme81s7+ozz5ZsAgCAo9rIEdprJL3FzN4k6RhJL5F0raQTzGy5OUrbLunhSY3dfbek3ZK09Xe2U9ANAObIpeFW23f3KyVdKUlmdr6k97n7O83sf0h6q8ZJx12SbpnhOCdKC8xGbaKCxkEoUZKWfzW5o2MffS5ss/SdH0xe8bunh22eP+lFE5evHhPHFdOiwQWJxTBZ1vLrP0wzlvxOS94mpcm/6dtlSctwc2kysiAtGBYGjptYlJrM0pRJsjbanh2M/8B8ZfKLOHz9SnER7ZIk45z/by8a86bYcE85Jv5M0p+a2T6Nr6ld186QAACzMr4PzYq/+myqG6vd/U5Jdzbf/1DSue0PCQAwS7XWcqzzWQEABofSVwAwIDV/HhoTGgAMzFqlJ+eY0ABgQNylVY7QZizKFieRVl8KCqtmbYJnvLo1bKKDxwex4m3Hhm2WTzo73mDAR5PHvbZcVoA4LHacvZbn/TrvLLZfUOQ33V68Ki5OXNCmYAx58eaowHdWgDheFY17dWvBEUHLcfbe/h8+h3HVesqxzuNOAMDg9OcIDQAwc+NQSJ3HMkxoADAwff8YmFJMaAAwIIcrhdSozuNOAMDgdHuEZpo6zeijOFIVpgJX4iGEnx6fvWMJxry6JX4/MDpu8rqi9FjbhYE7KoYaFtgt1WoCcvoiv+n22k5nFhQannZb6fZKCienY6jziKAN3V/OqvcaWp3PCgAQWpMVfx2NmR1jZneb2bfN7H4z+1Cz/NNm9iMzu7f5OqdZbmb212a2z8zuM7M/WLetXWb2/eZr19H65hoaAAxIBzdWH5B0gbs/Y2Yrkr5mZv+rWfcf3P3mIx7/RklnNV+vlvRJSa82s5MkXSVph8bnGe4xsz3u/mTUMRMaAAzMLE85urtLeqb5caX5ys5TXyTpM027r5vZCWZ2mqTzJd3m7k9IkpndJmmnpM9FG+KUIwCgVWY2MrN7JT2m8aR0V7Pq6ua04jVmdrg+0zZJD61rvr9ZFi0PMaEBwIAcrra/iQ/4PMXM9q77uuy3+nBfdfdzJG2XdK6Z/Z6kKyW9QtI/l3SSxh8S3SpOOQLAwGwk3JF43N13bOSB7v6Umd0haae7/1Wz+ICZ/TdJ72t+fljS6euabW+WPazxacf1y+/M+uvPhBbs3+xU71pUmHel3cKz0e0BS9ntAYeC/pNCsdG6knj3UddNqSwuPv2+bjuWXtRP0fbKivnG24v+IEq2laxsc2yZtm/hqElHt9EcNusbq83sZZKebyazYyW9TtJHzew0d3/EzEzSxZK+2zTZI+ndZnajxqGQnzePu1XSfzKzE5vHvV7jo7xQfyY0AEANTpN0g5mNNL6sdZO7f8nMvtJMdibpXkn/rnn8lyW9SdI+Sb+U9CeS5O5PmNlfSPpG87gPHw6IRJjQAGBgZpxyvE/SqyYsvyB4vEu6PFh3vaTrN9o3ExoADMlvwh3VYUIDgAFxbToU0ltMaAAwMLUeoXEfGgCgCt0foUVvDKLlyZTrQWw/SwivhSvjdyxhgjq6bUCSrUbLkyrmUWw/ifqXxPbbj5h3sK3C7XU3hqw6/fRj6KrafqT4kxJKPt0gHEPZEBZN1wdLNX8eGqccAWBgmNAAAAvvcOmrGjGhAcDA1JpyJBQCAKgCR2gAMCTONbSZ86Ug0pQdQ0Ypx7X4l7UWPuOkoPHS5O15UIBYkiwYW1Fx4oI2kuKnlLyYSxKLbaYc834KNliQFkx7mXMKtPXUZthPu//h9TqxOO+xkXJsTW8mNABAN2qd0LiGBgCoAkdoADAgxPYBANVwJjQAQA1qvQ+NCQ0ABsSJ7XegoDhxtM5HWQ53ckdZ+t2CzPGSJQWNowh+ULRYSmL7hQV7w0h/yxH8olh6UZw+qhIdt4m213Y0v7MCwF3F9lseQ6bXkf4OVDq3zEV/JjQAQCe4hgYAqAApRwBAJThCAwAsvJpLX1EpBABQhe6P0KZ9Y5A8PixoPMo6iSJVSZtg1VoSzwoThkGhYylJRhYWJ/aCtyvRU0r7aTEVWJrojLRaGLjlMYSpzcJ+pu2/WFepxKGkH7s+WPLCv9kFwClHABgYbqwGACw8V72hEK6hAQCqwBEaAAwK96EBACpBKAQAUIVar6H1fkLzLBofFAdO4+rh5qYvaJymrsMIfkHUfy2J+rdcnDjanhXt01gYZW+5+G6btxQUj6HFfkrG0Pqb8T68u+/DGFrS9dziXu+ERigEAFCF3h+hAQDaRSgEAFAFQiEAgCrUeg2NCQ0ABsRlw53QzOwYSV+VtLV5/M3ufpWZnSnpRkknS7pH0rvc/WDxSMKquMmOjyIt2fF0UBw4PwIP1ibpw1DWJFgXFmHW0YoTT95gSQHgNE1ZUAA4TFMWnArp8vRJUZKwo0LDre6GRT0ltYjjrnNumYuNpBwPSLrA3X9f0jmSdprZeZI+Kukad3+5pCclXTq7YQIA2uKb+Oqzo05oPvZM8+NK8+WSLpB0c7P8BkkXz2SEAID2NPehlX712YbuQzOzkZndK+kxSbdJ+oGkp9z9UPOQ/ZK2zWaIAIBWVXqItqEJzd1X3f0cSdslnSvpFRvtwMwuM7O9ZrZ39elnC4cJAEBuqkoh7v6UpDsk/aGkE8zscKhku6SHgza73X2Hu+8Yvfi4TQ0WALB5gz3laGYvM7MTmu+PlfQ6SQ9qPLG9tXnYLkm3zGqQAID2jOs5ln312UbuQztN0g1mNtJ4ArzJ3b9kZg9IutHM/lLStyRdN5MRZgWAgwyzBXF1SfLgJLAlHcW/wyzLHmwvi9lHNZCzqHbylsSCV1/6oiyI00fv2tq+PSBSEvUvvRbQ5h900S0K7XU/ow1ial0XJ9aAb6x29/skvWrC8h9qfD0NALAoXN2X+O8I1fYBAFWg9BUADEzfr4WVYkIDgKFhQgMALL7+x+9L9WdCi/Zv9k4iatN2AeCgUfaasLUgYZi2CfopSAumSraXtilIUwbKihNP/8dZlIxMB9FJk2rfWc9T66+FKc1lbqn0dUQoBABQhf4coQEAZs8HfB8aAKAynHIEANTBNvF1lC2bHWNmd5vZt83sfjP7ULP8TDO7y8z2mdnnzWxLs3xr8/O+Zv0Z67Z1ZbP8e2b2hqP1zYQGAGjTtB8KfamkJ5vl1zSPk5mdLekSSa+UtFPS3zQlGENMaAAwNDP8PLSCD4W+qPlZzfrXmpk1y2909wPu/iNJ+3SUcovdX0OLMrIlFymjCH6Sww0LF2dFg6NIf8mY09sQon6y7WXVmwvG0Oa59ZLixEmT+NeaVU6evp8ilV6TQAcqjO03R1L3SHq5pE8o/1DobZIekiR3P2RmP5d0crP86+s2e9QPkiYUAgBDsvnixKeY2d51P+92990v6MJ9VdI5zUeP/Z2m+FDozWBCA4CB2WQtx8fdfcfG+vGnzOwFHwrdHKWt/1DohyWdLml/86HRL5X0s3XLDws/SPowrqEBAFpT8KHQe5qf1az/irt7s/ySJgV5pqSzJN2d9c0RGgAMzWyvoU37odDXSfpbM9sn6QmNk41y9/vN7CZJD0g6JOny5lRmiAkNAIZmhpVCpv1QaHd/TtIfB9u6WtLVG+27/xNaVjk0+qWkwb/J2wvTj+NGkwUFiCXJg+2lT6cglWjJyfCSUGC4T9suglzQJlrVaUiMNOPi6mu1pzmMa94FmWel/xMaAKA9G7yfbBERCgEAVIEjNAAYFJvTh7DNHhMaAAxNpaccmdAAYGgqndC4hgYAqEJ/jtBKTum2mOOO4vzjzQUbzN4OBHH6tORM28WOS4oTt1nNt6hNy/ugxyq9jNFfPc2qz+V10M9dsWn9mdAAALO3+eLEvcWEBgAD09OD1U1jQgOAoal0QiMUAgCoAhMaAKAKnHIEgIHhGlofRb+VLMHTZqQ/66akan3Bx8gWffJsV9H4ojYFjUpuXcg+XKGrAFgfgmZ9GENHehvsm0tsv687Y3MWe0IDAEyHavsAAPQbR2gAMDSVHqExoQHAwBAKAQDUgQmth8KgTvLbitI9WegnWpd102YyMtNRkrDo9d/RH423ndhq++1rbanJHgfkfBEPPRZxzD212BMaAGB6lc6hTGgAMCDm9R4UMqEBwNBwYzUAoAqVHqFxYzUAoAocoQHAwHANbV6yU70FhWfDRiXnlNuO+pe8yjp6YRbdUtCZgtsQuvyLbnPX9fnXIM3/f8q+759J5lKceA59dqD/ExoAoD0Vpxy5hgYAqAJHaAAwNJUeoTGhAcDQMKEBAGpQ6zW0/kxo0R7O0nUFScKigsZhk5aTkSU6SkiFxZZrVPLX3off66KOu6ifil6Pi5jM7ClCIQCAKvTnCA0A0I2KDnDXY0IDgCGp+D40JjQAGBomNABAFSqd0AiFAACqcNQjNDM7XdJnJJ2q8by+292vNbOTJH1e0hmSfizpbe7+5OyGOoXWCxpPuzFI6kckvO2LBdEYOhqb9TzOb23u7x7E2Yv299SddPv/iHXfZWc2coR2SNJ73f1sSedJutzMzpb0AUm3u/tZkm5vfgYA9J1v4qvHjjqhufsj7v7N5vunJT0oaZukiyTd0DzsBkkXz2qQAICWNCnH0q8+m+oampmdIelVku6SdKq7P9Ks+qnGpyQntbnMzPaa2d7VZ57dxFABAIhteEIzs+MlfUHSe9z9F+vXuXt4MOruu919h7vvGB1/3KYGCwBowVBPOUqSma1oPJl91t2/2Cx+1MxOa9afJumx2QwRANCqSie0jaQcTdJ1kh5094+tW7VH0i5JH2n+vWUmI8xO2nZVHLjVZGSH+lwotiitV9Jm8tiK02vRcy15KWb7rWh70YqCNOX03TdjKHkttPe/ZCepxJbNY8x9vxZWaiM3Vr9G0rskfcfM7m2WfVDjiewmM7tU0k8kvW02QwQAtGqoE5q7f03xG7bXtjscAADKUPoKAIZkAa6FlWJCA4CBGfI1NABATZjQAAA14Aitj6LfSkmcP+2n3c21agEj82mTrJ8Wn2tpZL6raHw4vqSfkvh31E9plLwktl/WZuomrfbfpnn337akoP2fS/o3kv6+eegH3f3LTZsrJV0qaVXSv3f3W5vlOyVdK2kk6VPu/pGs78We0AAA05vtHHq4oP03zezFku4xs9uadde4+1+tf3BT7P4SSa+U9I8k/R8z+6fN6k9Iep2k/ZK+YWZ73P2BqGMmNAAYkhmnHJsav4803z9tZocL2kcuknSjux+Q9CMz2yfp3GbdPnf/oSSZ2Y3NY8MJjQ/4BIABsU1+TdXXCwvaS9K7zew+M7vezE5slm2T9NC6ZvubZdHyEBMaAGAapxz+BJXm67JJD5pQ0P6Tkn5X0jkaH8H957YHxilHABiazZ1yfNzdd2QPmFTQ3t0fXbf+v0r6UvPjw5JOX9d8e7NMyfKJ6pzQ2i5o3PYYpt5We5s66vbaTJwVJAlL0oclycji4rtLa1OPoSRJ2GZasGxsSZpyqlGNLaXpzDaLE3eTGFxq8W8y2zezMssuo4L2Znbaus/Q/CNJ322+3yPpv5vZxzQOhZwl6W6NX2pnmdmZGk9kl0h6R9Z3nRMaACA22zk0Kmj/djM7p+n9x5L+rSS5+/1mdpPGYY9Dki5391VJMrN3S7pV49j+9e5+f9YxExoADM1sU45RQfsvJ22ulnT1hOVfztodiVAIAKAKHKEBwJA4pa8AALVgQgMA1IAjtCHr6rffdtHggUTw09hzy5H5qF3WZmmpIBpfcLtBtB9K+sli6elzLRjDtNvK9LmOeJv7ZtMqndAIhQAAqsARGgAMDKccAQCLb8bV9ueJCQ0AhqbSCY1raACAKtR5hNZ2AeJoe31IP5YkwYLUXdbXoiYWSxKGUZusXba9UUmbknEH60ZBQWUpfmllv4eSQsOl25vW0gIeenRdnNjENTQAQC2Y0AAANTCvc0ZjQgOAIak45UgoBABQBY7QAGBgCIUAAOrAhNZDbcfz2+y/B2+B0th82CYYd/ZUo4h51k8QJe8qgr9UEGWX4gh8SQQ/i9OHEfy0n8nby2Lhyzb98ymJ4Kdtgv9dl4Kxta3r2PyRst/prPTgv6eZWOwJDQAwvUonNEIhAIAqcIQGAEPinHIEANSCCQ0AsOio5ThPJUnGPvyywoLGaaPp+ykJeqbxw2BxlnKMlmcpvjCxGPcTJROzYsIlBXujVGLWrqQ48XI6hsnrolRi2ibpZ9lWp25TklhcXprcT2Ye6b95GHWU5hyC/k9oAIB2UcsRAFCDWg9+mdAAYEgqLk7MhAYAA1PrZTturAYAVIEjNAAYGk45zlib8fySX1ZXhY6zq7FL0RiSNh0VSM5i6VE8P43Th0WDszaT+8ni3WGbgmh+1lcWc18ZTY6sZ3HtaHtZbH/L6NDUbaI4/dYkZr+SrGs1tq943Kstnlxam3OR83kURyYUAgBYfC5i+wCAOtR6hEYoBABQBY7QAGBoKj1CY0IDgAGhODEAoA7uhELmJtvvYWy/5Qr9bf7us3LyUT+jZHsFL8w0gh9W2y9pE48hju1PH5nP2qyMotj+9P2MtxdUp0+i8SVtSiL4W4M2W5YmL5ekrcG6aPnR1h2z9PzE5StBVf9Sqz758v9a8lESbbbJTHsbANX229P/CQ0A0CpOOQIA6sCEBgCoAUdoAIDF55LW6pzRuLEaAFCFox6hmdn1ki6U9Ji7/16z7CRJn5d0hqQfS3qbuz85u2EGStKMQaDI1pJttdjGk8Sir0QdxW1af0sSnYtICwBPnz4sKTS8XJBYDIv8ZinHLDUZFNndEiQZJWk5SPhlbaJkYlY0OEofHjs6GLZ5UbDu+NFzYZsXL8XrXrR0YOLy45biMawGL/Dn1lbCNs/5lonLn0/+wKJ1a0mSMRpb2mbKZGRU0Hmm6jxA29B/h5+WtPOIZR+QdLu7nyXp9uZnAMACMC//6rOjTmju/lVJTxyx+CJJNzTf3yDp4pbHBQCYlcM3V5d89VhpKORUd3+k+f6nkk6NHmhml0m6TJJGJ59Q2B0AoC19P9IqtekrMO7uSs7Iuvtud9/h7jtGxx+32e4AAJiodEJ71MxOk6Tm38faGxIAYGZ8k189Vjqh7ZG0q/l+l6Rb2hkOAGCWxtX2vfirzzYS2/+cpPMlnWJm+yVdJekjkm4ys0sl/UTS22Y5yE4k9UHt0OTo7uhAnKdf+cXkdQdfGr8gVoNqvnkKOHmBBSfKLdneUlgAOIvtB3H6IGYvxfH8KJovxXH6LIIfrYsKBo/HlhUNDiL4ScHeLUHUfjmJ4B87mlzkNysMHLV5URKZf3EQz89i+/9g+RfhujOWfzZx+Tlbt4Zt7j/4q4nLf7p6fNjmqdXJly+e8yTqH9wGkEX9V4P3/FlsP97W5L/vpXkUJ660HvJRJzR3f3uw6rUtjwUAgGKUvgKAgen7qcNSlL4CgCGZcSjEzE43szvM7AEzu9/MrmiWn2Rmt5nZ95t/T2yWm5n9tZntM7P7zOwP1m1rV/P475vZrqjPw5jQAGBQNnFT9caO7A5Jeq+7ny3pPEmXm9nZiitMvVHSWc3XZZI+Kf26xOJVkl4t6VxJVx2eBCNMaAAwMLMsfeXuj7j7N5vvn5b0oKRtiitMXSTpMz72dUknNLeDvUHSbe7+RFMr+Db9dhnGF1jsa2jR3s2KFgdTuC8nv6mgzaGkzaEXBWMYxW3CMSzHkSRLthets+RVaWGh4emLE3dVaDhLLIZtkoRhlpqM0oxRknG8vcnrskLDK0HyLUs5rgRFkKPlUpywGyUxuINJKvCnqy+ZuPzOX01OYErSqianGbPEYpREHXlSWNom77ss1RoVGl7NosJTmktx4s05xcz2rvt5t7vvnvRAMztD0qsk3c1bYDUAAAl7SURBVKW4wtQ2SQ+ta7a/WRYtDy32hAYAmN7mQiGPu/uOoz3IzI6X9AVJ73H3X9i625Lc3S17Z12IU44AMCQu2Vr510aY2YrGk9ln3f2LzeKowtTDkk5f13x7syxaHmJCA4ChmWEoxMaHYtdJetDdP7ZuVVRhao+kf92kHc+T9PPm1OStkl5vZic2YZDXN8tCnHIEALTpNZLeJek7ZnZvs+yDiitMfVnSmyTtk/RLSX8iSe7+hJn9haRvNI/7sLsf+VFmL8CEBgBDM8Mcirt/TQrqfE2oMNV8Ysvlwbaul3T9RvtmQgOAgam1Ukh/JrQo8JIk8ON4fvLLilYlVxO95JcfDS2Jv0djsKRgbxbbj4oDLyXbCyP4SZso5t52oeEonp+2CQsDTx/NH6+bvtBwFM/PIvhblybH3LMIfrQui6VHsoK9UZFfKS7am407Gl8UmZfiQr+ZLcEYVpP/L5aC2xfiPSCtTRlNsHnE9pnQAAALz1VttX1SjgCAKnCEBgADYur/B3WWYkIDgKFhQgMAVIEJrSXTh5OSbbX4SykZV1bktyDlGBUGjpZL0lKS1itJLEbrSgoNLydFg6PtdVVouCTJON5eQaHhgmRkNO6l5PcQpQWjAsSlnvfp/9tYTS7XRwWFS8adtYkSmFkh5lHB/wtZanKS8P+KWSEUAgBAv3HKEQAGhlAIAKAOTGgAgMW3sar5i4hraACAKnCEBgBD4qr2CK3/E1oWaQ0izEUx2OwWgGB72SeIW1RoOG0TxJfT2P7066LCwNm6kkLDKwX9ZJH5kjZRNL6kjSStBLHw9NaBoChuWrA3KoqbtFkqKHI7arkwblSYd6R43OG2kuLEbUf62xQVNI7MpThxpbH9/k9oAIBWkXIEANSh0gmNUAgAoAocoQHAkLiktTqP0JjQAGBQ6r0PrT8TWpRMbLsAcLS9JBkZFg2OmxQlFqPnkxUgTgsNB881295KVJw4G0NJoeEoLVjQJkslRmnGrE1RoeEkfbg1KIScJRajfqICxOM20/fTdvJv2oRfuq2CsbWd2sz2d1tK0qmbxoQGAKhCpRMaoRAAQBU4QgOAISEUAgCog0vBh6ouOiY0ABgarqEBANBf/TlCKyk0HMXpC6LxJUWDs7FF8fwsMh+2ScaWFxqe3G45icZHEfzlgkLDacHeqAhyEpPeMpocS8/aRPH8kmj+uK9oe5PHJsWx+ayfsE1a0DgoRp3EwkuKIGdR9uj1E/UjxfH8dNwlxYlbvKWgzTh/58WJuYYGAKhGpaccmdAAYGiY0AAAi6/e0leEQgAAVeAIDQCGxCWtcR/aTBUVGi5IOYaJwSRJGCcWS9q0W0w4SjKO10Xpw+kTi2nKMUh8bUkLDUdFg5OUY5AKzFOJ0ycwt5QkFgsKDadtSvoJ20z/fNIiyGnSst3UZCRKLJZsq6SgcZuJySzIPTOVnnLszYQGAOgIExoAYPF5tfehEQoBAFSBIzQAGBKXnOLEAIAqVHrKkQkNAIaGUEgLTHFGNSxOPH2h4ayYbxTpzyP4kw/P0zZRkdaCCH7eJolXB+uyNlHMPW0TrMvi7yXFiaMCwFEBYklaCbaXt4nXRWMoitO3fXtAFMFPIuZRpH9L0k9UTHjc1+TXcHbrQBS1bztO31Wkf1qdFyeuGEdoADAk7txYDQCoBKccAQA1cI7QAACLj2r7AAD02qaO0Mxsp6RrJY0kfcrdP1K+sWBxMuVG6cO0OHGUJBwlHxEfpQXTgsZRWnD6YsJZP1kx35WgOHCW9ipJH0aJwaiYcNYmSiVmbaLk4Xh70xdBLkkSHrP0/NRt8n4mP6cswbslaFPST5buazuxWFJoONpelsCMt9Xu6bdpk5FZknsmXNyHdiQzG0n6hKTXSdov6RtmtsfdH2hrcACAGaBSyG85V9I+d/+hJJnZjZIuksSEBgA95ZKcI7Tfsk3SQ+t+3i/p1ZsbDgBgptyrPUKbeSjEzC4zs71mtnf16Wdn3R0AYKA2c4T2sKTT1/28vVn2Au6+W9JuSdp65vY6j3MBYIFwyvG3fUPSWWZ2psYT2SWS3tHKqAAAs1PpKUfzTdxgZ2ZvkvRxjWP717v71Ud5/N9L+knz4ymSHi/uvA7sgzH2A/tAGu4++B13f1lXnZnZ/9Z4X5d63N13tjWeNm1qQttUx2Z73X3HXDrvCfbBGPuBfSCxD7B5VAoBAFSBCQ0AUIV5Tmi759h3X7APxtgP7AOJfYBNmts1NAAA2sQpRwBAFeYyoZnZTjP7npntM7MPzGMMXTOz683sMTP77rplJ5nZbWb2/ebfE+c5xlkzs9PN7A4ze8DM7jezK5rlg9kPZnaMmd1tZt9u9sGHmuVnmtldzd/E581sy7zHOmtmNjKzb5nZl5qfB7cP0K7OJ7R1VfrfKOlsSW83s7O7HsccfFrSkfdufEDS7e5+lqTbm59rdkjSe939bEnnSbq8+d0PaT8ckHSBu/++pHMk7TSz8yR9VNI17v5ySU9KunSOY+zKFZIeXPfzEPcBWjSPI7RfV+l394OSDlfpr5q7f1XSE0csvkjSDc33N0i6uNNBdczdH3H3bzbfP63xf2bbNKD94GPPND+uNF8u6QJJNzfLq94HkmRm2yW9WdKnmp9NA9sHaN88JrRJVfq3zWEcfXCquz/SfP9TSafOczBdMrMzJL1K0l0a2H5oTrXdK+kxSbdJ+oGkp9z98CdnDuFv4uOS3i/9+tM1T9bw9gFaRiikJ3wcNx1E5NTMjpf0BUnvcfdfrF83hP3g7qvufo7GBb3PlfSKOQ+pU2Z2oaTH3P2eeY8FddlMceJSG6rSPxCPmtlp7v6ImZ2m8Tv2qpnZisaT2Wfd/YvN4sHtB0ly96fM7A5JfyjpBDNbbo5Qav+beI2ktzS1YI+R9BJJ12pY+wAzMI8jtF9X6W9STJdI2jOHcfTBHkm7mu93SbpljmOZueY6yXWSHnT3j61bNZj9YGYvM7MTmu+PlfQ6ja8l3iHprc3Dqt4H7n6lu2939zM0/vv/iru/UwPaB5iNudxYPW2V/hqY2eckna9xletHJV0l6X9KuknSP9b4Uwje5u5HBkeqYWb/QtL/lfQd/ebayQc1vo42iP1gZv9M48DDSOM3lDe5+4fN7J9oHJA6SdK3JP0rdz8wv5F2w8zOl/Q+d79wqPsA7aFSCACgCoRCAABVYEIDAFSBCQ0AUAUmNABAFZjQAABVYEIDAFSBCQ0AUAUmNABAFf4/bCauF0CrxjYAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 720x432 with 2 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9xDqtxv1o6GH",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 218
        },
        "outputId": "a898412d-3ef4-49b6-c145-1cfa845eab4c"
      },
      "source": [
        "\"\"\"\"\"\"\"\"\"\"\"\n",
        "WELL SIMULATION\n",
        "\"\"\"\"\"\"\"\"\"\"\"\n",
        "\n",
        "# create well dataframe\n",
        "B = np.full(len(well_name), reservoir_input['B'])\n",
        "mu = np.full(len(well_name), reservoir_input['mu'])\n",
        "well_df = pd.DataFrame.from_dict(wells)\n",
        "well_df['well_mu'], well_df['well_B'], well_df['well_Gw'] = mu, B, Gw\n",
        "\n",
        "# print pressure in order (engineering notations: [0,0], [1,0], [2,0] as 1,2,3)\n",
        "p_sol = p_sol.T\n",
        "\n",
        "# FBHP and flow rate solution after pressure in each grid block solved\n",
        "well_report = solution_well2d(well_df, p_sol)\n",
        "well_report = well_df.drop(['well_mu', 'well_B', 'well_Gw'], axis=1)\n",
        "well_report"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>well_name</th>\n",
              "      <th>well_loc</th>\n",
              "      <th>well_rw</th>\n",
              "      <th>well_skin</th>\n",
              "      <th>well_condition</th>\n",
              "      <th>well_value</th>\n",
              "      <th>well_config</th>\n",
              "      <th>fbhp_sol</th>\n",
              "      <th>rate_sol</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>A</td>\n",
              "      <td>[10, 10]</td>\n",
              "      <td>3.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>constant_fbhp</td>\n",
              "      <td>2000.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>2000.000000</td>\n",
              "      <td>-2459.838800</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>B</td>\n",
              "      <td>[40, 40]</td>\n",
              "      <td>3.5</td>\n",
              "      <td>0.5</td>\n",
              "      <td>constant_rate</td>\n",
              "      <td>-600.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1829.524191</td>\n",
              "      <td>-600.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>C</td>\n",
              "      <td>[10, 40]</td>\n",
              "      <td>4.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>constant_rate</td>\n",
              "      <td>350.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1693.723877</td>\n",
              "      <td>350.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>D</td>\n",
              "      <td>[40, 10]</td>\n",
              "      <td>3.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>constant_fbhp</td>\n",
              "      <td>3000.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>3000.000000</td>\n",
              "      <td>-4645.809559</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>E</td>\n",
              "      <td>[25, 25]</td>\n",
              "      <td>2.5</td>\n",
              "      <td>1.0</td>\n",
              "      <td>constant_rate</td>\n",
              "      <td>-150.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>2776.645130</td>\n",
              "      <td>-150.000000</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "  well_name  well_loc  well_rw  ...  well_config     fbhp_sol     rate_sol\n",
              "0         A  [10, 10]      3.0  ...          0.0  2000.000000 -2459.838800\n",
              "1         B  [40, 40]      3.5  ...          0.0  1829.524191  -600.000000\n",
              "2         C  [10, 40]      4.0  ...          0.0  1693.723877   350.000000\n",
              "3         D  [40, 10]      3.0  ...          0.0  3000.000000 -4645.809559\n",
              "4         E  [25, 25]      2.5  ...          0.0  2776.645130  -150.000000\n",
              "\n",
              "[5 rows x 9 columns]"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 6
        }
      ]
    }
  ]
}