daviddoji/pbdl-book
1
0 Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

cleanup, working phiflow 2 version of DP NS opt

Browse Source
This commit is contained in:
NT committed 2021-03-10 12:14:33 +08:00
1 parent c443f2bfdf
commit 975d67d07a
1 file changed
+100 -115
Download Patch File Download Diff File Expand all files Collapse all files
diffphys-code-ns.ipynb
+100 -115
View File
@@ -8,9 +8,8 @@
"source": [
"# Differentiable Fluid Simulations\n",
"\n",
"<br>\n",
"\n",
"**TODO , NEXT finish opt case, diverges? see colab example , finish texts before end.**\n",
"<br>**TODO , NEXT finish opt case, diverges? see colab example , finish texts before end.**\n",
"<br>**TODO , fix bottom row, crop in display with [1:,0:]**\n",
"\n",
"<br><br>\n",
"\n",
@@ -70,22 +69,11 @@
},
{
"cell_type": "code",
"execution_count": 1,
"execution_count": null,
"metadata": {
"id": "da1uZcDXdVcF"
},
"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"
]
}
],
"outputs": [],
"source": [
"#!pip install --upgrade --quiet phiflow\n",
"#!pip install --upgrade --quiet git+https://github.com/tum-pbs/PhiFlow@develop\n",
@@ -154,7 +142,7 @@
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABHgAAAFtCAYAAAB80d3MAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8TklEQVR4nO3deZhcBZnv8d/b3dn3jZCQsCMCCgEjgrIJqMio4IKDjnMZ9QouuNyrznXUUcRxxg23kQsGwTDqqMgiMJdR2WRVIEBYQ0iAANnXTrqz9fbeP+pEm9j1nu46VafqdH8/z9NPuutX55w3p/u8fertU1Xm7gIAAAAAAEBxNdW7AAAAAAAAAGTDgAcAAAAAAKDgGPAAAAAAAAAUHAMeAAAAAACAgmPAAwAAAAAAUHAMeAAAAAAAAAqOAc8QZmbTzexOM2szs4vM7AIz+1kO273HzI6s9XZSajjczO6tZw3AUDPEe850M1tkZiPqWQcwlAzxnsN5DoaUIX6853qOsfu+zmOb6D8GPIOMmS0zs1P7efdzJa2XNN7dP13Dsv7MzN4qqc3dH06+foWZ/c7M1puZ93H/yWZ2nZltNbPnzey9A9jWfDPrMLP2Xh/NkuTuj0pqTeoBUKHB1HPMbISZXZ70mjYzW2hmbx7Ati4ws87des7+kuTuayTdrtI+AFChwdRzdlvuIDPbMZAHpJznYLAbjMe7mZ2dDGO2mtkzZnZ8P7dV1XOMZN9uT9azOuknY/u5eO77Gv3HgGdo20fSk+5e9oSjBj4s6ae9vu6UdJWkD5a5/8WSOiRNl/R3ki4xs8MGsL1vuvvYXh/dvbKfSzpvAOsCkE2j95wWSS9KOlHSBElflHSVme07gO39aree82yvjJ4D5KvRe05vF0t6oILtcZ4DlDT88W5mb5D0DUnvlzRO0gmSnu3rvmVU+xzjre4+VtIcSUdK+qd+LlfxvjazloEug4FhwDOImdk/mNndZvZtM9tkZs/t+mu0mc2XdI6kf0wmt381HTezt5nZE2bWamZ/MLNDktvfb2Y39rrfEjP7da+vXzSzOX2sb7ikkyXdses2d1/s7pdLeqKP+4+R9E5J/+zu7e5+t6QbJP19hbtkd3+QdEpelzMCg13Re467b3X3C9x9mbv3uPt/SXpO0qsq3ikvdZ+k/c1snyqtDxjSit5zei13tqRWSbcOdB+k+IM4z8EgMUiO969IutDd/5ScZ6xw9xUV7I6+VHyO4e6rJf1OpUGPJMnMjjGze5P99YiZnZTcPl+77WszazKzzyVXJG0ws6vMbHJy/33NzM3sg2b2gqTbkts/YKUrmTZZ6aqnfXpt283sw8n3otXMLjYz65V/KFm2zcyeNLOjkttnmtk1ZrYu+fn4xED3xWDAgGfwe42kxZKmSvqmpMvNzNz9H1Sa9O76y88tvRcys5dJ+oWkT0maJukmSTcmzewOSccnB/NMScMlHZsst7+ksZIe7aOWgyT1uPvyftb+Mkld7v50r9sekTSQK3g+amYbzexBM3tn7yBpqJ2SDh7A+gDEitxzXsLMpqvUh8o+MOvDW5Oe84SZfaR34O5dkpZKOqKSegD0qdA9x8zGS7pQ0v/u9//4pTjPwVBS2OPdSk+fnCtpmpktNbPlZvZDMxs1gP9/Tc4xzGyWpDcny8vM9pL0/yT9i6TJkj4j6Rozm1ZmX39c0pkqXQE9U9Imla5K7O1ESYdIepOZnSHp85LeodL34y6Vvj+9vUXSqyUdLundkt6U1HaWpAsk/Q9J4yW9TdIGM2uSdKNKjxX3knSKpE+Z2ZsGuj+KjgHP4Pe8u1+WXLJ7paQZKj3dKc3fSvp/7n6zu3dK+rakUZJem1wO2KbSlPcElSa+K83s5SodvHe5e08f65yYLNdfYyVt2e22zSpd0tgfP1Cp+e4h6Z8lzTez1+12n7akLgDVUeSe82dmNkylE5gr3f2pfi52lUonL9MkfUjSl8zsPbvdh54DVFfRe85XJV1e4SCa8xwMNUU+3qdLGibpXZKO11+eFvXFfi5fi3OM35hZm0pPT18r6cvJ7e+TdJO735RcaXSzpAWSTi+zng9L+oK7L3f3nSoNYN5lL3061gXJldLbk/v/m7svSgZT/yppzm5XH33d3Vvd/QWVXl9oTnL7/1RpuPSAlyx19+dVGgZNc/cL3b0j+b5eJunsAeyPQYHnwA1+q3d94u7bkqvb+vMCWjMlPd9r2R4ze1GliahUmnafJOnA5PNWlZrgsep1qeJuNqn/wxlJaldpMtvbePWzmbr7Q72+vMnMfq7SpPieXrePU6l2ANVR5J4jSUr+CvRTlV7/6/z+LufuT/b68l4z+75KJ3K9/ypFzwGqq7A9J3nax6kqPcgbMM5zMAQV9niXtD3599/dfZUkmdl3VBrwfCFt4RqdY5zp7reY2YmS/lOlK6NaVXqNnbPspS/SPkylQUtf9pF0nZn1HoR166XDtxd3u//37aXvwGUqfT92fZ9W98q26S/f59mSnilTw0wza+11W7NKVwcNKVzBg3JWqnSgSJKS5z3OlrTreaK7GuHxyed3qNQIT1T5Rrg0WdVeZfLdPS2pxcwO6nXbERrY0yV6c5Wah6Q/X344XKVLPQHUVyP0nF3bvVylk5J3Jn/pq9TuPadFpZPHRzKsE0B1NELPOUnSvpJeMLPVKj0N4p1m9lC0UIDzHKBvdT/e3X2TpOUqHad/vrnf/4M+VqkqnWO4+x2S5qt0ZZNUGsb81N0n9voY4+5fL7OKFyW9ebf7j9zt9YV8t/uft9v9R7n7vf0o90VJB5S5/bnd1jnO3ctddTRoMeBBOVdJ+hszOyV5qsKnJe2UtOvAu0PS6yWNSi4rvkvSaZKmSHq4rxW6e4ekW1RqlpJKXdHMRqp0AiIzG2nJiwG6+1ZJ10q60MzGJJcdn6Hk1ertLy/atW9f2zOzd5nZ2OQ5tW9U6XLDG3rd5URJtyWXEgKor7r3nMQlKl0C/dbkMuKXSHrOSX1tz8zOMLNJyTaOlvQJSdf3usvRkpYllxIDqK9G6DnzVHqgMif5uFSl173Y9VoTnOcA1dEIx7sk/UTSx81sDzObJOl/SfqvXstX7RzDzE6yMm/VXsb3JL3BzI6Q9DOVXu/nTWbWnPw/TrLSa/X05VJJX7PkKVZmNs1Kr7NTzqWS/smSd0Y2swlWem2d/vixpM+Y2auSfXFgst37JbWZ2f8xs1FJ3a8ws1f3c72DBgMe9MndF6t0ovDvktZLeqtKD3g6kvxplZ5CdVfy9RaV3ubvHn/pW3Tu7kd66btg7aPSJYu7rsrZrpf+pemjKj1Hdq1KlyB+xN133Xe2SpfxlXv1+U8mWaukb0n6kLv/oVf+dyo1GAB11gg9JzlBOE+lB1qrrfTuEO1m9ndJPlulp4g+VmZbZ6v0F702Sf8h6RvufmWvnJ4DNIhG6Dnuvs3dV+/6SLa3w93XJfflPAeogkY43hNflfSASs9SWKTS8OhrUk3OMWbrLwOsVEnf+Q9JX3L3F1X6o/rnJa1T6eqYz6r87OD7Kg2Xf2+l1/T5k0oviF1uW9ep9HbxvzSzLZIeV+lFnvtT569V2mf/qdK++I2kycn36S0qncM9p9L3+ceSJvRnvYOJ+cDfvh7IxMzukXS+u/c5ER/Aer4oaZ27/6iCZQ+X9CN3PzZLDQAaXxV7zvskHebu/1TBsnuo9BfCI919R5Y6ADQ2znOAoaNRzzHM7MeSfu3uv8tSF4qHAQ8AAAAAAEDB8RQtAAAAAACAgmPAAwAAAAAAUHAMeAAAAAAAAAquJc+NDfCt2gAMHuvdfVreG6XnAENWXXqORN8Bhip3t3psl54DDFl9nuvkOuApac5/kwDqrPv5+m2bngMMPfXsORJ9BxhqonfSzgM9Bxh6+j7X4SlaAAAAAAAABceABwAAAAAAoOAY8AAAAAAAABQcAx4AAAAAAICCY8ADAAAAAABQcAx4AAAAAAAACo4BDwAAAAAAQMEx4AEAAAAAACg4BjwAAAAAAAAFx4AHAAAAAACg4BjwAAAAAAAAFBwDHgAAAAAAgIJjwAMAAAAAAFBwDHgAAAAAAAAKjgEPAAAAAABAwTHgAQAAAAAAKDgGPAAAAAAAAAXHgAcAAAAAAKDgGPAAAAAAAAAUHAMeAAAAAACAgmPAAwAAAAAAUHAMeAAAAAAAAAqOAQ8AAAAAAEDBMeABAAAAAAAoOAY8AAAAAAAABceABwAAAAAAoOAY8AAAAAAAABQcAx4AAAAAAICCY8ADAAAAAABQcAx4AAAAAAAACo4BDwAAAAAAQMGlDnjMbKSZ3W9mj5jZE2b2leT2+Wb2nJktTD7m1LxaAIMePQdA3ug7APJEzwFQKy39uM9OSSe7e7uZDZN0t5n9d5J91t2vrl15AIYgeg6AvNF3AOSJngOgJlIHPO7uktqTL4clH17LogAMXfQcAHmj7wDIEz0HQK306zV4zKzZzBZKWivpZne/L4m+ZmaPmtl3zWxEmWXPNbMFZragOiUDGOzoOQDyRt8BkCd6DoBasNIAuZ93Npso6TpJH5e0QdJqScMlzZP0jLtfmLK8S80VFwugqLofdPe5A12KngOgMpX1HIm+A6AS3XJ3q2RJeg6AyvR9rjOgd9Fy91ZJt0s6zd1XeclOST+RdHRV6gSABD0HQN7oOwDyRM8BUE39eRetaclkWWY2StIbJD1lZjOS20zSmZIer12ZAIYKeg6AvNF3AOSJngOgVvrzLlozJF1pZs0qDYSucvf/MrPbzGyaJJO0UNKHa1cmgCGEngMgb/QdAHmi5wCoiQG9Bk/mjfEcUWCIqvz1MLKg5wBDVX16jkTfAYamyl+DJyt6DjBUVeE1eAAAAAAAANB4GPAAAAAAAAAUHAMeAAAAAACAgmPAAwAAAAAAUHAMeAAAAAAAAAqOAQ8AAAAAAEDBMeABAAAAAAAoOAY8AAAAAAAABceABwAAAAAAoOAY8AALine truncated
"image/png": "iVBORw0KGgoAAAANSUhEUgAABHgAAAFtCAYAAAB80d3MAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8TUlEQVR4nO3deZhcBZnv8d/b3dn3jZCQsCMCCgEjgrIJqMio4IKDjnMZ9QouuNyrznXUUcRxxg23kQsGwTDqqMgiMJdR2WRVIEBYQ0iAANnXTrqz9fbeP+pEm9j1nu46VafqdH8/z9NPuutX55w3p/u8fertU1Xm7gIAAAAAAEBxNdW7AAAAAAAAAGTDgAcAAAAAAKDgGPAAAAAAAAAUHAMeAAAAAACAgmPAAwAAAAAAUHAMeAAAAAAAAAqOAc8QZmbTzexOM2szs4vM7AIz+1kO273HzI6s9XZSajjczO6tZw3AUDPEe850M1tkZiPqWQcwlAzxnsN5DoaUIX6853qOsfu+zmOb6D8GPIOMmS0zs1P7efdzJa2XNN7dP13Dsv7MzN4qqc3dH06+foWZ/c7M1puZ93H/yWZ2nZltNbPnzey9A9jWfDPrMLP2Xh/NkuTuj0pqTeoBUKHB1HPMbISZXZ70mjYzW2hmbx7Ati4ws87des7+kuTuayTdrtI+AFChwdRzdlvuIDPbMZAHpJznYLAbjMe7mZ2dDGO2mtkzZnZ8P7dV1XOMZN9uT9azOuknY/u5eO77Gv3HgGdo20fSk+5e9oSjBj4s6ae9vu6UdJWkD5a5/8WSOiRNl/R3ki4xs8MGsL1vuvvYXh/dvbKfSzpvAOsCkE2j95wWSS9KOlHSBElflHSVme07gO39aree82yvjJ4D5KvRe05vF0t6oILtcZ4DlDT88W5mb5D0DUnvlzRO0gmSnu3rvmVU+xzjre4+VtIcSUdK+qd+LlfxvjazloEug4FhwDOImdk/mNndZvZtM9tkZs/t+mu0mc2XdI6kf0wmt381HTezt5nZE2bWamZ/MLNDktvfb2Y39rrfEjP7da+vXzSzOX2sb7ikkyXdses2d1/s7pdLeqKP+4+R9E5J/+zu7e5+t6QbJP19hbtkd3+QdEpelzMCg13Re467b3X3C9x9mbv3uPt/SXpO0qsq3ikvdZ+k/c1snyqtDxjSit5zei13tqRWSbcOdB+k+IM4z8EgMUiO969IutDd/5ScZ6xw9xUV7I6+VHyO4e6rJf1OpUGPJMnMjjGze5P99YiZnZTcPl+77WszazKzzyVXJG0ws6vMbHJy/33NzM3sg2b2gqTbkts/YKUrmTZZ6aqnfXpt283sw8n3otXMLjYz65V/KFm2zcyeNLOjkttnmtk1ZrYu+fn4xED3xWDAgGfwe42kxZKmSvqmpMvNzNz9H1Sa9O76y88tvRcys5dJ+oWkT0maJukmSTcmzewOSccnB/NMScMlHZsst7+ksZIe7aOWgyT1uPvyftb+Mkld7v50r9sekTSQK3g+amYbzexBM3tn7yBpqJ2SDh7A+gDEitxzXsLMpqvUh8o+MOvDW5Oe84SZfaR34O5dkpZKOqKSegD0qdA9x8zGS7pQ0v/u9//4pTjPwVBS2OPdSk+fnCtpmpktNbPlZvZDMxs1gP9/Tc4xzGyWpDcny8vM9pL0/yT9i6TJkj4j6Rozm1ZmX39c0pkqXQE9U9Imla5K7O1ESYdIepOZnSHp85LeodL34y6Vvj+9vUXSqyUdLundkt6U1HaWpAsk/Q9J4yW9TdIGM2uSdKNKjxX3knSKpE+Z2ZsGuj+KjgHP4Pe8u1+WXLJ7paQZKj3dKc3fSvp/7n6zu3dK+rakUZJem1wO2KbSlPcElSa+K83s5SodvHe5e08f65yYLNdfYyVt2e22zSpd0tgfP1Cp+e4h6Z8lzTez1+12n7akLgDVUeSe82dmNkylE5gr3f2pfi52lUonL9MkfUjSl8zsPbvdh54DVFfRe85XJV1e4SCa8xwMNUU+3qdLGibpXZKO11+eFvXFfi5fi3OM35hZm0pPT18r6cvJ7e+TdJO735RcaXSzpAWSTi+zng9L+oK7L3f3nSoNYN5lL3061gXJldLbk/v/m7svSgZT/yppzm5XH33d3Vvd/QWVXl9oTnL7/1RpuPSAlyx19+dVGgZNc/cL3b0j+b5eJunsAeyPQYHnwA1+q3d94u7bkqvb+vMCWjMlPd9r2R4ze1GliahUmnafJOnA5PNWlZrgsep1qeJuNqn/wxlJaldpMtvbePWzmbr7Q72+vMnMfq7SpPieXrePU6l2ANVR5J4jSUr+CvRTlV7/6/z+LufuT/b68l4z+75KJ3K9/ypFzwGqq7A9J3nax6kqPcgbMM5zMAQV9niXtD3599/dfZUkmdl3VBrwfCFt4RqdY5zp7reY2YmS/lOlK6NaVXqNnbPspS/SPkylQUtf9pF0nZn1HoR166XDtxd3u//37aXvwGUqfT92fZ9W98q26S/f59mSnilTw0wza+11W7NKVwcNKVzBg3JWqnSgSJKS5z3OlrTreaK7GuHxyed3qNQIT1T5Rrg0WdVeZfLdPS2pxcwO6nXbERrY0yV6c5Wah6Q/X344XKVLPQHUVyP0nF3bvVylk5J3Jn/pq9TuPadFpZPHRzKsE0B1NELPOUnSvpJeMLPVKj0N4p1m9lC0UIDzHKBvdT/e3X2TpOUqHad/vrnf/4M+VqkqnWO4+x2S5qt0ZZNUGsb81N0n9voY4+5fL7OKFyW9ebf7j9zt9YV8t/uft9v9R7n7vf0o90VJB5S5/bnd1jnO3ctddTRoMeBBOVdJ+hszOyV5qsKnJe2UtOvAu0PS6yWNSi4rvkvSaZKmSHq4rxW6e4ekW1RqlpJKXdHMRqp0AiIzG2nJiwG6+1ZJ10q60MzGJJcdn6Hk1ertLy/atW9f2zOzd5nZ2OQ5tW9U6XLDG3rd5URJtyWXEgKor7r3nMQlKl0C/dbkMuKXSHrOSX1tz8zOMLNJyTaOlvQJSdf3usvRkpYllxIDqK9G6DnzVHqgMif5uFSl173Y9VoTnOcA1dEIx7sk/UTSx81sDzObJOl/SfqvXstX7RzDzE6yMm/VXsb3JL3BzI6Q9DOVXu/nTWbWnPw/TrLSa/X05VJJX7PkKVZmNs1Kr7NTzqWS/smSd0Y2swlWem2d/vixpM+Y2auSfXFgst37JbWZ2f8xs1FJ3a8ws1f3c72DBgMe9MndF6t0ovDvktZLeqtKD3g6kvxplZ5CdVfy9RaV3ubvHn/pW3Tu7kd66btg7aPSJYu7rsrZrpf+pemjKj1Hdq1KlyB+xN133Xe2SpfxlXv1+U8mWaukb0n6kLv/oVf+dyo1GAB11gg9JzlBOE+lB1qrrfTuEO1m9ndJPlulp4g+VmZbZ6v0F702Sf8h6RvufmWvnJ4DNIhG6Dnuvs3dV+/6SLa3w93XJfflPAeogkY43hNflfSASs9SWKTS8OhrUk3OMWbrLwOsVEnf+Q9JX3L3F1X6o/rnJa1T6eqYz6r87OD7Kg2Xf2+l1/T5k0oviF1uW9ep9HbxvzSzLZIeV+lFnvtT569V2mf/qdK++I2kycn36S0qncM9p9L3+ceSJvRnvYOJ+cDfvh7IxMzukXS+u/c5ER/Aer4oaZ27/6iCZQ+X9CN3PzZLDQAaXxV7zvskHebu/1TBsnuo9BfCI919R5Y6ADQ2znOAoaNRzzHM7MeSfu3uv8tSF4qHAQ8AAAAAAEDB8RQtAAAAAACAgmPAAwAAAAAAUHAMeAAAAAAAAAquJc+NDfCt2gAMHuvdfVreG6XnAENWXXqORN8Bhip3t3psl54DDFl9nuvkOuApac5/kwDqrPv5+m2bngMMPfXsORJ9BxhqonfSzgM9Bxh6+j7X4SlaAAAAAAAABceABwAAAAAAoOAY8AAAAAAAABQcAx4AAAAAAICCY8ADAAAAAABQcAx4AAAAAAAACo4BDwAAAAAAQMEx4AEAAAAAACg4BjwAAAAAAAAFx4AHAAAAAACg4BjwAAAAAAAAFBwDHgAAAAAAgIJjwAMAAAAAAFBwDHgAAAAAAAAKjgEPAAAAAABAwTHgAQAAAAAAKDgGPAAAAAAAAAXHgAcAAAAAAKDgGPAAAAAAAAAUHAMeAAAAAACAgmPAAwAAAAAAUHAMeAAAAAAAAAqOAQ8AAAAAAEDBMeABAAAAAAAoOAY8AAAAAAAABceABwAAAAAAoOAY8AAAAAAAABQcAx4AAAAAAICCY8ADAAAAAABQcAx4AAAAAAAACo4BDwAAAAAAQMGlDnjMbKSZ3W9mj5jZE2b2leT2+Wb2nJktTD7m1LxaAIMePQdA3ug7APJEzwFQKy39uM9OSSe7e7uZDZN0t5n9d5J91t2vrl15AIYgeg6AvNF3AOSJngOgJlIHPO7uktqTL4clH17LogAMXfQcAHmj7wDIEz0HQK306zV4zKzZzBZKWivpZne/L4m+ZmaPmtl3zWxEmWXPNbMFZragOiUDGOzoOQDyRt8BkCd6DoBasNIAuZ93Npso6TpJH5e0QdJqScMlzZP0jLtfmLK8S80VFwugqLofdPe5A12KngOgMpX1HIm+A6AS3XJ3q2RJeg6AyvR9rjOgd9Fy91ZJt0s6zd1XeclOST+RdHRV6gSABD0HQN7oOwDyRM8BUE39eRetaclkWWY2StIbJD1lZjOS20zSmZIer12ZAIYKeg6AvNF3AOSJngOgVvrzLlozJF1pZs0qDYSucvf/MrPbzGyaJJO0UNKHa1cmgCGEngMgb/QdAHmi5wCoiQG9Bk/mjfEcUWCIqvz1MLKg5wBDVX16jkTfAYamyl+DJyt6DjBUVeE1eAAAAAAAANB4GPAAAAAAAAAUHAMeAAAAAACAgmPAAwAAAAAAUHAMeAAAAAAAAAqOAQ8AAAAAAEDBMeABAAAAAAAoOAY8AAAAAAAABceABwAAAAAAoOAY8AALine truncated
"text/plain": [
"<Figure size 1152x360 with 4 Axes>"
]
@@ -226,34 +214,18 @@
}
],
"source": [
"#smoke = initial_smoke = DOMAIN.grid(math.zeros(inflow_loc=4))\n",
"#velocity = initial_velocity = DOMAIN.staggered_grid(0) * math.ones(inflow_loc=4)\n",
"initial_smoke = DOMAIN.grid(math.zeros(inflow_loc=4))\n",
"initial_velocity = DOMAIN.staggered_grid(0) * math.ones(inflow_loc=4)\n",
"\n",
"print(\"Velocity tensor shape: \"+format(initial_velocity.shape))\n",
"\n",
"# old?\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",
"\n",
"def simulate(smoke: CenteredGrid, velocity: StaggeredGrid):\n",
" for _ in range(20):\n",
" smoke,velocity = step(smoke,velocity)\n",
" \n",
" loss = field.l2_loss(diffuse.explicit(smoke - field.stop_gradient(smoke.inflow_loc[-1]), 1, 1, 10))\n",
" \n",
" return loss, smoke, velocity\n",
"# smoke = advect.mac_cormack(smoke, velocity, dt=1) + INFLOW\n",
"# buoyancy_force = smoke * (0, 0.5) >> velocity\n",
"# velocity = advect.semi_lagrangian(velocity, velocity, dt=1) + buoyancy_force\n",
"# velocity, _, _, _ = fluid.make_incompressible(velocity, DOMAIN)\n",
"\n",
"\n"
" return loss, smoke, velocity\n"
]
},
{
@@ -262,7 +234,7 @@
"source": [
"Then, we use `field.functional_gradient()` to obtain the gradient with respect to the initial velocity. Since the velocity is the second argument of the `simulate()` function, we pass `wrt=[1]`. (Note that there's also a `field.spatal_gradient` with computes derivatives of tensors along spatial dimensions like `x,y`.)\n",
"\n",
"`functional_gradient` returns a gradient function, which we can evaluate with the initial states for smoke and velocity."
"`functional_gradient` returns a gradient function, which we can evaluate with the initial states for smoke and velocity. \n"
]
},
{
@@ -274,45 +246,35 @@
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 0.0000000e+00 0.0000000e+00 0.0000000e+00 ... 0.0000000e+00\n",
" 0.0000000e+00 0.0000000e+00]\n",
" [-1.2200639e-02 -1.9027891e-02 -2.3449246e-02 ... 0.0000000e+00\n",
" 0.0000000e+00 0.0000000e+00]\n",
" [-6.2169664e-02 -6.8051241e-02 -6.4626858e-02 ... -3.4051553e-43\n",
" 0.0000000e+00 0.0000000e+00]\n",
" ...\n",
" [ 2.0202854e-15 5.0669552e-14 3.2829501e-13 ... 0.0000000e+00\n",
" 0.0000000e+00 0.0000000e+00]\n",
" [ 1.0069069e-18 4.1775261e-17 4.3741229e-16 ... 0.0000000e+00\n",
" 0.0000000e+00 0.0000000e+00]\n",
" [ 1.9205262e-24 7.9680139e-23 8.3429929e-22 ... 0.0000000e+00\n",
" 0.0000000e+00 0.0000000e+00]]\n"
"Gradient info StaggeredGrid[(inflow_loc=4, x=32, y=40, vector=2), size=(32, 40) along vector, extrapolation=0]\n"
]
},
{
"data": {
"text/plain": [
"(x=33, y=40) float32 -1.1465684 < ... < 1.2370852"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sim_grad = field.functional_gradient(simulate, wrt=[1], get_output=False)\n",
"[velocity_grad] = sim_grad(initial_smoke, initial_velocity)\n",
"\n",
"for i in range(1):\n",
" sim_grad = field.functional_gradient(simulate, wrt=[1], get_output=True)\n",
" #[vals,velocity_grad]\n",
" [loss,s,v,velocity_grad] = sim_grad(initial_smoke, initial_velocity)\n",
"\n",
" print(format( velocity_grad.inflow_loc[0].vector[0].values.numpy('x,y') ))\n",
"\n",
"#v2 = initial_velocity + velocity_grad\n",
"#[loss2,s,v,velocity_grad] = sim_grad(initial_smoke, initial_velocity+velocity_grad)\n",
"\n"
"print(\"Gradient info \" + format(velocity_grad))\n",
"velocity_grad.inflow_loc[0].vector[0].values # one example, location 0, x component\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",
"The last two lines just print some information about the resulting gradient field. Naturally, it has the same shape as the velocity itself: it's a staggered grid with four inflow locations. The last line shows how to access the x-components of one of the gradients.\n",
"\n",
"Let's direcly take a look at it. The following images show the `x` components (via `vector[0].values`) of the staggered grid gradient in `grad`.\n",
"\n",
"**TODO , fix bottom row, crop in display with [1:,0:]**"
"We can directly use this to take a look at the content of the computed gradient. The following images show the `x` components (via `vector[0].values.numpy()`) of the staggered grid gradient in `velocity_grad`.\n"
]
},
{
@@ -324,12 +286,12 @@
"name": "stdout",
"output_type": "stream",
"text": [
"CenteredGrid[(x=33, y=40), size=(33.0, 40.0) along vector, extrapolation=boundary]\n"
"CenteredGrid[(x=33, y=40), size=(33.0, 40.0) along vector, extrapolation=0]\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABHUAAAFgCAYAAAA1sUWiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABYEUlEQVR4nO3de7hkVX3n//e3L4CICNgttlxsjSQGnYixg/ozMUTRYMaIM1EEjcEZkBjDJJlkZoKaqEHNYDKJSSaOsaMImqgwGENHSQiijLko0ihegBhbAtLIrbkIyM2G7++PvU939em9+1R11dp1e7+ep55TtWrX3mvXqfr06nXWXisyE0mSJEmSJE2XZeOugCRJkiRJkgZnp44kSZIkSdIUslNHkiRJkiRpCtmpI0mSJEmSNIXs1JEkSZIkSZpCdupIkiRJkiRNITt1JEmSJEmSekTEMRHxjYjYFBGnNTz/vIj4UkRsjYiXL3ruxIj4Zn07saf8mRHxtXqffxIRMWw97dSRJEmSJEmqRcRy4D3Ai4HDgRMi4vBFm30beC3wkUWvPQB4K/As4EjgrRGxf/30e4HXAYfVt2OGraudOpIkSZIkSdsdCWzKzGsy80HgY8CxvRtk5rWZ+VXg4UWv/Wngosy8PTPvAC4CjomINcC+mfmFzEzgQ8DLhq3oimF3IGk0jonILUPu43K4MDOH7u2VNLuGzRpzRtJSzBlJpY0gZ64E7u8pWp+Z63seHwRc3/N4M9XIm340vfag+ra5oXwodupIE2ILsHHIfQSsGkVdJM2uYbPGnJG0FHNGUmkjyJn7M3PdqOozTnbqSJNk2ZBXRD68eOSfJDUYJmvMGUn9MGcklVY2Z24ADul5fHBd1o8bgKMWvfaSuvzg3dxnKzt1pElip46kLvifLUmlmTOSSiubM5cBh0XEE6k6Xo4HXtXn3i8EfrdncuQXAW/MzNsj4q6IeDZwKfALwP8euO6LOFGyJEmSJElSLTO3AqdSddBcDZybmVdGxOkR8VKAiPixiNgMvAJ4X0RcWb/2duDtVB1DlwGn12UAbwDeD2wCvgX87bB1daSONCkihh+pI0lLMWsklWbOSCqtg5zJzAuACxaVvaXn/mXseDlV73ZnAmc2lG8EnjbKetqpI00SG0CSumDWSCrNnJFUmjkD2KkjTQ7/qiWpC2aNpNLMGUmlmTPb+C5IkiRJkiRNIUfqSJPE3mZJXTBrJJVmzkgqzZwB7NSRJovBJKkLZo2k0swZSaWZM4CdOtLk8LpQSV0waySVZs5IKs2c2cZOHWmSGEySumDWSCrNnJFUmjkDOFGyJEmSJEnSVHKkjjQpHEIoqQtmjaTSzBlJpZkz29ipI00Sg0lSF8waSaWZM5JKM2cAO3WkyWIwSeqCWSOpNHNGUmnmDOCcOpIkSZIkSVPJkTrSpPC6UEldMGsklWbOSCrNnNnGTh1pkhhMkrpg1kgqrXDORMQxwB8Dy4H3Z+YZi57/deBkYCtwK/CfM/O6opWS1C3bM4CdOtLksLdZUhfMGkmlFc6ZiFgOvAd4IbAZuCwiNmTmVT2bfRlYl5n3RsQvAb8HvLJYpSR1y/bMNnbqSJPEYJLUBbNGUmllc+ZIYFNmXgMQER8DjgW2depk5md7tv8C8PMlKyRpDGzPAE6ULEmSJGmyrIqIjT23UxY9fxBwfc/jzXVZm5OAvx11JSVpEjhSR5ok9jZL6oJZI6m04XJmS2auG0U1IuLngXXAT45if5ImiO0ZwE4daXJ4XaikLnSQNRGxH/B+4GlAUk1Q+vmiB5U0OcrnzA3AIT2PD67LFlUjjgbeDPxkZj5QskKSOub/nbaxU0eaJAaTpC6Uz5o/Bv4uM18eEXsAe5c+oKQJUzZnLgMOi4gnUnXmHA+8qneDiHgG8D7gmMy8pWRlJI2J/3cC7NSRJEkjFBGPBp4HvBYgMx8EHhxnnSTNlszcGhGnAhdSLWl+ZmZeGRGnAxszcwPw+8A+wP+NCIBvZ+ZLx1ZpSSrETh1pUjiEUFIXymfNE4FbgQ9GxNOBy4FfzczvlTyopAnSQZsmMy8ALlhU9pae+0cXrYCk8fL/Ttv4LkiTZNmy4W6S1I/hcmapVWlWAD8KvDcznwF8Dzit4zOUNG62ZySVZs4AjtSRJsuMBYykCVV2VZrNwObMvLR+fB526kjzxzaNpNLMGcBOHWlydLMizZnAS4BbMvNpDc8H1QSnPwPcC7w2M79UtFKSulU4azLzpoi4PiJ+KDO/AbwAuKrYASVNHi+LkFSaObON74I0X84CjtnF8y8GDqtvpwDv7aBOkmbPfwH+MiK+ChwB/O54qyNJkjSbHKkjTZLykwp+LiLW7mKTY4EPZWYCX4iI/SJiTWbeWLRikrpVPmuuAHZ1iZakWedf0CWVZs4AdupIk2M0QwhXRcTGnsfrM3P9AK8/CLi+5/HmusxOHWlWOFxZUmnmjKTSzJlt7NSRJsnwwbTUBKaSZCNIUnnmjKTSzBlgjufUiYhrI+Lo+v6bIuL9HR03IuKDEXFHRHyxi2OOQkScFRHvGOPx/ykinjGu49d1+JGI+Odx1qEDNwCH9Dw+uC7TbjBnBmPOQEQcGBFXR8Se46yHpoc5MxhzZm7aMxoxs2YwZo1tmi5NZKdORBwfEZdGxPci4pb6/hvqlXlGLjN/NzNPHnY/EbE2IjIidjUC6seBFwIHZ+aRwx5zGvT5vuzq9T8L3J2ZX64fPy0iLoyILRGRDdsfEBGfqD8/10XEqwY41lkR8WBE3NNzWw6QmV8F7qzrU8ayZcPdhrcB+IX6H9BnA9+d1fl0zJnZ0mXORMSeEfGBOl/ujogrIuLFAxzrbRHx/UU58ySAzLwZ+CzVROXljDdn5oY5M1u6bs/0vO6wiLg/Iv5igGONtz0D5kyHzJrZMo6sqT9DV9efoW9FxE/0eazxtmnMGWACO3Ui4jeollT+feBxwIHA64HnAnu0vGZ5ZxUc3hOAazPze4O+cHe/2DPg9cCHex5/HzgXOKll+/cAD1J9dl4NvDcinjrA8X4vM/fpuT3U89xfAr84wL76t3BdaMFwioiPAp8HfigiNkfESRHx+oh4fb3JBcA1wCbgz4E3FDnXMTNn2pkz2+wqZ1ZQzT31k8Cjgd8Czo1dT0K+2DmLcuaanufK5QwMnzXqiznTzpzZZqn2zIL3AJftxvHG054Bc6ZDZk07s2abXWZNRLwQeBfwn4BHAc+j+v9Av8bTpjFntpmos4mIRwOnA2/IzPMy8+6sfDkzX52ZD9TbnRUR742ICyLie8BPRcS/j4gvR8RdEXF9RLxt0b5fU/9V9baIePOi597W+9ePiHh2RPxzRNwZEV+JiKN6nrskIt4e1ZC2uyPi7yNiVf305+qfd9a9lM9ZdJyTgPcDz6mf/526/HURsSkibo+IDRHx+J7XZET8ckR8E/hmw3v2yoj4t4jYt3784oi4KSJWN2z7txFx6qKyr0TEf6zvPyUiLqrr8Y2IOG6nX9II1L+/90TEp+r38NKI+IGWbfcAng/8v4WyzPxGZn4AuLJh+0cCPwf8dmbek5n/SDX65DUjqv4lwAui1DDCwp06mXlCZq7JzJWZeXBmfiAz/ywz/6x+PjPzlzPzBzLz32XmxqX2OW3MGXOmYduBciYzv5eZb8vMazPz4cz8JPBvwDNHVP1LgSdFxBNGtL+d2QgqypwxZxq2HShnel53PHAncPFoa1+4PQPmTAfMGrOmYdvdyZrfAU7PzC/U7ZobMnNU0y+UbdMUzpmIOKb+3W6KiNMant8zIs6pn7806j/wRcSroxrJvXB7OCKOqJ+7pN7nwnOPHfptGHYHI/YcYE/g/D62fRXwTqrexH8Evgf8ArAf8O+BX4qIlwFExOHAe6n+Y/944DFUc4XsJCIOAj4FvAM4APhvwMcXfdFfRdWT+ViqHvD/Vpc/r/65X91L+fnefddfptcDn6+ff2tEPB/4n8BxwBrgOuBji6r1MuBZwOGL65uZ5wD/DPxJRDwG+ABwcmbe2nB6HwVO6DnXw6l6vz8VVWfIRcBH6vM6Hvg/9TYlHE8VIPtTjQp5Z8t2hwEPZ+bmPvf7g8DWzPzXnrKvAIOM1HlDHc6XR8TP9T5RB9z3gR8aYH/96WCkjgBzxpzZ2aA5s4OIOJAqe1r/Y9bgZ+ucuTIifqn3iczcWtf36btTnyX5l60umDPmzGID50z9n87TgV/fzbqNpz0D5kx3zBqzZrGBsiaqUVvrgNV1x8TmiPjTiHjEAHUbT5umcM7U7817gBdTfZZOaPj9ngTckZlPBt5NNeKJzPzLzDwiM4+g+h79W2Ze0fO6Vy88n5m3DPtWTFpqrqJavWfrQkFPr+99EfG8nm3Pz8x/qnsT78/MSzLza/Xjr1J9CX+y3vblwCcz83N1j/VvAw+31OHngQsy84J6XxcBG4Gf6dnmg5n5r5l5H9VQtiOGOOdXA2dm5pfqur2Rqjd6bc82/zMzb6+P1+SXqXpkLwH+pv6rcZNPAEf09JS+Gvir+rgvoRra+MHM3JrVNZgfB14xxLntyicy84v17/ovaX8P9wPuHmC/+wB3LSr7LtU/YP34E6owfCzV5+SsiHjuom3uruul6WTOmDOL7cdgObNNRKys9312Zv5Lny87F/hhYDXwOuAtEXHCom3Mmelmzpgzi+3H4DnzduADu9nhbHtmPpg1Zs1i+zFY1hwIrKT6nf9Evd9nUF1a3o9ZbtMcCWzKzGsy80GqzsNjF21zLHB2ff88qhGQi+eyOoGdOx5HatI6dW4DVkXP9Y+Z+f9l5n71c731vb73hRHxrIj4bETcGhHfperVXRja9/je7bO6JvO2ljo8AXhFHYZ3RsSdVBN0renZ5qae+/dSdSTsrsdT9TAv1O2eum4H9Wxz/eIX9crMO4H/CzwN+INdbHc3VU/68XXRCVShANV5P2vReb+a6trcXap7ZRcmxuprUi36fw/voP8OGYBLine truncated
"image/png": "iVBORw0KGgoAAAANSUhEUgAABHUAAAFgCAYAAAA1sUWiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABYE0lEQVR4nO3de7hkVX3n//e3L4CICNgttlxsjSQGnYixg/ozMUTRYMaIM1EEL8EZkBjDJJlkZoKaqCEmg8kkJpk4xo4iaLzAYAwdJSGIMuZiCI3iBYixJSCN3JqLgFwbvr8/9j7d1af37lPVVWvX7f16nnpO1apde699TtWnV69ae63ITCRJkiRJkjRdlo27ApIkSZIkSRqcnTqSJEmSJElTyE4dSZIkSZKkKWSnjiRJkiRJ0hSyU0eSJEmSJGkK2akjSZIkSZI0hezUkSRJkiRJ6hERx0TENyJiU0Sc1vD8CyLiSxGxNSJeuei5EyPim/XtxJ7yZ0fE1+p9/nFExLD1tFNHkiRJkiSpFhHLgfcCLwUOB06IiMMXbfZt4A3Axxa99gDgHcBzgCOBd0TE/vXT7wPeCBxW344Ztq526kiSJEmSJG13JLApM6/JzAeBTwDH9m6Qmddm5leBRxa99ieBizLz9sy8A7gIOCYi1gD7ZuY/ZWYCHwZeMWxFVwy7A0mjcUxEbhlyH5fDhZk5dG+vpNk1bNaYM5KWYs5IKm0EOXMlcH9P0frMXN/z+CDg+p7Hm6lG3vSj6bUH1bfNDeVDsVNHmhBbgI1D7iNg1SjqIml2DZs15oykpZgzkkobQc7cn5nrRlWfcbJTR5oky4a8IvKRxSP/JKnBMFljzkjqhzkjqbSyOXMDcEjP44Prsn7cABy16LWX1OUH7+Y+W9mpI00SO3UkdcH/bEkqzZyRVFrZnLkMOCwinkzV8XI88Jo+934h8Ds9kyO/BHhLZt4eEXdFxHOBS4GfBf73wHVfxImSJUmSJEmSapm5FTiVqoPmauDczLwyIk6PiJcDRMSPRMRm4FXA+yPiyvq1twO/RdUxdBlwel0G8GbgA8Am4FvAXw9bV0fqSJMiYviROpK0FLNGUmnmjKTSOsiZzLwAuGBR2dt77l/GjpdT9W53JnBmQ/lG4BmjrKedOtIksQEkqQtmjaTSzBlJpZkzgJ060uTwWy1JXTBrJJVmzkgqzZzZxt+CJEmSJEnSFHKkjjRJ7G2W1AWzRlJp5oyk0swZwE4dabIYTJK6YNZIKs2ckVSaOQPYqSNNDq8LldQFs0ZSaeaMpNLMmW3s1JEmicEkqQtmjaTSzBlJpZkzgBMlS5IkSZIkTSVH6kiTwiGEkrpg1kgqzZyRVJo5s42dOtIkMZgkdcGskVSaOSOpNHMGsFNHmiwGk6QumDWSSjNnJJVmzgDOqSNJkiRJkjSVHKkjTQqvC5XUBbNGUmnmjKTSzJlt7NSRJonBJKkLZo2k0grnTEQcA/wRsBz4QGaesej5XwFOBrYCtwL/OTOvK1opSd2yPQPYqSNNDnubJXXBrJFUWuGciYjlwHuBFwObgcsiYkNmXtWz2ZeBdZl5b0T8PPC7wKuLVUpSt2zPbGOnjjRJDCZJXTBrJJVWNmeOBDZl5jUAEfEJ4FhgW6dOZn6+Z/t/Al5XskKSxsD2DOBEyZIkSZImy6qI2NhzO2XR8wcB1/c83lyXtTkJ+OtRV1KSJoEjdaRJYm+zpC6YNZJKGy5ntmTmulFUIyJeB6wDfnwU+5M0QWzPAHbqSJPD60IldaGDrImI/YAPAM8AkmqC0i8WPaikyVE+Z24ADul5fHBdtqgacTTwNuDHM/OBkhWS1DH/77SNnTrSJDGYJHWhfNb8EfA3mfnKiNgD2Lv0ASVNmLI5cxlwWEQ8maoz53jgNb0bRMSzgPcDx2TmLSUrI2lM/L8TYKeOJEkaoYh4LPAC4A0Amfkg8OA46yRptmTm1og4FbiQaknzMzPzyog4HdiYmRuA3wP2Af5vRAB8OzNfPrZKS1IhdupIk8IhhJK6UD5rngzcCnwoIp4JXA78UmZ+r+RBJU2QDto0mXkBcMGisrf33D+6aAUkjZf/d9rG34I0SZYtG+4mSf0YLmeWWpVmBfDDwPsy81nA94DTOj5DSeNme0ZSaeYM4EgdabLMWMBImlBlV6XZDGzOzEvrx+dhp440f2zTSCrNnAHs1JEmRzcr0pwJvAy4JTOf0fB8UE1w+lPAvcAbMvNLRSslqVuFsyYzb4qI6yPiBzLzG8CLgKuKHVDS5PGyCEmlmTPb+FuQ5stZwDG7eP6lwGH17RTgfR3USdLs+S/ARyPiq8ARwO+MtzqSJEmzyZE60iQpP6ngFyJi7S42ORb4cGYm8E8RsV9ErMnMG4tWTFK3ymfNFcCuLtGSNOv8Bl1SaeYMYKeONDlGM4RwVURs7Hm8PjPXD/D6g4Drex5vrsvs1JFmhcOVJZVmzkgqzZzZxk4daZIMH0xLTWAqSTaCJJVnzkgqzZwB5nhOnYi4NiKOru+/NSI+0NFxIyI+FBF3RMQ/d3HMUYiIsyLiXWM8/j9ExLPGdfy6Dj8UEf84zjp04AbgkJ7HB9dl2g3mzGDMGYiIAyPi6ojYc5z10PQwZwZjzsxNe0YjZtYMxqyxTdOliezUiYjjI+LSiPheRNxS339zvTLPyGXm72TmycPuJyLWRkRGxK5GQP0o8GLg4Mw8cthjToM+fy+7ev1PA3dn5pfrx8+IiAsjYktEZMP2B0TEp+r3z3UR8ZoBjnVWRDwYEff03JYDZOZXgTvr+pSxbNlwt+FtAH62/gf0ucB3Z3U+HXNmtnSZMxGxZ0R8sM6XuyPiioh46QDHemdEPLQoZ54CkJk3A5+nmqi8nPHmzNwwZ2ZL1+2ZntcdFhH3R8SfD3Cs8bZnwJzpkFkzW8aRNfV76Or6PfStiPixPo813jaNOQNMYKdORPwq1ZLKvwc8ATgQeBPwfGCPltcs76yCw3sScG1mfm/QF+7uB3sGvAn4SM/jh4BzgZNatn8v8CDVe+e1wPsi4ukDHO93M3OfntvDPc99FPi5AfbVv4XrQguGU0R8HPgi8AMRsTkiToqIN0XEm+pNLgCuATYBfwa8uci5jpk5086c2WZXObOCau6pHwceC/w6cG7sehLyxc5ZlDPX9DxXLmdg+KxRX8yZdubMNku1Zxa8F7hsN443nvYMmDMdMmvamTXb7DJrIuLFwLuB/wQ8BngB1f8H+jWeNo05s81EnU1EPBY4HXhzZp6XmXdn5cuZ+drMfKDe7qyIeF9EXBAR3wN+IiL+fUR8OSLuiojrI+Kdi/b9+vpb1dsi4m2Lnntn77cfEfHciPjHiLgzIr4SEUf1PHdJRPxWVEPa7o6Iv42IVfXTX6h/3ln3Uj5v0XFOAj4APK9+/jfr8jdGxKaIuD0iNkTEE3tekxHxCxHxTeCbDb+zV0fEv0XEvvXjl0bETRGxumHbv46IUxeVfSUi/mN9/2kRcVFdj29ExHE7/ZFGoP77vTciPlP/Di+NiO9r2XYP4IXA/1soy8xvZOYHgSsbtn808DPAb2TmPZn591SjT14/oupfArwoSg0jLNypk5knZOaazFyZmQdn5gcz808z80/r5zMzfyEzvy8z/11mblxqn9PGnDFnGrYdKGcy83uZ+c7MvDYzH8nMTwP/Bjx7RNW/FHhKRDxpRPvbmY2goswZc6Zh24Fypud1xwN3AhePtvaF2zNgznTArDFrGrbdnaz5TeD0zPynul1zQ2aOavqFsm2awjkTEcfUf9tNEXFaw/N7RsQ59fOXRv0FX0S8NqqR3Au3RyLiiPq5S+p9Ljz3+KF/DcPuYMSeB+wJnN/Htq8BfpuqN/Hvge8BPwvsB/x74Ocj4hUAEXE48D6q/9g/EXgc1VwhO4mIg4DPAO8CDgD+G/DJRR/011D1ZD6eqgf8v9XlL6h/7lf3Un6xd9/1h+lNwBfr598RES8E/idwHLAGuA74xKJqvQJ4DnD44vpm5jnAPwJ/HBGPAz4InJyZtzac3seBE3rO9XCq3u/PRNUZchHwsfq8jgf+T71NCcdTBcj+VKNCfrtlu8OARzJzc5/7/X5ga2b+a0/ZV4BBRuq8uQ7nyyPiZ3qfqAPuIeAHBthffzoYqSPAnDFndjZozuwgIg6kyp7W/5g1+Ok6Z66MiJ/vfSIzt9b1febu1GdJfrPVBXPGnFls4Jyp/9N5OvAru1m38bRnwJzpjllj1iw2UNZENWprHbC67pjYHBF/EhGPGqBu42nTFM6Z+nfzXuClVO+lExr+vicBd2TmU4H3UI14IjM/mplHZOYRVJ+jf8vMK3pe99qF5zPzlmF/FZOWmquoVu/ZulDQ0+t7X0S8oGfb8zPzH+rexPsz85LM/Fr9+KtUH8Ifr7d9JfDpzPxC3WP9G8AjLXV4HXBBZl5Q7+siYCPwUz3bfCgz/zUz76MaynbEEOf8WuDMzPxSXbe3UPVGr+3Z5n9m5u318Zr8AlWP7CXAX9XfGjf5FHBET0/pa4G/qI/7MqqhjR/KzK1ZXYP5SeBVQ5zbrnwqM/+5/lt/lPbf4X7A3QPsdx/grkVl36X6B6wff0wVho+nep+cFRHPX7TN3XW9NJ3MGXNmsf0YLGe2iYiV9b7Pzsx/6fNl5wI/CKwG3gi8PSJOWLSNOTPdzBlzZrH9GDxnfgv44G52ONuemQ9mjVmz2H4MljUHAiup/uY/Vu/3WVSXlvdjlts0RwKbMvOazHyQqvPw2EXbHAucXd8/j2oE5OK5rE5g547HkZq0Tp3bgFXRc/1jZv5/mblf/Vxvfa/vfWFEPCciPh8Rt0bEd6l6dReG9j2xd/usrsm8raUOTwJeVYfhnRFxJ9UEXWt6trmp5/69VB0Ju+uJVD3MC3W7p67bQT3bXL/4Rb0y807g/wLPAH5/F9vdTdWTfnxddAJVKEB13s9ZdN6vpbo2d5fqXtmFibH6mlSL/n+Hd9B/hwzLine truncated
"text/plain": [
"<Figure size 1152x360 with 8 Axes>"
]
@@ -366,7 +328,9 @@
"\n",
"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-evalutes the gradient function, and iterates a few times to optimize the initial velocities with a learning rate (step size) `LR`."
"The following cell re-evalutes the gradient function, and iterates a few times to optimize the initial velocities with a learning rate (step size) `LR`.\n",
"\n",
"`field.functional_gradient` has a parameter `get_output` that determines whether the original results of the function (`simulate` here) are returned, or only the gradient. As it's interesting to track how the loss evolves over the course of the iterations, let's redefine the gradient function with `get_output=True`.\n"
]
},
{
@@ -378,22 +342,16 @@
"name": "stdout",
"output_type": "stream",
"text": [
"22.651642\n"
]
},
{
"ename": "RuntimeError",
"evalue": "you can only change requires_grad flags of leaf variables.",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mRuntimeError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-7-d1bfd8f0807f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 20\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 21\u001b[0m \u001b[0msim_grad\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfield\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfunctional_gradient\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msimulate\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mwrt\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mget_output\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[0;32m---> 22\u001b[0;31m \u001b[0;34m[\u001b[0m\u001b[0mloss\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvelocity_grad\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msim_grad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minitial_smoke\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minitial_velocity\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 23\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mloss\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 24\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/Dropbox/mbaDevelSelected/phiflow-v2/phi/field/_field_math.py\u001b[0m in \u001b[0;36mwrapper\u001b[0;34m(*fields)\u001b[0m\n\u001b[1;32m 244\u001b[0m \u001b[0mINPUT_FIELDS\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mextend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfields\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 245\u001b[0m \u001b[0mtensors\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mfield\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvalues\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mfield\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mfields\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 246\u001b[0;31m \u001b[0mresult_tensors\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtuple\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtensor_gradient\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mtensors\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[0m\u001b[1;32m 247\u001b[0m \u001b[0mproto_fields\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 248\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mget_output\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/Dropbox/mbaDevelSelected/phiflow-v2/phi/math/_functions.py\u001b[0m in \u001b[0;36mwrapper\u001b[0;34m(*values)\u001b[0m\n\u001b[1;32m 1232\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0m_\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn\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 1233\u001b[0m \u001b[0mARG_INDICES\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg_index\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1234\u001b[0;31m \u001b[0mresults_native\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgrad_native\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnatives\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[0m\u001b[1;32m 1235\u001b[0m \u001b[0mproto_tensors\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 1236\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mget_output\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/Dropbox/mbaDevelSelected/phiflow-v2/phi/math/_functions.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 1216\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgradf\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 1217\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgradf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbackend\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfunctional_gradient\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnative_function\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mshifted_wrt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mget_output\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mget_output\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1218\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgradf\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[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1219\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1220\u001b[0m \u001b[0mgrad_native\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mGradientFunction\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~/Dropbox/mbaDevelSelected/phiflow-v2/phi/torch/_torch_backend.py\u001b[0m in \u001b[0;36meval_grad\u001b[0;34m(*args)\u001b[0m\n\u001b[1;32m 548\u001b[0m \u001b[0mwrt_args\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0marg\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0marg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0menumerate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mwrt\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 549\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0marg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mwrt_args\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 550\u001b[0;31m \u001b[0marg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrequires_grad\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 551\u001b[0m \u001b[0moutput\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mf\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[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 552\u001b[0m \u001b[0mloss\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maux\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0moutput\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moutput\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moutput\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mtuple\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0moutput\u001b[0m\u001b[0;34m,\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[0;31mRuntimeError\u001b[0m: you can only change requires_grad flags of leaf variables."
"Optimization step 0, loss: 22.651701\n",
"Optimization step 1, loss: 22.533409\n",
"Optimization step 2, loss: 22.277126\n",
"Optimization step 3, loss: 22.404198\n",
"Optimization step 4, loss: 22.268917\n",
"Optimization step 5, loss: 22.391558\n",
"Optimization step 6, loss: 22.448061\n",
"Optimization step 7, loss: 22.387615\n",
"Optimization step 8, loss: 24.646418\n",
"Optimization step 9, loss: 26.083023\n"
]
}
],
@@ -401,26 +359,16 @@
"#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",
"sim_grad_wloss = field.functional_gradient(simulate, wrt=[1], get_output=True)\n",
"\n",
"for optim_step in range(10):\n",
" initial_smoke = DOMAIN.grid(math.zeros(inflow_loc=4))\n",
" #initial_velocity = DOMAIN.staggered_grid(0) * math.ones(inflow_loc=4)\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",
"\n",
" \n",
" # sim_grad__wloss returns: [loss,s,v,velocity _grad] \n",
" sim_grad = field.functional_gradient(simulate, wrt=[1], get_output=True)\n",
" [loss,s,v, velocity_grad] = sim_grad(initial_smoke, initial_velocity)\n",
" print(loss)\n",
" [loss,s,v, velocity_grad] = sim_grad_wloss(initial_smoke, initial_velocity)\n",
"\n",
" print('Optimization step %d, loss: %f' % (optim_step,loss))\n",
"\n",
" initial_velocity = initial_velocity - LR * velocity_grad\n"
]
@@ -431,18 +379,31 @@
"source": [
"The loss should have gone down significantly, from above 60 to below 40, and now we can also visualize the initial velocities that were obtained in the optimization.\n",
"\n",
"**TODO, test ini-vel images below**"
"The following images show the resulting three initial velocities (we're skipping the fourth set with `inflow_loc[0]` because it only contains zeros).\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 8,
"metadata": {},
"outputs": [],
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1cAAAFgCAYAAABXIoPMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABJ1klEQVR4nO3de7xcZX33/e+XEAgnTSAxRAKGFqzVVmPdDdoWy8NB0ac19Cki6oPhbijS3lr72ANRvMWi3g21Fm21thGQYFVAqBItihDhQVtFglLkICVSkIQA5oByChDyu/9Ya8+e2Xut7Fl71mFm9uf9es1rz7rmmjXXnL57XWut6xpHhAAAAAAAvdmt6QYAAAAAwDCgcwUAAAAAJaBzBQAAAAAloHMFAAAAACWgcwUAAAAAJaBzBQAAAAAl2L3pBgDT0fF2bO5xHTdLV0fE8aU0CMBQImsAVI2c6UTnCmjAZknrelyHpblltAXA8CJrAFStjpyxfbykj0uaIen8iFg57vY9JV0s6RWStkh6U0Tca3uRpDsl3ZVW/W5EnJHe5xWSLpK0l6SrJL0rSvgBYDpXQFN26/Gs3J07y2kHgOFG1gCoWoU5Y3uGpE9KOk7SBkk32V4TEXe0VVsuaVtEHGb7ZEnnSnpTetuPI2Jxxqo/JekPJd2opHN1vKSv9fZEGHMFNGe33Xq7AEA3yBoAVas2Z5ZIWh8R90TE05IukbR0XJ2lklan1y+XdIxt563Q9gJJz4mI76ZHqy6WdMIUnvkEpCYAAACAfnWQpPvbljekZZl1ImKHpJ9JOiC97VDbP7D9/9s+sq3+hknWOSWcFgg0wWaPMIDqkTUAqlZOzsy13T50a1VErOp1pZI2STokIrakY6y+bPslJaw3F50roCls8ACoQ8VZ08VA81dL+pikl0o6OSIur7RBAOrXe85sjoiRnNs2Sjq4bXlhWpZVZ4Pt3SU9V9KW9JS/pyQpIm62/WNJL0zrL5xknVPC1h3QhNG9PIyDAFClirOmbaD56yS9WNKbbb94XLWfSDpV0ucreIYAmlb9Ns1Nkg63fajtPSSdLGnNuDprJC1Lr58o6ZsREbbnpTkl278g6XBJ90TEJkk/t/3KdGzW2yRdWcbLwZErAAAwVa2B5pJke3SgeWsWr4i4N72NaQcBFBYRO2y/Q9LVSo6QXxgRt9s+R9K6iFgj6QJJn7W9XtJWJR0wSXq1pHNsPyNpp6QzImJretsfa2wq9q+phJkCJTpXQHM4+gSgDtVmTdZA8yOqfEAAfajibZqIuErJdOntZe9vu75d0hsz7neFpCty1rlO0q+U21I6V0Bz6FwBqEP/DjQHMCzYpmmhcwU0gRm8ANShnKzpdaA5gGHGNk0HOldAUwgiAHWoNmtaA82VdKpOlvSWKh8QQB9im6aFVwIAAExJ+mOdowPN75R02ehAc9tvkCTbv257g5LxEP9s+/bmWgwA1eLIFdAEDqEDqEMNWdPFQPOb1Pl7MgCGCds0HehcAU0hiADUgawBUDVypoXOFdAUgghAHcgaAFUjZ1p4JQAAAACgBBy5AprA+ckA6kDWAKgaOdOBzhXQFIIIQB3IGgBVI2da6FwBTWAvD4A6kDUAqkbOdOCVAAAAAIAScOQKaAp7eQDUgawBUDVypoXOFdAUgghAHcgaAFUjZ1roXAFN4PxkAHUgawBUjZzpwCsBAAAAACXgyBXQFPbyAKgDWQOgauRMC68E0ITRQ+i9XLp6GB9v+y7b622vyLh9T9uXprffaHtRWr7I9pO2b0kv/1TuCwCgFjVlDYBpjJzpwJEroCkVh4ntGZI+Kek4SRsk3WR7TUTc0VZtuaRtEXGY7ZMlnSvpTeltP46IxZU2EkD1hmzDBUAfImdaeCWA4bVE0vqIuCcinpZ0iaSl4+oslbQ6vX65pGNsu8Y2AgAADA2OXAFNqX4vz0GS7m9b3iDpiLw6EbHD9s8kHZDedqjtH0j6uaT3RcS3Km4vgCqwRxlA1ciZFjpXQBPKmbZ0ru11bcurImJVrytNbZJ0SERssf0KSV+2/ZKI+HlJ6wdQB6ZIBlA1cqYDnSugKb0H0eaIGNnF7RslHdy2vDAty6qzwfbukp4raUtEhKSnJCkibrb9Y0kvlLROAAYLGz0AqkbOtNC5AppQz16emyQdbvtQJZ2okyW9ZVydNZKWSfqOpBMlfTMiwvY8SVsj4lnbvyDpcEn3VN1gACVjjzKAqpEzHehcAUMqHUP1DklXS5oh6cKIuN32OZLWRcQaSRdI+qzt9ZK2KumASdKrJZ1j+xlJOyWdERFb638WAAAAg4POFdCUGvbyRMRVkq4aV/b+tuvbJb0x435XSLqi8gYCqB57lAFUjZxpoXMFNIUgAlAHsgZA1ciZFjpXQBM4PxlAHcgaAFUjZzrwSgAAAABACThyBTSFvTwA6kDWAKgaOdMyVK+E7a/ZXraL2//J9v/qcl3X2z6tvNa11hu2D+txHbfbPqqcFhV63Pm2b7D9qO2PlrTOebZ/ZHuvMtbXQzveafvcGh8wCaJeLmgMWVOtIc+a37V9aY0PSNYMKHKmWkOeM2zTNKjvn43te20f203diHhdRKxO73eq7W+Pu/2MiPhgFe2sU0S8JCKulyTbH7D9LzU99OmSNkt6TkT8WUnrXCHpooh4UpJsn2T7P2w/Yfv68ZVtL7Z9c3r7zbYXd/tA6WfpSduPpZdvtN38aUlvtf28Hp9P9wiivkLWTDRds8b2C21fafuntrfavtr2L3X7QOmG7Pa2rLlr9LaI+Iqkl9h+aUnPa3JkTd8gZyaarjnTzvbb0o5q1x1gtmn613A9G1TtBZLuiIgoekfbE05Btb2nkh+wbQ/SrZI+JmllRv09JF2Z1p8jabWkK9Pybv1uROybXl4zWphOSf41SW8rsC4A1Wg0ayTNVvID278kab6k7ynJniLe0ZY14ztmX1CyYQegOU3nzOj95kh6r6Tbi7ZDbNP0pYHqXI3uubH9t7a32f5v269ru/1626fZ/mVJ/yTpVWlv/pH09otsfyi9Psf2V9M9k9vS6wu7aMPz0z0F+7eVvdz2Ztsz0+U/sH1nut6rbb8gZ13PtX1x2ob7bL/P9m5tt/9hup5Hbd9h+9fS8nttH2v7eCVfyDelz/M/bb/R9s3jHufdtidsGNje3/YG27+bLu9re73tCV9G2xcpCY2/TB/rWNt72v6Y7QfSy8fScJHto9J1n2n7QUmfyXgJjpD0SERsGC2IiGsj4jJJD2TUP0rJOMGPRcRTEfH3kizp6KzXdwqul/R/l7SuybGXp2+RNdM7ayLiexFxQURsjYhnJJ0n6ZdsH5D1+k7B9SJrpj1yZnrnTJu/lvT3So6ilel6DVHO2D7e9l3pe7oi4/Y9bV+a3n6j7UVp+XFOznT6Yfr36Lb7XJ+u85b0UsqRvkFMzSMk3SVprqS/kXSBbbdXiIg7JZ0h6Ttpb352xnp2U/LleIGkQyQ9KekTkz14RDwg6TuSfr+t+C2SLo+IZ2wvVRIO/4+keZK+pWQvZZZ/kPRcSb8g6beV7GH4H5Jk+42SPpCWPUfSGyRtGdeWr0v635IuTZ/ny5TsbT00DeNRp0i6OOO5bJX0B5I+nX6gzpN0S0Rk1T1V0uck/U36WNdKOkvSKyUtlvQySUskva/tbgdK2l/Ja5y1l/ZXlbyX3XqJpFvH7WW6NS3v1ufS4P+G7ZeNu+1OJc+jepyfPAjImrG2TLesGe/Vkh6MiC2T1hzz1+kG6r974niSOyUtsv2cHtrUHbKm35EzY22Zdjlje4mkESWd56mYFts0tmdI+qSk10l6saQ3237xuGrLJW2LiMOUvPejY842KznC96tKOtSfHXe/t0bE4vTycG8vRGIQU/O+iPh0RDyr5LSwBUpO2ygkIrZExBUR8UREPCrpw0rCoBufl/RmSUpD8OS0TEoC8K8j4s6I2KEkKBaP39OTflBOlvSeiHg0Iu6V9FEloSFJpyn50t8UifURcV8Xz+spSZdK+n/Tx3mJpEWSvppT/xuSvihpraTXS3p7V69A4q2SzomIhyPip5L+qq39krRT0tnpUaYnM+4/W9KjBR5vX0k/G1f2M0n7FWjvIiXBeJ2kq23Pbrv9USX/GOrBBk+/I2t2/byGOWta0r3/n5T07gJ3O1PJBuZBklZJ+ortX2y7fbQts6fSpsLImn5Gzuz6eQ1tzqSv2T8qOYV4Z4F2trd3kabHNs0SSesj4p6IeFrSJZKWjquzVMl3SJIul3SMbUfED9KdCFJy6uVeo0ckqzKIqfng6JWIeCK9um/Rldje2/Y/p4eufy7pBkmz0w/7ZK5Qcnh+gZI9mjuV7M2Rkg/5x20/4uTQ/VYlp64dNG4dcyXNlNQeLve11TtY0o+LPq/UaklvSUPyFEmXpQGVZ5WkX1EyCLPIntnna2L7n9+2/NNIzvvNs03dd4wk6TEle7zaPUddhllE/HtEPJn+8/lrSY9IOrKtyn6a2HnD9EXWTG5Ys0ZSMvOXpG9I+seIyNtbP0FE3JhuYD4VyYQE/65kQ2/UaFseKdomDB1yZnLDmjN/rORsnO8WuE/LNNumOUjS/W3LGzTxM9iqk+4I+Jmk8ady/76k74/7/HwmPSXwf40/ajxVg9i56tZkAxT/TMlg5SMi4jlKAkVKQmPXK47YpuQf7puUHD6/JKJ1qtr9kt4eEbPbLntFxH+MW81mSc8oCa5Rh0ja2LaeX9TkJjzP9Iv6tJIv2Vs08RBoSxq8q5QcYv9jF5tS9QFNbH/7ecWTvQe3Snphgce7XdJLx334X6qpDQKVkva1r+uXJf3nFNdVDKfqDBOyZviyZnSQ+TckrYmIDxe5b4asrLk3In7e43onR9YMC3Jm+HLmGEm/Z/tBJ+O4fkPSR21PejpnjkHfpplre13bpdRJf9Kjnueq82jmWyM5XfDI9HJK1n2LGubUfEjSQufPJLefknOSH3EykPPsguvLine truncated
"text/plain": [
"<Figure size 864x360 with 6 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(initial_velocity.inflow_loc[i].vector[0].numpy('y,x'), origin='lower', cmap='bwr')\n",
"fig, axes = pylab.subplots(1, 3, figsize=(12, 5))\n",
"for i in range(INFLOW.shape.inflow_loc-1):\n",
" im = axes[i].imshow(initial_velocity.inflow_loc[i].vector[0].values.numpy('y,x'), origin='lower', cmap='bwr')\n",
" pylab.colorbar(im, ax=axes[i])\n",
" axes[i].set_title(f\"Initial velocity x for {math.to_int(INFLOW_LOCATION).inflow_loc[i].vector.unstack_spatial('x,y')}\")\n",
"pylab.tight_layout()"
@@ -450,13 +411,26 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 9,
"metadata": {},
"outputs": [],
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1YAAAFgCAYAAAC44OjyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABLaElEQVR4nO3de7xcdX3v/9eb3LgTQjDEBAyFWKoexHYX2p+15aeA2FvoOYqoR+M5cCht6Wl/9mKsVixqi22tttXWkwoSrRUo1hItNXL9eakiQblfSqRQAuESQuQaQsLn/LHWhsnea+09a8+stWbm+34+HvuRme98Z813Zc+89/quy2cUEZiZmZmZmdnM7db2AMzMzMzMzIadJ1ZmZmZmZmY98sTKzMzMzMysR55YmZmZmZmZ9cgTKzMzMzMzsx55YmVmZmZmZtaj2W0PwCxlJ0qxucdlXAfrIuLEvgzIzEaaM8fMmpJi3nhiZdaizcD6HpchWNiPsZjZ6HPmmFlTUswbT6zM2rZbj2fkPvdcf8ZhZmlw5phZUxLLG0+szNqWWOiYWcucOWbWlMTyxhMrszZJvYeOmVm3nDlm1pQE8yattTUzMzMzM6uBj1iZtS2xvTlm1jJnjpk1JbG88cTKrE0JHiY3sxY5c8ysKQnmjSdWZm1LLHTMrGXOHDNrSmJ5k9bampmZmZmZ1cBHrMzaltjeHDNrmTPHzJqSWN54YmXWpgTPPzazFjlzzKwpCeZNWmtrNoh22623HzOzKmrOHEknSrpD0gZJqwoenyfpwvzxayQty9uXSXpa0vX5z6f6v/Jm1qjEtnF8xMrMzMz6QtIs4JPA8cBG4FpJayPi1o5upwKPRsThkk4BPgK8OX/sBxFxVJNjNjPrF0+szNqU4GFyM2tR/ZlzNLAhIu7KXk4XACuAzonVCuAD+e2LgU9IUp2DMrMWJLiN44mVWdsSCx0za1m9mbMEuLfj/kbgmLI+EbFD0g+BA/LHDpX0feAx4H0R8Y06B2tmNUtsG8cTK7O2JRY6Ztay3jNnoaT1HfdXR8TqXhcKbAIOiYhHJP0E8M+SXh4Rj/Vh2WbWhsS2cTyxMjMzsyo2R8RYyWP3AQd33F+atxX12ShpNrAf8EhEBPAMQERcJ+kHwEuB9ZiZDQFPrMzalOD5x2bWovoz51pguaRDySZQpwBvndBnLbAS+DbwRuDKiAhJBwJbImKnpB8BlgN31TlYM6tRgts4nliZtS2x0DGzltWYOfk1U2cC64BZwHkRcYuks4H1EbEWOBf4nKQNwBayyRfAzwJnS3oWeA44IyK21DZYM6tfYts4nliZtSnBvTlm1qIGMiciLgUundD2/o7b24A3FTzvi8AXax2cmTUnwW2ctNbWzMzMzMysBj5iZda2xPbmmFnLnDlm1pTE8iattTUbRLvt1ttPFySdKOkOSRskrSp4fJ6kC/PHr5G0rN+raWYDooHMMTMDkssbH7Eya1MD5x9LmgV8Ejie7Ms6r5W0NiJu7eh2KvBoRBwu6RTgI8Cbax2YmTUvwWsezKwlCeaNJ1Zmbas/dI4GNkTEXQCSLgBWAJ0TqxXAB/LbFwOfkKT8e2XMbJQktqFjZi1KLG/SWluzNC0B7u24vzFvK+wTETuAHwIHNDI6MzMzsxHgI1ZmberPYfKFktZ33F8dEat7XaiZjaAET80xs5YkmDeeWJm1rffQ2RwRY1M8fh9wcMf9pXlbUZ+NkmYD+wGP9DowMxtAiW3omFmLEssbT6zM2lZ/6FwLLJd0KNkE6hTgrRP6rAVWAt8G3ghc6eurzEZUYhs6ZtaixPImrbU1S1B+zdSZwDrgNuCiiLhF0tmSfjnvdi5wgKQNwLuASSXZzczMzAbNTL9SRtLRkq7Pf26Q9Cu9jsVHrMza1ND5xxFxKXDphLb3d9zeBryp9oGYWbsSvObBzFoy+F8pczMwFhE7JC0GbpD05XyH9Ix4YmXWNm/kmFmTnDlm1pTB/kqZpzr67A70fAnEyKWrpH+VtHKKxz8l6Q+7XNbVkk7r3+ieX25IOrzHZdwi6dj+jKjS6y6S9HVJj0v6aJ+WeaCk2yXt0Y/l9TCO35T0kYZfNLlvJR81zpx6jXjm/JKkCxt+UWfOEHPe1GvE82ZUt3F6+koZScdIugW4CTijl6NVMCQTK0l3Szqum74R8YaIWJM/752Svjnh8TMi4oN1jLNJEfHyiLgaQNIHJP19Qy99OrAZ2DcifqdPy1wFnB8RTwNIOlnSv0l6StLVEztLOkrSdfnj10k6qtsXyt9LT0t6Iv/5WsfDfwe8TdKLelwfG3LOnMlSzRxJL5V0iaSHJW2RtE7Sj3b7QvnG67aOzLlj/LGI+DLwcklH9mm9bAg5byZLNW86SXpHPkntevI7ots4CyWt7/g5vZ8Lj4hrIuLlwE8C75G0ey/LG4qJlQ2UlwC3zqRinLIy3hPb5pFVo+sMzS3Ax4FzCvrPBS7J++8PrAEuydu79UsRsXf+c8J4Y36d0b8C76iwrN5577HZVFrNHGA+WdXMHwUWAd8ly6AqzuzInImTsi+Qbcw1x5ljVqbtvBl/3v7AHwC3VB0Ho7eNszkixjp+Jn5PZ5WvlBn/PU36SpmIuA14AnhFT6vby5PbML6HRtKfS3pU0n9IekPH41dLOk3SjwGfAn46n7VvzR8/X9KH8tv7S/pKvify0fz20i7G8OJ8j8CCjrZXSdosaU5+/39Kui1f7jpJLylZ1n6SPpuP4R5J75O0W8fj/ytfzuOSbpX043n73ZKOk3Qi2Yfvzfl63iDpTZKum/A675I0aWOgYt/zyQLi9/PXOk5ZpZWPS7o///l4HiRIOlbSRknvlvQA8JmC/4JjgK0RsXG8ISIuj4iLgPsL+h9Ldm3gxyPimYj4K0DAawv6zsTVwC/0aVnd8UbOQHPmpJ05EfHdiDg3IrZExLPAx4AflXRA0f/vDFyNM8dyzpu086bDnwB/RXb0rJ+uZvTy5vmvlFG2k/0Usp1hndaS/W6h4ytl8ufMBsjfw0cAd/e0ur08uUXHAHcAC4E/Bc6VpM4O+czzDODb+ax9fsFydiP7ILwEOAR4GvjEdC8eEfeTfd/Pf+tofitwcUQ8K2kFWRD8V+BA4BtkeyWL/DXZzPlHgJ8j25PwPyALBLKL7d4B7Av8MpNn2F8F/hi4MF/PV5K9gQ7Ng3fc24HPFrx+130j4p3A54E/zV/rcuC9wE8BRwGvJLuI8H0dTzsIWED2f1y0V/a/kP0uu/Vy4MYJe5NuzNu79fk85L8m6ZUTHruNbD2a4esdhoUz54WxpJY5E/0s8EBEVPkC7T/JN0q/pcnXjdwGLJO0bw9j6p4zZxg4b14YS3J5I+loYIxs4jwTSW3j9PiVMj9DVgnweuBLwK9HRE+T2WFNyHsi4u8iYifZqWCLyU7RqCQiHomIL0bEUxHxOPBhsg9+N/4BeAtAHnin5G2Qhd2fRMRt+S/8j4GjJu7RUVYi8hTgPRHxeETcDXyU7EMPcBrZB/zayGyIiHu6WK9ngAuB/56/zsuBZcBXeulb4m3A2RHxUEQ8DPxRx/gBngPOyo8uPV3w/PnA412+FsDeZBcddvohsE+F8S4jC8GrgHWS5nc8/jjZHwGzTs6cqddrlDPnefne/k+S/WHu1rvJNiqXAKuBL0s6rOPx8bHMn8mYbCQ5b6Zer5HNm/z/7G/ITh9+rtvnTRjvMhLbxomISyPipRFxWER8OG97f0SszW9vi4g3RcThEXF05BUEI+JzkV3Pd1RE/HhE/HOvYxnWidUD4zfihVKJe1ddiKQ9Jf2f/PD0Y8DXgfn5G3s6XyQ7BL+YbA/mc2R7bSB7Q/+lpK3KDs9vITtdbWKVkoXAHKAzSO7p6Hcw8IOq65VbA7w1D8S3k83gn+lD34lezOTxv7jj/sORnddb5lG6nxRBdv7rxD27+9JlcEXEtyLi6fwPzZ8AW4HXdHTZh8kTt3p57/EwcOZMb1QzB8gqewFfA/4mIsr2zk8S2YXRj+cbXmuAbwE/39FlfCxbq45pxpw5g855M71RzZtfJzsr5zsVnvM8b+O0b/hGXM10Fx/+DtkFycdExL5k4QFZQEy94IhHyf7IvpnsEPkFEc+fnnYv8KsRMb/jZ4+I+LcJi9kMPEsWUuMO4YWL7u4FDmN6k9Yz/1BuJ/tAvRX43BTr0nXfAvczefyd5w1P9zu4EXhphde7BThywmkRRzKzCzwhG1/nsn4MuGGGy6rOp+WMGmfO6GXO+IXkXwPWju8N7UFR5twdEY/1uNzuOHNGifNm9PLmdcCvSHogv27r/wE+KmnaUzhLeBunYcM34moeBJaqvGLcPmTnHG9VdpHmWRWX/w9k5wa/kRcOkUN2Xux78kPO4xdvvmnik/PD/BcBH5a0T34Y/V28UD3m08DvSvoJZQ6feKi9Yz2XqeOC0Nxnyc6nfjYivjn5aTPu2+kLwPuUfU/DQuD97Fr9ZjrfJduD9vyeLkmzlJW7nA3sJml35RfMkl14uRP438ouKj0zb78yf+47Jd1d9EKSDpH0aklz82X+HtketW91dPs5sqo5zUksdEacM2fEMkfZtU/rgG9FxKqJC1N2AXvhxpWk+ZJeny9vtqS3kW3cfrWjmzPHZsp5M2J5A7yTbPJzVP6znuz0w/fmz/U2zoAbvhFXcyXZkYwHJBVdjPZxYA+yvSrfYdc/dt1YCywnu5D5+T0AEfEl4CPABfnh95uBNxQvgt8EngTuAr5JFl7n5cv5R7Jzov+B7FS3fya7SHKif8z/fUTS9zraP0dWNrKbEKjSt9OHyD74N5J9udr38rauRMR24Hzy859zbyf7Y/C3ZHuYnib7/oXx/ieRhf1W4H8CJ+XtkJ1a0BkinfbJl/ko2R6zE4E3RH4Reh50P0922kBzEgudEefMGbHMAX6F7PtN/ode+G6YJyQdkj9+MDBxT/24OfnYHib7nf8mWV79e0eftwD/p9vLine truncated
"text/plain": [
"<Figure size 864x360 with 6 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(initial_velocity.inflow_loc[i].vector[1].numpy('y,x'), origin='lower', cmap='bwr')\n",
"fig, axes = pylab.subplots(1, 3, figsize=(12, 5))\n",
"for i in range(INFLOW.shape.inflow_loc-1):\n",
" im = axes[i].imshow(initial_velocity.inflow_loc[i].vector[1].values.numpy('y,x'), origin='lower', cmap='bwr')\n",
" pylab.colorbar(im, ax=axes[i])\n",
" axes[i].set_title(f\"Initial velocity y for {math.to_int(INFLOW_LOCATION).inflow_loc[i].vector.unstack_spatial('x,y')}\")\n",
"pylab.tight_layout()"
@@ -473,15 +447,26 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 10,
"metadata": {},
"outputs": [],
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA54AAAESCAYAAACVYgptAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAseElEQVR4nO3de3ycdZn38e810zRNk7bpmUI5FYqAHMpjBVRU5KAsIodd8LDqi3V1q6vu6uOuj+izj4qPuh7xsPqoCGxReUROLh520QooAitQpOVUgba00DM9JE3SNGlmfvtHhtdWvH5J7mR+kzuTz/v16qvpd2bu+3dP5prJ1Zlct4UQBAAAAABAKoXRXgAAAAAAoL7ReAIAAAAAkqLxBAAAAAAkReMJAAAAAEiKxhMAAAAAkBSNJwAAAAAgqQm13JmZce4WQFIIwUZ7DR5qFOhHjQL5Ro0C+ebVaE0bz37F2u8SyJXSaC9gENQoxjtqFMg3ahTIN79G+agtAAAAACApGk8AAAAAQFI0ngAAAACApGg8AQAAAABJ0XgCAAAAAJKi8QQAAAAAJEXjCQAAAABIisYTAAAAAJAUjScAAAAAICkaTwAAAABAUjSeAAAAAICkaDwBAAAAAEnReAIAAAAAkqLxBAAAAAAkReMJAAAAAEiKxhMAAAAAkBSNJwAAAAAgKRpPAAAAAEBSNJ4AAAAAgKRoPAEAAAAASdF4AgAAAACSGrTxNLNJZna/ma00s8fM7PJKvtTMnjazFZU/i5KvFsCfoEaBfKNGgXyjRoHamDCE6/RIOiOE0GlmDZLuNrP/qFz24RDCTemWB2AIqFEg36hRIN+oUaAGBm08QwhBUmflnw2VPyHlogAMHTUK5Bs1CuQbNQrUxpB+x9PMima2QtI2SctCCPdVLvqMmT1sZl8xs8ZUiwQwMGoUyDdqFMg3ahRIb0iNZwihFEJYJGm+pJPN7DhJH5V0tKSXSpoh6SPebc1siZktN7Pl1VkygBeiRoF8o0aBfKNGgfSs/9MFGW5g9nFJe0IIX9ovO13SP4YQzhvktkEqDmOZQD0pKYRgqbZOjQIjRY0C+UaNAvnm1+hQptrONrPWytdNks6W9Aczm1fJTNKFkh6t5nIBDA01CuQbNQrkGzUK1MZQptrOk3StmRXV36jeEEL4mZndYWazJZmkFZLek26ZAAZAjQL5Ro0C+UaNAjWQ+aO2I9oZHz8AlPojQiNBjQISNQrkHTUK5NswP2oLAAAAAMBI0HgCAAAAAJKi8QQAAAAAJEXjCQAAAABIisYTAAAAAJAUjScAAAAAICkaTwAAAABAUjSeAAAAAICkaDwBAAAAAEnReAIAAAAAkqLxBAAAAAAkNWG0F1D/LJLHev5yJA9VWAsAAAAA1B7veAIAAAAAkqLxBAAAAAAkReMJAAAAAEiKxhMAAAAAkBSNJwAAAAAgKabaVskB005z83fNOsPNj57S5+br9vjfkv/YutPN723/tpsH+dsHAAD7i02fZ5o8kAdmE928WGzJvK1SqdPNQ+jNvC1kxzueAAAAAICkaDwBAAAAAEnReAIAAAAAkqLxBAAAAAAkNWjjaWaTzOx+M1tpZo+Z2eWV/HAzu8/MVpvZjyz2m78AkqJGgXyjRoF8o0aB2rAQBp7aZmYmqTmE0GlmDZLulvQBSR+SdEsI4Xoz+7aklSGEbw2yrSAVq7T0tCZNPNDNP7XgXW7+j5/tcvPS+X9WlfUU77/fzdd9eq2bv+ueKW5+5+5vRvfBRK9aKSmEEBujmNl4rVEgHWp0bPC/RdNbjnfzFxdPd/OZhcluvr3sT798vHRXdEW7Oh+OXoZqokbHslj/fuK0t7n5JXPmu/nJM/Zk3vf9O/16v3HbBjdf2f4DN+dn5sH4NTroO56h3/PPvg2VP0HSGZJuquTXSrqwOgsFkAU1CuQbNQrkGzUK1MaQfsfTzIpmtkLSNknLJK2R1BZCeP5kkRskHZRkhQAGRY0C+UaNAvlGjQLpDanxDCGUQgiLJM2XdLKko4e6AzNbYmbLzWz58JYIYDDUKJBv1CiQb9QokF6mqbYhhDZJd0p6maRWM5tQuWi+pI2R21wZQlgcQlg8koUCGBw1CuQbNQrkGzUKpDOUqbazzay18nWTpLMlrVJ/UV5cudqlkm5NtEYAA6BGgXyjRoF8o0aB2hjKVNsT1P8L1UX1N6o3hBA+ZWYLJF0vaYakhyS9LYTQM8i2cjfpyzTBzZe/+m/d/MRbX+NvZ8Om6qzn0dVuXrr4PDePTbsNjY1u/qGL2qP7/tozn4lcUoreBsNR9Wl8dV2jo6Wx4QA3XzT5IjdvUZObb7Pn3HxNtz8Zc0/P+iGsDmlRo6PDP67pLS9287874EI3f9cx/nTKA1/rf0sLC+a4eXntNjff/pt9bi5JP3jsUDf/xpYH3Xx926/cnImZg6FG8yQ2pfYl0y518yuO9c+88LILd7h58ZXHuHlo9ifUDsS6/Em4pd+ucvP//LeZbv6hxzvc/MH2a918/NW0X6N+17WfEMLDkk5y8rXq/ww8gFFEjQL5Ro0C+UaNArWR6Xc8AQAAAADIisYTAAAAAJAUjScAAAAAICkaTwAAAABAUoNOta3qzkZp0pfZpOhlXzv2w27+tytf7ebFe//TzcMz/vS7rju3u/nGZ6ZF15TFi97pH1vpmV1uXpjpT96UpI992p/s94V1X3DzEPYOsjr4qjuNr5rqeRpfodDi5mdOebebX368P8158Vn+lNpyd9nN9271v9WbN0518x8+7dehJF21/df+ttrudfOgvui2MBBqdDRcMP0yN//Bxf6k5+a3RiZdzmit1pIys63+88Oem59y8zddf4ib/7zty5E9MGW+HzU6GmJngjh1mv86esdf+1NqJ77j5VVb02jp/Vf/dfeMa/wpuL9r/46b1+/rtF+jvOMJAAAAAEiKxhMAAAAAkBSNJwAAAAAgKRpPAAAAAEBSNJ4AAAAAgKTGxVTb6S0nRC/b8W8XuXnXj55081sfONzNG4v+/Vg0P2+I5DMae928c58/SSymoeBP2Hz5qzZFbzPxpfPcfOHfrHHzNbtuzbQmPI9pfNXh34UHTHuFm3/p8Fe5+ZsuftbNi0fPzrSafQ9u9rczNz5V279B/P8Dd/3Wf3646BfT3fzu9n+JbKl2z/tjEzVaDYXCZDd/w7S/d/ObPuK/PhXOXZxtxxNG7/6x9o5M1y/9dpWbv+0LB7n5DTuucPPxN2WeGk3Lv2tfMe19bn7H3+908wmXnFqV1YSZ/mtcNdkO/2wQWfXd+Ds3P+PrM9z8nvZvRrY01l+nmWoLAAAAABgFNJ4AAAAAgKRoPAEAAAAASdF4AgAAAACSovEEAAAAACSVbVRq7vlTuM5qPC96i95l/vTa79x9pJs3R6bXSv4U2diaXjpnR3RNnql9/oS0TV3Nbt7S0OfmbWsaovuYc/YUN3/XnBPd/GO7fu7mQf6+geHxa2j2VH/S5S3H+/kpb33O3/psf3rtvoe2uHmpy691i8xXLK/vdvNCi///fg0nHuBvSNLMN/r5HcdsdfPXf/VDbr6s/auRPZSi+wayOnraBW5+w3ufcfPCJWdk2r49udbNw6H+RFjblu11N8yZmen6w1F85TFu/r1dK938jm8e7+bP7X6gamsCikX/58Fvv8R//Ws4b5GbZ53Lal17/AtW+D+rV5MtnO/modmfzh0Tuy++fc86N1/0G/++LpV2Z9rvWME7ngAAAACApGg8AQAAAABJ0XgCAAAAAJKi8QQAAAAAJDVo42lmB5vZnWb2uJk9ZmYfqOSfNLONZrai8ufc9MsF8ELUKJBv1CiQb9QoUBsWwsAzp8xsnqR5IYTfm9kUSQ9KulDSGyV1hhC+NOSdmQXJn86a0o0v+Wj0stcsfNbNl/zyYDdvbvB79QVT/ON67QEdbn7Ugdvd/OFn57r5xII/bbK37O+3MzIFd4LFv9+nvdi/LybO8K+/8Optbr6h7fboPiBJJYUQIjNQs6uHGo1NrpWk+a3+pMufneRPgzzhY9PdPGza6ea7b9/l5j17/KHf5ZK/1mnzetzcMs4OL06L/3/gQBNvPW0/9qfdHnjjnW6+t3dTpu3XL2o0C7NJbn7dIn+q8ptuWuBvZ5v/uvjMp59y8/Zuf78LFvrTaye/7+X+ftf7j/u9y9a5+b6O+EOjcZafT3ytPylfLf5k+pgPv9mf+vmVZ/2HULkcmRI65lGj1WCRk1ucNe0Dbv6L2/wp8Cr4r1v7frbCzVfd7f9gubtnopu//GL/9VtvPtvPB3L9Mje+9yZ/TVMbe938mNP8NcWm2qrsTwR+3Tn+xP1ftX/NzcfOmSP8Gh30R6IQwmZJmytfd5jZKkn+rHIANUeNAvlGjQL5Ro0CtZHpdzzN7DBJJ0m6rxK938weNrNrzMx/mwFAzVCjQL5Ro0C+UaNAOkNuPM2sRdLNkj4YQtgt6VuSjpC0SP3/S/TlyO2WmNlyM1s+8uUCiKFGgXyjRoF8o0aBtIbUeJpZg/oL8boQwi2SFELYGkIohRDKkr4r6WTvtiGEK0MIi0MIi6u1aAB/jBoF8o0aBfKNGgXSG8pUW5N0taRVIYQr9svn7Xe1iyQ9Wv3lARgMNQrkGzUK5Bs1CtTGUOYtvkLS2yU9YmYrKtnHJL3FzBZJCpLWSXp3gvVlEpvONbfRnzYpSdeuPNzN15U3u3nr3ilu/rbD/H0vnOdP6WuY7E+pbZ3or/WXW/xfKyhFhtTG8gXN/n4l6acP+ffFX37X3/eJV/uTuDbojsgeBp6gjGEbMzUaUyz6dSVJXz/ypW5+wuf8CXTh0bVu/rt/bXLzvaUD3byp6E+OO/pQ/3Ef+2+8nrbItL9uf+rhpK4BJtaVt7hxw0n+tNvWi/wp2W+47a/c/KYdX3DzsTNFL7fGfI0ONHl6Yet5bh6bXqvOTje++f3tbn7Vav81aEPBf53+Yuchbn7O0xvcPHRLine truncated
"text/plain": [
"<Figure size 1152x360 with 4 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# re-run sim\n",
"smoke = initial_smoke # DOMAIN.grid(math.zeros(inflow_loc=4)) + INFLOW\n",
"smoke = initial_smoke \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",