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

smaller updates of DP part

Browse Source
This commit is contained in:
NT
2021-04-09 14:07:57 +08:00
parent 6d0651923e
commit 38da9fb388
4 changed files with 133 additions and 161 deletions
Download Patch File Download Diff File Expand all files Collapse all files

32
diffphys-code-ns.ipynb
View File

@@ -41,7 +41,7 @@
"source": [
"## Formulation\n",
"\n",
"As optimization objective we'll consider a more difficult variant of the previous example: the state of the observed density $d$ should match a given target after $n=20$ steps of simulation. In contrast to before, the marker $d$ cannot be modified in any way, but only the initial state of the velocity $\\mathbf{u}_0$ at $t=0$. This gives us a split between observable quantities for the loss formulation, and quantities that we can interact with during the optimization (or later on via NNs).\n",
"As optimization objective we'll consider a more difficult variant of the previous example: the state of the observed density $d$ should match a given target after $n=20$ steps of simulation. In contrast to before, the marker $d$ cannot be modified in any way, but only the initial state of the velocity $\\mathbf{u}_0$ at $t=0$. This gives us a split between observable quantities for the loss formulation and quantities that we can interact with during the optimization (or later on via NNs).\n",
"\n",
"With the notation from {doc}`overview-equations` this gives the minimization problem \n",
"\n",
@@ -70,7 +70,7 @@
},
{
"cell_type": "code",
"execution_count": 1,
"execution_count": null,
"metadata": {
"id": "da1uZcDXdVcF"
},
@@ -80,7 +80,7 @@
"#!pip install --upgrade --quiet git+https://github.com/tum-pbs/PhiFlow@develop\n",
"\n",
"from phi.torch.flow import * # The Dash GUI is not supported on Google Colab, ignore the warning\n",
"import pylab"
"import pylab\n"
]
},
{
@@ -240,7 +240,7 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 24,
"metadata": {},
"outputs": [
{
@@ -256,7 +256,7 @@
"(x=33, y=40) float32 -0.9229158759117126 < ... < 2.964700937271118"
]
},
"execution_count": 5,
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
@@ -266,7 +266,7 @@
"[velocity_grad] = sim_grad(initial_smoke, initial_velocity)\n",
"\n",
"print(\"Gradient info \" + format(velocity_grad))\n",
"velocity_grad.inflow_loc[0].vector[0].values # one example, location 0, x component\n"
"velocity_grad.inflow_loc[0].values.vector[0] # one example, location 0, x component\n"
]
},
{
@@ -280,14 +280,14 @@
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": 29,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CenteredGrid[(x=33, y=40), size=(33.0, 40.0) along vector, extrapolation=0]\n"
"x components: (x=33, y=40) float32 0.0 < ... < 0.0\n"
]
},
{
@@ -304,11 +304,11 @@
}
],
"source": [
"print(format( velocity_grad.inflow_loc[i].vector[0] ))\n",
"print(\"x components: \"+format( velocity_grad.inflow_loc[i].values.vector[0] ))\n",
"\n",
"fig, axes = pylab.subplots(1, 4, figsize=(16, 5))\n",
"for i in range(INFLOW.shape.inflow_loc):\n",
" im = axes[i].imshow(velocity_grad.inflow_loc[i].vector[0].values.numpy('y,x'), origin='lower', cmap='bwr')\n",
" im = axes[i].imshow(velocity_grad.inflow_loc[i].values.vector[0].numpy('y,x'), origin='lower', cmap='bwr')\n",
" pylab.colorbar(im, ax=axes[i])\n",
" axes[i].set_title(f\"Gradient for x vel - In {math.to_int(INFLOW_LOCATION).inflow_loc[i].vector.unstack_spatial('x,y')}\")\n",
"pylab.tight_layout()\n",
@@ -336,7 +336,7 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": 30,
"metadata": {},
"outputs": [
{
@@ -395,7 +395,7 @@
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": 32,
"metadata": {},
"outputs": [
{
@@ -414,7 +414,7 @@
"source": [
"fig, axes = pylab.subplots(1, 3, figsize=(12, 5))\n",
"for i in range(INFLOW.shape.inflow_loc-1):\n",
" im = axes[i].imshow(initial_velocity.inflow_loc[i].vector[0].values.numpy('y,x'), origin='lower', cmap='bwr')\n",
" im = axes[i].imshow(initial_velocity.inflow_loc[i].values.vector[0].numpy('y,x'), origin='lower', cmap='bwr')\n",
" pylab.colorbar(im, ax=axes[i])\n",
" axes[i].set_title(f\"Initial velocity x for {math.to_int(INFLOW_LOCATION).inflow_loc[i].vector.unstack_spatial('x,y')}\")\n",
"pylab.tight_layout()"
@@ -422,7 +422,7 @@
},
{
"cell_type": "code",
"execution_count": 9,
"execution_count": 33,
"metadata": {},
"outputs": [
{
@@ -441,7 +441,7 @@
"source": [
"fig, axes = pylab.subplots(1, 3, figsize=(12, 5))\n",
"for i in range(INFLOW.shape.inflow_loc-1):\n",
" im = axes[i].imshow(initial_velocity.inflow_loc[i].vector[1].values.numpy('y,x'), origin='lower', cmap='bwr')\n",
" im = axes[i].imshow(initial_velocity.inflow_loc[i].values.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()"
@@ -460,7 +460,7 @@
},
{
"cell_type": "code",
"execution_count": 10,
"execution_count": 34,
"metadata": {},
"outputs": [
{

18
diffphys-code-sol.ipynb
View File

@@ -239,11 +239,11 @@
"\n",
"In all cases we'll use fully convolutional networks, i.e. networks without any fully-connected layers. The\n",
"inputs are: \n",
"- 2 fields with x,y velociy\n",
"- 2 fields with x,y velocity\n",
"- plus the Reynolds number as constant channel.\n",
"\n",
"The output is : \n",
"- a 2 component field containing the x,y velocty\n",
"- a 2 component field containing the x,y velocity\n",
"\n",
"First, let's define a minimal network consisting only of three convolutional layers with ReLU activations (we're also using keras here for simplicity). The input channel dimension is defined via the `tensor_in`, then we'll go to 32 and 64 features, before reducing to 2 channels in the output. "
]
@@ -639,9 +639,9 @@
"source": [
"## Interleaving Simulation and Network\n",
"\n",
"Now comes the **most crucial** step in the whole setup: we define the chain of simulation steps and network evaluations to be used at training time. After all the work defining helper functions, it's acutally pretty simple: we loop over `msteps`, call the simulator via `KarmanFlow.step` for an input state, and afterwards evaluate the correction via `network(to_keras())`. The correction is then added to the last simultation state in the `prediction` list (we're actually simply overwriting the last simulated step `prediction[-1]` with `velocity + correction[-1]`.\n",
"Now comes the **most crucial** step in the whole setup: we define the chain of simulation steps and network evaluations to be used at training time. After all the work defining helper functions, it's acutally pretty simple: we loop over `msteps`, call the simulator via `KarmanFlow.step` for an input state, and afterwards evaluate the correction via `network(to_keras())`. The correction is then added to the last simulation state in the `prediction` list (we're actually simply overwriting the last simulated step `prediction[-1]` with `velocity + correction[-1]`.\n",
"\n",
"One other important things that's happening here is normalization: the inputs to the network are divided by the standard deviations in `dataset.dataStats`. This is slightly complicated as we have to append the scaling for the Reynolds numbers to the normalization for the velocity. After evaluating the `network`, we only have a velocity left, so we can simply multiply by the standard deviation again (`* dataset.dataStats['std'][1]`)."
"One other important thing that's happening here is normalization: the inputs to the network are divided by the standard deviations in `dataset.dataStats`. This is slightly complicated as we have to append the scaling for the Reynolds numbers to the normalization for the velocity. After evaluating the `network`, we only have a velocity left, so we can simply multiply by the standard deviation again (`* dataset.dataStats['std'][1]`)."
]
},
{
@@ -882,7 +882,7 @@
"\n",
"Let's see how well this works by applying it to a set of test data inputs with new Reynolds numbers that were not part of the training data.\n",
"\n",
"To keep things somewhat simple, we won't aim for a high performance version of our hybrid solver. For performance, please checck out the external code base: the network trained here should be directly useable in [this apply script](https://github.com/tum-pbs/Solver-in-the-Loop/blob/master/karman-2d/karman_apply.py).\n",
"To keep things somewhat simple, we won't aim for a high-performance version of our hybrid solver. For performance, please check out the external code base: the network trained here should be directly useable in [this apply script](https://github.com/tum-pbs/Solver-in-the-Loop/blob/master/karman-2d/karman_apply.py).\n",
"\n",
"---"
]
@@ -1112,9 +1112,9 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"They both start out with the same initial state at $t=0$ (a downsampled solution from the reference solution manifold), and at $t=20$ the solutions still share similarities. Over time, the source version strongly diffuses the strucutres in the flow and looses momentum. The flow behind the obstacles becomes straight, and lacks clear vortices. \n",
"They both start out with the same initial state at $t=0$ (a downsampled solution from the reference solution manifold), and at $t=20$ the solutions still share similarities. Over time, the source version strongly diffuses the structures in the flow and looses momentum. The flow behind the obstacles becomes straight, and lacks clear vortices. \n",
"\n",
"The version produced by the hybrid solver does much better. It preserves the vortex shedding even after more than one hundred updates. Note that both outputs were produced by the same underlying solver. The second version just profits from the learned corrector which has learned to revert the overly strong dissipation of the solver. However, it also produces some visible axis aligned structures, especially near the sides of the domain. These could be alleviated with improved training setups, e.g., more look-ahead, and a larger model.\n",
"The version produced by the hybrid solver does much better. It preserves the vortex shedding even after more than one hundred updates. Note that both outputs were produced by the same underlying solver. The second version just profits from the learned corrector which has learned to revert the overly strong dissipation of the solver. However, it also produces some visible axis-aligned structures, especially near the sides of the domain. These could be alleviated with improved training setups, e.g., more look-ahead, and a larger model.\n",
"\n",
"We can also compare and quantify how the models do w.r.t. reference data. The next cell plots one time step of the three versions: the reference data after 50 steps, and the re-simulated version of the source and our hybrid solver, together with a per-cell error of the two:"
]
@@ -1172,7 +1172,7 @@
"source": [
"This shows very clearly how the pure source simulation in the middle deviates from the reference on the left. Apart from the slight vertical streaks, the learned version stays much closer to the reference solution. \n",
"\n",
"The two per-cell error images on the right also illustrate this: the source version has much larger errors that show how it systematically underestimates the vortices that should form. The error for the learned version is much more evenly distributed, and significantly smaller.\n",
"The two per-cell error images on the right also illustrate this: the source version has much larger errors that show how it systematically underestimates the vortices that should form. The error for the learned version is much more evenly distributed and significantly smaller.\n",
"\n",
"Next, we can quantify these observations:"
]
@@ -1226,7 +1226,7 @@
"source": [
"The overall mean absolute error (MAE) is almost 50% larger for the regular simulation. The graph below shows this behavior over time: the error of the source version is clearly above the errors of the hybrid simulator.\n",
"\n",
"This concludes our evaluation. Note that the improved behavior of the hybrid solver is difficult to reliably measure with simple vector norms such as an MAE or $L^2$ norm. To achieve this, we'd need to employ other, domain specific metrics. In this case, metrics for fluids based on vorticity and turbulence properties of the flow would be applicable. However, in this text we instead want to focus on DL-related topics, and target another inverse problem with differentiable physics solvers in the next chapter."
"This concludes our evaluation. Note that the improved behavior of the hybrid solver is difficult to reliably measure with simple vector norms such as an MAE or $L^2$ norm. To achieve this, we'd need to employ other, domain-specific metrics. In this case, metrics for fluids based on vorticity and turbulence properties of the flow would be applicable. However, in this text, we instead want to focus on DL-related topics and target another inverse problem with differentiable physics solvers in the next chapter."
]
},
{

242
diffphys-control.ipynb
View File

File diff suppressed because one or more lines are too long

2
intro.md
View File

@@ -82,10 +82,10 @@ See also... Test link: {doc}`supervised`
## TODOs , include
- DP SoL, generate result
- DP control, show targets at bottom?
- update teaser image
- DP intro, check transpose of Jacobians in equations
- fix phiflow2 , diffphys-code-ns.ipynb
- outlook: include latent space physics & LSTMs, also ContConv benjamin (see below)