{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "83f86a7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import matplotlib.pyplot as plt\n",
    "from collections import OrderedDict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "38e97e1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Create PINN, following from tensorflow example\n",
    "\n",
    "class PINN(nn.Module):\n",
    "    def __init__(\n",
    "            self,\n",
    "            N_input,\n",
    "            N_hidden_arr,\n",
    "            N_output,\n",
    "            activation = nn.Tanh\n",
    "    ):\n",
    "        super(PINN, self).__init__()  # Create PINN object\n",
    "\n",
    "        # Create input layer w/ activation function\n",
    "        layers = [('Input', nn.Linear(N_input, N_hidden_arr[0]))]\n",
    "        layers.append(('Input activation', activation()))\n",
    "\n",
    "        # Create hidden layers\n",
    "        for i in range(len(N_hidden_arr)-1):\n",
    "            layers.append(\n",
    "                (\"Hidden %d\" % (i+1), nn.Linear(N_hidden_arr[i], N_hidden_arr[i+1]))\n",
    "            )\n",
    "            layers.append(('Hidden activation %d' % (i+1), activation()))\n",
    "        layers.append(('Output', nn.Linear(N_hidden_arr[-1], N_output)))\n",
    "        layerdict = OrderedDict(layers)\n",
    "        self.layers = nn.Sequential(layerdict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f0d875b6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sequential(\n",
      "  (Input): Linear(in_features=2, out_features=32, bias=True)\n",
      "  (Input activation): Tanh()\n",
      "  (Hidden 1): Linear(in_features=32, out_features=16, bias=True)\n",
      "  (Hidden activation 1): Tanh()\n",
      "  (Hidden 2): Linear(in_features=16, out_features=16, bias=True)\n",
      "  (Hidden activation 2): Tanh()\n",
      "  (Hidden 3): Linear(in_features=16, out_features=32, bias=True)\n",
      "  (Hidden activation 3): Tanh()\n",
      "  (Output): Linear(in_features=32, out_features=2, bias=True)\n",
      ")\n"
     ]
    }
   ],
   "source": [
    "test = PINN(2, [32,16,16,32], 2)\n",
    "print(test.layers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "b7d0f531",
   "metadata": {},
   "outputs": [],
   "source": [
    "## constant values\n",
    "u0 = 1  # defining velocity, density, and viscocity for fluid\n",
    "rho0 = 1\n",
    "nu0 = 0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "cb620c0e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([84, 2])\n",
      "tensor([[1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.],\n",
      "        [1., 0.]])\n",
      "torch.Size([84, 2])\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAATkxJREFUeJztnX9sVNeZ9x88OMQexk6j1HZspkaTkBYlW7wlgtA2C0va4Je0af5YlRrk0LLbtApoq6CWQn8Esv1Btq260XZpFuJSqghK2qg4q8avk25id9VCTF9Ss2mTzcK6eAZc7CAtY4QTZmKfV88Z23iGcZjxvc+598z9fiT3cM69/vLke6b3eZh77zmzlFKKAAAAAAA8osyrvxgAAAAAgEExAgAAAABPQTECAAAAAE9BMQIAAAAAT0ExAgAAAABPQTECAAAAAE9BMQIAAAAAT0ExAgAAAABPmU0WMDY2RgMDAxSJRGjWrFlehwMAAACAAuB1VS9cuED19fVUVlZmdzHChUg0GvU6DAAAAADMgEQiQfPmzbO7GOFvRCb+Y6qqqlzTTafT9Pzzz9Pdd99N5eXlrumCbOCzOeC1GeCzGeCz/T4PDw/rLxMm8rjVxcjErRkuRNwuRiorK7UmPuhywGdzwGszwGczwOfS8flqj1jgAVYAAAAAeAqKEQAAAAB4CooRAAAAAHgKihEAAAAAeAqKEQAAAAB4CooRAAAAAHgKihEAAAAAeAqKEQAAAAB4SuCKkXg8Tr/73e/yHuNxPj5TkskknT59Ou8xHufjftKV1LYxZkltxGxGGzGb0bYxZkltG2OOC+bCGaGK5Ne//rX62Mc+pm688UbFv37o0KGr/k5XV5f6y7/8S3XNNdeom266Sf34xz8u6u9MJpP67+LWCf39/Wru3Llq9uzZ6qWXXlKpVEq1t7frlvs8zsf5vGI5f/68uuOOO1QsFlPxeDzrGPd5nI/zeX7QNRnzVJ/9GrOktsmY3fIaPr+zNny2y2fJuG2cw37BXDjT/F10MdLR0aG++tWvql/84hcFFSN9fX2qsrJSbd68Wb366qvqBz/4gQqFQqqzs9N4MXL06FFtMmtxe/iZZ/QEcDt1nM8rlkQioT8YrKE/OD09Sr34om6njvN5ftA1GfOfjhzRPnPr15gltU3G7JbX8PmdteGzXT5Lxm3jHB4VzIXGipGsXy6gGNmyZYu69dZbs8bWrFmjVq1aZbwYYSaqPtaLVFToCeB2wnw+PlMmKlX9ASFSvx1vJz9IOZWt17qmYl447jO3fo5ZUttUzG56DZ+n14bP9vksGbeNc/iSYC70ZTFy5513qi984QtZY3v37lVVVVXT/s5bb72lA5/44aqP/65z587pr5Gc/nD1x6Zff/31egK45T6PO9XmCp7/j1Mx5Yf7PO5HXRMxT/XZ7zHbPoduew2f4XMp+Wy7HxUua0vmwokfztuFFCOzxouKGe/Cd+jQIbrvvvumPeeWW26hz3zmM7Rt27bJsY6ODrrnnntoZGSEKioqrvidHTt20COPPHLF+IEDB/TOggAAAADwP5zn165dqx+05V2Bp0UJfzOyYMEC9e1vfztr7Nlnn9W/OzIykvd38M1IaVTdNv7rxvY5xL/Y4TN8Ll0/Kkr4mxFf3qbJBc+M2Hc/0sb7vpLaNt5jh8/Ta8Nn+3yWjNvGOXwpaM+M8AOst912W9ZYS0uLJw+wGn+bpqtL5kltF3RNxiz+RLxlfkjGLPqWB3yGz5b6LBm3jXN4tBTeprlw4YL6/e9/r3/4L/j+97+v/zzxPvLWrVtVa2vrFa/2fulLX1Kvvfaa2rVrl2ev9mKdEW9inuqzX2OW1LZxXQb4/M7a8NkunyXjtnEO+0thnRFewIyFc3/Wr1+vj3O7fPnyK36nqalJL3rG5nm16BnD5k5Ue1MngOFxJ+bzB2K6CpXHZ/JhlNSV1J6qm+uzX2OW1DYVs5tew+fpteGzfT7natvoh5vakrlwJvnb0ds0phgeHqbq6uqrP41bJOl0Wr/Zs3r1aiovL3dNF2QDn80Br80An80An+33udD8Hbi9aQAAAADgL1CMAAAAAMBTUIwAAAAAwFNQjAAAAADAU1CMAAAAAMBTAleM8BO9p0+fznuMx/k4AAAAUMokfZYLA1WMsLnNzc20fPlySiQSWce4z+N8HAUJAACAUiXpw1wYqGLkwoULNDQ0RH19fbRixQo6c+yYHueW+zzOx/k8x3DF2dWVad1ESldSe2Aguw26H5IxS3kNn7OBz/K6tl47LJjDCyZzYaEoC3BzBVbJTZgmaWtTqqyMV5PLtNx3AyldSe22NpUKhzOr+4XDdsQsqS0cs4jX8PkKbfgsrGvrtcOiOYybyIWmNsqzsRhh9CZDRHobZp4AbrmvNyFyCi/bO/GBmfgJhTLjftQ1EHNq3GdufR+zpLaBmF33Gj7n1YbPlvosGbeFcxiXzIVF5u9A3aaZIHrxIj2ZM8b96MiIc/ETJ4jGxrLHRkeJTp70p66kto0xS2ojZjPaiNmMto0xS2pbGHNUMhcWSSCLkUQ4TK05Y9xPVFY6F1+wgKgsx9ZQiOjmm/2pK6ltY8yS2ojZjDZiNqNtY8yS2hbGnJDMhUUSuGKEnxRe0dJCfUQ0f3yMW+7zeO6TxUUzbx7Rnj2ZDwrD7e7dmXE/6iJmc9qIGTEj5tLVtizmhHQuLJJA7drL707zK0v8pHAsFqMX9u+n44ODtKi2lu5at25y/Ne//jXNc/oB4qed+Ss0rlzd+DBK6wpqp/v7qaO3l1Y3NVF5YyO5ioV+SMYs5jV8zgI+G9C19dphwRyeNpgLC83fsylARCIRqqmp0X/u7u6muro6Ot7RQQ2LF+s+v9LEx/k8x/AEuv1BlNSV1K6vJ+rtzbRuY6MfkjFLeQ2fs4HP8rq2XjssmMOIyVxYIIEqRrg66+zs1O9Oc7WXTqcnj0WjUV0Fsvl8HgAAAFCKVPswFwaqGGHY3OkMdnxrBgAAALCAap/lwsA9wAoAAAAAf4FiBAAAAACegmIEAAAAAJ6CYgQAAAAAnoJiBAAAAACeErhihBde4QVf8sHjfBwAAAAoZZI+y4WBKkbY3ObmZr3yXO5St9zncT6OggQAAECpkvRhLgxUMcILvAwNDemlbnmFuTPHjulxbrnP43ycz3MMV5xdXZnWTaR0JbUHBrLboPshGbOU1/A5G/gsr2vrtcOCObxgMhcWirKAZDLJ++fo1inxeFzFYjGtt7CiQrW3t+uW+zzOxx3T1qZUWRlv+pNpue8GUrqS2m1tKhUOa5+5tSJmSW3hmEW8hs9XaMNnYV1brx0WzWHcRC4sIn8Hrhhh4j09KkakKsYngFvu87hjEonLH5iJn1AoM+5HXQMxp8Z95tb3MUtqG4jZda/hc15t+Gypz5JxWziHcclcWGT+DtRtmgmiFy/Skzlj3I+OjDgXP3GCaGwse2x0NLPToh91JbVtjFlSGzGb0UbMZrRtjFlS28KYo5K5sEgCWYwkwmFqzRnjfqKy0rn4ggVEZTm2hkKZLZ/9qCupbWPMktqI2Yw2YjajbWPMktoWxpyQzIVFErhihJ8UXtHSQn1ENH98jFvu83juk8VFwxsM7dmT+aAw3O7e7XzbZyldxGxOGzEjZsRcutqWxZyQzoVFMovv1ZDPGR4e1rsL8mtGVVVVM9bhd6f5lSV+UjgWi9EL+/fT8cFBWlRbS3etWzc5ztsnO961kJ925q/QuHJ1cwdEKV1B7XR/P3X09tLqpiYqb2wkV7HQD8mYxbyGz1nAZwO6tl47LJjD0wZzYaH5ezYFiEgkQjU1NfrP3d3dVFdXR8c7Oqhh8WLd51ea+Dif5xieQIltmKV0JbXr64l6ezOt29joh2TMUl7D52zgs7yurdcOC+YwYjIXFkigihGuzjo7O/W701ztpdPpyWPRaFRXgWw+nwcAAACUItU+zIWBKkYYNnc6gx3fmgEAAAAsoNpnuTBwD7ACAAAAwF+gGAEAAACAp6AYAQAAAICnoBgBAAAAgKegGAEAAACApwSuGOGFV3jBl3zwOB8HAAAASpmkz3JhoIoRNre5uVmvPJe71C33eZyPoyABAABQqiR9mAsDVYzwAi9DQ0N6qVteYe7MsWN6nFvu8zgf5/McwxVnV1emdRMpXUntgYHsNuh+SMYs5TV8zgY+y+vaeu2wYA4vmMyFhaIsIJlM8v45unVKPB5XsVhM6y2sqFDt7e265T6P83HHtLUpVVbGm/5kWu67gZSupHZbm0qFw9pnbq2IWVJbOGYRr+HzFdrwWVjX1muHRXMYN5ELi8jfgStGmHhPj4oRqYrxCeCW+zzumETi8gdm4icUyoz7UddAzKlxn7n1fcyS2gZidt1r+JxXGz5b6rNk3BbOYVwyFxaZvwN1m2aC6MWL9GTOGPejIyPOxU+cIBobyx4bHc3stOhHXUltG2OW1EbMZrQRsxltG2OW1LYw5qhkLiySQBYjiXCYWnPGuJ+orHQuvmABUVmOraFQZstnP+pKatsYs6Q2YjajjZjNaNsYs6S2hTEnJHNhkQSuGOEnhVe0tFAfEc0fH+OW+zye+2Rx0fAGQ3v2ZD4oDLe7dzvf9llKFzGb00bMiBkxl662ZTEnpHNhkcziezXkc4aHh/XugvyaUVVV1Yx1+N1pfmWJnxSOxWL0wv79dHxwkBbV1tJd69ZNjvP2yY53LeSnnfkrNK5c3dwBUUpXUDvd308dvb20uqmJyhsbyVUs9EMyZjGv4XMW8NmArq3XDgvm8LTBXFho/p5NASISiVBNTY3+c3d3N9XV1dHxjg5qWLxY9/mVJj7O5zmGJ1BiG2YpXUnt+nqi3t5M6zY2+iEZs5TX8Dkb+Cyva+u1w4I5jJjMhQUSqGKEq7POzk797jRXe+l0evJYNBrVVSCbz+cBAAAApUi1D3NhoIoRhs2dzmDHt2YAAAAAC6j2WS4M3AOsAAAAAPAXKEYAAAAA4CkoRgAAAADgKShGAAAAAOApKEYAAAAA4CmBK0Z44RVe8CUfPM7HAQAAgFIm6bNcOKNiZNeuXTR//ny69tpraenSpXT06NF3PP+xxx6j9773vVRRUaHfYX7ooYforbfeItOwuc3NzXrludylbrnP43wcBQkAAIBSJenDXFh0MfLUU0/R5s2bafv27fTyyy/TokWLaNWqVTQ0NJT3/AMHDtDWrVv1+a+99hr96Ec/0hpf+cpXyDS8wAvHyUvd8gpzZ44d0+Pccp/H+Tif5xiuOLu6Mq2bSOlKag8MZLdB90MyZimv4XM28Fle19ZrhwVzeMFkLiwUVSRLlixRGzdunOyPjo6q+vp6tXPnzrzn87krV67MGtu8ebP60Ic+VPDfmUwmef8c3TolHo+rWCym9RZWVKj29nbdcp/H+bhj2tqUKivjTX8yLffdQEpXUrutTaXCYe0zt1bELKktHLOI1/D5Cm34LKxr67XDojmMm8iFReTvojbKS6VSVFlZSU8//TTdd999k+Pr16+n8+fP0zPPPJP3m5EHH3yQnn/+eVqyZImuuO655x5qbW2d9tuRS5cu6Z+pG+3w7Z1z58452ihvAq7+Vq9cSYMVFbR3717asGED1b75JnW8+KJem98RXMHfeivR2NjlMd5h8Q9/cLa/gpSugZjTc+bQr/bupY9u2EDlqZS/Y5bUNhCz617D57za8NlSnyXjtnAOz0jmwin5+4YbbrjqRnlFFSMDAwPU0NBAhw8fpmXLlk2Ob9myRa9l39PTk/f3/vmf/5m++MUv8rcw9Pbbb9PnP/95evzxx6f9e3bs2EGPPPJI3sKGiyEAAAAA+J+RkRFau3at97v28g6A3/72t+mHP/yhftj15MmT9IUvfIG+8Y1v0Ne//vW8v7Nt2zb9XEruNyN33303vhmxqOq28l83kto2/ksSPufVhs+W+iwZt4VzeMbQNyMFUcy9n0uXLqlQKKQOHTqUNX7//fere++9N+/vfPjDH1Zf/OIXs8aefPJJVVFRoZ83KdlnRkKhzL09bt28HymhK6nN933nzs3c9507146YJbWFYxbxGj5foQ2fhXVtvXZYNIdxnz0zMqMHWDdt2jTZ54KioaFh2gdYP/CBD6gtW7ZkjR04cEAXI2+//bbRYiSRSEyaz+2fjhzRE8Dt1HE+zzGs0dWVad1ESldQO3XqVOaCcuqUch0L/ZCMWcxr+JwFfDaga+u1w4I5TBjMhWLFyMGDB9WcOXPUvn371KuvvqoeeOABdd1116mzZ8/q462trWrr1q2T52/fvl1FIhH105/+VPX19annn39e3XTTTeqTn/yk6/8xV+P8+fPqjjvumKz6UqlU5oOeSk1WiXyczwPuMdVnIAu8NgN8NgN8lsFkLiw0fxf9zMiaNWvojTfeoIcffpjOnj1LTU1N1NnZSbW1tfp4PB6nsrLLy5d87Wtfo1mzZun2zJkz9O53v5s+/vGP07e+9S0yTXV1tY6V352eN28epdPpyWP8TAo/hBuJRPR5AAAAQClS7cNcOKMHWDdt2qR/pntgNesvmD1bL3jGP36AzZ3OYJ4UAAAAoNSp9lkuDNzeNAAAAADwFyhGAAAAAOApKEYAAAAA4CkoRgAAAADgKYErRnhJ2tPT7HjI4ya3TAYAAAC8wG+5MFDFCJvb3NxMy5cvp0QikXWM+zzOx1GQAAAAKFWSPsyFgSpG+J3qoaEhvXPwihUr9Lr8DLfc53E+zuc5hivOrq5M6yZSupLavK/C1DbofkjGLOU1fM4GPsvr2nrtsGAOL5jMhYWiLMDKvWnKyjJ7CHDr5r4HErqS2ry/RDicWd0vHLYjZklt4ZhFvIbPV2jDZ2FdW68dFs1h3Pa9aWwvRph4T4+KEen9cXgCuOU+jzuG1/Kf+MBM/PCmRk7X+JfSNRBzatxnbn0fs6S2gZhd9xo+59WGz5b6LBm3hXMYl8yFRebvQN2mmSB68SI9mTPG/ejIiHPxEyeyt3lmRkeJTp70p66kto0xS2ojZjPaiNmMto0xS2pbGHNUMhcWSSCLkUQ4TK05Y9xPVFY6F1+wgGjK3jyaUIjo5pv9qSupbWPMktqI2Yw2YjajbWPMktoWxpyQzIVFErhihJ8UXtHSQn1ENH98jFvu83juk8VFw2v679mT+aAw3O7enRn3oy5iNqeNmBEzYi5dbctiTkjnwiKZxfdqyOcMDw/rDX34NaOqqqoZ6/C70/zKEj8pHIvF6IX9++n44CAtqq2lu9atmxznHQsdbxTETzvzV2hcubq56ZCUrqB2ur+fOnp7aXVTE5U3NpKrWOiHZMxiXsPnLOCzAV1brx0WzOFpg7mw0Pw9o117bYW3RK6pqZncXbiuro6Od3RQw+LFus+vNPFxPs8xPIESOx9K6Upq19cT9fZmWrex0Q/JmKW8hs/ZwGd5XVuvHRbMYcRkLiyQQBUjXJ11dnbqd6e52kun05PHotGorgLZ/Om2VQYAAABsp9qHuTBQxQjD5k5nsONbMwAAAIAFVPssFwbuAVYAAAAA+AsUIwAAAADwFBQjAAAAAPAUFCMAAAAA8BQUIwAAAADwlMAVI7zwCi/4kg8e5+MAAABAKZP0WS4MVDHC5jY3N+uV53KXuuU+j/NxFCQAAABKlaQPc2GgihFe4GVoaEgvdcsrzJ05dkyPc8t9HufjfJ5juOLs6sq0biKlK6k9MJDdBt0PyZilvIbP2cBneV1brx0WzOEFk7mwUJQFJJNJ3j9Ht06Jx+MqFotpvYUVFaq9vV233OdxPu6Ytjalysp4059My303kNKV1G5rU6lwWPvMrRUxS2oLxyziNXy+Qhs+C+vaeu2waA7jJnJhEfk7cMUIE+/pUTEiVTE+Adxyn8cdk0hc/sBM/IRCmXE/6hqIOTXuM7e+j1lS20DMrnsNn/Nqw2dLfZaM28I5jEvmwiLzd6Bu00wQvXiRnswZ4350ZMS5+IkTRGNj2WOjo5mdFv2oK6ltY8yS2ojZjDZiNqNtY8yS2hbGHJXMhUUSyGIkEQ5Ta84Y9xOVlc7FFywgKsuxNRTKbPnsR11JbRtjltRGzGa0EbMZbRtjltS2MOaEZC4sksAVI/yk8IqWFuojovnjY9xyn8dznywuGt5gaM+ezAeF4Xb3bufbPkvpImZz2ogZMSPm0tW2LOaEdC4skll8r4Z8zvDwsN5dkF8zqqqqmrEOvzvNryzxk8KxWIxe2L+fjg8O0qLaWrpr3brJcd4+2fGuhfy0M3+FxpWrmzsgSukKaqf7+6mjt5dWNzVReWMjuYqFfkjGLOY1fM4CPhvQtfXaYcEcnjaYCwvN37MpQEQiEaqpqdF/7u7uprq6Ojre0UENixfrPr/SxMf5PMfwBEpswyylK6ldX0/U25tp3cZGPyRjlvIaPmcDn+V1bb12WDCHEZO5sEACVYxwddbZ2anfneZqL51OTx6LRqO6CmTz+TwAAACgFKn2YS4MVDHCsLnTGez41gwAAABgAdU+y4WBe4AVAAAAAP4CxQgAAAAAPAXFCAAAAAA8BcUIAAAAADwFxQgAAAAAPCVwxQgvvMILvuSDx/k4AAAAUMokfZYLA1WMsLnNzc165bncpW65z+N8HAUJAACAUiXpw1wYqGKEF3gZGhrSS93yCnNnjh3T49xyn8f5OJ/nGK44u7oyrZtI6UpqDwxkt0H3QzJmKa/hczbwWV7X1muHBXN4wWQuLBRlAclkkvfP0a1T4vG4isViWm9hRYVqb2/XLfd5nI87pq1NqbIy3vQn03LfDaR0JbXb2lQqHNY+c2tFzJLawjGLeA2fr9CGz8K6tl47LJrDuIlcWET+DlwxwsR7elSMSFWMTwC33OdxxyQSlz8wEz+hUGbcj7oGYk6N+8yt72OW1DYQs+tew+e82vDZUp8l47ZwDuOSubDI/B2o2zQTRC9epCdzxrgfHRlxLn7iBNHYWPbY6Ghmp0U/6kpq2xizpDZiNqONmM1o2xizpLaFMUclc2GRBLIYSYTD1Jozxv1EZaVz8QULiMpybA2FMls++1FXUtvGmCW1EbMZbcRsRtvGmCW1LYw5IZkLiyRwxQg/KbyipYX6iGj++Bi33Ofx3CeLi4Y3GNqzJ/NBYbjdvdv5ts9SuojZnDZiRsyIuXS1LYs5IZ0Li2QW36shnzM8PKx3F+TXjKqqqmasw+9O8ytL/KRwLBajF/bvp+ODg7SotpbuWrducpy3T3a8ayE/7cxfoXHl6uYOiFK6gtrp/n7q6O2l1U1NVN7YSK5ioR+SMYt5DZ+zgM8GdG29dlgwh6cN5sJC8/dsChCRSIRqamr0n7u7u6muro6Od3RQw+LFus+vNPFxPs8xPIES2zBL6Upq19cT9fZmWrex0Q/JmKW8hs/ZwGd5XVuvHRbMYcRkLiyQQBUjXJ11dnbqd6e52kun05PHotGorgLZfD4PAAAAKEWqfZgLA1WMMGzudAY7vjUDAAAAWEC1z3Jh4B5gBQAAAIC/QDECAAAAAE9BMQIAAAAAT0ExAgAAAABPQTECAAAAAE8JXDHCC6/wgi/54HE+DgAAAJQySZ/lwhkVI7t27aL58+fTtddeS0uXLqWjR4++4/nnz5+njRs30o033khz5syhW265hTo6Osg0bG5zc7NeeS53qVvu8zgfR0ECAACgVEn6MBcWXYw89dRTtHnzZtq+fTu9/PLLtGjRIlq1ahUNDQ3lPT+VStFHP/pROnXqFD399NP0+uuv0xNPPEENDQ1kGl7ghePkpW55hbkzx47pcW65z+N8nM9zDFecXV2Z1k2kdCW1Bway26D7IRmzlNfwORv4LK9r67XDgjm8YDIXFooqkiVLlqiNGzdO9kdHR1V9fb3auXNn3vMff/xxFYvFVCqVUjMlmUzy/jm6dUo8HtfxsN7CigrV3t6uW+7zOB93TFubUmVlvOlPpuW+G0jpSmq3talUOKx95taKmCW1hWMW8Ro+X6ENn4V1bb12WDSHcRO5sIj8XdRGefwtR2Vlpf6G47777pscX79+vb4V88wzz1zxO6tXr6brr79e/x4ff/e7301r166lL3/5yxSa2IEwh0uXLumfqRvt8BK1586dc7RR3gRc/a1euZIGKypo7969tGHDBqp9803qePFFvTa/I7iCv/VWorGxy2P83/mHPzjbX0FK10DM6Tlz6Fd799JHN2yg8lTK3zFLahuI2XWv4XNebfhsqc+ScVs4h2ckc+GU/H3DDTdcdaO8ooqRgYEBfXvl8OHDtGzZssnxLVu26LXse3p6rvid973vffoWzbp16+jBBx+kkydP6vbv//7v9a2efOzYsYMeeeSRK8YPHDigixoAAAAA+J+RkRH9BYTnu/aOjY3p3f/27NmjvwlZvHgxnTlzhr773e9OW4xs27ZNP5eS+83I3XffjW9GLKq6rfzXjaS2jf+ShM95teGzpT5Lxm3hHJ4x9M1IQRRz7+fSpUsqFAqpQ4cOZY3ff//96t577837O3/1V3+l7rrrrqyxjo4OfQ+J9Ur2mZFQKHNvj1s370dK6Epq833fuXMz933nzrUjZklt4ZhFvIbPV2jDZ2FdW68dFs1h3GfPjMzoAdZNmzZlPcDa0NAw7QOs27ZtU42Njfq8CR577DF14403Fvx3ulWMJBKJSfO5/dORI3oCuJ06zuc5hjW6ujKtm0jpCmqnTp3KXFBOnVKuY6EfkjGLeQ2fs4DPBnRtvXZYMIcJg7lQrBg5ePCgmjNnjtq3b5969dVX1QMPPKCuu+46dfbsWX28tbVVbd26dfJ8rq4ikYguYF5//XX1y1/+UtXU1KhvfvObrv/HXI3z58+rO+64Y7Lq4zd89Ac9lZqsEvk4nwfcY6rPQBZ4bQb4bAb4LIPJXFho/i76mZE1a9bQG2+8QQ8//DCdPXuWmpqaqLOzk2pra/XxeDxOZWWXly/hZz2ee+45euihh+j973+/fgD2C1/4gn6bxjTV1dU6Vn53et68eZROp7Pi5IdwI5GIPg8AAAAoRap9mAtn9ADrpk2b9E8+uru7rxjjN29eeukl8gNs7nQG86QAAAAApU61z3Jh4PamAQAAAIC/QDECAAAAAE9BMQIAAAAAT0ExAgAAAABPQTECAAAAAE8JXDHC6+Ofnmb7ZR7n4wAAAEApk/RZLgxUMcLmNjc30/LlyymRSGQd4z6P83EUJAAAAEqVpA9zYaCKEV7gZWhoiPr6+mjFihV6kyCGW+7zOB/n8xzDFWdXV6Z1EyldSW3e5GlqG3Q/JGOW8ho+ZwOf5XVtvXZYMIcXTObCQlEWYOVGeWVlmQ2NuHVzEyYJXUlt3uwqHM4sNRwO2xGzpLZwzCJew+crtOGzsK6t1w6L5jBu+0Z5thcjTLynR8WIVMX4BHDLfR53DG8sNPGBmfjhHRadbjgkpWsg5tS4z9z6PmZJbQMxu+41fM6rDZ8t9VkybgvnMC6ZC4vM34G6TTNB9OJFejJnjPvRkRHn4idOEI2NZY+NjhKdPOlPXUltG2OW1EbMZrQRsxltG2OW1LYw5qhkLiySQBYjiXCYWnPGuJ+orHQuvmAB0ZSNAjWhENHNN/tTV1LbxpgltRGzGW3EbEbbxpgltS2MOSGZC4skcMUIPym8oqWF+oho/vgYt9zn8dwni4uGNxjasyfzQWG43b07M+5HXcRsThsxI2bEXLralsWckM6FRTKL79WQzxkeHta7C/JrRlVVVTPW4Xen+ZUlflI4FovRC/v30/HBQVpUW0t3rVs3Oc7bJzvetZCfduav0LhydXMHRCldQe10fz919PbS6qYmKm9sJFex0A/JmMW8hs9ZwGcDurZeOyyYw9MGc2Gh+Xs2BYhIJEI1NTX6z93d3VRXV0fHOzqoYfFi3edXmvg4n+cYnkCJbZildCW16+uJenszrdvY6IdkzFJew+ds4LO8rq3XDgvmMGIyFxZIoIoRrs46Ozv1u9Nc7aXT6clj0WhUV4FsPp8HAAAAlCLVPsyFgSpGGDZ3OoMd35oBAAAALKDaZ7kwcA+wAgAAAMBfoBgBAAAAgKegGAEAAACAp6AYAQAAAICnBK4Y4Xed+R3rfPC4yS2TAQAAAC/wWy4MVDHC5jY3N+vFXnJXl+M+j/NxFCQAAABKlaQPc2GgihF+p3poaEivLseLupw5dkyPc8t9HufjfJ5juOLs6sq0biKlK6k9MJDdBt0PyZilvIbP2cBneV1brx0WzOEFk7mwUJQFFLoFcSHE43EVi8W03sLxbZO55T6P83HHtLVd3u6ZW+67gZSupHZbm0qFw5ltwMNhO2KW1BaOWcRr+HyFNnwW1rX12mHRHMZN5MIi8nfgihEm3tOjYkSqYnwCuOU+jzsmkbj8gZn4CYUy437UNRBzatxnbn0fs6S2gZhd9xo+59WGz5b6LBm3hXMYl8yFRebvQN2mmSB68SI9mTPG/ejIiHPxEyeIxsayx0ZHM5sb+VFXUtvGmCW1EbMZbcRsRtvGmCW1LYw5KpkLiySQxUgiHKbWnDHuJyornYsvWEBUlmMrb/fMuyz6UVdS28aYJbURsxltxGxG28aYJbUtjDkhmQuLJHDFCD8pvKKlhfqIaP74GLfc5/HcJ4uLhtf037Mn80FhuN292/lOi1K6iNmcNmJGzIi5dLUtizkhnQuLZBbfqyGfMzw8rDf04deMqqqqZqzD707zK0v8pHAsFqMX9u+n44ODtKi2lu5at25ynHcsdLxRED/tzF+hceXq5qZDUrqC2un+furo7aXVTU1U3thIrmKhH5Ixi3kNn7OAzwZ0bb12WDCHpw3mwkLzd6B27eUtkWtqavSfu7u7qa6ujo53dFDD4sW6z6808XE+zzE8gRI7H0rpSmrX1xP19mZat7HRD8mYpbyGz9nAZ3ldW68dFsxhxGQuLJBAFSNcnXV2dup3p7naS6fTk8ei0aiuAtn86bZVBgAAAGyn2oe5MFDFCMPmTmew41szAAAAgAVU+ywXBu4BVgAAAAD4CxQjAAAAAPAUFCMAAAAA8BQUIwAAAADwFBQjAAAAAPCUwBUjvPAKL/iSDx7n4wAAAEApk/RZLgxUMcLmNjc365Xncpe65T6P83EUJAAAAEqVpA9zYaCKEV7gZWhoSC91yyvMnTl2TI9zy30e5+N8nmO44uzqyrRuIqUrqT0wkN0G3Q/JmKW8hs/ZwGd5XVuvHRbM4QWTubBQlAUkk0neP0e3TonH4yoWi2m9hRUVqr29Xbfc53E+7pi2NqXKynjTn0zLfTeQ0pXUbmtTqXBY+8ytFTFLagvHLOI1fL5CGz4L69p67bBoDuMmcmER+TtwxQgT7+lRMSJVMT4B3HKfxx2TSFz+wEz8hEKZcT/qGog5Ne4zt76PWVLbQMyuew2f82rDZ0t9lozbwjmMS+bCIvN3oG7TTBC9eJGezBnjfnRkxLn4iRNEY2PZY6OjmZ0W/agrqW1jzJLaiNmMNmI2o21jzJLaFsYclcyFRRLIYiQRDlNrzhj3E5WVzsUXLCAqy7E1FMps+exHXUltG2OW1EbMZrQRsxltG2OW1LYw5oRkLiySwBUj/KTwipYW6iOi+eNj3HKfx3OfLC4a3mBoz57MB4Xhdvdu59s+S+kiZnPaiBkxI+bS1bYs5oR0LiySWXyvhnzO8PCw3l2QXzOqqqqasQ6/O82vLPGTwrFYjF7Yv5+ODw7Sotpaumvduslx3j7Z8a6F/LQzf4XGlaubOyBK6Qpqp/v7qaO3l1Y3NVF5YyO5ioV+SMYs5jV8zgI+G9C19dphwRyeNpgLC83fsylARCIRqqmp0X/u7u6muro6Ot7RQQ2LF+s+v9LEx/k8x/AESmzDLKUrqV1fT9Tbm2ndxkY/JGOW8ho+ZwOf5XVtvXZYMIcRk7mwQAJVjHB11tnZqd+d5movnU5PHotGo7oKZPP5PAAAAKAUqfZhLgxUMcKwudMZ7PjWDAAAAGAB1T7LhYF7gBUAAAAA/gLFCAAAAAA8BcUIAAAAADwFxQgAAAAAPAXFCAAAAAA8JXDFCC+8wgu+5IPH+TgAAABQyiR9lgtnVIzs2rWL5s+fT9deey0tXbqUjh49WtDvHTx4kGbNmkX33XcfeQGb29zcrFeey13qlvs8zsdRkAAAAChVkj7MhUUXI0899RRt3ryZtm/fTi+//DItWrSIVq1aRUNDQ+/4e6dOnaIvfvGLdOedd5JX8AIvHCcvdcsrzJ05dkyPc8t9HufjfJ5juOLs6sq0biKlK6k9MJDdBt0PyZilvIbP2cBneV1brx0WzOEFk7mwUFSRLFmyRG3cuHGyPzo6qurr69XOnTun/Z23335bffCDH1RtbW1q/fr16hOf+ERRf2cymeT9c3TrlHg8rmKxmNZbWFGh2tvbdct9HufjjmlrU6qsjDf9ybTcdwMpXUnttjaVCoe1z9xaEbOktnDMIl7D5yu04bOwrq3XDovmMG4iFxaRv4vaKC+VSlFlZSU9/fTTWbda1q9fT+fPn6dnnnkm7+/xtyj/+Z//SYcOHaJPf/rT+tz29vZp/55Lly7pn6kb7fAStefOnXO0Ud4EXP2tXrmSBisqaO/evbRhwwaqffNN6njxRb02vyO4gr/1VqKxsctjvMPiH/7gbH8FKV0DMafnzKFf7d1LH92wgcpTKX/HLKltIGbXvYbPebXhs6U+S8Zt4RyekcyFU/L3DTfccNWN8ooqRgYGBqihoYEOHz5My5YtmxzfsmWLXsu+p6fnit/5zW9+Q5/61Keot7dXB1RIMbJjxw565JFHrhg/cOCALoYAAAAA4H9GRkZo7dq13u7ay/ebWltb6YknntCFSKFs27ZNP5eS+83I3XffjW9GLKq6rfzXjaS2jf+ShM95teGzpT5Lxm3hHJ4x9M1IQRRz7+fSpUsqFAqpQ4cOZY3ff//96t57773i/N///vf6XhH/zsTPrFmz9A//+eTJk67ec/LVMyOhUObeHrdu3o+U0JXU5vu+c+dm7vvOnWtHzJLawjGLeA2fr9CGz8K6tl47LJrDuM+eGZnRA6ybNm3KeoC1oaEh7wOsb775pnrllVeyfvjh1ZUrV+o/c3Hj5n/M1UgkEpPmc/unI0f0BHA7dZzPcwxrdHVlWjeR0hXUTp06lbmgnDqlXMdCPyRjFvMaPmcBnw3o2nrtsGAOEwZzYaH5u+jbNHz7hB9Yvf3222nJkiX02GOP0cWLF+kzn/mMPn7//ffr50p27typ1yG57bbbsn7/uuuu023uuAkikQjV1NToP3d3d1NdXR0d7+jQX0dxn19p4uN8nmN4C2aJbZildCW1+WvE3l7nX1WWih+SMUt5DZ+zgc/yurZeOyyYw4jJXFggRRcja9asoTfeeIMefvhhOnv2LDU1NVFnZyfV1tbq4/F4nMrK/Lmwa3V1tY6Vn2WZN28epdPpyWP8TAo/hMvm83kAAABAKVLtw1w4owdYN23apH/ywVXVO7Fv3z7yEjZ3OoN5UgAAAIBSp9pnudCfX2EAAAAAIDCgGAEAAACAp6AYAQAAAICnoBgBAAAAgKegGAEAAACApwSuGOH18U9Ps/0yj/NxAAAAoJRJ+iwXBqoYYXObm5tp+fLllEgkso5xn8f5OAoSAAAApUrSh7kwUMUIL/AyNDREfX19eoU53iSI4Zb7PM7H+TzHcMXZ1ZVp3URKV1KbN3ma2gbdD8mYpbyGz9nAZ3ldW68dFszhBZO5sFCUBVi5UV5ZWWZDI27d3IRJQldSmze7Cocz+0uEw3bELKktHLOI1/D5Cm34LKxr67XDojmM275Rnu3FCBPv6VExIlUxPgHccp/HHcMbC018YCZ+eIdFpxsOSekaiDk17jO3vo9ZUttAzK57DZ/zasNnS32WjNvCOYxL5sIi83egbtNMEL14kZ7MGeN+dGTEufiJE0RjY9ljo6NEJ0/6U1dS28aYJbURsxltxGxG28aYJbUtjDkqmQuLJJDFSCIcptacMe4nKiudiy9YQJS7UWAoRHTzzf7UldS2MWZJbcRsRhsxm9G2MWZJbQtjTkjmwiIJXDHCTwqvaGmhPiKaPz7GLfd5PPfJ4qLhDYb27Ml8UBhud+92vu2zlC5iNqeNmBEzYi5dbctiTkjnwiKZxfdqyOcMDw/r3QX5NaOqqqoZ6/C70/zKEj8pHIvF6IX9++n44CAtqq2lu9atmxzn7ZMd71rITzvzV2hcubq5A6KUrqB2ur+fOnp7aXVTE5U3NpKrWOiHZMxiXsPnLOCzAV1brx0WzOFpg7mw0Pw9mwJEJBKhmpoa/efu7m6qq6uj4x0d1LB4se7zK018nM9zDE+gxDbMUrqS2vX1RL29mdZtbPRDMmYpr+FzNvBZXtfWa4cFcxgxmQsLJFDFCFdnnZ2d+t1prvbS6fTksWg0qqtANp/PAwAAAEqRah/mwkAVIwybO53Bjm/NAAAAABZQ7bNcGLgHWAEAAADgL1CMAAAAAMBTUIwAAAAAwFNQjAAAAADAU1CMAAAAAMBTAleM8MIrvOBLPnicjwMAAAClTNJnuTBQxQib29zcrFeey13qlvs8zsdRkAAAAChVkj7MhYEqRniBl6GhIb3ULa8wd+bYMT3OLfd5nI/zeY7hirOrK9O6iZSupPbAQHYbdD8kY5byGj5nA5/ldW29dlgwhxdM5sJCURaQTCZ5/xzdOiUej6tYLKb1FlZUqPb2dt1yn8f5uGPa2pQqK+NNfzIt991ASldSu61NpcJh7TO3VsQsqS0cs4jX8PkKbfgsrGvrtcOiOYybyIVF5O/AFSNMvKdHxYhUxfgEcMt9HndMInH5AzPxEwplxv2oayDm1LjP3Po+ZkltAzG77jV8zqsNny31WTJuC+cwLpkLi8zfgbpNM0H04kV6MmeM+9GREefiJ04QjY1lj42OZnZa9KOupLaNMUtqI2Yz2ojZjLaNMUtqWxhzVDIXFkkgi5FEOEytOWPcT1RWOhdfsICoLMfWUCiz5bMfdSW1bYxZUhsxm9FGzGa0bYxZUtvCmBOSubBIAleM8JPCK1paqI+I5o+Pcct9Hs99srhoeIOhPXsyHxSG2927nW/7LKWLmM1pI2bEjJhLV9uymBPSubBIZvG9GvI5w8PDendBfs2oqqpqxjr87jS/ssRPCsdiMXph/346PjhIi2pr6a516ybHeftkx7sW8tPO/BUaV65u7oAopSuone7vp47eXlrd1ETljY3kKhb6IRmzmNfwOQv4bEDX1muHBXN42mAuLDR/z6YAEYlEqKamRv+5u7ub6urq6HhHBzUsXqz7/EoTH+fzHMMTKLENs5SupHZ9PVFvb6Z1Gxv9kIxZymv4nA18lte19dphwRxGTObCAglUMcLVWWdnp353mqu9dDo9eSwajeoqkM3n8wAAAIBSpNqHuTBQxQjD5k5nsONbMwAAAIAFVPssFwbuAVYAAAAA+AsUIwAAAADwFBQjAAAAAPAUFCMAAAAA8JTAFSP8rjO/Y50PHje5ZTIAAADgBX7LhYEqRtjc5uZmvdhL7upy3OdxPo6CBAAAQKmS9GEuDFQxwu9UDw0N6dXleFGXM8eO6XFuuc/jfJzPcwxXnF1dmdZNpHQltQcGstug+yEZs5TX8Dkb+Cyva+u1w4I5vGAyFxaKsoBCtyAuhHg8rmKxmNZbOL5tMrfc53E+7pi2tsvbPXPLfTeQ0pXUbmtTqXA4sw14OGxHzJLawjGLeA2fr9CGz8K6tl47LJrDuIlcWET+DlwxwsR7elSMSFWMTwC33OdxxyQSlz8wEz+hUGbcj7oGYk6N+8yt72OW1DYQs+tew+e82vDZUp8l47ZwDuOSubDI/B2o2zQTRC9epCdzxrgfHRlxLn7iBNHYWPbY6GhmcyM/6kpq2xizpDZiNqONmM1o2xizpLaFMUclc2GRBLIYSYTD1Jozxv1EZaVz8QULiMpybOXtnnmXRT/qSmrbGLOkNmI2o42YzWjbGLOktoUxJyRzYZEErhjhJ4VXtLRQHxHNHx/jlvs8nvtkcdHwmv579mQ+KAy3u3c732lRShcxm9NGzIgZMZeutmUxJ6RzYZHM4ns15HOGh4f1hj78mlFVVdWMdfjdaX5liZ8UjsVi9ML+/XR8cJAW1dbSXevWTY7zjoWONwrip535KzSuXN3cdEhKV1A73d9PHb29tLqpicobG8lVLPRDMmYxr+FzFvDZgK6t1w4L5vC0wVxYaP4O1K69vCVyTU2N/nN3dzfV1dXR8Y4Oali8WPf5lSY+zuc5hidQYudDKV1J7fp6ot7eTOs2NvohGbOU1/A5G/gsr2vrtcOCOYyYzIUFEqhihKuzzs5O/e40V3vpdHryWDQa1VUgmz/dtsoAAACA7VT7MBcGqhhh2NzpDHZ8awYAAACwgGqf5cLAPcAKAAAAAH+BYgQAAAAAnoJiBAAAAACegmIEAAAAAJ6CYgQAAAAAnhK4YoQXXuEFX/LB43wcAAAAKGWSPsuFMypGdu3aRfPnz6drr72Wli5dSkePHp323CeeeILuvPNOete73qV/PvKRj7zj+ZKwuc3NzXrludylbrnP43wcBQkAAIBSJenDXFh0MfLUU0/R5s2bafv27fTyyy/TokWLaNWqVTQ0NJT3fF7NraWlhbq6uujIkSN6QZW7776bzpw5Q6bhBV44Tl7qlleYO3PsmB7nlvs8zsf5PMdwxdnVlWndREpXUntgILsNuh+SMUt5DZ+zgc/yurZeOyyYwwsmc2GhqCJZsmSJ2rhx42R/dHRU1dfXq507dxb0+2+//baKRCLqJz/5ScF/ZzKZ5P1zdOuUeDyuYrGY1ltYUaHa29t1y30e5+OOaWtTqqyMN/3JtNx3AyldSe22NpUKh7XP3FoRs6S2cMwiXsPnK7Ths7CurdcOi+YwbiIXFpG/i9ooL5VKUWVlJT399NN03333TY6vX7+ezp8/T88888xVNbjS4jXvf/7zn9PHPvaxvOdcunRJ/0zdaIe/UTl37pyjjfIm4Opv9cqVNFhRQXv37qUNGzZQ7ZtvUseLL+q1+R3BFfyttxKNjV0e4x0W//AHZ/srSOkaiDk9Zw79au9e+uiGDVSeSvk7ZkltAzG77jV8zqsNny31WTJuC+fwjGQunJK/b7jhhqtulFdUMTIwMEANDQ10+PBhWrZs2eT4li1b9Fr2PT09V9V48MEH6bnnnqM//vGP+pmTfOzYsYMeeeSRK8YPHDigiyEAAAAA+J+RkRFau3atv3btffTRR+ngwYP6OZLpChFm27Zt+rmU3G9G+FkTfDNiT9Vt5b9uJLVt/JckfM6rDZ8t9Vkybgvn8Iyhb0YKoph7P5cuXVKhUEgdOnQoa/z+++9X99577zv+7ne/+11VXV2tfve736lisfKZkVAoc2+PWzfvR0roSmrzfd+5czP3fefOtSNmSW3hmEW8hs9XaMNnYV1brx0WzWHcZ8+MzOgB1k2bNmU9wNrQ0PCOD7D+4z/+o6qqqlJHjhxRM8GtYiSRSEyaz+2fjhzRE8Dt1HE+zzGs0dWVad1ESldQO3XqVOaCcuqUch0L/ZCMWcxr+JwFfDaga+u1w4I5TBjMhWLFyMGDB9WcOXPUvn371KuvvqoeeOABdd1116mzZ8/q462trWrr1q2T5z/66KPqmmuuUU8//bT685//PPlz4cIF1/9jrsb58+fVHXfcMVn1pVKpzAc9lZqsEvk4nwfcY6rPQBZ4bQb4bAb4LIPJXFho/i76mZE1a9bQG2+8QQ8//DCdPXuWmpqaqLOzk2pra/XxeDxOZWWXly95/PHH9Vs4f/M3f5Olw+uU8IOqJqmurtax8hs98+bNo3Q6PXmMn0nhh3AjkYg+DwAAAChFqn2YC2f0AOumTZv0Tz744dSpnDp1ivwEmzudwTwpAAAAQKlT7bNcGLi9aQAAAADgL1CMAAAAAMBTUIwAAAAAwFNQjAAAAADAU1CMAAAAAMBTAleM8Pr4p6fZfpnH+TgAAABQyiR9lgsDVYywuc3NzbR8+XJKJBJZx7jP43wcBQkAAIBSJenDXBioYoQXeBkaGqK+vj5asWKF3iSI4Zb7PM7H+TzHcMXZ1ZVp3URKV1KbN3ma2gbdD8mYpbyGz9nAZ3ldW68dFszhBZO5sFCUBVi5UV5ZWWZDI27d3IRJQldSmze7CoczSw2Hw3bELKktHLOI1/D5Cm34LKxr67XDojmM275Rnu3FCBPv6VExIlUxPgHccp/HHcMbC018YCZ+eIdFpxsOSekaiDk17jO3vo9ZUttAzK57DZ/zasNnS32WjNvCOYxL5sIi83egbtNMEL14kZ7MGeN+dGTEufiJE0RjY9ljo6NEJ0/6U1dS28aYJbURsxltxGxG28aYJbUtjDkqmQuLJJDFSCIcptacMe4nKiudiy9YQDRlo0BNKER0883+1JXUtjFmSW3EbEYbMZvRtjFmSW0LY05I5sIiCVwxwk8Kr2hpoT4imj8+xi33eTz3yeKi4Q2G9uzJfFAYbnfvzoz7URcxm9NGzIgZMZeutmUxJ6RzYZHM4ns15HOGh4f17oL8mlFVVdWMdfjdaX5liZ8UjsVi9ML+/XR8cJAW1dbSXevWTY7z9smOdy3kp535KzSuXN3cAVFKV1A73d9PHb29tLqpicobG8lVLPRDMmYxr+FzFvDZgK6t1w4L5vC0wVxYaP6eTQEiEolQTU2N/nN3dzfV1dXR8Y4Oali8WPf5lSY+zuc5hidQYhtmKV1J7fp6ot7eTOs2NvohGbOU1/A5G/gsr2vrtcOCOYyYzIUFEqhihKuzzs5O/e40V3vpdHryWDQa1VUgm8/nAQAAAKVItQ9zYaCKEYbNnc5gx7dmAAAAAAuo9lkuDNwDrAAAAADwFyhGAAAAAOApKEYAAAAA4CkoRgAAAADgKShGAAAAAOApgStGeOEVXvAlHzzOxwEAAIBSJumzXBioYoTNbW5u1ivP5S51y30e5+MoSAAAAJQqSR/mwkAVI7zAy9DQkF7qlleYO3PsmB7nlvs8zsf5PMdwxdnVlWndREpXUntgILsNuh+SMUt5DZ+zgc/yurZeOyyYwwsmc2GhKAtIJpO8f45unRKPx1UsFtN6CysqVHt7u265z+N83DFtbUqVlfGmP5mW+24gpSup3damUuGw9plbK2KW1BaOWcRr+HyFNnwW1rX12mHRHMZN5MIi8nfgihEm3tOjYkSqYnwCuOU+jzsmkbj8gZn4CYUy437UNRBzatxnbn0fs6S2gZhd9xo+59WGz5b6LBm3hXMYl8yFRebvQN2mmSB68SI9mTPG/ejIiHPxEyeIxsayx0ZHMzst+lFXUtvGmCW1EbMZbcRsRtvGmCW1LYw5KpkLiySQxUgiHKbWnDHuJyornYsvWEBUlmNrKJTZ8tmPupLaNsYsqY2YzWgjZjPaNsYsqW1hzAnJXFgkgStG+EnhFS0t1EdE88fHuOU+j+c+WVw0vMHQnj2ZDwrD7e7dzrd9ltJFzOa0ETNiRsylq21ZzAnpXFgks/heDfmc4eFhvbsgv2ZUVVU1Yx1+d5pfWeInhWOxGL2wfz8dHxykRbW1dNe6dZPjvH2y410L+Wln/gqNK1c3d0CU0hXUTvf3U0dvL61uaqLyxkZyFQv9kIxZzGv4nAV8NqBr67XDgjk8bTAXFpq/Z1OAiEQiVFNTo//c3d1NdXV1dLyjgxoWL9Z9fqWJj/N5juEJlNiGWUpXUru+nqi3N9O6jY1+SMYs5TV8zgY+y+vaeu2wYA4jJnNhgQSqGOHqrLOzU787zdVeOp2ePBaNRnUVyObzeQAAAEApUu3DXBioYoRhc6cz2PGtGQAAAMACqn2WCwP3ACsAAAAA/AWKEQAAAAB4CooRAAAAAHgKihEAAAAAeAqKEQAAAAB4SuCKEV54hRd8yQeP83EAAACglEn6LBcGqhhhc5ubm/XKc7lL3XKfx/k4ChIAAAClStKHuTBQxQgv8DI0NKSXuuUV5s4cO6bHueU+j/NxPs8xXHF2dWVaN5HSldQeGMhug+6HZMxSXsPnbOCzvK6t1w4L5vCCyVxYKMoCkskk75+jW6fE43EVi8W03sKKCtXe3q5b7vM4H3dMW5tSZWW86U+m5b4bSOlKare1qVQ4rH3m1oqYJbWFYxbxGj5foQ2fhXVtvXZYNIdxE7mwiPwduGKEiff0qBiRqhifAG65z+OOSSQuf2AmfkKhzLgfdQ3EnBr3mVvfxyypbSBm172Gz3m14bOlPkvGbeEcxiVzYZH5O1C3aSaIXrxIT+aMcT86MuJc/MQJorGx7LHR0cxOi37UldS2MWZJbcRsRhsxm9G2MWZJbQtjjkrmwiIJZDGSCIepNWeM+4nKSufiCxYQleXYGgpltnz2o66kto0xS2ojZjPaiNmMto0xS2pbGHNCMhcWSeCKEX5SeEVLC/UR0fzxMW65z+O5TxYXDW8wtGdP5oPCcLt7t/Ntn6V0EbM5bcSMmBFz6WpbFnNCOhcWySy+V0M+Z3h4WO8uyK8ZVVVVzViH353mV5b4SeFYLEYv7N9PxwcHaVFtLd21bt3kOG+f7HjXQn7amb9C48rVzR0QpXQFtdP9/dTR20urm5qovLGRXMVCPyRjFvMaPmcBnw3o2nrtsGAOTxvMhYXm79kUICKRCNXU1Og/d3d3U11dHR3v6KCGxYt1n19p4uN8nmN4AiW2YZbSldSuryfq7c20bmOjH5IxS3kNn7OBz/K6tl47LJjDiMlcWCCBKka4Ouvs7NTvTnO1l06nJ49Fo1FdBbL5fB4AAABQilT7MBcGqhhh2NzpDHZ8awYAAACwgGqf5cLAPcAKAAAAAH+BYgQAAAAAnoJiBAAAAACegmIEAAAAAJ4SuGIkHo/T7373u7zHeJyPzxR+j5rf384Hj890O2YpXUltG2OW1EbMZrQRsxltG2OW1LYx5rhgLpwRM9n45l/+5V9UY2OjmjNnjlqyZInqucqmOj/72c/Ue9/7Xn3+bbfdpp599llPNsrr7+9Xc+fOVbNnz1YvvfSSSqVSmU2YUind53E+zucVy/nz59Udd9yRd7fDid0R+Tif5wddkzFP9dmvMUtqm4zZLa/h8ztrw2e7fJaM28Y57BfMhcZ27T148KC65ppr1N69e9Uf//hH9dnPflZdd911anBwMO/5v/3tb1UoFFLf+c531Kuvvqq+9rWvqfLycvXKK6+4/h9zNY4ePapNZi1uDz/zjJ4AbqeO83nFkkgkJrdj1h8cLtBefDGzK+KUcT7PD7omY/7TkSPaZ279GrOktsmY3fIaPr+zNny2y2fJuG2cw6OCudBYMcLfhGzcuHGyPzo6qurr69XOnTvznv/JT35S3XPPPVljS5cuVZ/73OeMFyPMRNXHepHxbZO5nTCfj8+UiUpVf0CI1G/H28kPUk5l67WuqZgXjvvMrZ9jltQ2FbObXsPn6bXhs30+S8Zt4xy+JJgLxYuRS5cu6W85Dh06lDV+//33q3vvvTfv70SjUfVP//RPWWMPP/ywev/73z/t3/PWW2/pwCd+uOrj/5hz587pr5Gc/nD1x6Zff/31egK45T6PO9XmCp7/j1Mx5Yf7PO5HXRMxT/XZ7zHbPoduew2f4XMp+Wy7HxUua0vmwokfztuFFCNFbZQ3MDBADQ0NdPjwYVq2bNnk+JYtW/TysT09PVf8zjXXXEM/+clPqKWlZXLshz/8IT3yyCM0ODiY9+/ZsWOHPp7LgQMHqNKDrY0BAAAAUDwjIyO0du1aOzfK27ZtG23evDlr1z9eL//uu+92tGvvBP/v//5fWvWpT1F5RQXt3buXNmzYQOk336TnDh6k2//P/3GkfebYMVq9ciWdmjLG2zJ3vPii3oTIb7omYh6c4nPtm2/6OmZJbRMxu+01fM6vDZ/t9Fkybhvn8P8J5sKp+bsg/HibJhc8M2Lf/Ugb7/tKatt4jx0+T68Nn+3zWTJuG+fwJZufGZl4gHXTpk1ZD7A2NDS84wOsH/vYx7LGli1b5skDrMbfpunqknlS2wVdkzGLPxFvmR+SMYu+5QGf4bOlPkvGbeMcHi2Ft2n41V5eL2Tfvn36Vd0HHnhAv9p79uxZfby1tVVt3bo169Ve/o/63ve+p1577TW1fft2z17txToj3sQ81We/xiypbeO6DPD5nbXhs10+S8Zt4xz2l8I6I8wPfvAD9Z73vEevN8LflEz9Omf58uVq/fr1Vyx6dsstt+jzb731Vs8WPWPY3Ilqb+oEMDzuxHz+QExXofL4TD6MkrqS2lN1c332a8yS2qZidtNr+Dy9Nny2z+dcbRv9cFNbMhfOJH8X9TaNV/ADMNXV1Vd9GrdY0uk0dXR00OrVq6m8vNw1XZANfDYHvDYDfDYDfLbf50Lzd+D2pgEAAACAv0AxAgAAAABPQTECAAAAAE9BMQIAAAAAT0ExAgAAAABPQTECAAAAAE9BMQIAAAAAT0ExAgAAAABPQTECAAAAAE+ZTRYwsUhswVsRF7Hq3MjIiNbF6n5ywGdzwGszwGczwGf7fZ7I21db7N2KYuTChQu6jUajXocCAAAAgBnkcV4Wfjqs2JtmbGyMBgYGKBKJ0KxZs1yt2LjASSQSru55A7KBz+aA12aAz2aAz/b7zCUGFyL19fVUVlZm9zcj/B8wb948MX02Hx90eeCzOeC1GeCzGeCz3T6/0zciE+ABVgAAAAB4CooRAAAAAHhKoIuROXPm0Pbt23UL5IDP5oDXZoDPZoDPwfHZigdYAQAAAFC6BPqbEQAAAAB4D4oRAAAAAHgKihEAAAAAeAqKEQAAAAB4SskXI7t27aL58+fTtddeS0uXLqWjR4++4/k///nP6X3ve58+/y/+4i+oo6PDWKxB8fmJJ56gO++8k971rnfpn4985CNXnRcw88/0BAcPHtQrGN93333iMQbR5/Pnz9PGjRvpxhtv1G8l3HLLLbh+CPj82GOP0Xvf+16qqKjQq4Y+9NBD9NZbbxmL10b+4z/+gz7+8Y/rVVD5GtDe3n7V3+nu7qYPfOAD+rN888030759+2SDVCXMwYMH1TXXXKP27t2r/vjHP6rPfvaz6rrrrlODg4N5z//tb3+rQqGQ+s53vqNeffVV9bWvfU2Vl5erV155xXjspezz2rVr1a5du9Tvf/979dprr6lPf/rTqrq6Wp0+fdp47KXu9QR/+tOfVENDg7rzzjvVJz7xCWPxBsXnS5cuqdtvv12tXr1a/eY3v9F+d3d3q97eXuOxl7LP+/fvV3PmzNEte/zcc8+pG2+8UT300EPGY7eJjo4O9dWvflX94he/4Ldn1aFDh97x/L6+PlVZWak2b96sc+EPfvADnRs7OzvFYizpYmTJkiVq48aNk/3R0VFVX1+vdu7cmff8T37yk+qee+7JGlu6dKn63Oc+Jx5rkHzO5e2331aRSET95Cc/EYwyuF6zvx/84AdVW1ubWr9+PYoRAZ8ff/xxFYvFVCqVMhhl8Hzmc1euXJk1xgnzQx/6kHispQIVUIxs2bJF3XrrrVlja9asUatWrRKLq2Rv06RSKTp27Ji+BTB1jxvuHzlyJO/v8PjU85lVq1ZNez6Ymc+58NbVvIX19ddfLxhpcL3+h3/4B6qpqaG//du/NRRp8Hz+t3/7N1q2bJm+TVNbW0u33XYbffvb36bR0VGDkZe+zx/84Af170zcyunr69O3wlavXm0s7iBwxINcaMVGeTPh3Llz+kLAF4apcP+//uu/8v7O2bNn857P48A9n3P58pe/rO9l5n74gXOvf/Ob39CPfvQj6u3tNRRlMH3mpPjiiy/SunXrdHI8efIkPfjgg7rI5pUtgTs+r127Vv/ehz/8Yb0b7Ntvv02f//zn6Stf+YqhqIPB2WlyIe/u++abb+rnddymZL8ZAXbw6KOP6gcrDx06pB9gA+7B23a3trbqB4ZvuOEGr8MpacbGxvS3T3v27KHFixfTmjVr6Ktf/Sr967/+q9ehlRT8UCV/4/TDH/6QXn75ZfrFL35Bzz77LH3jG9/wOjTgkJL9ZoQvvqFQiAYHB7PGuV9XV5f3d3i8mPPBzHye4Hvf+54uRv793/+d3v/+9wtHGjyv/+d//odOnTqln6KfmjSZ2bNn0+uvv0433XSTgchL/zPNb9CUl5fr35tg4cKF+l+YfDvimmuuEY87CD5//etf1wX23/3d3+k+v/F48eJFeuCBB3Txx7d5gHOmy4VVVVUi34owJTtz/H9+/hfKCy+8kHUh5j7f280Hj089n/nVr3417flgZj4z3/nOd/S/Zjo7O+n22283FG2wvOZX1F955RV9i2bi595776W//uu/1n/m1yKBO5/pD33oQ/rWzESxx/z3f/+3LlJQiLjnMz9flltwTBSA2GbNPTzJharEXxvj18D27dunX0964IEH9GtjZ8+e1cdbW1vV1q1bs17tnT17tvre976nXzndvn07Xu0V8PnRRx/Vr/M9/fTT6s9//vPkz4ULFzz8ryhNr3PB2zQyPsfjcf1G2KZNm9Trr7+ufvnLX6qamhr1zW9+08P/itLzma/J7PNPf/pT/frp888/r2666Sb9JiSYHr628lIK/MNp//vf/77+c39/vz7OHrPXua/2fulLX9K5kJdiwKu9DuH3o9/znvfo5Mevkb300kuTx5YvX64vzlP52c9+pm655RZ9Pr/a9Oyzz3oQdWn73NjYqP8PkfvDFxrg/md6KihG5Hw+fPiwXgqAkyu/5vutb31Lv1YN3PM5nU6rHTt26ALk2muvVdFoVD344IPqf//3fz2K3g66urryXnMnvOWWvc79naamJj0v/Hn+8Y9/LBrjLP4fue9dAAAAAAAC+swIAAAAAOwAxQgAAAAAPAXFCAAAAAA8BcUIAAAAADwFxQgAAAAAPAXFCAAAAAA8BcUIAAAAADwFxQgAAAAAPAXFCAAAAAA8BcUIAAAAADwFxQgAAAAAPAXFCAAAAADIS/4/+/Acq5nZccUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Generate data, plot for sake of testing\n",
    "\n",
    "# Create domain for collocation points\n",
    "d = 0.05\n",
    "x = torch.arange(0, 1+d, d)\n",
    "y = torch.arange(0, 1+d, d)\n",
    "\n",
    "# Collocation tensor\n",
    "X = torch.stack(torch.meshgrid(x, y)).reshape(2,-1).T\n",
    "\n",
    "# Get boundary conditions\n",
    "b_left = torch.stack(torch.meshgrid(x[0], y)).reshape(2,-1).T  # x = 0, y = var, left wall\n",
    "b_right = torch.stack(torch.meshgrid(x[-1],y)).reshape(2,-1).T\n",
    "b_bottom = torch.stack(torch.meshgrid(x, y[0])).reshape(2,-1).T\n",
    "b_top = torch.stack(torch.meshgrid(x, y[-1])).reshape(2,-1).T\n",
    "\n",
    "X_train = torch.cat([b_left, b_bottom, b_right, b_top])\n",
    "print(X_train.size())\n",
    "\n",
    "## Get solutions for each boundary condition\n",
    "walls_Y = torch.zeros(b_left.size())\n",
    "\n",
    "# Get BC for lid driven component\n",
    "wall_lid = torch.zeros(b_top.size()[0], b_top.size()[1])\n",
    "wall_lid[:,0] = u0\n",
    "print(wall_lid)\n",
    "\n",
    "train_Y = torch.cat([walls_Y, walls_Y, walls_Y, wall_lid])\n",
    "print(train_Y.size())\n",
    "\n",
    "# Try and plot for test\n",
    "plt.figure()\n",
    "plt.scatter(X[:,0].numpy(), X[:,1].numpy(), marker='.', c='r')\n",
    "plt.scatter(b_left[:,0].numpy(), b_left[:,1].numpy(), marker='x', c='k')\n",
    "plt.scatter(b_right[:,0].numpy(), b_right[:,1].numpy(), marker='x', c='k')\n",
    "plt.scatter(b_top[:,0].numpy(), b_top[:,1].numpy(), marker='x', c='k')\n",
    "plt.scatter(b_bottom[:,0].numpy(), b_bottom[:,1].numpy(), marker='x', c='k')\n",
    "plt.grid()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1368356",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pytorch_cpu",
   "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.10.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}