replaced phiflow >> op calls with .at()

This commit is contained in:
NT committed 2021-09-16 18:15:49 +02:00
1 parent bc8b2ddf73
commit 0bccd934cc
2 files changed
+289 -290

No files matched your search

+288 -288
View File
@@ -2,9 +2,6 @@
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "o4JZ84moBKMr"
},
"source": [
"# Differentiable Fluid Simulations\n",
"\n",
@@ -31,13 +28,13 @@
" \\frac{\\partial s}{\\partial{t}} + \\mathbf{u} \\cdot \\nabla s &= 0 \n",
"\\end{aligned}$$\n",
"\n"
]
],
"metadata": {
"id": "o4JZ84moBKMr"
}
},
{
"cell_type": "markdown",
"metadata": {
"id": "pwmMK-LW_trZ"
},
"source": [
"## Formulation\n",
"\n",
@@ -60,50 +57,50 @@
"Hence, the optimization can only change $\\mathbf{u}_0$ to align $y_{t_e,i}$ with the references $y^*_{t_e,i}$ as closely as possible.\n",
"\n",
"\n"
]
],
"metadata": {
"id": "pwmMK-LW_trZ"
}
},
{
"cell_type": "markdown",
"metadata": {
"id": "omVCYK5f_tra"
},
"source": [
"## Starting the Implementation\n",
"\n",
"First, let's get the loading of python modules out of the way. By importing `phi.torch.flow`, we get fluid simulation functions that work within pytorch graphs and can provide gradients (`phi.tf.flow` would be the alternative for tensorflow).\n",
"\n"
]
],
"metadata": {
"id": "omVCYK5f_tra"
}
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "da1uZcDXdVcF",
"outputId": "66973d88-ec2c-4137-ed5b-1696848be36f"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": []
}
],
"source": [
"#!pip install --upgrade --quiet phiflow\n",
"!pip install --upgrade --quiet git+https://github.com/tum-pbs/PhiFlow@develop\n",
"\n",
"from phi.torch.flow import * \n",
"import pylab # for visualizations later on\n"
]
],
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": []
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "da1uZcDXdVcF",
"outputId": "66973d88-ec2c-4137-ed5b-1696848be36f"
}
},
{
"cell_type": "markdown",
"metadata": {
"id": "rdSTbMoaS0Uz"
},
"source": [
"## Batched simulations\n",
"\n",
@@ -111,32 +108,14 @@
"\n",
"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. This leads to the creation of an `inflow_loc` dimensions in addition to the two spatial dimensions (the `vector` part).\n"
]
],
"metadata": {
"id": "rdSTbMoaS0Uz"
}
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "mphMP0sYIOz-",
"outputId": "0c3cea02-a1b2-4d88-e756-afbaef0544b0"
},
"outputs": [
{
"data": {
"text/plain": [
"(inflow_locᵇ=4, xˢ=32, yˢ=40)"
]
},
"execution_count": 2,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"# closed domain\n",
"INFLOW_LOCATION = tensor([(12, 4), (13, 6), (14, 5), (16, 5)], batch('inflow_loc'), channel('vector'))\n",
@@ -149,61 +128,59 @@
"#BND = extrapolation.BOUNDARY # open boundaries \n",
"\n",
"INFLOW.shape"
]
],
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"(inflow_locᵇ=4, xˢ=32, yˢ=40)"
]
},
"metadata": {
"tags": []
},
"execution_count": 2
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "mphMP0sYIOz-",
"outputId": "0c3cea02-a1b2-4d88-e756-afbaef0544b0"
}
},
{
"cell_type": "markdown",
"metadata": {
"id": "VSisaefX_trh"
},
"source": [
"The last statement verifies that our `INFLOW` grid likewise has an `inflow_loc` dimension in addition to the spatial dimensions `x` and `y`. 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.some_unknown_dim.exists` will give `False`. The $^b$ superscript indicates that `inflow_loc` is a batch dimension.\n",
"\n",
"Phiflow tensors are automatically broadcast to new dimensions via their names, and hence typically no-reshaping operations are required. E.g., you can easily add or multiply tensors with differing dimensions. Below we'll multiply a staggered grid with a tensor of ones along the `inflow_loc` dimension to get a staggered velocity that has `x,y,inflow_loc` as dimensions via `StaggeredGrid(...) * math.ones(batch(inflow_loc=4))`."
]
],
"metadata": {
"id": "VSisaefX_trh"
}
},
{
"cell_type": "markdown",
"metadata": {
"id": "SDmFUYlI_tri"
},
"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."
]
],
"metadata": {
"id": "SDmFUYlI_tri"
}
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 263
},
"id": "2PccI2r5_trj",
"outputId": "8dc97fa4-6c27-4406-e0c1-cdfb9c95b5e5"
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAr4AAADiCAYAAACoVMTQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxcVfk/8M8nadJ93yhtaSkUS2UptLKJiALKoqIIyirgUvSrX0QRvxUBRQFRqbgCwg9olU1WQUSwAgVLodACXeje0tI93dKmafY8vz/uDeTcM+1MJjPJ3NzP+/XKqzknz733ZPrM5Myds9DMICIiIiLS0RW1dwNERERERNqCOr4iIiIikgjq+IqIiIhIIqjjKyIiIiKJoI6viIiIiCSCOr4iIiIikgjq+IqIiEiskBxM8mWSFSQnk/wpyfvyeL0TSTaS3EXy1HxdpwXteYFkNckZbXS9G0huIbmxLa6XT+r47oWeWG37xJLcUN5yCskqkmvbuy2SOeUtV5CszefvXOhIriJ5cobhEwFsAdDLzK7MY7OaW29mPczsWQAgOYTkUyTXkzSSI6MHkDyZ5JskK0muJfmlTC4Uyc+mr4ubfm5mnwTwzUwbHj6f6sLzlJOcSfLYDI/dD8CVAMaa2T6ZXrNQJa7j29GeWCRvIbks/GOxmORXsrkoyXvC8x/YVNfSJ5bkTwfM21+RXENyJ8nVJK9uycVIDiT5AMkdJLeTvL/pZ2Z2CYDTWv8rSWt1tLxtQrIfyc0tuSlA8hKSDZGOzIlNPzezAwDc1LpfJ1FGAFho7bsLVyOAZwF8MdUPSY4F8ACAHwPoDeBwAHNacP6m/Gz6mtrK9v7NzHoAGADgRQCPZHjcfgC2mllZSy/IQEH1NQuqMQWo4J9YACoBfBbBk+piAL8jeVxLLkDyeAAHtKaRUlDikLd3AxhjZr0AHAfgApJnteD8jwPYiOAFeRCAW1rRVikMccjbJr8EsCiL878a6chMz+IciRC+UZgR3tzZTvJdkqeFP5uC4O/dD8M3EN6bK5KfI/lOeHdzOsmDw/pLSf6jWdwyko80K68hOS6TNprZJjO7DcAbewi5BsCfzexfZlZvZlvNbEWmj0G+mFk9gPsBDCU5EABI9iZ5N8kNJNcxGNpQHD620wDsGz7WU8L4Y8K7xuUk5zZ/Exc+3jeSfAXAbgCjSI4hOY3kNpJLmt/5ZvAp3Z9I/jO8iTeL5AHNfv7hZsduarpRQrKI5CQGn5ZsJfkwyX7pfv9Ed3w7whPLzH5iZovNrNHMZgH4L4CMPr4I29IJwB8A/G+mx0j76iB5u8TMKptVNQI4MFVsivZ/CsBwAFeZ2Q4zqzOztzI5VtpPR8jb8HzHATgEwL2Z/ebSCkcDWILgDuWvANxNkuGnOvcD+FX4BuI/zQ8ieRCABwFcAWAggGcA/INkKYCXAHws7DTtC6AU4d9MkqMA9AAwL0ftPyY87/ywQ3lfJh2zZgaFHb13Sd5KsnsuGhU+Dl8BsBXA9rB6CoB6BK/DRwD4FICvh4/tafjg7vMlJIcC+CeAGwD0A/ADAI81daJDFyH4FKcngM0IOs8PILhRcS6A2xjcEW9yLoDrAfQFsBzAjWFbewL4D4I3pPuG7Xs+POZ/AXwewMfDn20H8Kd0v3+iO76huD+xmrepK4CPAHinBYd9D8DLZpbz9khexT5vw3fquwCsBdAdwYtiJo5B8LtPDd/lv0Hy47lql+RVrPOWZDGAPwL4DoBs7kwfwWCC0FKS14Y3HmTPVpvZXWbWAGAqgCEABmdw3JcB/NPMpplZHYJPhLoCOM7MVgKoADAOwAkAngOwnuQYBB2o/5pZY47aPwxBB/CLAEaHbfhDhscuDts4BMAnAYwH8JtWtudLJMsBVAH4BoCzzaye5GAApwO4wswqwyENtyLojKZyIYBnzOyZ8KbbNACzw3M0mWJm74R3l08FsMrM7g3vfL8F4DEA5zSLf8LMXm92N7rpzepnAGw0s8lmVm1mFeFNPiAYivljM1trZjUAfgrg7HTPK3V84//Eau4OAHPD66VFcjiAywBcl4e2SH7FPm/N7GYEdwOOBPBXADsyPHQYgrsRLwLYB8BkAE+SHJCrtknexD1vLwcwy8xaMk6zycsI7hQPQtAROg/AVTlqV0f1/goCZrY7/LZHBsftC2B1s2MbAawBMDSsegnAiQjy5SUA0xHkysfDcq5UAbjXzJaa2S4EY7hPT3NMU5s3mtnCsGP5LoAfIv0QnHQeNrM+CJ5zCxB0poFgmFEJgA3hJyrlAP6MIFdTGQHgnKbYMP54BM/nJmsi8UdH4i9A8PrdpPlqEbvxwf/zcAB7Gh4yAsATzc65CEAD0rym6N1m5IlFEsjyiUUy1RPrwPD7cgRPqmOR2ycWAIDkrxG8qH6iBWPkfgvgZ2aWaYdDCkeHyNswV98i+WkEH3N9P4PDqhDcPbg7LD9E8scAPgrgyVy3UXIqtnkb3k2+HB90Flok7KA3mU/yZwg6vr/IQfPEtR7AoU0FBok2HMC6sOolBHNj9kfQGW3qiB2L4I5+rsyD+8lAa8avG3J0s9LMtpCcCGA2yQcQdFJrAAwI77imswbAX83sG3u7TCT+JTM7JYvmrsGe7zyvAfBVM3ulJSfUHd/srUfwbgPAHp9YJwL4WPj9S8jPO0qQvB7BGJxPmdnOFhx6EoBfk9zID9bme5Xk+blsnxSUgsnbiE7IfIJl9I8JUpSlYymEvD0KwR2theHr5e8AHBW+fhZncT4DwBy1TVwPAziD5EkkSxAsxVUDYGb485cAfAJAVzNbi2BuzKkA+gNo0XwBkl0AdA6LncNyk3sBXEpyFMluACYBeDrD836C5AgGhgO4GWne2DNYReWSTM5vZksQfDryQzPbAODfACaT7BUOGzpgL0PI7gPwWZKfZjABrguD5deG7SH+aQAHkbyIZEn49RGG4/TTeBrAEJJXkOxMsifJo8Of3QHgRpIjwt9/IMkz051QHd/sFcQTi+SPAJwP4GQz29rC3+EgBMurjMMH42k+C+CJFp5H4qPd8zZ8Ub2MZN/wRf0oAN/GBxMW0nkCQF+SF4cvumcjGP7Qonf9EivtnrcA/gVgJD54vbwuPPe4cOhGuvOeFo6lRDgM41roE4q8CDt1FyIYT7sFwd+1z5pZbfjzpQB2IcgThDeMVgJ4JZP/y4iq8FxAMC63qlk77gHwFwCzEHxiUYPgU4NMHIEgvyvDf+fv7dhwvHt/AK+1oO2/BjCR5CAEk91KASxEMEnsUbhDF95nZmsAnAngagQT19Yg+PQiZZ/SzCoQDE87F8Gb2I0IVkbpnCo+xbGnIPg/3AhgGYLnOhC8+XwKwL9JViD43Y9OdZ7mNNQhS2a2hGTTE2sogLcReWIxmLjz/hOL5EoAm7N8YjVZHP7bdKfgJgC1AJaHHxsCwE1mlnY9SIusyRcev8XMqlIfIXFXQHn7BQQf8ZYieCH8AzKc9GFm20h+DsBtCGbwLgZwppltaWH7JCYKIW/DyTPvD9UguQNAnZllupPVSQCmkOwBYBOCu2Zat7cZMxvZ7PspCFYaaP5zNvv+ksjPfhopP4G93MQxsyGR8oQ0zWsE0D8cS/plM3su2qY9XOcnAH6S5typjvsN9jKZjeQ0BBN9Xw+rjgLwtJktThUffXzCulkIJhY3+Vb4FY2bjuDmQvTYlHeEzezEFHVLAJyxh/hL9nY9M1uA4PkTPa4RwWPUokl/zHw4qOQbyaaJGTVo9sRqx/a8/8QyMy/pRICCzNu7EcwWLjOzjJZIk+QpwLxdgqBT/7CZfbU92yLSkanjKyIiIiKJoDG+IiIiIpII6viKiIiISCKo4ysiIiIiidCmqzqQ1IBiyVq62bP5oryV1lDeShwpbyWmtpjZwL0FtMNyZtms8y3S0hWJck15K9lQ3kocKW8lrhpWp4vQUAcRERERSQR1fEVEREQkEdTxFREREZFEUMdXRERERBJBHV8RERERSQR1fEVEREQkEdTxFREREZFEUMdXRERERBJBHV8RERERSQR1fEVEREQkEdTxFREREZFEUMdXRERERBJBHV8RERERSQR1fEVEREQkEdJ2fEl2Ifk6ybkk3yF5fVg/heS7JN8Ov8blv7kimVHeShwpbyWOlLcSJ50yiKkB8Ekz20WyBMAMkv8Kf3aVmT2av+aJZE15K3GkvJU4Ut5KbKTt+JqZAdgVFkvCL8tno0RaS3krcaS8lThS3kqcZDTGl2QxybcBlAGYZmazwh/dSHIeyVtJds5bK0WyoLyVOFLeShwpbyUuMur4mlmDmY0DMAzAUSQPAfAjAGMAfARAPwD/l+pYkhNJziY5O0dtFsmI8lbiSHkrcaS8lbhg8AlFCw4grwOw28xuaVZ3IoAfmNln0hxrQHE27ZTEa4CZMdujlbfSPpS3EkfKW4mrhjlmNmFvEZms6jCQZJ/w+64ATgGwmOSQsI4APg9gQesbLJIbyluJI+WtxJHyVuIkk1UdhgCYSrIYQUf5YTN7muQLJAcCIIC3AXwzj+0UaSnlrcSR8lbiSHkrsdHioQ6tupg+wpCste6jt9ZQ3kr2lLcSR8pbiascDHUQEREREekI1PEVERERkURQx1dEREREEkEdXxERERFJBHV8RURERCQR1PEVERERkURQx1dEREREEkEdXxERERFJBHV8RURERCQR1PEVERERkURQx1dEREREEqFTezdARPKvtGSQV7dP93FOuRENXszO2rVOuWL3Ci/GUN/K1om0BjOIsby3QiT3ijOI8V+3Ze90x1dEREREEkEdXxERERFJhLQdX5JdSL5Oci7Jd0heH9bvT3IWyeUk/0ayNP/NFcmM8lbiSHkrcaXclbig2d7HPpEkgO5mtotkCYAZAL4L4PsAHjezh0jeAWCumd2e5lyW2ZiVQuKOHysq6poixn3/YI3VXoTGQbZWA8wsk8FLine truncated
"text/plain": [
"<Figure size 720x432 with 4 Axes>"
]
},
"metadata": {
"needs_background": "light",
"tags": []
},
"output_type": "display_data"
}
],
"source": [
"smoke = CenteredGrid(0, extrapolation.BOUNDARY, x=32, y=40, bounds=Box[0:32, 0:40]) # sampled at cell centers\n",
"velocity = StaggeredGrid(0, BND, x=32, y=40, bounds=Box[0:32, 0:40]) # sampled in staggered form at face centers \n",
"\n",
"def step(smoke, velocity):\n",
" smoke = advect.mac_cormack(smoke, velocity, dt=1) + INFLOW\n",
" buoyancy_force = smoke * (0, 1) >> velocity\n",
" buoyancy_force = (smoke * (0, 1)).at(velocity)\n",
" velocity = advect.semi_lagrangian(velocity, velocity, dt=1) + buoyancy_force\n",
" velocity, _ = fluid.make_incompressible(velocity)\n",
" return smoke, velocity\n",
@@ -218,24 +195,44 @@
" axes[i].imshow(smoke_final.values.numpy('inflow_loc,y,x')[i,...], origin='lower', cmap='magma')\n",
" axes[i].set_title(f\"Inflow {INFLOW_LOCATION.numpy('inflow_loc,vector')[i]}\" + (\", Reference\" if i==3 else \"\"))\n",
"pylab.tight_layout()\n"
]
],
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": [
"<Figure size 720x432 with 4 Axes>"
],
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAr4AAADiCAYAAACoVMTQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxcVfk/8M8nadJ93yhtaSkUS2UptLKJiALKoqIIyirgUvSrX0QRvxUBRQFRqbgCwg9olU1WQUSwAgVLodACXeje0tI93dKmafY8vz/uDeTcM+1MJjPJ3NzP+/XKqzknz733ZPrM5Myds9DMICIiIiLS0RW1dwNERERERNqCOr4iIiIikgjq+IqIiIhIIqjjKyIiIiKJoI6viIiIiCSCOr4iIiIikgjq+IqIiEiskBxM8mWSFSQnk/wpyfvyeL0TSTaS3EXy1HxdpwXteYFkNckZbXS9G0huIbmxLa6XT+r47oWeWG37xJLcUN5yCskqkmvbuy2SOeUtV5CszefvXOhIriJ5cobhEwFsAdDLzK7MY7OaW29mPczsWQAgOYTkUyTXkzSSI6MHkDyZ5JskK0muJfmlTC4Uyc+mr4ubfm5mnwTwzUwbHj6f6sLzlJOcSfLYDI/dD8CVAMaa2T6ZXrNQJa7j29GeWCRvIbks/GOxmORXsrkoyXvC8x/YVNfSJ5bkTwfM21+RXENyJ8nVJK9uycVIDiT5AMkdJLeTvL/pZ2Z2CYDTWv8rSWt1tLxtQrIfyc0tuSlA8hKSDZGOzIlNPzezAwDc1LpfJ1FGAFho7bsLVyOAZwF8MdUPSY4F8ACAHwPoDeBwAHNacP6m/Gz6mtrK9v7NzHoAGADgRQCPZHjcfgC2mllZSy/IQEH1NQuqMQWo4J9YACoBfBbBk+piAL8jeVxLLkDyeAAHtKaRUlDikLd3AxhjZr0AHAfgApJnteD8jwPYiOAFeRCAW1rRVikMccjbJr8EsCiL878a6chMz+IciRC+UZgR3tzZTvJdkqeFP5uC4O/dD8M3EN6bK5KfI/lOeHdzOsmDw/pLSf6jWdwyko80K68hOS6TNprZJjO7DcAbewi5BsCfzexfZlZvZlvNbEWmj0G+mFk9gPsBDCU5EABI9iZ5N8kNJNcxGNpQHD620wDsGz7WU8L4Y8K7xuUk5zZ/Exc+3jeSfAXAbgCjSI4hOY3kNpJLmt/5ZvAp3Z9I/jO8iTeL5AHNfv7hZsduarpRQrKI5CQGn5ZsJfkwyX7pfv9Ed3w7whPLzH5iZovNrNHMZgH4L4CMPr4I29IJwB8A/G+mx0j76iB5u8TMKptVNQI4MFVsivZ/CsBwAFeZ2Q4zqzOztzI5VtpPR8jb8HzHATgEwL2Z/ebSCkcDWILgDuWvANxNkuGnOvcD+FX4BuI/zQ8ieRCABwFcAWAggGcA/INkKYCXAHws7DTtC6AU4d9MkqMA9AAwL0ftPyY87/ywQ3lfJh2zZgaFHb13Sd5KsnsuGhU+Dl8BsBXA9rB6CoB6BK/DRwD4FICvh4/tafjg7vMlJIcC+CeAGwD0A/ADAI81daJDFyH4FKcngM0IOs8PILhRcS6A2xjcEW9yLoDrAfQFsBzAjWFbewL4D4I3pPuG7Xs+POZ/AXwewMfDn20H8Kd0v3+iO76huD+xmrepK4CPAHinBYd9D8DLZpbz9khexT5vw3fquwCsBdAdwYtiJo5B8LtPDd/lv0Hy47lql+RVrPOWZDGAPwL4DoBs7kwfwWCC0FKS14Y3HmTPVpvZXWbWAGAqgCEABmdw3JcB/NPMpplZHYJPhLoCOM7MVgKoADAOwAkAngOwnuQYBB2o/5pZY47aPwxBB/CLAEaHbfhDhscuDts4BMAnAYwH8JtWtudLJMsBVAH4BoCzzaye5GAApwO4wswqwyENtyLojKZyIYBnzOyZ8KbbNACzw3M0mWJm74R3l08FsMrM7g3vfL8F4DEA5zSLf8LMXm92N7rpzepnAGw0s8lmVm1mFeFNPiAYivljM1trZjUAfgrg7HTPK3V84//Eau4OAHPD66VFcjiAywBcl4e2SH7FPm/N7GYEdwOOBPBXADsyPHQYgrsRLwLYB8BkAE+SHJCrtknexD1vLwcwy8xaMk6zycsI7hQPQtAROg/AVTlqV0f1/goCZrY7/LZHBsftC2B1s2MbAawBMDSsegnAiQjy5SUA0xHkysfDcq5UAbjXzJaa2S4EY7hPT3NMU5s3mtnCsGP5LoAfIv0QnHQeNrM+CJ5zCxB0poFgmFEJgA3hJyrlAP6MIFdTGQHgnKbYMP54BM/nJmsi8UdH4i9A8PrdpPlqEbvxwf/zcAB7Gh4yAsATzc65CEAD0rym6N1m5IlFEsjyiUUy1RPrwPD7cgRPqmOR2ycWAIDkrxG8qH6iBWPkfgvgZ2aWaYdDCkeHyNswV98i+WkEH3N9P4PDqhDcPbg7LD9E8scAPgrgyVy3UXIqtnkb3k2+HB90Flok7KA3mU/yZwg6vr/IQfPEtR7AoU0FBok2HMC6sOolBHNj9kfQGW3qiB2L4I5+rsyD+8lAa8avG3J0s9LMtpCcCGA2yQcQdFJrAAwI77imswbAX83sG3u7TCT+JTM7JYvmrsGe7zyvAfBVM3ulJSfUHd/srUfwbgPAHp9YJwL4WPj9S8jPO0qQvB7BGJxPmdnOFhx6EoBfk9zID9bme5Xk+blsnxSUgsnbiE7IfIJl9I8JUpSlYymEvD0KwR2theHr5e8AHBW+fhZncT4DwBy1TVwPAziD5EkkSxAsxVUDYGb485cAfAJAVzNbi2BuzKkA+gNo0XwBkl0AdA6LncNyk3sBXEpyFMluACYBeDrD836C5AgGhgO4GWne2DNYReWSTM5vZksQfDryQzPbAODfACaT7BUOGzpgL0PI7gPwWZKfZjABrguD5deG7SH+aQAHkbyIZEn49RGG4/TTeBrAEJJXkOxMsifJo8Of3QHgRpIjwt9/IMkz051QHd/sFcQTi+SPAJwP4GQz29rC3+EgBMurjMMH42k+C+CJFp5H4qPd8zZ8Ub2MZN/wRf0oAN/GBxMW0nkCQF+SF4cvumcjGP7Qonf9EivtnrcA/gVgJD54vbwuPPe4cOhGuvOeFo6lRDgM41roE4q8CDt1FyIYT7sFwd+1z5pZbfjzpQB2IcgThDeMVgJ4JZP/y4iq8FxAMC63qlk77gHwFwCzEHxiUYPgU4NMHIEgvyvDf+fv7dhwvHt/AK+1oO2/BjCR5CAEk91KASxEMEnsUbhDF95nZmsAnAngagQT19Yg+PQiZZ/SzCoQDE87F8Gb2I0IVkbpnCo+xbGnIPg/3AhgGYLnOhC8+XwKwL9JViD43Y9OdZ7mNNQhS2a2hGTTE2sogLcReWIxmLjz/hOL5EoAm7N8YjVZHP7bdKfgJgC1AJaHHxsCwE1mlnY9SIusyRcev8XMqlIfIXFXQHn7BQQf8ZYieCH8AzKc9GFm20h+DsBtCGbwLgZwppltaWH7JCYKIW/DyTPvD9UguQNAnZllupPVSQCmkOwBYBOCu2Zat7cZMxvZ7PspCFYaaP5zNvv+ksjPfhopP4G93MQxsyGR8oQ0zWsE0D8cS/plM3su2qY9XOcnAH6S5typjvsN9jKZjeQ0BBN9Xw+rjgLwtJktThUffXzCulkIJhY3+Vb4FY2bjuDmQvTYlHeEzezEFHVLAJyxh/hL9nY9M1uA4PkTPa4RwWPUokl/zHw4qOQbyaaJGTVo9sRqx/a8/8QyMy/pRICCzNu7EcwWLjOzjJZIk+QpwLxdgqBT/7CZfbU92yLSkanjKyIiIiKJoDG+IiIiIpII6viKiIiISCKo4ysiIiIiidCmqzqQ1IBiyVq62bP5oryV1lDeShwpbyWmtpjZwL0FtMNyZtms8y3S0hWJck15K9lQ3kocKW8lrhpWp4vQUAcRERERSQR1fEVEREQkEdTxFREREZFEUMdXRERERBJBHV8RERERSQR1fEVEREQkEdTxFREREZFEUMdXRERERBJBHV8RERERSQR1fEVEREQkEdTxFREREZFEUMdXRERERBJBHV8RERERSQR1fEVEREQkEdJ2fEl2Ifk6ybkk3yF5fVg/heS7JN8Ov8blv7kimVHeShwpbyWOlLcSJ50yiKkB8Ekz20WyBMAMkv8Kf3aVmT2av+aJZE15K3GkvJU4Ut5KbKTt+JqZAdgVFkvCL8tno0RaS3krcaS8lThS3kqcZDTGl2QxybcBlAGYZmazwh/dSHIeyVtJds5bK0WyoLyVOFLeShwpbyUuMur4mlmDmY0DMAzAUSQPAfAjAGMAfARAPwD/l+pYkhNJziY5O0dtFsmI8lbiSHkrcaS8lbhg8AlFCw4grwOw28xuaVZ3IoAfmNln0hxrQHE27ZTEa4CZMdujlbfSPpS3EkfKW4mrhjlmNmFvEZms6jCQZJ/w+64ATgGwmOSQsI4APg9gQesbLJIbyluJI+WtxJHyVuIkk1UdhgCYSrIYQUf5YTN7muQLJAcCIIC3AXwzj+0UaSnlrcSR8lbiSHkrsdHioQ6tupg+wpCste6jt9ZQ3kr2lLcSR8pbiascDHUQEREREekI1PEVERERkURQx1dEREREEkEdXxERERFJBHV8RURERCQR1PEVERERkURQx1dEREREEkEdXxERERFJBHV8RURERCQR1PEVERERkURQx1dEREREEqFTezdARPKvtGSQV7dP93FOuRENXszO2rVOuWL3Ci/GUN/K1om0BjOIsby3QiT3ijOI8V+3Ze90x1dEREREEkEdXxERERFJhLQdX5JdSL5Oci7Jd0heH9bvT3IWyeUk/0ayNP/NFcmM8lbiSHkrcaXclbig2d7HPpEkgO5mtotkCYAZAL4L4PsAHjezh0jeAWCumd2e5lyW2ZiVQuKOHysq6poixn3/YI3VXoTGQbZWA8wsk8FLine truncated
},
"metadata": {
"needs_background": "light",
"tags": []
}
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 263
},
"id": "2PccI2r5_trj",
"outputId": "8dc97fa4-6c27-4406-e0c1-cdfb9c95b5e5"
}
},
{
"cell_type": "markdown",
"metadata": {
"id": "V1tfVLqJ_trj"
},
"source": [
"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",
"(As a small side note: 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 PyTorch or TensorFlow tensor is also passed to the same operation. Thus, it is a good idea to verify that tensors are using the right backend once in a while, e.g., via `GRID.values.default_backend`.)\n"
]
],
"metadata": {
"id": "V1tfVLqJ_trj"
}
},
{
"cell_type": "markdown",
"metadata": {
"id": "vjIJq_Uz_trk"
},
"source": [
"## Gradients\n",
"\n",
@@ -250,27 +247,14 @@
"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 appropriate velocity tensor, we initialize a `StaggeredGrid` with a tensor of zeros along the `inflow_loc` batch dimension. As the staggered grid already has `y,x` and `vector` dimensions, this gives the desired four dimensions, as verified by the print statement below.\n",
"\n",
"Phiflow provides a unified API for gradients across different platforms by using functions that need to return a loss values, in addition to optional state values. It uses a loss function based interface, for which we define the `simulate` function below. `simulate` computes the $L^2$ error outlined above and returns the evolved `smoke` and `velocity` states after 20 simulation steps.\n"
]
],
"metadata": {
"id": "vjIJq_Uz_trk"
}
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "_j7fUGvQ_trl",
"outputId": "f06dea67-cc48-4af1-a963-46500f90d138"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Velocity dimensions: (inflow_locᵇ=4, xˢ=32, yˢ=40, vectorᵛ=2)\n"
]
}
],
"source": [
"initial_smoke = CenteredGrid(0, extrapolation.BOUNDARY, x=32, y=40, bounds=Box[0:32, 0:40])\n",
"initial_velocity = StaggeredGrid(math.zeros(batch(inflow_loc=4)), BND, x=32, y=40, bounds=Box[0:32, 0:40])\n",
@@ -285,84 +269,97 @@
" #loss = field.l2_loss(diffuse.explicit(smoke - field.stop_gradient(smoke.inflow_loc[-1]), 1, 1, 10))\n",
" \n",
" return loss, smoke, velocity\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "pjqeVedM_trm"
},
"source": [
"Phiflow's `field.functional_gradient()` function is the central function to compute gradients. Next, we'll use it 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]`. (Phiflow also has a `field.spatial_gradient` function which instead computes derivatives of tensors along spatial dimensions, like `x,y`.)\n",
"\n",
"`functional_gradient` generates a gradient function. As a demonstration, the next cell evaluates the gradient once with the initial states for smoke and velocity. The last statement prints a summary of a part of the resulting gradient tensor.\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
],
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"Velocity dimensions: (inflow_locᵇ=4, xˢ=32, yˢ=40, vectorᵛ=2)\n"
]
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "grYi0l-g_trn",
"outputId": "93dae39b-1b31-429a-cf81-f615e23bbf7d"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Some gradient info: StaggeredGrid[(inflow_locᵇ=4, xˢ=32, yˢ=40, vectorᵛ=2), size=(32, 40), extrapolation=0]\n",
"(xˢ=31, yˢ=40) float32 -17.366662979125977 < ... < 14.014090538024902\n"
]
}
"id": "_j7fUGvQ_trl",
"outputId": "f06dea67-cc48-4af1-a963-46500f90d138"
}
},
{
"cell_type": "markdown",
"source": [
"Phiflow's `field.functional_gradient()` function is the central function to compute gradients. Next, we'll use it 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]`. (Phiflow also has a `field.spatial_gradient` function which instead computes derivatives of tensors along spatial dimensions, like `x,y`.)\n",
"\n",
"`functional_gradient` generates a gradient function. As a demonstration, the next cell evaluates the gradient once with the initial states for smoke and velocity. The last statement prints a summary of a part of the resulting gradient tensor.\n"
],
"metadata": {
"id": "pjqeVedM_trm"
}
},
{
"cell_type": "code",
"execution_count": 5,
"source": [
"sim_grad = field.functional_gradient(simulate, wrt=[1], get_output=False)\n",
"(velocity_grad,) = sim_grad(initial_smoke, initial_velocity)\n",
"\n",
"print(\"Some gradient info: \" + format(velocity_grad))\n",
"print(format(velocity_grad.values.inflow_loc[0].vector[0])) # one example, location 0, x component, automatically prints size & content range\n"
]
],
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"Some gradient info: StaggeredGrid[(inflow_locᵇ=4, xˢ=32, yˢ=40, vectorᵛ=2), size=(32, 40), extrapolation=0]\n",
"(xˢ=31, yˢ=40) float32 -17.366662979125977 < ... < 14.014090538024902\n"
]
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "grYi0l-g_trn",
"outputId": "93dae39b-1b31-429a-cf81-f615e23bbf7d"
}
},
{
"cell_type": "markdown",
"metadata": {
"id": "-drJBAe2_tro"
},
"source": [
"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",
"We could use this to take a look at the content of the computed gradient with regular plotting functions, e.g., by converting the x component of one of the simulations to a numpy array via `velocity_grad.values.inflow_loc[0].vector[0].numpy('y,x')`. However, below we'll use phiflow's `view()` function instead. It automatically analyzes the grid content and provides UI buttons to choose different viewing modes. You can use them to show arrows, single components of the 2-dimensional velocity vectors, or their magnitudes. (Because of its interactive nature, the corresponding image won't show up outside of Jupyter, though.)\n"
]
],
"metadata": {
"id": "-drJBAe2_tro"
}
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "2LTHHjtZ_tro"
},
"outputs": [],
"source": [
"# neat phiflow helper function:\n",
"view(velocity_grad)\n"
]
],
"outputs": [],
"metadata": {
"id": "2LTHHjtZ_tro"
}
},
{
"cell_type": "markdown",
"metadata": {
"id": "J7IgFG4o_trp"
},
"source": [
"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 initial round inflow positions shown as positive and negative regions around the circular shape of the inflow. The ones for the larger distances on the left are also noticeably larger."
]
],
"metadata": {
"id": "J7IgFG4o_trp"
}
},
{
"cell_type": "markdown",
"metadata": {
"id": "3mpyowRYUSS4"
},
"source": [
"\n",
"\n",
@@ -375,22 +372,27 @@
"The following cell implements a simple steepest gradient descent optimization: it re-evaluates the gradient function, and iterates several times to optimize $\\mathbf{u}_0$ 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()` in our case) 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"
]
],
"metadata": {
"id": "3mpyowRYUSS4"
}
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "0U8h7_-b_trq",
"outputId": "11df5469-5c45-41b4-ba21-81114d2a3431"
},
"source": [
"sim_grad_wloss = field.functional_gradient(simulate, wrt=[1], get_output=True) # if we need outputs...\n",
"\n",
"LR = 1e-03 \n",
"for optim_step in range(80): \n",
" (loss, _smoke, _velocity), (velocity_grad,) = sim_grad_wloss(initial_smoke, initial_velocity)\n",
" initial_velocity = initial_velocity - LR * velocity_grad\n",
" if optim_step<3 or optim_step%10==9: print('Optimization step %d, loss: %f' % (optim_step, np.sum(loss.numpy()) ))\n"
],
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"name": "stdout",
"text": [
"Optimization step 0, loss: 1193.145020\n",
"Optimization step 1, loss: 1165.816650\n",
@@ -405,53 +407,28 @@
]
}
],
"source": [
"sim_grad_wloss = field.functional_gradient(simulate, wrt=[1], get_output=True) # if we need outputs...\n",
"\n",
"LR = 1e-03 \n",
"for optim_step in range(80): \n",
" (loss, _smoke, _velocity), (velocity_grad,) = sim_grad_wloss(initial_smoke, initial_velocity)\n",
" initial_velocity = initial_velocity - LR * velocity_grad\n",
" if optim_step<3 or optim_step%10==9: print('Optimization step %d, loss: %f' % (optim_step, np.sum(loss.numpy()) ))\n"
]
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "0U8h7_-b_trq",
"outputId": "11df5469-5c45-41b4-ba21-81114d2a3431"
}
},
{
"cell_type": "markdown",
"metadata": {
"id": "PI3xY-7m_trr"
},
"source": [
"The loss should have gone down significantly, from above 1100 to below 700, and now we can also visualize the initial velocities that were obtained in the optimization. \n",
"\n",
"The following images show the resulting three initial velocities in terms of their x (first set), and y components (second set of images). We're skipping the fourth set with `inflow_loc[0]` because it only contains zeros.\n"
]
],
"metadata": {
"id": "PI3xY-7m_trr"
}
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 316
},
"id": "FdL_orrX_trr",
"outputId": "b716256c-9b8a-40d9-fc01-d75e261a728c"
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAEXCAYAAACwK7IRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO29e7QlVX3v+/3uR+9+0C9oxBY4gtrRYM4I3LRoNBoOzzbRYM41CvF68Rw86B3hxjz0iMmIRBJzUJMQx/GVVhDMURHx1ceBkgbhoDcJ9kY5CCjStCDd0mB3Qz/o5977d/+o2rhq1tx71apVtVbNWd9Pjxq9ZtWcs+Zau741Z836/eaPZgYhhBBCCCFEwsiwGyCEEEIIIUST0ABZCCGEEEKIDjRAFkIIIYQQogMNkIUQQgghhOhAA2QhhBBCCCE6GBt2A4Sog/POO9127tzdU5m77vrxzWa2rqYmCSHmQHoVIizaoFkNkEWU7Ny5G3d+9x97KjM2+h9W1dQcIcQ8SK9ChEUbNKsBsogTAzAzM+xWCCGKIL0KERYt0KwGyCJSLHrxChEP0qsQYRG/ZuWkJ+JlZqa3rQAk15F8gORmkpd5jr+K5PdITpF8vXNsmuTd6bahom8pRBzUoFchRI1ErlnNIIs4MQAVh1EnOQrgowDOAbAVwCaSG8zs/o5sPwXwFgDv9FRxwMxOrbRRQsRADXoVQtRICzSrAbKIlFpe/5wOYLOZbQEAktcDOB/AMwNkM3s4PRbe47IQQyP+17VCxEX8mtUAWcRL7+JdRXKyI73ezNZ3pI8H8GhHeiuAl/ZQ/8K0/ikAV5rZV3ttoBDREnlnK0R0RK5ZDZBFnJTzsN1hZmtraM0szzWzbSSfB+BbJH9gZg/VeD4hwqAFHvFCREULNKsBsoiUWl7/bANwYkf6hHRfsRaZbUv/30LydgCnAdAAWYgWvK4VIi7i16wGyCJO6nm63QRgDcmTkQyMLwDw+0UKklwJYL+ZHSK5CsArAHyw6gYKESQtmI0SIipaoFkNkEWkGGjVitfMpkheCuBmAKMArjGz+0heAWDSzDaQfAmArwBYCeC1JN9nZi8G8MsA/jF13htBYoN8/xynEqJlVK9XIUSdxK9ZDZBFvNTwdGtmNwG4ydn33o7Pm5CYXrjl/gXAv6+8QULEQuSzUUJER+Sa1QBZxIkBmIl7jUYhokF6FSIsWqBZDZBFpMTvQCBEPEivQoRF/JrVAFnESQscCISIBulViLBogWZHht0AIWrDZnrbhBDDowa9klxH8gGSm0le5jn+KpLfIzlF8vXOsWmSd6fbhoq+pRDxEHkfqxlkESnxv/4RIh6q1yvJUQAfBXAOkqiXm0hucFaP+SmAtwB4p6eKA2Z2aqWNEiIa4u9jNUAWcdICBwIhoqEevZ4OYLOZbQEAktcDOB/AMwNkM3s4PRZ3Ty9E1bSgj9UAWURK/E+3QsRDKb2uIjnZkV5vZus70scDeLQjvRXAS3uof2Fa/xSSdcu/2msDhYiX+PtYDZBFvEQuXiGione97jCztXU0JeW5ZraN5PMAfIvkD8xMoeGFmCXyPlYDZBEnBjBy8QoRDfXodRuAEzvSJ6T7ijXJbFv6/xaStwM4DYAGyEIArehjtYqFiBQDrMdNCDEkatHrJgBrSJ5McgGACwAUWo2C5EqSE+nnVQBegQ7bZSFE/H2sZpBFvET+dCtEVFSsVzObInkpgJsBjAK4xszuI3kFgEkz20DyJQC+AmAlgNeSfJ+ZvRjALwP4x9R5bwSJDbIGyEJ0EnkfqwGyiJMWLGIuRDTUpFczuwnATc6+93Z83oTE9MIt9y8A/n3lDRIiFlrQx8rEQgghhBBCiA40gywixaJfo1GIeJBehQiL+DWrGWQRJ7Ovf3rZhBDDQXoVIixq0GyfoeEvIvlgul1UxVds7QCZ5Deq+hHnOccZJLfWWPcMyX0k19Vxjh7b8y2SB0l+Z9hteQZ1uNEgvVbenmtJHqjr+5ZCeo0Kabby9jxE8jDJ/zHstjxDhZrtCA3/agCnALiQ5ClOttnQ8J9zyh4N4HIkgYBOB3A5yZX9fr2oBsgkHyZ5dpG8ZvZqM7uu7jYVheRrSW5P/9Cz+84nuY3k8jmK/czMjjKzb6b5V5PcQPJnJI3kSc45/jZ9utpL8kck/++Sbb0mrf8Fs/vM7EwAby9TXz3EvwRN6EivXfX6QZKPktxD8hGSf9ZjG48l+TmSu0k+SfKzs8fM7C1IOqKGIL2GgDQ7v2Y76j2a5M97mTAi+RaS0+mAfHY7Y/a4mT0fwN8Ura9+KtfsM6HhzewwgNnQ8L84o9nDZnYPAHe0fR6AjWa2y8yeBLARQN8PNVENkEPGzP4ngG8BuAoASK4A8HEA/4+Z7S5YzQyAbwL4P+c4/jSA1wJYDuAiAB8m+fJe2knyNwA8v5cyQ0GvbEWNDEivVwN4kZktA/ByAG8i+R97aOaXAWwH8O8APAvA3/ZQdrBIr6JmBqTZWT4A4Iclmvmv6YB8dru9RB2DoZxmV5Gc7Ngu6ajRFxr++IKt6afsnEQ7QE6fxr6Tzpo+SfInJF/dcfx2km8tUM8bSU46+/6Y5Ib080R6jp+SfJzkJ0guKtnsPwTwapLnIRHx/zKzQgvbA4CZPW5mH0OyQL7v+OVm9iMzmzGzOwF8G8CvF62f5BiA/w7g/y1aZqjMWG+bGBrSq/f4A2b2dMeuGQAv8OV1IXkukihy7zKz3WZ2xMy+X7RtQ0F6DQpp1k866fQrAD5dso3h0Ltmd5jZ2o5t/bC/wnxEO0BOeSmABwCsAvBBAFeTZI91/E8ALyS5pmPf7+MXNjBXAvglAKci6byOB/BelMDMdgB4B4DPAngNEjHXQnqDeQmA+3oo9scA7khfcTQc04xUeEivDiQvI7kPyYzIEji2d/PwMiS/5XUkd5LcRPI3q25fdUivgSLNdsDEjvYjAC5FMsfaK6eR3EHyxyT/Ip2UaiiVa7af0PB9hZWfi9gHyI+Y2SfNbBrAdQBWAziulwrMbD+ArwG4EABSEb8IwIb0RnAJgD9ObV/2IrERuqCPNv8bEhOIfzazn/dRTzc+AeB/I4ky1RWSJwJ4G0remAaOXtmGiPTqYGZXAlgK4P8A8E8Air4KPgHAuQBuA/BsAH8H4GtMwiY3D+k1VKTZLH8I4E4zu6tE2TuQzDw/C4kJx4UA3lVh26qles2WDg2PZBxzLpMQ8SuR3PsKjW3mI/YB8vbZD6kIAeCoEvV8Dql4kTzZfjWt71gAiwHcRfIpkk8hsU86tnyTsR7AZwD8FsnC5g+9QPJDSIT4BrPC3i7/AOCKHmy1ho9e2YaG9OrBEr4P4ACA9xUsdgDAw2Z2dWpecT0SG71X1NHGSpBeQ0SaTSH5HCQD5D8vUz51TvtJagL5AwBXAHh9t3JDpULNmtkUkpn3m5HYb98wGxqe5O8AAMmXMFm15PeQhIK/Ly27C8BfIRlkb0IyVtnV79dr8PR9o9gI4FiSpyIR8R+n+3cg6YhebGZ9T+eTvBjJa4LXAJgE8CmSp6UenZVA8n1IvNd/08z29FD0LAC/QfKDHfv+leQ7zKzoa98BYoBplqmlRKNXhzEUd5C9B4lDbicNHlVKry0nBs2ejmQG/f7UymQRgEUktwM4Pp1l7wUD0Ku5ygCpXrNlQ8Onx64BcE2V7Yl9BrkSzOwIgC8C+BCAo5GIGWY2A+CTAK4i+SwAIHl86gDQE+nT54cA/BczO4TEBGInenwaJbkQwESanEjTs8feg+Tp/Gwz29ljE38JwK8isQM7Nd33WgBf6bGewWDQjFRLiUGvJEdIvi19ZUiSpwP4AwC3Fqz6KwBWMlk8f5TJovonAPj/emnfwJBeW00MmgXwDQAn4Rd95HsBfB/AqUUGxyRfTfK49POLAPwFEtOTZtICzWqADIDkK1NHmPn4HICzAXwxfRUwy7sBbAbwbyT3ALgFwAvnOM83OPdaph8DcL2ZfRtIXqsC+C8A/ojki4t/GxwAMPtdfpSmZ/kbJEs+beYv1lkstLaqmT1hZttnt3T3DjM7MG/BYSKbxihpkV5/F8BDAPYC+B9IVpD570UqTV8v/g6AdyKxW74MwPmpk1IzkV6jpQ2aNbNDTh+5G8CRjv6yG2cBuIfk00hmUb+MRq177CFyzbK4CapoEiRfhcRW5xCAN5pZ3wbpfbZnIxLP+e+a2VnDbAsArP2l423Tx97WU5mRcy6/y8zW1tQk0WIaqNerkdjxPWFmhZaOqxPpVTSNBmr2ASQreNxgZv95mG0B2qFZ2SAHipndgcTGqRGY2TnDbkMWC/KJVcRJA/V6MYCLh92OXyC9imbRQM16Z82HR/ya1QBZxEuANk9CtBbpVYiwiFyzhW2QU0eP75P8epo+meSdJDeT/EK6bp0QzcCQeNj2skWE9CqCouV6BaRZERgt0GwvM8jvQLI23bI0/QEAV5nZ9SQ/geR13cfnq2D5+GJ71oLlpRoqxCyb92/fYWZd1sEM02u2QvrW67KxxXbshPQq+mOL9FoU9bGiEaiPTSg0QCZ5AoDfBvB+AH/CZJG/M5EsGQYkEXT+El3E+6wFy/HhX/lPpRsrBAD89nf/2yOFMkYu3rmoSq/HTizHB3+5QWaqIkhef9f7pdcuqI8VTUJ9bELRGeR/APBfkYQ8BYBjADzVsRTLViTelTlIXoIkVCSOXbDMl0WI6pkNg9lOKtHrKulVDIp26xVQHytCowWa7TpAJvkaJEsB3UXyjF5PYGbrkYR2xJolqzOPGzPWPUhMkZ/fKqrHfRiyAkFs3DJFvtMgn7mKhOEZYfcWjTgV0fkW7nFvHZ59dM7ty1OkfV4if7r1UaVen+/otVD5EpopQlW6KnLuur5DVRTRmqvPInW4u8rcF4qce05aqFeg3X2s91xd9FekjiLfuy7K9lfd+tgiZdTHVksRJ71XAPgdkg8DuB7Line truncated
"text/plain": [
"<Figure size 720x288 with 6 Axes>"
]
},
"metadata": {
"needs_background": "light",
"tags": []
},
"output_type": "display_data"
}
],
"source": [
"fig, axes = pylab.subplots(1, 3, figsize=(10, 4))\n",
"for i in range(INFLOW.shape.get_size('inflow_loc')-1):\n",
@@ -459,34 +436,34 @@
" axes[i].set_title(f\"Ini. vel. X {INFLOW_LOCATION.numpy('inflow_loc,vector')[i]}\")\n",
" pylab.colorbar(im,ax=axes[i])\n",
"pylab.tight_layout()\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 318
},
"id": "mLhwxsZd_trs",
"outputId": "5bb43835-6e21-4d97-bb07-17b17215bbc8"
},
],
"outputs": [
{
"output_type": "display_data",
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAscAAAEYCAYAAACwdltJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO29e7glVXnu+77r1legbwgtjUC0DRGPabVF8xgNctE2FzHnKKKJwgkcwj7bnev2gHoeL0T3wZgE95O4jW0goNEAwahsg7KBwFZ3ItIgUS5BWgRtrvYV+rZ6Xb7zx6yFs0aNteaYNatqVo16f/3U07NqjlHjm7XqrTGq6vvGRzODEEIIIYQQAhgZtgFCCCGEEELUBQ2OhRBCCCGESNDgWAghhBBCiAQNjoUQQgghhEjQ4FgIIYQQQogEDY6FEEIIIYRI0OBYCCGEEEIMDZKbSD5AcivJiz3f/xHJ+0h+j+QtJI/r+u4ckg8myzmF2KN5joUQQgghxDAgOQrgBwDOALANwB0A3m5m93WVeR2A281sP8n/AOAUM3sbyVUAtgDYCMAA3Ang5Wa2axCb9ORYCCGEEEIMi5MBbDWzh8zsEICrAZzZXcDMbjWz/cnqtwGsSz6/AcBNZrYzGRDfBGDToAaNDboDIerOG95wsu3YsaevOnfe+YMbzWxggQkhwpFWhag/Jej0GAA/6VrfBuCVC+zuPABfW6DuMX0Z50GDYxE9O3bswe3f+XRfdcZGX7emJHOEEPMgrQpRf3Lq9ESSW7o2bTazzf22TfK30XGh+JV+6/aDBscifgzA7OywrRBC9EJaFaL+5NPpdjPbOM93jwI4tmt9XbItBcnTAbwfwK+Y2WRX3VOcurf1a5yLBseiBZg6XCEagbQqRP0pXKd3AFhP8gR0BrtnA3hHdwGSLwXwaQCbzOyprq9uBPBfSK5M1l8P4L2DGqTBsWgH6nCFaAbSqhD1p0Cdmtk0yXejM9AdBXCFmd1L8hIAW8zsegAfB7AcwD+QBIAfm9mbzGwnyT9BZ4ANAJeY2c5BbdLgWMSPAdCUhULUH2lViPpTgk7N7AYANzjbPtD1+fQF6l4B4Ioi7dHgWLQAvaoVohlIq0LUn/h1qsGxaAeRC1mIaJBWhag/ketUg2MRP4qAF6IZSKtC1J8W6FSDY9EC4n8FJEQcSKtC1J/4darBsYifFtzlChEF0qoQ9acFOtXgWLQAAy1uIQsRB9KqEPUnfp1qcCzaQeR3uUJEg7QqRP2JXKcaHIv4MQCzmjtViNojrQpRf1qgUw2ORQuIP3hAiDiQVoWoP/HrVINjET8tCB4QIgqkVSHqTwt0qsGxaAeRBw8IEQ3SqhD1J3KdanAsWkD8r4CEiANpVYj6E79ONTgW8dOC4AEhokBaFaL+tECnGhyLFhD/Xa4QcSCtClF/4tepBseiHUQuZCGiQVoVov5ErlMNjkX8GMDIhSxEFEirQtSfFuh0ZNgGDBuSXyN5TsltnEJyWwn7PZ6kkdxL8oKi95/DnitJHijjtw6GAdbnImqFdFqoPR8muS+xqWYPSKTVpiOtFmrPP5M8SPJbw7YlTfw6jXJwTPJhkqeHlDWzN5rZVWXbFALJo0huJ3mKs/0KklcvUHWFmW1Oyk6QvC45BubZ13tI3kPyGZI/IvmenLZ+INn/s8fZzM4F8MY8+yud2dn+FlE60umCOv1Dkg+RfJrkYyQv62cgS3Ipyf+W2LmH5DfmvjOzDwI4KXRflSOt1g5pdX6tdu1zguT9/Qzak0H+bDIYn1uevbEws1MBXBi6v0qJXKdRDo6bipk9CeAPAXyG5BIAIHkagF8H8J/62NW3APw2gCc83xHAuwCsBLAJwLtJnt2PnSSfD+CtAB7vp97QmJuwPGIhi+qoSKfXA3iZmR0O4MUAfhHA7/Wx780AVgH4heT/P+yj7vCQVkWBVKTVOd4D4Kc5zHzMzJZ3LbW4sViQFug0+sExyXNJfovkn5HclTwtfWPX97eRPD9gP28jucXZ9ockr08+L0ra+DHJJ0n+9ZwY+8HMPgfgAQCXJPU/DeD3zCxIdGZ2yMw+YWbfAjDj+f5PzewuM5s2swcAfAXAq/s085MALgJwqM96Q8I60870s4hKkU4z3//QzHbP/QQAswBeELJvkicCeBOAC8zsp2Y2Y2Z3htQdPtJq3ZFWvb/lBHQGz/9fv/Y1k/h1Gv3gOOGV6IhjDYA/BXA5Sfa5j/8O4OdJru/a9g4AX0g+XwrghQA2oNOJHQPgAzntvRDA7wC4GsA9ZrbQ65/cJMfgNQDu7aPOWwFMmtkNZdhUCi24y40E6bQLku8g+TSA7eg8Of50YNWTATwC4MPJK+Xvk/w/irStNKTVpiCtpvlLAO8DcCBH3eckg/8fJe5Tywq2rXhK0CnJTSQfILmV5MWe719L8i6S0yTf4nw3Q/LuZLm+iJ/YlsHxI2b2GTObAXAVgLUAjupnB2a2H52nrG8HgETQJwK4PrkoXADgD81sp5k9A+C/AOjLXaGrrW3oXAROB/Af8uwjkA+hcw78bUhhkoeh87t+v0SbykEdbhOQTtP7/0LiVvFCAH8N4MnAquvQccXYA+C5AN4N4CqSv1C0jaUgrTYBaTWB5G8CGDWzL+Wo/u/oDP7XAjgVwMsB/EWB5pVHgTolOYrOG+k3AngRgLeTfJFT7McAzsXPbp66OWBmG5LlTYP/uPYMjp/1E0oECQDLc+znC0iEjM4d7peT/R0JYCmAO0nuJrkbwNeT7Xm5F8AuMyvFr5fku9HxPf41M5sMrPYhAJ8zs4fLsKk84o+sjQTp1IOZPZi0898CqxwAMAXgI8kr4f8J4FYAry/JxAKRVhuCtAogecr7p+gvHuBZzOwJM7vPzGbN7EcA/h8ADXjLU7hOTwaw1cweMrND6DzhPzPVotnDZvY9dFzMSqdm0/jUnpsAHElyAzqCngty2Y5Oh3SSmT06LONCIfk7AC4G8NrkjjqU0wCsI/l/J+tHAriW5MfM7GNF21kYc6+ARFuIQqcOYwCeH1j2e55tzRhFSqtto+laXQ/geADfTLxKJgAcQfIJAK/K8SDJ0ISHlsXr9BgAP+la34aO604oixP/9WkAl5rZlwc1qP5/hBphZlMA/gHAx9GJAL8p2T4L4DMALiP5HAAgeQzJNwzDziSQYXGyOkFy8Zw/GMnfQuf11Blm9lCfuz4Nnde1G5LlMQC/i87rkHoTefCA+BmR6PT8LhtfBOC9AG4J3PU30HkF+V6SYyRfDeB1AG4s9heUhLTaGiLQ6j0AjsXP+sTz0XF/2oD0YG++/b6O5HHscCw6ftZfKeVHFE3/Ol1DckvXUuQ80seZ2UZ03j58gp0ZtQZCg+MuSL6G5N4exb6Ajt/SP5jZdNf2iwBsBfDtJIjmZgA/P087XyP5viJsnocH0LnrPgadDvEAgOOS7z4CYDWAO/izeRX/OmSnZrYjeQ30hJk9gU7k7i4z63XMhozJjzEiWqLTVwP4Psl9AG5IliBbkgHHmQB+FR2/488AeJeZ/Xuh1peCtBoTsWs1mfWpu0/cCWA2WffObOHwUgD/AmBf8v/3kdNFo1py6XS7mW3sWjZ37fBRdG4y5liXbAuzJnm7kDzwuw2d4zoQNPlsNRaSx6Ej2oMA3mNmnxmyPZejM//xU2YWNO1UFWw8cZ3dcXk/U1oCI7988Z3JnagQA1FDnX4QwB8BWARgWWAnXgnSqhgmNdTqTQBeBeA7ZnbaMG3ppmidspPg6AfovJ1+FMAdAN5hZpmZtEheCeCrZnZdsr4SwH4zmyS5BsC/AjjTzO7ry0AH+Rw3GDN7BMDingUrwszOA3DesO3wotevYkjUUKcfBvDhYdsxL9KqGBI11OoZw7ZhXgrUqZlNJ5ME3AhgFMAVZnYvyUsAbDGz60m+AsCX0Elg9hskP2xmJ6GT6OjTJGfR8Ya4dNCBMaDBsWgFBphevwpRf6RVIepP8TpNcifc4Gz7QNfnO9Bxt3Dr/QuA/61QY6DBsWgDBj2NEqIJSKtC1J8W6DQ4II/kKMnvkvxqsn4CyduTbCbXkJwoz0whBqQlQT7SqWg8LdCqdCoaT+Q67efJ8e8DuB/A4cn6xwBcZmZXJ7MdnAfgUwvt4IjxpXbUoiNyGSqEjycn92DP1P6F05a24C63i4F1unrxYnve8vpnMBXN4cd792HHwYO90wu3R6sD6xQAVowvsecuPrxXMSGCeezg09g9daD1fWrQ4JjkOgC/BuCjAP4omd/vVHTmlAM66SM/hB5iPmrREfjLF5+b11YhMvyne64MKGWl3LmS3ATgv6ITQPA3Znap8/0foTPv5TSAnwL4nSTgoxSK0unzli/Drb/2a2WZKVrI6/7pnwJLFq/VAJ2+FsAnALwEwNlzUfBlUZROAeC5iw/HZ1/2WyVZKtrIu+76fECpcvrUOhHqVvEJdNIazh2N1QB2d81JuA2d+f8ykLxgbtLnPVP7fUWEKJ+CEwswLBf8dwFsNLOXALgOnTSjZVKITrcfDM0mLkQJFKjVQJ3+GMC56My3WwW5dQqktbpr6kC5lgoxH5En6+k5OCb56+jMW3tnngbMbPPcpM9HjC/NswshBsPQiaztZ+lNSC74W81s7o7w2/BE2hZFkTpds3hRwdYJEUjxWg3R6cNm9j38bLBaGoPqFEhrdeX4kgKtEyKQcvrUWhHiVvFqAG8i+avozP93ODqvqFaQHEvudvvKZiKKpbcjXzzku/8s5c6131zw5wH4WtFGdCGdNhBGmqM0f19YuFb71WnZSKciApr5NLgfel6azey9ZrbOzI4HcDaAfzaz3wJwK4C3JMXOQVPygYt2MsQ88CR/G8BGAB8v6ue4SKciGoao1bKRTkU0RO5WMcg8xxcBuJrkR9Dxrby8GJOEKBhDnuCB7T1S0gblgid5OoD3A/gVMxuGM690KppD8VoN0mkNkE5Fc8iLine truncated
"text/plain": [
"<Figure size 720x288 with 6 Axes>"
]
],
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAEXCAYAAACwK7IRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO29e7QlVX3v+/3uR+9+0C9oxBY4gtrRYM4I3LRoNBoOzzbRYM41CvF68Rw86B3hxjz0iMmIRBJzUJMQx/GVVhDMURHx1ceBkgbhoDcJ9kY5CCjStCDd0mB3Qz/o5977d/+o2rhq1tx71apVtVbNWd9Pjxq9ZtWcs+Zau741Z836/eaPZgYhhBBCCCFEwsiwGyCEEEIIIUST0ABZCCGEEEKIDjRAFkIIIYQQogMNkIUQQgghhOhAA2QhhBBCCCE6GBt2A4Sog/POO9127tzdU5m77vrxzWa2rqYmCSHmQHoVIizaoFkNkEWU7Ny5G3d+9x97KjM2+h9W1dQcIcQ8SK9ChEUbNKsBsogTAzAzM+xWCCGKIL0KERYt0KwGyCJSLHrxChEP0qsQYRG/ZuWkJ+JlZqa3rQAk15F8gORmkpd5jr+K5PdITpF8vXNsmuTd6bahom8pRBzUoFchRI1ErlnNIIs4MQAVh1EnOQrgowDOAbAVwCaSG8zs/o5sPwXwFgDv9FRxwMxOrbRRQsRADXoVQtRICzSrAbKIlFpe/5wOYLOZbQEAktcDOB/AMwNkM3s4PRbe47IQQyP+17VCxEX8mtUAWcRL7+JdRXKyI73ezNZ3pI8H8GhHeiuAl/ZQ/8K0/ikAV5rZV3ttoBDREnlnK0R0RK5ZDZBFnJTzsN1hZmtraM0szzWzbSSfB+BbJH9gZg/VeD4hwqAFHvFCREULNKsBsoiUWl7/bANwYkf6hHRfsRaZbUv/30LydgCnAdAAWYgWvK4VIi7i16wGyCJO6nm63QRgDcmTkQyMLwDw+0UKklwJYL+ZHSK5CsArAHyw6gYKESQtmI0SIipaoFkNkEWkGGjVitfMpkheCuBmAKMArjGz+0heAWDSzDaQfAmArwBYCeC1JN9nZi8G8MsA/jF13htBYoN8/xynEqJlVK9XIUSdxK9ZDZBFvNTwdGtmNwG4ydn33o7Pm5CYXrjl/gXAv6+8QULEQuSzUUJER+Sa1QBZxIkBmIl7jUYhokF6FSIsWqBZDZBFpMTvQCBEPEivQoRF/JrVAFnESQscCISIBulViLBogWZHht0AIWrDZnrbhBDDowa9klxH8gGSm0le5jn+KpLfIzlF8vXOsWmSd6fbhoq+pRDxEHkfqxlkESnxv/4RIh6q1yvJUQAfBXAOkqiXm0hucFaP+SmAtwB4p6eKA2Z2aqWNEiIa4u9jNUAWcdICBwIhoqEevZ4OYLOZbQEAktcDOB/AMwNkM3s4PRZ3Ty9E1bSgj9UAWURK/E+3QsRDKb2uIjnZkV5vZus70scDeLQjvRXAS3uof2Fa/xSSdcu/2msDhYiX+PtYDZBFvEQuXiGione97jCztXU0JeW5ZraN5PMAfIvkD8xMoeGFmCXyPlYDZBEnBjBy8QoRDfXodRuAEzvSJ6T7ijXJbFv6/xaStwM4DYAGyEIArehjtYqFiBQDrMdNCDEkatHrJgBrSJ5McgGACwAUWo2C5EqSE+nnVQBegQ7bZSFE/H2sZpBFvET+dCtEVFSsVzObInkpgJsBjAK4xszuI3kFgEkz20DyJQC+AmAlgNeSfJ+ZvRjALwP4x9R5bwSJDbIGyEJ0EnkfqwGyiJMWLGIuRDTUpFczuwnATc6+93Z83oTE9MIt9y8A/n3lDRIiFlrQx8rEQgghhBBCiA40gywixaJfo1GIeJBehQiL+DWrGWQRJ7Ovf3rZhBDDQXoVIixq0GyfoeEvIvlgul1UxVds7QCZ5Deq+hHnOccZJLfWWPcMyX0k19Vxjh7b8y2SB0l+Z9hteQZ1uNEgvVbenmtJHqjr+5ZCeo0Kabby9jxE8jDJ/zHstjxDhZrtCA3/agCnALiQ5ClOttnQ8J9zyh4N4HIkgYBOB3A5yZX9fr2oBsgkHyZ5dpG8ZvZqM7uu7jYVheRrSW5P/9Cz+84nuY3k8jmK/czMjjKzb6b5V5PcQPJnJI3kSc45/jZ9utpL8kck/++Sbb0mrf8Fs/vM7EwAby9TXz3EvwRN6EivXfX6QZKPktxD8hGSf9ZjG48l+TmSu0k+SfKzs8fM7C1IOqKGIL2GgDQ7v2Y76j2a5M97mTAi+RaS0+mAfHY7Y/a4mT0fwN8Ura9+KtfsM6HhzewwgNnQ8L84o9nDZnYPAHe0fR6AjWa2y8yeBLARQN8PNVENkEPGzP4ngG8BuAoASK4A8HEA/4+Z7S5YzQyAbwL4P+c4/jSA1wJYDuAiAB8m+fJe2knyNwA8v5cyQ0GvbEWNDEivVwN4kZktA/ByAG8i+R97aOaXAWwH8O8APAvA3/ZQdrBIr6JmBqTZWT4A4Iclmvmv6YB8dru9RB2DoZxmV5Gc7Ngu6ajRFxr++IKt6afsnEQ7QE6fxr6Tzpo+SfInJF/dcfx2km8tUM8bSU46+/6Y5Ib080R6jp+SfJzkJ0guKtnsPwTwapLnIRHx/zKzQgvbA4CZPW5mH0OyQL7v+OVm9iMzmzGzOwF8G8CvF62f5BiA/w7g/y1aZqjMWG+bGBrSq/f4A2b2dMeuGQAv8OV1IXkukihy7zKz3WZ2xMy+X7RtQ0F6DQpp1k866fQrAD5dso3h0Ltmd5jZ2o5t/bC/wnxEO0BOeSmABwCsAvBBAFeTZI91/E8ALyS5pmPf7+MXNjBXAvglAKci6byOB/BelMDMdgB4B4DPAngNEjHXQnqDeQmA+3oo9scA7khfcTQc04xUeEivDiQvI7kPyYzIEji2d/PwMiS/5XUkd5LcRPI3q25fdUivgSLNdsDEjvYjAC5FMsfaK6eR3EHyxyT/Ip2UaiiVa7af0PB9hZWfi9gHyI+Y2SfNbBrAdQBWAziulwrMbD+ArwG4EABSEb8IwIb0RnAJgD9ObV/2IrERuqCPNv8bEhOIfzazn/dRTzc+AeB/I4ky1RWSJwJ4G0remAaOXtmGiPTqYGZXAlgK4P8A8E8Air4KPgHAuQBuA/BsAH8H4GtMwiY3D+k1VKTZLH8I4E4zu6tE2TuQzDw/C4kJx4UA3lVh26qles2WDg2PZBxzLpMQ8SuR3PsKjW3mI/YB8vbZD6kIAeCoEvV8Dql4kTzZfjWt71gAiwHcRfIpkk8hsU86tnyTsR7AZwD8FsnC5g+9QPJDSIT4BrPC3i7/AOCKHmy1ho9e2YaG9OrBEr4P4ACA9xUsdgDAw2Z2dWpecT0SG71X1NHGSpBeQ0SaTSH5HCQD5D8vUz51TvtJagL5AwBXAHh9t3JDpULNmtkUkpn3m5HYb98wGxqe5O8AAMmXMFm15PeQhIK/Ly27C8BfIRlkb0IyVtnV79dr8PR9o9gI4FiSpyIR8R+n+3cg6YhebGZ9T+eTvBjJa4LXAJgE8CmSp6UenZVA8n1IvNd/08z29FD0LAC/QfKDHfv+leQ7zKzoa98BYoBplqmlRKNXhzEUd5C9B4lDbicNHlVKry0nBs2ejmQG/f7UymQRgEUktwM4Pp1l7wUD0Ku5ygCpXrNlQ8Onx64BcE2V7Yl9BrkSzOwIgC8C+BCAo5GIGWY2A+CTAK4i+SwAIHl86gDQE+nT54cA/BczO4TEBGInenwaJbkQwESanEjTs8feg+Tp/Gwz29ljE38JwK8isQM7Nd33WgBf6bGewWDQjFRLiUGvJEdIvi19ZUiSpwP4AwC3Fqz6KwBWMlk8f5TJovonAPj/emnfwJBeW00MmgXwDQAn4Rd95HsBfB/AqUUGxyRfTfK49POLAPwFEtOTZtICzWqADIDkK1NHmPn4HICzAXwxfRUwy7sBbAbwbyT3ALgFwAvnOM83OPdaph8DcL2ZfRtIXqsC+C8A/ojki4t/GxwAMPtdfpSmZ/kbJEs+beYv1lkstLaqmT1hZttnt3T3DjM7MG/BYSKbxihpkV5/F8BDAPYC+B9IVpD570UqTV8v/g6AdyKxW74MwPmpk1IzkV6jpQ2aNbNDTh+5G8CRjv6yG2cBuIfk00hmUb+MRq177CFyzbK4CapoEiRfhcRW5xCAN5pZ3wbpfbZnIxLP+e+a2VnDbAsArP2l423Tx97WU5mRcy6/y8zW1tQk0WIaqNerkdjxPWFmhZaOqxPpVTSNBmr2ASQreNxgZv95mG0B2qFZ2SAHipndgcTGqRGY2TnDbkMWC/KJVcRJA/V6MYCLh92OXyC9imbRQM16Z82HR/ya1QBZxEuANk9CtBbpVYiwiFyzhW2QU0eP75P8epo+meSdJDeT/EK6bp0QzcCQeNj2skWE9CqCouV6BaRZERgt0GwvM8jvQLI23bI0/QEAV5nZ9SQ/geR13cfnq2D5+GJ71oLlpRoqxCyb92/fYWZd1sEM02u2QvrW67KxxXbshPQq+mOL9FoU9bGiEaiPTSg0QCZ5AoDfBvB+AH/CZJG/M5EsGQYkEXT+El3E+6wFy/HhX/lPpRsrBAD89nf/2yOFMkYu3rmoSq/HTizHB3+5QWaqIkhef9f7pdcuqI8VTUJ9bELRGeR/APBfkYQ8BYBjADzVsRTLViTelTlIXoIkVCSOXbDMl0WI6pkNg9lOKtHrKulVDIp26xVQHytCowWa7TpAJvkaJEsB3UXyjF5PYGbrkYR2xJolqzOPGzPWPUhMkZ/fKqrHfRiyAkFs3DJFvtMgn7mKhOEZYfcWjTgV0fkW7nFvHZ59dM7ty1OkfV4if7r1UaVen+/otVD5EpopQlW6KnLuur5DVRTRmqvPInW4u8rcF4qce05aqFeg3X2s91xd9FekjiLfuy7K9lfd+tgiZdTHVksRJ71XAPgdkg8DuB7Line truncated
},
"metadata": {
"needs_background": "light",
"tags": []
},
"output_type": "display_data"
}
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 316
},
"id": "FdL_orrX_trr",
"outputId": "b716256c-9b8a-40d9-fc01-d75e261a728c"
}
},
{
"cell_type": "code",
"execution_count": 9,
"source": [
"fig, axes = pylab.subplots(1, 3, figsize=(10, 4))\n",
"for i in range(INFLOW.shape.get_size('inflow_loc')-1):\n",
@@ -494,47 +471,47 @@
" axes[i].set_title(f\"Ini. vel. Y {INFLOW_LOCATION.numpy('inflow_loc,vector')[i]}\")\n",
" pylab.colorbar(im,ax=axes[i])\n",
"pylab.tight_layout()\n"
]
],
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": [
"<Figure size 720x288 with 6 Axes>"
],
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAscAAAEYCAYAAACwdltJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO29e7glVXnu+77r1legbwgtjUC0DRGPabVF8xgNctE2FzHnKKKJwgkcwj7bnev2gHoeL0T3wZgE95O4jW0goNEAwahsg7KBwFZ3ItIgUS5BWgRtrvYV+rZ6Xb7zx6yFs0aNteaYNatqVo16f/3U07NqjlHjm7XqrTGq6vvGRzODEEIIIYQQAhgZtgFCCCGEEELUBQ2OhRBCCCGESNDgWAghhBBCiAQNjoUQQgghhEjQ4FgIIYQQQogEDY6FEEIIIYRI0OBYCCGEEEIMDZKbSD5AcivJiz3f/xHJ+0h+j+QtJI/r+u4ckg8myzmF2KN5joUQQgghxDAgOQrgBwDOALANwB0A3m5m93WVeR2A281sP8n/AOAUM3sbyVUAtgDYCMAA3Ang5Wa2axCb9ORYCCGEEEIMi5MBbDWzh8zsEICrAZzZXcDMbjWz/cnqtwGsSz6/AcBNZrYzGRDfBGDToAaNDboDIerOG95wsu3YsaevOnfe+YMbzWxggQkhwpFWhag/Jej0GAA/6VrfBuCVC+zuPABfW6DuMX0Z50GDYxE9O3bswe3f+XRfdcZGX7emJHOEEPMgrQpRf3Lq9ESSW7o2bTazzf22TfK30XGh+JV+6/aDBscifgzA7OywrRBC9EJaFaL+5NPpdjPbOM93jwI4tmt9XbItBcnTAbwfwK+Y2WRX3VOcurf1a5yLBseiBZg6XCEagbQqRP0pXKd3AFhP8gR0BrtnA3hHdwGSLwXwaQCbzOyprq9uBPBfSK5M1l8P4L2DGqTBsWgH6nCFaAbSqhD1p0Cdmtk0yXejM9AdBXCFmd1L8hIAW8zsegAfB7AcwD+QBIAfm9mbzGwnyT9BZ4ANAJeY2c5BbdLgWMSPAdCUhULUH2lViPpTgk7N7AYANzjbPtD1+fQF6l4B4Ioi7dHgWLQAvaoVohlIq0LUn/h1qsGxaAeRC1mIaJBWhag/ketUg2MRP4qAF6IZSKtC1J8W6FSDY9EC4n8FJEQcSKtC1J/4darBsYifFtzlChEF0qoQ9acFOtXgWLQAAy1uIQsRB9KqEPUnfp1qcCzaQeR3uUJEg7QqRP2JXKcaHIv4MQCzmjtViNojrQpRf1qgUw2ORQuIP3hAiDiQVoWoP/HrVINjET8tCB4QIgqkVSHqTwt0qsGxaAeRBw8IEQ3SqhD1J3KdanAsWkD8r4CEiANpVYj6E79ONTgW8dOC4AEhokBaFaL+tECnGhyLFhD/Xa4QcSCtClF/4tepBseiHUQuZCGiQVoVov5ErlMNjkX8GMDIhSxEFEirQtSfFuh0ZNgGDBuSXyN5TsltnEJyWwn7PZ6kkdxL8oKi95/DnitJHijjtw6GAdbnImqFdFqoPR8muS+xqWYPSKTVpiOtFmrPP5M8SPJbw7YlTfw6jXJwTPJhkqeHlDWzN5rZVWXbFALJo0huJ3mKs/0KklcvUHWFmW1Oyk6QvC45BubZ13tI3kPyGZI/IvmenLZ+INn/s8fZzM4F8MY8+yud2dn+FlE60umCOv1Dkg+RfJrkYyQv62cgS3Ipyf+W2LmH5DfmvjOzDwI4KXRflSOt1g5pdX6tdu1zguT9/Qzak0H+bDIYn1uevbEws1MBXBi6v0qJXKdRDo6bipk9CeAPAXyG5BIAIHkagF8H8J/62NW3APw2gCc83xHAuwCsBLAJwLtJnt2PnSSfD+CtAB7vp97QmJuwPGIhi+qoSKfXA3iZmR0O4MUAfhHA7/Wx780AVgH4heT/P+yj7vCQVkWBVKTVOd4D4Kc5zHzMzJZ3LbW4sViQFug0+sExyXNJfovkn5HclTwtfWPX97eRPD9gP28jucXZ9ockr08+L0ra+DHJJ0n+9ZwY+8HMPgfgAQCXJPU/DeD3zCxIdGZ2yMw+YWbfAjDj+f5PzewuM5s2swcAfAXAq/s085MALgJwqM96Q8I60870s4hKkU4z3//QzHbP/QQAswBeELJvkicCeBOAC8zsp2Y2Y2Z3htQdPtJq3ZFWvb/lBHQGz/9fv/Y1k/h1Gv3gOOGV6IhjDYA/BXA5Sfa5j/8O4OdJru/a9g4AX0g+XwrghQA2oNOJHQPgAzntvRDA7wC4GsA9ZrbQ65/cJMfgNQDu7aPOWwFMmtkNZdhUCi24y40E6bQLku8g+TSA7eg8Of50YNWTATwC4MPJK+Xvk/w/irStNKTVpiCtpvlLAO8DcCBH3eckg/8fJe5Tywq2rXhK0CnJTSQfILmV5MWe719L8i6S0yTf4nw3Q/LuZLm+iJ/YlsHxI2b2GTObAXAVgLUAjupnB2a2H52nrG8HgETQJwK4PrkoXADgD81sp5k9A+C/AOjLXaGrrW3oXAROB/Af8uwjkA+hcw78bUhhkoeh87t+v0SbykEdbhOQTtP7/0LiVvFCAH8N4MnAquvQccXYA+C5AN4N4CqSv1C0jaUgrTYBaTWB5G8CGDWzL+Wo/u/oDP7XAjgVwMsB/EWB5pVHgTolOYrOG+k3AngRgLeTfJFT7McAzsXPbp66OWBmG5LlTYP/uPYMjp/1E0oECQDLc+znC0iEjM4d7peT/R0JYCmAO0nuJrkbwNeT7Xm5F8AuMyvFr5fku9HxPf41M5sMrPYhAJ8zs4fLsKk84o+sjQTp1IOZPZi0898CqxwAMAXgI8kr4f8J4FYAry/JxAKRVhuCtAogecr7p+gvHuBZzOwJM7vPzGbN7EcA/h8ADXjLU7hOTwaw1cweMrND6DzhPzPVotnDZvY9dFzMSqdm0/jUnpsAHElyAzqCngty2Y5Oh3SSmT06LONCIfk7AC4G8NrkjjqU0wCsI/l/J+tHAriW5MfM7GNF21kYc6+ARFuIQqcOYwCeH1j2e55tzRhFSqtto+laXQ/geADfTLxKJgAcQfIJAK/K8SDJ0ISHlsXr9BgAP+la34aO604oixP/9WkAl5rZlwc1qP5/hBphZlMA/gHAx9GJAL8p2T4L4DMALiP5HAAgeQzJNwzDziSQYXGyOkFy8Zw/GMnfQuf11Blm9lCfuz4Nnde1G5LlMQC/i87rkHoTefCA+BmR6PT8LhtfBOC9AG4J3PU30HkF+V6SYyRfDeB1AG4s9heUhLTaGiLQ6j0AjsXP+sTz0XF/2oD0YG++/b6O5HHscCw6ftZfKeVHFE3/Ol1DckvXUuQ80seZ2UZ03j58gp0ZtQZCg+MuSL6G5N4exb6Ajt/SP5jZdNf2iwBsBfDtJIjmZgA/P087XyP5viJsnocH0LnrPgadDvEAgOOS7z4CYDWAO/izeRX/OmSnZrYjeQ30hJk9gU7k7i4z63XMhozJjzEiWqLTVwP4Psl9AG5IliBbkgHHmQB+FR2/488AeJeZ/Xuh1peCtBoTsWs1mfWpu0/cCWA2WffObOHwUgD/AmBf8v/3kdNFo1py6XS7mW3sWjZ37fBRdG4y5liXbAuzJnm7kDzwuw2d4zoQNPlsNRaSx6Ej2oMA3mNmnxmyPZejM//xU2YWNO1UFWw8cZ3dcXk/U1oCI7988Z3JnagQA1FDnX4QwB8BWARgWWAnXgnSqhgmNdTqTQBeBeA7ZnbaMG3ppmidspPg6AfovJ1+FMAdAN5hZpmZtEheCeCrZnZdsr4SwH4zmyS5BsC/AjjTzO7ry0AH+Rw3GDN7BMDingUrwszOA3DesO3wotevYkjUUKcfBvDhYdsxL9KqGBI11OoZw7ZhXgrUqZlNJ5ME3AhgFMAVZnYvyUsAbDGz60m+AsCX0Elg9hskP2xmJ6GT6OjTJGfR8Ya4dNCBMaDBsWgFBphevwpRf6RVIepP8TpNcifc4Gz7QNfnO9Bxt3Dr/QuA/61QY6DBsWgDBj2NEqIJSKtC1J8W6DQ4II/kKMnvkvxqsn4CyduTbCbXkJwoz0whBqQlQT7SqWg8LdCqdCoaT+Q67efJ8e8DuB/A4cn6xwBcZmZXJ7MdnAfgUwvt4IjxpXbUoiNyGSqEjycn92DP1P6F05a24C63i4F1unrxYnve8vpnMBXN4cd792HHwYO90wu3R6sD6xQAVowvsecuPrxXMSGCeezg09g9daD1fWrQ4JjkOgC/BuCjAP4omd/vVHTmlAM66SM/hB5iPmrREfjLF5+b11YhMvyne64MKGWl3LmS3ATgv6ITQPA3Znap8/0foTPv5TSAnwL4nSTgoxSK0unzli/Drb/2a2WZKVrI6/7pnwJLFq/VAJ2+FsAnALwEwNlzUfBlUZROAeC5iw/HZ1/2WyVZKtrIu+76fECpcvrUOhHqVvEJdNIazh2N1QB2d81JuA2d+f8ykLxgbtLnPVP7fUWEKJ+CEwswLBf8dwFsNLOXALgOnTSjZVKITrcfDM0mLkQJFKjVQJ3+GMC56My3WwW5dQqktbpr6kC5lgoxH5En6+k5OCb56+jMW3tnngbMbPPcpM9HjC/NswshBsPQiaztZ+lNSC74W81s7o7w2/BE2hZFkTpds3hRwdYJEUjxWg3R6cNm9j38bLBaGoPqFEhrdeX4kgKtEyKQcvrUWhHiVvFqAG8i+avozP93ODqvqFaQHEvudvvKZiKKpbcjXzzku/8s5c6131zw5wH4WtFGdCGdNhBGmqM0f19YuFb71WnZSKciApr5NLgfel6azey9ZrbOzI4HcDaAfzaz3wJwK4C3JMXOQVPygYt2MsQ88CR/G8BGAB8v6ue4SKciGoao1bKRTkU0RO5WMcg8xxcBuJrkR9Dxrby8GJOEKBhDnuCB7T1S0gblgid5OoD3A/gVMxuGM690KppD8VoN0mkNkE5Fc8iLine truncated
},
"metadata": {
"needs_background": "light",
"tags": []
}
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 318
},
"id": "mLhwxsZd_trs",
"outputId": "5bb43835-6e21-4d97-bb07-17b17215bbc8"
}
},
{
"cell_type": "markdown",
"metadata": {
"id": "3mitRMPx_trt"
},
"source": [
"## Re-simulation\n",
"\n",
"We can also visualize how the full simulation over the course of 20 steps turns out, given the new initial velocity conditions for each of the inflow locations. This is what happened internally at optimization time for every gradient calculation, and what was measured by our loss function. Hence, it's good to get an intuition for which solutions the optimization has found.\n",
"\n",
"Below, we re-run the forward simulation with the new initial conditions from `initial_velocity`:\n"
]
],
"metadata": {
"id": "3mitRMPx_trt"
}
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 263
},
"id": "OINe99pf_trt",
"outputId": "a7c9ade3-3685-40e6-a5f1-97a5084deec9"
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAr4AAADiCAYAAACoVMTQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxcVfk/8M8nadJ93yhtaSkUS2UptLKJiALKoqIIyirgUvSrX0QRvxUBRQFRqbgCwg9olU1WQUSwAgVLodACXeje0tI93dKmafY8vz/uDeTcM+1MJjPJ3NzP+/XKqzknz733ZPrM5Myds9DMICIiIiLS0RW1dwNERERERNqCOr4iIiIikgjq+IqIiIhIIqjjKyIiIiKJoI6viIiIiCSCOr4iIiIikgjq+IqIiEiskBxM8mWSFSQnk/wpyfvyeL0TSTaS3EXy1HxdpwXteYFkNckZbXS9G0huIbmxLa6XT+r47oWeWG37xJLcUN5yCskqkmvbuy2SOeUtV5CszefvXOhIriJ5cobhEwFsAdDLzK7MY7OaW29mPczsWQAgOYTkUyTXkzSSI6MHkDyZ5JskK0muJfmlTC4Uyc+mr4ubfm5mnwTwzUwbHj6f6sLzlJOcSfLYDI/dD8CVAMaa2T6ZXrNQJa7j29GeWCRvIbks/GOxmORXsrkoyXvC8x/YVNfSJ5bkTwfM21+RXENyJ8nVJK9uycVIDiT5AMkdJLeTvL/pZ2Z2CYDTWv8rSWt1tLxtQrIfyc0tuSlA8hKSDZGOzIlNPzezAwDc1LpfJ1FGAFho7bsLVyOAZwF8MdUPSY4F8ACAHwPoDeBwAHNacP6m/Gz6mtrK9v7NzHoAGADgRQCPZHjcfgC2mllZSy/IQEH1NQuqMQWo4J9YACoBfBbBk+piAL8jeVxLLkDyeAAHtKaRUlDikLd3AxhjZr0AHAfgApJnteD8jwPYiOAFeRCAW1rRVikMccjbJr8EsCiL878a6chMz+IciRC+UZgR3tzZTvJdkqeFP5uC4O/dD8M3EN6bK5KfI/lOeHdzOsmDw/pLSf6jWdwyko80K68hOS6TNprZJjO7DcAbewi5BsCfzexfZlZvZlvNbEWmj0G+mFk9gPsBDCU5EABI9iZ5N8kNJNcxGNpQHD620wDsGz7WU8L4Y8K7xuUk5zZ/Exc+3jeSfAXAbgCjSI4hOY3kNpJLmt/5ZvAp3Z9I/jO8iTeL5AHNfv7hZsduarpRQrKI5CQGn5ZsJfkwyX7pfv9Ed3w7whPLzH5iZovNrNHMZgH4L4CMPr4I29IJwB8A/G+mx0j76iB5u8TMKptVNQI4MFVsivZ/CsBwAFeZ2Q4zqzOztzI5VtpPR8jb8HzHATgEwL2Z/ebSCkcDWILgDuWvANxNkuGnOvcD+FX4BuI/zQ8ieRCABwFcAWAggGcA/INkKYCXAHws7DTtC6AU4d9MkqMA9AAwL0ftPyY87/ywQ3lfJh2zZgaFHb13Sd5KsnsuGhU+Dl8BsBXA9rB6CoB6BK/DRwD4FICvh4/tafjg7vMlJIcC+CeAGwD0A/ADAI81daJDFyH4FKcngM0IOs8PILhRcS6A2xjcEW9yLoDrAfQFsBzAjWFbewL4D4I3pPuG7Xs+POZ/AXwewMfDn20H8Kd0v3+iO76huD+xmrepK4CPAHinBYd9D8DLZpbz9khexT5vw3fquwCsBdAdwYtiJo5B8LtPDd/lv0Hy47lql+RVrPOWZDGAPwL4DoBs7kwfwWCC0FKS14Y3HmTPVpvZXWbWAGAqgCEABmdw3JcB/NPMpplZHYJPhLoCOM7MVgKoADAOwAkAngOwnuQYBB2o/5pZY47aPwxBB/CLAEaHbfhDhscuDts4BMAnAYwH8JtWtudLJMsBVAH4BoCzzaye5GAApwO4wswqwyENtyLojKZyIYBnzOyZ8KbbNACzw3M0mWJm74R3l08FsMrM7g3vfL8F4DEA5zSLf8LMXm92N7rpzepnAGw0s8lmVm1mFeFNPiAYivljM1trZjUAfgrg7HTPK3V84//Eau4OAHPD66VFcjiAywBcl4e2SH7FPm/N7GYEdwOOBPBXADsyPHQYgrsRLwLYB8BkAE+SHJCrtknexD1vLwcwy8xaMk6zycsI7hQPQtAROg/AVTlqV0f1/goCZrY7/LZHBsftC2B1s2MbAawBMDSsegnAiQjy5SUA0xHkysfDcq5UAbjXzJaa2S4EY7hPT3NMU5s3mtnCsGP5LoAfIv0QnHQeNrM+CJ5zCxB0poFgmFEJgA3hJyrlAP6MIFdTGQHgnKbYMP54BM/nJmsi8UdH4i9A8PrdpPlqEbvxwf/zcAB7Gh4yAsATzc65CEAD0rym6N1m5IlFEsjyiUUy1RPrwPD7cgRPqmOR2ycWAIDkrxG8qH6iBWPkfgvgZ2aWaYdDCkeHyNswV98i+WkEH3N9P4PDqhDcPbg7LD9E8scAPgrgyVy3UXIqtnkb3k2+HB90Flok7KA3mU/yZwg6vr/IQfPEtR7AoU0FBok2HMC6sOolBHNj9kfQGW3qiB2L4I5+rsyD+8lAa8avG3J0s9LMtpCcCGA2yQcQdFJrAAwI77imswbAX83sG3u7TCT+JTM7JYvmrsGe7zyvAfBVM3ulJSfUHd/srUfwbgPAHp9YJwL4WPj9S8jPO0qQvB7BGJxPmdnOFhx6EoBfk9zID9bme5Xk+blsnxSUgsnbiE7IfIJl9I8JUpSlYymEvD0KwR2theHr5e8AHBW+fhZncT4DwBy1TVwPAziD5EkkSxAsxVUDYGb485cAfAJAVzNbi2BuzKkA+gNo0XwBkl0AdA6LncNyk3sBXEpyFMluACYBeDrD836C5AgGhgO4GWne2DNYReWSTM5vZksQfDryQzPbAODfACaT7BUOGzpgL0PI7gPwWZKfZjABrguD5deG7SH+aQAHkbyIZEn49RGG4/TTeBrAEJJXkOxMsifJo8Of3QHgRpIjwt9/IMkz051QHd/sFcQTi+SPAJwP4GQz29rC3+EgBMurjMMH42k+C+CJFp5H4qPd8zZ8Ub2MZN/wRf0oAN/GBxMW0nkCQF+SF4cvumcjGP7Qonf9EivtnrcA/gVgJD54vbwuPPe4cOhGuvOeFo6lRDgM41roE4q8CDt1FyIYT7sFwd+1z5pZbfjzpQB2IcgThDeMVgJ4JZP/y4iq8FxAMC63qlk77gHwFwCzEHxiUYPgU4NMHIEgvyvDf+fv7dhwvHt/AK+1oO2/BjCR5CAEk91KASxEMEnsUbhDF95nZmsAnAngagQT19Yg+PQiZZ/SzCoQDE87F8Gb2I0IVkbpnCo+xbGnIPg/3AhgGYLnOhC8+XwKwL9JViD43Y9OdZ7mNNQhS2a2hGTTE2sogLcReWIxmLjz/hOL5EoAm7N8YjVZHP7bdKfgJgC1AJaHHxsCwE1mlnY9SIusyRcev8XMqlIfIXFXQHn7BQQf8ZYieCH8AzKc9GFm20h+DsBtCGbwLgZwppltaWH7JCYKIW/DyTPvD9UguQNAnZllupPVSQCmkOwBYBOCu2Zat7cZMxvZ7PspCFYaaP5zNvv+ksjPfhopP4G93MQxsyGR8oQ0zWsE0D8cS/plM3su2qY9XOcnAH6S5typjvsN9jKZjeQ0BBN9Xw+rjgLwtJktThUffXzCulkIJhY3+Vb4FY2bjuDmQvTYlHeEzezEFHVLAJyxh/hL9nY9M1uA4PkTPa4RwWPUokl/zHw4qOQbyaaJGTVo9sRqx/a8/8QyMy/pRICCzNu7EcwWLjOzjJZIk+QpwLxdgqBT/7CZfbU92yLSkanjKyIiIiKJoDG+IiIiIpII6viKiIiISCKo4ysiIiIiidCmqzqQ1IBiyVq62bP5oryV1lDeShwpbyWmtpjZwL0FtMNyZtms8y3S0hWJck15K9lQ3kocKW8lrhpWp4vQUAcRERERSQR1fEVEREQkEdTxFREREZFEUMdXRERERBJBHV8RERERSQR1fEVEREQkEdTxFREREZFEUMdXRERERBJBHV8RERERSQR1fEVEREQkEdTxFREREZFEUMdXRERERBJBHV8RERERSQR1fEVEREQkEdJ2fEl2Ifk6ybkk3yF5fVg/heS7JN8Ov8blv7kimVHeShwpbyWOlLcSJ50yiKkB8Ekz20WyBMAMkv8Kf3aVmT2av+aJZE15K3GkvJU4Ut5KbKTt+JqZAdgVFkvCL8tno0RaS3krcaS8lThS3kqcZDTGl2QxybcBlAGYZmazwh/dSHIeyVtJds5bK0WyoLyVOFLeShwpbyUuMur4mlmDmY0DMAzAUSQPAfAjAGMAfARAPwD/l+pYkhNJziY5O0dtFsmI8lbiSHkrcaS8lbhg8AlFCw4grwOw28xuaVZ3IoAfmNln0hxrQHE27ZTEa4CZMdujlbfSPpS3EkfKW4mrhjlmNmFvEZms6jCQZJ/w+64ATgGwmOSQsI4APg9gQesbLJIbyluJI+WtxJHyVuIkk1UdhgCYSrIYQUf5YTN7muQLJAcCIIC3AXwzj+0UaSnlrcSR8lbiSHkrsdHioQ6tupg+wpCste6jt9ZQ3kr2lLcSR8pbiascDHUQEREREekI1PEVERERkURQx1dEREREEkEdXxERERFJBHV8RURERCQR1PEVERERkURQx1dEREREEkEdXxERERFJBHV8RURERCQR1PEVERERkURQx1dEREREEqFTezdARPKvtGSQV7dP93FOuRENXszO2rVOuWL3Ci/GUN/K1om0BjOIsby3QiT3ijOI8V+3Ze90x1dEREREEkEdXxERERFJhLQdX5JdSL5Oci7Jd0heH9bvT3IWyeUk/0ayNP/NFcmM8lbiSHkrcaXclbig2d7HPpEkgO5mtotkCYAZAL4L4PsAHjezh0jeAWCumd2e5lyW2ZiVQuKOHysq6poixn3/YI3VXoTGQbZWA8wsk8FLine truncated
"text/plain": [
"<Figure size 720x432 with 4 Axes>"
]
},
"metadata": {
"needs_background": "light",
"tags": []
},
"output_type": "display_data"
}
],
"source": [
"smoke = initial_smoke \n",
"velocity = initial_velocity\n",
@@ -547,45 +524,45 @@
" axes[i].imshow(smoke_final.values.numpy('inflow_loc,y,x')[i,...], origin='lower', cmap='magma')\n",
" axes[i].set_title(f\"Inflow {INFLOW_LOCATION.numpy('inflow_loc,vector')[i]}\" + (\", Reference\" if i==3 else \"\"))\n",
"pylab.tight_layout()\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "ootUIrQq_tru"
},
"source": [
"Naturally, the image on the right is the same (this is the reference), and the other three simulations now exhibit a shift towards the right. As the differences are a bit subtle, let's visualize the difference between the target configuration and the different final states.\n",
"\n",
"The following images contain the difference between the evolved simulated and target density. Hence, dark regions indicate where the target should be, but isn't. The top row shows the original states with the initial velocity being zero, while the bottom row shows the versions after the optimization has tuned the initial velocities. Hence, in each column you can compare before (top) and after (bottom): \n"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 410
},
"id": "7lm9Rrel640-",
"outputId": "6ef7914a-8f63-49af-9a1a-fde1aed37682"
},
],
"outputs": [
{
"output_type": "display_data",
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlEAAAF1CAYAAADIqb9jAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOy9eZycVZX//zm9L+nurJ2dLJAQQpAtRHBhERBQERdEcBlUHB1nHHXG78z4nfl+3b7OjMv35zYuY0YYcRvhy6CgosiqIAiEnSSEhOz70kt6766q8/ujntTznJuq7uru2m7V551XvXJP3We5VX0/9Zzn3vOcK6oKQgghhBAyPqqK3QBCCCGEEB+hE0UIIYQQMgHoRBFCCCGETAA6UYQQQgghE4BOFCGEEELIBKATRQghhBAyAehEjYGIXCgiuyP2ehG5MCiLiPyniHSKyOPBex8RkQMi0isiM0Y57vtEJB5sd0reP8gYiMjLIjIsIj8udltIaVNBmrhfRAZF5OFit4WUNpWgCRGpD9oxIiJfKGZbSomSd6KCTvS8iPSLyH4R+a6ITC1We1T1VFV9MDBfA+BSAAtUdY2I1AL4KoDXq+oUVT0yxuEeDbbbCAAiskpE7haRwyJiEngFHfhGEdkhIj0i8oyIXDGRzyAi94mIikhN5HOdCOBfJnI8UlioiRAR+bGI7BORoyLykoh8cDxtF5GlIvKrQFOHReTLkc/1OgB/MZ7jkeJATRyPiCwLbgKyvjEWkc8GTlJv5LU0+ExDqjoFwE+yPV4lUNJOlIh8EsCXAPwdgDYA5wJYBOAeEanLsE9NuvfzxCIA21W1L7BnA2gAsH6CxxsBcCuAG9LU1QDYBeACJL+L/wXgVhFZPJ4TiMi7AdROsH2kyFATx/GvABaraiuANwP4goicnc2Bg+/rHgD3A5gDYAEAjsR6BjWRkW8DeGICx78lcNqOvbZO4BgVQ8k6USLSCuBzAP5aVX+rqiOquh3ANQAWA3hPsN1nReS24I70KID3icgSEflDcHd5r4h8O1tvXEQaReQHwdDrBgDnOPXbReQSEbkBwPcBnBd46/8FYFOwWZeI3D/ez6yqm1T1RqQRl6r2qepnVXW7qiZU9VcAtgHI6oIRtL0NwGcA/P1420aKDzWRtn69qg4dM4PXiVke/n0A9qrqVwN9Darqc+NtIyke1ETG9l0LoAvAfeM9PhkfJetEAXgVkt767dE3VbUXwF1IDo8e4yoAtwGYiuRQ408BPA5gBoDPAnjvOM77GSR/hE8EcBmA69NtFHTiv0A41HodgFOD6qnBVEDeEJHZAJZjfHcz/wLguwD256VRJN9QE2kQke+ISD+AFwHsQ/K7yIZzAWwXkd8EUyMPishp+WgjyRvUhEPgWH4ewN9O8BBXikiHJOO6PpLDppUlpexEzQRwWFVjaer2BfXHeFRVf6GqCQCzkLwr+LSqDqvqwwDuHMd5rwHwz6raoaq7AHxzgu3PG8Gc+k8A3KyqL2a5z2oArwbwb/lsG8kr1EQaVPUvAbQAeC2SF9Oh0fdIsQDAtUh+nnkAfg3gjkxTQKQkoSaO5/8AuFFVd4+55fHcCuAUJL+fPwfwaRG5LpeNKzdK2Yk6DGBmhrnruUH9MXZFyvMAdKhqf4b6sZjnbL9jHPvmHRGpAvAjAMMAPjqOfb4D4OMZfmyIH1ATGVDVeHAhXAAg27vnAQAPq+pvVHUYwP9FclSi6E8GkqyhJiKIyBkALgHwtYnsr6obVHVvoKdHAHwDwNW5bGO5UcpO1KNI3lG+LfqmiEwBcAXsXG/0CYV9AKaLSFPkvYXjOO8+Z/sTxrFvXhERAXAjkoGJb1fVkSx3bQWwGsAtIrIfYbDhbhF5be5bSvIENTE2Ncg+Juo52O+J+Ac1YbkQyViwncFv/f8A8HYReWqCx1MAkpumlScl60SpajeSAYP/JiKXi0ht8CTarQB2Izkak26/HQDWAfisiNSJyHkArhzHqW8F8D9FZJqILADw15P4GONCkjQAqAvsBhGpj2zyXSTvkq9U1YFxHLobyTunM4LXG4L3zwbw2KQbTgoCNWE1ISLtInKtiEwRkWoRuQzAdcg+mPbHAM4NAoCrAXwCyZGLjbn/JCQfUBPHXSfWInkTcey3/t+RnKa+LMtjXxV8JhGRNQA+BuCOXH+GcqJknSgAUNUvA/hHJIfZjyJ5wd8F4OLIEznpeDeA8wAcAfAFALcgEicRPCWRaQTmc0gOzW4D8DtkEGE2BAGr/ziOXRYhOcVwLFh8AMGTHCKyCMCHkRTGfglzeLx7rINqkv3HXgAOBVUHgmkM4gnURKgJJO+SP4LkxbITye/kE6qaVWyLqm5C8umtfw/2vwrAm6kJv6AmQk2oar/zW98LYFBVD6U/1HFcC2ALgB4APwTwJVW9eRxtqzhEtfxHs0XkFgAvqupnit2WY4jIewF8D8nYpvOOJVIrYns2AZgP4FZV/UAx20LyDzWRVXvuQfIJvsdV9eJitoXkH2pizLbUAziAZJ7BL6vq54rVllKiLJ0oETkHQAeSdwmvB/ALJDvg00VtGCFFgpogxEJNkFxQ0tN5k2AOgAeRHMr8JoCPUBh+E8Q7bBKRLSLyqTT1fysiG0TkOUkua7OoGO0sYaiJMoOamDTURJlRDE2U5UgUKS+CoN+XkEyctxvJpwuvU9UNkW0uAvCYqvZLMkHchar6zqI0mJA8Q00QYimWJsp1JIqUF2sAbFHVrUHQ78+QDAJOoaoPRHK+/AnJfEGElCvUBCGWomiCThTxgfmwie12B+9l4gYAv8lriwgpLtQEIZaiaKKQK1mjtaZJZ9W3FfKUpAAcGurG0Vj/qAnZLrtsjR450p227sknX1oPYDDy1lpVXTuRtojIe5BMLHrBRPYvNG21TTqbmig7Dgx1o3uEmpgIbbVN2l5HTZQbB4fLUxMFdaJm1bfhy6fcUMhTkgLw9xtvHHObI4e78Nij30lbV1N/6aCqrh5l9z2w2YEXBO8ZROQSAP8E4IIx8sOUDLPr2/CtVe8rdjNIjvnoCz8YcxtqIj3tdW345qr3F7sZJMd87IX/HHMbHzVRUCeKVDAKIB6f6N5PAFgmIkuQFMW1AN4V3UBEzkQyn8rlqnpwEi0lpDBQE4RYPNQEnShSIBRIJCa2p2pMRD4K4G4A1QBuUtX1IvJ5AOuCDNVfATAFwP8TEQDYqapvzk3bCckH1AQhFv80QSeKFAYFEItNfHfVuwDc5bz36Uj5kgkfnJBiQE0QYvFQE3SiSGFQnZQ4CCk7qAlCLB5qgk4UKRAK0YkN0xJSnlAThFj80wSdKFIYFEBswgGDhJQf1AQhFg81QSeKFAYPh2kJySvUBCEWDzVBJ4oUCAXifomDkPxCTRBi8U8TdKJIYVAACS52TUgKaoIQi4eaoBNFCoR/w7SE5BdqghCLf5qgE0UKgypk4ploCSk/qAlCLB5qgk4UKQyTTKJGSNlBTRBi8VATdKJIgVDv5roJyS/UBCEW/zRBJ4oUBg/vMAjJK9QEIRYPNUEnihQG1cmszk1I+UFNEGLxUBN0okjh8GyYlpC8Q00QYvFME3SiSGHwMBMtIXmFmiDE4qEm6ESRAuHfMC0h+YWaIMTinyboRJHC4OHCkoTkFWqCEIuHmqATRQqDh+n8Cckr1AQhFg81QSeKFAj17g6DkPxCTRBi8U8TdKJIYfDw0VVC8go1QYjFQ03QiSKFwcO5bkLyCjVBiMVDTdCJIoXDs7luQvIONUGIxTNN0IkihUH9m+smJK9QE4RYPNQEnShSGBRAPFHsVhBSOlAThFg81ASdKFIg/Fudm5D8Qk0QYvFPE1XFbkA5oZCMr4rnWMBgulcWiMjlIrJJRLaIyKfS1J8vIk+JSExErs518yuZppqYec2f0pd6zWu2r2n1w+ZVJZp6KWBeFQ81UZZE+3y6l6sDaiKCh5qgE0UKgyoQS6R/jYGIVAP4NoArAKwEcJ2IrHQ22wngfQB+muOWE5IfqAlCLB5qYkwnSkQaRORxEXlWRNaLyOeC95eIyGOBx3eLiNTlqlGkPNF4Iu0rC9YA2KKqW1V1GMDPAFxljq26XVWfA5D3CXVqguQKaoIQi2+ayGYkagjA61T1dABnALhcRM4F8CUAX1PVkwB0ArghV40iZcixdP7pXmMzH8CuiL07eK9YUBNk8lAThFg81MSYgeWqqgB6A7M2eCmA1wF4V/D+zQA+C+C7uW9ifqlywpVqJfMfq3uk2ti7BuzXd2Qo83lm1Ft7XoOd422u8euxznFzbJg2PTNFZF3EXquqawvQqglRbpporLZ9b8X8Q8aefp7t91Unz02VtbPX1A0/ecDY256flio/dXiaqRtM2Hu46XUxYzdXUxMRm5ooIaqd68Sc5v5UubF2xNTtPtpi7K19janygUGrgaYae9wFjdREhJLURFZP5wVzjU8COAnJOceXAXSp6rG/cLHvgogHaGZxHFbV1aPsugfAwoi9IHivaFATJBdQE4RYfNNEVoHlqhpX1TOQbNQaACuyPYGIfEhE1onIuqOx/rF3IOWJYsIBgwCeALAsiK+oA3AtgDvz2dyxyJUmukeoiYqFmkjB6wQB4KUmxpUnSlW7ROQBAOcBmCoiNcFdRkaPLxhuWwsAJzbPLfpTnC3OtNnpC+30w7RT7PDpoWfDebhfvLzA1D2434r98fjvUuV4ws7tTa1bZOylah8aWDO9NVVe0WLb2FJrO5D4+DCsTjz/h6rGROSjAO4GUA3gJlVdLyKfB7BOVe8UkXMA/BzANABXisjnVPXUXDV/lLZNShPLpxRGEzVV9jQnTu9MlRde02DqZMVpduch25d1/c5U+cj9g6Zuy745xt4/GB57RO3cecewvYdbd6TW2LMbwxjkVW12iqS1xurUS6iJ6D4pTSwrgetEbZX9zZ3TZH/rl67sMLZErqRHd9vY+Ya+pozn6XWLine truncated
"text/plain": [
"<Figure size 720x432 with 12 Axes>"
]
"<Figure size 720x432 with 4 Axes>"
],
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAr4AAADiCAYAAACoVMTQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxcVfk/8M8nadJ93yhtaSkUS2UptLKJiALKoqIIyirgUvSrX0QRvxUBRQFRqbgCwg9olU1WQUSwAgVLodACXeje0tI93dKmafY8vz/uDeTcM+1MJjPJ3NzP+/XKqzknz733ZPrM5Myds9DMICIiIiLS0RW1dwNERERERNqCOr4iIiIikgjq+IqIiIhIIqjjKyIiIiKJoI6viIiIiCSCOr4iIiIikgjq+IqIiEiskBxM8mWSFSQnk/wpyfvyeL0TSTaS3EXy1HxdpwXteYFkNckZbXS9G0huIbmxLa6XT+r47oWeWG37xJLcUN5yCskqkmvbuy2SOeUtV5CszefvXOhIriJ5cobhEwFsAdDLzK7MY7OaW29mPczsWQAgOYTkUyTXkzSSI6MHkDyZ5JskK0muJfmlTC4Uyc+mr4ubfm5mnwTwzUwbHj6f6sLzlJOcSfLYDI/dD8CVAMaa2T6ZXrNQJa7j29GeWCRvIbks/GOxmORXsrkoyXvC8x/YVNfSJ5bkTwfM21+RXENyJ8nVJK9uycVIDiT5AMkdJLeTvL/pZ2Z2CYDTWv8rSWt1tLxtQrIfyc0tuSlA8hKSDZGOzIlNPzezAwDc1LpfJ1FGAFho7bsLVyOAZwF8MdUPSY4F8ACAHwPoDeBwAHNacP6m/Gz6mtrK9v7NzHoAGADgRQCPZHjcfgC2mllZSy/IQEH1NQuqMQWo4J9YACoBfBbBk+piAL8jeVxLLkDyeAAHtKaRUlDikLd3AxhjZr0AHAfgApJnteD8jwPYiOAFeRCAW1rRVikMccjbJr8EsCiL878a6chMz+IciRC+UZgR3tzZTvJdkqeFP5uC4O/dD8M3EN6bK5KfI/lOeHdzOsmDw/pLSf6jWdwyko80K68hOS6TNprZJjO7DcAbewi5BsCfzexfZlZvZlvNbEWmj0G+mFk9gPsBDCU5EABI9iZ5N8kNJNcxGNpQHD620wDsGz7WU8L4Y8K7xuUk5zZ/Exc+3jeSfAXAbgCjSI4hOY3kNpJLmt/5ZvAp3Z9I/jO8iTeL5AHNfv7hZsduarpRQrKI5CQGn5ZsJfkwyX7pfv9Ed3w7whPLzH5iZovNrNHMZgH4L4CMPr4I29IJwB8A/G+mx0j76iB5u8TMKptVNQI4MFVsivZ/CsBwAFeZ2Q4zqzOztzI5VtpPR8jb8HzHATgEwL2Z/ebSCkcDWILgDuWvANxNkuGnOvcD+FX4BuI/zQ8ieRCABwFcAWAggGcA/INkKYCXAHws7DTtC6AU4d9MkqMA9AAwL0ftPyY87/ywQ3lfJh2zZgaFHb13Sd5KsnsuGhU+Dl8BsBXA9rB6CoB6BK/DRwD4FICvh4/tafjg7vMlJIcC+CeAGwD0A/ADAI81daJDFyH4FKcngM0IOs8PILhRcS6A2xjcEW9yLoDrAfQFsBzAjWFbewL4D4I3pPuG7Xs+POZ/AXwewMfDn20H8Kd0v3+iO76huD+xmrepK4CPAHinBYd9D8DLZpbz9khexT5vw3fquwCsBdAdwYtiJo5B8LtPDd/lv0Hy47lql+RVrPOWZDGAPwL4DoBs7kwfwWCC0FKS14Y3HmTPVpvZXWbWAGAqgCEABmdw3JcB/NPMpplZHYJPhLoCOM7MVgKoADAOwAkAngOwnuQYBB2o/5pZY47aPwxBB/CLAEaHbfhDhscuDts4BMAnAYwH8JtWtudLJMsBVAH4BoCzzaye5GAApwO4wswqwyENtyLojKZyIYBnzOyZ8KbbNACzw3M0mWJm74R3l08FsMrM7g3vfL8F4DEA5zSLf8LMXm92N7rpzepnAGw0s8lmVm1mFeFNPiAYivljM1trZjUAfgrg7HTPK3V84//Eau4OAHPD66VFcjiAywBcl4e2SH7FPm/N7GYEdwOOBPBXADsyPHQYgrsRLwLYB8BkAE+SHJCrtknexD1vLwcwy8xaMk6zycsI7hQPQtAROg/AVTlqV0f1/goCZrY7/LZHBsftC2B1s2MbAawBMDSsegnAiQjy5SUA0xHkysfDcq5UAbjXzJaa2S4EY7hPT3NMU5s3mtnCsGP5LoAfIv0QnHQeNrM+CJ5zCxB0poFgmFEJgA3hJyrlAP6MIFdTGQHgnKbYMP54BM/nJmsi8UdH4i9A8PrdpPlqEbvxwf/zcAB7Gh4yAsATzc65CEAD0rym6N1m5IlFEsjyiUUy1RPrwPD7cgRPqmOR2ycWAIDkrxG8qH6iBWPkfgvgZ2aWaYdDCkeHyNswV98i+WkEH3N9P4PDqhDcPbg7LD9E8scAPgrgyVy3UXIqtnkb3k2+HB90Flok7KA3mU/yZwg6vr/IQfPEtR7AoU0FBok2HMC6sOolBHNj9kfQGW3qiB2L4I5+rsyD+8lAa8avG3J0s9LMtpCcCGA2yQcQdFJrAAwI77imswbAX83sG3u7TCT+JTM7JYvmrsGe7zyvAfBVM3ulJSfUHd/srUfwbgPAHp9YJwL4WPj9S8jPO0qQvB7BGJxPmdnOFhx6EoBfk9zID9bme5Xk+blsnxSUgsnbiE7IfIJl9I8JUpSlYymEvD0KwR2theHr5e8AHBW+fhZncT4DwBy1TVwPAziD5EkkSxAsxVUDYGb485cAfAJAVzNbi2BuzKkA+gNo0XwBkl0AdA6LncNyk3sBXEpyFMluACYBeDrD836C5AgGhgO4GWne2DNYReWSTM5vZksQfDryQzPbAODfACaT7BUOGzpgL0PI7gPwWZKfZjABrguD5deG7SH+aQAHkbyIZEn49RGG4/TTeBrAEJJXkOxMsifJo8Of3QHgRpIjwt9/IMkz051QHd/sFcQTi+SPAJwP4GQz29rC3+EgBMurjMMH42k+C+CJFp5H4qPd8zZ8Ub2MZN/wRf0oAN/GBxMW0nkCQF+SF4cvumcjGP7Qonf9EivtnrcA/gVgJD54vbwuPPe4cOhGuvOeFo6lRDgM41roE4q8CDt1FyIYT7sFwd+1z5pZbfjzpQB2IcgThDeMVgJ4JZP/y4iq8FxAMC63qlk77gHwFwCzEHxiUYPgU4NMHIEgvyvDf+fv7dhwvHt/AK+1oO2/BjCR5CAEk91KASxEMEnsUbhDF95nZmsAnAngagQT19Yg+PQiZZ/SzCoQDE87F8Gb2I0IVkbpnCo+xbGnIPg/3AhgGYLnOhC8+XwKwL9JViD43Y9OdZ7mNNQhS2a2hGTTE2sogLcReWIxmLjz/hOL5EoAm7N8YjVZHP7bdKfgJgC1AJaHHxsCwE1mlnY9SIusyRcev8XMqlIfIXFXQHn7BQQf8ZYieCH8AzKc9GFm20h+DsBtCGbwLgZwppltaWH7JCYKIW/DyTPvD9UguQNAnZllupPVSQCmkOwBYBOCu2Zat7cZMxvZ7PspCFYaaP5zNvv+ksjPfhopP4G93MQxsyGR8oQ0zWsE0D8cS/plM3su2qY9XOcnAH6S5typjvsN9jKZjeQ0BBN9Xw+rjgLwtJktThUffXzCulkIJhY3+Vb4FY2bjuDmQvTYlHeEzezEFHVLAJyxh/hL9nY9M1uA4PkTPa4RwWPUokl/zHw4qOQbyaaJGTVo9sRqx/a8/8QyMy/pRICCzNu7EcwWLjOzjJZIk+QpwLxdgqBT/7CZfbU92yLSkanjKyIiIiKJoDG+IiIiIpII6viKiIiISCKo4ysiIiIiidCmqzqQ1IBiyVq62bP5oryV1lDeShwpbyWmtpjZwL0FtMNyZtms8y3S0hWJck15K9lQ3kocKW8lrhpWp4vQUAcRERERSQR1fEVEREQkEdTxFREREZFEUMdXRERERBJBHV8RERERSQR1fEVEREQkEdTxFREREZFEUMdXRERERBJBHV8RERERSQR1fEVEREQkEdTxFREREZFEUMdXRERERBJBHV8RERERSQR1fEVEREQkEdJ2fEl2Ifk6ybkk3yF5fVg/heS7JN8Ov8blv7kimVHeShwpbyWOlLcSJ50yiKkB8Ekz20WyBMAMkv8Kf3aVmT2av+aJZE15K3GkvJU4Ut5KbKTt+JqZAdgVFkvCL8tno0RaS3krcaS8lThS3kqcZDTGl2QxybcBlAGYZmazwh/dSHIeyVtJds5bK0WyoLyVOFLeShwpbyUuMur4mlmDmY0DMAzAUSQPAfAjAGMAfARAPwD/l+pYkhNJziY5O0dtFsmI8lbiSHkrcaS8lbhg8AlFCw4grwOw28xuaVZ3IoAfmNln0hxrQHE27ZTEa4CZMdujlbfSPpS3EkfKW4mrhjlmNmFvEZms6jCQZJ/w+64ATgGwmOSQsI4APg9gQesbLJIbyluJI+WtxJHyVuIkk1UdhgCYSrIYQUf5YTN7muQLJAcCIIC3AXwzj+0UaSnlrcSR8lbiSHkrsdHioQ6tupg+wpCste6jt9ZQ3kr2lLcSR8pbiascDHUQEREREekI1PEVERERkURQx1dEREREEkEdXxERERFJBHV8RURERCQR1PEVERERkURQx1dEREREEkEdXxERERFJBHV8RURERCQR1PEVERERkURQx1dEREREEqFTezdARPKvtGSQV7dP93FOuRENXszO2rVOuWL3Ci/GUN/K1om0BjOIsby3QiT3ijOI8V+3Ze90x1dEREREEkEdXxERERFJhLQdX5JdSL5Oci7Jd0heH9bvT3IWyeUk/0ayNP/NFcmM8lbiSHkrcaXclbig2d7HPpEkgO5mtotkCYAZAL4L4PsAHjezh0jeAWCumd2e5lyW2ZiVQuKOHysq6poixn3/YI3VXoTGQbZWA8wsk8FLine truncated
},
"metadata": {
"needs_background": "light",
"tags": []
},
"output_type": "display_data"
}
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 263
},
"id": "OINe99pf_trt",
"outputId": "a7c9ade3-3685-40e6-a5f1-97a5084deec9"
}
},
{
"cell_type": "markdown",
"source": [
"Naturally, the image on the right is the same (this is the reference), and the other three simulations now exhibit a shift towards the right. As the differences are a bit subtle, let's visualize the difference between the target configuration and the different final states.\n",
"\n",
"The following images contain the difference between the evolved simulated and target density. Hence, dark regions indicate where the target should be, but isn't. The top row shows the original states with the initial velocity being zero, while the bottom row shows the versions after the optimization has tuned the initial velocities. Hence, in each column you can compare before (top) and after (bottom): \n"
],
"metadata": {
"id": "ootUIrQq_tru"
}
},
{
"cell_type": "code",
"execution_count": 14,
"source": [
"fig, axes = pylab.subplots(2, 3, figsize=(10, 6))\n",
"for i in range(INFLOW.shape.get_size('inflow_loc')-1):\n",
@@ -596,24 +573,44 @@
" axes[1,i].imshow(smoke.values.numpy('inflow_loc,y,x')[i,...] - smoke_final.values.numpy('inflow_loc,y,x')[3,...], origin='lower', cmap='magma')\n",
" axes[1,i].set_title(f\"Result {INFLOW_LOCATION.numpy('inflow_loc,vector')[i]}\")\n",
" pylab.colorbar(im,ax=axes[1,i])"
]
],
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": [
"<Figure size 720x432 with 12 Axes>"
],
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlEAAAF1CAYAAADIqb9jAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOy9eZycVZX//zm9L+nurJ2dLJAQQpAtRHBhERBQERdEcBlUHB1nHHXG78z4nfl+3b7OjMv35zYuY0YYcRvhy6CgosiqIAiEnSSEhOz70kt6766q8/ujntTznJuq7uru2m7V551XvXJP3We5VX0/9Zzn3vOcK6oKQgghhBAyPqqK3QBCCCGEEB+hE0UIIYQQMgHoRBFCCCGETAA6UYQQQgghE4BOFCGEEELIBKATRQghhBAyAehEjYGIXCgiuyP2ehG5MCiLiPyniHSKyOPBex8RkQMi0isiM0Y57vtEJB5sd0reP8gYiMjLIjIsIj8udltIaVNBmrhfRAZF5OFit4WUNpWgCRGpD9oxIiJfKGZbSomSd6KCTvS8iPSLyH4R+a6ITC1We1T1VFV9MDBfA+BSAAtUdY2I1AL4KoDXq+oUVT0yxuEeDbbbCAAiskpE7haRwyJiEngFHfhGEdkhIj0i8oyIXDGRzyAi94mIikhN5HOdCOBfJnI8UlioiRAR+bGI7BORoyLykoh8cDxtF5GlIvKrQFOHReTLkc/1OgB/MZ7jkeJATRyPiCwLbgKyvjEWkc8GTlJv5LU0+ExDqjoFwE+yPV4lUNJOlIh8EsCXAPwdgDYA5wJYBOAeEanLsE9NuvfzxCIA21W1L7BnA2gAsH6CxxsBcCuAG9LU1QDYBeACJL+L/wXgVhFZPJ4TiMi7AdROsH2kyFATx/GvABaraiuANwP4goicnc2Bg+/rHgD3A5gDYAEAjsR6BjWRkW8DeGICx78lcNqOvbZO4BgVQ8k6USLSCuBzAP5aVX+rqiOquh3ANQAWA3hPsN1nReS24I70KID3icgSEflDcHd5r4h8O1tvXEQaReQHwdDrBgDnOPXbReQSEbkBwPcBnBd46/8FYFOwWZeI3D/ez6yqm1T1RqQRl6r2qepnVXW7qiZU9VcAtgHI6oIRtL0NwGcA/P1420aKDzWRtn69qg4dM4PXiVke/n0A9qrqVwN9Darqc+NtIyke1ETG9l0LoAvAfeM9PhkfJetEAXgVkt767dE3VbUXwF1IDo8e4yoAtwGYiuRQ408BPA5gBoDPAnjvOM77GSR/hE8EcBmA69NtFHTiv0A41HodgFOD6qnBVEDeEJHZAJZjfHcz/wLguwD256VRJN9QE2kQke+ISD+AFwHsQ/K7yIZzAWwXkd8EUyMPishp+WgjyRvUhEPgWH4ewN9O8BBXikiHJOO6PpLDppUlpexEzQRwWFVjaer2BfXHeFRVf6GqCQCzkLwr+LSqDqvqwwDuHMd5rwHwz6raoaq7AHxzgu3PG8Gc+k8A3KyqL2a5z2oArwbwb/lsG8kr1EQaVPUvAbQAeC2SF9Oh0fdIsQDAtUh+nnkAfg3gjkxTQKQkoSaO5/8AuFFVd4+55fHcCuAUJL+fPwfwaRG5LpeNKzdK2Yk6DGBmhrnruUH9MXZFyvMAdKhqf4b6sZjnbL9jHPvmHRGpAvAjAMMAPjqOfb4D4OMZfmyIH1ATGVDVeHAhXAAg27vnAQAPq+pvVHUYwP9FclSi6E8GkqyhJiKIyBkALgHwtYnsr6obVHVvoKdHAHwDwNW5bGO5UcpO1KNI3lG+LfqmiEwBcAXsXG/0CYV9AKaLSFPkvYXjOO8+Z/sTxrFvXhERAXAjkoGJb1fVkSx3bQWwGsAtIrIfYbDhbhF5be5bSvIENTE2Ncg+Juo52O+J+Ac1YbkQyViwncFv/f8A8HYReWqCx1MAkpumlScl60SpajeSAYP/JiKXi0ht8CTarQB2Izkak26/HQDWAfisiNSJyHkArhzHqW8F8D9FZJqILADw15P4GONCkjQAqAvsBhGpj2zyXSTvkq9U1YFxHLobyTunM4LXG4L3zwbw2KQbTgoCNWE1ISLtInKtiEwRkWoRuQzAdcg+mPbHAM4NAoCrAXwCyZGLjbn/JCQfUBPHXSfWInkTcey3/t+RnKa+LMtjXxV8JhGRNQA+BuCOXH+GcqJknSgAUNUvA/hHJIfZjyJ5wd8F4OLIEznpeDeA8wAcAfAFALcgEicRPCWRaQTmc0gOzW4D8DtkEGE2BAGr/ziOXRYhOcVwLFh8AMGTHCKyCMCHkRTGfglzeLx7rINqkv3HXgAOBVUHgmkM4gnURKgJJO+SP4LkxbITye/kE6qaVWyLqm5C8umtfw/2vwrAm6kJv6AmQk2oar/zW98LYFBVD6U/1HFcC2ALgB4APwTwJVW9eRxtqzhEtfxHs0XkFgAvqupnit2WY4jIewF8D8nYpvOOJVIrYns2AZgP4FZV/UAx20LyDzWRVXvuQfIJvsdV9eJitoXkH2pizLbUAziAZJ7BL6vq54rVllKiLJ0oETkHQAeSdwmvB/ALJDvg00VtGCFFgpogxEJNkFxQ0tN5k2AOgAeRHMr8JoCPUBh+E8Q7bBKRLSLyqTT1fysiG0TkOUkua7OoGO0sYaiJMoOamDTURJlRDE2U5UgUKS+CoN+XkEyctxvJpwuvU9UNkW0uAvCYqvZLMkHchar6zqI0mJA8Q00QYimWJsp1JIqUF2sAbFHVrUHQ78+QDAJOoaoPRHK+/AnJfEGElCvUBCGWomiCThTxgfmwie12B+9l4gYAv8lriwgpLtQEIZaiaKKQK1mjtaZJZ9W3FfKUpAAcGurG0Vj/qAnZLrtsjR450p227sknX1oPYDDy1lpVXTuRtojIe5BMLHrBRPYvNG21TTqbmig7Dgx1o3uEmpgIbbVN2l5HTZQbB4fLUxMFdaJm1bfhy6fcUMhTkgLw9xtvHHObI4e78Nij30lbV1N/6aCqrh5l9z2w2YEXBO8ZROQSAP8E4IIx8sOUDLPr2/CtVe8rdjNIjvnoCz8YcxtqIj3tdW345qr3F7sZJMd87IX/HHMbHzVRUCeKVDAKIB6f6N5PAFgmIkuQFMW1AN4V3UBEzkQyn8rlqnpwEi0lpDBQE4RYPNQEnShSIBRIJCa2p2pMRD4K4G4A1QBuUtX1IvJ5AOuCDNVfATAFwP8TEQDYqapvzk3bCckH1AQhFv80QSeKFAYFEItNfHfVuwDc5bz36Uj5kgkfnJBiQE0QYvFQE3SiSGFQnZQ4CCk7qAlCLB5qgk4UKRAK0YkN0xJSnlAThFj80wSdKFIYFEBswgGDhJQf1AQhFg81QSeKFAYPh2kJySvUBCEWDzVBJ4oUCAXifomDkPxCTRBi8U8TdKJIYVAACS52TUgKaoIQi4eaoBNFCoR/w7SE5BdqghCLf5qgE0UKgypk4ploCSk/qAlCLB5qgk4UKQyTTKJGSNlBTRBi8VATdKJIgVDv5roJyS/UBCEW/zRBJ4oUBg/vMAjJK9QEIRYPNUEnihQG1cmszk1I+UFNEGLxUBN0okjh8GyYlpC8Q00QYvFME3SiSGHwMBMtIXmFmiDE4qEm6ESRAuHfMC0h+YWaIMTinyboRJHC4OHCkoTkFWqCEIuHmqATRQqDh+n8Cckr1AQhFg81QSeKFAj17g6DkPxCTRBi8U8TdKJIYfDw0VVC8go1QYjFQ03QiSKFwcO5bkLyCjVBiMVDTdCJIoXDs7luQvIONUGIxTNN0IkihUH9m+smJK9QE4RYPNQEnShSGBRAPFHsVhBSOlAThFg81ASdKFIg/Fudm5D8Qk0QYvFPE1XFbkA5oZCMr4rnWMBgulcWiMjlIrJJRLaIyKfS1J8vIk+JSExErs518yuZppqYec2f0pd6zWu2r2n1w+ZVJZp6KWBeFQ81UZZE+3y6l6sDaiKCh5qgE0UKgyoQS6R/jYGIVAP4NoArAKwEcJ2IrHQ22wngfQB+muOWE5IfqAlCLB5qYkwnSkQaRORxEXlWRNaLyOeC95eIyGOBx3eLiNTlqlGkPNF4Iu0rC9YA2KKqW1V1GMDPAFxljq26XVWfA5D3CXVqguQKaoIQi2+ayGYkagjA61T1dABnALhcRM4F8CUAX1PVkwB0ArghV40iZcixdP7pXmMzH8CuiL07eK9YUBNk8lAThFg81MSYgeWqqgB6A7M2eCmA1wF4V/D+zQA+C+C7uW9ifqlywpVqJfMfq3uk2ti7BuzXd2Qo83lm1Ft7XoOd422u8euxznFzbJg2PTNFZF3EXquqawvQqglRbpporLZ9b8X8Q8aefp7t91Unz02VtbPX1A0/ecDY256flio/dXiaqRtM2Hu46XUxYzdXUxMRm5ooIaqd68Sc5v5UubF2xNTtPtpi7K19janygUGrgaYae9wFjdREhJLURFZP5wVzjU8COAnJOceXAXSp6rG/cLHvgogHaGZxHFbV1aPsugfAwoi9IHivaFATJBdQE4RYfNNEVoHlqhpX1TOQbNQaACuyPYGIfEhE1onIuqOx/rF3IOWJYsIBgwCeALAsiK+oA3AtgDvz2dyxyJUmukeoiYqFmkjB6wQB4KUmxpUnSlW7ROQBAOcBmCoiNcFdRkaPLxhuWwsAJzbPLfpTnC3OtNnpC+30w7RT7PDpoWfDebhfvLzA1D2434r98fjvUuV4ws7tTa1bZOylah8aWDO9NVVe0WLb2FJrO5D4+DCsTjz/h6rGROSjAO4GUA3gJlVdLyKfB7BOVe8UkXMA/BzANABXisjnVPXUXDV/lLZNShPLpxRGEzVV9jQnTu9MlRde02DqZMVpduch25d1/c5U+cj9g6Zuy745xt4/GB57RO3cecewvYdbd6TW2LMbwxjkVW12iqS1xurUS6iJ6D4pTSwrgetEbZX9zZ3TZH/rl67sMLZErqRHd9vY+Ya+pozn6XWLine truncated
},
"metadata": {
"needs_background": "light",
"tags": []
}
}
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 410
},
"id": "7lm9Rrel640-",
"outputId": "6ef7914a-8f63-49af-9a1a-fde1aed37682"
}
},
{
"cell_type": "markdown",
"metadata": {
"id": "cuGIjhGx_trw"
},
"source": [
"These difference images clearly show that the optimization managed to align the upper region of the plumes very well. Each original image (at the top) shows a clear misalignment in terms of a black halo, while the states after optimization largely overlap the target smoke configuration of the reference, and exhibit differences closer to zero for the front of each smoke cloud.\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. E.g., this is noticeable at the stems of the smoke plumes, which still show a black halo after the optimization. The optimization was not able to shift the inflow position, and hence needs to focus on aligning the upper regions of the plumes.\n"
]
],
"metadata": {
"id": "cuGIjhGx_trw"
}
},
{
"cell_type": "markdown",
"metadata": {
"id": "jHdQureE_trx"
},
"source": [
"## Conclusions\n",
"\n",
@@ -621,13 +618,13 @@
"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 jointly solved a single inverse problem. There were three cases to solve as a mini-batch, of course, but nonetheless the setup still represents a direct optimization. Thus, in line with the PINN example of {doc}`physicalloss-code` we've not really dealt with a _machine learning_ task here. However, that will be the topic of the next chapters.\n"
]
],
"metadata": {
"id": "jHdQureE_trx"
}
},
{
"cell_type": "markdown",
"metadata": {
"id": "ffaDxZ3x_try"
},
"source": [
"## Next steps\n",
"\n",
@@ -636,7 +633,10 @@
"- Modify the setup of the simulation to differ more strongly across the four instances, run longer, or use a finer spatial discretization (i.e. larger grid size). Note that this will make the optimization problem tougher, and hence might not converge directly with this simple setup.\n",
"\n",
"- As a larger change, add a multi-resolution optimization to handle cases with larger differences. I.e., first solve with a coarse discretization, and then uses this solution as an initial guess for a finer discretization. \n"
]
],
"metadata": {
"id": "ffaDxZ3x_try"
}
}
],
"metadata": {
@@ -665,4 +665,4 @@
},
"nbformat": 4,
"nbformat_minor": 1
}
}
+1 -2
View File
@@ -174,8 +174,7 @@
"source": [
"A lot has happened in this `step()` call: we've advected the smoke field, added an upwards force via a Boussinesq model, advected the velocity field, and finally made it divergence free via a pressure solve.\n",
"\n",
"The Boussinesq model uses a multiplication by a tuple `(0, buoyancy_factor)` to turn the smoke field into a staggered, 2 component force field, sampled at the locations of the velocity components via the `>>` operator. \n",
"This operator makes sure the individual force components are correctly interpolated for the velocity components of the staggered velocity. `>>` could be rewritten via the sampling function `at()` as `(smoke*(0,buoyancy_factor)).at(velocity)`. However, `>>` also directly ensure the boundary conditions of the original grid are kept. Hence, it internally also does `StaggeredGrid(..., extrapolation.ZERO,...)` for the force grid. Long story short: the `>>` operator does the same thing, and typically produces shorter and more readable code.\n",
"The Boussinesq model uses a multiplication by a tuple `(0, buoyancy_factor)` to turn the smoke field into a staggered, 2 component force field, sampled at the locations of the velocity components via the `at()` function. This function makes sure the individual force components are correctly interpolated for the velocity components of the staggered velocity. Note that this also directly ensure the boundary conditions of the original grid are kept. It internally also does `StaggeredGrid(..., extrapolation.ZERO,...)` for the resulting force grid. \n",
"\n",
"The pressure projection step in `make_incompressible` is typically the computationally most expensive step in the sequence above. It solves a Poisson equation for the boundary conditions of the domain, and updates the velocity field with the gradient of the computed pressure.\n",
"\n",