From 1e074ef720af8c5775e1dfa3b84e7d00463093c2 Mon Sep 17 00:00:00 2001 From: Sichao He <1310722434@qq.com> Date: Wed, 18 Dec 2024 20:35:36 +0800 Subject: [PATCH] Update Navier_Stokes_inverse.ipynb --- .../Navier_Stokes_inverse.ipynb | 395 +++++++++++++++++- 1 file changed, 392 insertions(+), 3 deletions(-) diff --git a/docs/unit-examples-inverse/Navier_Stokes_inverse.ipynb b/docs/unit-examples-inverse/Navier_Stokes_inverse.ipynb index 076bd79..92d461f 100644 --- a/docs/unit-examples-inverse/Navier_Stokes_inverse.ipynb +++ b/docs/unit-examples-inverse/Navier_Stokes_inverse.ipynb @@ -173,13 +173,18 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 22, "metadata": {}, "outputs": [], "source": [ + "import re\n", + "\n", "import brainstate as bst\n", "import brainunit as u\n", "import numpy as np\n", + "from scipy.io import loadmat\n", + "import jax\n", + "import matplotlib.pyplot as plt\n", "\n", "import pinnx\n" ] @@ -193,7 +198,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ @@ -202,7 +207,391 @@ "unit_of_t = u.second\n", "unit_of_rho = u.kilogram / u.meter**3\n", "unit_of_c1 = u.UNITLESS\n", - "unit_of_c2 = u.meter2 / u.second" + "unit_of_c2 = u.meter2 / u.second\n", + "unit_of_u = u.meter / u.second\n", + "unit_of_v = u.meter / u.second\n", + "unit_of_p = u.pascal" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define Navier Stokes Equations:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "def Navier_Stokes_Equation(x, y):\n", + " jacobian = net.jacobian(x)\n", + " x_hessian = net.hessian(x, y=['u', 'v'], xi=['x'], xj=['x'])\n", + " y_hessian = net.hessian(x, y=['u', 'v'], xi=['y'], xj=['y'])\n", + "\n", + " u = y['u']\n", + " v = y['v']\n", + " p = y['p']\n", + " du_x = jacobian['u']['x']\n", + " du_y = jacobian['u']['y']\n", + " du_t = jacobian['u']['t']\n", + " dv_x = jacobian['v']['x']\n", + " dv_y = jacobian['v']['y']\n", + " dv_t = jacobian['v']['t']\n", + " dp_x = jacobian['p']['x']\n", + " dp_y = jacobian['p']['y']\n", + " du_xx = x_hessian['u']['x']['x']\n", + " du_yy = y_hessian['u']['y']['y']\n", + " dv_xx = x_hessian['v']['x']['x']\n", + " dv_yy = y_hessian['v']['y']['y']\n", + " continuity = du_x + dv_y\n", + " x_momentum = unit_of_rho * du_t + unit_of_rho * C1.value * (u * du_x + v * du_y) + dp_x - unit_of_rho * C2.value * (du_xx + du_yy)\n", + " y_momentum = unit_of_rho * dv_t + unit_of_rho * C1.value * (u * dv_x + v * dv_y) + dp_y - unit_of_rho * C2.value * (dv_xx + dv_yy)\n", + " return [continuity, x_momentum, y_momentum]\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define the neural network model:" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "net = pinnx.nn.Model(\n", + " pinnx.nn.DictToArray(x=unit_of_x, y=unit_of_y, t=unit_of_t),\n", + " pinnx.nn.FNN([3] + [50] * 6 + [3], \"tanh\"),\n", + " pinnx.nn.ArrayToDict(u=unit_of_u, v=unit_of_v, p=unit_of_p),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define Spatio-temporal domain" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "# Rectangular\n", + "Lx_min, Lx_max = 1.0, 8.0\n", + "Ly_min, Ly_max = -2.0, 2.0\n", + "# Spatial domain: X × Y = [1, 8] × [−2, 2]\n", + "space_domain = pinnx.geometry.Rectangle([Lx_min, Ly_min], [Lx_max, Ly_max])\n", + "# Time domain: T = [0, 7]\n", + "time_domain = pinnx.geometry.TimeDomain(0, 7)\n", + "# Spatio-temporal domain\n", + "geomtime = pinnx.geometry.GeometryXTime(space_domain, time_domain)\n", + "geomtime = geomtime.to_dict_point(x=unit_of_x, y=unit_of_y, t=unit_of_t)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define the true values of the $\\lambda_1, \\lambda_2$ and the function to load training data." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "# true values\n", + "C1true = 1.0\n", + "C2true = 0.01 * unit_of_c2\n", + "\n", + "\n", + "# Load training data\n", + "def load_training_data(num):\n", + " data = loadmat(\"../dataset/cylinder_nektar_wake.mat\")\n", + " U_star = data[\"U_star\"] # N x 2 x T\n", + " P_star = data[\"p_star\"] # N x T\n", + " t_star = data[\"t\"] # T x 1\n", + " X_star = data[\"X_star\"] # N x 2\n", + " N = X_star.shape[0]\n", + " T = t_star.shape[0]\n", + " # Rearrange Problem\n", + " XX = np.tile(X_star[:, 0:1], (1, T)) # N x T\n", + " YY = np.tile(X_star[:, 1:2], (1, T)) # N x T\n", + " TT = np.tile(t_star, (1, N)).T # N x T\n", + " UU = U_star[:, 0, :] # N x T\n", + " VV = U_star[:, 1, :] # N x T\n", + " PP = P_star # N x T\n", + " x = XX.flatten()[:, None] # NT x 1\n", + " y = YY.flatten()[:, None] # NT x 1\n", + " t = TT.flatten()[:, None] # NT x 1\n", + " u = UU.flatten()[:, None] # NT x 1\n", + " v = VV.flatten()[:, None] # NT x 1\n", + " p = PP.flatten()[:, None] # NT x 1\n", + " # training domain: X × Y = [1, 8] × [−2, 2] and T = [0, 7]\n", + " data1 = np.concatenate([x, y, t, u, v, p], 1)\n", + " data2 = data1[:, :][data1[:, 2] <= 7]\n", + " data3 = data2[:, :][data2[:, 0] >= 1]\n", + " data4 = data3[:, :][data3[:, 0] <= 8]\n", + " data5 = data4[:, :][data4[:, 1] >= -2]\n", + " data_domain = data5[:, :][data5[:, 1] <= 2]\n", + " # choose number of training points: num =7000\n", + " idx = np.random.choice(data_domain.shape[0], num, replace=False)\n", + " x_train = data_domain[idx, 0]\n", + " y_train = data_domain[idx, 1]\n", + " t_train = data_domain[idx, 2]\n", + " u_train = data_domain[idx, 3]\n", + " v_train = data_domain[idx, 4]\n", + " p_train = data_domain[idx, 5]\n", + " # return [x_train, y_train, t_train, u_train, v_train, p_train]\n", + " return [x_train * unit_of_x, y_train * unit_of_y, t_train * unit_of_t,\n", + " u_train * unit_of_u, v_train * unit_of_v, p_train * unit_of_p]\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define the Parameters to be identified:" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "C1 = bst.ParamState(0.0)\n", + "C2 = bst.ParamState(0.0 * unit_of_c2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Get the training data:" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "[ob_x, ob_y, ob_t, ob_u, ob_v, ob_p] = load_training_data(num=7000)\n", + "ob_xyt = {\"x\": ob_x, \"y\": ob_y, \"t\": ob_t}\n", + "ob_yv = {\"u\": ob_u, \"v\": ob_v, }\n", + "observe_bc = pinnx.icbc.PointSetBC(ob_xyt, ob_yv)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define the problem and solve the inverse problem:" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "problem = pinnx.problem.TimePDE(\n", + " geomtime,\n", + " Navier_Stokes_Equation,\n", + " [observe_bc],\n", + " net,\n", + " num_domain=700,\n", + " num_boundary=200,\n", + " num_initial=100,\n", + " anchors=ob_xyt,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setup the model trainer:" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "model = pinnx.Trainer(problem, external_trainable_variables=[C1, C2])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define the callbacks for storing results:" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [], + "source": [ + "fnamevar = \"variables.dat\"\n", + "variable = pinnx.callbacks.VariableValue([C1, C2], period=100, filename=fnamevar)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Compile, train and save trainer:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model.compile(bst.optim.Adam(1e-3)).train(iterations=10000, callbacks=[variable],\n", + " display_every=1000, disregard_previous_best=True)\n", + "model.saveplot(issave=True, isplot=True)\n", + "model.compile(bst.optim.Adam(1e-4)).train(iterations=10000, callbacks=[variable],\n", + " display_every=1000, disregard_previous_best=True)\n", + "model.saveplot(issave=True, isplot=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Calculate the mean residual error:" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Mean residual: [0.04728882 * becquerel, 0.04445616 * kilogram / klitre * (meter / second) / second, 0.0641908 * kilogram / klitre * (meter / second) / second]\n" + ] + } + ], + "source": [ + "f = model.predict(ob_xyt, operator=Navier_Stokes_Equation)\n", + "print(\"Mean residual:\", jax.tree.map(lambda x: u.math.mean(u.math.abs(x)), f))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot Variables:" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/var/folders/x0/_jqxxbbn0rsdn6b4h6fxbrjr0000gn/T/ipykernel_24545/4284654551.py:6: DeprecationWarning: string or file could not be read to its end due to unmatched data; this will raise a ValueError in the future.\n", + " np.fromstring(\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAHHCAYAAABTMjf2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/GU6VOAAAACXBIWXMAAA9hAAAPYQGoP6dpAABDrUlEQVR4nO3dfVxUdd7/8feAMNwJSCDiTd5Shrel4mJbalqoXVaWrbnuhrprW2J3xpZZqdnPZbfa1kourb2xutzurCxXSzNT28y8jVXTLArTVMCbBAERZL6/P1xGJ26cwYEBzuv5eMxD5pwP53zmMBzenvOdc2zGGCMAAAAL8vN1AwAAAL5CEAIAAJZFEAIAAJZFEAIAAJZFEAIAAJZFEAIAAJZFEAIAAJZFEAIAAJZFEAIAAJZFEALQ5Kxdu1Y2m01r1671+HvHjx+vsLAwt2ptNptmzZrl8ToANBwEIQD14oYbblBISIhOnDhRbc24ceMUGBioo0eP1mNnAKyMIASgXowbN04nT57UkiVLqpxfXFys9957T8OGDdNFF110Qeu6+uqrdfLkSV199dUXtBwATR9BCEC9uOGGG9S8eXO9+uqrVc5/7733VFRUpHHjxtV6HSUlJXI4HPLz81NQUJD8/NjFAagZewkA9SI4OFg333yzVq9erby8vErzX331VTVv3lw///nPlZaWph49eigsLEzh4eEaPny4/vOf/7jUV4wDev311/Xoo4+qTZs2CgkJUUFBQZVjhP7973/r1ltv1cUXXyy73a527drp/vvv18mTJ6vs97vvvlNycrJCQ0PVunVrzZ49W8aY877OAwcOaOLEiYqNjZXdble3bt30j3/8o1Ld888/r27duikkJEQtWrRQ3759qw2JAOpOM183AMA6xo0bp5dffllvvvmmpkyZ4px+7NgxrVy5UmPHjtWhQ4f07rvv6tZbb1XHjh2Vm5urF154QQMHDtSuXbvUunVrl2U+8cQTCgwMVFpamk6dOqXAwMAq17148WIVFxfrrrvu0kUXXaRNmzbp+eef1w8//KDFixe71JaXl2vYsGH62c9+pieffFIrVqzQzJkzdfr0ac2ePbva15ebm6uf/exnstlsmjJlimJiYvTBBx/oN7/5jQoKCnTfffdJkv7617/qnnvu0ejRo3XvvfeqpKRE27dv18aNG/XLX/6yllsXQK0YAKgnp0+fNnFxcSYpKcll+oIFC4wks3LlSlNSUmLKy8td5mdnZxu73W5mz57tnLZmzRojyXTq1MkUFxe71FfMW7NmjXPaT2uMMSY9Pd3YbDbz/fffO6elpKQYSebuu+92TnM4HOb66683gYGB5vDhw87pkszMmTOdz3/zm9+YuLg4c+TIEZf13HbbbSYiIsLZw4033mi6detW3WYCUI84NQag3vj7++u2227Thg0btHfvXuf0V199VbGxsRoyZIjsdrtzbE95ebmOHj2qsLAwXXrppdq2bVulZaakpCg4OPi86z63pqioSEeOHNGAAQNkjNEXX3xRqf7cI1YVR3hKS0v10UcfVbl8Y4zefvttjRw5UsYYHTlyxPlITk5Wfn6+s//IyEj98MMP2rx583n7BlC3CEIA6lXFYOiK8TA//PCD/v3vf+u2226Tv7+/HA6H/vKXvyg+Pl52u13R0dGKiYnR9u3blZ+fX2l5HTt2dGu9+/bt0/jx4xUVFaWwsDDFxMRo4MCBklRpuX5+furUqZPLtEsuuUSSXALcuQ4fPqzjx4/rxRdfVExMjMtjwoQJkuQcG/XQQw8pLCxMiYmJio+PV2pqqtavX+/W6wDgXYwRAlCv+vTpo65du+q1117T9OnT9dprr8kY4wxIf/jDH/TYY49p4sSJeuKJJxQVFSU/Pz/dd999cjgclZbnztGg8vJyXXvttTp27Jgeeughde3aVaGhoTpw4IDGjx9f5XI9VbGMX/3qV0pJSamypmfPnpKkyy67THv27NGyZcu0YsUKvf322/rf//1fzZgxQ48//vgF9wLAfQQhAPVu3Lhxeuyxx7R9+3a9+uqrio+PV79+/SRJb731lgYPHqy///3vLt9z/PhxRUdH12p9O3bs0Ndff62XX35Zt99+u3P6qlWrqqx3OBz67rvvnEeBJOnrr7+WJHXo0KHK74mJiVHz5s1VXl6uoUOHnren0NBQjRkzRmPGjFFpaaluvvlmzZkzRw8//LCCgoI8eHUALgSnxgDUu4qjPzNmzFBmZqbLtYP8/f0rfUx98eLFOnDgQK3X5+/vL0kuyzXG6Nlnn632e+bNm+dSO2/ePAUEBGjIkCHVruOWW27R22+/rZ07d1aaf/jwYefXP71ydmBgoBISEmSMUVlZmXsvCoBXcEQIQL3r2LGjBgwYoPfee0+SXILQ//zP/2j27NmaMGGCBgwYoB07duif//xnpTE7nujatas6d+6stLQ0HThwQOHh4Xr77bf1448/VlkfFBSkFStWKCUlRf3799cHH3yg5cuXa/r06YqJial2PX/84x+1Zs0a9e/fX5MmTVJCQoKOHTumbdu26aOPPtKxY8ckSdddd51atWqlK6+8UrGxsdq9e7fmzZun66+/Xs2bN6/16wTgOY4IAfCJivCTmJioLl26OKdPnz5dDzzwgFauXKl7771X27Zt0/Lly9WuXbtarysgIED/+te/1Lt3b6Wnp+vxxx9XfHy8XnnllSrr/f39tWLFCuXk5Oj3v/+9Nm/erJkzZ+qJJ56ocT2xsbHatGmTJkyYoHfeeUdTpkzRs88+q2PHjulPf/qTs+53v/udCgsL9cwzzyg1NVXvvvuu7rnnHi1atKjWrxFA7djMT49BAwAAWARHhAAAgGURhAAAgGURhAAAgGURhAAAgGURhAAAgGURhAAAgGVxQcUaOBwOHTx4UM2bN5fNZvN1OwAAwA3GGJ04cUKtW7eWn1/Nx3wIQjU4ePDgBV3EDQAA+M7+/fvVtm3bGmsIQjWouNT9/v37FR4e7uNuAACAOwoKCtSuXTu3bllDEKpBxemw8PBwghAAAI2MO8NaGCwNAAAsiyAEAAAsiyAEAAAsiyAEAAAsiyAEAAAsiyAEAAAsiyAEAAAsiyAEAAAsq8kHoWXLlunSSy9VfHy8/va3v/m6HQAA0IA06StLnz59WlOnTtWaNWsUERGhPn36aNSoUbrooot83RoAAGgAmvQRoU2bNqlbt25q06aNwsLCNHz4cH344Ye+bgsAADQQDToIffLJJxo5cqRat24tm82md999t1JNRkaGOnTooKCgIPXv31+bNm1yzjt48KDatGnjfN6mTRsdOHCgPloHAACNQIM+NVZUVKRevXpp4sSJuvnmmyvNf+ONNzR16lQtWLBA/fv319y5c5WcnKw9e/aoZcuWPujYQ0VF1c/z95eCgtyr9fOTgoNrV1tcLBlTda3NJoWE1K725EnJ4ai+j9DQ2tWWlEjl5d6pDQk507cknTolnT7tndrg4DPbWZJKS6WyMu/UBgWdeV94WltWdqa+Ona71KyZ57WnT5/ZFtUJDJQCAjyvLS8/87OrTkDAmXpPax2OM+81b9Q2a3ZmW0hnfieKi71T68nvPfuIqmvZR3he2xD2Eb5kGglJZsmSJS7TEhMTTWpqqvN5eXm5ad26tUlPTzfGGLN+/Xpz0003Oeffe++95p///Ge16ygpKTH5+fnOx/79+40kk5+f790XU+HMLqPqx4gRrrUhIdXXDhzoWhsdXX1t376ute3bV1+bkOBam5BQfW379q61fftWXxsd7Vo7cGD1tSEhrrUjRtS83c41enTNtYWFZ2tTUmquzcs7Wzt5cs212dlna9PSaq7dufNs7cyZNddu2nS29skna65ds+Zs7bx5NdcuW3a2duHCmmvffPNs7Ztv1ly7cOHZ2mXLaq6dN+9s7Zo1Ndc++eTZ2k2baq6dOfNs7c6dNdempZ2tzc6uuXby5LO1eXk116aknK0tLKy5dvRo46KmWvYRZx7sI84+GvM+wsvy8/ONu3+/G/SpsZqUlpZq69atGjp0qHOan5+fhg4dqg0bNkiSEhMTtXPnTh04cECFhYX64IMPlJycXO0y09PTFRER4Xy0a9euzl8HAADwHZsxxvi6CXfYbDYtWbJEN910k6Sz438+++wzJSUlOesefPBBrVu3Ths3bpQkLV26VGlpaXI4HHrwwQd1xx13VLuOU6dO6dQ5h/EKCgrUrl075efnKzw83PsvisPentdy2Nvz2oZw2JtTY+7VcmrsLPYRntc25n2ElxUUFCgiIsKtv98NeoyQN9xwww264YYb3Kq12+2yV+yw6sO5v5C+qj13x+TN2nN3pN6sPXfH781au/3sHytv1gYGnv3j6qvagAD3dzae1DZrdnaH581af3/338Oe1Pr51U2tzVY3tVLDqGUfcQb7CM9rPfm996FGe2osOjpa/v7+ys3NdZmem5urVq1a+agrAADQmDTaIBQYGKg+ffpo9erVzmkOh0OrV692OVUGAABQnQZ9zKqwsFBZWVnO59nZ2crMzFRUVJQuvvhiTZ06VSkpKerbt68SExM1d+5cFRUVacKECT7sGgAANBYNOght2bJFgwcPdj6fOnWqJCklJUUvvfSSxowZo8OHD2vGjBnKyclR7969tWLFCsXGxl7QejMyMpSRkaHymgbRAQCARq/RfGrMFzwZdQ4AABoGT/5+N9oxQgAAABeKIAQAACyLIAQAACyLIAQAACyLIAQAACyLIAQAACyLIFSFjIwMJSQkqF+/fr5uBQAA1CGuI1QDriMEAEDjw3WEAAAA3EAQAgAAlkUQAgAAlkUQAgAAlkUQAgAAlkUQAgAAlkUQqgLXEQIAwBq4jlANuI4QAACND9cRAgAAcANBCAAAWBZBCAAAWBZBCAAAWBZBCAAAWBZBCAAAWBZBCAAAWBZBCAAAWBZBqApcWRoAAGvgytI14MrSAAA0PlxZGgAAwA0EIQAAYFkEIQAAYFkEIQAAYFkEIQAAYFkEIQAAYFkEIQAAYFkEIQAAYFkEIQAAYFkEoSpwiw0AAKyBW2zUgFtsAADQ+HCLDQAAADcQhAAAgGURhAAAgGURhAAAgGURhAAAgGURhAAAgGURhAAAgGURhAAAgGURhAAAgGURhAAAgGURhAAAgGURhAAAgGURhKrA3ecBALAG7j5fA+4+DwBA48Pd5wEAANxAEAIAAJZFEAIAAJZFEAIAAJZFEAIAAJZFEAIAAJZFEAIAAJZFEAIAAJZFEAIAAJZFEAIAAJZFEAIAAJZFEAIAAJZFEAIAAJZFEAIAAJZFEAIAAJZFEAIAAJZFEKpCRkaGEhIS1K9fP1+3AgAA6pDNGGN83URDVVBQoIiICOXn5ys8PNzX7QAAADd48vebI0IAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEIAAMCyCEJVyMjIUEJCgvr16+frVgAAQB2yGWOMr5toqAoKChQREaH8/HyFh4f7uh0AAOAGT/5+c0QIAABYFkEIAABYFkEIAABYFkEIAABYFkEIAABYFkEIAABYFkEIAABYFkEIAABYFkEIAABYFkEIAABYFkEIAABYFkEIAABYFkEIAABYFkEIAABYFkEIAABYFkEIAABYFkEIAABYFkEIAABYFkEIAABYFkEIAABYFkEIAABYVjNfNwAAQGNRXl6usrIyX7cBSQEBAfL397/g5RCEAAA4D2OMcnJydPz4cV+3gnNERkaqVatWstlstV4GQQgAgPOoCEEtW7ZUSEjIBf3hxYUzxqi4uFh5eXmSpLi4uFoviyAEAEANysvLnSHooosu8nU7+K/g4GBJUl5enlq2bFnr02QMlgYAoAYVY4JCQkJ83Al+quJnciHjtghCAAC4gdNhDY83fiYEIQAAYFkEIQAAUK29e/fKZrMpMzPT163UCYIQAABNWE5Oju6++2516tRJdrtd7dq108iRI7V69WpJ0osvvqhBgwYpPDxcNputzi4RMGjQIN133311suwLQRACAKCJ2rt3r/r06aOPP/5YTz31lHbs2KEVK1Zo8ODBSk1NlSQVFxdr2LBhmj59uo+79Q2CEAAATdTkyZNls9m0adMm3XLLLbrkkkvUrVs3TZ06VZ9//rkk6b777tO0adP0s5/9rMZlfffddxo8eLBCQkLUq1cvbdiwwTnv6NGjGjt2rNq0aaOQkBD16NFDr732mnP++PHjtW7dOj377LOy2Wyy2Wzau3dvnbxmT3EdIQAAPGWMVFxc/+sNCZHc/KTUsWPHtGLFCs2ZM0ehoaGV5kdGRnq06kceeURPP/204uPj9cgjj2js2LHKyspSs2bNVFJSoj59+uihhx5SeHi4li9frl//+tfq3LmzEhMT9eyzz+rrr79W9+7dNXv2bElSTEyMR+uvK14LQrm5uXrhhRc0Y8YMby0SAICGqbhYCgur//UWFkpVhJqqZGVlyRijrl27emXVaWlpuv766yVJjz/+uLp166asrCx17dpVbdq0UVpamrP27rvv1sqVK/Xmm28qMTFRERERCgwMVEhIiFq1auWVfrzFa6fGcnJy9Pjjj3trcQAA4AIYY7y6vJ49ezq/rrilRcUtLsrLy/XEE0+oR48eioqKUlhYmFauXKl9+/Z5tYe64PYRoe3bt9c4f8+ePRfcDAAAjUJIyJmjM75Yr5vi4+Nls9n01VdfeWXVAQEBzq8rLmTocDgkSU899ZSeffZZzZ07Vz169FBoaKjuu+8+lZaWemXddcntINS7d2/ZbLYqE2bF9IZ61c1Ro0Zp7dq1GjJkiN566y1ftwMAaOxsNrdPUflKVFSUkpOTlZGRoXvuuafSOKHjx497PE6oOuvXr9eNN96oX/3qV5LOBKSvv/5aCQkJzprAwECVl5d7ZX3e5PapsaioKP31r39VdnZ2pcd3332nZcuW1WWfF+Tee+/VK6+84us2AACoVxkZGSovL1diYqLefvttffPNN9q9e7eee+45JSUlSToztCUzM1NZWVmSpB07digzM1PHjh1zez3x8fFatWqVPvvsM+3evVu/+93vlJub61LToUMHbdy4UXv37tWRI0ecR5N8ze0g1KdPHx08eFDt27ev8tGmTRuvn4/0lkGDBql58+a+bgMAgHrVqVMnbdu2TYMHD9YDDzyg7t2769prr9Xq1as1f/58SdKCBQt0+eWXa9KkSZKkq6++WpdffrmWLl3q9noeffRRXXHFFUpOTtagQYPUqlUr3XTTTS41aWlp8vf3V0JCgmJiYhrM+CGbcTO9LFmyREVFRc7DXj/1448/aunSpUpJSfGogU8++URPPfWUtm7dqkOHDmnJkiWVNl5GRoaeeuop5eTkqFevXnr++eeVmJjo0XrWrl2refPmeXRqrKCgQBEREcrPz1d4eLhH6wMANA0lJSXKzs5Wx44dFRQU5Ot2cI7qfjae/P12e4zQqFGjapzfokULj0OQJBUVFalXr16aOHGibr755krz33jjDU2dOlULFixQ//79NXfuXCUnJ2vPnj1q2bKlpDPjl06fPl3pez/88EO1bt3a454AAIA1XNB1hNavX6++ffvKbrfXehnDhw/X8OHDq53/zDPPaNKkSZowYYKkM4fwli9frn/84x+aNm2aJHntRnCnTp3SqVOnnM8LCgq8slwAANAwXdB1hIYPH64DBw54q5dKSktLtXXrVg0dOtQ5zc/PT0OHDnW5tLe3pKenKyIiwvlo166d19cBAAAajgsKQnU9OPrIkSMqLy9XbGysy/TY2Fjl5OS4vZyhQ4fq1ltv1fvvv6+2bdtWG6Iefvhh5efnOx/79++/oP4BAEDDZol7jX300Udu1dnt9gs6zQcAABqXCzoi9MILL1Q6WuNN0dHR8vf3r3Qtgtzc3AZ3rxIAAND4XFAQ+uUvf1nlHW29JTAwUH369NHq1aud0xwOh1avXu28EBQAAEBteXxq7JtvvlF8fLzXGigsLHRezVKSsrOzlZmZqaioKF188cWaOnWqUlJS1LdvXyUmJmru3LkqKipyfooMAACgtjwKQlu3btUNN9zg1U+KbdmyRYMHD3Y+nzp1qiQpJSVFL730ksaMGaPDhw9rxowZysnJUe/evbVixYo6PSWXkZHhvCw5AABouty+svTHH3+sW265Rc8995x+/etf13VfDQJXlgYAWP3K0mvXrtXgwYP1448/eu0mrd7ijStLuzVG6J133tH//M//KD093TIhCACApiAnJ0d33323OnXqJLvdrnbt2mnkyJFavXq1jh07prvvvluXXnqpgoODdfHFF+uee+5Rfn6+1/vo0KGD5s6d6/XlXii3To2NGTNGs2bN0p133lnX/QAAAC/Zu3evrrzySkVGRuqpp55Sjx49VFZWppUrVyo1NVVvvfWWDh48qKeffloJCQn6/vvvdeedd+rgwYMe3ZuzMXPriFDnzp314Ycf6uTJk3XdDwAA8JLJkyfLZrNp06ZNuuWWW3TJJZeoW7dumjp1qj7//HN1795db7/9tkaOHKnOnTvrmmuu0Zw5c/Svf/2r0j08t27dqr59+yokJEQDBgzQnj17nPO+/fZb3XjjjYqNjVVYWJj69evncg2/QYMG6fvvv9f9998vm80mm81Wb9vgfNwKQp9++qmKi4s1atQolZWV1XVPAAA0aMYYFRUV1fvDkzs6HDt2TCtWrFBqamqVl7qpbrxPxbiaZs1cTxo98sgj+vOf/6wtW7aoWbNmmjhxonNeYWGhRowYodWrV+uLL77QsGHDNHLkSO3bt0/SmSE2bdu21ezZs3Xo0CEdOnTI7ddR19w6NRYdHa01a9bohhtu0C9+8QstWbKkrvsCAKDBKi4uVlhYWL2vt7Cw0O3r92VlZckYo65du7q9/CNHjuiJJ57QHXfcUWnenDlzNHDgQEnStGnTdP3116ukpERBQUHq1auXevXq5ax94okntGTJEi1dulRTpkxRVFSU/P391bx58wZ3QWS3L6gYFhamDz74QH5+F3QNxkYhIyNDCQkJ6tevn69bAQCgVjy9H2hBQYGuv/56JSQkaNasWZXm9+zZ0/l1XFycJCkvL0/SmYCWlpamyy67TJGRkQoLC9Pu3budR4QaMo+uI2S327V48eK66qXBSE1NVWpqqvPjdwAAnCskJESFhYU+Wa+74uPjZbPZ9NVXX5239sSJExo2bJiaN2+uJUuWKCAgoFLNudMqxvg4HA5JUlpamlatWqWnn35aXbp0UXBwsEaPHq3S0lK3+/UVj68sbYUjQgAA1MRms9XpLaa8ISoqSsnJycrIyNA999xTqd/jx48rMjJSBQUFSk5Olt1u19KlS2t1raT169dr/PjxGjVqlKQzR4j27t3rUhMYGNggL1TscarZtm2bduzY4Xz+3nvv6aabbtL06dMbRfIDAMAqKu6SkJiYqLffflvffPONdu/ereeee05JSUkqKCjQddddp6KiIv39739XQUGBcnJylJOT41FoiY+P1zvvvKPMzEz95z//0S9/+Uvn0aIKHTp00CeffKIDBw7oyJEj3n6pteZxEPrd736nr7/+WpL03Xff6bbbblNISIgWL16sBx980OsNAgCA2unUqZO2bdumwYMH64EHHlD37t117bXXavXq1Zo/f762bdumjRs3aseOHerSpYvi4uKcj/3797u9nmeeeUYtWrTQgAEDNHLkSCUnJ+uKK65wqZk9e7b27t2rzp07KyYmxtsvtdbcvsVGhYiICG3btk2dO3fWn/70J3388cdauXKl1q9fr9tuu82jDdfQcYsNAIDVb7HRkNXbLTbOZYxxHu766KOPNGLECElSu3btGtShLgAAgPPxOAj17dtX/+///T/93//9n9atW6frr79ekpSdnV2nd4QHAADwNo+D0Ny5c7Vt2zZNmTJFjzzyiLp06SJJeuuttzRgwACvNwgAAFBXPP74fM+ePV0+NVbhqaeekr+/v1ea8rWMjAznSHsAANB0ee2iQEFBQVVegKkxSk1N1a5du7R582ZftwIAAOoQV0cEAACWRRACAACWRRACAACWRRACAACW5VEQOnTokBYtWqT333+/0n3FioqKNHv2bK82BwAAUJfcDkKbN29WQkKCUlNTNXr0aHXr1k1ffvmlc35hYaEef/zxOmkSAAC4z2az1fiYNWtWvfaTlZWlCRMmqG3btrLb7erYsaPGjh2rLVu2OGvmzJmjAQMGKCQkRJGRkfXWm9tBaPr06Ro1apR+/PFH5ebm6tprr9XAgQP1xRdf1GV/AADAQ4cOHXI+5s6dq/DwcJdpaWlpzlpjjE6fPl1nvWzZskV9+vTR119/rRdeeEG7du3SkiVL1LVrVz3wwAPOutLSUt16662666676qyXKhk3tWjRwuzZs8dlWnp6umnRooXZtGmTycnJMX5+fu4urkGbN2+eueyyy8wll1xiJJn8/HxftwQA8JGTJ0+aXbt2mZMnT/q6lVpZuHChiYiIcD5fs2aNkWTef/99c8UVV5iAgACzZs0ak5KSYm688UaX77333nvNwIEDnc/Ly8vNH/7wB9OhQwcTFBRkevbsaRYvXlztuh0Oh+nWrZvp06ePKS8vrzT/xx9/PG+/NanuZ5Ofn+/232+PrixdUlLi8nzatGlq1qyZrrvuOv3jH//wXjrzsdTUVKWmpjrvXgsAQJWKiqqf5+8vnXu3+ppq/fyk4OCaa0NDPe+vBtOmTdPTTz+tTp06qUWLFm59T3p6uhYtWqQFCxYoPj5en3zyiX71q18pJiZGAwcOrFSfmZmpL7/8Uq+++qr8/CqfhKrPU2DVcTsIde/eXZ999pl69uzpMj0tLU0Oh0Njx471enMAADRoYWHVzxsxQlq+/Ozzli2l4uKqawcOlNauPfu8QwfpyBHXGmNq22WVZs+erWuvvdbt+lOnTukPf/iDPvroIyUlJUmSOnXqpE8//VQvvPBClUHom2++kSR17drVO03XAbeD0O23365169bpzjvvrDTvwQcflDFGCxYs8GpzAACgbvTt29ej+qysLBUXF1cKT6Wlpbr88sur/B7j5fBWF9wOQr/97W/129/+ttr5Dz30kB566CGvNAUAQKNQWFj9vJ/eiDwvr/ran5422ru31i25K/Qnp9r8/PwqBZeysjLn14X/fa3Lly9XmzZtXOrsdnuV67jkkkskSV999VW1YcnXPL77PAAA+C9Pxu3UVa2XxMTEaOfOnS7TMjMznTdUT0hIkN1u1759+6o8DVaV3r17KyEhQX/+8581ZsyYSuOEjh8/7vNxQh5fWfqzzz6riz4AAIAPXXPNNdqyZYteeeUVffPNN5o5c6ZLMGrevLnS0tJ0//336+WXX9a3336rbdu26fnnn9fLL79c5TJtNpsWLlyor7/+WldddZXef/99fffdd9q+fbvmzJmjG2+80Vm7b98+ZWZmat++fSovL1dmZqYyMzOdR6LqikdHhN5//31NmDBBubm5ddUPAADwgeTkZD322GN68MEHVVJSookTJ+r222/Xjh07nDVPPPGEYmJilJ6eru+++06RkZG64oorNH369GqXm5iYqC1btmjOnDmaNGmSjhw5ori4OA0YMEBz58511s2YMcMlUFWcSluzZo0GDRrk9ddbwWbcHMm0aNEiTZ48We+8846GDh1aZw01JBUfn8/Pz1d4eLiv2wEA+EBJSYmys7PVsWNHBZ37cXj4XHU/G0/+frt1amzu3Ln67W9/q0WLFlkmBAEAgKbPrVNjU6dO1XPPPacbbrihrvsBAACoN24dEbryyiv1v//7vzp69Ghd9wMAAFBv3ApCq1atUseOHXXttdeqoKCgrnvyuYyMDCUkJKhfv36+bgUAANQht4JQUFCQli5dqoSEBA0bNqyue/K51NRU7dq1S5s3b/Z1KwAAoA65fR0hf39/LVq0SImJiXXZDwAAQL3x+IKK537mHwAAoDHzOAgBAAA0FbW611hJSYm2b9+uvLw8ORwOl3l8xB4AADQWHgehFStW6Pbbb9eRI0cqzbPZbCovL/dKYwAAAHXN41Njd999t2699VYdOnRIDofD5UEIAgDA92w2W42PWbNm1Ws/WVlZmjBhgtq2bSu73a6OHTtq7Nix2rJliyRp7969+s1vfqOOHTsqODhYnTt31syZM1VaWlrnvXl8RCg3N1dTp05VbGxsXfQDAAAu0KFDh5xfv/HGG5oxY4b27NnjnBYWFub82hij8vJyNWtWq9Ey57VlyxYNGTJE3bt31wsvvKCuXbvqxIkTeu+99/TAAw9o3bp1+uqrr+RwOPTCCy+oS5cu2rlzpyZNmqSioiI9/fTTddJXBY9f9ejRo7V27Vp17ty5LvoBAKDRKCoqqnaev7+/y41Aa6r18/NTcHBwjbWhoaFu99WqVSvn1xEREbLZbM5pa9eu1eDBg/X+++/r0Ucf1Y4dO/Thhx/qpZde0vHjx/Xuu+86v/e+++5TZmam1q5dK0lyOBz605/+pBdffFE5OTm65JJL9Nhjj2n06NFV9mGM0fjx4xUfH69///vf8vM7eyKqd+/euvfeeyVJw4YNc7lOYadOnbRnzx7Nnz+/4QWhefPm6dZbb9W///1v9ejRQwEBAS7z77nnHq81BwBAQ3bukZWfGjFihJYvX+583rJlSxUXF1dZO3DgQGfYkKQOHTpUGotrjLmwZn9i2rRpevrpp9WpUye1aNHCre9JT0/XokWLtGDBAsXHx+uTTz7Rr371K8XExGjgwIGV6jMzM/Xll1/q1VdfdQlBFSIjI6tdV35+vqKiotx+PbXlcRB67bXX9OGHHyooKEhr166VzWZzzrPZbAQhAAAagdmzZ+vaa691u/7UqVP6wx/+oI8++khJSUmSzhy5+fTTT/XCCy9UGYS++eYbSVLXrl096i0rK0vPP/98nR8NkmoRhB555BE9/vjjmjZtWpXpDgAAqygsLKx2nr+/v8vzvLy8amt/+vd07969F9SXO/r27etRfVZWloqLiyuFp9LSUl1++eVVfk9tjmIdOHBAw4YN06233qpJkyZ5/P2e8jgIlZaWasyYMYQgAIDleTJup65qa+un6/Dz86sUXMrKypxfV4S+5cuXq02bNi51dru9ynVccsklkqSvvvqq2rB0roMHD2rw4MEaMGCAXnzxxfO/CC/wOM2kpKTojTfeqIteGgzuPg8AsJqYmBiXT5tJZ8b4VEhISJDdbte+ffvUpUsXl0e7du2qXGbv3r2VkJCgP//5z5UuwCxJx48fd3594MABDRo0SH369NHChQvr7YCLx0eEysvL9eSTT2rlypXq2bNnpcHSzzzzjNea85XU1FSlpqaqoKBAERERvm4HAIA6d8011+ipp57SK6+8oqSkJC1atEg7d+50Hslp3ry50tLSdP/998vhcOjnP/+58vPztX79eoWHhyslJaXSMm02mxYuXKihQ4fqqquu0iOPPKKuXbuqsLBQ//rXv/Thhx9q3bp1zhDUvn17Pf300zp8+LBzGed+Aq4ueByEduzY4dwoO3fudJl37sBpAADQeCQnJ+uxxx7Tgw8+qJKSEk2cOFG33367duzY4ax54oknFBMTo/T0dH333XeKjIzUFVdcoenTp1e73MTERG3ZskVz5szRpEmTdOTIEcXFxWnAgAHOG7mvWrVKWVlZysrKUtu2bV2+39uflvspm6nrNTRiFUeE8vPzFR4e7ut2AAA+UFJSouzsbHXs2NHlukDwvep+Np78/WbEMwAAsCyCEAAAsCyCEAAAsCyCEAAAsCyCEAAAbuCzRQ2PN34mBCEAAGpQcb286m6YCt+p+Jn89JqGnvD4OkIAAFiJv7+/IiMjnfcKCwkJ4bp5PmaMUXFxsfLy8hQZGVnpvm6eIAgBAHAeFVc3runGqah/kZGRF3zlaYIQAADnYbPZFBcXp5YtW7rciBS+ExAQcEFHgioQhAAAcJO/v79X/vii4WCwNAAAsCyCEAAAsCyCEAAAsCyCEAAAsCyCUBUyMjKUkJCgfv36+boVAABQh2yGa4ZXq6CgQBEREcrPz1d4eLiv2wEAAG7w5O83R4QAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYQAAIBlEYSqkJGRoYSEBPXr18/XrQAAgDpkM8YYXzfRUBUUFCgiIkL5+fkKDw/3dTsAAMANnvz95ogQAACwLIIQAACwLIIQAACwLIIQAACwLIIQAACwLIIQAACwLIIQAACwLIIQAACwLIIQAACwLIIQAACwrGa+bgAAADQixkglJVJhoXTiRPX/1jTv3H8DAqQDB3z2cghCAAA0ZeXlZwKHO6HE3X/Ly73XX0CA95ZVCwQhAAAaCmOkU6e8F1pOnJBOnqy7fkNCpObNpbCwC/vXGMlmq7s+a0AQAgCgthwOqajIO6eIKv49fbpuevX3v7Cw8tNpoaFnltnIEYQAANbhzaMthYVnQlBdCQ6ufUip6l+73WdHXRoyghAAoGFyOKTiYu+FlhMnpLKyuunVz6924aS6f0NDpWb8ia4PbGUAgHeUlbk/bsXdoy3G1E2vQUEXPq7l3H+Dgzna0kgRhICmZMcO6W9/q7sxBrAud47OlJbWzbptNu+FloqjLT7+pBIaDoIQ0JQ8+qi0dKmvu4DVBQZe2CDcn/4bEsLRFtQZghDQlFQM3Bw1SurRw7e9oOkJDT1/aAkLOxOEgEaCIAQ0RbfeKo0d6+suAKDB415jAADAsghCQFNSV5+wAYAmiiAEAAAsiyAENEV8wgYA3EIQApoSTo0BgEcIQgAAwLIIQkBTxKkxAHALQQgAAFgWQQhoShgjBAAeIQgBAADLIggBTRFjhADALQQhAABgWQQhoClhjBAAeIQgBAAALKvJB6H9+/dr0KBBSkhIUM+ePbV48WJftwTUPcYIAYBbmvm6gbrWrFkzzZ07V71791ZOTo769OmjESNGKDQ01NetAd7HqTEA8EiTD0JxcXGKi4uTJLVq1UrR0dE6duwYQQgAAPj+1Ngnn3yikSNHqnXr1rLZbHr33Xcr1WRkZKhDhw4KCgpS//79tWnTplqta+vWrSovL1e7du0usGuggePUGAC4xedHhIqKitSrVy9NnDhRN998c6X5b7zxhqZOnaoFCxaof//+mjt3rpKTk7Vnzx61bNlSktS7d2+dPn260vd++OGHat26tSTp2LFjuv322/XXv/61bl+QB4qKiqqd5+/vr6CgILdq/fz8FBwcXKva4uJimWpOp9hsNoWEhNSq9uTJk3I4HNX2ce4ROU9qS0pKVF5e7pXakJAQ2f4bGE6dOlXle6g2tcHBwfLzO/N/jNLSUpWVlXmlNigoSP7+/jXX/vf1BpWXy/+/k8rKylRaWlrtcu12u5o1a+Zx7enTp3Xq1KlqawMDAxUQEOBxbXl5uUpKSqqtDQgIUGBgoMe1DodDJ0+e9Epts2bNZLfbJUnGGBUXF3ul1pPfe/YRVdeyjzjPPqKK2oawj/Ap04BIMkuWLHGZlpiYaFJTU53Py8vLTevWrU16errbyy0pKTFXXXWVeeWVV85bl5+f73zs37/fSDL5+fkevQ53Sar2MWLECJfakJCQamsHDhzoUhsdHV1tbd++fV1q27dvX21tQkKCS21CQkK1te3bt3ep7du3b7W10dHRLrUDBw6stjYkJMSldsSIETVut3ONHj26xtrCwkJnbUpKSo21eXl5ztrJkyfXWJudne2sTUtLq7F2586dztqZM2fWWLtp0yZn7ZNPPllj7ZpZs5y18+bNq7F22bJlztqFCxfWWPvmm286a998880aaxcuXOisXbZsWY218+bNc9auWbOmxtonn3zSWbtp06Yaa2fOnOms3blzZ421aWlpztrs7OwaaydPnuyszcvLq7E2JSXFWVtYWFhj7ejRo13ewzXVso8482Afcfbh0T5izRpnbUPYR3hbfn6+kdz7++3zU2M1KS0t1datWzV06FDnND8/Pw0dOlQbNmxwaxnGGI0fP17XXHONfv3rX9dYm56eroiICOeDU2gAADRtNmMazsdMbDablixZoptuukmSdPDgQbVp00afffaZkpKSnHUPPvig1q1bp40bN553mZ9++qmuvvpq9ezZ0znt//7v/9SjR49KtadOnXI5jFdQUKB27dopPz9f4eHhF/DKqsZhb89rOex9nsPe110nffaZgt54Q/6/+IWkhnHYm1NjnBpjH3EWp8bOqMtTYwUFBYqIiHDr77fPxwjVtZ///Oc1/gKdy263O3dY9cGTT67VVe25OyZv1p67I/Vm7bk7fm/WevKz96Q2MDDQ+ce1Xmr/u2Nz/qszf+Dd3dl4UtusWTPnDs+btf7+/m6/hz2p9fPzq5Nam81WJ7US+4ja1LKP8Ly2IewjfKlBnxqLjo6Wv7+/cnNzXabn5uaqVatWPuoKaMAazgFeAGgUGnQQCgwMVJ8+fbR69WrnNIfDodWrV7ucKgMAAKgNnx+zKiwsVFZWlvN5dna2MjMzFRUVpYsvvlhTp05VSkqK+vbtq8TERM2dO1dFRUWaMGGCD7sGGjiuIwQAbvF5ENqyZYsGDx7sfD516lRJUkpKil566SWNGTNGhw8f1owZM5STk6PevXtrxYoVio2N9VXLQMPFqTEA8IjPg9CgQYOq/ZRBhSlTpmjKlCn11NGZK1lnZGTU+GkCAADQ+DXoMUK+kpqaql27dmnz5s2+bgWoHU6NAYBbCEIAAMCyCEJAU8IYIQDwCEEIAABYFkEIaIoYIwQAbiEIAQAAyyIIVSEjI0MJCQnq16+fr1sBPMMYIQDwCEGoCnx8Ho0ep8YAwC0EIQAAYFkEIQAAYFkEIaApYYwQAHiEIAQ0RYwRAgC3EIQAAIBlEYSApoRTYwDgEYJQFbiOEAAA1kAQqgLXEUKjxxghAHALQQgAAFgWQQhoShgjBAAeIQgBTRGnxgDALQQhAABgWQQhAABgWQQhoClhjBAAeIQgBDRFjBECALcQhAAAgGURhKrAlaXRaHFqDAA8QhCqAleWBgDAGpr5ugEAdYAxQgAaEGOMjh07pry8POXm5ro8Tp8+rT/96U8+640gBAAAPFZeXq4jR464hJqqgk7F9NOnT1e5nODgYP3xj3+UzUf/gSMIAU0JY4QAXIDS0lIdPny4yjDz07Bz5MgRORwOj5YfGRmp2NjYSg+HwyF/f/86elU1IwgBTRGnxgD818mTJ906apObm6sff/zRo2XbbDZFR0erZcuWVQacikfLli3VsmVL2e32OnqVtUcQAgCgETHGqLCw0K2jNrm5uTpx4oRHy/f3968UbKoLOtHR0WrWrHFHicbdPQAATYAxRj/++KNbR23y8vJ08uRJj5YfGBh43iM2FV9HRUXJz886HyonCAFNCWOEgAajvLxcR48ePe8Rm4rnZWVlHi0/NDT0vEdsKh7h4eE+G4zc0BGEgKaIHR5QJ8rKypwhpqajNrUdTBwREXHeIzYVj9DQ0Dp6ldZCEAIAWFpJSYlbY21yc3N17Ngxj5fvyWDioKCgOniFqAlBqAoZGRnKyMhQeXm5r1sBPMOpMcBlMPH5jtrk5eWpoKDAo+X7+/srJibGrcHEMTExjX4wcVPHT6cKqampSk1NVUFBgSIiInzdDgBYnjFGx48fd+uoTW5urseDiQMCAtweTHzRRRdZajBxU0cQApoixgihEagYTOzOUZu8vDyVlpZ6tPyQkBC3BxNHREQwmNiiCEIAAK8pKyurdGXi6oLO4cOHPR5MHB4e7vZg4rCwsDp6lWhKCEJAU2KMDkoqy82Vvv/e192gCanpo+Dnhp2jR496vOyLLrrIrcHEsbGxDCaG1xGEgCZk+sGDSpeklBRftwIL8/Pz82gwcUBAgK9bhoURhIAmZFNxsSSpmb+/mvHHBV5ks9kUFRXl9mBiX91AE/AUQQhogl75/e81Nj3d120AQIPH5/8AAIBlEYSAJsR5OUU+BgwAbiEIAQAAyyIIAU0Qx4MAwD0EIaAJMdxrDAA8QhACmiLGCAGAWwhCVcjIyFBCQoL69evn61aAWuGeSQDgHoJQFVJTU7Vr1y5t3rzZ160AAIA6RBACmhBGCAGAZwhCQBPEqTEAcA9BCAAAWBZBCGhC+Pg8AHiGIAQ0QTY/frUBwB3sLQEAgGU183UDDVnFaYaCggIfdwK45/R/37PFp07xvgVgWRX7P3eGC9gMgwqq9cMPP6hdu3a+bgMAANTC/v371bZt2xprCEI1cDgcOnjwoJo3b+71jyMXFBSoXbt22r9/v8LDw7267KaGbeU+tpX72FaeYXu5j23lvrraVsYYnThxQq1bt5bfecZMcmqsBn5+fudNkhcqPDycXxQ3sa3cx7ZyH9vKM2wv97Gt3FcX2yoiIsKtOgZLAwAAyyIIAQAAyyII+YjdbtfMmTNlt9t93UqDx7ZyH9vKfWwrz7C93Me2cl9D2FYMlgYAAJbFESEAAGBZBCEAAGBZBCEAAGBZBCEAAGBZBCEfycjIUIcOHRQUFKT+/ftr06ZNvm7Jp2bNmiWbzeby6Nq1q3N+SUmJUlNTddFFFyksLEy33HKLcnNzfdhx/frkk080cuRItW7dWjabTe+++67LfGOMZsyYobi4OAUHB2vo0KH65ptvXGqOHTumcePGKTw8XJGRkfrNb36jwsLCenwV9eN822r8+PGV3mvDhg1zqbHCtkpPT1e/fv3UvHlztWzZUjfddJP27NnjUuPO792+fft0/fXXKyQkRC1bttTvf/97nT59uj5fSr1wZ3sNGjSo0nvrzjvvdKmxwvaaP3++evbs6bxIYlJSkj744APn/Ib2viII+cAbb7yhqVOnaubMmdq2bZt69eql5ORk5eXl+bo1n+rWrZsOHTrkfHz66afOeffff7/+9a9/afHixVq3bp0OHjyom2++2Yfd1q+ioiL16tVLGRkZVc5/8skn9dxzz2nBggXauHGjQkNDlZycrJKSEmfNuHHj9OWXX2rVqlVatmyZPvnkE91xxx319RLqzfm2lSQNGzbM5b322muvucy3wrZat26dUlNT9fnnn2vVqlUqKyvTddddp6KiImfN+X7vysvLdf3116u0tFSfffaZXn75Zb300kuaMWOGL15SnXJne0nSpEmTXN5bTz75pHOeVbZX27Zt9cc//lFbt27Vli1bdM011+jGG2/Ul19+KakBvq8M6l1iYqJJTU11Pi8vLzetW7c26enpPuzKt2bOnGl69epV5bzjx4+bgIAAs3jxYue03bt3G0lmw4YN9dRhwyHJLFmyxPnc4XCYVq1amaeeeso57fjx48Zut5vXXnvNGGPMrl27jCSzefNmZ80HH3xgbDabOXDgQL31Xt9+uq2MMSYlJcXceOON1X6PVbdVXl6ekWTWrVtnjHHv9+799983fn5+Jicnx1kzf/58Ex4ebk6dOlW/L6Ce/XR7GWPMwIEDzb333lvt91h5e7Vo0cL87W9/a5DvK44I1bPS0lJt3bpVQ4cOdU7z8/PT0KFDtWHDBh925nvffPONWrdurU6dOmncuHHat2+fJGnr1q0qKytz2WZdu3bVxRdfbPltJknZ2dnKyclx2T4RERHq37+/c/ts2LBBkZGR6tu3r7Nm6NCh8vPz08aNG+u9Z19bu3atWrZsqUsvvVR33XWXjh496pxn1W2Vn58vSYqKipLk3u/dhg0b1KNHD8XGxjprkpOTVVBQ4Pzff1P10+1V4Z///Keio6PVvXt3PfzwwyouLnbOs+L2Ki8v1+uvv66ioiIlJSU1yPcVN12tZ0eOHFF5ebnLD1iSYmNj9dVXX/moK9/r37+/XnrpJV166aU6dOiQHn/8cV111VXauXOncnJyFBgYqMjISJfviY2NVU5Ojm8abkAqtkFV76mKeTk5OWrZsqXL/GbNmikqKspy23DYsGG6+eab1bFjR3377beaPn26hg8frg0bNsjf39+S28rhcOi+++7TlVdeqe7du0uSW793OTk5Vb7vKuY1VVVtL0n65S9/qfbt26t169bavn27HnroIe3Zs0fvvPOOJGttrx07digpKUklJSUKCwvTkiVLlJCQoMzMzAb3viIIoUEYPny48+uePXuqf//+at++vd58800FBwf7sDM0Nbfddpvz6x49eqhnz57q3Lmz1q5dqyFDhviwM99JTU3Vzp07XcbloXrVba9zx5H16NFDcXFxGjJkiL799lt17ty5vtv0qUsvvVSZmZnKz8/XW2+9pZSUFK1bt87XbVWJU2P1LDo6Wv7+/pVGyOfm5qpVq1Y+6qrhiYyM1CWXXKKsrCy1atVKpaWlOn78uEsN2+yMim1Q03uqVatWlQbjnz59WseOHbP8NuzUqZOio6OVlZUlyXrbasqUKVq2bJnWrFmjtm3bOqe783vXqlWrKt93FfOaouq2V1X69+8vSS7vLatsr8DAQHXp0kV9+vRRenq6evXqpWeffbZBvq8IQvUsMDBQffr00erVq53THA6HVq9eraSkJB921rAUFhbq22+/VVxcnPr06aOAgACXbbZnzx7t27ePbSapY8eOatWqlcv2KSgo0MaNG53bJykpScePH9fWrVudNR9//LEcDodzZ21VP/zwg44ePaq4uDhJ1tlWxhhNmTJFS5Ys0ccff6yOHTu6zHfn9y4pKUk7duxwCY6rVq1SeHi4EhIS6ueF1JPzba+qZGZmSpLLe8sq2+unHA6HTp061TDfV14ffo3zev31143dbjcvvfSS2bVrl7njjjtMZGSkywh5q3nggQfM2rVrTXZ2tlm/fr0ZOnSoiY6ONnl5ecYYY+68805z8cUXm48//ths2bLFJCUlmaSkJB93XX9OnDhhvvjiC/PFF18YSeaZZ54xX3zxhfn++++NMcb88Y9/NJGRkea9994z27dvNzfeeKPp2LGjOXnypHMZw4YNM5dffrnZuHGj+fTTT018fLwZO3asr15SnalpW504ccKkpaWZDRs2mOzsbPPRRx+ZK664wsTHx5uSkhLnMqywre666y4TERFh1q5daw4dOuR8FBcXO2vO93t3+vRp0717d3PdddeZzMxMs2LFChMTE2MefvhhX7ykOnW+7ZWVlWVmz55ttmzZYrKzs817771nOnXqZK6++mrnMqyyvaZNm2bWrVtnsrOzzfbt2820adOMzWYzH374oTGm4b2vCEI+8vzzz5uLL77YBAYGmsTERPP555/7uiWfGjNmjImLizOBgYGmTZs2ZsyYMSYrK8s5/+TJk2by5MmmRYsWJiQkxIwaNcocOnTIhx3XrzVr1hhJlR4pKSnGmDMfoX/sscdMbGyssdvtZsiQIWbPnj0uyzh69KgZO3asCQsLM+Hh4WbChAnmxIkTPng1daumbVVcXGyuu+46ExMTYwICAkz79u3NpEmTKv0nxArbqqptJMksXLjQWePO793evXvN8OHDTXBwsImOjjYPPPCAKSsrq+dXU/fOt7327dtnrr76ahMVFWXsdrvp0qWL+f3vf2/y8/NdlmOF7TVx4kTTvn17ExgYaGJiYsyQIUOcIciYhve+shljjPePMwEAADR8jBECAACWRRACAACWRRACAACWRRACAACWRRACAACWRRACAACWRRACAACWRRACgPOw2Wx69913fd0GgDpAEALQoI0fP142m63SY9iwYb5uDUAT0MzXDQDA+QwbNkwLFy50mWa3233UDYCmhCNCABo8u92uVq1auTxatGgh6cxpq/nz52v48OEKDg5Wp06d9NZbb7l8/44dO3TNNdcoODhYF110ke644w4VFha61PzjH/9Qt27dZLfbFRcXpylTprjMP3LkiEaNGqWQkBDFx8dr6dKlznk//vijxo0bp5iYGAUHBys+Pr5ScAPQMBGEADR6jz32mG655Rb95z//0bhx43Tbbbdp9+7dkqSioiIlJyerRYsW2rx5sxYvXqyPPvrIJejMnz9fqampuuOOO7Rjxw4tXbpUXbp0cVnH448/rl/84hfavn27RowYoXHjxunYsWPO9e/atUsffPCBdu/erfnz5ys6Orr+NgCA2quTW7kCgJekpKQYf39/Exoa6vKYM2eOMebMXcHvvPNOl+/p37+/ueuuu4wxxrz44oumRYsWprCw0Dl/+fLlxs/Pz3nX+datW5tHHnmk2h4kmUcffdT5vLCw0EgyH3zwgTHGmJEjR5oJEyZ45wUDqFeMEQLQ4A0ePFjz5893mRYVFeX8OikpyWVeUlKSMjMzJUm7d+9Wr169FBoa6px/5ZVXyuFwaM+ePbLZbDp48KCGDBlSYw89e/Z0fh0aGqrw8HDl5eVJku666y7dcsst2rZtm6677jrddNNNGjBgQK1eK4D6RRAC0OCFhoZWOlXlLcHBwW7VBQQEuDy32WxyOBySpOHDh+v777/X+++/r1WrVmnIkCFKTU3V008/7fV+AXgXY4QANHqff/55peeXXXaZJOmyyy7Tf/7zHxUVFTnnr1+/Xn5+frr00kvVvHlzdejQQatXr76gHmJiYpSSkqJFixZp7ty5evHFFy9oeQDqB0eEADR4p06dUk5Ojsu0Zs2aOQckL168WH379tXPf/5z/fOf/9SmTZv097//XZI0btw4zZw5UykpKZo1a5YOHz6su+++W7/+9a8VGxsrSZo1a5buvPNOtWzZUsOHD9eJEye0fv163X333W71N2PGDPXp00fdunXTqVOntGzZMmcQA9CwEYQANHgrVqxQXFycy7RLL71UX331laQzn+h6/fXXNXnyZMXFxem1115TQkKCJCkkJEQrV67Uvffeq379+ikkJES33HKLnnnmGeeyUlJSVFJSor/85S9KS0tTdHS0Ro8e7XZ/gYGBevjhh7V3714FBwfrqquu0uuvv+6FVw6grtmMMcbXTQBAbdlsNi1ZskQ33XSTr1sB0AgxRggAAFgWQQgAAFgWY4QANGqc3QdwITgiBAAALIsgBAAALIsgBAAALIsgBAAALIsgBAAALIsgBAAALIsgBAAALIsgBAAALIsgBAAALOv/A/QvZaA8KZfpAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# reopen saved data using callbacks in fnamevar\n", + "lines = open(fnamevar, \"r\").readlines()\n", + "# read output data in fnamevar\n", + "Chat = np.array(\n", + " [\n", + " np.fromstring(\n", + " min(re.findall(re.escape(\"[\") + \"(.*?)\" + re.escape(\"]\"), line), key=len),\n", + " sep=\",\",\n", + " )\n", + " for line in lines\n", + " ]\n", + ")\n", + "l, c = Chat.shape\n", + "plt.semilogy(range(0, l * 100, 100), Chat[:, 0], \"r-\")\n", + "plt.semilogy(range(0, l * 100, 100), Chat[:, 1], \"k-\")\n", + "plt.semilogy(range(0, l * 100, 100), np.ones(Chat[:, 0].shape) * C1true, \"r--\")\n", + "plt.semilogy(range(0, l * 100, 100), np.ones(Chat[:, 1].shape) * C2true, \"k--\")\n", + "plt.legend([\"C1hat\", \"C2hat\", \"True C1\", \"True C2\"], loc=\"right\")\n", + "plt.xlabel(\"Epochs\")\n", + "plt.title(\"Variables\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the velocity distribution of the flow field:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for t in range(0, 8):\n", + " t = t * unit_of_t\n", + " [ob_x, ob_y, ob_t, ob_u, ob_v, ob_p] = load_training_data(num=140000)\n", + "\n", + " xyt_pred = {\"x\": ob_x, \"y\": ob_y, \"t\": t * np.ones((len(ob_x),))}\n", + " uvp_pred = model.predict(xyt_pred)\n", + "\n", + " x_pred, y_pred, t_pred = xyt_pred['x'], xyt_pred['y'], xyt_pred['t']\n", + " u_pred, v_pred, p_pred = uvp_pred['u'], uvp_pred['v'], uvp_pred['p']\n", + " x_true = ob_x[ob_t == t]\n", + " y_true = ob_y[ob_t == t]\n", + " u_true = ob_u[ob_t == t]\n", + " fig, ax = plt.subplots(2, 1)\n", + " cntr0 = ax[0].tricontourf(x_pred.mantissa, y_pred.mantissa, u_pred.mantissa, levels=80, cmap=\"rainbow\")\n", + " cb0 = plt.colorbar(cntr0, ax=ax[0])\n", + " cntr1 = ax[1].tricontourf(x_true.mantissa, y_true.mantissa, u_true.mantissa, levels=80, cmap=\"rainbow\")\n", + " cb1 = plt.colorbar(cntr1, ax=ax[1])\n", + " ax[0].set_title(\"u-PINN \" + \"(t=\" + str(t) + \")\", fontsize=9.5)\n", + " ax[0].axis(\"scaled\")\n", + " ax[0].set_xlabel(\"X\", fontsize=7.5, family=\"Arial\")\n", + " ax[0].set_ylabel(\"Y\", fontsize=7.5, family=\"Arial\")\n", + " ax[1].set_title(\"u-Reference solution \" + \"(t=\" + str(t) + \")\", fontsize=9.5)\n", + " ax[1].axis(\"scaled\")\n", + " ax[1].set_xlabel(\"X\", fontsize=7.5, family=\"Arial\")\n", + " ax[1].set_ylabel(\"Y\", fontsize=7.5, family=\"Arial\")\n", + " fig.tight_layout()\n", + " plt.show()" ] } ],