dp ns code update

diffphys-code-ns.ipynb

@@ -47,6 +47,10 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
+      "/Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v2/phi/physics/_boundaries.py:384: SyntaxWarning: \"is\" with a literal. Did you mean \"==\"?\n",
+      "  return self.velocity is 0 and self.angular_velocity is 0\n",
+      "/Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v2/phi/physics/_boundaries.py:384: SyntaxWarning: \"is\" with a literal. Did you mean \"==\"?\n",
+      "  return self.velocity is 0 and self.angular_velocity is 0\n",
       "/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",
@@ -178,20 +182,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Velocity tensor shape: (inflow_loc=4, x=32, y=40, vector=2)\n",
-      "WARNING:tensorflow:From /Users/thuerey/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/ops/math_grad.py:297: setdiff1d (from tensorflow.python.ops.array_ops) is deprecated and will be removed after 2018-11-30.\n",
-      "Instructions for updating:\n",
-      "This op will be removed after the deprecation date. Please switch to tf.sets.difference().\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "smoke = DOMAIN.grid(math.zeros(inflow_loc=4))\n",
     "velocity = initial_velocity = DOMAIN.staggered_grid(0) * math.ones(inflow_loc=4)\n",
@@ -216,22 +209,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAABHYAAAFgCAYAAADehv6hAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABa40lEQVR4nO3de7xkVX3n/c+XqyYKDXSLLWgaA06CiSEPHSQXlQFEMo8TnAkiymgzAyGMMZNMEgcMGTVEMjgmIU9mjEkrCF6BYNROQiSAEjXxQqMEQaO0RENjc2mgESVCGn7PH3uf7upz6nTXOXVq17l83q9Xvbr22mvvWnUu395n1Vprp6qQJEmSJEnSwrPbuBsgSZIkSZKk2bFjR5IkSZIkaYGyY0eSJEmSJGmBsmNHkiRJkiRpgbJjR5IkSZIkaYGyY0eSJEmSJGmBsmNHkiRJkiRpkiQnJvlqkg1Jzu2z/4VJvpBka5KTJ+1bk+T29rGmp/zIJF9qz/lHSTJsO+3YkSRJkiRJ6pFkd+DtwM8ChwOvTHL4pGr/DJwOfGDSsfsDbwKeDxwFvCnJfu3udwC/ABzWPk4ctq127EiSJEmSJO3oKGBDVd1RVY8BlwMn9Vaoqm9U1S3AE5OOfQlwbVU9UFUPAtcCJyZZCexTVZ+tqgLeA7xs2IbuMewJJM3OiUltHvIcN8E1VTV0D6+kxc28kdSVYfPGrJE0qDnIm9uA7/UUra2qtT3bBwF39mxvpBmBM4h+xx7UPjb2KR+KHTvSmGwG1g95jsDyuWiLpMXNvJHUlWHzxqyRNKg5yJvvVdXquWrPODkVSxqn3XYb7iFJgzJvJHVlxFkzwGKmeye5ot3/uSSreva9oS3/apKXzN2bljQWo82bu4Bn9mwf3JYNc+xd7fPZnHNaXqlJ4+QfWpK6Yt5I6soIs2bAxUzPAB6sqkOBi4C3tsceDpwKPJdmsdI/bs8naaEa7bXNjcBhSQ5JshdNfqwbsGXXACck2a9dNPkEmqmmm4BvJzm6vRvWa4CPzvyN78grNUmSJEkLxS4XM223L2ufXwUc1/4BdRJweVU9WlX/BGxozydJU1TVVuB1NJ00XwGurKrbkpyf5OcAkvxEko3Ay4E/TXJbe+wDwO/QdA7dCJzflgG8FngXTQZ9HfjrYdvqGjvSuCR+Ci6pG+aNpK4MnzfLk/QumzGbxUy31amqrUkeAg5oyz876dihFy2VNCYdXN9U1dXA1ZPK3tjz/EZ2nFrVW+8S4JI+5euBH5nLdtqxI42Tf2hJ6op5I6krw+XN5sWymKmkDnh9A9ixI42Pn6BL6op5I6kro8+bQRYznaizMckewL7A/QMeK2mh8PpmG78KkiRJkhaKQRYzXQesaZ+fDHy8qqotP7W9a9YhwGHA5ztqtySNjCN2pHGyh1lSV8wbSV0ZYd60a+ZMLGa6O3DJxGKmwPqqWgdcDLw3yQbgAZrOH9p6VwJfBrYCv1RVj4+ssZJGz+sbwI4dabwMIkldMW8kdWX8i5l+j+YONf2OvQC4YKQNlNQdr28AO3ak8XFOqKSumDeSumLeSOqKebONHTvSOBlEkrpi3kjqinkjqSvmDeDiyZIkSZIkSQuWI3akcXHooKSumDeSumLeSOqKebONHTvSOBlEkrpi3kjqinkjqSvmDWDHjjReBpGkrpg3krpi3kjqinkDuMaOJEmSJEnSguWIHWlcnBMqqSvmjaSumDeSumLebGPHjjROBpGkrpg3krpi3kjqinkD2LEjjY89zJK6Yt5I6op5I6kr5s02duxI42QQSeqKeSOpK+aNpK6YN4CLJ0uSJEmSJC1YjtiRxskeZkldMW8kdcW8kdQV8wawY0caH+eESuqKeSOpK+aNpK6YN9vYsSONk0EkqSvmjaSumDeSumLeAK6xIy1qSU5M8tUkG5Kc22f/C5N8IcnWJCdP2rcmye3tY013rZYkSZIkDcoRO9K4jHjoYJLdgbcDLwY2AjcmWVdVX+6p9s/A6cBvTDp2f+BNwGqggJvaYx8cWYMljY5DlSV1xbyR1BXzZhs7dqRxGm0QHQVsqKo7AJJcDpwEbOvYqapvtPuemHTsS4Brq+qBdv+1wInAB0fZYEkj5IWPpK6YN5K6Yt4AduxI4zXaIDoIuLNneyPw/CGOPWiO2iVpHLzwkdQV80ZSV8wbwI4daXzmZujg8iTre7bXVtXaYU8qaZFxqLKkrpg3krpi3mzjV0Fa2DZX1eqeR2+nzl3AM3u2D27LBjHMsZKWoAEWa987yRXt/s8lWdWz7w1t+VeTvKTThkuSJC1wjtiRxmm0Pcw3AoclOYSmU+ZU4FUDHnsN8LtJ9mu3TwDeMPdNlNSZ8S/WfgbwYFUdmuRU4K3AK5IcTpNPzwWeAVyX5DlV9fjIGixptPwEXVJXzBvAETvS+EwMHRzmsRNVtRV4HU0nzVeAK6vqtiTnJ/m5pgn5iSQbgZcDf5rktvbYB4DfoekcuhE4f2IhZUkL0Ijzhp7F2qvqMWBisfZeJwGXtc+vAo5Lkrb88qp6tKr+CdjQnk/SQjRs3kjSoMybbRyxI43TiAOlqq4Grp5U9sae5zfSTLPqd+wlwCUjbaCk7ox2Ta9BFmvfVqeqtiZ5CDigLf/spGNdrF1ayBbZH0yS5jHzBrBjR5IkDWZzVa0edyMkSZK0Izt2pHGyh1lSV0abN4MsuD5RZ2OSPYB9gfsHPFbSQuL1jaSumDfAIl9jJ8k3khzfPv/NJO/q6HWT5N1JHkzy+S5ec1yS/Nck9yT5TpIDklSSQ0f8mi9J8pFRvsaA7fhQkp8d4gSjXvNCHTJvRm+J583vJ/mvQ5xg1HmzbbH2JHvRLIa8blKddcCa9vnJwMerqtryU9u7Zh0CHAYs6p/lYZg1o7fEs2a4a5vmJF7bLBLmzegt8bwZ7tqmOYl50xrbu0lyanu70+8mubd9/tp2IcU5V1W/W1VnDnueJKvaX7idjXb6GZo7gxxcVSNfALJt0yeSPJLkHycCeJbnOj3JpwesuyfwB8AJVfWUqrp/tq87QxcAF/a043eSfCnJ1iRv7tPOVyX5Zvuz9pEk+w/yIj3f6+/0PP5nT5W3Am8Z6p3YsdMJ82bumDfT502S/zfJp5NsSXJ3kncleeqgL9R+r7/bkze9F9C/B/xm22kyOyPMm0EWawcuBg5IsgH4NeDc9tjbgCuBLwMfA35pod4Ry6yZO2bNzq9teupdMpM/BDu5tgGvbTpg3swd82aXf0utSPKBJA+l6Wx7/6AvNPJrGxh53iQ5MclXk2xIcm6f/XsnuaLd/7kkq9ry05Lc3PN4IskR7b4b2nNO7HvaUF8DxtSxk+TXgf8PeBvwdOBA4Gzgp4G+39g0t1JdKH4A+EZVfXemB+4i5KbzQeCLNItQngdclWTFLM4zUwcCTwJu6+C1gOYuTsC+VdW70OYG4H8Af9Wn/nOBPwVeTdPeR4A/nuHLLmvD9ilV9TsThVX1eWCfJLNbc8IRO50wb6Zn3uzcTPOGZmrRW2hu2f3DNAsAv22GL/tjPXmz7QK6qjYB/wj83PSH7kQHeVNVV1fVc6rqB6vqgrbsjVW1rn3+vap6eVUdWlVHVdUdPcde0B73b6rqr2f1HsfMrJmeWbNzs8iaieN+BvjBWb7saK5tmoZ5bTNi5s30zJudm2Xe/DlwN/As4Gk0HTIzMZprGxh53rS/N28HfhY4HHhlksMnVTsDeLCqDgUuoukgp6reX1VHVNURNH+L/lNV3dxz3GkT+6vq3ll/DVqdp2eSfYHzgddW1VVV9XA1vlhVp1XVo229S5O8I8nVSb4L/Ns0n4Z+Mcm3k9w5uUcxyavTjMy4P8l5k/a9Ocn7eraPTvL3aT5Z/Yckx/Tsu6Htufy7JA8n+Zsky9vdn2z/3dL2Ov7kpNc5A3gX8JPt/t9uy3+h7cV7IMm6JM/oOaaS/FKS24Hb+3zN/jrJ6yaV/UOS/5jkOcD/A7ypqv6lqj4EfAn4+V18KwaSZgjmbyS5pe2lvSLJk9rX/WrP1+LjfY7dN8l7ktzXfl9+K8lu7b5vJjmyfX5a+zV4brt9RqYfHvizwN/2FlTVZe0fAg/3qX8a8BdV9cmq+g7wP4H/mBl8ir4LNwD/7xydS3PMvDFv2n2d5E1VfaCqPlZVj1TVg8A7aS6y58oNmDfzkllj1rT7Osma9nx7AP8H+OWZfwV26QbMmnnLvDFv2n2d5E2SE2jWwXt9VT1UVf9aVV+c3VejrxuY33lzFLChqu6oqseAy4GTJtU5CbisfX4VcFwyZeTcK9tjR2Yc3eI/CewNfHSAuq+iGSr2VODTwHeB1wDLaH4A/muSlwGk6Tl7B01v2DNoelz73sY5yUE0PZJvAfYHfgP4UHbsmX0V8J9peiX3ausAvLD9d+KTjs/0nruqLqbpMf9Mu/9NSY4F/hdwCrAS+CZTv7Evo7k17OQeQGh6kV/Z0/7DaXqy/wp4LnBHVfX+Iv5DWz5XTgFOBA4BngecXlVf63mNZVV1bJ/j/g/NJ9jPBl5E8737z+2+vwWOaZ+/CLiD7V/bFzEpcHr8KNtDcBDPpfl6AFBVXwceA54zg3N8M8nGNHN9l0/a9xXgx2Zwrh05YmfUzBvzBrrLm8leyMw/hftkmmlcf552KG8P82b+MmvMGug2a/478MmqumWGx00Y3bUNmDWjZd6YN9Bd3hzd1r+s7fC7McmLZnA8jPLaBobNm+VJ1vc8zpp09oOAO3u2N7ZlfetUMzX9IZrfn16voPk57PXuNNOw/mefjqAZG0d6Lqe5ZerWiYKe3t5/SfLCnrofraq/q6onqhnCfUNVfandvoXmizPxg3Uy8JftyIxHaUZmPDFNG/4TcHU1w8afqKprgfXAv+up8+6q+lpV/QvN3P8jhnjPpwGXVNUX2ra9gaYXelVPnf9VVQ+0rzfZh4EjkvxAz/n+vD3XU2h+eHo9RBPgc+WPqupbVfUA8BcM8LVLine truncated
-      "text/plain": [
-       "<Figure size 1152x360 with 8 Axes>"
-      ]
-     },
-     "metadata": {
-      "needs_background": "light"
-     },
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "fig, axes = pylab.subplots(1, 4, figsize=(16, 5))\n",
     "for i in range(INFLOW.shape.inflow_loc):\n",
@@ -258,26 +238,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "65.591805\n",
-      "65.37787\n",
-      "65.1976\n",
-      "64.12548\n",
-      "64.667114\n",
-      "64.404465\n",
-      "74.69762\n",
-      "76.985146\n",
-      "78.58485\n",
-      "82.51491\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "initial_velocity = DOMAIN.staggered_grid(0) * math.ones(inflow_loc=4)\n",
     "LR = 1e-02\n",
@@ -303,7 +266,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The loss should have gone down significantly, from above 60 to below 40, and now we can visualize how well the reconstruction for the last frame turned out for simulations 0 to 2."
+    "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**"
    ]
   },
   {
@@ -314,9 +279,23 @@
    "source": [
     "fig, axes = pylab.subplots(1, 4, figsize=(16, 5))\n",
     "for i in range(INFLOW.shape.inflow_loc):\n",
-    "  im = axes[i].imshow(grad.inflow_loc[i].vector[0].numpy('y,x'), origin='lower', cmap='bwr')\n",
+    "  im = axes[i].imshow(initial_velocity.inflow_loc[i].vector[0].numpy('y,x'), origin='lower', cmap='bwr')\n",
     "  pylab.colorbar(im, ax=axes[i])\n",
-    "  axes[i].set_title(f\"Gradient for v0 - Inflow {math.to_int(INFLOW_LOCATION).inflow_loc[i].vector.unstack_spatial('x,y')}\")\n",
+    "  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()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axes = pylab.subplots(1, 4, figsize=(16, 5))\n",
+    "for i in range(INFLOW.shape.inflow_loc):\n",
+    "  im = axes[i].imshow(initial_velocity.inflow_loc[i].vector[1].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()"
    ]
   },
@@ -324,15 +303,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Naturally, the image on the right is the same (this is simulation (3), our reference), and the other three simulations now exhibit a noticeable shift towards the right. \n",
+    "## Re-simulation\n",
     "\n",
-    "**XXX TODO** This is especially visible for the far right simulation (3), which was moving straight up in its original state, and now has an off-center plume shape that tries to match the reference on the left. \n",
-    "\n",
-    "Note that all three simulations need to \"work\" with a fixed inflow, hence they cannot simply \"produce\" marker density out of the blue to match the target. Also each simulation needs to take into account how the non-linear model equations change the state of the system over the course of 20 time steps. So the optimization goal is quite difficult, and it is not possible to exactly satisfy the constraints to match the reference simulation in this scenario. This is exactly what makes it an interesting test case, though.\n",
-    "\n",
-    "**XXX TODO**  Now we can also have a look at the now-optimized initial velocity field. Here we see the y-components (from `data[...,0]`).\n",
-    "\n",
-    "... re-simulate "
+    "Most importantly, we can also visualize how the full simulation over the course of 20 steps turns out given the new initial conditions for each of the inflow locations. This is what happened internally at optimization time, and was measured by our loss function, but it's good to get an intuition for which solution the optimization has found.\n"
    ]
   },
   {
@@ -359,359 +332,34 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "This example illustrates how the differentiable physics approach can easily be extended towards significanlty more complex PDEs. Above, we've optimized for a mini-batch of 20 steps of a full Navier-Stokes solver.\n",
+    "Naturally, the image on the right is the same (this is simulation (3), our reference), and the other three simulations now exhibit a noticeable shift towards the right. This is especially visible for the far left simulation (0), which was strongly curved in its original state, and now has a shape that tries to match the straighter reference on the right. \n",
     "\n",
-    "This is a powerful basis to bring NNs into the picture! As you might have noticed, our degrees of freedom were still a regular grid, and we've solved a single inverse problem. There was more than one case to solve, but nonetheless no question of generalization and no neural network structure.\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "<br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>\n",
-    "\n",
-    "\n",
-    "--- \n",
-    "\n",
-    "## END UP TO HERE , REMOVE REST?\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "xxx\n",
-    "\n",
-    "#use keras GD?\n",
-    "\n",
-    "# \"reset\" necessary???\n",
-    "#smoke = DOMAIN.grid(math.zeros(inflow_loc=4))\n",
-    "if 0:\n",
-    "    velocity = initial_velocity = DOMAIN.staggered_grid(0) * math.ones(inflow_loc=4)\n",
-    "    smoke = DOMAIN.grid(math.zeros(inflow_loc=4)) + INFLOW\n",
-    "\n",
-    "    with math.record_gradients(initial_velocity.values):\n",
-    "        buoyancy_force = smoke * (0, 1) >> velocity\n",
-    "        velocity, _, _, _ = fluid.make_incompressible(velocity + buoyancy_force, DOMAIN)\n",
-    "        for _ in range(2):\n",
-    "            smoke,velocity = step(smoke,velocity)\n",
-    "\n",
-    "        loss = field.l2_loss(smoke - field.stop_gradient(smoke.inflow_loc[-1]))\n",
-    "        print(format(loss))    \n",
-    "        grad = math.gradients(loss, initial_velocity.values)\n",
-    "        print(type(grad))\n",
-    "\n",
-    "    print(\"ok\")\n",
-    "# ---\n",
-    "\n",
-    "# https://www.tensorflow.org/api_docs/python/tf/keras/optimizers/Optimizer\n",
-    "opt = tf.keras.optimizers.SGD(learning_rate=1e-02)\n",
-    "#opt = tf.keras.optimizers.Adam(learning_rate=0.1)\n",
-    "    \n",
-    "velocity = initial_velocity = DOMAIN.staggered_grid(0) * math.ones(inflow_loc=4)\n",
-    "smoke = DOMAIN.grid(math.zeros(inflow_loc=4)) + INFLOW\n",
-    "\n",
-    "#with math.record_gradients(initial_velocity.values): # PHI\n",
-    "with tf.GradientTape() as tape: # same\n",
-    "    buoyancy_force = smoke * (0, 1) >> velocity\n",
-    "    velocity, _, _, _ = fluid.make_incompressible(velocity + buoyancy_force, DOMAIN)\n",
-    "    for _ in range(2):\n",
-    "        smoke,velocity = step(smoke,velocity)\n",
-    "\n",
-    "    loss = field.l2_loss(smoke - field.stop_gradient(smoke.inflow_loc[-1]))\n",
-    "    #print(type(loss.native()))\n",
-    "    print(\"Loss \"+format(loss))\n",
-    "    \n",
-    "    #grad = math.gradients(loss, initial_velocity.values) # PHI\n",
-    "    \n",
-    "    #print(type(grad.vector[0]))\n",
-    "    print(type(grad.vector[0].native()))\n",
-    "    #print(format(grad.vector[0].shape))\n",
-    "    print(type(initial_velocity.vector[0].values.native()))\n",
-    "    #print(format(initial_velocity.vector[0].shape))\n",
-    "\n",
-    "    grads = tape.gradient(loss.native(), [initial_velocity.vector[0].values.native()] ) # works\n",
-    "    \n",
-    "    # a = grad.vector[0].native() + initial_velocity.vector[0].values.native()  # works PHI\n",
-    "    #opt.apply_gradients(zip([grad.vector[0].native()], [initial_velocity.vector[0].values.native()] ))\n",
-    "\n",
-    "    opt.apply_gradients(zip(grads, [initial_velocity.vector[0].values.native()] ))\n",
-    "    # -> same error as opt.minimize?\n",
-    "    \n",
-    "    # minimize:\n",
-    "    #lossL = lambda: (field.l2_loss(smoke - field.stop_gradient(smoke.inflow_loc[-1]))).native()\n",
-    "    #opt.minimize(lossL, [initial_velocity.vector[0].values.native()])\n",
-    "\n",
-    "\n",
-    "\n",
-    "#initial_velocity = initial_velocity - 1e-02 * grad\n",
-    "#step_count = opt.minimize(loss, [var1]).numpy()\n",
-    "#step_count = opt.minimize(loss, [initial_velocity.values.native()]).numpy()\n",
-    "\n",
-    "#print(step_count)\n",
+    "Note that all three simulations need to \"work\" with a fixed inflow, hence they cannot simply \"produce\" marker density out of the blue to match the target. Also each simulation needs to take into account how the non-linear model equations change the state of the system over the course of 20 time steps. So the optimization goal is quite difficult, and it is not possible to exactly satisfy the constraints to match the reference simulation in this scenario. This is exactly what makes it an interesting test case, though.\n",
     "\n"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
-   "source": [
-    "..."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "NlJMJikaHOL6"
-   },
-   "outputs": [],
    "source": [
     "\n",
+    "## Conclusions\n",
     "\n",
+    "This example illustrates how the differentiable physics approach can easily be extended towards significanlty more \n",
+    "complex PDEs. Above, we've optimized for a mini-batch of 20 steps of a full Navier-Stokes solver.\n",
     "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "\n",
-    "xxx UP TO HERE xxx\n",
-    "\n",
-    "\n",
-    "world = World()\n",
-    "fluid = world.add(Fluid(Domain([40, 32], boundaries=CLOSED), buoyancy_factor=0.05, batch_size=4), physics=IncompressibleFlow())\n",
-    "world.add(Inflow(Sphere(center=centers, radius=3), rate=0.2));\n",
-    "fluid.velocity = variable(fluid.velocity)  # create TensorFlow variable\n",
-    "initial_state = fluid.state  # Remember the state at t=0 for later visualization\n",
-    "session.initialize_variables()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "id": "vSdGIEzCgq5-"
-   },
-   "source": [
-    "The simulation now contains variables in the initial state.\n",
-    "Since all later states depend on the value of the variable, the `step` method cannot directly compute concrete state values.\n",
-    "Instead, `world.step` will extend the TensorFlow graph by the operations needed to perform the step.\n",
-    "\n",
-    "To execute the graph with actual data, we can use `session.run`, just like with regular TensorFlow 1.x. While `run` would usually be used to infer predictions from a learning model, it now executes the graph of simulation steps."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "wSrIezfWHjcQ"
-   },
-   "outputs": [],
-   "source": [
-    "world.step()\n",
-    "pylab.imshow(np.concatenate(session.run(fluid.density).data[...,0], axis=1), origin='lower', cmap='magma')"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "id": "iJc6UdYHhtOH"
-   },
-   "source": [
-    "Let's build a graph for the full simulation."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "b9xHtdDQRrjL"
-   },
-   "outputs": [],
-   "source": [
-    "for step in range(20):\n",
-    "  print('Building graph for time step %d' % step)\n",
-    "  world.step(dt=1.5)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "id": "2VQ92g2rs6wM"
-   },
-   "source": [
-    "When calling `session.run` now, the full simulation is evaluated using TensorFlow operations.\n",
-    "This will take advantage of your GPU, if available.\n",
-    "If you compile Φ<sub>Flow</sub> with [CUDA support](https://github.com/tum-pbs/PhiFlow/blob/master/documentation/Installation_Instructions.md), the TensorFlow graph will also use optimized operators for efficient simulation and training runs.\n",
-    "\n",
-    "The `session.run()` call of the following code block will now retrieve the final fluid density, and for that it actually needs to process all 20 simulation steps of the TF graph that we just constructed."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "TA6Ibs-mXsTc"
-   },
-   "outputs": [],
-   "source": [
-    "print('Computing frames...')\n",
-    "original_state = np.concatenate(session.run(fluid.density).data[...,0], axis=1)\n",
-    "pylab.imshow(original_state, origin='lower', cmap='magma')"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "id": "IClfRMfoyGUa"
-   },
-   "source": [
-    "Next, we define the *loss* function. This is the value we want to decrease via optimization.\n",
-    "For this example, we want the marker densities of all final simulation states to match the left-most one, called `target`, in terms of an $L^2$ norm.\n",
-    "\n",
-    "For the optimizer, we choose again gradient descent for this example."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "7KPpyIwjYETi"
-   },
-   "outputs": [],
-   "source": [
-    "target = session.run(fluid.density).data[0,...]\n",
-    "loss = math.l2_loss(fluid.density.data[1:,...] - target)\n",
-    "optim = tf.train.GradientDescentOptimizer(learning_rate=0.01).minimize(loss)\n",
-    "session.initialize_variables()\n",
-    "\n",
-    "print('Initial loss: %f' % session.run(loss))"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "id": "AALD66-N0U5F"
-   },
-   "source": [
-    "With the loss and optimizer set up, all that's left is to run the actual optimization."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "pvvF6xqmaRLX"
-   },
-   "outputs": [],
-   "source": [
-    "for optim_step in range(100):\n",
-    "    _, loss_value = session.run([optim, loss])\n",
-    "    if optim_step<3 or optim_step%10==0: \n",
-    "        print('Optimization step %d, loss %f %s' % (optim_step, loss_value, '' if optim_step else '. The first step sets up the adjoint graph.'))\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, of course, but nonetheless no question of generalization and no neural network structure. These will be topics of the next chapters.\n"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The loss should have gone down significantly, from above 60 to below 40, and now we can visualize how well the reconstruction for the last frame turned out for simulations 1 to 3."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "SQBtCmhZaYYj"
-   },
-   "outputs": [],
-   "source": [
-    "optimized_state = np.concatenate(session.run(fluid.density).data[...,0], axis=1)\n",
-    "pylab.title(\"Original state at t=20 (top), and optimized version (bottom):\")\n",
-    "pylab.imshow( np.concatenate([optimized_state, original_state],axis=0), origin='lower', cmap='magma')"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "id": "HP7aDQfpKifp"
-   },
-   "source": [
-    "Naturally, the image on the left is the same (this is simulation (0), our reference), and the other three simulations now exhibit a noticeable curved shape towards the left. This is especially visible for the far right simulation (3), which was moving straight up in its original state, and now has an off-center plume shape that tries to match the reference on the left. \n",
+    "## Next Steps\n",
     "\n",
-    "Note that this simulation (like both others) needs to \"work\" with a fixed inflow, hence it cannot simply \"produce\" marker density out of the blue to match the target. Also it 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 tough, and it is not possible to exactly satisfy the constraints to match simulation (0) in this scenario. This is exactly what makes it an interesting test case, though.\n",
-    "\n",
-    "Now we can also have a look at the now-optimized initial velocity field. Here we see the y-components (from `data[...,0]`)."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "i7ZahlUudex8"
-   },
-   "outputs": [],
-   "source": [
-    "optimized_velocity_field = session.run(initial_state.velocity).at_centers()\n",
-    "\n",
-    "pylab.title('Initial y-velocity (optimized)')\n",
-    "pylab.imshow(np.concatenate(optimized_velocity_field.data[...,0], axis=1), origin='lower')"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "This nicely shows how the initial perturbation of the velocity field grows for the larger inflow deviations towards the right:\n",
-    "\n",
-    "For comlpeteness, here are the x-components:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "Pqw5BDxmdkut"
-   },
-   "outputs": [],
-   "source": [
-    "pylab.title('Initial x-velocity (optimized)')\n",
-    "pylab.imshow(np.concatenate(optimized_velocity_field.data[...,1], axis=1), origin='lower')"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "id": "ooqVxCPM8PXl"
-   },
-   "source": [
-    "This example illustrates how the differentiable physics approach can easily be extended towards significanlty more complex PDEs. Above, we've optimized for a mini-batch of 20 steps of a full Navier-Stokes solver.\n",
-    "\n",
-    "This is a powerful basis to bring NNs into the picture! As you might have noticed, our degrees of freedom were still a regular grid, and we've solved a single inverse problem. There was more than one case to solve, but nonetheless no question of generalization and no neural network structure."
+    "Based on this code example, we can recommend experimenting with the following:\n",
+    "- 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). This will make the optimization problem tougher, and might not converge directly with this simple setup.\n",
+    "- For more complex cases, add a multi-resolution optimization, that first solves with a coarse discretization, and then uses this solution as an initial guess for a finer discretization. (Ideally, implement it with a loop that continuously refines the discretization.)\n"
    ]
   },
   {

diffphys-discuss.md

@@ -1,12 +1,13 @@
 Discussion
 =======================
 
-**TODO** discussion of DP only
+The training via differentiable physics as described so far allows us
+to integrate full numerical simulations and the training of deep neural networks
+interacting with these simulations. While we've only hinted at what could be
+achieved via DP approaches it is nonetheless a good time to summarize the pros and cons.
 
 ## Summary
 
-Thus, to summarize ...
-
 ✅ Pro: 
 - uses physical model and numerical methods for discretization
 - efficiency of selected methods carries over to training
@@ -14,4 +15,6 @@ Thus, to summarize ...
 
 ❌ Con: 
 - not compatible with all simulators (need to provide gradients)
-- ...
+- require more heavy machinery (in terms of framework support) than previously discussed methods
+
+Especially the last point is one that is bound to strongly improve in a fairly short time, but for now it's important to keep in mind that not every simulator is suitable for DP training out of the box. Hence, in this book we'll focus on examples using phiflow, which was designed for interfacing with deep learning frameworks. Next we can target more a complex scenarios involving fluids.