{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 30,
   "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": 31,
   "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)\n",
    "\n",
    "    \n",
    "    def forward(self, x):\n",
    "        y = self.layers(x)\n",
    "        return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "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": 33,
   "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": 34,
   "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": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVXBJREFUeJztvX9sG+ed5/+26F8UTTENspQikWuDsdMa6da+deEf7ebIS9qE57Zp/jg0KxuKd3O32UVsXFGj6ybbbpxcd5tct9grrpf27HC9KQx7nW5Q04etTm0vEXlo7Mj5KpEvrXM9+yY2Z6y1WQOn0cBKTUX6fP8gxZj6kZAaPkM+5vsFEAM+JN/z9vsZzOdjzojPEhEREEIIIYQ0iLZGGyCEEEJIa8NmhBBCCCENhc0IIYQQQhoKmxFCCCGENBQ2I4QQQghpKGxGCCGEENJQ2IwQQgghpKGwGSGEEEJIQ1naaAPVMD09jdHRUQSDQSxZsqTRdgghhBBSBSICx3HQ3d2NtraFv//QohkZHR1FNBpttA1CCCGELALTNBGJRBZ8XYtmJBgMAij+Yzo6OuqmOzk5iZ/97Ge47777sGzZsrrpkkqYs3cwa29gzt7AnL1BZc7j4+OIRqPlOr4QWjQjM5dmOjo66t6MtLe3o6Ojgwe6QpizdzBrb2DO3sCcvcGLnD/sFgvewEoIIYSQhsJmhBBCCCENhc0IIYQQQhoKmxFCCCGENBQ2I4QQQghpKGxGCCGEENJQ2IwQQgghpKGwGSGEEEJIQ2m5ZiSXy+H111+f97XXX38duVxu0dq2bcOyrHlfsywLtm03la5KbR09q9SmZ2+06dkbbR09q9TW0bPKWrgopEay2ax8/vOfl9tvv10AyPHjxz/0M4ODg/Iv/sW/kOXLl8sdd9whf//3f1/TPm3bFgBi23atdiu4ePGirFq1SpYuXSqvvfaaFAoFSafTUigU5LXXXpOlS5fKqlWr5OLFizVrj42NydatWyUWi0kul6t4LZfLSSwWk61bt8rY2FhT6Hrp+cacm9WzSm0vPdcra+b8wdrMWa+cVfrWcQ5V1sLZVFu/a25G+vv75etf/7r8+Mc/rqoZMQxD2tvbZe/evXL27Fn53ve+Jz6fTwYGBqreZ72akdOnT8vSpUsFgCxdulROnjgh6XRaTp44UTF++vTpmrVN05RYLCYAigfO0JDIK69IbmioYtw0zabQ9dLzO6dOSTqdlndOnWpazyq1vfRcr6yZ8wdrM2e9clbpW8c5VFkLZ6OsGan4cBXNyL59++Suu+6qGHvooYfk/vvvr3o/9WpGRKTc9QGQoN8v6XRagn5/OfzXXntt0doznSoAiQHyamlbPpBmdbaN1vXK8/pSzutLOTerZ5XaXnmuZ9bMeWFt5qxfzip96ziHKmvhjTRNM3L33XfLl7/85YqxQ4cOSUdHx4Kf+e1vfyu2bZcfpmkKALl69aoUCgXXj5MnTkjQ75dbb71V0um03HrrrRL0++XkiROutd85dUrW+/3iv+Gx3u+Xd06dakpdLzzfmHOze9Z9DuudNXNmzjdTzrrnUW9tlbVw5nH16tWqmpElpaZiUSxZsgTHjx/Hgw8+uOB77rzzTvzxH/8xnnjiifJYf38/Pve5z2FiYgJ+v3/OZ5566ik8/fTTc8aPHj2K9vb2xdolhBBCiIdMTExgx44dsG0bHR0dC7+xuu9A5gdVfDOybt06+da3vlUx9pOf/EQAyMTExLyf4TcjN0fXreP/bnSfQ/6PnTkz55s3j5v5m5GmvEwzG94zot/1SB2v+6rU1vEaO3NeWJs565ezSt86zmHL3TOyb98++fjHP14x1tvb25AbWD3/a5rBQTV3atdB10vPyu+I1ywPlZ6V/pUHc2bOHnvW4dyh4xzeFH9N4ziOvPnmm/Lmm28KAPnbv/1befPNN8t/j/z4449LX19f+f0zf9r753/+5/L222/Lc88917A/7eXvjDTG8405N6tnldo6/i4Dc/5gbeasV84qfes4hzfF74wMDg4KSl8R3fjYtWuXiIjs2rVL4vH4nM9s3LhRli9fLrFYrGE/eiZSnISZbu/GCRApdotuwh8bG1uwQzVNc1EHo0pdldo36s7OuVk9q9T2ynM9s2bOC2szZ/1ynq09Gx3yqKe2ylp4I9XWb1d/TeMV4+PjCIVCH343bo1MTk6iv78f27dvx7Jly+qmSyphzt7BrL2BOXsDc/YGlTlXW79bbm0aQgghhDQXbEYIIYQQ0lDYjBBCCCGkobAZIYQQQkhDYTNCCCGEkIbScs2IbduwLGve1yzLgm3bHjsihBBCvKXZamFLNSO2bSOZTCIej8M0zYrXTNNEPB5HMplkQ0IIIeSmpRlrYUs1I47jIJ/PwzAMJBIJXBoeBgBcGh5GIpGAYRjI5/NwHMf9ziwLGBwsbuuJKl2V2qOjldt6omMeKj2rypo5V8Kc1esCep47NJhDT2thtdTlJ9YUU89fYFW5CFOZVEqkrU0EKG5TKfeaKnVVaqdSUggEir/uFwjo4VmltmLPSrJmznO0mbNi3ZK2ducOjebQk1ooHi2U5xX1bEZEpLjIECD+0gT4/X6JAcVFiNximu8fMDMPn6843oy6HngulHIu+P3N71mltgee6541c55Xmzl741mrc4eGc6i0Fpaotn631GWaGaLXruHwrLHDAKITE+7Fz50Dpqcrx6amgPPnm1NXpbaOnlVq07M32vTsjbaOnlVqa+hZaS2skZZsRsxAAH2zxvoAmO3t7sXXrQPaZsXq8wFr1zanrkptHT2r1KZnb7Tp2RttHT2r1NbQs9JaWCMt14yYpolEby8MAGtKY2sAGAASvb1z7iyumUgEOHiweKAAxe2BA8XxZtRVqa2jZ5Xa9OyNNj17o62jZ5XamnlWXgtrpKVW7bUsC/F4HIZhIBaL4eUjR3DmyhVs6OzEvTt3lsez2Swibg8gyyp+hbZ2bX0ORtW6CrUnL15E/8gItm/ciGWrV9dNF4CWeaj0rCxr5lwBc/ZAF5qeOzSYQy9rYbX1e6mrvWhGMBhEOBwGAGQyGXR1deFMfz96Nm1CJpNBIpFAOBxGMBh0v7NIpP4Hokpdldrd3cDISHFbb3TMQ6VnVVkz50qYs3pdQM9zhwZz6GktrJKWakZCoRAGBgbgOA4ikQgmJyfLr0WjUWSzWQSDQYRCoQa6JIQQQtTRjLWwpZoRoDgJCwXs+tIMIYQQogHNVgtb7gZWQgghhDQXbEYIIYQQ0lDYjBBCCCGkobAZIYQQQkhDYTNCCCGEkIbScs2IbduwFlh+2bIs2LbtsSNCCCHEW5qtFrZUM2LbNpLJJOLx+JyfujVNE/F4HMlkkg0JIYSQm5ZmrIUt1Yw4joN8Pg/DMJBIJHBpeBgAcGl4GIlEAoZhIJ/Pw3Ec9zuzLGBwsLitJ6p0VWqPjlZu64mOeaj0rCpr5lwJc1avC+h57tBgDj2thdUiGmDbtgAQ27Zda+VyOYnFYgJA1vv9kk6nZb3fLwAkFotJLpdzbziVEmlrEwGK21TKvaZKXZXaqZQUAgFJp9NSCAT08KxSW7FnJVkz5znazFmxbklbu3OHRnPoSS2U6ut3yzUjIiK5oSGJAeIvTYDf75cYILmhIffipvn+ATPz8PmK482o64HnQinngt/f/J5Vanvgue5ZM+d5tZmzN561OndoOIdKa2GJaut3S12mmSF67RoOzxo7DCA6MeFe/Nw5YHq6cmxqqrjSYjPqqtTW0bNKbXr2RpuevdHW0bNKbQ09K62FNdKSzYgZCKBv1lgfALO93b34unVA26xYfb7iks/NqKtSW0fPKrXp2RttevZGW0fPKrU19Ky0FtZIyzUjpmki0dsLA8Ca0tgaAAaARG/vnDuLayYSAQ4eLB4oQHF74ID7ZZ9V6arU1tGzSm169kabnr3R1tGzSm3NPCuvhTWyRETE0z0ugvHxcYRCIdi2jY6OjkXrWJaFeDwOwzAQi8Xw8pEjOHPlCjZ0duLenTvL49ls1v2qhZZV/Apt7dr6HIyqdRVqT168iP6REWzfuBHLVq+umy4ALfNQ6VlZ1sy5AubsgS40PXdoMIde1sJq6/dSV3vRjGAwiHA4DADIZDLo6urCmf5+9GzahEwmg0QigXA4jGAw6H5nkUj9D0SVuiq1u7uBkZHitt7omIdKz6qyZs6VMGf1uoCe5w4N5tDTWlglLdWMhEIhDAwMwHEcRCIRTE5Oll+LRqPIZrMIBoMIhUINdEkIIYSooxlrYUs1I0BxEhYK2PWlGUIIIUQDmq0WttwNrIQQQghpLtiMEEIIIaShsBkhhBBCSENhM0IIIYSQhsJmhBBCCCENpeWaEdu2YS2w/LJlWbBt22NHhBBCiLc0Wy1sqWbEtm0kk0nE4/E5P3Vrmibi8TiSySQbEkIIITctzVgLW6oZcRwH+XwehmEgkUjg0vAwAODS8DASiQQMw0A+n4fjOO53ZlnA4GBxW09U6arUHh2t3NYTHfNQ6VlV1sy5EuasXhfQ89yhwRx6WgurRTTAtm0BILZtu9bK5XISi8UEgKz3+yWdTst6v18ASCwWk1wu595wKiXS1iYCFLeplHtNlboqtVMpKQQCkk6npRAI6OFZpbZiz0qyZs5ztJmzYt2StnbnDo3m0JNaKNXX75ZrRkREckNDEgPEX5oAv98vMUByQ0PuxU3z/QNm5uHzFcebUdcDz4VSzgW/v/k9q9T2wHPds2bO82ozZ288a3Xu0HAOldbCEtXW75a6TDND9No1HJ41dhhAdGLCvfi5c8D0dOXY1FRxpcVm1FWpraNnldr07I02PXujraNnldoaelZaC2ukJZsRMxBA36yxPgBme7t78XXrgLZZsfp8xSWfm1FXpbaOnlVq07M32vTsjbaOnlVqa+hZaS2skZZrRkzTRKK3FwaANaWxNQAMAIne3jl3FtdMJAIcPFg8UIDi9sAB98s+q9JVqa2jZ5Xa9OyNNj17o62jZ5XamnlWXgtrZImIiKd7XATj4+MIhUKwbRsdHR2L1rEsC/F4HIZhIBaL4eUjR3DmyhVs6OzEvTt3lsez2az7VQstq/gV2tq19TkYVesq1J68eBH9IyPYvnEjlq1eXTddAFrmodKzsqyZcwXM2QNdaHru0GAOvayF1dbvpa72ohnBYBDhcBgAkMlk0NXVhTP9/ejZtAmZTAaJRALhcBjBYND9ziKR+h+IKnVVand3AyMjxW290TEPlZ5VZc2cK2HO6nUBPc8dGsyhp7WwSlqqGQmFQhgYGIDjOIhEIpicnCy/Fo1Gkc1mEQwGEQqFGuiSEEIIUUcz1sKWakaA4iQsFLDrSzOEEEKIBjRbLWy5G1gJIYQQ0lywGSGEEEJIQ2EzQgghhJCGwmaEEEIIIQ2FzQghhBBCGkrLNSO2bcNaYPlly7Jg27bHjgghhBBvabZauKhm5LnnnsOaNWuwcuVKbNmyBadPn/7A93/3u9/FRz/6Ufj9fkSjUXzlK1/Bb3/720UZdoNt20gmk4jH43N+6tY0TcTjcSSTSTYkhBBCblqasRbW3Iy8+OKL2Lt3L/bv34833ngDGzZswP333498Pj/v+48ePYrHH38c+/fvx9tvv42/+7u/w4svvoi/+Iu/cG2+VhzHQT6fh2EYSCQSuDQ8DAC4NDyMRCIBwzCQz+fhOI77nVkWMDhY3NYTVboqtUdHK7f1RMc8VHpWlTVzroQ5q9cF9Dx3aDCHntbCapEa2bx5s+zevbv8fGpqSrq7u+WZZ56Z9/27d++We+65p2Js79698ulPf7rqfdq2LQDEtu1a7c4hl8tJLBYTALLe75d0Oi3r/X4BILFYTHK5nOt9SCol0tYmAhS3qZR7TZW6KrVTKSkEApJOp6UQCOjhWaW2Ys9KsmbOc7SZs2LdkrZ25w6N5tCTWijV1++aFsorFApob2/HSy+9hAcffLA8vmvXLoyNjeHEiRNzPnP06FE89thj+NnPfobNmzfDMAx87nOfQ19f34Lfjly/fh3Xr18vPx8fH0c0GsXVq1ddLZQ3w6XhYWy/5x5c8ftx6NAhPPLII+h89130v/IKejZtcic+OgrcdRcwPf3+mM8H/PKX7tZXUKWrUrukO7liBX5+6BA++8gjWFYoNLdnldoeeK571sx5Xm3mrFD3Bm2tzh0azqHSWlhifHwct91224culFdTMzI6Ooqenh6cPHkS27ZtK4/v27cP2WwWQ0ND837uP//n/4yvfvWrEBG89957+LM/+zP84Ac/WHA/Tz31FJ5++uk540ePHkV7e3u1dgkhhBDSQCYmJrBjx47Gr9qbyWTwrW99C9///vexZcsWnD9/Hl/+8pfxzW9+E3/5l38572eeeOIJ7N27t/x85puR++67j9+MaNR1a/m/G5XaOv5PkjnPq82cFereoK3VuUPDOfTqm5GqqOXaz/Xr18Xn88nx48crxh9++GF54IEH5v3MH/zBH8hXv/rVirHDhw+L3++Xqampqvar5T0jPl/x2p7PV9/rkSp0VWqnUlJYtap43XfVKj08q9RW7FlJ1sx5jjZzVqxb0tbu3KHRHDbbPSOLuoF1z5495edTU1PS09Oz4A2sv//7vy/79u2rGDt69Kj4/X557733qtpnvZoR0zTL4cdiMXnn1ClJp9PyzqlTFeOmabraT2lnIoODxW09UaWrULtw4ULxhHLhQl11RUTLPFR6VpY1c66AOXugK5qeOzSYQy9robJm5NixY7JixQp54YUX5OzZs/Loo4/KLbfcIpcvXxYRkb6+Pnn88cfL79+/f78Eg0H5h3/4BzEMQ372s5/JHXfcIV/60pfq/o/5MMbGxmTr1q3lrq9QKBQP9EKh3CVu3bpVxsbGXO2HVHJjzkQtzNobmLM3MGc1eFkLq63fNd8z8tBDD+E3v/kNnnzySVy+fBkbN27EwMAAOjs7AQC5XA5tbe//fMk3vvENLFmyBN/4xjdw6dIl/M7v/A6+8IUv4K//+q9r3bVrQqEQBgYG4DgOIpEIJicny69Fo1Fks1kEg0GEQiHPvRFCCCFe0Iy1cFE3sO7Zswd79uyZ97VMJlO5g6VLsX//fuzfv38xu6o7oVBowYAjkYjHbgghhBDvabZa2HJr0xBCCCGkuWAzQgghhJCGwmaEEEIIIQ2FzQghhBBCGkrLNSO2bcNaYMVDy7I8XTKZEEIIaQTNVgtbqhmxbRvJZBLxeBymaVa8Zpom4vE4kskkGxJCCCE3Lc1YC1uqGXEcB/l8HoZhIJFI4NLwMIDi7/MnEgkYhoF8Pg/HcdzvzLKAwcHitp6o0lWpPTpaua0nOuah0rOqrJlzJcxZvS6g57lDgzn0tBZWi+ufV/MALdemaWsrriHQ1lbfdQ9U6KrUTqWkEAgUf90vENDDs0ptxZ6VZM2c52gzZ8W6JW3tzh0azaH2a9M0gno2IyIiuaEhiQHiL02A3++XGCC5oSH34qb5/gEz8/D53K9ToErXA8+FUs4Fv7/5PavU9sBz3bNmzvNqM2dvPGt17tBwDpXWwhLV1u+WukwzQ/TaNRyeNXYYQHRiwr34uXOVyzwDwNQUcP58c+qq1NbRs0ptevZGm5690dbRs0ptDT0rrYU10pLNiBkIoG/WWB8As73dvfi6dUDbrFh9PmDt2ubUVamto2eV2vTsjTY9e6Oto2eV2hp6VloLa6TlmhHTNJHo7YUBYE1pbA0AA0Cit3fOncU1E4kABw8WDxSguD1woDjejLoqtXX0rFKbnr3RpmdvtHX0rFJbM8/Ka2GNLBER8XSPi2B8fByhUAi2baOjo2PROpZlIR6PwzAMxGIxvHzkCM5cuYINnZ24d+fO8ng2m3W/UJBlFb9CW7u2Pgejal2F2pMXL6J/ZATbN27EstWr66YLQMs8VHpWljVzroA5e6ALTc8dGsyhl7Ww2vq9qFV7dSUYDCIcDgMori7c1dWFM/396Nm0CZlMBolEAuFwGMFg0P3OIpH6H4gqdVVqd3cDIyPFbb3RMQ+VnlVlzZwrYc7qdQE9zx0azKGntbBKWqoZCYVCGBgYgOM4iEQimJycLL8WjUaRzWYRDAYXXFaZEEII0Z1mrIUt1YwAxUlYKGDXl2YIIYQQDWi2WthyN7ASQgghpLlgM0IIIYSQhsJmhBBCCCENhc0IIYQQQhoKmxFCCCGENJSWa0Zs24a1wPLLlmXBtm2PHRFCCCHe0my1sKWaEdu2kUwmEY/H5/zUrWmaiMfjSCaTbEgIIYTctDRjLWypZsRxHOTzeRiGgUQigUvDwwCAS8PDSCQSMAwD+XwejuO435llAYODxW09UaWrUnt0tHJbT3TMQ6VnVVkz50qYs3pdQM9zhwZz6GktrBbRANu2BYDYtu1aK5fLSSwWEwCy3u+XdDot6/1+ASCxWExyuZx7w6mUSFubCFDcplLuNVXqqtROpaQQCEg6nZZCIKCHZ5Xaij0ryZo5z9Fmzop1S9ranTs0mkNPaqFUX79brhkREckNDUkMEH9pAvx+v8QAyQ0NuRc3zfcPmJmHz1ccb0ZdDzwXSjkX/P7m96xS2wPPdc+aOc+rzZy98azVuUPDOVRaC0tUW79b6jLNDNFr13B41thhANGJCffi584B09OVY1NTxZUWm1FXpbaOnlVq07M32vTsjbaOnlVqa+hZaS2skZZsRsxAAH2zxvoAmO3t7sXXrQPaZsXq8xWXfG5GXZXaOnpWqU3P3mjTszfaOnpWqa2hZ6W1sEZarhkxTROJ3l4YANaUxtYAMAAkenvn3FlcM5EIcPBg8UABitsDB9wv+6xKV6W2jp5VatOzN9r07I22jp5VamvmWXktrJElIiKe7nERjI+PIxQKwbZtdHR0LFrHsizE43EYhoFYLIaXjxzBmStXsKGzE/fu3Fkez2az7lcttKziV2hr19bnYFStq1B78uJF9I+MYPvGjVi2enXddAFomYdKz8qyZs4VMGcPdKHpuUODOfSyFlZbv5e62otmBINBhMNhAEAmk0FXVxfO9PejZ9MmZDIZJBIJhMNhBINB9zuLROp/IKrUVand3Q2MjBS39UbHPFR6VpU1c66EOavXBfQ8d2gwh57WwippqWYkFAphYGAAjuMgEolgcnKy/Fo0GkU2m0UwGEQoFGqgS0IIIUQdzVgLW6oZAYqTsFDAri/NEEIIIRrQbLWw5W5gJYQQQkhzwWaEEEIIIQ2FzQghhBBCGgqbEUIIIYQ0FDYjhBBCCGkoLdeM2LYNa4Hlly3Lgm3bHjsihBBCvKXZamFLNSO2bSOZTCIej8/5qVvTNBGPx5FMJtmQEEIIuWlpxlrYUs2I4zjI5/MwDAOJRAKXhocBAJeGh5FIJGAYBvL5PBzHcb8zywIGB4vbeqJKV6X26Gjltp7omIdKz6qyZs6VMGf1uoCe5w4N5tDTWlgtogG2bQsAsW3btVYul5NYLCYAZL3fL+l0Wtb7/QJAYrGY5HI594ZTKZG2NhGguE2l3Guq1FWpnUpJIRCQdDothUBAD88qtRV7VpI1c56jzZwV65a0tTt3aDSHntRCqb5+t1wzIiKSGxqSGCD+0gT4/X6JAZIbGnIvbprvHzAzD5+vON6Muh54LpRyLvj9ze9ZpbYHnuueNXOeV5s5e+NZq3OHhnOotBaWqLZ+t9Rlmhmi167h8KyxwwCiExPuxc+dA6anK8empoorLTajrkptHT2r1KZnb7Tp2RttHT2r1NbQs9JaWCMt2YyYgQD6Zo31ATDb292Lr1sHtM2K1ecrLvncjLoqtXX0rFKbnr3RpmdvtHX0rFJbQ89Ka2GNtFwzYpomEr29MACsKY2tAWAASPT2zrmzuGYiEeDgweKBAhS3Bw64X/ZZla5KbR09q9SmZ2+06dkbbR09q9TWzLPyWlgjS0REPN3jIhgfH0coFIJt2+jo6Fi0jmVZiMfjMAwDsVgMLx85gjNXrmBDZyfu3bmzPJ7NZt2vWmhZxa/Q1q6tz8GoWleh9uTFi+gfGcH2jRuxbPXquukC0DIPlZ6VZc2cK2DOHuhC03OHBnPoZS2stn4vdbUXzQgGgwiHwwCATCaDrq4unOnvR8+mTchkMkgkEgiHwwgGg+53FonU/0BUqatSu7sbGBkpbuuNjnmo9Kwqa+ZcCXNWrwvoee7QYA49rYVV0lLNSCgUwsDAABzHQSQSweTkZPm1aDSKbDaLYDCIUCjUQJeEEEKIOpqxFrZUMwIUJ2GhgF1fmiGEEEI0oNlqYcvdwEoIIYSQ5oLNCCGEEEIaCpsRQgghhDQUNiOEEEIIaShsRgghhBDSUFquGbFtG9YCyy9blgXbtj12RAghhHhLs9XCRTUjzz33HNasWYOVK1diy5YtOH369Ae+f2xsDLt378btt9+OFStW4M4770R/f/+iDLvBtm0kk0nE4/E5P3Vrmibi8TiSySQbEkIIITctzVgLa25GXnzxRezduxf79+/HG2+8gQ0bNuD+++9HPp+f9/2FQgGf/exnceHCBbz00kv49a9/jeeffx49PT2uzdeK4zjI5/MwDAOJRAKXhocBAJeGh5FIJGAYBvL5PBzHcb8zywIGB4vbeqJKV6X26Gjltp7omIdKz6qyZs6VMGf1uoCe5w4N5tDTWlgtUiObN2+W3bt3l59PTU1Jd3e3PPPMM/O+/wc/+IHEYjEpFAq17qqMbdsCQGzbXrTGDLlcTmKxmACQ9X6/pNNpWe/3CwCJxWKSy+Vc70NSKZG2NhGguE2l3Guq1FWpnUpJIRCQdDothUBAD88qtRV7VpI1c56jzZwV65a0tTt3aDSHntRCqb5+17RQXqFQQHt7O1566SU8+OCD5fFdu3ZhbGwMJ06cmPOZ7du349Zbb0V7eztOnDiB3/md38GOHTvwta99Db6ZFQhncf36dVy/fr38fHx8HNFoFFevXnW1UN4Ml4aHsf2ee3DF78ehQ4fwyCOPoPPdd9H/yivo2bTJnfjoKHDXXcD09PtjPh/wy1+6W19Bla5K7ZLu5IoV+PmhQ/jsI49gWaHQ3J5Vanvgue5ZM+d5tZmzQt0btLU6d2g4h0prYYnx8XHcdtttH7pQXk3NyOjoKHp6enDy5Els27atPL5v3z5ks1kMDQ3N+czHPvYxXLhwATt37sRjjz2G8+fP47HHHsO///f/Hvv37593P0899RSefvrpOeNHjx5Fe3t7tXYJIYQQ0kAmJiawY8eOxq/aOz09jXA4jIMHD8Ln82HTpk24dOkS/uZv/mbBZuSJJ57A3r17y89nvhm57777+M2IRl23lv+7Uamt4/8kmfO82sxZoe4N2lqdOzScQ6++GamKWq79XL9+XXw+nxw/frxi/OGHH5YHHnhg3s/8y3/5L+Xee++tGOvv7xcAcv369ar2q+U9Iz5f8dqez1ff65EqdFVqp1JSWLWqeN131So9PKvUVuxZSdbMeY42c1asW9LW7tyh0Rw22z0ji7qBdc+ePeXnU1NT0tPTs+ANrE888YSsXr1apqamymPf/e535fbbb696n/VqRkzTLIcfi8XknVOnJJ1OyzunTlWMm6bpaj+lnYkMDha39USVrkLtwoULxRPKhQt11RURLfNQ6VlZ1sy5Aubsga5oeu7QYA69rIXKmpFjx47JihUr5IUXXpCzZ8/Ko48+KrfccotcvnxZRET6+vrk8ccfL78/l8tJMBiUPXv2yK9//Wv5p3/6JwmHw/JXf/VXdf/HfBhjY2OydevWctdXKBSKB3qhUO4St27dKmNjY672Qyq5MWeiFmbtDczZG5izGryshdXW75rvGXnooYfwm9/8Bk8++SQuX76MjRs3YmBgAJ2dnQCAXC6Htrb3f74kGo3ipz/9Kb7yla/gE5/4BHp6evDlL38ZX/va12rdtWtCoRAGBgbgOA4ikQgmJycrfGazWQSDQYRCIc+9EUIIIV7QjLVwUTew7tmzB3v27Jn3tUwmM2ds27ZteO211xazq7oTCoUWDDgSiXjshhBCCPGeZquFLbc2DSGEEEKaCzYjhBBCCGkobEYIIYQQ0lDYjBBCCCGkobAZIYQQQkhDablmxLZtWAssv2xZFmzb9tgRIYQQ4i3NVgtbqhmxbRvJZBLxeBymaVa8Zpom4vE4kskkGxJCCCE3Lc1YC1uqGXEcB/l8HoZhIJFI4NLwMIDiYkGJRAKGYSCfz8NxHPc7syxgcLC4rSeqdFVqj45WbuuJjnmo9Kwqa+ZcCXNWrwvoee7QYA49rYXV4vq3Xj1Ay4Xy2tqKCxq1tdV3ESYVuiq1UykpBALFnxoOBPTwrFJbsWclWTPnOdrMWbFuSVu7c4dGc6j9QnmNoJ7NiIhIbmhIYoD4SxPg9/slBkhuaMi9uGm+f8DMPHw+94smqdL1wHOhlHPB729+zyq1PfBc96yZ87zazNkbz1qdOzScQ6W1sES19bulLtPMEL12DYdnjR0GEJ2YcC9+7hwwPV05NjUFnD/fnLoqtXX0rFKbnr3RpmdvtHX0rFJbQ89Ka2GNtGQzYgYC6Js11gfAbG93L75uHdA2K1afD1i7tjl1VWrr6FmlNj17o03P3mjr6FmltoaeldbCGmm5ZsQ0TSR6e2EAWFMaWwPAAJDo7Z1zZ3HNRCLAwYPFAwUobg8cKI43o65KbR09q9SmZ2+06dkbbR09q9TWzLPyWlgjS0REPN3jIhgfH0coFIJt2+jo6Fi0jmVZiMfjMAwDsVgMLx85gjNXrmBDZyfu3bmzPJ7NZt2vWmhZxa/Q1q6tz8GoWleh9uTFi+gfGcH2jRuxbPXquukC0DIPlZ6VZc2cK2DOHuhC03OHBnPoZS2stn4vdbUXzQgGgwiHwwCATCaDrq4unOnvR8+mTchkMkgkEgiHwwgGg+53FonU/0BUqatSu7sbGBkpbuuNjnmo9Kwqa+ZcCXNWrwvoee7QYA49rYVV0lLNSCgUwsDAABzHQSQSweTkZPm1aDSKbDaLYDCIUCjUQJeEEEKIOpqxFrZUMwIUJ2GhgF1fmiGEEEI0oNlqYcvdwEoIIYSQ5oLNCCGEEEIaCpsRQgghhDQUNiOEEEIIaSgt14zYtg1rgRUPLcvydMlkQgghpBE0Wy1sqWbEtm0kk0nE4/E5vy5nmibi8TiSySQbEkIIITctzVgLW6oZcRwH+XwehmEgkUjg0vAwAODS8DASiQQMw0A+n4fjOO53ZlnA4GBxW09U6arUHh2t3NYTHfNQ6VlV1sy5EuasXhfQ89yhwRx6WgurpW7rBCuk2iWIqyGXy0ksFhMAsr60bPJ6v18ASCwWk1wu595wKvX+cs9tbcXn9UCVrkrtVEoKgUBxGfBAQA/PKrUVe1aSNXOeo82cFeuWtLU7d2g0h57UQqm+frdcMyIikhsakhgg/tIE+P1+iQGSGxpyL26a7x8wMw+frzjejLoeeC6Uci74/c3vWaW2B57rnjVznlebOXvjWatzh4ZzqLQWlqi2frfUZZoZoteu4fCsscMAohMT7sXPnQOmpyvHpqaKixs1o65KbR09q9SmZ2+06dkbbR09q9TW0LPSWlgjLdmMmIEA+maN9QEw29vdi69bB7TNitXnK66y2Iy6KrV19KxSm5690aZnb7R19KxSW0PPSmthjbRcM2KaJhK9vTAArCmNrQFgAEj09s65s7hmIhHg4MHigQIUtwcOuF9pUZWuSm0dPavUpmdvtOnZG20dPavU1syz8lpYI0tERDzd4yIYHx9HKBSCbdvo6OhYtI5lWYjH4zAMA7FYDC8fOYIzV65gQ2cn7t25szyezWbdLxRkWcWv0Naure9y0qp0FWpPXryI/pERbN+4EctWr66bLgAt81DpWVnWzLkC5uyBLjQ9d2gwh17Wwmrrd0ut2hsMBhEOhwEAmUwGXV1dONPfj55Nm5DJZJBIJBAOhxEMBt3vLBKp/4GoUleldnc3MDJS3NYbHfNQ6VlV1sy5EuasXhfQ89yhwRx6WgurpKWakVAohIGBATiOg0gkgsnJyfJr0WgU2WwWwWBwwWWVCSGEEN1pxlrYUs0IUJyEhQJ2fWmGEEII0YBmq4UtdwMrIYQQQpoLNiOEEEIIaShsRgghhBDSUNiMEEIIIaShsBkhhBBCSENpuWbEtm1YCyy/bFkWbNv22BEhhBDiLc1WC1uqGbFtG8lkEvF4fM5P3ZqmiXg8jmQyyYaEEELITUsz1sKWakYcx0E+n4dhGEgkErg0PAwAuDQ8jEQiAcMwkM/n4TiO+51ZFjA4WNzWE1W6KrVHRyu39UTHPFR6VpU1c66EOavXBfQ8d2gwh57WwmoRDbBtWwCIbduutXK5nMRiMQEg6/1+SafTst7vFwASi8Ukl8u5N5xKibS1iQDFbSrlXlOlrkrtVEoKgYCk02kpBAJ6eFaprdizkqyZ8xxt5qxYt6St3blDozn0pBZK9fW75ZoREZHc0JDEAPGXJsDv90sMkNzQkHtx03z/gJl5+HzF8WbU9cBzoZRzwe9vfs8qtT3wXPesmfO82szZG89anTs0nEOltbBEtfW7pS7TzBC9dg2HZ40dBhCdmHAvfu4cMD1dOTY1VVxpsRl1VWrr6FmlNj17o03P3mjr6FmltoaeldbCGmnJZsQMBNA3a6wPgNne7l583TqgbVasPl9xyedm1FWpraNnldr07I02PXujraNnldoaelZaC2uk5ZoR0zSR6O2FAWBNaWwNAANAord3zp3FNROJAAcPFg8UoLg9cMD9ss+qdFVq6+hZpTY9e6NNz95o6+hZpbZmnpXXwhpZIiLi6R4Xwfj4OEKhEGzbRkdHx6J1LMtCPB6HYRiIxWJ4+cgRnLlyBRs6O3Hvzp3l8Ww2637VQssqfoW2dm19DkbVugq1Jy9eRP/ICLZv3Ihlq1fXTReAlnmo9Kwsa+ZcAXP2QBeanjs0mEMva2G19Xupq71oRjAYRDgcBgBkMhl0dXXhTH8/ejZtQiaTQSKRQDgcRjAYdL+zSKT+B6JKXZXa3d3AyEhxW290zEOlZ1VZM+dKmLN6XUDPc4cGc+hpLaySlmpGQqEQBgYG4DgOIpEIJicny69Fo1Fks1kEg0GEQqEGuiSEEELU0Yy1sKWaEaA4CQsF7PrSDCGEEKIBzVYLW+4GVkIIIYQ0F2xGCCGEENJQ2IwQQgghpKGwGSGEEEJIQ2EzQgghhJCG0nLNiG3bsBZYftmyLNi27bEjQgghxFuarRYuqhl57rnnsGbNGqxcuRJbtmzB6dOnq/rcsWPHsGTJEjz44IOL2a1rbNtGMplEPB6f81O3pmkiHo8jmUyyISGEEHLT0oy1sOZm5MUXX8TevXuxf/9+vPHGG9iwYQPuv/9+5PP5D/zchQsX8NWvfhV33333os26xXEc5PN5GIaBRCKBS8PDAIBLw8NIJBIwDAP5fB6O47jfmWUBg4PFbT1RpatSe3S0cltPdMxDpWdVWTPnSpizel1Az3OHBnPoaS2sFqmRzZs3y+7du8vPp6ampLu7W5555pkFP/Pee+/Jpz71KUmlUrJr1y754he/WNM+bdsWAGLbdq1255DL5SQWiwkAWe/3SzqdlvV+vwCQWCwmuVzO9T4klRJpaxMBittUyr2mSl2V2qmUFAIBSafTUggE9PCsUluxZyVZM+c52sxZsW5JW7tzh0Zz6EktlOrrd00L5RUKBbS3t+Oll16quNSya9cujI2N4cSJE/N+bv/+/fhf/+t/4fjx4/ijP/ojjI2NIZ1OL7if69ev4/r16+Xn4+PjiEajuHr1qquF8ma4NDyM7ffcgyt+Pw4dOoRHHnkEne++i/5XXkHPpk3uxEdHgbvuAqan3x/z+YBf/tLd+gqqdFVql3QnV6zAzw8dwmcfeQTLCoXm9qxS2wPPdc+aOc+rzZwV6t6grdW5Q8M5VFoLS4yPj+O222770IXyampGRkdH0dPTg5MnT2Lbtm3l8X379iGbzWJoaGjOZ37xi1/gD//wDzEyMoLbbrutqmbkqaeewtNPPz1n/OjRo2hvb6/WLiGEEEIayMTEBHbs2NHYVXsdx0FfXx+ef/553HbbbVV/7oknnsDevXvLz2e+Gbnvvvv4zYhGXbeW/7tRqa3j/ySZ87zazFmh7g3aWp07NJxDr74ZqYparv1cv35dfD6fHD9+vGL84YcflgceeGDO+998800BID6fr/xYsmSJLFmyRHw+n5w/f76q/Wp5z4jPV7y25/PV93qkCl2V2qmUFFatKl73XbVKD88qtRV7VpI1c56jzZwV65a0tTt3aDSHzXbPyKJuYN2zZ0/5+dTUlPT09Mx7A+u7774rb731VsXji1/8otxzzz3y1ltvyfXr16vaZ72aEdM0y+HHYjF559QpSafT8s6pUxXjpmm62k9pZyKDg8VtPVGlq1C7cOFC8YRy4UJddUVEyzxUelaWNXOugDl7oCuanjs0mEMva2G19bvmyzR79+7Frl278MlPfhKbN2/Gd7/7XVy7dg1//Md/DAB4+OGH0dPTg2eeeQYrV67Exz/+8YrP33LLLQAwZ9wLgsEgwuEwACCTyaCrqwtn+vvRs2kTMpkMEokEwuEwgsGg+51FIsVHvVGlq1K7uxsYGXH/VeV86JiHSs+qsmbOlTBn9bqAnucODebQ01pYJTU3Iw899BB+85vf4Mknn8Tly5exceNGDAwMoLOzEwCQy+XQ1tacP+waCoUwMDAAx3EQiUQwOTlZfi0ajSKbzSIYDCIUCjXQJSGEEKKOZqyFi7qBdc+ePdizZ8+8r2UymQ/87AsvvLCYXdaNUCi0YMARVd0sIYQQ0kQ0Wy1szq8wCCGEENIysBkhhBBCSENhM0IIIYSQhsJmhBBCCCENhc0IIYQQQhpKyzUjtm3DWmD5ZcuyYNu2x44IIYQQb2m2WthSzYht20gmk4jH4zBNs+I10zQRj8eRTCbZkBBCCLlpacZa2FLNiOM4yOfzMAwDiUQCl4aHARQXC0okEjAMA/l8Ho7juN+ZZQGDg8VtPVGlq1J7dLRyW090zEOlZ1VZM+dKmLN6XUDPc4cGc+hpLawW1z887wFaLpTX1lZc0Kitrb6LMKnQVamdSkkhECiuLxEI6OFZpbZiz0qyZs5ztJmzYt2StnbnDo3mUPuF8hpBPZsREZHc0JDEAPGXJsDv90sMkNzQkHtx03z/gJl5+HzuF01SpeuB50Ip54Lf3/yeVWp74LnuWTPnebWZszeetTp3aDiHSmthiWrrd0tdppkheu0aDs8aOwwgOjHhXvzcOWB6unJsago4f745dVVq6+hZpTY9e6NNz95o6+hZpbaGnpXWwhppyWbEDATQN2usD4DZ3u5efN06YPZCgT4fsHZtc+qq1NbRs0ptevZGm5690dbRs0ptDT0rrYU10nLNiGmaSPT2wgCwpjS2BoABINHbO+fO4pqJRICDB4sHClDcHjjgftlnVboqtXX0rFKbnr3RpmdvtHX0rFJbM8/Ka2GNLBER8XSPi2B8fByhUAi2baOjo2PROpZlIR6PwzAMxGIxvHzkCM5cuYINnZ24d+fO8ng2m3W/aqFlFb9CW7u2Pgejal2F2pMXL6J/ZATbN27EstWr66YLQMs8VHpWljVzroA5e6ALTc8dGsyhl7Ww2vq91NVeNCMYDCIcDgMAMpkMurq6cKa/Hz2bNiGTySCRSCAcDiMYDLrfWSRS/wNRpa5K7e5uYGSkuK03Ouah0rOqrJlzJcxZvS6g57lDgzn0tBZWSUs1I6FQCAMDA3AcB5FIBJOTk+XXotEostksgsEgQqFQA10SQggh6mjGWthSzQhQnISFAnZ9aYYQQgjRgGarhS13AyshhBBCmgs2I4QQQghpKGxGCCGEENJQ2IwQQgghpKGwGSGEEEJIQ2m5ZsS2bVgLLL9sWRZs2/bYESGEEOItzVYLW6oZsW0byWQS8Xh8zk/dmqaJeDyOZDLJhoQQQshNSzPWwpZqRhzHQT6fh2EYSCQSuDQ8DAC4NDyMRCIBwzCQz+fhOI77nVkWMDhY3NYTVboqtUdHK7f1RMc8VHpWlTVzroQ5q9cF9Dx3aDCHntbCahENsG1bAIht2661crmcxGIxASDr/X5Jp9Oy3u8XABKLxSSXy7k3nEqJtLWJAMVtKuVeU6WuSu1USgqBgKTTaSkEAnp4Vqmt2LOSrJnzHG3mrFi3pK3duUOjOfSkFkr19bvlmhERkdzQkMQA8ZcmwO/3SwyQ3NCQe3HTfP+AmXn4fMXxZtT1wHOhlHPB729+zyq1PfBc96yZ87zazNkbz1qdOzScQ6W1sES19bulLtPMEL12DYdnjR0GEJ2YcC9+7hwwPV05NjVVXGmxGXVVauvoWaU2PXujTc/eaOvoWaW2hp6V1sIaaclmxAwE0DdrrA+A2d7uXnzdOqBtVqw+X3HJ52bUVamto2eV2vTsjTY9e6Oto2eV2hp6VloLa6TlmhHTNJHo7YUBYE1pbA0AA0Cit3fOncU1E4kABw8WDxSguD1wwP2yz6p0VWrr6FmlNj17o03P3mjr6FmltmaeldfCGlkiIuLpHhfB+Pg4QqEQbNtGR0fHonUsy0I8HodhGIjFYnj5yBGcuXIFGzo7ce/OneXxbDbrftVCyyp+hbZ2bX0ORtW6CrUnL15E/8gItm/ciGWrV9dNF4CWeaj0rCxr5lwBc/ZAF5qeOzSYQy9rYbX1e6mrvWhGMBhEOBwGAGQyGXR1deFMfz96Nm1CJpNBIpFAOBxGMBh0v7NIpP4Hokpdldrd3cDISHFbb3TMQ6VnVVkz50qYs3pdQM9zhwZz6GktrJKWakZCoRAGBgbgOA4ikQgmJyfLr0WjUWSzWQSDQYRCoQa6JIQQQtTRjLWwpZoRoDgJCwXs+tIMIYQQogHNVgtb7gZWQgghhDQXbEYIIYQQ0lDYjBBCCCGkobAZIYQQQkhDablmxLZtWAuseGhZlqdLJhNCCCGNoNlqYUs1I7ZtI5lMIh6Pz/l1OdM0EY/HkUwm2ZAQQgi5aWnGWthSzYjjOMjn8zAMA4lEApeGhwEAl4aHkUgkYBgG8vk8HMdxvzPLAgYHi9t6okpXpfboaOW2nuiYh0rPqrJmzpUwZ/W6gJ7nDg3m0NNaWC11WydYIdUuQVwNuVxOYrGYAJD1pWWT1/v9AkBisZjkcjn3hlOp95d7bmsrPq8HqnRVaqdSUggEisuABwJ6eFaprdizkqyZ8xxt5qxYt6St3blDozn0pBZK9fW75ZoREZHc0JDEAPGXJsDv90sMkNzQkHtx03z/gJl5+HzF8WbU9cBzoZRzwe9vfs8qtT3wXPesmfO82szZG89anTs0nEOltbBEtfW7pS7TzBC9dg2HZ40dBhCdmHAvfu4cMD1dOTY1VVzcqBl1VWrr6FmlNj17o03P3mjr6FmltoaeldbCGmnJZsQMBNA3a6wPgNne7l583TqgbVasPl9xlcVm1FWpraNnldr07I02PXujraNnldoaelZaC2uk5ZoR0zSR6O2FAWBNaWwNAANAord3zp3FNROJAAcPFg8UoLg9cMD9SouqdFVq6+hZpTY9e6NNz95o6+hZpbZmnpXXwhpZIiLi6R4Xwfj4OEKhEGzbRkdHx6J1LMtCPB6HYRiIxWJ4+cgRnLlyBRs6O3Hvzp3l8Ww2636hIMsqfoW2dm19l5NWpatQe/LiRfSPjGD7xo1Ytnp13XQBaJmHSs/KsmbOFTBnD3Sh6blDgzn0shZWW79batXeYDCIcDgMAMhkMujq6sKZ/n70bNqETCaDRCKBcDiMYDDofmeRSP0PRJW6KrW7u4GRkeK23uiYh0rPqrJmzpUwZ/W6gJ7nDg3m0NNaWCUt1YyEQiEMDAzAcRxEIhFMTk6WX4tGo8hmswgGgwsuq0wIIYToTjPWwpZqRoDiJCwUsOtLM4QQQogGNFstbLkbWAkhhBDSXLAZIYQQQkhDYTNCCCGEkIbCZoQQQgghDYXNCCGEEEIaSss1I7Ztw1pg+WXLsmDbtseOCCGEEG9ptlq4qGbkueeew5o1a7By5Ups2bIFp0+fXvC9zz//PO6++2585CMfwUc+8hF85jOf+cD3q8S2bSSTScTj8Tk/dWuaJuLxOJLJJBsSQgghNy3NWAtrbkZefPFF7N27F/v378cbb7yBDRs24P7770c+n5/3/ZlMBr29vRgcHMSpU6cQjUZx33334dKlS67N14rjOMjn8zAMA4lEApeGhwEAl4aHkUgkYBgG8vk8HMdxvzPLAgYHi9t6okpXpfboaOW2nuiYh0rPqrJmzpUwZ/W6gJ7nDg3m0NNaWC1SI5s3b5bdu3eXn09NTUl3d7c888wzVX3+vffek2AwKD/84Q+r3qdt2wJAbNuu1e4ccrmcxGIxASDr/X5Jp9Oy3u8XABKLxSSXy7neh6RSIm1tIkBxm0q511Spq1I7lZJCICDpdFoKgYAenlVqK/asJGvmPEebOSvWLWlrd+7QaA49qYVSff2uaaG8QqGA9vZ2vPTSS3jwwQfL47t27cLY2BhOnDjxoRqO4yAcDuMf//Ef8fnPf37e91y/fh3Xr18vPx8fH0c0GsXVq1ddLZQ3w6XhYWy/5x5c8ftx6NAhPPLII+h89130v/IKejZtcic+OgrcdRcwPf3+mM8H/PKX7tZXUKWrUrukO7liBX5+6BA++8gjWFYoNLdnldoeeK571sx5Xm3mrFD3Bm2tzh0azqHSWlhifHwct91224culFdTMzI6Ooqenh6cPHkS27ZtK4/v27cP2WwWQ0NDH6rx2GOP4ac//Sl+9atfYeXKlfO+56mnnsLTTz89Z/zo0aNob2+v1i4hhBBCGsjExAR27NjRXKv2Pvvsszh27BgymcyCjQgAPPHEE9i7d2/5+cw3I/fddx+/GdGo69byfzcqtXX8nyRznlebOSvUvUFbq3OHhnPo1TcjVVHLtZ/r16+Lz+eT48ePV4w//PDD8sADD3zgZ//mb/5GQqGQvP7667XsUkQ0vWfE5yte2/P56ns9UoWuSu1USgqrVhWv+65apYdnldqKPSvJmjnP0WbOinVL2tqdOzSaw2a7Z2RRN7Du2bOn/Hxqakp6eno+8AbW//gf/6N0dHTIqVOnat2diNSvGTFNsxx+LBaTd06dknQ6Le+cOlUxbpqmq/2UdiYyOFjc1hNVugq1CxcuFE8oFy7UVVdEtMxDpWdlWTPnCpizB7qi6blDgzn0shYqa0aOHTsmK1askBdeeEHOnj0rjz76qNxyyy1y+fJlERHp6+uTxx9/vPz+Z599VpYvXy4vvfSS/PM//3P54ThO3f8xH8bY2Jhs3bq13PUVCoXigV4olLvErVu3ytjYmKv9kEpuzJmohVl7A3P2BuasBi9rYbX1u+Z7Rh566CH85je/wZNPPonLly9j48aNGBgYQGdnJwAgl8uhre39ny/5wQ9+gEKhgH/zb/5Nhc7+/fvx1FNP1bp7V4RCIQwMDMBxHEQiEUxOTpZfi0ajyGazCAaDCIVCnvoihBBCvKIZa+GibmDds2cP9uzZM+9rmUym4vmFCxcWswtlhEKhBQOORCIeuyGEEEK8p9lqYcutTUMIIYSQ5oLNCCGEEEIaCpsRQgghhDQUNiOEEEIIaShsRgghhBDSUFquGbFtG9YCyy9blgXbtj12RAghhHhLs9XClmpGbNtGMplEPB6HaZoVr5mmiXg8jmQyyYaEEELITUsz1sKWakYcx0E+n4dhGEgkErg0PAyguFhQIpGAYRjI5/NwHMf9ziwLGBwsbuuJKl2V2qOjldt6omMeKj2rypo5V8Kc1esCep47NJhDT2thtbj+rVcP0HKhvLa24oJGbW31XYRJha5K7VRKCoFA8aeGAwE9PKvUVuxZSdbMeY42c1asW9LW7tyh0Rxqv1BeI6hnMyIikhsakhgg/tIE+P1+iQGSGxpyL26a7x8wMw+fz/2iSap0PfBcKOVc8Pub37NKbQ881z1r5jyvNnP2xrNW5w4N51BpLSxRbf1uqcs0M0SvXcPhWWOHAUQnJtyLnzsHTE9Xjk1NAefPN6euSm0dPavUpmdvtOnZG20dPavU1tCz0lpYIy3ZjJiBAPpmjfUBMNvb3YuvWwe0zYrV5wPWrm1OXZXaOnpWqU3P3mjTszfaOnpWqa2hZ6W1sEZarhkxTROJ3l4YANaUxtYAMAAkenvn3FlcM5EIcPBg8UABitsDB4rjzairUltHzyq16dkbbXr2RltHzyq1NfOsvBbWyBIREU/3uAjGx8cRCoVg2zY6OjoWrWNZFuLxOAzDQCwWw8tHjuDMlSvY0NmJe3fuLI9ns1n3qxZaVvErtLVr63MwqtZVqD158SL6R0awfeNGLFu9um66ALTMQ6VnZVkz5wqYswe60PTcocEcelkLq63fS13tRTOCwSDC4TAAIJPJoKurC2f6+9GzaRMymQwSiQTC4TCCwaD7nUUi9T8QVeqq1O7uBkZGitt6o2MeKj2rypo5V8Kc1esCep47NJhDT2thlbRUMxIKhTAwMADHcRCJRDA5OVl+LRqNIpvNIhgMIhQKNdAlIYQQoo5mrIUt1YwAxUlYKGDXl2YIIYQQDWi2WthyN7ASQgghpLlgM0IIIYSQhsJmhBBCCCENhc0IIYQQQhoKmxFCCCGENJSWa0Zs24a1wPLLlmXBtm2PHRFCCCHe0my1sKWaEdu2kUwmEY/H5/zUrWmaiMfjSCaTbEgIIYTctDRjLWypZsRxHOTzeRiGgUQigUvDwwCAS8PDSCQSMAwD+XwejuO435llAYODxW09UaWrUnt0tHJbT3TMQ6VnVVkz50qYs3pdQM9zhwZz6GktrBbRANu2BYDYtu1aK5fLSSwWEwCy3u+XdDot6/1+ASCxWExyuZx7w6mUSFubCFDcplLuNVXqqtROpaQQCEg6nZZCIKCHZ5Xaij0ryZo5z9Fmzop1S9ranTs0mkNPaqFUX79brhkREckNDUkMEH9pAvx+v8QAyQ0NuRc3zfcPmJmHz1ccb0ZdDzwXSjkX/P7m96xS2wPPdc+aOc+rzZy98azVuUPDOVRaC0tUW79b6jLNDNFr13B41thhANGJCffi584B09OVY1NTxZUWm1FXpbaOnlVq07M32vTsjbaOnlVqa+hZaS2skZZsRsxAAH2zxvoAmO3t7sXXrQPaZsXq8xWXfG5GXZXaOnpWqU3P3mjTszfaOnpWqa2hZ6W1sEZarhkxTROJ3l4YANaUxtYAMAAkenvn3FlcM5EIcPBg8UABitsDB9wv+6xKV6W2jp5VatOzN9r07I22jp5VamvmWXktrJElIiKe7nERjI+PIxQKwbZtdHR0LFrHsizE43EYhoFYLIaXjxzBmStXsKGzE/fu3Fkez2az7lcttKziV2hr19bnYFStq1B78uJF9I+MYPvGjVi2enXddAFomYdKz8qyZs4VMGcPdKHpuUODOfSyFlZbv5e62otmBINBhMNhAEAmk0FXVxfO9PejZ9MmZDIZJBIJhMNhBINB9zuLROp/IKrUVand3Q2MjBS39UbHPFR6VpU1c66EOavXBfQ8d2gwh57WwippqWYkFAphYGAAjuMgEolgcnKy/Fo0GkU2m0UwGEQoFGqgS0IIIUQdzVgLW6oZAYqTsFDAri/NEEIIIRrQbLWw5W5gJYQQQkhzwWaEEEIIIQ2FzQghhBBCGgqbEUIIIYQ0FDYjhBBCCGkoLdeM2LYNa4Hlly3Lgm3bHjsihBBCvKXZamFLNSO2bSOZTCIej8/5qVvTNBGPx5FMJtmQEEIIuWlpxlrYUs2I4zjI5/MwDAOJRAKXhocBAJeGh5FIJGAYBvL5PBzHcb8zywIGB4vbeqJKV6X26Gjltp7omIdKz6qyZs6VMGf1uoCe5w4N5tDTWlgtogG2bQsAsW3btVYul5NYLCYAZL3fL+l0Wtb7/QJAYrGY5HI594ZTKZG2NhGguE2l3Guq1FWpnUpJIRCQdDothUBAD88qtRV7VpI1c56jzZwV65a0tTt3aDSHntRCqb5+t1wzIiKSGxqSGCD+0gT4/X6JAZIbGnIvbprvHzAzD5+vON6Muh54LpRyLvj9ze9ZpbYHnuueNXOeV5s5e+NZq3OHhnOotBaWqLZ+t9Rlmhmi167h8KyxwwCiExPuxc+dA6anK8empoorLTajrkptHT2r1KZnb7Tp2RttHT2r1NbQs9JaWCMt2YyYgQD6Zo31ATDb292Lr1sHtM2K1ecrLvncjLoqtXX0rFKbnr3RpmdvtHX0rFJbQ89Ka2GNtFwzYpomEr29MACsKY2tAWAASPT2zrmzuGYiEeDgweKBAhS3Bw64X/ZZla5KbR09q9SmZ2+06dkbbR09q9TWzLPyWlgjS0REPN3jIhgfH0coFIJt2+jo6Fi0jmVZiMfjMAwDsVgMLx85gjNXrmBDZyfu3bmzPJ7NZt2vWmhZxa/Q1q6tz8GoWleh9uTFi+gfGcH2jRuxbPXquukC0DIPlZ6VZc2cK2DOHuhC03OHBnPoZS2stn4vdbUXzQgGgwiHwwCATCaDrq4unOnvR8+mTchkMkgkEgiHwwgGg+53FonU/0BUqatSu7sbGBkpbuuNjnmo9Kwqa+ZcCXNWrwvoee7QYA49rYVV0lLNSCgUwsDAABzHQSQSweTkZPm1aDSKbDaLYDCIUCjUQJeEEEKIOpqxFrZUMwIUJ2GhgF1fmiGEEEI0oNlqYcvdwEoIIYSQ5oLNCCGEEEIaCpsRQgghhDQUNiOEEEIIaSgt14zkcjm8/vrr8772+uuvI5fLLVrbtm1YC6ymaFnWopdjVqWrUltHzyq16dkbbXr2RltHzyq1dfSsshYuisUsfPNf/st/kdWrV8uKFStk8+bNMvQhi+r86Ec/ko9+9KOyYsUK+fjHPy4/+clPatpfvRbKu3jxoqxatUqWLl0qr732mhQKheIiTIWCvPbaa7J06VJZtWqVXLx4sWbtsbEx2bp167yrHc6sjrh161YZGxtrCl0vPd+Yc7N6Vqntped6Zc2cP1ibOeuVs0rfOs6hylo4G2Wr9h47dkyWL18uhw4dkl/96lfyJ3/yJ3LLLbfIlStX5n3/q6++Kj6fT7797W/L2bNn5Rvf+IYsW7ZM3nrrrar3Wa9m5PTp07J06VIBIEuXLpWTJ05IOp2WkydOVIyfPn26Zm3TNMvLMcdiseKqh6+8UlwV8YZxs8ZVFlXpeun5nVOnJJ1OyzunTjWtZ5XaXnquV9bM+YO1mbNeOav0reMcqqyFs1HWjGzevFl2795dfj41NSXd3d3yzDPPzPv+L33pS/K5z32uYmzLli3yp3/6p1Xvs17NiIiUuz4AEiwtmxz0+8vhv/baa4vWnulUAUgMkFdL2/KBNKuzbbSuV57Xl3JeX8q5WT2r1PbKcz2zZs4LazNn/XJW6VvHOVRZC29ESTNy/fp18fl8cvz48Yrxhx9+WB544IF5PxONRuU//af/VDH25JNPyic+8YkF9/Pb3/5WbNsuP0zTFABy9epVKRQKrh8nT5yQoN8vt956q6TTabn11lsl6PfLyRMnXGu/c+qUrPf7xX/DY73fL++cOtWUul54vjHnZves+xzWO2vmzJxvppx1z6Pe2ipr4czj6tWrVTUjNS2UNzo6ip6eHpw8eRLbtm0rj+/btw/ZbBZDQ0NzPrN8+XL88Ic/RG9vb3ns+9//Pp5++mlcuXJl3v089dRTePrpp+eMHz16FO0NWNqYEEIIIbUzMTGBHTt26LlQ3hNPPIG9e/eWn4+PjyMajeK+++5ztWrvDP/ff//vuP8P/xDL/H4cOnQIjzzyCCbffRc/PXYMn/zX/9qV9qXhYWy/5x5cuGFsDYD+V15Bz6ZNTaerUntG98oNOXe++25Te1ap7YXnemfNnOfXZs565nyj9oUbxtZAjzzqra2yFs4wPj5e3RuruDpTxqvLNLPhPSP6XY/U8bqvSm0dr7Ez54W1mbN+Oav0reMcan3PiEjxBtY9e/aUn09NTUlPT88H3sD6+c9/vmJs27ZtDbmB1fO/phkcVHOndh10vfSs/I54zfJQ6VnpX3kwZ+bssWcdzh06zuFN8dc0x44dkxUrVsgLL7wgZ8+elUcffVRuueUWuXz5soiI9PX1yeOPP15+/6uvvipLly6V73znO/L222/L/v37G/anvfydkcZ4vjHnZvWsUlvH32Vgzh+szZz1ylmlbx3n8Kb4nRERke9973vyu7/7u7J8+XLZvHlzxdc58Xhcdu3aVfH+H/3oR3LnnXfK8uXL5a677mrYj56JFCdhptu7cQJEit2im/DHxsYW7FBN01zUwahSV6X2jbqzc25Wzyq1vfJcz6yZ88LazFm/nGdrz0aHPOqprbIW3ki19bumv6ZpFOPj4wiFQh96N26tTE5Oor+/H9u3b8eyZcvqpksqYc7eway9gTl7A3P2BpU5V1u/W25tGkIIIYQ0F2xGCCGEENJQ2IwQQgghpKGwGSGEEEJIQ2EzQgghhJCGwmaEEEIIIQ2FzQghhBBCGgqbEUIIIYQ0FDYjhBBCCGkoSxttoBpmfiS26qWIq2RychITExMYHx/nr/sphDl7B7P2BubsDczZG1TmPFO3P+zH3rVoRhzHAQBEo9EGOyGEEEJIrTiOg1AotODrWqxNMz09jdHRUQSDQSxZsqRuuuPj44hGozBNs65r3pBKmLN3MGtvYM7ewJy9QWXOIgLHcdDd3Y22toXvDNHim5G2tjZEIhFl+h0dHTzQPYA5ewez9gbm7A3M2RtU5fxB34jMwBtYCSGEENJQ2IwQQgghpKG0dDOyYsUK7N+/HytWrGi0lZsa5uwdzNobmLM3MGdvaIactbiBlRBCCCE3Ly39zQghhBBCGg+bEUIIIYQ0FDYjhBBCCGkobEYIIYQQ0lBu+mbkueeew5o1a7By5Ups2bIFp0+f/sD3/+M//iM+9rGPYeXKlfi93/s99Pf3e+RUb2rJ+fnnn8fdd9+Nj3zkI/jIRz6Cz3zmMx86L+R9aj2mZzh27BiWLFmCBx98UK3Bm4Racx4bG8Pu3btx++23Y8WKFbjzzjt5/qiCWnP+7ne/i49+9KPw+/2IRqP4yle+gt/+9rceudWT//k//ye+8IUvoLu7G0uWLEE6nf7Qz2QyGfz+7/8+VqxYgbVr1+KFF15Qa1JuYo4dOybLly+XQ4cOya9+9Sv5kz/5E7nlllvkypUr877/1VdfFZ/PJ9/+9rfl7Nmz8o1vfEOWLVsmb731lsfO9aLWnHfs2CHPPfecvPnmm/L222/LH/3RH0koFBLLsjx2rh+1Zj3DO++8Iz09PXL33XfLF7/4RW/MakytOV+/fl0++clPyvbt2+UXv/iFvPPOO5LJZGRkZMRj53pRa85HjhyRFStWyJEjR+Sdd96Rn/70p3L77bfLV77yFY+d60V/f798/etflx//+McCQI4fP/6B7zcMQ9rb22Xv3r1y9uxZ+d73vic+n08GBgaUebypm5HNmzfL7t27y8+npqaku7tbnnnmmXnf/6UvfUk+97nPVYxt2bJF/vRP/1SpT92pNefZvPfeexIMBuWHP/yhKos3DYvJ+r333pNPfepTkkqlZNeuXWxGqqDWnH/wgx9ILBaTQqHglcWbglpz3r17t9xzzz0VY3v37pVPf/rTSn3eTFTTjOzbt0/uuuuuirGHHnpI7r//fmW+btrLNIVCAcPDw/jMZz5THmtra8NnPvMZnDp1at7PnDp1quL9AHD//fcv+H6yuJxnMzExgcnJSdx6662qbN4ULDbr//Af/gPC4TD+7b/9t17Y1J7F5Pzf/tt/w7Zt27B79250dnbi4x//OL71rW9hamrKK9vasZicP/WpT2F4eLh8KccwDPT392P79u2eeG4VGlELtVgobzFcvXoVU1NT6OzsrBjv7OzE//7f/3vez1y+fHne91++fFmZT91ZTM6z+drXvobu7u45Bz+pZDFZ/+IXv8Df/d3fYWRkxAOHNweLydkwDLzyyivYuXMn+vv7cf78eTz22GOYnJzE/v37vbCtHYvJeceOHbh69Sr+4A/+ACKC9957D3/2Z3+Gv/iLv/DCcsuwUC0cHx/Hu+++C7/fX/d93rTfjBA9ePbZZ3Hs2DEcP34cK1eubLSdmwrHcdDX14fnn38et912W6Pt3NRMT08jHA7j4MGD2LRpEx566CF8/etfx3/9r/+10dZuKjKZDL71rW/h+9//Pt544w38+Mc/xk9+8hN885vfbLQ14pKb9puR2267DT6fD1euXKkYv3LlCrq6uub9TFdXV03vJ4vLeYbvfOc7ePbZZ/E//sf/wCc+8QmVNm8Kas36//7f/4sLFy7gC1/4QnlsenoaALB06VL8+te/xh133KHWtIYs5pi+/fbbsWzZMvh8vvLY+vXrcfnyZRQKBSxfvlypZx1ZTM5/+Zd/ib6+Pvy7f/fvAAC/93u/h2vXruHRRx/F17/+dbS18f/X9WChWtjR0aHkWxHgJv5mZPny5di0aRNefvnl8tj09DRefvllbNu2bd7PbNu2reL9APDzn/98wfeTxeUMAN/+9rfxzW9+EwMDA/jkJz/phVXtqTXrj33sY3jrrbcwMjJSfjzwwAP4V//qX2FkZATRaNRL+9qwmGP605/+NM6fP19u9gDg//yf/4Pbb7+djcgCLCbniYmJOQ3HTAMoXGatbjSkFiq7NbYJOHbsmKxYsUJeeOEFOXv2rDz66KNyyy23yOXLl0VEpK+vTx5//PHy+1999VVZunSpfOc735G3335b9u/fzz/trYJac3722Wdl+fLl8tJLL8k///M/lx+O4zTqn6ANtWY9G/41TXXUmnMul5NgMCh79uyRX//61/JP//RPEg6H5a/+6q8a9U/Qglpz3r9/vwSDQfmHf/gHMQxDfvazn8kdd9whX/rSlxr1T9ACx3HkzTfflDfffFMAyN/+7d/Km2++KRcvXhQRkccff1z6+vrK75/5094///M/l7fffluee+45/mmvW773ve/J7/7u78ry5ctl8+bN8tprr5Vfi8fjsmvXror3/+hHP5I777xTli9fLnfddZf85Cc/8dixntSS8+rVqwXAnMf+/fu9N64htR7TN8JmpHpqzfnkyZOyZcsWWbFihcRiMfnrv/5ree+99zx2rR+15Dw5OSlPPfWU3HHHHbJy5UqJRqPy2GOPyf/7f//Pe+MaMTg4OO85dybbXbt2STwen/OZjRs3yvLlyyUWi8nf//3fK/W4RITfbRFCCCGkcdy094wQQgghRA/YjBBCCCGkobAZIYQQQkhDYTNCCCGEkIbCZoQQQgghDYXNCCGEEEIaCpsRQgghhDQUNiOEEEIIaShsRgghhBDSUNiMEEIIIaShsBkhhBBCSENhM0IIIYSQhvL/A/vwHKusV1D2AAAAAElFTkSuQmCC",
      "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": 35,
   "id": "d1368356",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Testing input/output for model\n",
    "test = PINN(\n",
    "    N_input=2,\n",
    "    N_hidden_arr=[32,16,16,32],\n",
    "    N_output=2\n",
    ")\n",
    " # output and input are the same size/dimensions, means things are mappable???\n",
    "\n",
    "\n",
    "# Plan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e52e8c17",
   "metadata": {},
   "outputs": [],
   "source": [
    "## create neural network\n",
    "\n",
    "class Net:\n",
    "    def __init__(self):\n",
    "        # Create PINN, send off to device if CUDA is available\n",
    "        device = torch.device(\"cpu\")\n",
    "        print(\"Using CPU\")\n",
    "        self.model = PINN(\n",
    "            N_input=2,\n",
    "            N_hidden_arr=[32,16,16,32],\n",
    "            N_output = 2\n",
    "        ).to(device)\n",
    "\n",
    "        # Create data\n",
    "        self.d = 0.05\n",
    "        x = torch.arange(0, 1+d, d)\n",
    "        y = torch.arange(0, 1+d, d)\n",
    "\n",
    "        self.X = torch.stack(torch.meshgrid(x,y)).reshape(2,-1).T\n",
    "\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",
    "        self.X_train = torch.cat([b_left, b_bottom, b_right, b_top])\n",
    "\n",
    "        # Get wall velocities in form of [u,v]\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",
    "\n",
    "        self.data_train = torch.cat([walls_Y, walls_Y, walls_Y, wall_lid])\n",
    "\n",
    "        # Send data off to devices\n",
    "        self.X = self.X.to(device)\n",
    "        self.X_train = self.X_train.to(device)\n",
    "        self.data_train = self.data_train.to(device)\n",
    "        self.X.requires_grad = True\n",
    "        self.X_train.requires_grad = True\n",
    "\n",
    "\n",
    "        self.criterion = nn.MSELoss()\n",
    "        self.iter = 1\n",
    "\n",
    "        # Set up optimiser for LBFGS\n",
    "        self.optimiser = torch.optim.LBFGS(\n",
    "            params=self.model.parameters(),\n",
    "            lr=1.0,\n",
    "            max_iter = 75*10**3,\n",
    "            max_eval = 75*10**3,\n",
    "            history_size=50,\n",
    "            tolerance_change=1e-15,\n",
    "            tolerance_grad=1e-15,\n",
    "            line_search_fn=\"strong_wolfe\"\n",
    "        )\n",
    "\n",
    "        self.adam = torch.optim.Adam(self.model.parameters())\n",
    "    \n",
    "    def loss_function(self):\n",
    "        self.optimiser.zero_grad()\n",
    "        self.adam.zero_grad()\n",
    "\n",
    "        def get_gradient(output, input):\n",
    "            return torch.autograd.grad(outputs=output, \n",
    "                                       inputs=input,\n",
    "                                       grad_outputs=torch.ones_like(output),\n",
    "                                       create_graph=True,\n",
    "                                       retain_graph=True)[0]\n",
    "\n",
    "        # Get data loss\n",
    "        output = self.model(self.X_train)\n",
    "        psi_pred, p_pred = output[:,0], output[:,1]\n",
    "        psipred_dX = torch.autograd.grad(  # get velocities from psi \n",
    "            outputs=psi_pred,\n",
    "            inputs = self.X_train,\n",
    "            grad_outputs = torch.ones_like(psi_pred),\n",
    "            retain_graph=True, \n",
    "            create_graph=True,\n",
    "        )[0]\n",
    "\n",
    "        u_pred = psipred_dX[:,1]  # get each component from grad\n",
    "        v_pred = -psipred_dX[:,0]\n",
    "        vel_pred = torch.column_stack([u_pred, v_pred])\n",
    "\n",
    "        data_loss = self.criterion(vel_pred, self.data_train)  # compute loss from known values\n",
    "\n",
    "        # Get physics loss\n",
    "        output_phys = self.model(self.X)\n",
    "        psi, p = output_phys[:,0], output_phys[:,1]\n",
    "        psi_dX, p_dX = map(lambda x: get_gradient(x, self.X), [psi, p])\n",
    "        self.u, self.v = psi_dX[:,1], -1*psi_dX[:,0]\n",
    "        \n",
    "        p_x, p_y = p_dX[:,0], p_dX[:,1]\n",
    "\n",
    "        u_dX, v_dX = map(lambda x: get_gradient(x, self.X), [self.u,self.v])\n",
    "        u_dXX, v_dXX = map(lambda x: get_gradient(x, self.X), [u_dX, v_dX])\n",
    "\n",
    "        u_x, u_y = u_dX[:,0], u_dX[:,1]\n",
    "        v_x, v_y = v_dX[:,0], v_dX[:,1]\n",
    "\n",
    "        u_xx, u_yy = u_dXX[:,0], u_dXX[:,1]\n",
    "        v_xx, v_yy = v_dXX[:,0], v_dXX[:,1]\n",
    "\n",
    "        continuity_loss = self.criterion(u_x+v_y, torch.zeros(len(u_x)))\n",
    "        x_mom_loss = self.criterion(self.u*u_x+self.v*u_y+p_x/rho0, nu0*(u_xx+u_yy))\n",
    "        y_mom_loss = self.criterion(self.u*v_x+self.v*v_y+p_y/rho0, nu0*(v_xx+v_yy))\n",
    "\n",
    "        phys_loss = continuity_loss+x_mom_loss+y_mom_loss\n",
    "\n",
    "        # Put more weight on data loss?\n",
    "        # Change scheme as data goes on\n",
    "        if self.iter >= 1000:\n",
    "            loss = data_loss\n",
    "        elif self.iter > 1000 and self.iter < 2500:\n",
    "            loss = 0.7*data_loss+0.3*phys_loss\n",
    "        else:\n",
    "            loss = data_loss+phys_loss\n",
    "        loss.backward()\n",
    "\n",
    "        if self.iter % 100 == 0:\n",
    "            print(f\"Iteration {self.iter:5}, Total Loss {loss:.9f}\")\n",
    "\n",
    "            print((vel_pred).size())\n",
    "        self.iter+= 1\n",
    "\n",
    "        return loss\n",
    "        \n",
    "    def train(self):\n",
    "        self.model.train()\n",
    "        for i in range(1000):\n",
    "            self.adam.step(self.loss_function)\n",
    "        self.optimiser.step(self.loss_function)\n",
    "\n",
    "    def eval(self):\n",
    "        self.model.eval()\n",
    "\n",
    "        \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "id": "38c6e9e1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using CPU\n",
      "Iteration   100, Total Loss 0.026230508\n",
      "torch.Size([84, 2])\n",
      "Iteration   200, Total Loss 0.014797823\n",
      "torch.Size([84, 2])\n",
      "Iteration   300, Total Loss 0.013963315\n",
      "torch.Size([84, 2])\n",
      "Iteration   400, Total Loss 0.013744197\n",
      "torch.Size([84, 2])\n",
      "Iteration   500, Total Loss 0.013631235\n",
      "torch.Size([84, 2])\n",
      "Iteration   600, Total Loss 0.013560196\n",
      "torch.Size([84, 2])\n",
      "Iteration   700, Total Loss 0.013503134\n",
      "torch.Size([84, 2])\n",
      "Iteration   800, Total Loss 0.013437574\n",
      "torch.Size([84, 2])\n",
      "Iteration   900, Total Loss 0.013379758\n",
      "torch.Size([84, 2])\n",
      "Iteration  1000, Total Loss 0.012725575\n",
      "torch.Size([84, 2])\n",
      "Iteration  1100, Total Loss 0.008411493\n",
      "torch.Size([84, 2])\n"
     ]
    }
   ],
   "source": [
    "test = Net()\n",
    "test.train()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78e07c9a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pytorch",
   "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
}