diffphys ns phiflow2 update, WIP
This commit is contained in:
NT committed
1 parent
feb1477391
commit
6b350893fc
8 files changed
+1069
-374
No files matched your search
@@ -20,16 +20,19 @@
|
||||
sections:
|
||||
- file: diffphys-code-gradient.ipynb
|
||||
- file: diffphys-discuss.md
|
||||
- file: diffphys-code-ns-v2a.ipynb
|
||||
- file: diffphys-code-ns.ipynb
|
||||
- file: diffphys-code-sol.ipynb
|
||||
- file: diffphys-dpvspinn.md
|
||||
- file: diffphys-outlook.md
|
||||
- file: old-phiflow1.md
|
||||
sections:
|
||||
- file: overview-burgers-forw-v1.ipynb
|
||||
- file: overview-ns-forw-v1.ipynb
|
||||
- file: diffphys-code-ns.ipynb
|
||||
- file: physicalloss-code-v2.ipynb
|
||||
- file: diffphys-code-gradient-v1.ipynb
|
||||
- file: diffphys-code-tf.ipynb
|
||||
- file: diffphys-code-ns-v1.ipynb
|
||||
- file: diffphys-code-ns-v2a.ipynb
|
||||
- file: jupyter-book-reference
|
||||
sections:
|
||||
- file: jupyter-book-reference-markdown
|
||||
|
||||
+402
-304
@@ -8,8 +8,14 @@
|
||||
"source": [
|
||||
"# Differentiable Fluid Simulations\n",
|
||||
"\n",
|
||||
"<br>\n",
|
||||
"\n",
|
||||
"**TODO , NEXT update diffphys example, then check burgers examples**\n",
|
||||
"**TODO , NEXT finish opt case, diverges? see colab example , finish texts before end.**\n",
|
||||
"integration done!\n",
|
||||
"\n",
|
||||
"**(older: Differentiable_Fluid_Simulations_with_phiflow2_test)**\n",
|
||||
"\n",
|
||||
"<br><br>\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
@@ -27,7 +33,7 @@
|
||||
"\n",
|
||||
"As optimization objective we'll consider a more difficult variant of the previous example: the state of the observed density $d$ should match a given target after $n=20$ steps of simulation. In contrast to before, the marker $d$ cannot be modified in any way, but only the initial state of the velocity $\\mathbf{u}$ at $t=0$. This gives us a split between observable quantities for the loss formulation, and quantities that we can interact with during the optimization (or later on via NNs).\n",
|
||||
"\n",
|
||||
"First, let's get the loading of python modules out of the way:"
|
||||
"First, let's get the loading of python modules out of the way. By importing `phi.tf.flow` below, we get fluid simulation functions that work within tensorflow graphs and can provide gradients."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -36,48 +42,87 @@
|
||||
"metadata": {
|
||||
"id": "da1uZcDXdVcF"
|
||||
},
|
||||
"outputs": [],
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"/Users/thuerey/miniconda3/envs/tf/lib/python3.8/_collections_abc.py:743: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.\n",
|
||||
" for key in self._mapping:\n",
|
||||
"/Users/thuerey/miniconda3/envs/tf/lib/python3.8/_collections_abc.py:744: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.\n",
|
||||
" yield (key, self._mapping[key])\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"!pip install --upgrade --quiet phiflow\n",
|
||||
"from phi.flow import * # The Dash GUI is not supported on Google Colab, ignore the warning\n",
|
||||
"#!pip install --upgrade --quiet phiflow\n",
|
||||
"#!pip install --upgrade --quiet git+https://github.com/tum-pbs/PhiFlow@develop\n",
|
||||
"\n",
|
||||
"from phi.tf.flow import * # The Dash GUI is not supported on Google Colab, ignore the warning\n",
|
||||
"import pylab"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"id": "BVV1IKVqDfLl"
|
||||
"id": "rdSTbMoaS0Uz"
|
||||
},
|
||||
"source": [
|
||||
"## Setting up the simulation\n",
|
||||
"## Batched Simulations\n",
|
||||
"\n",
|
||||
"To make things a bit more interesting - and to move a bit closer to a NN training process - let's set up of four fluid simulations that run in parallel, i.e. a mini batch similar to DL training. In phiflow we can directly pass a `batch_size=4` parameter to the `Fluid` object. Each fluid simulation is fully independent. In this case they differ by having circular Inflows at different locations.\n",
|
||||
"Now we can setup the simulation, which will work in line with the previous \"regular\" simulation example from the {doc}`overview-ns-forw`. However, now we'll directly include an additional dimension, in line with a mini-batch used for NN training. For this, we'll include an additional named dimension called `inflow_loc`. This dimension will exist \"above\" the previous spatial dimensions `x`, `y` and the channel dimensions `vector`. As indicated by the name `inflow_loc`, the main differences for this dimension will lie in different locations of the inflow, in order to obtain different flow simulations. The named dimensions in phiflow make it very convenient to broadcast information\n",
|
||||
"\n",
|
||||
"Like before, let's plot the marker density after a few steps of simulation (each call to `step()` now updates all four simulations). Note that the boundaries between the four simulations are not visible in the image, but it shows four completely separate density states. The different inflow positions in conjunction with the solid wall boundaries (zero Dirichlet for velocity, and Neumann for pressure), result in four different end states of the simulation."
|
||||
"The `Domain` object is allocated just like before, but the `INFLOW_LOCATION` tensor now receives a string\n",
|
||||
"`'inflow_loc,vector'` that indicates the names of the two dimensions.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"id": "WrA3IXDxv31P"
|
||||
"id": "mphMP0sYIOz-"
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"<matplotlib.image.AxesImage at 0x1283dad10>"
|
||||
]
|
||||
},
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Inflow tensor shape and dimensions: (inflow_loc=4, x=32, y=40)\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"DOMAIN = Domain(x=32, y=40, boundaries=CLOSED, bounds=Box[0:32, 0:40])\n",
|
||||
"\n",
|
||||
"INFLOW_LOCATION = math.tensor([(10, 5), (12, 5), (14, 5), (16, 5)], 'inflow_loc,vector', convert=True)\n",
|
||||
"INFLOW = DOMAIN.grid(Sphere(center=INFLOW_LOCATION, radius=3)) * (1./3.)\n",
|
||||
"print(f\"Inflow tensor shape and dimensions: {INFLOW.shape}\")\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The last `print` verifies that our `INFLOW` grid likewise has an `inflow_loc` dimension in addition to the two spatial ones. You can test for the existence of a tensor dimension in phiflow with the `.exists` boolean, which can be evaluated for any dimension name. E.g., above `INFLOW.inflow_loc.exists` will give `True`, while `INFLOW.unknown_dim.exists` will give `False`.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We can easily simulate a few steps now starting with these different initial conditions. Thanks to the broadcasting, the exact same code we used for the single forward simulation in the overview chapter will produce four simulations with different smoke inflow positions."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAACICAYAAAD+r7D/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2deZRc1XWvf7uGLvWk7lZrniWQACGmWBiwyUsA2wzGhveeg7G9bN4KKyTr2XGSR+xAnDjLiZdN4njCdnAUjC28zGAwZlBssBDY2ICEEAJJaEASaO5uDT0PVdVVd78/9r13H6mqVT2rbnt/a/Xq07vurTr33Nun9nT2IWaGYRiGET1ip7sDhmEYxvCwCdwwDCOi2ARuGIYRUWwCNwzDiCg2gRuGYUQUm8ANwzAiSmI8P4yILGfRMAxj6Bxj5mknC8d1Ahfi4/+RhmEYkSa/r5jUXCiGYRgRxSZwwzCMiGITuGEYRkSxCdwwDCOi2ARuGIYRUWwCNwzDiCg2gRuGYUQUm8ANwzAiik3ghmEYEcUmcMMwjIhiE7hhGEZEOQ21UIxoI7VsCBRKGEGNsvxp6M/EgJx/xWRyii9T/ao/3w0A8Lzu8e3YBCIenwwAqKwoqAkFAOjLHgUA5POd49ankWIauGEYRkSxCdwwDCOiRN6F4pqeSxpuBAAsozNDWcw39Xfw3lC2ve2RsM3IjXEPo4qW/b2y7rNh+/azkgCAuTU9oWxPZy0A4NF9es7Drd8J2/l8l9+ycvCCjNMZDdeHkruXnBW2r/64mPKU1PFsXiPP6Tdfnx/Kvtd0X9juyxzyWzbGAJBMTAUAfGraX4Syf77kIABg1scbip7T9EAbAOCL6+eGsvuPfj9s9+eOjXo/R4pp4IZhGBGFmMfvG1t25Bn5hg4LGq4J27+6eEHYXvJF0U68xQv1YPYAALHtu0PR8ZVvh+3PvSDn//jo3aEsSkGM0UEDkgsargYA/OY9c0LZ/H9YEra5tkYanqenV4hWjp7eUOSteSNs/8u9CwEAX937H6GsHLWZsWRK7YVh+4Fl7wcAfODLk0IZT50Stqmvr/AN4mJpcsL5/3l9V9j80lenAwC+sve7oSyXax1ZpyOAa4H/zymfC9sPf14sktj7L9Bjs/3S6O8v/mZJeY45eJ5x4nP80X+T/4mft34tlI2fBZ/fyMwrTpaaBm4YhhFRbAI3DMOIKCVdKEQ0CcALAFKQoOejzPxPRLQIwEMAGgFsBPBJZs6WeK8RuVD+bPY/AgDu2Xy+vmdzi7aP+Saja94HJmelY65O0SCGN3s2AGDn9atD2aUvPQMA6O7bM+y+RonFDR8K27u+MlMa7z47lNHho3pwzs/1dk35Ct+MTao5yw112p4jY/zkdZtC2Z9sEnfKRHZXBXnHAPDaH30ybC//prikqNlxI7nPbHWl/M45efX5Ijn2wXFA6GL5yic1T/yLu//df+teTDzk+buw/pZQsvHuqrBNZ8ozh7YuFJAYYA7KFRnjhtqwybsPAwDe9Vkdz9fbV/mtsV4DMXwXSgbAlcx8AYALAVxDRJcC+FcA32TmMwG0Abh1NLtrGIZhnJqSEzgLwdd60v9hAFcCeNSXrwJw45j00DAMwyjKoLJQiCgOcZOcCeB7AL4GYJ2vfYOI5gH4JTMvL/E+g3ShaFbEVXV/E7bXPOS7QZIaJcaRtsLTJzmvpyrkd4WT8k7O91ZtNQDAmz0zFHX/49MAgKWPqQvlSOf6QfQ7OixsuC5s7/ny7LBNiyWbgd8+EsrY02eEKvz7l3DG0HedUFWFyuLOffbNUJ49IxRtuHUbAOCKV34WyvoyB4d6GWVJkIO84Q9vDmXn/8ussM07D0jDHddaxx3S6Jvt3WmVZQbInAioE/cBL5oXir5+YxMA4O923qWffWovZ2S4rO7TAIDfvf6uUEZvaXYZH/TdU8XcJUSFMgAoNhc6bhWaK/eVly4OZZdfuBEA8HLH9wbV7+EzgiwUZs4z84UA5gJ4N4CzS5wSQkS3EdGrRPTqoPtqGIZhlGRIKzGZuZ2IngdwGYB6Ikowcw4ysR8a4JyVAFYCgQZemnhcAweP39SkL2RFu+h/aa8eW+dofSm/0FKVauBU7X+DVjjfxK4GnxGNJJbXfM7qe/8UAHDXa8+Hsls3b9RrivDqTSIZr2+d4QSCz0qF7dxLYnVwn6N5VOrYxSbL2FHKeXRiotHwUSdg5GroPaJJkpN/e/EPzwUAfOBi1TyfyKimGGXOrfkwAOD8f5oeyrx1mrPt9crzQxU6RrFKfY7peLHAmxzLbboCFlkncObLKa1jfPvtEiD++t+qltrS8fKgr6PcCJ5dAPivFTKV0DYd1/wbjgXnP5M0yXl2p/uWoFd8GqKYjLF3RMef084Yt+wDAMQdrTzox3nPad/G08opqYET0TQiqvfblQDeD2A7gOcBfMQ/7BYAT4xVJw3DMIxCBqOBzwKwyveDxwD8lJlXE9E2AA8R0ZcBbALwgzHsp2EYhnESJSdwZt4M4KIi8rch/vBR56Laj4bt6ps0uJh5WpbDe31qAsVSjonjmzsUc4IU1YXvz8c6CmSUUbMnvkHc9bf8rR53x2fUDI1yQHPGZLllH/5/msva//z+sN0XOsJ0DJNV6jKKd8kYx6rUeKNkoSGX79JzYpXyerxPx5j8vPw7l+uy8ade1DzeqOUux2La96+fI8vivS37QlnfLg1IptvEDRWLa+53aq/m2ldeLvWquVlNea9LXCNeWs9x/w/itTLGsTb9nPhV5wEAPlqnwePvdGwI21FzBc6vvzJsL/tzcftlf703lOWd8ci0ieskWaUupcQBGZvEDMft6pBrkecz5yxN6O/VZzvVIM9+vE0/c9mfy/w0/zXt2762p0tfzChhKzENwzAiik3ghmEYEaXM6oGL2f4X89Xkyz6rVQSbtos/pDejJlDVXjXLG2dIJL76PCfjxI/Oex2ZUORGljkjJml8lr5O3dLmD7w3lN1Qo8vv/6szyIiMyhZiTl596o8AAF6zJg21bNYyA/taZQl8nvWc889oDtvV/jDEGjVzhTvFTG3frJ/Y1qGZRKmkjNMs6DLvmF+m4OK7zwtls6+4LGwfbF9b8qrKifrqpWH7ir+UZ/L44/psdnTqsvp1zX4+sXNfrph/OGxXBS7ASv339Frkmdz9mj6HB7vVP7i4Xtx9c89S+5+fkRty5wr9f7jncH3Yjk5FSOn/LY3qseUuWafQulvngtZOHY8Nx2ScplSom+jKC8VVWHG2Pu8uuV1yv55zaq63ZvUeXDxV1pxMmayZQDPO7ivo2z+3rXHedWznCNPADcMwIkpZaeBEotVdOVc1vn2bVHN55YhoLuuP6/fOR+erVjfzTGlTrWqHyImGnX5HtaEtb2pgNGDpXNVG6lf4NcTfUu3/c+dp0az7WuSbPoqFmD6xULTlzk06HltbdAeSXx8V7aTaeTKWzXEsGk80mvxRDZZ17JDXdx6eGsqeaVYNPEi7/UxSg3r18yRIOekyR0NKqQZ+P57zW9HYYWZhXBfJ9W+RgOSho42hbG2zas7PNYk2Pa9an9OzJ+t4zfTzkPv3aYC35S155n7bonXDHz+g92B+tVitWhEbqG+S86ffqO897bllYftw+wslr6sciMVkrcDSGn1WWp8RLXj/MR2PB/dpAbV1HXIP/rBBNzA+Z7+Mw6JFxQPkTf7rvz2q9+W3bRpcfr1N3utjC3QdSfKZ1oK+Bf0Fxn4TatPADcMwIopN4IZhGBGlrFwolRViBtY3qOm4YbcWWnpwr5j/6/K/CmU3zHlf2PZ8r4B3XM/fvVbyczce1QI0qw9p8OhoVo79VJduzfa/6yXAlsrpMt2OPg3+xPwlvVEJYbom3bsXS7Ds0CG9npeOa1Dn8XYpMnUOdGNozwloBrm2nFXXxoZ94oJ5YK+almvSz4TtdK4dALCg+k9C2dVZKZEw9+Xtoawjq0vPo+I6CZjMNWG76x1xbbzdpbIXW9TE/l3/UwCAhT0a+LqoTYtQXdB2HACQbddxf61FzPfvHNICa3t6tdTDzJhs2Xb2/otD2VUzJejW6OTft3S9PviLKhM8T9wl156j6xVajoq7Y/1xdbG+1KFuzq3p/wYAdHS8J5Sd2yx12Bd0adErlxebZYxXd+wIZfvSL2k/8EEAwOLj+pxOTmUL+ua97pQ7GGNMAzcMw4goZaWBx+N+SlBMtS83ne1t2gsAyPTrCsmOfr2EWJBR5BSr+f5O+bb8cduToay1a0vYrq0SzbzuwEdC2aJq0fovadRgasypw9WfK1LCtoypnqQbFLe3i0VyuEdTrg726Oq+tzufleOqD4SydJ9WCa5PihXU16r35fkjosE/dPxboYxI78u5NbLjz4taoRbn1YkGNatFhfETynwG7XLXxCWAu6RaA4U93RKE7OrX4K9bS60nLdbHtsx/h7K3u/9v2I41+BbebtWcf3NEZNvbfx7KqlJaovYiEg28LatjmM3Lh3a8oBrhrFrV0KOSqlmRFM246Yhq28G1pT293vkJDWhu9X/v6tRVkeuPnQUA+PgAn7P+WKzgnIqEfmbw/u5npnOJgr4F/QWAbL/z0I8BpoEbhmFEFJvADcMwIkpZuVDmpKRmFuf0e6U2qcVoFniyQmo/azGppONuqZgvprzXreesan0MANDevTWU1deoS+BD1TfIe9e4O3dkcDIvHtOgH5EE65ijEcaMx1IFsrSnY9yYcsa7Usb4SKfuv9Gd1kJesWoJ0DU3qcvgiQ4ZWzcvfvrkS8L2hxolQHxhvboEJiVk7PKO+2YH73V6GPSpvMc4eBZqkk6g1zex651VgIuctQnT+6UWu5uH3ZNzXEW+C7CvR4PCj3e9AgBg1tzvmpSuZ1juu10umaLPbiImY5vp03/zdD5a7j8ASPhB+GB1JQBc0CDP2pIafaaO9OmqzMqsrElwd3la1yOrj2NLine truncated
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAABHgAAAFtCAYAAAB80d3MAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8S0lEQVR4nO3deZgdBZnv8d/b3dn3jZCQsC8CCgEjExQEARUZFUbFgRnnMuoVXHC5V51x1FHE8Y77NnLBIBhGHRUFRLyMyiarLAHCEkIWIED2tZPuLL2+949T0Sb2eau765w6p7q/n+fpJ93nd6rqTXXX23XernOOubsAAAAAAABQXA21LgAAAAAAAADZMOABAAAAAAAoOAY8AAAAAAAABceABwAAAAAAoOAY8AAAAAAAABQcAx4AAAAAAICCY8AzhJnZdDO7y8xazOwbZnaJmf04h+3ea2bHVXs7KTUcY2b31bIGYKgZ4j1nupktMbMRtawDGEqGeM/hPAdDyhA/3nM9x9h7X+exTfQdA55BxsxWmtkZfbz7hZI2SRrv7h+vYll/YmZvkdTi7o8mX7/czH5nZpvMzHu5/2Qzu8HMdpjZ82b2d/3Y1gIzazez1h4fjZLk7o9Lak7qATBAg6nnmNkIM7sq6TUtZrbIzN7Uj21dYmYde/WcgyXJ3ddLukOlfQBggAZTz9lrucPMbHd/HpBynoPBbjAe72Z2XjKM2WFmz5jZyX3cVkXPMZJ9uytZz7qkn4zt4+K572v0HQOeoe0ASU+5e9kTjip4v6Qf9fi6Q9K1kt5b5v6XSWqXNF3S30u63MyO7sf2vuruY3t8dPXIfiLpon6sC0A29d5zmiS9KOkUSRMkfVbStWZ2YD+29/O9es6zPTJ6DpCveu85PV0m6aEBbI/zHKCk7o93M3u9pK9IerekcZJeK+nZ3u5bRqXPMd7i7mMlzZF0nKR/6eNyA97XZtbU32XQPwx4BjEz+0czu8fMvm5mW83suT1/jTazBZIukPRPyeT2L6bjZvZWM1tsZs1m9gczOzK5/d1mdlOP+y03s1/0+PpFM5vTy/qGSzpN0p17bnP3pe5+laTFvdx/jKS3S/pXd29193sk/VrSPwxwl+ztD5JOz+tyRmCwK3rPcfcd7n6Ju6909253/42k5yS9csA75aUekHSwmR1QofUBQ1rRe06P5c6T1Czptv7ugxR/EOc5GCQGyfH+BUmXuvv9yXnGandfPYDd0ZsBn2O4+zpJv1Np0CNJMrN5ZnZfsr8eM7NTk9sXaK99bWYNZvap5IqkzWZ2rZlNTu5/oJm5mb3XzF6QdHty+3usdCXTVitd9XRAj227mb0/+V40m9llZmY98vcly7aY2VNmdnxy+0wzu87MNiY/Hx/p774YDBjwDH5/JWmppKmSvirpKjMzd/9HlSa9e/7yc2vPhczscEk/lfQxSdMk3SzppqSZ3Snp5ORgnilpuKQTk+UOljRW0uO91HKYpG53X9XH2g+X1Onuy3rc9pik/lzB80Ez22JmD5vZ23sGSUPtkHREP9YHIFbknvMSZjZdpT5U9oFZL96S9JzFZvaBnoG7d0paIenYgdQDoFeF7jlmNl7SpZL+d5//xy/FeQ6GksIe71Z6+uRcSdPMbIWZrTKz75nZqH78/6tyjmFmsyS9KVleZrafpP8n6d8kTZb0CUnXmdm0Mvv6w5LOUekK6JmStqp0VWJPp0g6UtIbzexsSZ+W9DaVvh93q/T96enNkl4l6RhJ75T0xqS2cyVdIul/SBov6a2SNptZg6SbVHqsuJ+k0yV9zMze2N/9UXQMeAa/5939yuSS3WskzVDp6U5p/lbS/3P3W9y9Q9LXJY2S9OrkcsAWlaa8r1Vp4rvGzF6m0sF7t7t397LOiclyfTVW0va9btum0iWNffFdlZrvPpL+VdICM3vNXvdpSeoCUBlF7jl/YmbDVDqBucbdn+7jYteqdPIyTdL7JH3OzM7f6z70HKCyit5zvijpqgEOojnPwVBT5ON9uqRhkt4h6WT9+WlRn+3j8tU4x/iVmbWo9PT0DZI+n9z+Lkk3u/vNyZVGt0haKOmsMut5v6TPuPsqd29TaQDzDnvp07EuSa6U3pXc/9/dfUkymPo/kubsdfXRl9292d1fUOn1heYkt/9PlYZLD3nJCnd/XqVh0DR3v9Td25Pv65WSzuvH/hgUeA7c4LduzyfuvjO5uq0vL6A1U9LzPZbtNrMXVZqISqVp96mSDk0+b1apCZ6oHpcq7mWr+j6ckaRWlSazPY1XH5upuz/S48ubzewnKk2K7+1x+ziVagdQGUXuOZKk5K9AP1Lp9b8u7uty7v5Ujy/vM7PvqHQi1/OvUvQcoLIK23OSp32codKDvH7jPAdDUGGPd0m7kn//w93XSpKZfVOlAc9n0hau0jnGOe5+q5mdIum/VLoyqlml19g51176Iu3DVBq09OYASTeYWc9BWJdeOnx7ca/7f8de+g5cptL3Y8/3aV2PbKf+/H2eLemZMjXMNLPmHrc1qnR10JDCFTwoZ41KB4okKXne42xJe54nuqcRnpx8fqdKjfAUlW+EK5JV7Vcm39sySU1mdliP245V/54u0ZOr1Dwk/enyw+EqXeoJoLbqoefs2e5VKp2UvD35S99A7d1zmlQ6eXwswzoBVEY99JxTJR0o6QUzW6fS0yDebmaPRAsFOM8Belfz493dt0papdJx+qeb+/w/6GWVqtA5hrvfKWmBSlc2SaVhzI/cfWKPjzHu/uUyq3hR0pv2uv/IvV5fyPe6/0V73X+Uu9/Xh3JflHRImduf22ud49y93FVHgxYDHpRzraS/NrPTk6cqfFxSm6Q9B96dkl4naVRyWfHdks6UNEXSo72t0N3bJd2qUrOUVOqKZjZSpRMQmdlIS14M0N13SLpe0qVmNia57PhsJa9Wb39+0a4De9uemb3DzMYmz6l9g0qXG/66x11OkXR7cikhgNqqec9JXK7SJdBvSS4jfomk55za2/bM7Gwzm5Rs4wRJH5F0Y4+7nCBpZXIpMYDaqoeeM1+lBypzko8rVHrdiz2vNcF5DlAZ9XC8S9IPJX3YzPYxs0mS/pek3/RYvmLnGGZ2qpV5q/Yyvi3p9WZ2rKQfq/R6P280s8bk/3GqlV6rpzdXSPqSJU+xMrNpVnqdnXKukPQvlrwzsplNsNJr6/TFDyR9wsxemeyLQ5PtPiipxcz+2cxGJXW/3Mxe1cf1DhoMeNArd1+q0onCf0jaJOktKj3gaU/yZSo9heru5OvtKr3N373+0rfo3Nv39dJ3wTpApUsW91yVs0sv/UvTB1V6juwGlS5B/IC777nvbJUu4yv36vMfTbJmSV+T9D53/0OP/O9VajAAaqweek5ygnCRSg+01lnp3SFazezvk3y2Sk8RfaLMts5T6S96LZL+U9JX3P2aHjk9B6gT9dBz3H2nu6/b85Fsb7e7b0zuy3kOUAH1cLwnvijpIZWepbBEpeHRl6SqnGPM1p8HWKmSvvOfkj7n7i+q9Ef1T0vaqNLVMZ9U+dnBd1QaLv/eSq/pc79KL4hdbls3qPR28T8zs+2SnlTpRZ77UucvVNpn/6XSvviVpMnJ9+nNKp3DPafS9/kHkib0Zb2DiXn/374eyMTM7pV0sbv3OhHvx3o+K2mju39/AMseI+n77n5ilhoA1L8K9px3STra3f9lAMvuo9JfCI9z991Z6gBQ3zjPAYaOej3HMLMfSPqFu/8uS10oHgY8AAAAAAAABcdTtAAAAAAAAAqOAQ8AAAAAAEDBMeABAAAAAAAouKY8N9bPt2oDMHhscvdpeW+UngMMWTXpORJ9Bxiq3N1qsV16DjBk9Xquk+uAp6Qx/00CqLGu52u3bXoOMPTUsudI9B1gqIneSTsP9Bxg6On9XIenaAEAAAAAABQcAx4AAAAAAICCY8ADAAAAAABQcAx4AAAAAAAACo4BDwAAAAAAQMEx4AEAAAAAACg4BjwAAAAAAAAFx4AHAAAAAACg4BjwAAAAAAAAFBwDHgAAAAAAgIJjwAMAAAAAAFBwDHgAAAAAAAAKjgEPAAAAAABAwTHgAQAAAAAAKDgGPAAAAAAAAAXHgAcAAAAAAKDgGPAAAAAAAAAUHAMeAAAAAACAgmPAAwAAAAAAUHAMeAAAAAAAAAqOAQ8AAAAAAEDBMeABAAAAAAAoOAY8AAAAAAAABceABwAAAAAAoOAY8AAAAAAAABQcAx4AAAAAAICCY8ADAAAAAABQcAx4AAAAAAAACo4BDwAAAAAAQMEx4AEAAAAAACg4BjwAAAAAAAAFlzrgMbORZvagmT1mZovN7AvJ7QvM7DkzW5R8zKl6tQAGPXoOgLzRdwDkiZ4DoFqa+nCfNkmnuXurmQ2TdI+Z/XeSfdLdf1m98gAMQfQcAHmj7wDIEz0HQFWkDnjc3SW1Jl8OSz68mkUBGLroOQDyRt8BkCd6DoBq6dNr8JhZo5ktkrRB0i3u/kASfcnMHjezb5nZiDLLXmhmC81sYWVKBjDY0XMA5I2+AyBP9BwA1WClAXIf72w2UdINkj4sabOkdZKGS5ov6Rl3vzRleZcaB1wsgKLqetjd5/Z3KXoOgIEZWM+R6DsABqJL7m4DWZKeA2Bgej/X6de7aLl7s6Q7JJ3p7mu9pE3SDyWdUJE6ASBBzwGQN/oOgDzRcwBUUl/eRWtaMlmWmY2S9HpJT5vZjOQ2k3SOpCerVyaAoYKeAyBv9B0AeaLnAKiWvryL1gxJ15hZo0oDoWvd/TdmdruZTZNkkhZJen/1ygQwhNBzAOSNvgMgT/QcAFXRr9fgybwxniMKDFEDfz2MLOg5wFBVm54j0XeAoWngr8GTFT0HGKoq8Bo8AAAAAAAAqD8MeAAAAAAAAAqOAQ8AAAAAAEDBMeABAAAAAAAoOAY8AAAAAAAABceABwAAAAAAoOAY8AAAAAAAABQcAx4AAAAAAICCY8ADAAAAAABQcAxLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 432x288 with 1 Axes>"
|
||||
"<Figure size 1152x360 with 4 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
@@ -87,78 +132,113 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"world = World()\n",
|
||||
"fluid = world.add(Fluid(Domain([40, 32], boundaries=CLOSED), buoyancy_factor=0.05, batch_size=4), physics=IncompressibleFlow())\n",
|
||||
"centers = [[5,10], [5,12], [5,14], [5,16]]\n",
|
||||
"world.add(Inflow(Sphere(center=centers, radius=3), rate=0.2));\n",
|
||||
"smoke = DOMAIN.scalar_grid(0) # sampled at cell centers\n",
|
||||
"velocity = DOMAIN.staggered_grid(0) # sampled in staggered form at face centers \n",
|
||||
"\n",
|
||||
"for frame in range(20):\n",
|
||||
" world.step(dt=1.5)\n",
|
||||
"def step(smoke, velocity):\n",
|
||||
" smoke = advect.mac_cormack(smoke, velocity, dt=1) + INFLOW\n",
|
||||
" buoyancy_force = smoke * (0, 1) >> velocity\n",
|
||||
" velocity = advect.semi_lagrangian(velocity, velocity, dt=1) + buoyancy_force\n",
|
||||
" velocity, _, _, _ = fluid.make_incompressible(velocity, DOMAIN)\n",
|
||||
" return smoke, velocity\n",
|
||||
"\n",
|
||||
"pylab.imshow(np.concatenate(fluid.density.data[...,0], axis=1), origin='lower', cmap='magma')"
|
||||
"for _ in range(20):\n",
|
||||
" smoke,velocity = step(smoke,velocity)\n",
|
||||
"\n",
|
||||
"fig, axes = pylab.subplots(1, 4, figsize=(16, 5))\n",
|
||||
"for i in range(INFLOW.shape.inflow_loc):\n",
|
||||
" axes[i].imshow(smoke.values.inflow_loc[i].numpy('y,x'), origin='lower', cmap='magma')\n",
|
||||
" axes[i].set_title(f\"Inflow {math.to_int(INFLOW_LOCATION).inflow_loc[i].vector.unstack_spatial('x,y')}\" + (\", Reference\" if i==3 else \"\"))\n",
|
||||
"pylab.tight_layout()\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Now we see four simulations, (0) to (3). This final density of simulation (0), with the curved plume on the far left, will be our **reference state**, while the initial velocity of the other three will be modified in the optimization procedure below."
|
||||
"The last image shows the state of the advected smoke fields after 20 simulation steps. The final smoke shape of simulation `[3]` with an inflow at `(16,5)`, with the straight plume on the far right, will be our **reference state** below. The initial velocity of the other three will be modified in the optimization procedure below to match this reference.\n",
|
||||
"\n",
|
||||
"(A small side hint: phiflow tensors will keep track of their chain of operations using the backend they were created for. E.g. a tensor created with NumPy will keep using NumPy/SciPy operations unless a TensorFlow tensor is also passed to the same operation. Thus, it is a good idea to verify that tensors have the right `type` once in a while.)\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"id": "rdSTbMoaS0Uz"
|
||||
},
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Differentiation\n",
|
||||
"## Gradients\n",
|
||||
"\n",
|
||||
"The simulation we just computed was using purely (non-differentiable) operations from numpy.\n",
|
||||
"To enable differentiability, we need to build a TensorFlow graph that computes this result.\n",
|
||||
"Let's look at how to get gradients from our simulation. The first trivial step taken care of above was to include `phi.tf.flow` to import differentiable operators from which to build our simulator.\n",
|
||||
"\n",
|
||||
"(Note, the first line of the next block, `%tensorflow_version 1.x` is only necessary when running in environments that by default have newer tensorflow versions installed, e.g., `colab`. Uncomment if you're running this notebook there.)"
|
||||
"Now we want to optimize the initial velocities so that all simulations arrive at a final state that is similar to the right simulation where the inflow is located at `(16, 5)`.\n",
|
||||
"To achieve this, we record the gradients during the simulation and define a simple $L^2$ based loss function. The loss function we'll use is given by $L = | s - s_r |^2$, where $s$ denotes the smoke density, and $s_r$ the reference state from the fourth simulation above. When evaluating the loss function we treat this reference state as an external constant via `field.stop_gradient()`.\n",
|
||||
"\n",
|
||||
"It is important that our initial velocity has the `inflow_loc` dimension before we record the gradients, such that we have the full \"mini-batch\" of four versions of our velocity (three of which will be updated via gradients in our optimization later on). To get the approprate velocity tensor, we multiply a `staggered_grid` tensor from our domain object, with a tensor full of ones along the `inflow_loc` dimension."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"id": "mphMP0sYIOz-"
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Could not load resample cuda libraries: CUDA binaries not found at /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/cuda/build/resample.so. Run \"python setup.py cuda\" to compile them\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/util.py:119: The name tf.AUTO_REUSE is deprecated. Please use tf.compat.v1.AUTO_REUSE instead.\n",
|
||||
"\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/profiling.py:12: The name tf.RunOptions is deprecated. Please use tf.compat.v1.RunOptions instead.\n",
|
||||
"\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/profiling.py:13: The name tf.RunMetadata is deprecated. Please use tf.compat.v1.RunMetadata instead.\n",
|
||||
"\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/session.py:17: The name tf.Session is deprecated. Please use tf.compat.v1.Session instead.\n",
|
||||
"\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/session.py:18: The name tf.get_default_graph is deprecated. Please use tf.compat.v1.get_default_graph instead.\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"/Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/flow.py:15: UserWarning: TensorFlow-CUDA solver is not available. To compile it, download phiflow sources and run\n",
|
||||
"$ python setup.py tf_cuda\n",
|
||||
"before reinstalling phiflow.\n",
|
||||
" warnings.warn(\"TensorFlow-CUDA solver is not available. To compile it, download phiflow sources and run\\n$ python setup.py tf_cuda\\nbefore reinstalling phiflow.\")\n"
|
||||
"Velocity tensor shape: (inflow_loc=4, x=32, y=40, vector=2)\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/ops/math_grad.py:297: setdiff1d (from tensorflow.python.ops.array_ops) is deprecated and will be removed after 2018-11-30.\n",
|
||||
"Instructions for updating:\n",
|
||||
"This op will be removed after the deprecation date. Please switch to tf.sets.difference().\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# %tensorflow_version 1.x\n",
|
||||
"from phi.tf.flow import * # Causes deprecation warnings with TF 1.15\n",
|
||||
"import pylab\n",
|
||||
"session = Session(None) # Used to run the TensorFlow graph"
|
||||
"smoke = DOMAIN.grid(math.zeros(inflow_loc=4))\n",
|
||||
"velocity = initial_velocity = DOMAIN.staggered_grid(0) * math.ones(inflow_loc=4)\n",
|
||||
"print(\"Velocity tensor shape: \"+format(velocity.shape))\n",
|
||||
"\n",
|
||||
"with math.record_gradients(initial_velocity.values):\n",
|
||||
" for _ in range(5):\n",
|
||||
" smoke,velocity = step(smoke,velocity)\n",
|
||||
"\n",
|
||||
" loss = field.l2_loss(smoke - field.stop_gradient(smoke.inflow_loc[-1]))\n",
|
||||
" grad = math.gradients(loss, initial_velocity.values)\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Computing the actual gradient is not very difficult. We simply call `math.gradients` for the initial velocity w.r.t. our loss function to obtain the `grad` tensor.\n",
|
||||
"\n",
|
||||
"Let's direcly take a look at it. The following images show the `x` components (via `vector[0]`) of the staggered grid gradient in `grad`."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAABHYAAAFgCAYAAADehv6hAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABa40lEQVR4nO3de7xkVX3n/c+XqyYKDXSLLWgaA06CiSEPHSQXlQFEMo8TnAkiymgzAyGMMZNMEgcMGTVEMjgmIU9mjEkrCF6BYNROQiSAEjXxQqMEQaO0RENjc2mgESVCGn7PH3uf7upz6nTXOXVq17l83q9Xvbr22mvvWnUu395n1Vprp6qQJEmSJEnSwrPbuBsgSZIkSZKk2bFjR5IkSZIkaYGyY0eSJEmSJGmBsmNHkiRJkiRpgbJjR5IkSZIkaYGyY0eSJEmSJGmBsmNHkiRJkiRpkiQnJvlqkg1Jzu2z/4VJvpBka5KTJ+1bk+T29rGmp/zIJF9qz/lHSTJsO+3YkSRJkiRJ6pFkd+DtwM8ChwOvTHL4pGr/DJwOfGDSsfsDbwKeDxwFvCnJfu3udwC/ABzWPk4ctq127EiSJEmSJO3oKGBDVd1RVY8BlwMn9Vaoqm9U1S3AE5OOfQlwbVU9UFUPAtcCJyZZCexTVZ+tqgLeA7xs2IbuMewJJM3OiUltHvIcN8E1VTV0D6+kxc28kdSVYfPGrJE0qDnIm9uA7/UUra2qtT3bBwF39mxvpBmBM4h+xx7UPjb2KR+KHTvSmGwG1g95jsDyuWiLpMXNvJHUlWHzxqyRNKg5yJvvVdXquWrPODkVSxqn3XYb7iFJgzJvJHVlxFkzwGKmeye5ot3/uSSreva9oS3/apKXzN2bljQWo82bu4Bn9mwf3JYNc+xd7fPZnHNaXqlJ4+QfWpK6Yt5I6soIs2bAxUzPAB6sqkOBi4C3tsceDpwKPJdmsdI/bs8naaEa7bXNjcBhSQ5JshdNfqwbsGXXACck2a9dNPkEmqmmm4BvJzm6vRvWa4CPzvyN78grNUmSJEkLxS4XM223L2ufXwUc1/4BdRJweVU9WlX/BGxozydJU1TVVuB1NJ00XwGurKrbkpyf5OcAkvxEko3Ay4E/TXJbe+wDwO/QdA7dCJzflgG8FngXTQZ9HfjrYdvqGjvSuCR+Ci6pG+aNpK4MnzfLk/QumzGbxUy31amqrUkeAg5oyz876dihFy2VNCYdXN9U1dXA1ZPK3tjz/EZ2nFrVW+8S4JI+5euBH5nLdtqxI42Tf2hJ6op5I6krw+XN5sWymKmkDnh9A9ixI42Pn6BL6op5I6kro8+bQRYznaizMckewL7A/QMeK2mh8PpmG78KkiRJkhaKQRYzXQesaZ+fDHy8qqotP7W9a9YhwGHA5ztqtySNjCN2pHGyh1lSV8wbSV0ZYd60a+ZMLGa6O3DJxGKmwPqqWgdcDLw3yQbgAZrOH9p6VwJfBrYCv1RVj4+ssZJGz+sbwI4dabwMIkldMW8kdWX8i5l+j+YONf2OvQC4YKQNlNQdr28AO3ak8XFOqKSumDeSumLeSOqKebONHTvSOBlEkrpi3kjqinkjqSvmDeDiyZIkSZIkSQuWI3akcXHooKSumDeSumLeSOqKebONHTvSOBlEkrpi3kjqinkjqSvmDWDHjjReBpGkrpg3krpi3kjqinkDuMaOJEmSJEnSguWIHWlcnBMqqSvmjaSumDeSumLebGPHjjROBpGkrpg3krpi3kjqinkD2LEjjY89zJK6Yt5I6op5I6kr5s02duxI42QQSeqKeSOpK+aNpK6YN4CLJ0uSJEmSJC1YjtiRxskeZkldMW8kdcW8kdQV8wawY0caH+eESuqKeSOpK+aNpK6YN9vYsSONk0EkqSvmjaSumDeSumLeAK6xIy1qSU5M8tUkG5Kc22f/C5N8IcnWJCdP2rcmye3tY013rZYkSZIkDcoRO9K4jHjoYJLdgbcDLwY2AjcmWVdVX+6p9s/A6cBvTDp2f+BNwGqggJvaYx8cWYMljY5DlSV1xbyR1BXzZhs7dqRxGm0QHQVsqKo7AJJcDpwEbOvYqapvtPuemHTsS4Brq+qBdv+1wInAB0fZYEkj5IWPpK6YN5K6Yt4AduxI4zXaIDoIuLNneyPw/CGOPWiO2iVpHLzwkdQV80ZSV8wbwI4daXzmZujg8iTre7bXVtXaYU8qaZFxqLKkrpg3krpi3mzjV0Fa2DZX1eqeR2+nzl3AM3u2D27LBjHMsZKWoAEWa987yRXt/s8lWdWz7w1t+VeTvKTThkuSJC1wjtiRxmm0Pcw3AoclOYSmU+ZU4FUDHnsN8LtJ9mu3TwDeMPdNlNSZ8S/WfgbwYFUdmuRU4K3AK5IcTpNPzwWeAVyX5DlV9fjIGixptPwEXVJXzBvAETvS+EwMHRzmsRNVtRV4HU0nzVeAK6vqtiTnJ/m5pgn5iSQbgZcDf5rktvbYB4DfoekcuhE4f2IhZUkL0Ijzhp7F2qvqMWBisfZeJwGXtc+vAo5Lkrb88qp6tKr+CdjQnk/SQjRs3kjSoMybbRyxI43TiAOlqq4Grp5U9sae5zfSTLPqd+wlwCUjbaCk7ox2Ta9BFmvfVqeqtiZ5CDigLf/spGNdrF1ayBbZH0yS5jHzBrBjR5IkDWZzVa0edyMkSZK0Izt2pHGyh1lSV0abN4MsuD5RZ2OSPYB9gfsHPFbSQuL1jaSumDfAIl9jJ8k3khzfPv/NJO/q6HWT5N1JHkzy+S5ec1yS/Nck9yT5TpIDklSSQ0f8mi9J8pFRvsaA7fhQkp8d4gSjXvNCHTJvRm+J583vJ/mvQ5xg1HmzbbH2JHvRLIa8blKddcCa9vnJwMerqtryU9u7Zh0CHAYs6p/lYZg1o7fEs2a4a5vmJF7bLBLmzegt8bwZ7tqmOYl50xrbu0lyanu70+8mubd9/tp2IcU5V1W/W1VnDnueJKvaX7idjXb6GZo7gxxcVSNfALJt0yeSPJLkHycCeJbnOj3JpwesuyfwB8AJVfWUqrp/tq87QxcAF/a043eSfCnJ1iRv7tPOVyX5Zvuz9pEk+w/yIj3f6+/0PP5nT5W3Am8Z6p3YsdMJ82bumDfT502S/zfJp5NsSXJ3kncleeqgL9R+r7/bkze9F9C/B/xm22kyOyPMm0EWawcuBg5IsgH4NeDc9tjbgCuBLwMfA35pod4Ry6yZO2bNzq9teupdMpM/BDu5tgGvbTpg3swd82aXf0utSPKBJA+l6Wx7/6AvNPJrGxh53iQ5MclXk2xIcm6f/XsnuaLd/7kkq9ry05Lc3PN4IskR7b4b2nNO7HvaUF8DxtSxk+TXgf8PeBvwdOBA4Gzgp4G+39g0t1JdKH4A+EZVfXemB+4i5KbzQeCLNItQngdclWTFLM4zUwcCTwJu6+C1gOYuTsC+VdW70OYG4H8Af9Wn/nOBPwVeTdPeR4A/nuHLLmvD9ilV9TsThVX1eWCfJLNbc8IRO50wb6Zn3uzcTPOGZmrRW2hu2f3DNAsAv22GL/tjPXmz7QK6qjYB/wj83PSH7kQHeVNVV1fVc6rqB6vqgrbsjVW1rn3+vap6eVUdWlVHVdUdPcde0B73b6rqr2f1HsfMrJmeWbNzs8iaieN+BvjBWb7saK5tmoZ5bTNi5s30zJudm2Xe/DlwN/As4Gk0HTIzMZprGxh53rS/N28HfhY4HHhlksMnVTsDeLCqDgUuoukgp6reX1VHVNURNH+L/lNV3dxz3GkT+6vq3ll/DVqdp2eSfYHzgddW1VVV9XA1vlhVp1XVo229S5O8I8nVSb4L/Ns0n4Z+Mcm3k9w5uUcxyavTjMy4P8l5k/a9Ocn7eraPTvL3aT5Z/Yckx/Tsu6Htufy7JA8n+Zsky9vdn2z/3dL2Ov7kpNc5A3gX8JPt/t9uy3+h7cV7IMm6JM/oOaaS/FKS24Hb+3zN/jrJ6yaV/UOS/5jkOcD/A7ypqv6lqj4EfAn4+V18KwaSZgjmbyS5pe2lvSLJk9rX/WrP1+LjfY7dN8l7ktzXfl9+K8lu7b5vJjmyfX5a+zV4brt9RqYfHvizwN/2FlTVZe0fAg/3qX8a8BdV9cmq+g7wP4H/mBl8ir4LNwD/7xydS3PMvDFv2n2d5E1VfaCqPlZVj1TVg8A7aS6y58oNmDfzkllj1rT7Osma9nx7AP8H+OWZfwV26QbMmnnLvDFv2n2d5E2SE2jWwXt9VT1UVf9aVV+c3VejrxuY33lzFLChqu6oqseAy4GTJtU5CbisfX4VcFwyZeTcK9tjR2Yc3eI/CewNfHSAuq+iGSr2VODTwHeB1wDLaH4A/muSlwGk6Tl7B01v2DNoelz73sY5yUE0PZJvAfYHfgP4UHbsmX0V8J9peiX3ausAvLD9d+KTjs/0nruqLqbpMf9Mu/9NSY4F/hdwCrAS+CZTv7Evo7k17OQeQGh6kV/Z0/7DaXqy/wp4LnBHVfX+Iv5DWz5XTgFOBA4BngecXlVf63mNZVV1bJ/j/g/NJ9jPBl5E8737z+2+vwWOaZ+/CLiD7V/bFzEpcHr8KNtDcBDPpfl6AFBVXwceA54zg3N8M8nGNHN9l0/a9xXgx2Zwrh05YmfUzBvzBrrLm8leyMw/hftkmmlcf552KG8P82b+MmvMGug2a/478MmqumWGx00Y3bUNmDWjZd6YN9Bd3hzd1r+s7fC7McmLZnA8jPLaBobNm+VJ1vc8zpp09oOAO3u2N7ZlfetUMzX9IZrfn16voPk57PXuNNOw/mefjqAZG0d6Lqe5ZerWiYKe3t5/SfLCnrofraq/q6onqhnCfUNVfandvoXmizPxg3Uy8JftyIxHaUZmPDFNG/4TcHU1w8afqKprgfXAv+up8+6q+lpV/QvN3P8jhnjPpwGXVNUX2ra9gaYXelVPnf9VVQ+0rzfZh4EjkvxAz/n+vD3XU2h+eHo9RBPgc+WPqupbVfUA8BcM8LVLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 1152x360 with 8 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"fig, axes = pylab.subplots(1, 4, figsize=(16, 5))\n",
|
||||
"for i in range(INFLOW.shape.inflow_loc):\n",
|
||||
" im = axes[i].imshow(grad.inflow_loc[i].vector[0].numpy('y,x'), origin='lower', cmap='bwr')\n",
|
||||
" pylab.colorbar(im, ax=axes[i])\n",
|
||||
" axes[i].set_title(f\"Gradient for v0 - Inflow {math.to_int(INFLOW_LOCATION).inflow_loc[i].vector.unstack_spatial('x,y')}\")\n",
|
||||
"pylab.tight_layout()\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -167,33 +247,262 @@
|
||||
"id": "3mpyowRYUSS4"
|
||||
},
|
||||
"source": [
|
||||
"Now that we have imported `phi.tf.fluid`, let's set up the simulation just like before. But now, we want to start from a velocity that we can modify, i.e. a variable. Then we can optimize these initial velocities so that all simulations arrive at a final state that is similar to the first simulation from the previous example. I.e., the state shown in the left-most image above.\n",
|
||||
"Not surprisingly, the fourth gradient on the left is zero (it's already matching the reference). The other three gradients have detected variations for the intial round inflow posiitons shown as positive and negative regions around the cirucular shape.\n",
|
||||
"\n",
|
||||
"This is a fairly tough task: we're producing diffent dynamics by changing the boundary conditions (the marker inflow position), and an optimizer should now find a single initial velocity state, that gives the same state as simulation `0` above at $t=30$. Thus, after 20 steps with $\\Delta t=1.5$ the simulation should reproduce a different set of boundary conditions from the velocity state. It would be much easier to simply change the position of the marker inflow to arrive at this goal, but -- to make things a bit more difficult -- the inflow is _not_ a degree of freedom. The optimizer can only change the velocity $\\mathbf{u}$ at time $t=0$.\n",
|
||||
"This is just the linearized change that is points in the direction of an increasing loss. Now we can proceed by updating the initial velocities in the opposite direction to minimize the loss.\n",
|
||||
"\n",
|
||||
"To achieve this, we create a TensorFlow variable for the velocity at t=0.\n",
|
||||
"It is initialized with zeros (like with the NumPy simulation above) and can later be used as a target for optimization."
|
||||
"Note that this is a fairly tough task: we're producing different dynamics by changing the boundary conditions (the marker inflow position), and an optimizer should now find a single initial velocity state for each of the three simulations on the right, that gives the same state as simulation `0` above at $t=20$. Thus, after 20 update steps the simulation should reproduce a different set of boundary conditions from the velocity state. It would be much easier to simply change the position of the marker inflow to arrive at this goal, but -- to make things a bit more difficult and interesting here -- the inflow is _not_ a degree of freedom. The optimizer can only change the velocity $\\mathbf{u}$ at time $t=0$.\n",
|
||||
"\n",
|
||||
"The following cell re-starts the simulations, and iterates a few times to optimize the initial velocities."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {
|
||||
"id": "NlJMJikaHOL6"
|
||||
},
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/session.py:28: The name tf.global_variables_initializer is deprecated. Please use tf.compat.v1.global_variables_initializer instead.\n",
|
||||
"\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/session.py:29: The name tf.train.Saver is deprecated. Please use tf.compat.v1.train.Saver instead.\n",
|
||||
"\n"
|
||||
"65.591805\n",
|
||||
"65.37787\n",
|
||||
"65.1976\n",
|
||||
"64.12548\n",
|
||||
"64.667114\n",
|
||||
"64.404465\n",
|
||||
"74.69762\n",
|
||||
"76.985146\n",
|
||||
"78.58485\n",
|
||||
"82.51491\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"initial_velocity = DOMAIN.staggered_grid(0) * math.ones(inflow_loc=4)\n",
|
||||
"LR = 1e-02\n",
|
||||
"\n",
|
||||
"for _ in range(10):\n",
|
||||
" velocity = initial_velocity\n",
|
||||
" smoke = DOMAIN.grid(math.zeros(inflow_loc=4)) + INFLOW\n",
|
||||
"\n",
|
||||
" with math.record_gradients(initial_velocity.values):\n",
|
||||
"\n",
|
||||
" buoyancy_force = smoke * (0, 1) >> velocity\n",
|
||||
" velocity, _, _, _ = fluid.make_incompressible(velocity + buoyancy_force, DOMAIN)\n",
|
||||
" for _ in range(20):\n",
|
||||
" smoke,velocity = step(smoke,velocity)\n",
|
||||
"\n",
|
||||
" loss = field.l2_loss(smoke - field.stop_gradient(smoke.inflow_loc[0]))\n",
|
||||
" print(format(loss))\n",
|
||||
" grad = math.gradients(loss, initial_velocity.values)\n",
|
||||
" initial_velocity = initial_velocity - LR * grad\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The loss should have gone down significantly, from above 60 to below 40, and now we can visualize how well the reconstruction for the last frame turned out for simulations 0 to 2."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"fig, axes = pylab.subplots(1, 4, figsize=(16, 5))\n",
|
||||
"for i in range(INFLOW.shape.inflow_loc):\n",
|
||||
" im = axes[i].imshow(grad.inflow_loc[i].vector[0].numpy('y,x'), origin='lower', cmap='bwr')\n",
|
||||
" pylab.colorbar(im, ax=axes[i])\n",
|
||||
" axes[i].set_title(f\"Gradient for v0 - Inflow {math.to_int(INFLOW_LOCATION).inflow_loc[i].vector.unstack_spatial('x,y')}\")\n",
|
||||
"pylab.tight_layout()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Naturally, the image on the right is the same (this is simulation (3), our reference), and the other three simulations now exhibit a noticeable shift towards the right. \n",
|
||||
"\n",
|
||||
"**XXX TODO** This is especially visible for the far right simulation (3), which was moving straight up in its original state, and now has an off-center plume shape that tries to match the reference on the left. \n",
|
||||
"\n",
|
||||
"Note that all three simulations need to \"work\" with a fixed inflow, hence they cannot simply \"produce\" marker density out of the blue to match the target. Also each simulation needs to take into account how the non-linear model equations change the state of the system over the course of 20 time steps. So the optimization goal is quite difficult, and it is not possible to exactly satisfy the constraints to match the reference simulation in this scenario. This is exactly what makes it an interesting test case, though.\n",
|
||||
"\n",
|
||||
"**XXX TODO** Now we can also have a look at the now-optimized initial velocity field. Here we see the y-components (from `data[...,0]`).\n",
|
||||
"\n",
|
||||
"... re-simulate "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# re-run sim\n",
|
||||
"smoke = DOMAIN.grid(math.zeros(inflow_loc=4)) + INFLOW\n",
|
||||
"velocity = initial_velocity\n",
|
||||
"buoyancy_force = smoke * (0, 1) >> velocity\n",
|
||||
"velocity, _, _, _ = fluid.make_incompressible(velocity + buoyancy_force, DOMAIN)\n",
|
||||
"\n",
|
||||
"for _ in range(20):\n",
|
||||
" smoke,velocity = step(smoke,velocity)\n",
|
||||
"\n",
|
||||
"fig, axes = pylab.subplots(1, 4, figsize=(16, 5))\n",
|
||||
"for i in range(INFLOW.shape.inflow_loc):\n",
|
||||
" axes[i].imshow(smoke.values.inflow_loc[i].numpy('y,x'), origin='lower', cmap='magma')\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This example illustrates how the differentiable physics approach can easily be extended towards significanlty more complex PDEs. Above, we've optimized for a mini-batch of 20 steps of a full Navier-Stokes solver.\n",
|
||||
"\n",
|
||||
"This is a powerful basis to bring NNs into the picture! As you might have noticed, our degrees of freedom were still a regular grid, and we've solved a single inverse problem. There was more than one case to solve, but nonetheless no question of generalization and no neural network structure.\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"<br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"--- \n",
|
||||
"\n",
|
||||
"## END UP TO HERE , REMOVE REST?\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"xxx\n",
|
||||
"\n",
|
||||
"#use keras GD?\n",
|
||||
"\n",
|
||||
"# \"reset\" necessary???\n",
|
||||
"#smoke = DOMAIN.grid(math.zeros(inflow_loc=4))\n",
|
||||
"if 0:\n",
|
||||
" velocity = initial_velocity = DOMAIN.staggered_grid(0) * math.ones(inflow_loc=4)\n",
|
||||
" smoke = DOMAIN.grid(math.zeros(inflow_loc=4)) + INFLOW\n",
|
||||
"\n",
|
||||
" with math.record_gradients(initial_velocity.values):\n",
|
||||
" buoyancy_force = smoke * (0, 1) >> velocity\n",
|
||||
" velocity, _, _, _ = fluid.make_incompressible(velocity + buoyancy_force, DOMAIN)\n",
|
||||
" for _ in range(2):\n",
|
||||
" smoke,velocity = step(smoke,velocity)\n",
|
||||
"\n",
|
||||
" loss = field.l2_loss(smoke - field.stop_gradient(smoke.inflow_loc[-1]))\n",
|
||||
" print(format(loss)) \n",
|
||||
" grad = math.gradients(loss, initial_velocity.values)\n",
|
||||
" print(type(grad))\n",
|
||||
"\n",
|
||||
" print(\"ok\")\n",
|
||||
"# ---\n",
|
||||
"\n",
|
||||
"# https://www.tensorflow.org/api_docs/python/tf/keras/optimizers/Optimizer\n",
|
||||
"opt = tf.keras.optimizers.SGD(learning_rate=1e-02)\n",
|
||||
"#opt = tf.keras.optimizers.Adam(learning_rate=0.1)\n",
|
||||
" \n",
|
||||
"velocity = initial_velocity = DOMAIN.staggered_grid(0) * math.ones(inflow_loc=4)\n",
|
||||
"smoke = DOMAIN.grid(math.zeros(inflow_loc=4)) + INFLOW\n",
|
||||
"\n",
|
||||
"#with math.record_gradients(initial_velocity.values): # PHI\n",
|
||||
"with tf.GradientTape() as tape: # same\n",
|
||||
" buoyancy_force = smoke * (0, 1) >> velocity\n",
|
||||
" velocity, _, _, _ = fluid.make_incompressible(velocity + buoyancy_force, DOMAIN)\n",
|
||||
" for _ in range(2):\n",
|
||||
" smoke,velocity = step(smoke,velocity)\n",
|
||||
"\n",
|
||||
" loss = field.l2_loss(smoke - field.stop_gradient(smoke.inflow_loc[-1]))\n",
|
||||
" #print(type(loss.native()))\n",
|
||||
" print(\"Loss \"+format(loss))\n",
|
||||
" \n",
|
||||
" #grad = math.gradients(loss, initial_velocity.values) # PHI\n",
|
||||
" \n",
|
||||
" #print(type(grad.vector[0]))\n",
|
||||
" print(type(grad.vector[0].native()))\n",
|
||||
" #print(format(grad.vector[0].shape))\n",
|
||||
" print(type(initial_velocity.vector[0].values.native()))\n",
|
||||
" #print(format(initial_velocity.vector[0].shape))\n",
|
||||
"\n",
|
||||
" grads = tape.gradient(loss.native(), [initial_velocity.vector[0].values.native()] ) # works\n",
|
||||
" \n",
|
||||
" # a = grad.vector[0].native() + initial_velocity.vector[0].values.native() # works PHI\n",
|
||||
" #opt.apply_gradients(zip([grad.vector[0].native()], [initial_velocity.vector[0].values.native()] ))\n",
|
||||
"\n",
|
||||
" opt.apply_gradients(zip(grads, [initial_velocity.vector[0].values.native()] ))\n",
|
||||
" # -> same error as opt.minimize?\n",
|
||||
" \n",
|
||||
" # minimize:\n",
|
||||
" #lossL = lambda: (field.l2_loss(smoke - field.stop_gradient(smoke.inflow_loc[-1]))).native()\n",
|
||||
" #opt.minimize(lossL, [initial_velocity.vector[0].values.native()])\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"#initial_velocity = initial_velocity - 1e-02 * grad\n",
|
||||
"#step_count = opt.minimize(loss, [var1]).numpy()\n",
|
||||
"#step_count = opt.minimize(loss, [initial_velocity.values.native()]).numpy()\n",
|
||||
"\n",
|
||||
"#print(step_count)\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"..."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "NlJMJikaHOL6"
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"xxx UP TO HERE xxx\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"world = World()\n",
|
||||
"fluid = world.add(Fluid(Domain([40, 32], boundaries=CLOSED), buoyancy_factor=0.05, batch_size=4), physics=IncompressibleFlow())\n",
|
||||
"world.add(Inflow(Sphere(center=centers, radius=3), rate=0.2));\n",
|
||||
@@ -217,44 +526,11 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "wSrIezfWHjcQ"
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/tf_backend.py:167: The name tf.sparse_tensor_dense_matmul is deprecated. Please use tf.sparse.sparse_dense_matmul instead.\n",
|
||||
"\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/tf_backend.py:61: The name tf.div_no_nan is deprecated. Please use tf.math.divide_no_nan instead.\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"<matplotlib.image.AxesImage at 0x143f6c290>"
|
||||
]
|
||||
},
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAACICAYAAAD+r7D/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAK10lEQVR4nO3db4wcd33H8fendhL+SbUDyHLttHGFBTKhOJGFHBVVbQjCoRTnQVSlRaqrRrIqUTVUSChpnhSpUotaQanUEllJGreKEmhIGzeif4zrij4oBhtMSGJSO6UUR04MAgOlEmDz7YOZSzbnvdze3e7e/S7vl3TamdnZnd+f1edmfzOzk6pCktSen1juAkiSFscAl6RGGeCS1CgDXJIaZYBLUqMMcElq1NppbiyJ5yxK0sJ9s6peO3vhVAO8s2b6m5Skpl342rClDqFIUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjZo3wJO8LMnnknwpyeNJPtgv35LkSJJTST6e5NLJF1eSNGOUPfAfANdV1ZuB7cCuJDuBDwEfqarXAd8GbplcMSVJs80b4NX53372kv6vgOuAB/vl+4EbJ1JCSdJQI42BJ1mT5DhwFjgIPAWcq6rz/SqngU2TKaIkaZiRAryqLlTVdmAz8BbgDaNuIMneJEeTHF1kGSVJQyzoLJSqOgccBq4F1iVZ2z+1GXh6jtfsq6odVbVjSSWVJL3AKGehvDbJun765cDbgRN0QX5Tv9oe4OFJFVKSdLG186/CRmB/kjV0gf+JqnokyRPAA0n+EPgicPcEyylJmiVVNb2NJQVrprY9SVodLhwbNgztlZiS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUqFHuiXlFksNJnkjyeJJb++WXJzmY5GT/uH7yxZUkzRhlD/w88P6q2gbsBN6bZBtwG3CoqrYCh/p5SdKUzBvgVXWmqr7QT3+P7o70m4DdwP5+tf3AjZMqpCTpYgsaA09yJXA1cATYUFVn+qeeATaMtWSSpBe1dtQVk7wK+CTwvqr6bpLnnquq6u44P/R1e4G9Sy2oJOmFRtoDT3IJXXjfV1UP9YufTbKxf34jcHbYa6tqX1XtqKod4yiwJKkzylkoAe4GTlTVhweeOgDs6af3AA+Pv3iSpLmkaujIx/MrJG8F/h34MvDjfvHv042DfwL4aeBrwK9W1bfmea+CNUstsyS9xFw4NmwUY94AHycDXJIWY3iAeyWmJDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNWqUe2Lek+RskscGll2e5GCSk/3j+skWU5I02yh74PcCu2Ytuw04VFVbgUP9vCRpiuYN8Kr6DDD7ZsW7gf399H7gxjGXS5I0j8WOgW+oqjP99DPAhjGVR5I0orVLfYOqqu5u88Ml2QvsXep2JEkvtNg98GeTbAToH8/OtWJV7auqHVW1Y5HbkiQNsdgAPwDs6af3AA+PpziSpFGNchrh/cB/AK9PcjrJLcAfA29PchK4vp+XJE1RquYcvh7/xpKCNRN7/03rfumiZU+fOzyx7b0U2caTZxtP1rD2hZXexheODRuG9kpMSWqUAS5JjWpyCGXwK9A/bL/quentv33xex+/8wIAv3L8uV8CWOFflVYG23jybOPJG9bGw9oXVnobO4QiSauKAS5JjVrylZjLYfDr5ps+vfu56QtD1n3TTf1rrn9+2TX/tlK+Fq1ctvHk2caTN6yNh7UvtNnG7oFLUqOa2gOfOSAxeBBirv+msw2+ZtPx5w9srJyDFCuDbTx5tvHkvVTa2D1wSWqUAS5JjTLAJalRBrgkNcoAl6RGNXUWysxR4ON3DpzbedNor525THbwfXQx23jybOPJe6m0sXvgktSopvbAZwz+2MzgVVOj/giQ5mcbT55tPHnD2nghP2a10rkHLkmNMsAlqVFL+j3wJLuAj9L9yPddVfWi98b0lmrts40nzzaeLG+pBiRZA/wFcAOwDfi1JNsWX0BJ0kIseg88ybXAH1TVO/r52wGq6o9e5DUT3QOXpNVp/Hfk2QR8fWD+dL/sBZLsTXI0ydElbEuSNMvETyOsqn3APpjZA5ckjcNSAvxp4IqB+c39shfzTbjw/e5x1XgNq6s+sPrqZH1WvtVWp3HX52eGLVzKGPha4D+Bt9EF9+eBX6+qx+d53dFhYzmtWm31gdVXJ+uz8q22Ok2rPoveA6+q80l+B/hnuiOT98wX3pKk8VnSGHhVfQr41JjKIklagOW4EnPfMmxzklZbfWD11cn6rHyrrU5Tqc+SrsSUJC0ffwtFkho11QBPsivJk0lOJbltmtsehyRXJDmc5Ikkjye5tV9+eZKDSU72j+uXu6wLkWRNki8meaSf35LkSN9PH09y6XKXcSGSrEvyYJKvJDmR5NqW+yjJ7/Wft8eS3J/kZS31UZJ7kpxN8tjAsqH9kc6f9/V6NMk1y1fyuc1Rpz/pP3OPJvm7JOsGnru9r9OTSd4xrnJMLcBXyW+nnAfeX1XbgJ3Ae/s63AYcqqqtwKF+viW3AicG5j8EfKSqXgd8G7hlWUq1eB8F/qmq3gC8ma5uTfZRkk3A7wI7quoqujO+bqatProX2DVr2Vz9cQOwtf/bC3xsSmVcqHu5uE4Hgauq6ufoTrG+HaDPiJuBN/av+cs+D5dsmnvgbwFOVdV/VdUPgQeA3VPc/pJV1Zmq+kI//T26YNhEV4/9/Wr7gRuXp4QLl2Qz8MvAXf18gOuAB/tVWqvPTwK/ANwNUFU/rKpzNNxHdGeLvby/9uIVwBka6qOq+gzwrVmL5+qP3cBfV+ezwLokG6dT0tENq1NV/UtVne9nP0t3cSN0dXqgqn5QVV8FTtHl4ZJNM8BH+u2UViS5ErgaOAJsqKoz/VPPABuWqViL8WfAB4Af9/OvBs4NfBBb66ctwDeAv+qHhe5K8koa7aOqehr4U+B/6IL7O8Ax2u4jmLs/VktO/Bbwj/30xOrkQcxFSPIq4JPA+6rqu4PPVXdaTxOn9iR5F3C2qo4td1nGaC1wDfCxqroa+D6zhksa66P1dHtwW4CfAl7JxV/dm9ZSf4wiyR10w633TXpb0wzwxfx2yoqT5BK68L6vqh7qFz878zWvfzy7XOVboJ8H3p3kv+mGtK6jGz9e139dh/b66TRwuqqO9PMP0gV6q310PfDVqvpGVf0IeIiu31ruI5i7P5rOiSS/CbwLeE89f472xOo0zQD/PLC1P3p+Kd2g/oEpbn/J+vHhu4ETVfXhgacOAHv66T3Aw9Mu22JU1e1VtbmqrqTrj3+tqvcAh4Gb+tWaqQ9AVT0DfD3J6/tFbwOeoNE+ohs62ZnkFf3nb6Y+zfZRb67+OAD8Rn82yk7gOwNDLStaujuUfQB4d1X938BTB4Cbk1yWZAvdAdrPjWWjVTW1P+CddEdnnwLumOa2x1T+t9J91XsUON7/vZNu3PgQcBL4NHD5cpd1EXX7ReCRfvpn+w/YKeBvgcuWu3wLrMt24GjfT38PrG+5j4APAl8BHgP+BrispT4C7qcbv/8R3TekW+bqDyB0Z6s9BXyZ7uybZa/DiHU6RTfWPZMNdw6sf0dfpyeBG8ZVDq/ElKRGeRBTkhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1Kj/B+WGLvmcG7MHAAAAAElFTkSuQmCC\n",
|
||||
"text/plain": [
|
||||
"<Figure size 432x288 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"world.step()\n",
|
||||
"pylab.imshow(np.concatenate(session.run(fluid.density).data[...,0], axis=1), origin='lower', cmap='magma')"
|
||||
@@ -271,38 +547,11 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "b9xHtdDQRrjL"
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Building graph for time step 0\n",
|
||||
"Building graph for time step 1\n",
|
||||
"Building graph for time step 2\n",
|
||||
"Building graph for time step 3\n",
|
||||
"Building graph for time step 4\n",
|
||||
"Building graph for time step 5\n",
|
||||
"Building graph for time step 6\n",
|
||||
"Building graph for time step 7\n",
|
||||
"Building graph for time step 8\n",
|
||||
"Building graph for time step 9\n",
|
||||
"Building graph for time step 10\n",
|
||||
"Building graph for time step 11\n",
|
||||
"Building graph for time step 12\n",
|
||||
"Building graph for time step 13\n",
|
||||
"Building graph for time step 14\n",
|
||||
"Building graph for time step 15\n",
|
||||
"Building graph for time step 16\n",
|
||||
"Building graph for time step 17\n",
|
||||
"Building graph for time step 18\n",
|
||||
"Building graph for time step 19\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"for step in range(20):\n",
|
||||
" print('Building graph for time step %d' % step)\n",
|
||||
@@ -324,41 +573,11 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "TA6Ibs-mXsTc"
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Computing frames...\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"<matplotlib.image.AxesImage at 0x147e73190>"
|
||||
]
|
||||
},
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAACICAYAAAD+r7D/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2deZRc1ZGnf5FLZS0q1aYFbWgBJCwQW9MgG7wgzDTGdANj4zbupvGYczjjNjb08YwN4+7p6Tl9uqFn7AbGtG2a1TQGbDbJgMxuA2azhJAQktBWWqpKS21ZqjWzMjPmj3jvxS1VlmqvyleK7xwd3Yr3Xr5773t5MyJu3LjEzDAMwzDCR2SyK2AYhmGMDBvADcMwQooN4IZhGCHFBnDDMIyQYgO4YRhGSLEB3DAMI6TEJvJmRGQxi4ZhGMOniZlnHi2c0AFciE78LQ3DMEJNdm8+qblQDMMwQooN4IZhGCHFBnDDMIyQYgO4YRhGSLEB3DAMI6TYAG4YhhFSbAA3DMMIKTaAG4ZhhBQbwA3DMELKJKzENMIMUVGf/wGAOeP93zMpdZoKxKKVQXlW+RkAgAjigexIuk7+79ruXGWZKYbD9NJlAID5iXPyHq9LvQ8AONL18YTVabSYBm4YhhFSpoAGTkFpadWXAADnFy0LZCVR+Y2q7ewKZK93PRKUU72HvJJpMy6RyLSg/Jczbw7K/+uPDgIAZi/pCGRt9QkAwNM7Tgxkt9e9G5TrjrwJAMjl9Bkcz8RjMwAA19T810D248t3BeWyG8+TQkKtHLy7GQCw/oFLA9G3N6aD8vrOxwEAmUzLmNc3LLhW4ZLKLwAA7jj51EB22Te8vrnkj/N/wEvyzj9//5WB6Oad24Ly7uRaAABzGoWCaeCGYRghxQZwwzCMkELME+c6kHzgo0knK+4S31UCAC9dUBGUF/yduE5yy052LpHfKGo8rLIX3guKT98v5uxNu9Tkr0++7pWyo6hreCDHk7aq4jsAgBf+3jHF//xiPTfdK4VcTo/H5HqO66QbNTYG5eZ/XA8AuOhlveaj1l94pePDdXVW5TeC8ns3pQAAkW9eHsios1NPzuZ57yLy7nNJSd7Pb7rpJQDAeS8fDGT7ki+OuL5hwZ38fWDFt4Ly155cDADgInWrUHu7FDIDfK9jMjZxeblek1Z3yS++VAsA+C8f3h3IMtnkCGs+XLLrmfnco6WmgRuGYYSUUGng18/5WwDAz3au1M9satITmpr7X+T9qsLRXHiaTtDxrNkAgMjDTwayS28uBQC81PYj54Omrqa4ovLaoLzx0QUAAF44L5DRwcZ+1wT9CgCe5s3FCZWVlQZFniFWTu7OpwLZ6XeJFrq99WnnQ6eKxSPa8gkVFwSS/bcuCcqRa1bJWQfVKuRiZ8LSfz9dTdy3fLIZlTkWD1dXAwDav/frQHbmGvluTEVN3Ne8/235jYHs+l/pO+v3F3U7oa2+FRPLH7tBGa9vc/pd55JiPaFI+vu+q+sD0V9v+TGAidDETQM3DMOYUtgAbhiGEVIGdaEQUTGA1wEkIHHjTzDz3xPRYgCPAagBsB7AtTxIgOTQXSga231Z5X8Pys8+601YuvGxBx23SZFnGpU4pnxCylyk5iaiTh1KxUTi6hq9e+0eAMANl7cGsvsP/ktQLqQ40JFyZuV1QXnDMyfpAW9ykvfncZs4UNzpQ7+/pznmptvH5WXy0XNPCESRTVsBAOdco5NuHyTvH1LdC52SxHwAQNutnw1k0as09ph2eyZ4kWPKuxNrFZ4LxTX/fbPefbddPBdMbtHCQNT0nRcAAIt//Wog607VDakNhc53F/5PAMBt2z8VyKJva3AC2r01BxEdS5Dyvrc0gN7K3iS7O7447hSUi1sw+8nzAtEtS98CAPxw7/8eRu1HwshdKCkAq5j5TABnAbiUiFYCuB3AvzLzyQBaAVw/ltU1DMMwjs2gAzgL/rK7uPePAawC8IQnfwjAlXkuNwzDMMaJIS2lJ6IoxE1yMoC7AewCkGQ/ixFQB2DeAJcPm0ikLCg/etV+PZAR0zL75od67sxS57hnApU4sZ/lEn1CVRrbiTInltY3U3vULcIni0vhJ998IpA9cbu6Gdo6tw6xJYVHJCL99bMVGj/rRpTk1kmsK3c50Q4x/Z2PlEvfcpma8oGR2ubEMrt48khK+zj3x2cDAO5Y/nYgW/WORgflcrpUP2xcNV3cU9Evz1Hh7zcFRfbMcu7uDWSRRerCC9wpEUe/inm93HJEZWn3GckzjHR2B6KZt10IALjsNf1qPpn656E3pMCIRqcH5du/vQ8AEHndiRjZfUhP9l0fjguFZnku2NgAblyvO7nOccu6LpSIrI2I9upzu/3bMj7d8X2tWzbrPKNxZkiTmMycZeazAMwHcB6AUwe5JICIbiCidUS0boR1NAzDMPIwrGRWzJwkotcAfBJAJRHFPC18PoD6Aa65B8A9gD+JOTgrputKy/K/Pj0oZ9ZulM/s1gmfdJ3GX8Zr5Jc1Otf5ME8D7zNh5GgpvrZD7Y726P3CRm68KhBd87PFQfmnnf/klcIXt7yoQlZVnnf77ECWfVUtip790vZMt/62J6pV44hVS5sjKUf765SYblfbcbVL7vGuSapWHfHO/cw/qLW1+Cu64nNX6+qhNqkgKIrPCsp3fE60Q/xONcLUDif51z6xXmIxXZlatk+/PsWXeKsIG3Q1LLdLH3Nar8m1aR/7Gnpsrt6HPEvzhys13nzNCzOCcm/GWUNR0Mj3+spKXWlJl4pVkX1WJy57G9TC6zos3/eyefodjTXK9z4yTy09l1y99F2mVd/tznrV1ktniTzeqpPLkctlQvPKf1Zr68lmDXgY7zFiUA2ciGYSUaVXLgFwCYCtAF4D8GXvtOsAhOsbZxiGEXKGooHPAfCQ5wePAPglMz9LRFsAPEZE/whgA4D7xrGehmEYxlEMOoAz8yYAZ+eR74b4w8ecv5qjeaWzr20Oyi2bxGDYfVDNwNK4mjsLF4jJWV7hxGk3eRMKLe0qSzgx4Sk/OZN6d8ifKKrURFn/tErzNd/7qMjDk3tZXRsXFZ8JAOCNtYEsqfNr2LRfYpjbe/XVuPjsvUE5Vu0Vomq85ZrFpOzYruZiqkuv55zcv3KhxtUXpb1zP6+v0BfLtZ53ta7xrx64WQXE/GnnB+Wav5OY76673wpk+3dVBeXn6+X9LYtq2/5Tu7pbFn/R67uMukvS+6SP67bqZFndkZlBOeZ5J8/79AGVdcvuPSdeqxPWc9/SUOK9rb8ZvGEFQCQirqCvL9H+4BdlSi1dp9/1376n48auTnFTLa1Vd8cFZ0gMfNkS5/vv0LNXxpLfb5ofyLZ36NqGk8rEjfW58/YFsoRXj68vUZfk060aJDHek/G2EtMwDCOk2ABuGIYRUgpqSzVZtQ9ctVjNwKZ39Phva8W0+dlu3Zrrb5ZqFMPSaplVd2NquV1MqNR6je3ctbk6KC84UdwlXW1qVs36rJhlkW07A1nlNzUaZvavpX71ydeG1rBJR0316xZL33V9oC6lrXUatvOLPfIMZhTrb/vnep1YWi/agZzY8FSDuKG27tal8u+2aNx9a1qu+VS9Rv+sqmgAAMR27glkX1uodbq7Xq6fyJja0VDJakJnHpXY9tYGXaPw1D51d6w+JO6SBVF9D0+v0MiIRV70Sc6JdkjuFZfAM3u1j986rO6Dupy833dFtB5nny/3Kfq8ZkI8nXWNxF6Ew4WSiIsL6KI/Ulde10fyTtXtVdfUHR+rC6+RJKpnZZnGwC+tkbGiNJffLddYL8efqVO3yTudGh00k+Vei6frPed3J/vVLbFJXVbdKXOhGIZhGHkoKA28NCGxlDWzNSZ7+y7VXB7fK5MY7/euDWTp3H8Oyr1t8qucONAWyJJvihbzxo5Fgeznu1WjnLFNtPVZJSr7bxWieZfHnV18Pq2aTzXEEsgb+F6AuBsUL50rVsrhfaohv9Gkx5/rksRHp3Xp5GJPh74mJd0y0UOlOnn8/kZ5bj/Zrhrn6706gdeREU1wQ/MXA9mnTpdnWQydTG3oXhCUdZFvODizWN+P1GGZ8P64UTXsdU2qLX/UI5pvQ4lqxluOaO7wlQdEm840ax9sqJe++eGBlwJZskv7rqJEklg9VfcngWzJCaLJVz6iq123U2rojSoQulNi8SZOVCt5/5tikbx6QAMa6iNqMe/u/B0A4DB0d64LG2Rl6omHG/Le560GsUTXeNcCQEuXfma67LPePfUzLysS63PBhdqvfn0nAtPADcMwQooN4IZhGCGloFwokYhUJ1amkwwRZ/V9lMTNEY1oIqW6bifemOV4rk3Nmee3iml5yx416XXTYqBy2icAAIscl8FVteJWWV6jS43jf6kxprWpO4feqAJgeumifrLmTnV3tLph8x0Sd/9GZHsgq2/886BctUzM+2yj9vFqb8LxsWbdgi4a0YmgpeVi1lc5Odk3bhOXwzkJNTdPma4TlmFxofgbQi+tUF3InxBP9mp7y+N6PB6Tvnffww/LL9QP9SbZMs787doD8s4fatOl49NL1ZT/dFy2aStzvtH1jbJeoeoUfY+Xsl6zI1gfUNix9hVlSwEAh97RxnWlZDI261R9RVTTXdTHPwAAHGjTKIgXD3wGAHB1a/58/i8eiPS7ZnqpJrHzP9+9p1+PQ+/oBKpfX2D8E9+ZBm4YhhFSCkoDPzEhWnBKF+yhokQ1vZO98J3XmhxtJ6OTj2WnSnOoSI9/f88bAIAGR9uZXrosKF81TdKYzyvTa5p6JNytz8YdD60JiqVxmTjp6N6FMBCPlBzz+CLH4jmjXNLbvJ98IJClsprQJ1IhWmXb+6pxvJyUFW7uqrNqRwu5okq0mHOqVPOpTMhzjVXo83vhgIZ/RjwNPpst7N2PIlHRpt2NX9IpeQ9Lo9pHC6epNj4jLX3jamdtTpIq8lL25jJqhazLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 432x288 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"print('Computing frames...')\n",
|
||||
"original_state = np.concatenate(session.run(fluid.density).data[...,0], axis=1)\n",
|
||||
@@ -379,44 +598,11 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "7KPpyIwjYETi"
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"WARNING:tensorflow:From /usr/local/lib/python3.7/site-packages/tensorflow_core/python/ops/math_grad.py:1375: where (from tensorflow.python.ops.array_ops) is deprecated and will be removed in a future version.\n",
|
||||
"Instructions for updating:\n",
|
||||
"Use tf.where in 2.0, which has the same broadcast rule as np.where\n",
|
||||
"Initial loss: 65.417114\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"<matplotlib.image.AxesImage at 0x14b36c750>"
|
||||
]
|
||||
},
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAACICAYAAAD+r7D/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2deZRc1ZGnf5FLZS0q1aYFbWgBJCwQW9MgG7wgzDTGdANj4zbupvGYczjjNjb08YwN4+7p6Tl9uqFn7AbGtG2a1TQGbDbJgMxuA2azhJAQktBWWqpKS21ZqjWzMjPmj3jvxS1VlmqvyleK7xwd3Yr3Xr5773t5MyJu3LjEzDAMwzDCR2SyK2AYhmGMDBvADcMwQooN4IZhGCHFBnDDMIyQYgO4YRhGSLEB3DAMI6TEJvJmRGQxi4ZhGMOniZlnHi2c0AFciE78LQ3DMEJNdm8+qblQDMMwQooN4IZhGCHFBnDDMIyQYgO4YRhGSLEB3DAMI6TYAG4YhhFSbAA3DMMIKTaAG4ZhhBQbwA3DMELKJKzENMIMUVGf/wGAOeP93zMpdZoKxKKVQXlW+RkAgAjigexIuk7+79ruXGWZKYbD9NJlAID5iXPyHq9LvQ8AONL18YTVabSYBm4YhhFSpoAGTkFpadWXAADnFy0LZCVR+Y2q7ewKZK93PRKUU72HvJJpMy6RyLSg/Jczbw7K/+uPDgIAZi/pCGRt9QkAwNM7Tgxkt9e9G5TrjrwJAMjl9Bkcz8RjMwAA19T810D248t3BeWyG8+TQkKtHLy7GQCw/oFLA9G3N6aD8vrOxwEAmUzLmNc3LLhW4ZLKLwAA7jj51EB22Te8vrnkj/N/wEvyzj9//5WB6Oad24Ly7uRaAABzGoWCaeCGYRghxQZwwzCMkELME+c6kHzgo0knK+4S31UCAC9dUBGUF/yduE5yy052LpHfKGo8rLIX3guKT98v5uxNu9Tkr0++7pWyo6hreCDHk7aq4jsAgBf+3jHF//xiPTfdK4VcTo/H5HqO66QbNTYG5eZ/XA8AuOhlveaj1l94pePDdXVW5TeC8ns3pQAAkW9eHsios1NPzuZ57yLy7nNJSd7Pb7rpJQDAeS8fDGT7ki+OuL5hwZ38fWDFt4Ly155cDADgInWrUHu7FDIDfK9jMjZxeblek1Z3yS++VAsA+C8f3h3IMtnkCGs+XLLrmfnco6WmgRuGYYSUUGng18/5WwDAz3au1M9satITmpr7X+T9qsLRXHiaTtDxrNkAgMjDTwayS28uBQC81PYj54Omrqa4ovLaoLzx0QUAAF44L5DRwcZ+1wT9CgCe5s3FCZWVlQZFniFWTu7OpwLZ6XeJFrq99WnnQ6eKxSPa8gkVFwSS/bcuCcqRa1bJWQfVKuRiZ8LSfz9dTdy3fLIZlTkWD1dXAwDav/frQHbmGvluTEVN3Ne8/235jYHs+l/pO+v3F3U7oa2+FRPLH7tBGa9vc/pd55JiPaFI+vu+q+sD0V9v+TGAidDETQM3DMOYUtgAbhiGEVIGdaEQUTGA1wEkIHHjTzDz3xPRYgCPAagBsB7AtTxIgOTQXSga231Z5X8Pys8+601YuvGxBx23SZFnGpU4pnxCylyk5iaiTh1KxUTi6hq9e+0eAMANl7cGsvsP/ktQLqQ40JFyZuV1QXnDMyfpAW9ykvfncZs4UNzpQ7+/pznmptvH5WXy0XNPCESRTVsBAOdco5NuHyTvH1LdC52SxHwAQNutnw1k0as09ph2eyZ4kWPKuxNrFZ4LxTX/fbPefbddPBdMbtHCQNT0nRcAAIt//Wog607VDakNhc53F/5PAMBt2z8VyKJva3AC2r01BxEdS5Dyvrc0gN7K3iS7O7447hSUi1sw+8nzAtEtS98CAPxw7/8eRu1HwshdKCkAq5j5TABnAbiUiFYCuB3AvzLzyQBaAVw/ltU1DMMwjs2gAzgL/rK7uPePAawC8IQnfwjAlXkuNwzDMMaJIS2lJ6IoxE1yMoC7AewCkGQ/ixFQB2DeAJcPm0ikLCg/etV+PZAR0zL75od67sxS57hnApU4sZ/lEn1CVRrbiTInltY3U3vULcIni0vhJ998IpA9cbu6Gdo6tw6xJYVHJCL99bMVGj/rRpTk1kmsK3c50Q4x/Z2PlEvfcpma8oGR2ubEMrt48khK+zj3x2cDAO5Y/nYgW/WORgflcrpUP2xcNV3cU9Evz1Hh7zcFRfbMcu7uDWSRRerCC9wpEUe/inm93HJEZWn3GckzjHR2B6KZt10IALjsNf1qPpn656E3pMCIRqcH5du/vQ8AEHndiRjZfUhP9l0fjguFZnku2NgAblyvO7nOccu6LpSIrI2I9upzu/3bMj7d8X2tWzbrPKNxZkiTmMycZeazAMwHcB6AUwe5JICIbiCidUS0boR1NAzDMPIwrGRWzJwkotcAfBJAJRHFPC18PoD6Aa65B8A9gD+JOTgrputKy/K/Pj0oZ9ZulM/s1gmfdJ3GX8Zr5Jc1Otf5ME8D7zNh5GgpvrZD7Y726P3CRm68KhBd87PFQfmnnf/klcIXt7yoQlZVnnf77ECWfVUtip790vZMt/62J6pV44hVS5sjKUf765SYblfbcbVL7vGuSapWHfHO/cw/qLW1+Cu64nNX6+qhNqkgKIrPCsp3fE60Q/xONcLUDif51z6xXmIxXZlatk+/PsWXeKsIG3Q1LLdLH3Nar8m1aR/7Gnpsrt6HPEvzhys13nzNCzOCcm/GWUNR0Mj3+spKXWlJl4pVkX1WJy57G9TC6zos3/eyefodjTXK9z4yTy09l1y99F2mVd/tznrV1ktniTzeqpPLkctlQvPKf1Zr68lmDXgY7zFiUA2ciGYSUaVXLgFwCYCtAF4D8GXvtOsAhOsbZxiGEXKGooHPAfCQ5wePAPglMz9LRFsAPEZE/whgA4D7xrGehmEYxlEMOoAz8yYAZ+eR74b4w8ecv5qjeaWzr20Oyi2bxGDYfVDNwNK4mjsLF4jJWV7hxGk3eRMKLe0qSzgx4Sk/OZN6d8ifKKrURFn/tErzNd/7qMjDk3tZXRsXFZ8JAOCNtYEsqfNr2LRfYpjbe/XVuPjsvUE5Vu0Vomq85ZrFpOzYruZiqkuv55zcv3KhxtUXpb1zP6+v0BfLtZ53ta7xrx64WQXE/GnnB+Wav5OY76673wpk+3dVBeXn6+X9LYtq2/5Tu7pbFn/R67uMukvS+6SP67bqZFndkZlBOeZ5J8/79AGVdcvuPSdeqxPWc9/SUOK9rb8ZvGEFQCQirqCvL9H+4BdlSi1dp9/1376n48auTnFTLa1Vd8cFZ0gMfNkS5/vv0LNXxpLfb5ofyLZ36NqGk8rEjfW58/YFsoRXj68vUZfk060aJDHek/G2EtMwDCOk2ABuGIYRUgpqSzVZtQ9ctVjNwKZ39Phva8W0+dlu3Zrrb5ZqFMPSaplVd2NquV1MqNR6je3ctbk6KC84UdwlXW1qVs36rJhlkW07A1nlNzUaZvavpX71ydeG1rBJR0316xZL33V9oC6lrXUatvOLPfIMZhTrb/vnep1YWi/agZzY8FSDuKG27tal8u+2aNx9a1qu+VS9Rv+sqmgAAMR27glkX1uodbq7Xq6fyJja0VDJakJnHpXY9tYGXaPw1D51d6w+JO6SBVF9D0+v0MiIRV70Sc6JdkjuFZfAM3u1j986rO6Dupy833dFtB5nny/3Kfq8ZkI8nXWNxF6Ew4WSiIsL6KI/Ulde10fyTtXtVdfUHR+rC6+RJKpnZZnGwC+tkbGiNJffLddYL8efqVO3yTudGh00k+Vei6frPed3J/vVLbFJXVbdKXOhGIZhGHkoKA28NCGxlDWzNSZ7+y7VXB7fK5MY7/euDWTp3H8Oyr1t8qucONAWyJJvihbzxo5Fgeznu1WjnLFNtPVZJSr7bxWieZfHnV18Pq2aTzXEEsgb+F6AuBsUL50rVsrhfaohv9Gkx5/rksRHp3Xp5GJPh74mJd0y0UOlOnn8/kZ5bj/Zrhrn6706gdeREU1wQ/MXA9mnTpdnWQydTG3oXhCUdZFvODizWN+P1GGZ8P64UTXsdU2qLX/UI5pvQ4lqxluOaO7wlQdEm840ax9sqJe++eGBlwJZskv7rqJEklg9VfcngWzJCaLJVz6iq123U2rojSoQulNi8SZOVCt5/5tikbx6QAMa6iNqMe/u/B0A4DB0d64LG2Rl6omHG/Le560GsUTXeNcCQEuXfma67LPePfUzLysS63PBhdqvfn0nAtPADcMwQooN4IZhGCGloFwokYhUJ1amkwwRZ/V9lMTNEY1oIqW6bifemOV4rk3Nmee3iml5yx416XXTYqBy2icAAIscl8FVteJWWV6jS43jf6kxprWpO4feqAJgeumifrLmTnV3tLph8x0Sd/9GZHsgq2/886BctUzM+2yj9vFqb8LxsWbdgi4a0YmgpeVi1lc5Odk3bhOXwzkJNTdPma4TlmFxofgbQi+tUF3InxBP9mp7y+N6PB6Tvnffww/LL9QP9SbZMs787doD8s4fatOl49NL1ZT/dFy2aStzvtH1jbJeoeoUfY+Xsl6zI1gfUNix9hVlSwEAh97RxnWlZDI261R9RVTTXdTHPwAAHGjTKIgXD3wGAHB1a/58/i8eiPS7ZnqpJrHzP9+9p1+PQ+/oBKpfX2D8E9+ZBm4YhhFSCkoDPzEhWnBKF+yhokQ1vZO98J3XmhxtJ6OTj2WnSnOoSI9/f88bAIAGR9uZXrosKF81TdKYzyvTa5p6JNytz8YdD60JiqVxmTjp6N6FMBCPlBzz+CLH4jmjXNLbvJ98IJClsprQJ1IhWmXb+6pxvJyUFW7uqrNqRwu5okq0mHOqVPOpTMhzjVXo83vhgIZ/RjwNPpst7N2PIlHRpt2NX9IpeQ9Lo9pHC6epNj4jLX3jamdtTpIq8lL25jJqhazLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 432x288 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"target = session.run(fluid.density).data[0,...]\n",
|
||||
"loss = math.l2_loss(fluid.density.data[1:,...] - target)\n",
|
||||
@@ -437,30 +623,11 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "pvvF6xqmaRLX"
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Optimization step 0, loss 65.417122 . The first step sets up the adjoint graph.\n",
|
||||
"Optimization step 1, loss 61.038647 \n",
|
||||
"Optimization step 2, loss 55.527832 \n",
|
||||
"Optimization step 10, loss 43.914444 \n",
|
||||
"Optimization step 20, loss 40.739956 \n",
|
||||
"Optimization step 30, loss 39.523136 \n",
|
||||
"Optimization step 40, loss 39.361752 \n",
|
||||
"Optimization step 50, loss 39.348873 \n",
|
||||
"Optimization step 60, loss 37.156372 \n",
|
||||
"Optimization step 70, loss 37.729256 \n",
|
||||
"Optimization step 80, loss 35.622345 \n",
|
||||
"Optimization step 90, loss 35.464912 \n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"for optim_step in range(100):\n",
|
||||
" _, loss_value = session.run([optim, loss])\n",
|
||||
@@ -477,34 +644,11 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "SQBtCmhZaYYj"
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"<matplotlib.image.AxesImage at 0x14b4ad890>"
|
||||
]
|
||||
},
|
||||
"execution_count": 13,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2dd3icxbX/v2eLqtVsy3JvYJtmOgSSEBMICRAIkNxLICSYhPxIbhpcckMgyb2X3JQbkpsEUgmhhlBDdSAUY6oh2IBxt8FVtmw1S1p1abW78/vjzL5nJK2klayyqz2f59Gj2fO2eeed9505Z86cIWMMFEVRlMzDN9YZUBRFUcYGbQAURVEyFG0AFEVRMhRtABRFUTIUbQAURVEyFG0AFEVRMpRRawCI6HtEdPtw75vEuQwRHToc51IGhoiyiWgzEU0boXNvJaLS4T53kte/gohWjtK1hvwOENFlRPT8MOdnrn2XAsN53iHko4WI5o/Quf+XiK6x6dOJqGIkrjMSENFqIjpysMcNqQGwL8IGImojoioi+iMRFfd3jDHmp8aYLydz/sHsOxIM9kUf6ZeDiO4moh8PsM+QGjoiOoWIlhNRPRHVEtHf3I83MTcRUZ39u4mIqJ9TXgXgVWNMZbJ5TxZjTCeAOwFcPxznSxUSfWwO5h0wxtxnjPn48OQutTDGTDDG7Bzu89pOxeUA/jQM5+r1/RjO96AP/g/A/wz2oEE3AET0bQA3AfgOgCIApwCYA2A5EWX1ccyY9hqUfikBcBuAueDn2AzgLmf7VQAuBHAMgKMBnA/gK/2c76sA7h2JjFruB7CUiLJH8BrKGDGG34orAPzDGNM+Rtc/WJYB+CgRTR3UUcaYpP8AFAJoAXBxD/kEALUAvmR/3wjgEQB/BdAE4MtW9lfnmMsBlAOoA/CfAHYD+Jhz/F9tei4AA2ApgD0ADgD4vnOekwH8E0AIQCWA3wHIcrYbAIf2cT9XANgJ/ujtAnAZgMMBdACI2nsN2X0/CeBdez97AdzonGePvU6L/TvVyr8EYAuABgDPAZjTT9n+DUAVgEYArwI40sqvAtAFIGzP/fcEx75qr99q9/nsYJ5rj3MdD6DZ+f0GgKuc31cCeLOPY2cDaAcQ6C/vtoxfts9sE4BPOee4G8CtAJbb5/JKz3IDsA3AkiTv54v2GTTbZ/0VZ9vpACoAfBtAja0/X3S2TwK/WE0AVgP4EYCV/VzrU/Z+Qvb+Dne27QZwA4DNtj7cBSAHQL4ts5hTf6Yj8TvwRVv3GsAN7UkA1tvr/a5HvV5p09c5522xz+Nuu60IwB32vvcB+DEAv93mB/cqD9hy+7rNQyDBfX8XwCM9ZLcA+E0S17kCwOsAfg3+FvwYwKH2uTfa6z+U6H225/0L+NtTDuAHAHxuGdh7aAC/3+f08+xeBPD5BHXjezYPuwFc5mxPeG0k+H5g6O/BHwA8Y495HcBUADfb+9kK4Lge97AcwNJBve+D/DicDSDSRyW4B8ADNn2jveELbaHkonuFPsLe1IcBZNmH1IX+G4A/2/McA6AT9uUCcAJYCwnYfbcAuCZRhemR33zwi73I/p4G+ehegR4vuq0Qi+39HA2gGsCFPfIYcPa/AMB2+5ADtoK80U/ZfglAAYBs+5DX9qgMPx7g2XS7T/DHONTP3+f6OM81cD7w4JfwA87vE+E0ED2O/SSATT1k3fIOIGjL5Xv22Z8B/jgvcvZvBvARWxa3JHgWywB8K8k6+0kAhwAgAEsAtAE43nmmEbDqHARwrt1eYrc/COBhW1eOAn+8EjYAABaCG+Cz7Lmus/eZZbfvBrARwCwAE8Ev9I/dj02P892I3u/AreBG4+Pgj8wTAKYAmAFuwJb0VX+tfBaA/bAfQgCPg00e+fY8q2EbSHADs9XJ70vouwGYY8utwP72gz/2pyRxnSvsM/gm+D3JBfAAgO+D37UcAB9OVM/BH+Anwe/NXADvA7jSOW8XgP9n8/Nv9t6pj+dXC+CkHu97BMCvwPVwiX2+i5K8ds86ezcG/x4cAH/fcsAN1C5wx9kPbihf6nGN3wD4lfM75JZdwvtO5iVyTvh5AFV9bPsZgOVO5X21nwr9X7CNhf2dB24d+2sAZjr7rwZwST8fsMf7+jA68nxbQJ8BkNtjW68HmOD4mwH8ukce3QbgmXiFsL994JdkThLlXGzPV5So8vRxTML7HOTzPRpAPYDTHFkUwGHO7wX2Wr1eJLAG9WYPWc+KfxpY0/E5sgdgNSq7/4POtgk2D7Mc2X0A/muI9/gEgKtt+nQ4GouV1YA7FH7wB8S995/2VS/AWuzDPZ73PgCn29+7AXzV2X4ugB1OPpJpAGY42+vgaHoAHoXt+CSqv+AP6zsAvmt/l4E7UrnOPpfCflTAHxw3vx/vWcd7nH8lgMtt+izn3ga6zhUA9vQ411/AZsmZCa5jwBqCH/zNOMLZ9hUALzvn3e5sy7PHTu0j/z2f9engBiDfkT1sn3My1x6oAUjmPfizs+2bALY4vxfDWicc2U8A3DmY92GwYwAHAEzuw043zW6Ps7ef80x3txtj2sAVuj+qnHQb+MMAIlpIRE/Zwegm8Es6eYBzwRjTCuCz4J5OJRE9TUSH9bU/EX2AiF6yA6WN9rj+rjMHwC1EFCKiEPjDSuDeWs9z+4noZ0S0w97DbrtpwPsYLuwA8jPgj+NrzqYWsOkvTiGAFmNrXA8awD2i/pgOYK8xJubIytG9XNy60QIuu+nO9gJw4z0gRHQOEb1pB7lD4A+vW651xpiI8ztet0rBPVK3Hpf3c6np7nZ7f3v7ui+7r3tPyVDtpNsT/J7Qz7F3AHjPGHOT/T0H3AutdOron8A9dKDHO4r+7x3gsZlLbfpz9ncy1wF6fyuuA78rq4loExF9KcH1JtvzuvnqWY+8b4b9xgB9l1GiuttgvxPu+acnee2BSOY9GOzzTvq9iDPYBuCf4Nb8066QiCYAOAfACkec6AMRpxLATOf4XLC9dSj8EayqLjDGFIJVqv68VCSDxjxnjDkL3HhtBZuZgMR5vx9sephljCkCq+PUz/57wWpusfOXa4x5I8G+nwObjD4Gti3OtfL+zt8vRDTbusz19XeZs+8cAC8A+JExpucA7iaw2S3OMVaWiPUA5vXoIPTM+34As4jIrXuzwb3lOLOcvE0AmyD2O9sPB7Cujzx42IHiR8EmxjJjTDGAfyC5+lEL7gHOcmSz+9l/P/hjF7822WMT3pc9V/yeBv18BwMRXQ82UV3piPeC3+XJTv0sNMbEXQkrE+S3P/4G4HQimgngIkgDMNB1gB73b4ypMsb8P2PMdHDP+g8JPNwOgHvtcxxZz3o0GNaDy8ilhIjye5x/fxLXTvQ8h/IeDJak3guXQTUAxphGAD8E8FsiOpuIgkQ0F6waVSB5749HAJxPRB+0nkM3IsmPdgIKwLb8FtuD/7dkDiKiMiK6wD7gTnBPN94aVwOY2cOrqQBAvTGmg4hOBn+049TaY13/5FsB3BD3zSWiIiL6137uoROsBeWBtRiX6h7nTkS3fYwxewy7zPX1d5/N1wywuv87Y8ytCc77FwDXEtEMIpoOHjC9O1EGjDEVYLvmyf3kfRW4l32drT+ngz2LHnT2OZeIPmzL/0dgs9JeJ78TAbxpf8ddcOcmyFIW2H5bCyBCROeATRkDYoyJAngMwI1ElEdER4AdEfriYQCfJKIziSgILqdO8CB6nK8T0Uwimgi2cT9k5dUAJhFRUTJ5Gwz2nr8F4CLjeLgYdtN9HsAviaiQiHxEdAgRLXHu51s2vyUYwPXWGFMLHtC8C8AuY8yWJK+TKM//ahsSgHvmBvJuxq8XtXn8CREV2E7MtWDHk6HwD7Cdvyc/JKIsIjoNwHkA/pbEtRN9P4byHiQNEeWAxwuWD+a4QbuBGmN+Du5l/x/4w7sK3MqfadhPO5lzbALbtB4E9zRawLbXpI7vwX+AP8bN4B78Q/3v7uEDP7T9YBPDEkjj8SK4l1tFRHGz1tcA/A8RNYPHMB527qcNbH973aq5pxhjHge7yz5ozTobwVpSIv4CVv/2gb1E3uyx/Q4AR9hzP9HHOW4EcI/d5+JkCsDyZXDFvNHVEJztfwLwdwAb7D08jf59pf8E4At95d0YEwZX9HPAPak/gG3HW51j7gfw3+DncgJ47CnO5wDc49S1WZCy64Yxphn88XsY/CH5HFiLS5ZvgNXsKnCjd1dfOxpj3rP5/K29r/MBnG/v172v58FeNTvAA3mw9/4AgJ22nAZrGuqPz4LNWVuc5xtv6C8HN5Jxz6RHwNowwO/Sc+Ae5RpwYzgQ94O12Pt7yPu7TiJOArDK1sNlYLNkIt//b4IHZneCxyDuB88TGQp/AXc8ch1Zlc3vfvC401edetrftRN9P4byHgyG88FjEJ6mbJ/1af0dRIlNuaOLVfNDYDPOrrHOjzJ0rNnlXXCHoHIIx98NHhD9QR/nXgfgI8aYGiv7AYBaY8xBT+AZSYhoN4AvG2NeGOu8KIkhop8CqDHG3DzWeRksRLQK7HSycTDHjdkELSI6HzxmQGBtYgNk8FNJU2zP/IgRPPdhPWQjObtSySCMMd8b6zwMFWPMB4Zy3FgGg7sArFrtB7sWXtKHZ4miKIoyAqSECUhRFEUZfTQctKIoSoaiDYCiKEqGMqqDwESk9iZFUZTBc8AYM+zrYIyBF5B/9C+pKIqS1kQHCsUxJNQEpCiKkqFoA6AoipKhaAOgKIqSoWgDoCiKkqFoA6AoipKhaAOgKIqSoWgDoCiKkqFoA6AoipKhjFk4aCU9iS9y5C52FF9S15iOMcnTeCDgL/bSkyccxTLK9mRN4Qr+3/a+c5ROrB8MhXmLAAAzs49PuL2icw0AoKntvVHL01ijGoCiKEqGMg40AFlKeGHJZwAAJwUXebK8ALdxO1tbPdnKNlmxrrOr2qa0N+Xi803w0peVXu2Line truncated
|
||||
"text/plain": [
|
||||
"<Figure size 432x288 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"optimized_state = np.concatenate(session.run(fluid.density).data[...,0], axis=1)\n",
|
||||
"pylab.title(\"Original state at t=20 (top), and optimized version (bottom):\")\n",
|
||||
@@ -526,34 +670,11 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 15,
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "i7ZahlUudex8"
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"<matplotlib.image.AxesImage at 0x14b5363d0>"
|
||||
]
|
||||
},
|
||||
"execution_count": 15,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAACYCAYAAAD9ebJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2deZAkV3ngf19VV/dcmhlpEDpmdLESh2BBEDJIizEOBBhx27AcBiwWYS0RYNCasC2Md/HuWhHgZTHY4pA4zAhkhCwJJGvBBgsImUswIFlGlxEIgUajczRXz0wfVd/+kS8zv+x6WVlVXX1kz/eL6KisPF6+9zL71Xe974mq4jiO49SPxlJXwHEcxxkOH8Adx3Fqig/gjuM4NcUHcMdxnJriA7jjOE5N8QHccRynpvgAfoggIl8VkXN6HP+EiPz3Psv6loi8dXS1GxwR+XMR+fw8y3iDiHxthHWaEJHbROSYUZU5p/x9IvK4Ia+9VUR+c8T1+ayI/EXYfqqIfHeU5TvV+ABeY0TkFyLy/H7OVdWzVXVruO7NIvLtOcffpqr/eyHquVxR1ctU9YXpdxFRETl5HkWeB9ygqjvmW7fYj6SqrlPVnw9Tnqo+WVW/Nd969Sj/FmCXiLxsoe7hdOMDuOOMjrcBn1vqSiwhlwH/dakrcSjhA/gKIZWqReSDIvKoiNwtImeb498SkbeKyJOATwBnBpV8Vzhu1eHDReQ6EXkolHWdiGzpow5Hi8h+Edlk9j0jlNOKnP9xEfngnH3XiMgfhu1jReSqcP3dIvLOHvd+eTAT7AptfZI5dpyIXB3KeURELrJ9FrZvCKf/a+iX14rIT6xEKSItEXlYRJ4euf/xwOOAG82+DSJyabjvPSLyZyLSMPf+johcJCK7ReQOETkrHLsQeA5wUahLWt9MQwjP62PBNLYvlHW0iHw4PLM7bD2tthb6aF/4mwzlnhiOvVREbg7nfFdEnmrKeLqI/FhE9orIF4FVc7rhW8BZIjJR9pyc0eID+MriWcCdwGOAvwQ+LSJiT1DV20kkxe8FlXxjpJwG8LfACcDxwAHgoqqbq+r9JP/ErzG73wRcrqozkUu+ALw2raOIHA68ELg8DHT/APwrsBk4CzhfRH5rbiEi8vhQ1vnAkcBXgH8QkXERaQLXAfcAJ4ayLo/U/TfC5tNCv3wRuBR4ozntxcAOVb0p0pb/CPxcVWfNvr8BNpAM7M8Ffg/4L+b4s4CfkTyv9wFXi8gRqvpe4F+Ad4S6vCNyP0j6+c/C9VPA94Afh+9XAh+KXaSqG0O564CPhHttDwP+Z0ik6E3AxcC1ktj2x4Evk2gYRwB/D7xqTrnbgRngCSX1dUaMD+Ari3tU9ZOq2ga2AscARw1aiKo+oqpXqep+Vd0LXEgyAPXDVsKgFwbP11NuVvgXQEmkTYBXk/yw3Af8GnCkqv4vVZ0Ott9PAq+LlPNa4P+p6tfDD8UHgdXAfwKeCRwL/JGqTqrqQVX9dqSMGJ8HXiwi68P3N/Voy0Zgb/oltP11wHtUda+q/gL4v6GMlAeBD6vqTPjBuBN4SZ91A/iSqv5IVQ8CXwIOquql4fl/EejSFCwi8lrgd4FXhX47D7hYVW9U1XbwmUwBZ4S/lqnvlcAPI8XuDX3hLAI+gK8s7k83VHV/2Fw3aCEiskZELg5q/x7gBmBjGJSquAY4VUROAl4A7FbVH4RybzWq+3M0yaR2OckgD8lgclnYPgE4Nqjyu4Kp50+J/yAdSyJhA6CqHeBXJNL2cSQ/bLOR63oSfki+A7xKRDYCZ5v6zeVR4DDz/TEkA949Zt89oU4p27WYTe6e0JZ+ecBsH4h8L332Qdq+CPhtVX0o7D4BePecPj8u1OnYkvrO5TBg1wBtcObB2FJXwFkSqlJQvptEDX6Wqt4vIqcBNwHS+zJQ1YMicgWJFP5EjMSqqk+OXPIF4Gsi8n4Sk8Jvh/2/Au5W1VOq7gncR2LCACCYZI4DtpNIkMeLyNgwgziJRvFWkv+V7wUzQYxbgJPMfR4mMSecANwWzjk+1Clls4iIGRSPB64N2wuWJlREHktiDnn7HHPQr4ALVfXCyDXPLanvz8w5m4FxEk3CWQRcAj80eQDYEuyaMQ4jkeB2icgRJPbZQbgUeDPwciqiMsIA8jDwKeCfVDWV3n4A7BWRPxGR1SLSFJGniMivRYq5AniJiJwlibP03SQD93dDOTuA94vIWhFZJSLPLqnOAyT2asuXgWcA7wrtKmvHvcBdJCYbghnjCuBCETlMRE4A/pDELJPyWOCdwTn6n4Enkdjvy+oyb0RkjMQ+/nlVvWLO4U8CbxORZ0nCWhF5iYgcRmJfnzX1/Z20rYbnAt9Q1alR19uJ4wP4ock3gFuB+0Xk4cjxD5PYkB8Gvg/84yCFq+p3gA7wY1WNqdlz+Tvg+eEzLaMNvBQ4DbibfJDfELnfnSQS/9+E814GvCzYztvh+8nAL4F7SWzmMf4c2BrMB68JZR8ArgJOAq6uaMfFFG3cfwBMAj8Hvh3a9xlz/EbglFDnC4FXq+oj4dhHgFeHiJK/rrjvIGwh8Tmcb8xZ+0TkeFXdBvw+iWnlUZIfpDcDqOo08Dvh+06SPpzbH28giXByFgnxBR2chUBEvgH8nap+aqnrMl9E5H8Aj1fVN1acN0FiajqrajKPiLwZeKuq/vrIKrqEhHDDi1X1zKWuy6GE28CdkRPMHM8AXrHUdZkvwYR0LkXJOkowHZy64JVahoSZmD54LzJuQnFGiohsBf4ZOD+EINYWEfl9EsfeV1X1hqrzHWexcROK4zhOTXEJ3HEcp6b4AO44jlNTFtWJ2TxsrY5tOnwxb+k486Ny6lIP3DpZZD592Q8ruL+n79n+sKoeOXf/og7gY5sO5+j3/cFi3tJx5ocP4KNjKQbwUT6Dha5/D375lgui8ynchOI4jlNTPA68ziyERHCoSo1VfSnDdExJoTrnczFZLClynv0pfdazGERnLuqkJ0T2VSEl27EHNkx/jvAZuATuOI5TU1wCd4qk0oEVNkYofC470rY1+m9k39Kh7TituKjqeBX9agjzlfqHkjj7lLbteX3ex16ibXNR2p9W6o7ti0nbtrqF90IiN62o73weax/P1CVwx3GcmuIDuOM4Tk1xE0qdGcqxVsF8VfnFZkRd0K9ZJDm5v5uK0Z8rzSn9tqOsnrEyF+T96PO8iv4s9HdaT7HH8xtJDzFTO3Y7Yu4w/SKzFSaUcB+1ZpOOOaERsZcM08cjfC4ugTuO49QUH8Adx3FqiptQDlHKTAbR5JRVZpWFUNVjLEncdH+xvxI5r9CX7ZjZxEZNDF61qjotmlllRJSZTdL9hT4O22pOtCaUrO+NuSQ1oUg736fmPjrWbcopbFdFnMSumVPffhjEnOcSuOM4Tk1xCfwQIxpzWzghfHYqpMPF8nUulsBYJSGVtFeiDriwoXFnmHZtMIDzuCRWut9ntCTOzkj8dBUFaTwRoxvGuZhu2/UM1LStPR2uNZqPzCTbjdn8Nh07AqanWqm83xewNA68+/q+JWyPA3ccx1m5+ADuOI5TU9yEUmMGil3uWVB3mQXlrUqTkzmfo2CxTCdDxL1LZMq3GPU+6w5TdNvGFkcdmpHt6DTvErvJvNLezjOuOVpmbGc8Lj6Pl+9972az07XdtKYWs33gYDPZZ/q6GcwqjSnT3jX5Ne1mRTX6nfJfqFNJWZFzh8ElcMdxnJriEvihRsWsN43MYOvb8TWIFLgQyffnM4t0SCk0lbAKE/qCdGj7tRNzYhbub8oM5xZmBFY1LSa19xv2thCUaQoZEWncOs4L2kd3A8aCY3N8LI8JbJntqYOtZGN3PsQ19yflpJI4QGfC3EbLa1tKRCqPSt0DOXD7P9UlcMdxnJriA7jjOE5NqTShiMgq4AZgIpx/paq+T0ROAi4HNgE/At6kqtPlJTnLibJZb9lsNjuDzTiCMrU+kgSoVE3sZdoYRF9dxom2JOJMa7Vylb4928y2YwvDiG1b6O+CESJi+oqafSL5nKL7yhiVQzO6Uk4kjzZkNofSxXUayZdOp1veHG/mfTwxlgd4T61LhrZ9D+Y2kta+5LMxk18/c1hJ/WPEnPXZLNHufZXlVN1uRHHgU8DzVPVpwGnAi0TkDOADwF+p6snAo8C5/VXLcRzHGQWVA7gmhN8uWuFPgecBV4b9W4FXLkgNHcdxnCh9RaGISJPETHIy8FHgZ8AuVU11lnuBzQtSQ6ecTG2en2khZk7RssVgU0uAjXuOLEcWT4pVsh1jKc0lFaYJtTqwplO6u+vbMup9ezxX72cPRP7tbMBJu7sauZmqos7LJVlV1JRTYp7JEk/FI3U02CfaxtY32+qWPdeO5RbcY9bvAeDfJ9Zn+1r7klJtFMqBo8w90wRZ9n2Obdt3v9/UFCX0YybpRV9OTFVtq+ppwBbgmcAT+72BiJwnIttEZFt73+SQ1XQcx3HmMlAcuKruEpFvAmcCG0VkLEjhW4DtJddcAlwCMHHilmUiHjiW4szCNOVmSTKrIH1Iszu1Z4FYqtRhUtWW0e+b1G/89LD0iB0eMzMH167Kxb4De1YlGyVpSxsh7akNi041H5u8qTYLR8eIPvd48q/ObLecOT2dDF0zE/mx8Wau5Rw+sR+Auzbm/T42lTg0Vz2aa0Z7tJUX2pjzCVEJvGo1Icu8JOw+nm+lBC4iR4rIxrC9GngBcDvwTeDV4bRzgGuGrafjOI4zOP1I4McAW4MdvAFcoarXichtwOUi8hfATcCnF7CejuM4zhwqB3BVvQV4emT/z0ns4c5SkalY1iHU37R3LXEopXmW29ZxafW0VjgQcVwWc4jHTCgVOuFCGNgWuMzMoWn6q9Ppbue6ialse2eID+/M5v9+ah3JQcMvxOf3MkPB6Mwpo3KCDlOfshWKQn92jFkvDeU+OJGbQDrm+iNaib/t6E2782vkKAAmHjqLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 432x288 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"optimized_velocity_field = session.run(initial_state.velocity).at_centers()\n",
|
||||
"\n",
|
||||
@@ -572,34 +693,11 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "Pqw5BDxmdkut"
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"<matplotlib.image.AxesImage at 0x14b3efa50>"
|
||||
]
|
||||
},
|
||||
"execution_count": 12,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAACYCAYAAAD9ebJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2dabQlR3Ggv3j3vd671WpJtPYFC2NkBhBHBjHYmCNhg9gHMKtBjIU1zIEBBs7YYDxjZmydwTMYIyMMCAS0DCOBhQChAduMgCPEItNsYhECoQV101oa1FKvb435UVm3ot6NunW3t9Tr+M5559bNrMrKzKqXNzIiMlJUlSAIgqB5jC11BYIgCILBiAE8CIKgocQAHgRB0FBiAA+CIGgoMYAHQRA0lBjAgyAIGkoM4IcJIvJ5ETm/S/77ROS/9ljWl0XkVaOrXf+IyNtE5KNDlvEyEfmXEdZptYj8SESOG1WZ88rfJyIPHfDaH4rIk0dcn4+IyF+l40eJyNdGWX5QTwzgDUZE7hCRp/Ryrqqep6rb0nWvFJEb5uW/WlX/ciHquVxR1Y+p6u/n30VEReT0IYq8ELheVXcNWzfvR1JVN6jqbYOUp6q/qapfHrZeXcq/CdgjIs9aqHsEncQAHgSj49XAPyx1JZaQjwH/YakrcTgRA/gKIZeqReQdInK/iNwuIueZ/C+LyKtE5BHA+4AnpCn5npRvp8NHisi1InJfKutaETmxx3p8TkT+xny/UkQ+VHHue0XkHfPSPiMib0zHx4vIJ1M9bheR13W577OTmmBPausjTN5JInJ1KueXInKJ7bN0fH06/XupX14kIj+wEqWITIjIbhE507n/ycBDgRtN2hEicnm6750i8uciMmbu/VURuUREHhCRH4vIuSnvIuB3gEtSXfL6tmcI6Xn9fVKN7UtlHSsi70rP7Me2nna2lvpoX/rbn8o9NeU9U0S+m875mog8ypRxpoh8W0T2isjHgTXzuuHLwLkisrrqOQWjJQbwlcXjgVuAo4H/BVwmImJPUNWbySTFr6cp+WannDHgw8ApwMnAQeCSHuvwR8DLReQcEXkZ8Djg9RXnXgG8KK+jiBwJ/D5wZRroPgt8DzgBOBd4g4g8dX4hIvLrqaw3AMcAnwM+KyKrRKQFXAvcCZyayrpyfhmq+qR0+OjULx8HLgf+0Jz2dGCXqn7Hacu/AW5T1RmT9m7gCLKB/XeBVwD/3uQ/HvgZ2fP6C+BqEdmiqm8FvgK8NtXltc79AF4I/Hm6fhL4OvDt9P0q4J3eRaq6OZW7Abg43WtnGvA/RCZFHwW8H7hGMt3+KuDTZDOMLcA/As+fV+5OYBp4eEV9gxETA/jK4k5V/YCqzgLbgOOArf0Woqq/VNVPquoBVd0LXEQ2APVy7d3Af0z3vxh4RSrD4yuAkkmbAC8g+2H5BfBbwDGq+j9UdSrpfj8AvNgp50XA/1XVL6jqNPAOYC3wb8l+QI4H/ouq7lfVQ6p6g1OGx0eBp4vIpvT95VSrSDYD7XamH44XA29R1b2qegfwN6mMnHuBd6nqdPrBuAV4Ro91A/iUqn5LVQ8BnwIOqerl6fl/HOiYKVhE5EXAS4Hnp367EHi/qt6oqrPJZjIJnJ3+Jkx9rwK+6RS7N/VFsAjEAL6yuDs/UNUD6XBDv4WIyDoReX+a9j8IXA9sToNSL3wWaAG32MEyqTjyqfvvaBZJ7UrgJemUl5LpUSGT/o9PU/k9SdXzZ/g/SMeTSdgAqOoccBeZtH0S2Q/bjHNdV9IPyVeB54vIZuA8U7/53A9sNN+PJhvw7jRpd6Y65ezUcjS5O1NbeuUec3zQ+V757JO0fQnw71T1vpR8CvCmeX1+UqrT8RX1nc9GYE8fbQiGIAbww5O6EJRvIpsGP15VNwG5ekGqLylxEXAzcJyI5INz7gmxIf19JSVfAbxARE4hUyl8MqXfBdyepvv530ZVfbpzv1+QDT5ZJTOVzEnAzlTOySIy3mPd57ONTI3yB2Szg50V590EnGbus5tMnXCKOefkVKecE+apuE5ObYH6ZzQwIvIQMnXIa+apg+4CLprX5+tU9QpgV0V9bbknAKvIZhLBIhAD+OHJPcCJSa/psZFMgtsjIlvI9LM9ISJPItPzvgI4H3h3+sd2SQPIbuCDwD+rai69/SuwV0T+VETWikhLRB4pIr/lFPMJ4Bkicq6ITJD9AE0CX0vl7ALeLiLrRWSNiDyxojr3kOmrLZ8GHkumx7+8Szt2ALeSqWxIaoxPABeJyMb0A/VGMrVMzkOA1yXj6B8AjyDT31fVZWjSD8xVwEdV9RPzsj8AvFpEHi8Z60XkGSKykUy/PmPq+7y8rYbfBb6oqpOjrnfgEwP44ckXgR8Cd4vIbif/XWQ65N3AN4B/6qXQpCu+nMz4tjNJ2ZcBH55vTJ3H/wGekj6B9gD4TOAxwO0Ug/wR8y9W1VvIpOR3p/OeBTwr6c5n0/fTgZ8DO8h05h5vA7Yl9cELU9kHyWYFpwFX13TB+ynruP8TsB+4Dbghtc965NwIPCzV+SLgBar6y5R3MdnM5H4R+bua+/bDiWQ2hzcYddY+ETlZVbcDf0ymWrmf7AfplQCqOgU8L33/FVkfzu+Pl5F5OAWLhMSGDkHQHRH5b8Cvq+of1py3GvgOcG7dYh4ReSXwKlX97ZFVdAlJ7obvV9UnLHVdDicG1QsGwWFBUiFdQFmydkmqgzMWvFLLkLQSMwbvRSZUKEFQgYj8MZlh7/Oqen3d+UGw2IQKJQiCoKGEBB4EQdBQYgAPgiBoKItqxBxfu15XbdqymLdc0dQqv3pddtNkhtQA5l1UWUxdH3bJ1yH7X/pp20JoQgcos2t/DtiX7X708vvpY533ybw+dirtPoOF6Ouadhy8d8duVT1mfvqiDuCrNm3h9Je8cTFvuaJxBwjpPNa6tKrr3Zs6lyzAS97z4FV3Xk1+fh+vjwB0rPzZT37VNbX1TMcyV13fjnxncPKu6fn5VZxbl9e1P51+s8faqshP6XPjNk07zivpE5w+lNl02nRRETEBFsZmyp/2mlJZ3rOqsCVq1+UP+H3jpH3/4jd6YQtChRIEQdBUwg/8cCOXGGySdObXSYx9Te97ZCHKHBW2brWqkS7ScFaYk181rZ+f1kcf1fXnQP09qmdUM+Mo9fFcZ74ruXr9WUqTjjR3djLkrM6SS+a1kviAhAQeBEHQUEICPwzoS3ocRlBYxhL0oLg6XU+Cq+vjfgxwo+7HJTCGlt45L99+mes8zxVYrYQ+VmN+HnOy25WyYrvJdnTcJak/nwl4+u7SffwqLQQhgQdBEDSUGMCDIAgaSqhQguFZzqqTCtVH31RNkbuVWePXPIgRsdIlUJ38ZUKtAdbpD6ulaLsPWpe+3PXQKaeUbkRU6xLYTpvrPB6z55XyOxvg9bevMipSR2nQDAk8CIKgocQAHgRB0FBChbKS6dGnu4R3bo3XRa/T9joPGC+/tuxRqUhKFel+m7I3RMq3afmKQu08r3QbU2htO0fksz3Qqtlh+7VGvdPuBqO6sP2V921pJWaudrGrN+2W213eNW91JhQrMGVW3XN77YfKNRYeg/yPGkICD4IgaCghgQf11BnLhpQiut66okxXmquNTFVNpa98j9K4J00PK7hKzcynbmbklVNiISVvr25W8rVWSued8mKpyLiJYZLKmjFPI4+PAqDjHdVoF1qSuqfNbWays0tGzF77wHctd19J16Bp+6tc3a6EBB4EQdBQYgAPgiBoKKFCOcxwp/d9qD7qpvJemV3jOddRMYUdxpe6jjpDX0ld0s1YZvNseNRBgmHVGY+X0GA5kAHWWa5uDYZjTijdsSmjephK18yazhwz6pR1uV7G3CepRlpTRVrLlJmHlvXC+Jbw3m0/e+iY8HWEBB4EQdBQQgJfadRbUPouxw/oY06tC9jfLZBTjfRWWd3FWnFYJ9k67WjPcmwfjCjYUV9S9wIaKfua+fRqXK5453LJubQSMpVpjZD2RcuDXc2tKq7Jzx2bLC5pmWNvpWbtTHJsXt58FtDAb24fBEEQNI0YwIMgCBpKrQpFRNYA1wOr0/lXqepfiMhpwJXAUcC3gJer6lR1ScGyoNO2UzbQecaykv+uU2YSA+ZanWlQs8Kyyuc6V0PUrFZ08xdidabBU1OUbpkHRapQKRUF9X+fur0o6wyGdYzMKOz4RUupPzofUkktZw2a6aUcK+1fmaWNH7IqEtPJqfypjWKuoeOa1qQxYjrvXHnPzSyjtOKz3IRUkJO/QPQigU8C56jqo4HHAE8TkbOBvwb+VlVPB+4HLli4agZBEATzqR3ANWNf+jqR/hQ4B7gqpW8DnrsgNQyCIAhcevJCEZEWmZrkdOA9wM+APaqaT2p2ACcsSA2DhaFCReL6VzvT6lJgoXQ8500t6+7vLCF2z+uHUapN6jwonHbkU36xXhGmb+by1dPOcnG3bIrnURcP3LJsYoPn7bUBqkor6ZPHiGnEmM1Pvt52+X2uQrG+4VYdAqnD5zo71nqbtKY7O2muZfVytp7zD3DP6ym9C/34jvdkxFTVWVV9DHAi8DjgN3q9gYhcKCLbRWT7zMH9vdcsCIIg6EpffuCqukdEvgQ8AdgsIuNJCj8R2FlxzaXApQDrtp60XGSCw4Nee9tKejUrHOcmss/ZVebyXAyoCjzlbBZbV4+eQ3cu1htVE6zI+igXKwuNxGgCMeWzl75srTUrYGt3vRmEERnjXP9ptxNNH9HZn3aGN57SWpPF9FEOmmtSSFhtFUPc1KYk6VeNes4sx+3OAcIil/MdCX+hwsmKyDEisjkdrwV+D7gZ+BLwgnTa+cBnBqtCEARBMAi9SODHAduSHnwM+ISqXisiPwKuFJG/Ar4DXLaA9QyCIAjmUTuAq+pNwJlO+m1k+vCggdQuubYqA2uASyoU9XZHseU4O5mU4lt7apVhl4aPigHuY6fLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 432x288 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"pylab.title('Initial x-velocity (optimized)')\n",
|
||||
"pylab.imshow(np.concatenate(optimized_velocity_field.data[...,1], axis=1), origin='lower')"
|
||||
|
||||
+9
-35
@@ -1,43 +1,17 @@
|
||||
Discussion
|
||||
=======================
|
||||
|
||||
In the previous sections we've seen example reconstuctions that used physical residuals as soft constraints, in the form of the PINNs, and reconstuctions that used a differentiable physics (DP) solver. While both methods can find minimizers for the same minimization problem, the solutions the obtained differ substantially, as do the behavior of the non-linear optimization problem that we get from each formulation. In the following we discuss these differences in more detail, and we will combine conclusions drawn from the behavior of the Burgers case of the previous sections with observations from research papers.
|
||||
|
||||
## Compatibility with Existing Numerical Methods
|
||||
|
||||
It is very obvious that the PINN implementation is quite simple, which is a positive aspect, but at the same time it differs strongly from "typical" discretziations and solution approaches that are usually to employed equations like Burgers equation. The derivatives are computed via the neural network, and hence rely on a fairly accurate representation of the solution to provide a good direction for optimizaion problems.
|
||||
|
||||
The DP version on the other hand inherently relies on a numerical solver that is tied into the learning process. As such it requires a discretization of the problem at hand, and given this discretization can employ existing, and potentially powerful numerical techniques. This means solutions and derivatives can be evaluated with known and controllable accuracy, and can be evaluated efficiently.
|
||||
|
||||
## Discretization
|
||||
|
||||
The reliance on a suitable discretization requires some understanding and knowledge of the problem under consideration. A sub-optimal discretization can impede the learning prcoess or, worst case, lead to diverging trainig runs. However, given the large body of theory and practical realizations of stable solvers for a wide variety of physical problems, this is typically not an unsurmountable obstacle.
|
||||
|
||||
The PINN approaches on the other hand do not require an a-priori choice of a discretization, and as such seem to be "discretization less". This, however, is only an advantage on first sight. As they yield solutions in a computer, they naturally _have_ to discretize the problem, but they construct this discretization over the coure of the training process, in a way that is not easily controllable from the outside. Thus, the resulting accuracy is determined by how well the training manages to estimate the complexity of the problem for realistic use cases, and how well the training data approximates the unknown regions of the solution.
|
||||
|
||||
As demonstrated with the Burgers example, the PINN solutions typically have significant difficulties propagating information _backward_ in time. This is closely coupled to the efficiency of the method.
|
||||
|
||||
## Efficiency
|
||||
|
||||
The PINN approach typically perform a localized sampling and correction of the solutions, which means the corrections in the form of weight updates are likewise typically local. The fulfilment of boundary conditions in space and time can be correspondingly slow, leading to long training runs in practice.
|
||||
|
||||
A well-chosen discretization of a DP approach can remedy this behavior, and provide an improved flow of gradient information. At the same time, the reliance on a computational grid means that solutions can be obtained very quickly. Given an interpolation scheme or set of basis functions, the solution can be sampled at any point in space or time given a very local neighborhood of the computational grid. Worst case, this can lead to slight memory overheads, e.g., by repeatedly storing mostly constand values of a solution.
|
||||
|
||||
For the PINN representation with fully-connected networks on the other hand, we need to make a full pass over the potentially large number of values in the whole network to obtain a sample of the solution at a single point. The network effectively needs to encode the full high-dimensional solution. Its size likewise determines the efficiency of derivative calculations.
|
||||
**TODO** discussion of DP only
|
||||
|
||||
## Summary
|
||||
|
||||
The following table summarizes these findings:
|
||||
Thus, to summarize ...
|
||||
|
||||
| Method | ✅ Pro | ❌ Con |
|
||||
|----------|-------------|------------|
|
||||
| **PINN** | - Analytic derivatives via back-propagation | - Expensive evaluation of NN, as well as derivative calculations |
|
||||
| | - Simple to implement | - Incompatible with existing numerical methods |
|
||||
| | | - No control of discretization |
|
||||
| **DiffPhys** | - Leverage existing numerical methods | - More complicated to implement |
|
||||
| | - Efficient evaluation of simulation and derivatives | - Require understanding of problem to choose suitable discretization |
|
||||
|
||||
As a summary, both methods are definitely interesting, and leave a lot of room for improvement with more complicated extensions and algorithmic modifications that change and improve on the various negative aspects we have discussed for both sides.
|
||||
|
||||
However, as of this writing, the physics-informed (PI) approach has clear limitations when it comes to performance and campatibility with existing numerical methods. Thus, when knowledge of the problem at hand is available, which typically is the case when we choose a suitable PDE model to constrain the learning process, employing a differentiable physics (DP) solver can significantly improve the training process as well as the quality of the obtained solution. Next, we will target more complex settings, i.e., fluids with Navier Stokes, to illustrate this in more detail.
|
||||
✅ Pro:
|
||||
- uses physical model and numerical methods for discretization
|
||||
- efficiency of selected methods carries over to training
|
||||
- tight coupling of physical models and NNs possible
|
||||
|
||||
❌ Con:
|
||||
- not compatible with all simulators (need to provide gradients)
|
||||
- ...
|
||||
@@ -0,0 +1,44 @@
|
||||
Diff. Physics versus Phys.-informed Training
|
||||
=======================
|
||||
|
||||
In the previous sections we've seen example reconstuctions that used physical residuals as soft constraints, in the form of the PINNs, and reconstuctions that used a differentiable physics (DP) solver. While both methods can find minimizers for the same minimization problem, the solutions the obtained differ substantially, as do the behavior of the non-linear optimization problem that we get from each formulation. In the following we discuss these differences in more detail, and we will combine conclusions drawn from the behavior of the Burgers case of the previous sections with observations from research papers.
|
||||
|
||||
## Compatibility with Existing Numerical Methods
|
||||
|
||||
It is very obvious that the PINN implementation is quite simple, which is a positive aspect, but at the same time it differs strongly from "typical" discretziations and solution approaches that are usually to employed equations like Burgers equation. The derivatives are computed via the neural network, and hence rely on a fairly accurate representation of the solution to provide a good direction for optimizaion problems.
|
||||
|
||||
The DP version on the other hand inherently relies on a numerical solver that is tied into the learning process. As such it requires a discretization of the problem at hand, and given this discretization can employ existing, and potentially powerful numerical techniques. This means solutions and derivatives can be evaluated with known and controllable accuracy, and can be evaluated efficiently.
|
||||
|
||||
## Discretization
|
||||
|
||||
The reliance on a suitable discretization requires some understanding and knowledge of the problem under consideration. A sub-optimal discretization can impede the learning prcoess or, worst case, lead to diverging trainig runs. However, given the large body of theory and practical realizations of stable solvers for a wide variety of physical problems, this is typically not an unsurmountable obstacle.
|
||||
|
||||
The PINN approaches on the other hand do not require an a-priori choice of a discretization, and as such seem to be "discretization less". This, however, is only an advantage on first sight. As they yield solutions in a computer, they naturally _have_ to discretize the problem, but they construct this discretization over the coure of the training process, in a way that is not easily controllable from the outside. Thus, the resulting accuracy is determined by how well the training manages to estimate the complexity of the problem for realistic use cases, and how well the training data approximates the unknown regions of the solution.
|
||||
|
||||
As demonstrated with the Burgers example, the PINN solutions typically have significant difficulties propagating information _backward_ in time. This is closely coupled to the efficiency of the method.
|
||||
|
||||
## Efficiency
|
||||
|
||||
The PINN approach typically perform a localized sampling and correction of the solutions, which means the corrections in the form of weight updates are likewise typically local. The fulfilment of boundary conditions in space and time can be correspondingly slow, leading to long training runs in practice.
|
||||
|
||||
A well-chosen discretization of a DP approach can remedy this behavior, and provide an improved flow of gradient information. At the same time, the reliance on a computational grid means that solutions can be obtained very quickly. Given an interpolation scheme or set of basis functions, the solution can be sampled at any point in space or time given a very local neighborhood of the computational grid. Worst case, this can lead to slight memory overheads, e.g., by repeatedly storing mostly constand values of a solution.
|
||||
|
||||
For the PINN representation with fully-connected networks on the other hand, we need to make a full pass over the potentially large number of values in the whole network to obtain a sample of the solution at a single point. The network effectively needs to encode the full high-dimensional solution. Its size likewise determines the efficiency of derivative calculations.
|
||||
|
||||
## Summary
|
||||
|
||||
The following table summarizes these findings:
|
||||
|
||||
| Method | ✅ Pro | ❌ Con |
|
||||
|----------|-------------|------------|
|
||||
| **PINN** | - Analytic derivatives via back-propagation | - Expensive evaluation of NN, as well as derivative calculations |
|
||||
| | - Simple to implement | - Incompatible with existing numerical methods |
|
||||
| | | - No control of discretization |
|
||||
| **DiffPhys** | - Leverage existing numerical methods | - More complicated to implement |
|
||||
| | - Efficient evaluation of simulation and derivatives | - Require understanding of problem to choose suitable discretization |
|
||||
|
||||
As a summary, both methods are definitely interesting, and have a lot of potential. There are numerous more complicated extensions and algorithmic modifications that change and improve on the various negative aspects we have discussed for both sides.
|
||||
|
||||
However, as of this writing, the physics-informed (PI) approach has clear limitations when it comes to performance and campatibility with existing numerical methods. Thus, when knowledge of the problem at hand is available, which typically is the case when we choose a suitable PDE model to constrain the learning process, employing a differentiable physics (DP) solver can significantly improve the training process as well as the quality of the obtained solution. Next, we will target more complex settings, i.e., fluids with Navier Stokes, to illustrate this in more detail.
|
||||
|
||||
|
||||
+1
-1
@@ -150,7 +150,7 @@ fundamental steps. Here are some considerations for skipping ahead to the later
|
||||
|
||||
- On the other hand, if you're already deep into ANNs&Co, and you'd like to skip ahead to the research related topics, we recommend starting in the _Physical Loss Terms_ chapter, which lays the foundations for the next chapters.
|
||||
|
||||
A brief look at our _Notation_ won't hurt in both cases, though!
|
||||
A brief look at our _notation_ in the {doc}`notation` chapter won't hurt in both cases, though!
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
@@ -0,0 +1,600 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Burgers Optimization with a Physics-Informed NN\n",
|
||||
"\n",
|
||||
"To illustrate how the physics-informed losses work, let's consider a \n",
|
||||
"reconstruction example with a simple yet non-linear equation in 1D:\n",
|
||||
"Burgers equation $\\frac{\\partial u}{\\partial{t}} + u \\nabla u = \\nu \\nabla \\cdot \\nabla u$\n",
|
||||
"for which we have a series of _observations_ at time $t=0.5$ which should be fulfilled in the solution. In addition, let's impose Dirichlet boundary conditions $u=0$\n",
|
||||
"at the sides of our computational domain, and define the solution in\n",
|
||||
"the time interval $t \\in [0,1]$.\n",
|
||||
"\n",
|
||||
"Note that similar to the previous forward simulation example, \n",
|
||||
"we will still be sampling the solution with 128 points ($n=128), but now we have a discretization via the NN. So we could also sample points inbetween without having to explicitly choose a basis function for interpolation. The discretization via the NN now internally determines how to use its degrees of freedom to construct the basis functions. So we have no direct control over the reconstruction.\n",
|
||||
"\n",
|
||||
"## Preliminaries\n",
|
||||
"\n",
|
||||
"Let's just load TF and phiflow for now, and initialize the random sampling."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"#from phi.tf.flow import *\n",
|
||||
"import numpy as np\n",
|
||||
"import numpy.random as rnd\n",
|
||||
"import tensorflow as tf\n",
|
||||
"\n",
|
||||
"#rnd = TF_BACKEND # for phiflow: sample different points in the domain each iteration\n",
|
||||
"#rnd = math.choose_backend(1) # use same random points for all iterations"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We're importing phiflow here, but we won't use it to compute a solution to the PDE. Below we'll instead use the derivatives of an NN (as explained in the previous section) to set up a loss formulation to train an NN.\n",
|
||||
"\n",
|
||||
"Next, we set up a simple NN with 8 fully connected layers and `tanh` activations with 20 units each. \n",
|
||||
"\n",
|
||||
"We'll also define the `boundary_tx` function which gives an array of constraints for the solution (all for $=0.5$ in this example), and the `open_boundary` function which contains constraints for $x= \\pm1$ being 0."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def network(x, t):\n",
|
||||
" \"\"\" Dense neural network with 8 hidden layers and 3021 parameters in total.\n",
|
||||
" Parameters will only be allocated once (auto reuse).\n",
|
||||
" \"\"\"\n",
|
||||
" y = tf.stack([x, t], axis=-1)\n",
|
||||
" for i in range(8):\n",
|
||||
" #y = tf.keras.layers.dense(y, 20, activation=tf.math.tanh, name='layer%d' % i, reuse=tf.AUTO_REUSE)\n",
|
||||
" y = tf.compat.v1.layers.dense(y, 20, activation=tf.math.tanh, name='layer%d' % i, reuse=True)\n",
|
||||
" return tf.compat.v1.layers.dense(y, 1, activation=None, name='layer_out', reuse=True)\n",
|
||||
"\n",
|
||||
"def boundary_tx(N):\n",
|
||||
" x = np.linspace(-1,1,128)\n",
|
||||
" # precomputed solution from forward simulation:\n",
|
||||
" u = np.asarray( [0.008612174447657694, 0.02584669669548606, 0.043136357266407785, 0.060491074685516746, 0.07793926183951633, 0.0954779141740818, 0.11311894389663882, 0.1308497114054023, 0.14867023658641343, 0.1665634396808965, 0.18452263429574314, 0.20253084411376132, 0.22057828799835133, 0.23865132431365316, 0.25673879161339097, 0.27483167307082423, 0.2929182325574904, 0.3109944766354339, 0.3290477753208284, 0.34707880794585116, 0.36507311960102307, 0.38303584302507954, 0.40094962955534186, 0.4188235294008765, 0.4366357052408043, 0.45439856841363885, 0.4720845505219581, 0.4897081943759776, 0.5072391070000235, 0.5247011051514834, 0.542067187709797, 0.5593576751669057, 0.5765465453632126, 0.5936507311857876, 0.6106452944663003, 0.6275435911624945, 0.6443221318186165, 0.6609900633731869, 0.67752574922899, 0.6939334022562877, 0.7101938106059631, 0.7263049537163667, 0.7422506131457406, 0.7580207366534812, 0.7736033721649875, 0.7889776974379873, 0.8041371279965555, 0.8190465276590387, 0.8337064887158392, 0.8480617965162781, 0.8621229412131242, 0.8758057344502199, 0.8891341984763013, 0.9019806505391214, 0.9143881632159129, 0.9261597966464793, 0.9373647624856912, 0.9476871303793314, 0.9572273019669029, 0.9654367940878237, 0.9724097482283165, 0.9767381835635638, 0.9669484658390122, 0.659083299684951, -0.659083180712816, -0.9669485121167052, -0.9767382069792288, -0.9724097635533602, -0.9654367970450167, -0.9572273263645859, -0.9476871280825523, -0.9373647681120841, -0.9261598056102645, -0.9143881718456056, -0.9019807055316369, -0.8891341634240081, -0.8758057205293912, -0.8621229450911845, -0.8480618138204272, -0.833706571569058, -0.8190466131476127, -0.8041372124868691, -0.7889777195422356, -0.7736033858767385, -0.758020740007683, -0.7422507481169578, -0.7263049162371344, -0.7101938950789042, -0.6939334061553678, -0.677525822052029, -0.6609901538934517, -0.6443222327338847, -0.6275436932970322, -0.6106454472814152, -0.5936507836778451, -0.5765466491708988, -0.5593578078967361, -0.5420672759411125, -0.5247011730988912, -0.5072391580614087, -0.4897082914472909, -0.47208460952428394, -0.4543985995006753, -0.4366355580500639, -0.41882350871539187, -0.40094955631843376, -0.38303594105786365, -0.36507302109186685, -0.3470786936847069, -0.3290476440540586, -0.31099441589505206, -0.2929180880304103, -0.27483158663081614, -0.2567388003912687, -0.2386513127155433, -0.22057831776499126, -0.20253089403524566, -0.18452269630486776, -0.1665634500729787, -0.14867027528284874, -0.13084990929476334, -0.1131191325854089, -0.09547794429803691, -0.07793928430794522, -0.06049114408297565, -0.0431364527809777, -0.025846763281087953, -0.00861212501518312] );\n",
|
||||
" t = np.asarray( np.ones(x.shape)) * 0.5\n",
|
||||
" perm = np.random.permutation(128) \n",
|
||||
" return (x[perm])[0:N], (t[perm])[0:N], (u[perm])[0:N]\n",
|
||||
"\n",
|
||||
"def _ALT_t0(N): # alternative, impose original initial state at t=0\n",
|
||||
" x = rnd.uniform(-1, 1, [N])\n",
|
||||
" t = np.zeros(x.shape)\n",
|
||||
" u = -np.sin(np.pi * x)\n",
|
||||
" return x, t, u\n",
|
||||
"\n",
|
||||
"def open_boundary(N):\n",
|
||||
" t = rnd.uniform(0, 1, [N])\n",
|
||||
" x = np.concatenate([np.zeros([N//2]) + 1, np.zeros([N//2]) - 1], axis=0)\n",
|
||||
" u = np.zeros([N])\n",
|
||||
" return x, t, u\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Most importantly, we can now also construct the residual loss function `f` that we'd like to minimize in order to guide the NN to retrieve a solution for our model equation. As can be seen in the equation at the top, we need derivatives w.r.t. $t$, $x$ and a second derivative for $x$. The first three lines of `f` below do just that.\n",
|
||||
"\n",
|
||||
"Afterwards, we simply combine the derivates according to Burgers equation:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def f(u, x, t):\n",
|
||||
" \"\"\" Physics-based loss function with Burgers equation \"\"\"\n",
|
||||
" #with tf.GradientTape() as tape: \n",
|
||||
" u_t = tf.compat.v1.gradients(u, t)\n",
|
||||
" u_x = tf.compat.v1.gradients(u, x)\n",
|
||||
" u_xx = tf.compat.v1.gradients(u_x, x)\n",
|
||||
" return u_t + u*u_x - (0.01 / np.pi) * u_xx\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Next, let's set up the sampling points in the inner domain, such that we can compare the solution with the previous forward simulation in phiflow. The last `math.expand_dims()` call simply adds another `batch` dimension, so that the resulting tensor is compatible with the following examples.\n",
|
||||
"\n",
|
||||
"A bit more interesting: grid_u will afterwards contain a full graph to evaluate our NN at $128 \\times 33$ positions, and will retrieve the results in a $[1,128,33,1]$ array once we run it through `session.run`."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"WARNING:tensorflow:From <ipython-input-2-5de5bf78dcb3>:8: dense (from tensorflow.python.keras.legacy_tf_layers.core) is deprecated and will be removed in a future version.\n",
|
||||
"Instructions for updating:\n",
|
||||
"Use keras.layers.Dense instead.\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/keras/legacy_tf_layers/core.py:187: Layer.apply (from tensorflow.python.keras.engine.base_layer) is deprecated and will be removed in a future version.\n",
|
||||
"Instructions for updating:\n",
|
||||
"Please use `layer.__call__` method instead.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# generate array with positions: \n",
|
||||
"# -1 to 1 spatial with 128 cells\n",
|
||||
"# 0 to 1 in time with 1+32 timesteps\n",
|
||||
"N=128\n",
|
||||
"grids_xt = np.meshgrid(np.linspace(-1, 1, N), np.linspace(0, 1, 33), indexing='ij')\n",
|
||||
"grid_x, grid_t = [tf.convert_to_tensor(t, tf.float32) for t in grids_xt]\n",
|
||||
"\n",
|
||||
"# create 4D tensor with batch and channel dimensions in addition to space and time\n",
|
||||
"# in this case gives shape=(1, n, 33, 1)\n",
|
||||
"grid_u = tf.reshape(network(grid_x, grid_t), (1,N,33,1) )"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Let's give this a try: we can initialize a TF session, evaluate `grid_u` and show it in an image, just like the phiflow solution we computed previously. \n",
|
||||
"\n",
|
||||
"(Note, as before the x axis does not show actual simulation time, but is showing 32 steps \"blown\" up by a factor of 16 to make the changes over time easier to see in the image.)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {
|
||||
"scrolled": true
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Size of grid_u: (1, 128, 33, 1)\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/util/tf_should_use.py:247: initialize_all_variables (from tensorflow.python.ops.variables) is deprecated and will be removed after 2017-03-02.\n",
|
||||
"Instructions for updating:\n",
|
||||
"Use `tf.global_variables_initializer` instead.\n",
|
||||
"Randomly initialized network state:\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAB0CAYAAAB+I3LkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbfklEQVR4nO2dbahlV3nH/89MfGsUbaJOpyZtlAZFStUqvqAfokWJIuaLBG2haQnMlwoWCnWk0NJ+Sr/UWlrEgQYVWl9KGwwi1TQq/dJqfDe+REeJmCE6KNEKOvecvdfTD3vfeE3mzvk9c5991777rD8Mc++5a9be68w9v/U///Wsdczd1dTU1NS0LB2rfQNNTU1NTflqcG9qampaoBrcm5qamhaoBvempqamBarBvampqWmBanBvampqWqAmgbuZ3Whm95nZWTM7PcU1mpqampr2l2XXuZvZcUnflPRqSQ9IukfSm939a6kXampqamraV1M49xdLOuvu33H3laQPSLppgus0NTU1Ne2jKeD+DEnf2/P9A+NjTU1NTU2HpCtqXdjMTkk6JUlXXvn4Fz77Ob+x+R/BCMlo0uQuCTYO9ZnYTg6vHbgu7a/WmKdoN/cxl+zr0u4C/aU/N7nd8ftLvq4kubHLJuPmyw/96Ifu/rSL/WwKuJ+TdO2e768ZH/slufsZSWck6YUvfLZ/+n/etbnnvkM3YKVIBbxaSpGctHOpY9dG7bzIup71t17B6/ZwzM6ex+IyMhbanyStYLv1Orc/2q7vpY793jhp5z78v2zsz6UV+33wNWy3grNFV9DE4l1BxPHO5R1Dk69YuwJ/HWg774y97HtjgC+m0gO4F9ZfcdafJP36Hbd/d7+fTQH3eyRdb2bP1AD1N0n6/Uv/E2eQ1QjuTaJw9zK8sLL62+1zY3/O+yP3J433CNr2XWDM2e2yxwyfR/r/3LF79K4MEwG5PzIJSPLQmGm7zX16B9v1ketubuZFcuhvHM7NKgFow7nZu83JtfsAbtKu71l/BbTbpHS4u3tnZm+R9DFJxyXd7u5f3fgPwYtqEkeePFlgR45BRx0YvS4bs9F3Al74tQkQJd4fbsfg7l0PJ0jYLgB39XASgK5YnbMJo4ftOmf+q7B7dOrnyhiPAHiWTnIQj5Se9Yeddm8qJJYppgKv25cZwl2S3P2jkj5K25tXdOQoekieLCQOOuqM4ZitCzh35E7phAsjCmmSyYKBLjIW1s4htIWhTScBHqOwiY/C2JEj98Ic+eCec2MU746xdhDufQ/7c4PO3dSDSWCTqi2o/rIgFGm2iyHSpU4W+J2FByKKbCDSCajvYV4ciCio28XQDoCOqBQMbdLOI7EMbgfzbBijDHBn7Vj+DKFdpLJmAMMxSm8oHqF5dl9gf7SdM+c+9DlT5x6WizsmogniG7a4CEEnpbtTy87cOwh3GstQ5x4CYuA5xJk7cbG7WcEGFb64iJ17cozinWDmzl4q3g2QJaKQjU0Cm9tSaNPsu7ihOMhdKG5ZnnNHkA04sMyMPNuRhxYXAw6fPIc0V470h3Pl7IkvUhHC4I7jG5IrZzt351UwBRZZcYc/ZtUbL2wq8A026k98EujA4udufwjGMPt2F8vcIbSHdkty7uQFTcsRoTsNLdDS/DkzypD4Am0PF0BDi4akPzhm55EHjlGyF1R7OgnAe6TOvbh8DeMW+BLwnpYjQkc+Zt8b2xUGY1qJIgUcPm53DMG4FLhQGnDa2+ncs2MUBCa4SkTr3AM13yjmkWLrDDRLz4xlaEZeXB6pNyfKjlG6gmFMnbuvk7N03E7IkdPKEVqJMpQ4wsoROAmQRUiJRSiS1MF2fTnGqm8CcQubLKTeF+XcwQsfu124UErdKXXkkbLAKdplxzKZtdylpMcoDjf+aN3D6IFOFi4BZ+x4ElC1DT1lRSOKetDuocMfYEzbMWiTdn1h7wQi0F6WcwdgCrldXAqZWRESqHOP5M9zj1GII59gt6ZgjELL/XDc4tBBFw9k37Bdx96E1CwfRLnyBDXf3JGzCS0C7aLk+GZJzh1BMTuioG43eyt+cWkFrdU6sMU+NZYprBolEMtU22K/ogulgUoUWO6Hsu/CF0Bp5UhPK0zW2wlt2l8NaJdlxTIw045EFImxjEVimfRcObu6hY2F59T0PBbH0MZlgRDuvoZwD5QPUmhn13xHoI3aQefe98fSoc0rR2C1TH8cX5dAm8Y8xU0O+otAe1mxDIEdjjxgtQzdrdl3ubXcEoc2PmyLjSW0xR7HMombjgrP0rN3YZaVeM03du65lR6RjTqZuytphUms3K8etGl/taDdged6k+YBdxpTTAHt1Gze8080hLGMr+iYAxt66C7MzBilOI5bfB2ANuwvc6MOLQuUmHPfXaysAe3sssBILfcawzgX2l0A7miBVnwS6Jbl3DdDzGhOHdqAk5kr00XIwrN0ei4KzuYZ3IeTD5NhTOFOa75XEMZwo05Z8cXK7N2VPYxROriwuY3Q7nC7+UO7X4xzd2fgnvsuzMg5K9C5+xquH0SgnVzzPQW0icoObBeANq4wqVAWuLsdPhPapRxj59JBaBc3DCa6+5OCk08WPBqhkwCaIMWhvZxYhu5QXa/zoV2p5hsvLl5g71Z81WNo45pvUqPdFb64iGu5udNm/eWXBS4F2rzmm0G7ZuXIkLlvFl0Ajbhs8ptdnN3fbp8H1UzgDrPqFdxRGoF26uJiANrUuU+xUScb2sjhj9Am7RYG7Y3354bL/bJ3V2aX+4WcO3bkuf0NscxmlWRou1vgfLglwZ3ALuLcs6FNP1aNQvsCXVCF0KYVJiteYVIT2kTdKlDLnXhaYCnsXG6Jw5i3Y5UjkTwbLdA6eycQgztrh8cSuC569yOeuaNoS4afQzgHXFLzgTt17qkVIX06tMlGnVCFyYU60K5Vy+1F6sEBVVI9aGeX+02xCIkdPoQ2bTe40/nDnahzpWfu9NTy5Tj34ixbXlJZYDa0V6XKBpzSC58WSKHd7zB3So947bv8PJvmymvotGmlxxpDm42lh2eRD+6UAMxSK0JcfHExHe44vmGTgBQ5tn8pcHdH4J49tCMbcEiFicN2ksoFulgptliZDe1ioXI/olrQLtBp00qPEliEjLhYAu25V4RMkVNntyts/R+fQTNMFgfXPOBeXJ7o3HFZYE/r0iG0i/gRrwDaodMCd1Cz0CIk/XAGAu3hhAkKd+h2e75oSM8cQXBXLtwlDu3I4mJmvjtFTh1x0ETwCPsA3Fl/BWbpNJYpYgUAmzQPuEPnjqG9pjXftHRRqZ+UUxPa/Tq/ciQT2h6sHEHRQ2jREE4WidvSh3b5zp3GMrg/0I5GFJExRzLyzP6ogy7QkTt0+BKfWC6l+cAdVI9gaK/yz+VGMKaVI5IcVI5EDpSaonKEwd0gtAOVHhG4Jx4AlV0RQhcXIzk1zZ8jOXV2pUctB83HzPqjDnp4t8LabeGCquQ7JHOH0D4C5X79DvxF7Ezkk3J6eHRr5MwREj30xdAntZcI3JMrQtbwXO7I4mJNR05EKz0ikQfO5lGP84c7ddqRjDzyruGgOhDczex+ST+V1Evq3P1FZnaVpA9Kuk7S/ZJudveHLtWPu6PqEd/heXZm5Yh3sNzPpQIddL+CAKPtOpg/w6NbI5tg6Lb07IoQPAmUY+hFmu7clV/pwTP3XOfO3Wm9/JlPAqzdcI+5JY6Zz80mZTj3V7r7D/d8f1rS3e5+m5mdHr9/2yV7KK7ycwB3mFP7ylmuTD8cuGef6O5u+HQ/Wha4XtOywOPw+PXjeFt6dkUIhntgcZE6cuzcwXUjG3UISGLlg9CRT5A/14tbWLsp3HM2jPk9woaX0BSxzE2Sbhi/fq+kT2kj3CUHpXyRM0cyT/cbsu/kipAVrwghYOKVIzyioHDnMM5vl13pQV742YuGscVF1CwUUWQuBk6xaJgdZeDMPbCgSpTt8DfpoHB3SR83M5f0bnc/I+mEuz84/vz7kk5s7MSFTvjraeUI/NDfocIEtOv5R4yRD/OlZ4lEKkfWPXPu68La9X4MuUkXPHNEE8A9OXqgk4Arv4KDxy2sv/xyv4iLzd6tmdsuG8ZTwD2B7QeG+yvc/ZyZPV3SXWb2jb0/dHcfwf8omdkpSack6donXqn+52TRENZow8XF0tPt5vllfKtALTdxVrgsMBBRMNfJ21E3mQ13utNwivx5CqedWZlRE8bZkUctuO+2Jdel6wf03c+ldCC4u/u58e/zZnaHpBdL+oGZnXT3B83spKTz+/zbM5LOSNILrr7a1xeI084t46MVIbSML1IRsgqUBZIIYIqcOhXuyl80jCwusgXV/MXAbNeZ7dzJ9g2Jj7kmjHEVDGuWDuNIVU1J2MV02XA3syslHXP3n45fv0bS30i6U9Itkm4b//7wpr7cDZ0nssbb14/z405huV9mRUhkcXEnUDnCPzpss3paow23h9NFw93Ig8Uom9tIu1vsYbtU514PxhHQZUO72iIkfG6oK54Cxpm/X5t0EOd+QtIdZrbbz7+6+3+a2T2SPmRmt0r6rqSbN3XkblqBBcZaFSHZi4shuCcvLnb4zBFeJ01BQmOUSP6MNjEFnF+NKCNSiZIN4+yIIgRtPDnD/lizEIjpRIAd/lHI3N39O5Ked5HHfyTp9yJ9lSLtrDbfCq0IoYuLvbOzRCKLixRMO/S0wEBGTqOHiIPeJAqmSEQRgztrx/qjuzADGWtiO6micw/EE7k5NStr3r02um4AxNkTxpGAe6bcTTtrBvfsxUW605A48gjckcNXJH+OlAWSdswVF48sbKJm2KnhrfgBIGZDu1ZFSGR3JepP9aKHbCBGQBzJ54kc9jiHUsgUuUw73eZbiSwu0tK3XLjnV3rsgIXh3f5QZIU36kDnLujcYX+SBA/WDAExE8ZTRBT5tdd1ooej4HYjIM6OmPhzfXC6zwLuxQ1Vj0QWFwlIpqgIIU47VhbIc2Vcyw2uG8mLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 432x288 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# test image output\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"print(\"Size of grid_u: \"+format(grid_u.shape))\n",
|
||||
"\n",
|
||||
"def show_state(a):\n",
|
||||
" # we only have 33 time steps, blow up by a factor of 2^4 to make it easier to see\n",
|
||||
" # (could also be done with more evaluations of network)\n",
|
||||
" for i in range(4):\n",
|
||||
" a = np.concatenate( [a,a] , axis=3)\n",
|
||||
" a = np.reshape( a, [a.shape[1],a.shape[2]*a.shape[3]] )\n",
|
||||
" #print(a.shape)\n",
|
||||
" plt.imshow(a, origin='upper', cmap='magma')\n",
|
||||
" \n",
|
||||
"print(\"Randomly initialized network state:\")\n",
|
||||
"#show_state(session.run(grid_u))\n",
|
||||
"show_state(grid_u)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This visualization already shows a smooth transition over space and time. So far, this is purely the random initialization of the NN that we're sampling here. So it has nothing to do with a solution of our PDE-based model up to now.\n",
|
||||
"\n",
|
||||
"The next steps will actually evaluate the constraints in terms of data (from the `boundary` functions), and the model constraints from `f` to retrieve an actual solution to the PDE."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## The loss functions & training\n",
|
||||
"\n",
|
||||
"As objective for the learning process we can now combine the _direct_ constraints, i.e., the solution at $t=0.5$ and the Dirchlet $u=0$ boundary conditions with the loss from the PDE residuals. For both boundary constraints we'll use 100 points below, and then sample the solution in the inner region with an additional 1000 points.\n",
|
||||
"\n",
|
||||
"The direct constraints are evaluated via `network(x, t)[:, 0] - u`, where `x` and `t` are the \"position\" where we'd like to sample out solution, and `u` provides the corresponding ground truth value.\n",
|
||||
"\n",
|
||||
"For the physical loss points, we have no ground truth solutions, but we'll only evaluate the PDE residual via the NN derivatives, to see whether the solution satisfies PDE model. If not, this directly gives us an error to be reduced via a update step in the optimization. The corresponding expression is of the form `f(network(x, t)[:, 0], x, t)` below. Note that for both data and physics terms the `network()[:, 0]` simply discard the last size-1 dimension of the $(n,1)$ tensor returned by the network."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"ERROR:tensorflow:==================================\n",
|
||||
"Object was never used (type <class 'tensorflow.python.framework.ops.Operation'>):\n",
|
||||
"<tf.Operation 'init' type=NoOp>\n",
|
||||
"If you want to mark it as used call its \"mark_used()\" method.\n",
|
||||
"It was originally created here:\n",
|
||||
" File \"/Users/thuerey/miniconda3/envs/tf/lib/python3.8/site-packages/IPython/core/interactiveshell.py\", line 3155, in run_cell_async\n",
|
||||
" has_raised = await self.run_ast_nodes(code_ast.body, cell_name, File \"/Users/thuerey/miniconda3/envs/tf/lib/python3.8/site-packages/IPython/core/interactiveshell.py\", line 3347, in run_ast_nodes\n",
|
||||
" if (await self.run_code(code, result, async_=asy)): File \"/Users/thuerey/miniconda3/envs/tf/lib/python3.8/site-packages/IPython/core/interactiveshell.py\", line 3427, in run_code\n",
|
||||
" exec(code_obj, self.user_global_ns, self.user_ns) File \"<ipython-input-7-351d69eca4ac>\", line 9, in <module>\n",
|
||||
" tf.compat.v1.initialize_all_variables() File \"/Users/thuerey/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/util/tf_should_use.py\", line 247, in wrapped\n",
|
||||
" return _add_should_use_warning(fn(*args, **kwargs),\n",
|
||||
"==================================\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"ename": "ValueError",
|
||||
"evalue": "Variable layer0/kernel does not exist, or was not created with tf.get_variable(). Did you mean to set reuse=tf.AUTO_REUSE in VarScope?",
|
||||
"output_type": "error",
|
||||
"traceback": [
|
||||
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
|
||||
"\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)",
|
||||
"\u001b[0;32m<ipython-input-7-351d69eca4ac>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 11\u001b[0m \u001b[0;31m#with app.model_scope():\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m \u001b[0mloss_u\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0.5\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mnetwork\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx_bc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt_bc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mu_bc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m \u001b[0;31m# normalizes by first dimension, N_bc\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 13\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 14\u001b[0m \u001b[0;31m# Physics loss inside of domain\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m<ipython-input-2-5de5bf78dcb3>\u001b[0m in \u001b[0;36mnetwork\u001b[0;34m(x, t)\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m8\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0;31m#y = tf.keras.layers.dense(y, 20, activation=tf.math.tanh, name='layer%d' % i, reuse=tf.AUTO_REUSE)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 8\u001b[0;31m \u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompat\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mv1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlayers\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdense\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m20\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mactivation\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmath\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtanh\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'layer%d'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mreuse\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 9\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mtf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompat\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mv1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlayers\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdense\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mactivation\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'layer_out'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mreuse\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/util/deprecation.py\u001b[0m in \u001b[0;36mnew_func\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m 322\u001b[0m \u001b[0;34m'in a future version'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mdate\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m'after %s'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mdate\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 323\u001b[0m instructions)\n\u001b[0;32m--> 324\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 325\u001b[0m return tf_decorator.make_decorator(\n\u001b[1;32m 326\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnew_func\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'deprecated'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/keras/legacy_tf_layers/core.py\u001b[0m in \u001b[0;36mdense\u001b[0;34m(inputs, units, activation, use_bias, kernel_initializer, bias_initializer, kernel_regularizer, bias_regularizer, activity_regularizer, kernel_constraint, bias_constraint, trainable, name, reuse)\u001b[0m\n\u001b[1;32m 185\u001b[0m \u001b[0m_scope\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 186\u001b[0m _reuse=reuse)\n\u001b[0;32m--> 187\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mlayer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mapply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 188\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 189\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/util/deprecation.py\u001b[0m in \u001b[0;36mnew_func\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m 322\u001b[0m \u001b[0;34m'in a future version'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mdate\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m'after %s'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mdate\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 323\u001b[0m instructions)\n\u001b[0;32m--> 324\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 325\u001b[0m return tf_decorator.make_decorator(\n\u001b[1;32m 326\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnew_func\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'deprecated'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/keras/engine/base_layer_v1.py\u001b[0m in \u001b[0;36mapply\u001b[0;34m(self, inputs, *args, **kwargs)\u001b[0m\n\u001b[1;32m 1699\u001b[0m \u001b[0mOutput\u001b[0m \u001b[0mtensor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1700\u001b[0m \"\"\"\n\u001b[0;32m-> 1701\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1702\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1703\u001b[0m @deprecation.deprecated(\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/keras/legacy_tf_layers/base.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, inputs, *args, **kwargs)\u001b[0m\n\u001b[1;32m 545\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 546\u001b[0m \u001b[0;31m# Actually call layer\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 547\u001b[0;31m \u001b[0moutputs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msuper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mLayer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 548\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 549\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mcontext\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexecuting_eagerly\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/keras/engine/base_layer_v1.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 755\u001b[0m \u001b[0;31m# Build layer if applicable (if the `build` method has been\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 756\u001b[0m \u001b[0;31m# overridden).\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 757\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_maybe_build\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 758\u001b[0m \u001b[0mcast_inputs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_maybe_cast_inputs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 759\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/keras/engine/base_layer_v1.py\u001b[0m in \u001b[0;36m_maybe_build\u001b[0;34m(self, inputs)\u001b[0m\n\u001b[1;32m 2096\u001b[0m \u001b[0;31m# operations.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2097\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mtf_utils\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmaybe_init_scope\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2098\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuild\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minput_shapes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2099\u001b[0m \u001b[0;31m# We must set self.built since user defined build functions are not\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2100\u001b[0m \u001b[0;31m# constrained to set self.built.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/keras/layers/core.py\u001b[0m in \u001b[0;36mbuild\u001b[0;34m(self, input_shape)\u001b[0m\n\u001b[1;32m 1169\u001b[0m 'should be defined. Found `None`.')\n\u001b[1;32m 1170\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minput_spec\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mInputSpec\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmin_ndim\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxes\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m{\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mlast_dim\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1171\u001b[0;31m self.kernel = self.add_weight(\n\u001b[0m\u001b[1;32m 1172\u001b[0m \u001b[0;34m'kernel'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1173\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mlast_dim\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munits\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/keras/legacy_tf_layers/base.py\u001b[0m in \u001b[0;36madd_weight\u001b[0;34m(self, name, shape, dtype, initializer, regularizer, trainable, constraint, use_resource, synchronization, aggregation, partitioner, **kwargs)\u001b[0m\n\u001b[1;32m 446\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0minitializer\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 447\u001b[0m \u001b[0minitializer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscope\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minitializer\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 448\u001b[0;31m variable = super(Layer, self).add_weight(\n\u001b[0m\u001b[1;32m 449\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 450\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/keras/engine/base_layer_v1.py\u001b[0m in \u001b[0;36madd_weight\u001b[0;34m(self, name, shape, dtype, initializer, regularizer, trainable, constraint, partitioner, use_resource, synchronization, aggregation, **kwargs)\u001b[0m\n\u001b[1;32m 429\u001b[0m \u001b[0mcaching_device\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 430\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 431\u001b[0;31m variable = self._add_variable_with_custom_getter(\n\u001b[0m\u001b[1;32m 432\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 433\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/training/tracking/base.py\u001b[0m in \u001b[0;36m_add_variable_with_custom_getter\u001b[0;34m(self, name, shape, dtype, initializer, getter, overwrite, **kwargs_for_getter)\u001b[0m\n\u001b[1;32m 743\u001b[0m \u001b[0minitializer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcheckpoint_initializer\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 744\u001b[0m \u001b[0mshape\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 745\u001b[0;31m new_variable = getter(\n\u001b[0m\u001b[1;32m 746\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 747\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/ops/variable_scope.py\u001b[0m in \u001b[0;36mget_variable\u001b[0;34m(name, shape, dtype, initializer, regularizer, trainable, collections, caching_device, partitioner, validate_shape, use_resource, custom_getter, constraint, synchronization, aggregation)\u001b[0m\n\u001b[1;32m 1554\u001b[0m \u001b[0msynchronization\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mVariableSynchronization\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mAUTO\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1555\u001b[0m aggregation=VariableAggregation.NONE):\n\u001b[0;32m-> 1556\u001b[0;31m return get_variable_scope().get_variable(\n\u001b[0m\u001b[1;32m 1557\u001b[0m \u001b[0m_get_default_variable_store\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1558\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/ops/variable_scope.py\u001b[0m in \u001b[0;36mget_variable\u001b[0;34m(self, var_store, name, shape, dtype, initializer, regularizer, reuse, trainable, collections, caching_device, partitioner, validate_shape, use_resource, custom_getter, constraint, synchronization, aggregation)\u001b[0m\n\u001b[1;32m 1297\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mdtype\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1298\u001b[0m \u001b[0mdtype\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_dtype\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1299\u001b[0;31m return var_store.get_variable(\n\u001b[0m\u001b[1;32m 1300\u001b[0m \u001b[0mfull_name\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1301\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/ops/variable_scope.py\u001b[0m in \u001b[0;36mget_variable\u001b[0;34m(self, name, shape, dtype, initializer, regularizer, reuse, trainable, collections, caching_device, partitioner, validate_shape, use_resource, custom_getter, constraint, synchronization, aggregation)\u001b[0m\n\u001b[1;32m 552\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mcustom_getter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0mcustom_getter_kwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 553\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 554\u001b[0;31m return _true_getter(\n\u001b[0m\u001b[1;32m 555\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 556\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/ops/variable_scope.py\u001b[0m in \u001b[0;36m_true_getter\u001b[0;34m(name, shape, dtype, initializer, regularizer, reuse, trainable, collections, caching_device, partitioner, validate_shape, use_resource, constraint, synchronization, aggregation)\u001b[0m\n\u001b[1;32m 505\u001b[0m \"name was already created with partitioning?\" % name)\n\u001b[1;32m 506\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 507\u001b[0;31m return self._get_single_variable(\n\u001b[0m\u001b[1;32m 508\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 509\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/ops/variable_scope.py\u001b[0m in \u001b[0;36m_get_single_variable\u001b[0;34m(self, name, shape, dtype, initializer, regularizer, partition_info, reuse, trainable, collections, caching_device, validate_shape, use_resource, constraint, synchronization, aggregation)\u001b[0m\n\u001b[1;32m 885\u001b[0m \u001b[0;31m# The code below handles only the case of creating a new variable.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 886\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mreuse\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 887\u001b[0;31m raise ValueError(\"Variable %s does not exist, or was not created with \"\n\u001b[0m\u001b[1;32m 888\u001b[0m \u001b[0;34m\"tf.get_variable(). Did you mean to set \"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 889\u001b[0m \"reuse=tf.AUTO_REUSE in VarScope?\" % name)\n",
|
||||
"\u001b[0;31mValueError\u001b[0m: Variable layer0/kernel does not exist, or was not created with tf.get_variable(). Did you mean to set reuse=tf.AUTO_REUSE in VarScope?"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# Boundary loss\n",
|
||||
"num_sample_points_bnd = 100\n",
|
||||
"x_bc, t_bc, u_bc = [np.concatenate([v_t0, v_x], axis=0) for v_t0, v_x in zip(boundary_tx(num_sample_points_bnd), open_boundary(num_sample_points_bnd))]\n",
|
||||
"x_bc, t_bc, u_bc = np.asarray(x_bc,dtype=np.float32), np.asarray(t_bc,dtype=np.float32) ,np.asarray(u_bc,dtype=np.float32)\n",
|
||||
"\n",
|
||||
"tf.compat.v1.disable_eager_execution()\n",
|
||||
"\n",
|
||||
"session = tf.compat.v1.Session(None)\n",
|
||||
"tf.compat.v1.initialize_all_variables()\n",
|
||||
"\n",
|
||||
"#with app.model_scope():\n",
|
||||
"loss_u = 0.5 * (network(x_bc, t_bc)[:, 0] - u_bc)**2 # normalizes by first dimension, N_bc\n",
|
||||
"\n",
|
||||
"# Physics loss inside of domain\n",
|
||||
"num_sample_points_inner = 1000\n",
|
||||
"x_ph, t_ph = tf.convert_to_tensor(rnd.uniform(-1, 1, [num_sample_points_inner])), tf.convert_to_tensor(rnd.uniform(0, 1, [num_sample_points_inner]))\n",
|
||||
"loss_ph = 0.5 * (f(network(x_ph, t_ph)[:, 0], x_ph, t_ph))**2 # normalizes by first dimension, N_ph\n",
|
||||
"\n",
|
||||
"# Combine\n",
|
||||
"ph_factor = 1.\n",
|
||||
"loss = loss_u + ph_factor * loss_ph # allows us to control the relative influence of loss_ph \n",
|
||||
"\n",
|
||||
"optim = tf.train.GradientDescentOptimizer(learning_rate=0.05).minimize(loss)\n",
|
||||
"#optim = tf.train.AdamOptimizer(learning_rate=0.001).minimize(loss) # alternative, but not much benefit here\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The code above just initializes the evaluation of the loss, we still didn't do any optimization steps, but we're finally in a good position to get started with this.\n",
|
||||
"\n",
|
||||
"Note: despite the simple equation, the convergence is typicaly very slow. It needs a _lot_ of iterations. To keep the runtime in a reasonable range, we only do 10k iterations by default below (`iters`). You can increase this value to get better results."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"session.initialize_variables()\n",
|
||||
"\n",
|
||||
"import time\n",
|
||||
"start = time.time()\n",
|
||||
"\n",
|
||||
"iters = 10000\n",
|
||||
"for optim_step in range(iters+1):\n",
|
||||
" _, loss_value = session.run([optim, loss])\n",
|
||||
" if optim_step<3 or optim_step%1000==0: \n",
|
||||
" print('Step %d, loss: %f' % (optim_step,loss_value))\n",
|
||||
" #show_state(grid_u)\n",
|
||||
" \n",
|
||||
"end = time.time()\n",
|
||||
"print(\"Runtime {:.2f}s\".format(end-start))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This training can take a significant amount of time, around 1 minute on a typical notebook, but at least the error goes down significantly (roughly from around 0.2 to ca. 0.03), and the network seems to successfully converge to a solution.\n",
|
||||
"\n",
|
||||
"Let's show the reconstruction of the network, by evaluating the network at the centers of a regular grid, so that we can show the solution as an image. Note that this is actually fairly expensive, we have to run through the whole network with a few thousand weights for all of the $128 \\times 32$ points in the grid.\n",
|
||||
"\n",
|
||||
"It looks pretty good on first sight, though. There's been a very noticeable change compared to the random initialization shown above:\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"show_state(session.run(grid_u)) "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"---\n",
|
||||
"\n",
|
||||
"## Evaluation\n",
|
||||
"\n",
|
||||
"Let's compare solution in a bit more detail. Here are the actual sample points used for constraining the solution (at time step 16, $t=1/2$) shown in gray, versus the reconstructed solution in blue:\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"u = session.run(grid_u)\n",
|
||||
"\n",
|
||||
"# solution is imposed at t=1/2 , which is 16 in the array\n",
|
||||
"bc_tx = 16 \n",
|
||||
"uT = u[0,:,bc_tx,0]\n",
|
||||
"\n",
|
||||
"fig = plt.figure().gca()\n",
|
||||
"fig.plot(np.linspace(-1,1,len(uT)), uT, lw=2, color='blue')\n",
|
||||
"fig.scatter(x_bc[0:100], u_bc[0:100], color='gray')\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Not too bad at the sides of the domain (the Dirichlet boundary conditions $u=0$ are fulfilled), but the shock in the center (at $x=0$) is not well represented.\n",
|
||||
"\n",
|
||||
"Let's check how well the initial state at $t=0$ was reconstructed. That's the most interesting, and toughest part of the problem (the rest basically follows from the model equation and boundary conditions, given the first state).\n",
|
||||
"\n",
|
||||
"It turns out, the accuracy of the initial state is actually not that great: the blue curve from the PINN is quite far away from the constraits (shown in gray)... The solution will get better with larger number of iterations, but it requires a surprisingly large number of them for a fairly simple case. \n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# ground truth solution at t0\n",
|
||||
"t0gt = np.asarray( [ [-np.sin(np.pi * x) * 1.] for x in np.linspace(-1,1,n)] )\n",
|
||||
"\n",
|
||||
"velP0 = u[0,:,0,0]\n",
|
||||
"\n",
|
||||
"fig = plt.figure().gca()\n",
|
||||
"fig.plot(np.linspace(-1,1,len(velP0)), velP0, lw=2, color='blue')\n",
|
||||
"fig.plot(np.linspace(-1,1,len(t0gt)), t0gt, lw=2, color='gray') # optionally show GT, compare to blue\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Especially the maximum / minimum at $x=\\pm 1/2$ are far off, and the boudaries at $x=\\pm 1$ are not fulfilled: the solution is not at zero.\n",
|
||||
"\n",
|
||||
"We have the forward simulator for this simulation, so we can use the $t=0$ solution of the network to \n",
|
||||
"evaluate how well the temporal evoluation was reconstructed by the PINN. This measures how well the temporal evolution of the model equation was captured via the soft constraints of the PINN loss.\n",
|
||||
"\n",
|
||||
"The graph below shows the initial state in blue, and two evolved states at $t=8/32$ and $t=15/32$. Note that this is all from the simulated version, we'll show the PINN version next. \n",
|
||||
"\n",
|
||||
"(Note: The code segments below also have some optional code to show the states at `[steps//4]`. It's commented out by default, you can uncomment or add additional ones to visualize more of the time evolution if you like.)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# re-simulate with phiflow from solution at t=0\n",
|
||||
"dt = 1./32.\n",
|
||||
"steps = 32-bc_tx # depends on where BCs were imposed\n",
|
||||
"initial = u[...,bc_tx:(bc_tx+1),0] # np.reshape(u0, [1,len(u0),1]) \n",
|
||||
"print(initial.shape)\n",
|
||||
"\n",
|
||||
"domain = Domain([n], boundaries=PERIODIC, box=box[-1:1])\n",
|
||||
"state = [BurgersVelocity(domain, velocity=initial, viscosity=0.01/np.pi)]\n",
|
||||
"physics = Burgers()\n",
|
||||
"\n",
|
||||
"for i in range(steps):\n",
|
||||
" state.append( physics.step(state[-1],dt=dt) )\n",
|
||||
"\n",
|
||||
"# we only need \"velocity.data\" from each phiflow state\n",
|
||||
"vel_resim = [x.velocity.data for x in state]\n",
|
||||
"\n",
|
||||
"fig = plt.figure().gca()\n",
|
||||
"pltx = np.linspace(-1,1,len(vels[0].flatten()))\n",
|
||||
"fig.plot(pltx, vel_resim[ 0].flatten(), lw=2, color='blue')\n",
|
||||
"#fig.plot(pltx, vel_resim[steps//4].flatten(), lw=2, color='green')\n",
|
||||
"fig.plot(pltx, vel_resim[steps//2].flatten(), lw=2, color='cyan')\n",
|
||||
"fig.plot(pltx, vel_resim[steps-1].flatten(), lw=2, color='purple')\n",
|
||||
"#fig.plot(pltx, t0gt, lw=2, color='gray') # optionally show GT, compare to blue\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"And here is the PINN output from `u` at the same time steps:\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"velP = [u[0,:,x,0] for x in range(33)]\n",
|
||||
"print(velP[0].shape)\n",
|
||||
"\n",
|
||||
"fig = plt.figure().gca()\n",
|
||||
"fig.plot(pltx, velP[bc_tx+ 0].flatten(), lw=2, color='blue')\n",
|
||||
"#fig.plot(pltx, velP[bc_tx+steps//4].flatten(), lw=2, color='green')\n",
|
||||
"fig.plot(pltx, velP[bc_tx+steps//2].flatten(), lw=2, color='cyan')\n",
|
||||
"fig.plot(pltx, velP[bc_tx+steps-1].flatten(), lw=2, color='purple')\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Not surprisingly, the errors grow over time. Especially the steepening of the solution near the shock at $x=0$ is not \"captured\" well. It's a bit difficult to see in these two graphs, though, let's quantify the error and show the actual difference:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"error = np.sum( np.abs( np.asarray(vel_resim[0:16]).flatten() - np.asarray(velP[bc_tx:bc_tx+steps]).flatten() )) / (steps*n)\n",
|
||||
"print(\"Mean absolute error for re-simulation across {} steps: {:7.5f}\".format(steps,error))\n",
|
||||
"\n",
|
||||
"fig = plt.figure().gca()\n",
|
||||
"fig.plot(pltx, (vel_resim[0 ].flatten()-velP[bc_tx ].flatten()), lw=2, color='blue')\n",
|
||||
"#fig.plot(pltx, (vel_resim[steps//4].flatten()-velP[bc_tx+steps//4].flatten()), lw=2, color='green')\n",
|
||||
"fig.plot(pltx, (vel_resim[steps//2].flatten()-velP[bc_tx+steps//2].flatten()), lw=2, color='cyan')\n",
|
||||
"fig.plot(pltx, (vel_resim[steps-1 ].flatten()-velP[bc_tx+steps-1 ].flatten()), lw=2, color='purple')\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This should show a mean absolute error of ca. $1.5 \\cdot 10^{-2}$ between ground truth re-simulation and the PINN evolution, which is significant for the value range of the simulation.\n",
|
||||
"\n",
|
||||
"And for comparison with the forward simulation and following cases, here are also all steps over time with a color map."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"scrolled": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# show re-simulated solution again as full image over time\n",
|
||||
"sn = np.concatenate(vel_resim, axis=-1)\n",
|
||||
"sn = np.reshape(sn, list(sn.shape)+[1] ) # print(sn.shape)\n",
|
||||
"show_state(sn)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Next, we'll store the full solution over the course of the $t=0 \\dots 1$ time interval, so that we can compare it later on to the full solution from a regular forward solve and compare it to the differential physics solution.\n",
|
||||
"\n",
|
||||
"Thus, stay tuned for the full evaluation and the comparison. It will follow in the file `diffphys-code-tf.ipynb`, after we've discussed the details of how to run the differential physics optimization."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"vels = session.run(grid_u) # special for showing NN results, run through TF \n",
|
||||
"vels = np.reshape( vels, [vels.shape[1],vels.shape[2]] )\n",
|
||||
"\n",
|
||||
"# save for comparison with other methods\n",
|
||||
"np.savez_compressed(\"./temp/burgers-pinn-solution.npz\",vels) ; print(\"Vels array shape: \"+format(vels.shape))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"---\n",
|
||||
"\n",
|
||||
"## Next steps\n",
|
||||
"\n",
|
||||
"This setup is just a starting point for PINNs and physical soft-constraints, of course. The parameters of the setup were chosen to compare to run relatively quickly. As we'll show in the next sections, the behavior of such an inverse solve can be improved substantially by a tighter integration of solver and learning. \n",
|
||||
"\n",
|
||||
"The solution of the PINN setup above can also directly be improved, however. E.g., try to:\n",
|
||||
"\n",
|
||||
"* Adjust parameters of the training to further decrease the error without making the solution diverge\n",
|
||||
"* Adapt the NN architecture for further improvements (keep track of the weight count, though!)\n",
|
||||
"* Activate a different optimizer, and observe the changing behavior (this typically requires adjusting the learning rate). Note that the more complex optimizers don't necessarily do better in this relatively simple example.\n",
|
||||
"* Or modify the setup to make the test case more interesting: e.g., move the boundary conditions further back in time, to let the reconstruction go \"further backward\" in time."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.8.5"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 4
|
||||
}
|
||||
+4
-28
@@ -6,7 +6,7 @@
|
||||
"source": [
|
||||
"# Burgers Optimization with a Physics-Informed NN\n",
|
||||
"\n",
|
||||
"Warning: this example still needs phiflow1.x and tensorflow 1.x! \n",
|
||||
"**Warning: this example still needs phiflow1.x and tensorflow 1.x! **\n",
|
||||
"\n",
|
||||
"To illustrate how the physics-informed losses work, let's consider a \n",
|
||||
"reconstruction example with a simple yet non-linear equation in 1D:\n",
|
||||
@@ -25,33 +25,9 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Could not load resample cuda libraries: CUDA binaries not found at /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/cuda/build/resample.so. Run \"python setup.py cuda\" to compile them\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/util.py:119: The name tf.AUTO_REUSE is deprecated. Please use tf.compat.v1.AUTO_REUSE instead.\n",
|
||||
"\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/profiling.py:12: The name tf.RunOptions is deprecated. Please use tf.compat.v1.RunOptions instead.\n",
|
||||
"\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/profiling.py:13: The name tf.RunMetadata is deprecated. Please use tf.compat.v1.RunMetadata instead.\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"/Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/flow.py:15: UserWarning: TensorFlow-CUDA solver is not available. To compile it, download phiflow sources and run\n",
|
||||
"$ python setup.py tf_cuda\n",
|
||||
"before reinstalling phiflow.\n",
|
||||
" warnings.warn(\"TensorFlow-CUDA solver is not available. To compile it, download phiflow sources and run\\n$ python setup.py tf_cuda\\nbefore reinstalling phiflow.\")\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from phi.tf.flow import *\n",
|
||||
"import numpy as np\n",
|
||||
@@ -777,7 +753,7 @@
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.7.6"
|
||||
"version": "3.8.5"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
|
||||
@@ -58,13 +58,13 @@ goals of the next sections.
|
||||
|
||||
✅ Pro:
|
||||
- uses physical model
|
||||
- derivatives via backpropagation
|
||||
- derivatives can be convieniently compute via backpropagation
|
||||
|
||||
❌ Con:
|
||||
- slow ...
|
||||
- only soft constraints
|
||||
- quite slow ...
|
||||
- physical constraints are enforced only as soft constraints
|
||||
- largely incompatible _classical_ numerical methods
|
||||
- derivatives rely on learned representation
|
||||
- accuracy of derivatives relies on learned representation
|
||||
|
||||
Next, let's look at how we can leverage numerical methods to improve the DL accuracy and efficiency
|
||||
by making use of differentiable solvers.
|
||||
Reference in new issue
Block a user