update dp burgers for phiflow2
This commit is contained in:
NT committed
1 parent
e99b2a0b4a
commit
feb1477391
7 files changed
+783
-151
No files matched your search
@@ -19,7 +19,6 @@
|
||||
- file: diffphys
|
||||
sections:
|
||||
- file: diffphys-code-gradient.ipynb
|
||||
- file: diffphys-code-tf.ipynb
|
||||
- file: diffphys-discuss.md
|
||||
- file: diffphys-code-ns-v2a.ipynb
|
||||
- file: diffphys-code-sol.ipynb
|
||||
@@ -29,6 +28,8 @@
|
||||
- file: overview-burgers-forw-v1.ipynb
|
||||
- file: overview-ns-forw-v1.ipynb
|
||||
- file: diffphys-code-ns.ipynb
|
||||
- file: diffphys-code-gradient-v1.ipynb
|
||||
- file: diffphys-code-tf.ipynb
|
||||
- file: jupyter-book-reference
|
||||
sections:
|
||||
- file: jupyter-book-reference-markdown
|
||||
|
||||
@@ -0,0 +1,285 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Burgers Optimization with a (Manual) Differentiable Physics Gradient\n",
|
||||
"\n",
|
||||
"To illustrate the process of compute gradients in a differentiable physics setting, we first demonstrate a process that isn't recommended for more complex real-world cases, but more clearly shows the way we can get gradient information from a physical simulation in a DL framework. We'll target tensorflow in the following.\n",
|
||||
"\n",
|
||||
"We target the same reconstruction task like for the PINN example. This has some immediate implications: the evolution of the system is now fully determined by our PDE formulatio. Hence, the only real unknown is the initial state! We will still need to re-compute all the states betwwen the initial and target state many times, just now we won't need an NN for this step, we can rely on our discretized model. Also, as we choose an initial discretization for the DP approach, the unknown initial state consists of the sampling points of the involved physical fields, and we can simply represent these unknowns as floating point variables. Hence, even for the initial state we do not need to set up an NN. Thus, our Burgers reconstruction problem reduces to a gradient-based opitmization without any NN when solving it with DP. Nonetheless, it's a very good starting point to illustrate the process.\n",
|
||||
"\n",
|
||||
"First, we'll set up our discretized simulation. Here we can employ the phiflow solver, as shown in the overview section on _Burgers forward simulations_. phiflow directly gives us a computational graph in tensorflow. So, as a first step, let's set up our grid with 128 points, and construct a tensorflow graph for 32 steps of the simulation. (As usual, you can ignore the warnings when loading TF and phiflow.)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Each velocity is a phiflow grid like this one: Grid[128(1), size=[2.], <class 'numpy.float32'>]\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"from phi.tf.flow import *\n",
|
||||
"\n",
|
||||
"# run with phiflow\n",
|
||||
"n = 128\n",
|
||||
"dt = 1./32.\n",
|
||||
"initial = np.zeros([n,1]) # start from 0\n",
|
||||
"\n",
|
||||
"domain = Domain([n], boundaries=PERIODIC, box=box[-1:1])\n",
|
||||
"state0 = BurgersVelocity(domain, velocity=initial, viscosity=0.01/np.pi)\n",
|
||||
"\n",
|
||||
"# start with a state0, to be modified\n",
|
||||
"state_in = state0.copied_with(velocity=placeholder)\n",
|
||||
"states = [state_in]\n",
|
||||
"\n",
|
||||
"for i in range(32):\n",
|
||||
" states.append( Burgers().step(states[-1],dt=dt) )\n",
|
||||
"\n",
|
||||
"print(\"Note: each velocity is a phiflow grid like this one \" + format(states[-1].velocity) )"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"**A lot** has happened in these first steps, but we didnt run a single calculation for our simulation. In contrast to the first example, which used numpy as phiflow backend, we now have used the tensorflow backend (the trick here was importing `phi.tf.flow` instad of `phi.flow`). This gave us a tensorflow graph, but didn't run any real simulation steps. That will only happen once we run a tensorflow session and pass some data through this graph. So, in a way _nothing_ has happened so far in terms of our simulation! \n",
|
||||
"\n",
|
||||
"Next, let's set up the loss, such that we can start the optimization. That's actually pretty simple, we want the solution at $t=0.5$, i.e. step number 16, to respect the data in terms of an L2 norm. Let's also directly initialize the TF variables, and see what loss we get from the initial zero'ed initialization:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/session.py:17: The name tf.Session is deprecated. Please use tf.compat.v1.Session instead.\n",
|
||||
"\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/session.py:18: The name tf.get_default_graph is deprecated. Please use tf.compat.v1.get_default_graph instead.\n",
|
||||
"\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/session.py:28: The name tf.global_variables_initializer is deprecated. Please use tf.compat.v1.global_variables_initializer instead.\n",
|
||||
"\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/session.py:29: The name tf.train.Saver is deprecated. Please use tf.compat.v1.train.Saver instead.\n",
|
||||
"\n",
|
||||
"Initial loss 0.3829152584075928\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"sess = Session(None) \n",
|
||||
"\n",
|
||||
"# enforce constraints at time t=0.5 (step 16)\n",
|
||||
"solution_t16 = np.asarray( [0.008612174447657694, 0.02584669669548606, 0.043136357266407785, 0.060491074685516746, 0.07793926183951633, 0.0954779141740818, 0.11311894389663882, 0.1308497114054023, 0.14867023658641343, 0.1665634396808965, 0.18452263429574314, 0.20253084411376132, 0.22057828799835133, 0.23865132431365316, 0.25673879161339097, 0.27483167307082423, 0.2929182325574904, 0.3109944766354339, 0.3290477753208284, 0.34707880794585116, 0.36507311960102307, 0.38303584302507954, 0.40094962955534186, 0.4188235294008765, 0.4366357052408043, 0.45439856841363885, 0.4720845505219581, 0.4897081943759776, 0.5072391070000235, 0.5247011051514834, 0.542067187709797, 0.5593576751669057, 0.5765465453632126, 0.5936507311857876, 0.6106452944663003, 0.6275435911624945, 0.6443221318186165, 0.6609900633731869, 0.67752574922899, 0.6939334022562877, 0.7101938106059631, 0.7263049537163667, 0.7422506131457406, 0.7580207366534812, 0.7736033721649875, 0.7889776974379873, 0.8041371279965555, 0.8190465276590387, 0.8337064887158392, 0.8480617965162781, 0.8621229412131242, 0.8758057344502199, 0.8891341984763013, 0.9019806505391214, 0.9143881632159129, 0.9261597966464793, 0.9373647624856912, 0.9476871303793314, 0.9572273019669029, 0.9654367940878237, 0.9724097482283165, 0.9767381835635638, 0.9669484658390122, 0.659083299684951, -0.659083180712816, -0.9669485121167052, -0.9767382069792288, -0.9724097635533602, -0.9654367970450167, -0.9572273263645859, -0.9476871280825523, -0.9373647681120841, -0.9261598056102645, -0.9143881718456056, -0.9019807055316369, -0.8891341634240081, -0.8758057205293912, -0.8621229450911845, -0.8480618138204272, -0.833706571569058, -0.8190466131476127, -0.8041372124868691, -0.7889777195422356, -0.7736033858767385, -0.758020740007683, -0.7422507481169578, -0.7263049162371344, -0.7101938950789042, -0.6939334061553678, -0.677525822052029, -0.6609901538934517, -0.6443222327338847, -0.6275436932970322, -0.6106454472814152, -0.5936507836778451, -0.5765466491708988, -0.5593578078967361, -0.5420672759411125, -0.5247011730988912, -0.5072391580614087, -0.4897082914472909, -0.47208460952428394, -0.4543985995006753, -0.4366355580500639, -0.41882350871539187, -0.40094955631843376, -0.38303594105786365, -0.36507302109186685, -0.3470786936847069, -0.3290476440540586, -0.31099441589505206, -0.2929180880304103, -0.27483158663081614, -0.2567388003912687, -0.2386513127155433, -0.22057831776499126, -0.20253089403524566, -0.18452269630486776, -0.1665634500729787, -0.14867027528284874, -0.13084990929476334, -0.1131191325854089, -0.09547794429803691, -0.07793928430794522, -0.06049114408297565, -0.0431364527809777, -0.025846763281087953, -0.00861212501518312] )\n",
|
||||
"target = state0.copied_with(velocity=np.reshape(solution_t16,[n,1]))\n",
|
||||
"\n",
|
||||
"loss = math.sum( (states[16].velocity.data - target.velocity.data)**2 ) / n\n",
|
||||
"sess.initialize_variables()\n",
|
||||
"\n",
|
||||
"# compute initial loss\n",
|
||||
"s,l = sess.run([states[-1],loss], feed_dict={state_in: state0})\n",
|
||||
"print(\"Initial loss \"+format(l))\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Not surprisingly, because we're starting from zero, there's a significant intial error of ca. 0.38 for the 16th simulation step.\n",
|
||||
"\n",
|
||||
"Note that because we're only constraining timestep 16, we could actually omit steps 17 to 31 in this setup. They don't have any degrees of freedom and are not constrained in any way. For fairness and for direct comparison with the previous PINN case, let's include them.\n",
|
||||
"\n",
|
||||
"Now we have our simulation graph in TF, we can use TF to give us a gradient for the initial state for the loss. All we need to do is run `tf.gradients(loss, [state_in.velocity.data]`, which will give us a direction for each velocity variable. Based on a linear approximation, this direction tells us how to change each of them to increase the loss function (gradients _always_ point \"upwards\"). In the following code snipped, we're additionally saving all these gradients in a list called `grads`, such that we can visualize them later on. (Normally, we could discard each gradient after performing an update step.)\n",
|
||||
"\n",
|
||||
"Based on the gradient, we can now take a step in the opposite direction to bring the loss down (instead of increasing it). Below we're using a learning rate `LR=5` for this step. Afterwards, we're re-evaluating the loss for the updated state to check how we did. As these updates aren't exactly fast, we're only computing five of them here:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"WARNING:tensorflow:From /usr/local/lib/python3.7/site-packages/tensorflow_core/python/ops/math_grad.py:1375: where (from tensorflow.python.ops.array_ops) is deprecated and will be removed in a future version.\n",
|
||||
"Instructions for updating:\n",
|
||||
"Use tf.where in 2.0, which has the same broadcast rule as np.where\n",
|
||||
"Step 0, loss = 0.3829152584075928\n",
|
||||
"Step 1, loss = 0.3334265947341919\n",
|
||||
"Step 2, loss = 0.2905231714248657\n",
|
||||
"Step 3, loss = 0.25346457958221436\n",
|
||||
"Step 4, loss = 0.2215319573879242\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"LR = 5.\n",
|
||||
"stateN = state0\n",
|
||||
"grads = []\n",
|
||||
"for i in range(5):\n",
|
||||
" grads.append( tf.gradients(loss, [state_in.velocity.data] ) ) # use TF\n",
|
||||
" grads.append( gradients(loss, state_in.velocity) ) # phiflow wrapper also works\n",
|
||||
"\n",
|
||||
" state_updated = state_in.copied_with( velocity=(state_in.velocity - LR*grads[-1]) )\n",
|
||||
" stateN,l = sess.run([state_updated,loss], feed_dict={state_in: stateN})\n",
|
||||
"\n",
|
||||
" #s,l = sess.run([states[-1],loss], feed_dict={state_in: stateN})\n",
|
||||
" print( \"Step {}, loss = {}\".format(i,l) )\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"One of the reasons for the bad performance here is that we're querying states via `sess.run()` multiple times, instead of computing the update in one go. Nonetheless, this simple example already clearly shows that the loss is nicely going down: from above 0.38 to ca. 0.22. This is a good sign!\n",
|
||||
"\n",
|
||||
"Before improving on the performance, let's visualize the reconstruction.\n",
|
||||
"\n",
|
||||
"Starting from the zero state `state0` (shown in blue), the first gradient is shown as a red line here. If you compare it with the solution it points in the opposite direction, as expected. The solution is much larger in magnitude, so we omit it here (see the next graph).\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[<matplotlib.lines.Line2D at 0x13fa0ae90>]"
|
||||
]
|
||||
},
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAD4CAYAAAAUymoqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de9zUc/r48dfV3UkpKreWQiFs2E3uDYtO0kkqhApltVrUOsR3ZX2t0+5+sb+VczY5xC6VsN1YkpKz6o6UIt2FVZu6VXIu5fr9cb1nG7e5T83c85nD9Xw85jGf81wz3c01n/dRVBXnnHOuJupEHYBzzrns48nDOedcjXnycM45V2OePJxzztWYJw/nnHM1VjfqANJht9120zZt2kQdhnPOZZUFCxZ8qqqFifblRfJo06YNJSUlUYfhnHNZRUQ+qmifF1s555yrsZQkDxHpLSLLRKRURMYm2N9ARKaE/XNFpE3Y3kJEXhCRL0XkjnLnHC4ii8M5t4mIhO3NRWSmiCwPz81S8R6cc85VX9LJQ0QKgDuBPkB7YIiItC932Ahgo6ruD4wDbgzbvwWuAi5LcOnxwLlAu/DoHbaPBWapajtgVlh3zjmXRqm48+gElKrqSlXdAkwGBpQ7ZgAwKSxPA44TEVHVr1T1FSyJ/JeI7AE0VdU31MZPeRAYmOBak+K2O+ecS5NUJI9WwMdx66vCtoTHqOpWYBPQooprrqrgmi1VdU1Y/gRomegCIjJSREpEpKSsrKw678M551w1ZXWFebgrSTiyo6pOUNUiVS0qLEzY0sw559wOSkXyWA3sFbfeOmxLeIyI1AV2AdZXcc3WFVxzbSjWihVvrdvhyJ1zzu2QVPTzmA+0E5G22Bf8YGBouWOKgeHA68AgYLZWMha8qq4Rkc9F5EhgLjAMuL3ctW4Iz9NT8B6cc+Vt2QIvvQQvvwxz58L330OjRtC6NXTrBl27QovKSp9dLks6eajqVhEZDcwACoD7VHWJiFwHlKhqMXAv8JCIlAIbsAQDgIh8CDQF6ovIQKCnqi4FLgAeAHYCngkPsKQxVURGAB8BpyX7HpxzcVThqafgkktgxYrEx9x5J9SpA/36wejRcNxxtu7yhuTDZFBFRUXqPcydq4YvvoAzz4TiYlvff3848UQ4+mjYeWf46itYuhRmz4ZXXoHvvrPjDj0UbrgB+vQB65LlcoCILFDVooT7PHk45wD4/HPo2xdefRWaNoVrr4VRo6BevcTHr10LEyfC+PGwOlRJdukCt99uycRlvcqSh99nOucscfTubYljr73gzTfh4osrThwALVvClVdCaSn89a/QvDm8+CJ07AiXX253KS5nefJwLt+pwogR8PrrsPfeMGcO7Ldf9c9v2BDGjLH6kVGjYNs2uOkmOOwwmD+/1sJ20fLk4Vy+e+ghmDYNmjSBWbNg33137Dq77gp33GEtsw45BJYvh6OOgj/+EbZuTW3MLnKePJzLZx9+aK2lwOoq9t8/+Wv+4hd2x3HJJXYXctVVVheycmXy13YZw5OHc/lKFc4+21pYnXIKDBuWums3bAg33wwzZ0KrVvDaa/Dzn8Mjj6TuNVykPHk4l6+Ki62Cu7AQ/va32mli26MHLFoEp54KX34JQ4fChRdaB0SX1Tx5OJePtm6FK66w5T/8oXZ7ijdvDlOmwF13Weut22+33umrVlV5qstcnjycy0eTJsG771rl+MiRtf96InD++TbUSevW1rKrY0frbOiykicP5/LNN9/A1Vfb8h//CPXrp++1jzjC+pD06AFlZXD88XDbbVb/4rKKJw/n8s1991mP8A4d4PTT0//6hYXw7LPw+9/bYIsXXQTnnbd9qBOXFTx5OJdPtm2DW26x5SuvjG4ww4IC+NOf4OGHoUEDmDABevaE9ZXN1OAyiScP5/LJU0/ZcCJt2sDADJjBecgQG/b9Jz+xnu1HHGF1MS7jefJwLp+MG2fPF10EdVMxnU8KdOpknQo7drQhTo480oq1XEbz5OFcvliwwPp1NGkC55wTdTQ/1Lq13YEMGmSDNJ5wAtx6q1ekZzBPHs7li1tvtedzz7Uh1zNN48bWH+Sqq6wi/eKL4be/tXoal3E8eTiXDzZuhEcfteVRo6KNpTJ16sB111lFev36NmPhKafA119HHZkrJyXJQ0R6i8gyESkVkbEJ9jcQkSlh/1wRaRO374qwfZmI9ArbDhSRhXGPz0Xk4rDvGhFZHbevbyreg3M57eGH4dtvrX/Fjo6am05Dhti4WLvuCtOn2zS3n34adVQuTtLJQ0QKgDuBPkB7YIiItC932Ahgo6ruD4wDbgzntsfmMz8Y6A3cJSIFqrpMVTuoagfgcOBr4Im4642L7VfVfyX7HpzLeffea8+//nW0cdRE5842OdXee8Mbb8Avf1nxnOou7VJx59EJKFXVlaq6BZgMDCh3zABgUlieBhwnIhK2T1bVzar6AVAarhfvOGCFqn6Uglidyz9vvglvvWVjTGVC89yaaN/eEsdhh22fH2TevKijcqQmebQCPo5bXxW2JTxGVbcCm4AW1Tx3MFB+HOfRIrJIRO4TkWaJghKRkSJSIiIlZWVlNXk/zuWW2F3HWWdZh7xss8ce1kqsVy8b0qRrV3jyyaijynsZXWEuIvWB/sCjcZvHA/sBHYA1wF8TnauqE1S1SFWLCgsLaz1W5zLSN9/AP/5hyyNGRBtLMpo0sYTxq1/Zexo4EO65J+qo8loqksdqYK+49dZhW8JjRKQusAuwvhrn9gHeVNW1sQ2qulZVt6nq98A9/LiYyzkX8+STsGkTFBXBoYdGHU1y6tWzu6g//MGa8o4cCf/3f94XJCKpSB7zgXYi0jbcKQwGissdUwwMD8uDgNmqqmH74NAaqy3QDogv0BxCuSIrEdkjbvUk4J0UvAfnclPsruPMM6ONI1VE4NprbW4QERtc8dJLLZm4tEp6fAJV3Soio4EZQAFwn6ouEZHrgBJVLQbuBR4SkVJgA5ZgCMdNBZYCW4FRqroNQEQaA8cDvyn3kjeJSAdAgQ8T7HfOgQ0y+Mwz1nciitFza9P559sEVmeeaUOufPqp3ZXUqxd1ZHlDNA9u+YqKirSkpCTqMJxLr7vvti/ZXr1yd6yomTPhpJPgq6+gXz/rod6oUdRR5QwRWaCqRYn2ZXSFuXMuCbEiqzPOiDaO2nT88TBrljVDfuopS5SffRZ1VHnBk4dzuejDD+GVV2CnnbKvb0dNHXGEvdfWre25SxdYsybqqHKeJw/nctEjoZ3JgAHWzDXX/fSn1hv9wANh0SI45hhLoK7WePJwLhdNmWLPQ4ZEG0c67b03vPwyHH44rFxpw5u8/37UUeUsTx7O5Zply+Dtt2GXXawOIJ8UFlodyNFHw8cfWwJZvDjqqHKSJw/nck3srmPgwOwcjiRZu+wCM2bYCMJr19pwJt7aMuU8eTiXa6ZOtedc69tRE40bW+/6E0+EDRuge3erTHcp48nDuVyyZIk9mjWzOTDyWcOG8NhjlkS/+AJ69rR+IS4lPHk4l0tidx0nn2wz8eW7evWsv0tsQMV+/aC4/OhJbkd48nAuV6hur+/I5yKr8goKYOJEmw99yxZLrJMnRx1V1vPk4VyuWLTIWlrttht06xZ1NJmlTh249VYYOxa2bYOhQ+G++6KOKqt58nAuV8SKrE45BeomPeZp7hGxIdz/9Ce7SxsxAu68M+qospYnD+dygRdZVd/vf28j8QKMHm13JK7GPHk4lwvefBNWrICWLa1jnKvcxRfDHXdsX7755mjjyUKePJzLBbEiq0GDrILYVW3UKPjb32z50kvhppuijSfLePJwLtupesfAHTVypE0iJQKXXw5//nPUEWUNTx7OZbt582wE2T33tDGdXM2ccw7cf78lkCuvhOuuizqirJCS5CEivUVkmYiUisjYBPsbiMiUsH+uiLSJ23dF2L5MRHrFbf9QRBaLyEIRKYnb3lxEZorI8vDcLBXvwbmsFRt+/bTTrEmqq7nhw+HBB+3zu/pq+MMf7I7OVSjpvzQRKQDuBPoA7YEhItK+3GEjgI2quj8wDrgxnNsem8/8YKA3cFe4Xkw3Ve1QbhrEscAsVW0HzArrzuWnbdu2d3gbOjTaWLLdmWdab/SCArj+ersL8QRSoVT8TOkElKrqSlXdAkwGBpQ7ZgAwKSxPA44TEQnbJ6vqZlX9ACgN16tM/LUmATk+TZpzlXjhBRs5dv/9oSjhVNOuJgYPtmRct671Cbn8ck8gFUhF8mgFfBy3vipsS3iMqm4FNgEtqjhXgedEZIGIjIw7pqWqxuaY/ARomSgoERkpIiUiUlJWVlbzd+VcNnj4YXs+4wwrs3fJGzTIGiDUrQt/+Yu1xPIE8iOZXEB6jKp2xIrDRonIjxqvq6piSeZHVHWCqhapalFhYWEth+pcBL791kaNhfyaMTAdTjrJPtt69axD4WWXeQIpJxXJYzWwV9x667At4TEiUhfYBVhf2bmqGnteBzzB9uKstSKyR7jWHsC6FLwH57LP00/D55/btKsHHhh1NLmnf394/HFLIDffDL/7nSeQOKlIHvOBdiLSVkTqYxXg5cc8LgaGh+VBwOxw11AMDA6tsdoC7YB5ItJYRJoAiEhjoCfwToJrDQemp+A9OJd9YkVWXlFee/r1g2nTLIH8v/9nAyt6AgFSkDxCHcZoYAbwLjBVVZeIyHUi0j8cdi/QQkRKgTGEFlKqugSYCiwFngVGqeo2rB7jFRF5G5gHPK2qz4Zr3QAcLyLLgR5h3bn88tlnduch4h0Da1v//tvrQG66ycbG8gSCaB58CEVFRVricxi7XHL//da5rVs3mD076mjyw+OPW1+abdusGe/11+d8IwURWVCuq8R/ZXKFuXOuIl5klX6xSaQKCmxY92uuiTqiSHnycC7brFljdxv16tncHS59Bg2yHv0FBTaMybXXRh1RZDx5OJdLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 432x288 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"\n",
|
||||
"fig = plt.figure().gca()\n",
|
||||
"pltx = np.linspace(-1,1,n)\n",
|
||||
"\n",
|
||||
"# first gradient\n",
|
||||
"[grad0,s0] = sess.run([grads[0],state0], feed_dict={state_in: stateN}) #.velocity.data\n",
|
||||
"fig.plot(pltx, grad0[0].flatten() , lw=2, color='red') \n",
|
||||
"\n",
|
||||
"fig.plot(pltx, state0.velocity.data.flatten(), lw=2, color='mediumblue')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"So the gradient looks good (feel free to plot the others from the five iterations above).\n",
|
||||
"\n",
|
||||
"How well does the 16th state of the simulation actually match the target after the 5 update steps? This is what the loss measures, after all. The next graph shows the constraints (i.e. the solution we'd like to obtain) in green, and the reconstructed state after the initial state was updated 16 times by the solver:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[<matplotlib.lines.Line2D at 0x13fbd43d0>]"
|
||||
]
|
||||
},
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3xUVfr48c+T3miht9CbSyeEAAoiqAio+911LbsqKqxlpQnSxIrKoihNXRWVddfefruLgAIqWCCUIFW61FBCCKGkZ2bO74+ZxElIaHOTm0me9+uVV2ZumXlyZ3Kfe8655xwxxqCUUqryCrA7AKWUUvbSRKCUUpWcJgKllKrkNBEopVQlp4lAKaUquSC7A7gctWrVMk2bNrU7DKWU8ivr168/YYypXXS5XyaCpk2bkpiYaHcYSinlV0TkQHHLtWpIKaUqOU0ESilVyWkiUEqpSs6SRCAi80XkuIhsLWG9iMhcEdkjIptFpKvXuqEistvzM9SKeJRSSl08q0oE7wIDz7P+BqCV5+d+4HUAEYkGngJ6AHHAUyJSw6KYlFJKXQRLEoEx5gfg5Hk2uRn4t3FbDVQXkfrA9cAyY8xJY0wasIzzJxSllFIWK6s2gobAIa/nSZ5lJS0/h4jcLyKJIpKYkpJSaoEqpVRl4zf9CIwx84B5ALGxsTp2tipTuc5c1h9ez5ZjW9iRsoMcRw43XXET/Vv0JyjAb/6NlCpWWX2DDwONvZ438iw7DFxdZPmKMopJqfPKyM3g21+/5etdX/PT/p/IyMsotP7r3V9TL6oeMwbNoFdML5uiVMp3ZZUIFgAjRORj3A3Dp40xR0VkCTDNq4H4OmByGcWk1DmMMaxLWsfHmz9m6e6lZDmyCta1rtWaHo160LZOWzJzM/lg0wfsT9vPJ5s/0USg/JoliUBEPsJ9ZV9LRJJw3wkUDGCMeQNYDAwC9gCZwL2edSdF5FlgneelphpjztforFSpyHHk8J9t/+Hd9e+yO3V3wfIuDbowpO0Qrm1xLQ2rFW6+almzJfd+cS9pWWllHa5SlrIkERhj7rjAegM8XMK6+cB8K+JQ6lKdzDzJB5s+4L0N75GamQpA7cja/KnDn7i1/a00rt64xH2rh1cH0ESg/J62cqlK6djZY7yx9g0+2/IZ2Y5sAH5X93cMjx3ODa1vIDgw+IKvER0eDWgiUP5PE4GqVFIyUnhjzRt8uOlDcp25APRt1pe/dv8r8Y3jEZGLfq0a4e6mLU0Eyt9pIlCVQmpmKvPWzuP9je8XlAAGtR7EiJ4jaFO7zWW9ZkRwBCGBIWQ7ssnKyyI8ONzKkJUqM5oIVIWW48jhn+v/yetrXic9Nx2AAS0HMKbXGNrVaefTa4sINcJrkJyeTFpWmiYC5bc0EagKyRjDwp0LmfHDDA6fOQzAVU2vYuyVY+lYr6Nl7+OdCBpUbWDZ6ypVljQRqApn/eH1TFsxjY1HNwLQplYbJl89mauaXmX5e2k7gaoINBGoCuPgqYPM+GEGi3ctBqBWRC3GXjmWW9rfQmBAYKm8Z40wdyI4maXdX5T/0kSg/F52XjZvrH2DN9e+Sa4zl9CgUIbHDuf+uPuJCokq1ffOLxGcyjpVqu+jVGnSRKD82vK9y5n67VQOnj4IwM3tbubRqx4ts/p6rRpSFYEmAuWXDp8+zLPLn2XZnmUAtKrZiqkDphLXOK5M49Dexaoi0ESg/EquM5f5ifN5JeEVsh3ZRAZHMqrXKIZ2HXpRvYGtVtC7OFsTgfJfmgiU31h7aC2PL3ucX0/+CsCgNoOYcvUU6lWpZ1tMWjWkKgJNBKrcO5tzlhd+eIGPNn0EQNMaTXmm/zNc2fRKmyPTqiFVMWgiUOXakt1LePqbpzmecZzggGAe7PEgD/V4iNCgULtDA3TgOVUxaCJQ5VJyejLPfPsMS3YvAaBL/S5Mu34arWu1tjmywrRqSFUEmghUueIyLj7d8inTv5/O2ZyzRAZH8uhVj/KXzn8ptU5hvtCB51RFYNUMZQOBOUAg8LYxZnqR9bOAfp6nEUAdY0x1zzonsMWz7qAx5iYrYlL+J+l0EhO/nsjqQ6sBuKb5NTwz4JlyPYaPDjynKgKfE4GIBAKvAdcCScA6EVlgjNmWv40x5hGv7UcCXbxeIssY09nXOJT/Msbw0aaPmP79dDLyMogOj+bJ/k8ypM2QS5ofwC468Jzyd1aUCOKAPcaYvQCeCepvBraVsP0duOc0VorDpw8zeelkVh5YCbjnCHh6wNPUjKhpb2CXQNsJlL+zIhE0BA55PU8CehS3oYg0AZoB33ktDhORRMABTDfG/LeEfe8H7geIiYmxIGxlJ2MMH2/+mOnfTyc9N50a4TV4pv8zDG472O7QLpkOPKf8XVk3Ft8OfG6McXota2KMOSwizYHvRGSLMebXojsaY+YB8wBiY2NN2YSrSsORM0eYvGQyPx34CYCBrQbyzIBnqBVZy+bILo8OPKf8nRWJ4DDQ2Ot5I8+y4twOPOy9wBhz2PN7r4iswN1+cE4iUP7PGMNnWz/j+eXPk56bTvWw6jw94Gm/aQsoiVYNKX9nRSJYB7QSkWa4E8DtwJ+LbiQibYEaQILXshpApjEmR0RqAb2BFy2ISZUzyenJTPp6Ej/s/wGA61pdx7MDnvXbUoA37V2s/J3PicAY4xCREcAS3LePzjfG/CIiU4FEY8wCz6a3Ax8bY7yrddoBb4qICwjA3UZQUiOz8lOLdy7miWVPcCr7FNXCqvF0/6e5se2Nfl0K8KYDzyl/Z0kbgTFmMbC4yLInizx/upj9VgEdrIhBlT9nss/wzHfP8N9t7vb/Pk378MLAF6gTVcfmyKylVUPK32nPYlUqEg4mMP6r8Rw9e5SwoDAeu/ox/tzpzxWmFOBNE4Hyd5oIlKVyHDm89ONLzF8/H4CO9Toyc9BMmkU3szmy0qOJQPk7TQTKMtuOb2PsorHsTt1NoAQyoucIHurxkC0TxpQlTQTK32kiUD5zupy8ve5tZq2cRZ4rj2Y1mjFz0Ew61u9od2hlIiI4gqCAILId2eQ4csrNENlKXSxNBMonR84cYdzicaxNWgvAXZ3vYmLfiZVq8DURITggGIfLgcPlIBRNBMq/aCJQl23xzsVMWTqFMzlnqB1ZmxcGvkDfZn3tDssWQYFB4ACHy2F3KEpdMk0E6pJl5Gbw7HfP8tnWzwDo36I/f7/+7341UJzVggLc/0qaCJQ/0kSgLsnmY5t5ZNEj7E/bT2hQKI/1fYy/dP5Lhbwt9FJoIlD+TBOBuigu4+KtdW8x86eZOFwO2tRqw+whs8vd1JF2KUgETk0Eyv9oIlAXdOzsMcZ/NZ5VB1cBMLTrUCb2mah3x3gpSARGE4HyP5oI1Hkt27OMyUsmk5aVRnR4NC/e8CL9mve78I6VjJYIlD/TRKCKlZWXxbQV0/hw04eAe5ygF294kdqRtW2OrHzSNgLlzzQRqHPsSNnB6C9Hs+fkHkICQ5jQZwJDuw4lQALsDq3c0kSg/JkmAlXAGMOHmz7kueXPkevMpWV0S2YPmU27Ou3sDq3cy08ETpfzAlsqVf5oIlCAe8joyUsm8/XurwG4reNtPNHviUrVQ9gX+Ykgz5VncyRKXTpNBIqNRzcy+svRJJ1JIiokiueve54hbYfYHZZfCQrUqiHlvyyp9BWRgSKyU0T2iMikYtbfIyIpIrLR8zPca91QEdnt+RlqRTzq4riMi3lr53HbR7eRdCaJDnU78OXdX2oSuAyBEghoIlD+yecSgYgEAq8B1wJJwDoRWVDMlJOfGGNGFNk3GngKiAUMsN6zr47nW8pSM1MZ/9V4vt/3PQD3dbuP8X3GExIYYnNk/il/qG1NBMofWVE1FAfsMcbsBRCRj4GbgYuZe/h6YJkx5qRn32XAQOAjC+JSJVh9cDVjF48lOT2Z6mHVmXHDDK5pcY3dYfk1bSxW/syKqqGGwCGv50meZUX9UUQ2i8jnItL4EvdFRO4XkUQRSUxJSbEg7MrH6XIyZ+Uc7vrsLpLTk+neqDsLhy7UJGCBwAB31ZA2Fit/VFY3hn8JNDXGdASWAf+61BcwxswzxsQaY2Jr19ZOTZfq2Nlj3PXpXcxNmIsxhhHxI3j/1vepX6W+3aFVCMEBnqoh7Vms/JAVVUOHgcZezxt5lhUwxqR6PX0beNFr36uL7LvCgpiUl+V7lzPhqwmczDpJ7cjazBw8k14xvewOq0LJLxE4jVYNKf9jRYlgHdBKRJqJSAhwO7DAewMR8b7svAnY7nm8BLhORGqISA3gOs8yZYFcZy7TVkxj+P8bzsmsk1zV9CoW3r1Qk0ApKOhH4NSqIeV/fC4RGGMcIjIC9wk8EJhvjPlFRKYCicaYBcAoEbkJcAAngXs8+54UkWdxJxOAqfkNx8o3h08fZuSXI9l0bBOBEsi4q8bx1+5/1WEiSkl+1ZCWCJQ/sqRDmTFmMbC4yLInvR5PBiaXsO98YL4VcSi37379jke/epTT2adpUKUBc26cQ9cGXe0Oq0LLrxrSNgLlj7RncQXicDmY+dNM3lz7JgDXNL+GGTfMoHp4dZsjq/jySwR615DyR5oIKohjZ48xZtEY1iWtI1ACefSqRxnefbhWBZWRgsZi7Ueg/JAmggpg5YGVjFk4hpNZJ6kbVZc5Q+bQvVF3u8OqVHTQOeXPNBH4MafLyWurX2PuqrkYDL2b9GbmoJnUiqxld2iVTv4QE1oiUP5IE4GfOpFxgnGLx/HTgZ8QhNG9RvNw/MMFVRSqbOmgc8qfaSLwQ4lJiYxaOIrk9GSiw6OZNXgWVza90u6wKrX8Yai1akj5I00EfsQYw1vr3uKlH1/CaZzENoxlzpA51KtSz+7QKj0ddE75M00EfuJU1ikmfD2Bb3/9FoAH4h5g7JVjC05Ayl4FcxZrPwLlh/Qs4gc2H93MiC9HcPjMYaqFVWPGDTPo36K/3WEpLwWJwGgLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 432x288 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"fig = plt.figure().gca()\n",
|
||||
"\n",
|
||||
"# target constraint at t=0.5\n",
|
||||
"fig.plot(pltx, target.velocity.data.flatten(), lw=2, color='forestgreen') \n",
|
||||
"\n",
|
||||
"# constrained state of simulation\n",
|
||||
"contrained_state = sess.run(states[16], feed_dict={state_in: stateN}).velocity.data\n",
|
||||
"fig.plot(pltx, contrained_state.flatten(), lw=2, color='mediumblue')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This seems to be going in the right direction! It's not there yet, but we've also only computed 5 steps so far. The two peaks with a positive velocity on the left side of the shock and the negative peak on the right side are starting to show.\n",
|
||||
"\n",
|
||||
"This is a good indicator that the backpropagation of gradients through all of our 16 simulated steps is behaving correctly, and driving the solution in the right direction. This hints at how powerful this setup is: the gradient that we obtain from each of the simulation steps (and each operation within them) can easily be chained together into more complex sequences. Here, we're backpropagation through all 16 steps of the simulation, and we could easily enlarge this \"look-ahead\" of the optimization without any additional effort.\n",
|
||||
"\n",
|
||||
"A simple direct reconstruction problem like this one is always a good initial test for a DP solver, e.g., before moving to more complex setups like coupling it with an NN. If the direct optimization does not converge, there's probably still something fundamentally wrong, and there's no point involving an NN. \n",
|
||||
"\n",
|
||||
"Our Burgers example here looks good - before including an NN, let's address the performance issue so: here we've updated our solutoin via a gradient that is \"manually\" computed with multiple TF session run() calls. TF's own optimizers do this internally, more flexibly and more efficiently. Let's start over, and switch to a `tf.train.Optimizer()` instead."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.8.5"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 4
|
||||
}
|
||||
+463
-123
@@ -4,13 +4,50 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Burgers Optimization with a (Manual) Differentiable Physics Gradient\n",
|
||||
"# Burgers Optimization with a Differentiable Physics Gradient\n",
|
||||
"\n",
|
||||
"To illustrate the process of compute gradients in a differentiable physics setting, we first demonstrate a process that isn't recommended for more complex real-world cases, but more clearly shows the way we can get gradient information from a physical simulation in a DL framework. We'll target tensorflow in the following.\n",
|
||||
"To illustrate the process of computing gradients in a _differentiable physics_ setting, we target the same reconstruction task like for the PINN example. This has some immediate implications: the evolution of the system is now fully determined by our PDE formulation. Hence, the only real unknown is the initial state! We will still need to re-compute all the states betwwen the initial and target state many times, just now we won't need an NN for this step. Instead, we can rely on our discretized model. \n",
|
||||
"\n",
|
||||
"We target the same reconstruction task like for the PINN example. This has some immediate implications: the evolution of the system is now fully determined by our PDE formulatio. Hence, the only real unknown is the initial state! We will still need to re-compute all the states betwwen the initial and target state many times, just now we won't need an NN for this step, we can rely on our discretized model. Also, as we choose an initial discretization for the DP approach, the unknown initial state consists of the sampling points of the involved physical fields, and we can simply represent these unknowns as floating point variables. Hence, even for the initial state we do not need to set up an NN. Thus, our Burgers reconstruction problem reduces to a gradient-based opitmization without any NN when solving it with DP. Nonetheless, it's a very good starting point to illustrate the process.\n",
|
||||
"Also, as we choose an initial discretization for the DP approach, the unknown initial state consists of the sampling points of the involved physical fields, and we can simply represent these unknowns as floating point variables. Hence, even for the initial state we do not need to set up an NN. Thus, our Burgers reconstruction problem reduces to a gradient-based opitmization without any NN when solving it with DP. Nonetheless, it's a very good starting point to illustrate the process.\n",
|
||||
"\n",
|
||||
"First, we'll set up our discretized simulation. Here we can employ the phiflow solver, as shown in the overview section on _Burgers forward simulations_. phiflow directly gives us a computational graph in tensorflow. So, as a first step, let's set up our grid with 128 points, and construct a tensorflow graph for 32 steps of the simulation. (As usual, you can ignore the warnings when loading TF and phiflow.)"
|
||||
"First, we'll set up our discretized simulation. Here we can employ phiflow, as shown in the overview section on _Burgers forward simulations_. \n",
|
||||
"\n",
|
||||
"## Initialization\n",
|
||||
"\n",
|
||||
"phiflow directly gives us a sequence of differentiable operations, provided that we don't use the _numpy_ backend.\n",
|
||||
"The important step here is to include `phi.tf.flow` instad of `phi.flow` (for _pytorch_ you could use `phi.torch.flow`).\n",
|
||||
"\n",
|
||||
"So, as a first step, let's set up some constants, and intialize our domain with 128 points. We'll directly allocate a `velocity` field, and our constraint at $t=0.5$ (step 16), now as a `scalar_grid` in phiflow.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"#!pip install --upgrade --quiet git+https://github.com/tum-pbs/PhiFlow@develop\n",
|
||||
"#!pip install --upgrade --quiet phiflow\n",
|
||||
"from phi.tf.flow import *\n",
|
||||
"\n",
|
||||
"N = 128\n",
|
||||
"STEPS = 32\n",
|
||||
"DT = 1./STEPS\n",
|
||||
"NU = 0.01/np.pi\n",
|
||||
"\n",
|
||||
"DOMAIN = Domain(x=N, boundaries=PERIODIC, bounds=Box[-1:1])\n",
|
||||
"\n",
|
||||
"# initialization of velocities\n",
|
||||
"velocity = DOMAIN.scalar_grid() # start from zero\n",
|
||||
"\n",
|
||||
"SOLUTION_T16 = DOMAIN.scalar_grid(np.asarray( [0.008612174447657694, 0.02584669669548606, 0.043136357266407785, 0.060491074685516746, 0.07793926183951633, 0.0954779141740818, 0.11311894389663882, 0.1308497114054023, 0.14867023658641343, 0.1665634396808965, 0.18452263429574314, 0.20253084411376132, 0.22057828799835133, 0.23865132431365316, 0.25673879161339097, 0.27483167307082423, 0.2929182325574904, 0.3109944766354339, 0.3290477753208284, 0.34707880794585116, 0.36507311960102307, 0.38303584302507954, 0.40094962955534186, 0.4188235294008765, 0.4366357052408043, 0.45439856841363885, 0.4720845505219581, 0.4897081943759776, 0.5072391070000235, 0.5247011051514834, 0.542067187709797, 0.5593576751669057, 0.5765465453632126, 0.5936507311857876, 0.6106452944663003, 0.6275435911624945, 0.6443221318186165, 0.6609900633731869, 0.67752574922899, 0.6939334022562877, 0.7101938106059631, 0.7263049537163667, 0.7422506131457406, 0.7580207366534812, 0.7736033721649875, 0.7889776974379873, 0.8041371279965555, 0.8190465276590387, 0.8337064887158392, 0.8480617965162781, 0.8621229412131242, 0.8758057344502199, 0.8891341984763013, 0.9019806505391214, 0.9143881632159129, 0.9261597966464793, 0.9373647624856912, 0.9476871303793314, 0.9572273019669029, 0.9654367940878237, 0.9724097482283165, 0.9767381835635638, 0.9669484658390122, 0.659083299684951, -0.659083180712816, -0.9669485121167052, -0.9767382069792288, -0.9724097635533602, -0.9654367970450167, -0.9572273263645859, -0.9476871280825523, -0.9373647681120841, -0.9261598056102645, -0.9143881718456056, -0.9019807055316369, -0.8891341634240081, -0.8758057205293912, -0.8621229450911845, -0.8480618138204272, -0.833706571569058, -0.8190466131476127, -0.8041372124868691, -0.7889777195422356, -0.7736033858767385, -0.758020740007683, -0.7422507481169578, -0.7263049162371344, -0.7101938950789042, -0.6939334061553678, -0.677525822052029, -0.6609901538934517, -0.6443222327338847, -0.6275436932970322, -0.6106454472814152, -0.5936507836778451, -0.5765466491708988, -0.5593578078967361, -0.5420672759411125, -0.5247011730988912, -0.5072391580614087, -0.4897082914472909, -0.47208460952428394, -0.4543985995006753, -0.4366355580500639, -0.41882350871539187, -0.40094955631843376, -0.38303594105786365, -0.36507302109186685, -0.3470786936847069, -0.3290476440540586, -0.31099441589505206, -0.2929180880304103, -0.27483158663081614, -0.2567388003912687, -0.2386513127155433, -0.22057831776499126, -0.20253089403524566, -0.18452269630486776, -0.1665634500729787, -0.14867027528284874, -0.13084990929476334, -0.1131191325854089, -0.09547794429803691, -0.07793928430794522, -0.06049114408297565, -0.0431364527809777, -0.025846763281087953, -0.00861212501518312] ))\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We can verify that the fields of our simulation are now backed by TensorFlow."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -19,41 +56,29 @@
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Each velocity is a phiflow grid like this one: Grid[128(1), size=[2.], <class 'numpy.float32'>]\n"
|
||||
]
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"tensorflow.python.framework.ops.EagerTensor"
|
||||
]
|
||||
},
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"from phi.tf.flow import *\n",
|
||||
"\n",
|
||||
"# run with phiflow\n",
|
||||
"n = 128\n",
|
||||
"dt = 1./32.\n",
|
||||
"initial = np.zeros([n,1]) # start from 0\n",
|
||||
"\n",
|
||||
"domain = Domain([n], boundaries=PERIODIC, box=box[-1:1])\n",
|
||||
"state0 = BurgersVelocity(domain, velocity=initial, viscosity=0.01/np.pi)\n",
|
||||
"\n",
|
||||
"# start with a state0, to be modified\n",
|
||||
"state_in = state0.copied_with(velocity=placeholder)\n",
|
||||
"states = [state_in]\n",
|
||||
"\n",
|
||||
"for i in range(32):\n",
|
||||
" states.append( Burgers().step(states[-1],dt=dt) )\n",
|
||||
"\n",
|
||||
"print(\"Note: each velocity is a phiflow grid like this one \" + format(states[-1].velocity) )"
|
||||
"type(velocity.values.native())"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"**A lot** has happened in these first steps, but we didnt run a single calculation for our simulation. In contrast to the first example, which used numpy as phiflow backend, we now have used the tensorflow backend (the trick here was importing `phi.tf.flow` instad of `phi.flow`). This gave us a tensorflow graph, but didn't run any real simulation steps. That will only happen once we run a tensorflow session and pass some data through this graph. So, in a way _nothing_ has happened so far in terms of our simulation! \n",
|
||||
"## Gradients\n",
|
||||
"\n",
|
||||
"Next, let's set up the loss, such that we can start the optimization. That's actually pretty simple, we want the solution at $t=0.5$, i.e. step number 16, to respect the data in terms of an L2 norm. Let's also directly initialize the TF variables, and see what loss we get from the initial zero'ed initialization:"
|
||||
"Now we can use the `record_gradients` function of phiflow to trigger the generation of a gradient tape to compute gradients of a simulation.\n",
|
||||
"\n",
|
||||
"For that we need a loss function: we want the solution at $t=0.5$ to respect the data in terms of an L2 norm. Thus we can simply compute an $L^2$ difference between step number 16 and our constraint array as `loss`. Afterwards, we evaluate the gradient of our initial velocity state `velocity` with respect to this loss."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -65,110 +90,63 @@
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/session.py:17: The name tf.Session is deprecated. Please use tf.compat.v1.Session instead.\n",
|
||||
"\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/session.py:18: The name tf.get_default_graph is deprecated. Please use tf.compat.v1.get_default_graph instead.\n",
|
||||
"\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/session.py:28: The name tf.global_variables_initializer is deprecated. Please use tf.compat.v1.global_variables_initializer instead.\n",
|
||||
"\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/Dropbox/mbaDevelSelected/phiflow-v1.5/phi/tf/session.py:29: The name tf.train.Saver is deprecated. Please use tf.compat.v1.train.Saver instead.\n",
|
||||
"\n",
|
||||
"Initial loss 0.3829152584075928\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",
|
||||
"Loss: 0.382915\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"sess = Session(None) \n",
|
||||
"velocities = [velocity]\n",
|
||||
"\n",
|
||||
"# enforce constraints at time t=0.5 (step 16)\n",
|
||||
"solution_t16 = np.asarray( [0.008612174447657694, 0.02584669669548606, 0.043136357266407785, 0.060491074685516746, 0.07793926183951633, 0.0954779141740818, 0.11311894389663882, 0.1308497114054023, 0.14867023658641343, 0.1665634396808965, 0.18452263429574314, 0.20253084411376132, 0.22057828799835133, 0.23865132431365316, 0.25673879161339097, 0.27483167307082423, 0.2929182325574904, 0.3109944766354339, 0.3290477753208284, 0.34707880794585116, 0.36507311960102307, 0.38303584302507954, 0.40094962955534186, 0.4188235294008765, 0.4366357052408043, 0.45439856841363885, 0.4720845505219581, 0.4897081943759776, 0.5072391070000235, 0.5247011051514834, 0.542067187709797, 0.5593576751669057, 0.5765465453632126, 0.5936507311857876, 0.6106452944663003, 0.6275435911624945, 0.6443221318186165, 0.6609900633731869, 0.67752574922899, 0.6939334022562877, 0.7101938106059631, 0.7263049537163667, 0.7422506131457406, 0.7580207366534812, 0.7736033721649875, 0.7889776974379873, 0.8041371279965555, 0.8190465276590387, 0.8337064887158392, 0.8480617965162781, 0.8621229412131242, 0.8758057344502199, 0.8891341984763013, 0.9019806505391214, 0.9143881632159129, 0.9261597966464793, 0.9373647624856912, 0.9476871303793314, 0.9572273019669029, 0.9654367940878237, 0.9724097482283165, 0.9767381835635638, 0.9669484658390122, 0.659083299684951, -0.659083180712816, -0.9669485121167052, -0.9767382069792288, -0.9724097635533602, -0.9654367970450167, -0.9572273263645859, -0.9476871280825523, -0.9373647681120841, -0.9261598056102645, -0.9143881718456056, -0.9019807055316369, -0.8891341634240081, -0.8758057205293912, -0.8621229450911845, -0.8480618138204272, -0.833706571569058, -0.8190466131476127, -0.8041372124868691, -0.7889777195422356, -0.7736033858767385, -0.758020740007683, -0.7422507481169578, -0.7263049162371344, -0.7101938950789042, -0.6939334061553678, -0.677525822052029, -0.6609901538934517, -0.6443222327338847, -0.6275436932970322, -0.6106454472814152, -0.5936507836778451, -0.5765466491708988, -0.5593578078967361, -0.5420672759411125, -0.5247011730988912, -0.5072391580614087, -0.4897082914472909, -0.47208460952428394, -0.4543985995006753, -0.4366355580500639, -0.41882350871539187, -0.40094955631843376, -0.38303594105786365, -0.36507302109186685, -0.3470786936847069, -0.3290476440540586, -0.31099441589505206, -0.2929180880304103, -0.27483158663081614, -0.2567388003912687, -0.2386513127155433, -0.22057831776499126, -0.20253089403524566, -0.18452269630486776, -0.1665634500729787, -0.14867027528284874, -0.13084990929476334, -0.1131191325854089, -0.09547794429803691, -0.07793928430794522, -0.06049114408297565, -0.0431364527809777, -0.025846763281087953, -0.00861212501518312] )\n",
|
||||
"target = state0.copied_with(velocity=np.reshape(solution_t16,[n,1]))\n",
|
||||
"with math.record_gradients(velocity.values):\n",
|
||||
"\n",
|
||||
"loss = math.sum( (states[16].velocity.data - target.velocity.data)**2 ) / n\n",
|
||||
"sess.initialize_variables()\n",
|
||||
" for time_step in range(STEPS):\n",
|
||||
" v1 = diffuse.explicit(1.0*velocities[-1], NU, DT)\n",
|
||||
" v2 = advect.semi_lagrangian(v1, v1, DT)\n",
|
||||
" velocities.append(v2)\n",
|
||||
"\n",
|
||||
"# compute initial loss\n",
|
||||
"s,l = sess.run([states[-1],loss], feed_dict={state_in: state0})\n",
|
||||
"print(\"Initial loss \"+format(l))\n"
|
||||
" loss = field.l2_loss(velocities[16] - SOLUTION_T16)*2./N # MSE\n",
|
||||
"\n",
|
||||
" grad = math.gradients(loss, velocity.values)\n",
|
||||
"\n",
|
||||
"print('Loss: %f' % (loss))\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Not surprisingly, because we're starting from zero, there's a significant intial error of ca. 0.38 for the 16th simulation step.\n",
|
||||
"Because we're only constraining timestep 16, we could actually omit steps 17 to 31 in this setup. They don't have any degrees of freedom and are not constrained in any way. For fairness and for direct comparison with the previous PINN case, we include them.\n",
|
||||
"\n",
|
||||
"Note that because we're only constraining timestep 16, we could actually omit steps 17 to 31 in this setup. They don't have any degrees of freedom and are not constrained in any way. For fairness and for direct comparison with the previous PINN case, let's include them.\n",
|
||||
"Note that we've done a lot of calcuations here: first the 32 steps of our simulation, and then another 16 steps backwards from the loss. There were recorded by the gradient tape, and used to backpropagate the loss to the initial state of the simulation.\n",
|
||||
"\n",
|
||||
"Now we have our simulation graph in TF, we can use TF to give us a gradient for the initial state for the loss. All we need to do is run `tf.gradients(loss, [state_in.velocity.data]`, which will give us a direction for each velocity variable. Based on a linear approximation, this direction tells us how to change each of them to increase the loss function (gradients _always_ point \"upwards\"). In the following code snipped, we're additionally saving all these gradients in a list called `grads`, such that we can visualize them later on. (Normally, we could discard each gradient after performing an update step.)\n",
|
||||
"Not surprisingly, because we're starting from zero, there's also a significant intial error of ca. 0.38 for the 16th simulation step.\n",
|
||||
"\n",
|
||||
"Based on the gradient, we can now take a step in the opposite direction to bring the loss down (instead of increasing it). Below we're using a learning rate `LR=5` for this step. Afterwards, we're re-evaluating the loss for the updated state to check how we did. As these updates aren't exactly fast, we're only computing five of them here:"
|
||||
"So what do we get as a gradient here? It has the same dimensions as the velocity, and we can easily visualize it:\n",
|
||||
"Starting from the zero state for `velocity` (shown in blue), the first gradient is shown as a green line below. If you compare it with the solution it points in the opposite direction, as expected. The solution is much larger in magnitude, so we omit it here (see the next graph).\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"WARNING:tensorflow:From /usr/local/lib/python3.7/site-packages/tensorflow_core/python/ops/math_grad.py:1375: where (from tensorflow.python.ops.array_ops) is deprecated and will be removed in a future version.\n",
|
||||
"Instructions for updating:\n",
|
||||
"Use tf.where in 2.0, which has the same broadcast rule as np.where\n",
|
||||
"Step 0, loss = 0.3829152584075928\n",
|
||||
"Step 1, loss = 0.3334265947341919\n",
|
||||
"Step 2, loss = 0.2905231714248657\n",
|
||||
"Step 3, loss = 0.25346457958221436\n",
|
||||
"Step 4, loss = 0.2215319573879242\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"LR = 5.\n",
|
||||
"stateN = state0\n",
|
||||
"grads = []\n",
|
||||
"for i in range(5):\n",
|
||||
" grads.append( tf.gradients(loss, [state_in.velocity.data] ) ) # use TF\n",
|
||||
" grads.append( gradients(loss, state_in.velocity) ) # phiflow wrapper also works\n",
|
||||
"\n",
|
||||
" state_updated = state_in.copied_with( velocity=(state_in.velocity - LR*grads[-1]) )\n",
|
||||
" stateN,l = sess.run([state_updated,loss], feed_dict={state_in: stateN})\n",
|
||||
"\n",
|
||||
" #s,l = sess.run([states[-1],loss], feed_dict={state_in: stateN})\n",
|
||||
" print( \"Step {}, loss = {}\".format(i,l) )\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"One of the reasons for the bad performance here is that we're querying states via `sess.run()` multiple times, instead of computing the update in one go. Nonetheless, this simple example already clearly shows that the loss is nicely going down: from above 0.38 to ca. 0.22. This is a good sign!\n",
|
||||
"\n",
|
||||
"Before improving on the performance, let's visualize the reconstruction.\n",
|
||||
"\n",
|
||||
"Starting from the zero state `state0` (shown in blue), the first gradient is shown as a red line here. If you compare it with the solution it points in the opposite direction, as expected. The solution is much larger in magnitude, so we omit it here (see the next graph).\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[<matplotlib.lines.Line2D at 0x13fa0ae90>]"
|
||||
"[<matplotlib.lines.Line2D at 0x7fead6acbb80>]"
|
||||
]
|
||||
},
|
||||
"execution_count": 5,
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAD4CAYAAAAUymoqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de9zUc/r48dfV3UkpKreWQiFs2E3uDYtO0kkqhApltVrUOsR3ZX2t0+5+sb+VczY5xC6VsN1YkpKz6o6UIt2FVZu6VXIu5fr9cb1nG7e5T83c85nD9Xw85jGf81wz3c01n/dRVBXnnHOuJupEHYBzzrns48nDOedcjXnycM45V2OePJxzztWYJw/nnHM1VjfqANJht9120zZt2kQdhnPOZZUFCxZ8qqqFifblRfJo06YNJSUlUYfhnHNZRUQ+qmifF1s555yrsZQkDxHpLSLLRKRURMYm2N9ARKaE/XNFpE3Y3kJEXhCRL0XkjnLnHC4ii8M5t4mIhO3NRWSmiCwPz81S8R6cc85VX9LJQ0QKgDuBPkB7YIiItC932Ahgo6ruD4wDbgzbvwWuAi5LcOnxwLlAu/DoHbaPBWapajtgVlh3zjmXRqm48+gElKrqSlXdAkwGBpQ7ZgAwKSxPA44TEVHVr1T1FSyJ/JeI7AE0VdU31MZPeRAYmOBak+K2O+ecS5NUJI9WwMdx66vCtoTHqOpWYBPQooprrqrgmi1VdU1Y/gRomegCIjJSREpEpKSsrKw678M551w1ZXWFebgrSTiyo6pOUNUiVS0qLEzY0sw559wOSkXyWA3sFbfeOmxLeIyI1AV2AdZXcc3WFVxzbSjWihVvrdvhyJ1zzu2QVPTzmA+0E5G22Bf8YGBouWOKgeHA68AgYLZWMha8qq4Rkc9F5EhgLjAMuL3ctW4Iz9NT8B6cc+Vt2QIvvQQvvwxz58L330OjRtC6NXTrBl27QovKSp9dLks6eajqVhEZDcwACoD7VHWJiFwHlKhqMXAv8JCIlAIbsAQDgIh8CDQF6ovIQKCnqi4FLgAeAHYCngkPsKQxVURGAB8BpyX7HpxzcVThqafgkktgxYrEx9x5J9SpA/36wejRcNxxtu7yhuTDZFBFRUXqPcydq4YvvoAzz4TiYlvff3848UQ4+mjYeWf46itYuhRmz4ZXXoHvvrPjDj0UbrgB+vQB65LlcoCILFDVooT7PHk45wD4/HPo2xdefRWaNoVrr4VRo6BevcTHr10LEyfC+PGwOlRJdukCt99uycRlvcqSh99nOucscfTubYljr73gzTfh4osrThwALVvClVdCaSn89a/QvDm8+CJ07AiXX253KS5nefJwLt+pwogR8PrrsPfeMGcO7Ldf9c9v2BDGjLH6kVGjYNs2uOkmOOwwmD+/1sJ20fLk4Vy+e+ghmDYNmjSBWbNg33137Dq77gp33GEtsw45BJYvh6OOgj/+EbZuTW3MLnKePJzLZx9+aK2lwOoq9t8/+Wv+4hd2x3HJJXYXctVVVheycmXy13YZw5OHc/lKFc4+21pYnXIKDBuWums3bAg33wwzZ0KrVvDaa/Dzn8Mjj6TuNVykPHk4l6+Ki62Cu7AQ/va32mli26MHLFoEp54KX34JQ4fChRdaB0SX1Tx5OJePtm6FK66w5T/8oXZ7ijdvDlOmwF13Weut22+33umrVlV5qstcnjycy0eTJsG771rl+MiRtf96InD++TbUSevW1rKrY0frbOiykicP5/LNN9/A1Vfb8h//CPXrp++1jzjC+pD06AFlZXD88XDbbVb/4rKKJw/n8s1991mP8A4d4PTT0//6hYXw7LPw+9/bYIsXXQTnnbd9qBOXFTx5OJdPtm2DW26x5SuvjG4ww4IC+NOf4OGHoUEDmDABevaE9ZXN1OAyiScP5/LJU0/ZcCJt2sDADJjBecgQG/b9Jz+xnu1HHGF1MS7jefJwLp+MG2fPF10EdVMxnU8KdOpknQo7drQhTo480oq1XEbz5OFcvliwwPp1NGkC55wTdTQ/1Lq13YEMGmSDNJ5wAtx6q1ekZzBPHs7li1tvtedzz7Uh1zNN48bWH+Sqq6wi/eKL4be/tXoal3E8eTiXDzZuhEcfteVRo6KNpTJ16sB111lFev36NmPhKafA119HHZkrJyXJQ0R6i8gyESkVkbEJ9jcQkSlh/1wRaRO374qwfZmI9ArbDhSRhXGPz0Xk4rDvGhFZHbevbyreg3M57eGH4dtvrX/Fjo6am05Dhti4WLvuCtOn2zS3n34adVQuTtLJQ0QKgDuBPkB7YIiItC932Ahgo6ruD4wDbgzntsfmMz8Y6A3cJSIFqrpMVTuoagfgcOBr4Im4642L7VfVfyX7HpzLeffea8+//nW0cdRE5842OdXee8Mbb8Avf1nxnOou7VJx59EJKFXVlaq6BZgMDCh3zABgUlieBhwnIhK2T1bVzar6AVAarhfvOGCFqn6Uglidyz9vvglvvWVjTGVC89yaaN/eEsdhh22fH2TevKijcqQmebQCPo5bXxW2JTxGVbcCm4AW1Tx3MFB+HOfRIrJIRO4TkWaJghKRkSJSIiIlZWVlNXk/zuWW2F3HWWdZh7xss8ce1kqsVy8b0qRrV3jyyaijynsZXWEuIvWB/sCjcZvHA/sBHYA1wF8TnauqE1S1SFWLCgsLaz1W5zLSN9/AP/5hyyNGRBtLMpo0sYTxq1/Zexo4EO65J+qo8loqksdqYK+49dZhW8JjRKQusAuwvhrn9gHeVNW1sQ2qulZVt6nq98A9/LiYyzkX8+STsGkTFBXBoYdGHU1y6tWzu6g//MGa8o4cCf/3f94XJCKpSB7zgXYi0jbcKQwGissdUwwMD8uDgNmqqmH74NAaqy3QDogv0BxCuSIrEdkjbvUk4J0UvAfnclPsruPMM6ONI1VE4NprbW4QERtc8dJLLZm4tEp6fAJV3Soio4EZQAFwn6ouEZHrgBJVLQbuBR4SkVJgA5ZgCMdNBZYCW4FRqroNQEQaA8cDvyn3kjeJSAdAgQ8T7HfOgQ0y+Mwz1nciitFza9P559sEVmeeaUOufPqp3ZXUqxd1ZHlDNA9u+YqKirSkpCTqMJxLr7vvti/ZXr1yd6yomTPhpJPgq6+gXz/rod6oUdRR5QwRWaCqRYn2ZXSFuXMuCbEiqzPOiDaO2nT88TBrljVDfuopS5SffRZ1VHnBk4dzuejDD+GVV2CnnbKvb0dNHXGEvdfWre25SxdYsybqqHKeJw/nctEjoZ3JgAHWzDXX/fSn1hv9wANh0SI45hhLoK7WePJwLhdNmWLPQ4ZEG0c67b03vPwyHH44rFxpw5u8/37UUeUsTx7O5Zply+Dtt2GXXawOIJ8UFlodyNFHw8cfWwJZvDjqqHKSJw/nck3srmPgwOwcjiRZu+wCM2bYCMJr19pwJt7aMuU8eTiXa6ZOtedc69tRE40bW+/6E0+EDRuge3erTHcp48nDuVyyZIk9mjWzOTDyWcOG8NhjlkS/+AJ69rR+IS4lPHk4l0tidx0nn2wz8eW7evWsv0tsQMV+/aC4/OhJbkd48nAuV6hur+/I5yKr8goKYOJEmw99yxZLrJMnRx1V1vPk4VyuWLTIWlrttht06xZ1NJmlTh249VYYOxa2bYOhQ+G++6KOKqt58nAuV8SKrE45BeomPeZp7hGxIdz/9Ce7SxsxAu68M+qospYnD+dygRdZVd/vf28j8QKMHm13JK7GPHk4lwvefBNWrICWLa1jnKvcxRfDHXdsX7755mjjyUKePJzLBbEiq0GDrILYVW3UKPjb32z50kvhppuijSfLePJwLtupesfAHTVypE0iJQKXXw5//nPUEWUNTx7OZbt582wE2T33tDGdXM2ccw7cf78lkCuvhOuuizqirJCS5CEivUVkmYiUisjYBPsbiMiUsH+uiLSJ23dF2L5MRHrFbf9QRBaLyEIRKYnb3lxEZorI8vDcLBXvwbmsFRt+/bTTrEmqq7nhw+HBB+3zu/pq+MMf7I7OVSjpvzQRKQDuBPoA7YEhItK+3GEjgI2quj8wDrgxnNsem8/8YKA3cFe4Xkw3Ve1QbhrEscAsVW0HzArrzuWnbdu2d3gbOjTaWLLdmWdab/SCArj+ersL8QRSoVT8TOkElKrqSlXdAkwGBpQ7ZgAwKSxPA44TEQnbJ6vqZlX9ACgN16tM/LUmATk+TZpzlXjhBRs5dv/9oSjhVNOuJgYPtmRct671Cbn8ck8gFUhF8mgFfBy3vipsS3iMqm4FNgEtqjhXgedEZIGIjIw7pqWqxuaY/ARomSgoERkpIiUiUlJWVlbzd+VcNnj4YXs+4wwrs3fJGzTIGiDUrQt/+Yu1xPIE8iOZXEB6jKp2xIrDRonIjxqvq6piSeZHVHWCqhapalFhYWEth+pcBL791kaNhfyaMTAdTjrJPtt69axD4WWXeQIpJxXJYzWwV9x667At4TEiUhfYBVhf2bmqGnteBzzB9uKstSKyR7jWHsC6FLwH57LP00/D55/btKsHHhh1NLmnf394/HFLIDffDL/7nSeQOKlIHvOBdiLSVkTqYxXg5cc8LgaGh+VBwOxw11AMDA6tsdoC7YB5ItJYRJoAiEhjoCfwToJrDQemp+A9OJd9YkVWXlFee/r1g2nTLIH8v/9nAyt6AgFSkDxCHcZoYAbwLjBVVZeIyHUi0j8cdi/QQkRKgTGEFlKqugSYCiwFngVGqeo2rB7jFRF5G5gHPK2qz4Zr3QAcLyLLgR5h3bn88tlnduch4h0Da1v//tvrQG66ycbG8gSCaB58CEVFRVricxi7XHL//da5rVs3mD076mjyw+OPW1+abdusGe/11+d8IwURWVCuq8R/ZXKFuXOuIl5klX6xSaQKCmxY92uuiTqiSHnycC7brFljdxv16tncHS59Bg2yHv0FBTaMybXXRh1RZDx5OJdLine truncated
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD4CAYAAAAZ1BptAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwxUlEQVR4nO3dd3wUdf7H8dcnnd4FpAuIdBICiCgo0kWKgKJUBRJE5U7Penqnh6enng1EJSSgiAUR4QhIEQEVlZaQUASBiKKAQKSXJJDk+/sjk/xW3CWBLZPdfJ6Pxz6yM/Od2TeTJZ+d73d2RowxKKWUUs4E2R1AKaVU8aVFQimllEtaJJRSSrmkRUIppZRLWiSUUkq5FGJ3AE+qWrWqqV+/vt0xlFLKryQnJ/9ujKnmbFlAFYn69euTlJRkdwyllPIrIrLX1TLtblJKKeWSFgmllFIuaZFQSinlkhYJpZRSLmmRUEop5ZIWCaWUUi5pkVBKKeWSR74nISK9gMlAMJBgjHnhguXhwHtAW+AIcIcx5mcRqQLMA9oB7xpj7ndY50ugJpBhzephjDnsibxKFQdJB5J4atVTlA4tzRVlrqDlFS3p3bg3V1W6yu5oShVwu0iISDDwJtAd2AdsFJFEY8x2h2ZjgGPGmEYiMhR4EbgDyAT+AbSwHhcaZozRb8epgHMu5xwjFozgh99/+OOCpdCkShNub347d7a4k6bVmtoTUCmLJ7qb2gNpxpg9xphzwByg/wVt+gOzrOfzgJtFRIwxZ4wx35BXLJQqMV5b+xo//P4DV1e5mrmD5zKl1xSGNBtChfAK7Dyyk2e/fpZmbzWj44yOfLDlA7Kys+yOrEooT3Q31QJ+dZjeB3Rw1cYYky0iJ4AqwO+FbPsdEckBPgX+bZzcRk9EYoAYgLp1617WP0ApX/rlxC9M+noSAFN7T6V7w+4APNDhAbJzs/ny5y/5cOuHfLrjU9btW8e6fet46POHGBc1jti2sdSpUMfO+KqEKc4D18OMMS2BG6zHCGeNjDHTjTHRxpjoatWcXp9KqWLlkRWPcPb8WYY0G1JQIPKFBIXQ7apuzOw/kwMPHSCubxwtr2jJ4TOHeW7NczSY3ICRC0ayI32HTelVSeOJIrEfcPxoU9ua57SNiIQAFcgbwHbJGLPf+nkK+JC8bi2l/N5nuz4D4KXuL120XZmwMsS0jWHz+M18Pfpr7mh+BwCzt8ym+VvNuf2T29lyaIvX86qSzRNFYiPQWEQaiEgYMBRIvKBNIjDKej4YWOWs6yifiISISFXreSjQF9jmgaxK2epczjnOnD9DsARTr0K9Iq0jItxQ7wbmDJ7D7gd2E9s2lpCgED7Z/gmtp7VmwJwBbD201cvJVUnldpEwxmQD9wPLgR3AXGPM9yIySUT6Wc1mAFVEJA14CHg8f30R+Rl4FRgtIvtEpBkQDiwXkS1AKnlHIvHuZlXKbsczjwNQMaIiInLJ6zeo1IBpfaex5y97mNh+IhEhESzcuZDW01ozcsFIfj7+s2cDqxJPLvKB3u9ER0cbvZ+EKs52HdlFk6lNaFipIWkT09ze3qHTh3huzXNMS5rG+dzzhAaFcm/0vTzV+SmqldExOlU0IpJsjIl2tqw4D1wrFXDyjyQqlarkke1VL1udKb2n8MP9PzCs5TCyc7OZsmEKjd5oxCvfvcK5nHMeeR1VcmmRUMqHjmUcA/K6mzzpqkpX8f5t75MSm0LvRr05mXWSh1c8TIu3WhQMlCt1ObRIKOVDBUcSEZ45krhQ6xqtWTJsCUvuWkKTKk3YfXQ3fT/qS58P+vz5291KFYEWCaV86Fimd44kLtS7cW+23LuFV3u8Svnw8ixNW0rLt1vyt+V/40TmCa++tgosWiSU8iFvH0k4CgsO48GOD7L7gd2MixpHTm4Or657lcZvNGb25tkE0kkrynu0SCjlQ94ak7iYK8pcwfRbp5Mck8z1da8n/Ww6I/83km6zu7Hz950+y6H8kxYJpXzI02c3XYrImpF8Pfpr3u3/LlVKVWHVT6toNa0VT69+msxsvcamck6LhFI+5KsxCVdEhFFtRrHz/p3c0+YezuWcY9LXk2j1diu+2POFLZlU8aZFQikf8uWYxMVUKV2FGf1n8NXor2hatSm7j+6m++zuDJs/jMNn9N5e6v9pkVDKh+w+krhQ53qdSR2fyvNdnyciJIIPt35Iszeb8dHWj3RgWwFaJJTyKTvHJFwJCw7jiRue4PsJ33Nzg5s5knGEu+bfxYCPB3Dg1AG74ymbaZFQyofsOLupqK6qdBUrRqwg/tZ4yoeXJ3FnIs3ebMbMlJl6VFGCaZFQykeMMX+4CmxxJCKMjRrL9xO+p+/VfTmRdYIxiWPo+X5PvcJsCaVFQikfOX3uNDkmh9KhpQkLDrM7zkXVLl+bxKGJvD/wfSqXqsyKPSto8VYL3tzwJrkm1+54yoe0SCjlI8XlzKaiEhGGtRrG9gnbGdJsCGfOn+H+pfdz47s3suvILrvjKR/RIqGUjxT3riZXqpetztwhc5k3ZB7Vy1RnzS9raDOtDW+sf0OPKkoALRJK+Uj+6a/F6cymSzGo2SC237ed4a2Gk5GdwcRlE+k+uzu/nPjF7mjKi7RIKOUj/nok4ahyqcrMHjibT2//lKqlq7Lqp1W0fLsls1Jn6RlQAUqLhFI+kn/6q7+MSVzMbU1vY9u92+jfpD8ns04yeuFoBn48UL+tHYC0SCjlI4FwJOGoetnqLLhjAe/2f5fy4eVZuHMhLd5qwYIdC+yOpjxIi4RSPlIwJhEARxL58i8YuPXerXRt0JX0s+ncNvc2Ri4YWVAUlX/TIqGUjwTakYSjuhXqsmLECqb0mkJESASzt8ym5dstWblnpd3RlJu0SCjlI/5+dlNhgiSIBzo8QGpsKu1rtWffyX10n92dR1c8SlZ2lt3x1GXySJEQkV4islNE0kTkcSfLw0XkY2v5ehGpb82vIiKrReS0iEy9YJ22IrLVWmeKiIgnsipll0A+knDUpGoTvr3nW57p8gxBEsR/v/svHWd0ZEf6DrujqcvgdpEQkWDgTaA30Ay4U0SaXdBsDHDMGNMIeA140ZqfCfwDeNjJpt8GxgGNrUcvd7MqZadAOrupMCFBITx949OsuXsNV1W6ipSDKbSd3pa3N76tp8r6GU8cSbQH0owxe4wx54A5QP8L2vQHZlnP5wE3i4gYY84YY74hr1gUEJGaQHljzDqT9456DxjggaxK2aakHEk46linI6mxqYxqPYqM7AwmLJlAvzn99FRZP+KJIlEL+NVhep81z2kbY0w2cAKoUsg29xWyTQBEJEZEkkQkKT09/RKjK+U7gT4m4Uq58HK8O+Bd5gyaQ8WIiizetZhWb7diWdoyu6OpIvD7gWtjzHRjTLQxJrpatWp2x1HKpZJ4JOHojhZ3sGX8FrrU68KhM4fo/UFv/rL0L2Scz7A7mroITxSJ/UAdh+na1jynbUQkBKgAHClkm7UL2aZSfuN8znlOnztNkARRLqyc3XFsU6dCHVaOXMkLN79ASFAIUzZMoX1Ce7Ye2mp3NOWCJ4rERqCxiDQQkTBgKJB4QZtEYJT1fDCwylxk9MoY8xtwUkSutc5qGgks9EBWpWxxIusEkHcUUdJP1AsOCuax6x9j3Zh1XF3larYd3ka7+HZMXjdZrypbDLldJKwxhvuB5cAOYK4x5nsRmSQi/axmM4AqIpIGPAQUnCYrIj8DrwKjRWSfw5lRE4AEIA34EVjqblal7FKSzmwqqrZXtmVTzCZiomLIysnir8v/yi0f3qKD2sWMBNLpaNHR0SYpKcnuGEr9ycb9G2mf0J62NduSFKPv0Qv974f/MSZxDEczjlKjbA3eH/g+N191s92xSgwRSTbGRDtb5vcD10r5g4K70pWwM5uKasA1A9g8fjOd63Xm4OmDdJ/dnSdXPkl2brbd0Uo8LRJK+UD+mET58PI2Jym+apevzaqRq3imyzOICM9/8zxd3u3C3uN77Y5WommRUMoHzuWcAyA8ONzmJMVbcFAwT9/4NKtGrqJWuVp89+t3tIlrw/wd8+2OVmJpkVDKB87nnAcgNDjU5iT+oUv9LqSOT6Xv1X05nnmcQXMHMeGzCfqdChtokVDKB/KPJEKDtEgUVdXSVUkcmsjkXpMJCw7j7aS36ZDQQS8U6GNaJJTygfO5eUcSYcFhNifxLyLCxA4TWTtmLY0rN2br4a1Ex0czM2WmXijQR7RIKOUDBd1NeiRxWaJqRpEck8yIViM4e/4sYxLHMGz+ME5mnbQ7WsDTIqGUD+QfSeiYxOUrF16O9wa+x6wBsygTWoaPtn1EZFwkG/dvtDtaQNMioZQP5I9JaHeT+0a2Hsmm2E20qdGGPcf2cN3M63jlu1e0+8lLtEgo5QPa3eRZV1e5mnVj1jGx/USyc7N5eMXD9J/Tn6MZR+2OFnC0SCjlA9rd5HnhIeFM7j2Z/93xPypGVGTRrkVExUWxft96u6MFFC0SSvmAdjd5T/9r+rMpZhPtrmzH3hN7ueGdG3h93eva/eQhWiSU8gHtbvKuBpUa8M093zCx/UTO557nweUPMmjuoIJrZqnLp0VCKR/Q7ibvCwsOY3LvycwbMo/y4eVZ8MMCouKiSDqgV911hxYJpXxAu5t8Z1CzQWyK2URUzSh+Ov4TnWZ24s0Nb2r302XSIqGUDxQcSWh3k080rNyQb+/5lgnREziXc477l97PHfPu4ETmCbuj+R0tEkr5gF7gz/ciQiJ485Y3mTNoDuXCyvHJ9k+Ijo8m9WCq3dH8ihYJpXxAL/Bnnzta3EFyTDKtq7cm7Wga1yZcS1xSnHY/FZEWCaV8QC/wZ6/GVRqzdszagvtpj/9sPMPmD+NU1im7oxV7WiSU8gHtbrJfqdBSxN0axwe3fVBw7afo+Gi2HNpid7RiTYuEUj6g3U3Fx10t7yIpJokWV7Rg15FddEjoQMKmBO1+ckGLhFI+oN1Line truncated
|
||||
"text/plain": [
|
||||
"<Figure size 432x288 with 1 Axes>"
|
||||
]
|
||||
@@ -183,42 +161,96 @@
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"\n",
|
||||
"fig = plt.figure().gca()\n",
|
||||
"pltx = np.linspace(-1,1,n)\n",
|
||||
"pltx = np.linspace(-1,1,N)\n",
|
||||
"\n",
|
||||
"# first gradient\n",
|
||||
"[grad0,s0] = sess.run([grads[0],state0], feed_dict={state_in: stateN}) #.velocity.data\n",
|
||||
"fig.plot(pltx, grad0[0].flatten() , lw=2, color='red') \n",
|
||||
"fig.plot(pltx, grad.numpy('x') , lw=2, color='green') \n",
|
||||
"\n",
|
||||
"fig.plot(pltx, state0.velocity.data.flatten(), lw=2, color='mediumblue')"
|
||||
"fig.plot(pltx, velocity.values.numpy('x'), lw=2, color='mediumblue')\n",
|
||||
"\n",
|
||||
"# some (optional) other fields to plot:\n",
|
||||
"#fig.plot(pltx, (velocities[16]).values.numpy('x') , lw=2, color='cyan') \n",
|
||||
"#fig.plot(pltx, (SOLUTION_T16).values.numpy('x') , lw=2, color='red') \n",
|
||||
"#fig.plot(pltx, (velocities[16] - SOLUTION_T16).values.numpy('x') , lw=2, color='blue') \n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"So the gradient looks good (feel free to plot the others from the five iterations above).\n",
|
||||
"Now we have our simulation graph in TF, we can use TF to give us a gradient for the initial state for the loss. All we need to do is run `tf.gradients(loss, [state_in.velocity.data]`, which will give us a \n",
|
||||
"\n",
|
||||
"How well does the 16th state of the simulation actually match the target after the 5 update steps? This is what the loss measures, after all. The next graph shows the constraints (i.e. the solution we'd like to obtain) in green, and the reconstructed state after the initial state was updated 16 times by the solver:"
|
||||
"Thus now we have \"search direction\" for each velocity variable. Based on a linear approximation, the gradient tells us how to change each of them to increase the loss function (gradients _always_ point \"upwards\"). In the following code block, we're additionally saving all these gradients in a list called `grads`, such that we can visualize them later on. (Normally, we could discard each gradient after performing an update step.)\n",
|
||||
"\n",
|
||||
"Based on the gradient, we can now take a step in the opposite direction to bring the loss down (instead of increasing it). Below we're using a learning rate `LR=5` for this step. Afterwards, we're re-evaluating the loss for the updated state to check how we did. "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Optimization step 0, loss: 0.382915\n",
|
||||
"Optimization step 1, loss: 0.326966\n",
|
||||
"Optimization step 2, loss: 0.281119\n",
|
||||
"Optimization step 3, loss: 0.242881\n",
|
||||
"Optimization step 4, loss: 0.210734\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"LR = 5.\n",
|
||||
"\n",
|
||||
"grads=[]\n",
|
||||
"for optim_step in range(5):\n",
|
||||
" velocities = [velocity]\n",
|
||||
" with math.record_gradients(velocity.values):\n",
|
||||
" for time_step in range(STEPS):\n",
|
||||
" v1 = diffuse.explicit(1.0*velocities[-1], NU, DT)\n",
|
||||
" v2 = advect.semi_lagrangian(v1, v1, DT)\n",
|
||||
" velocities.append(v2)\n",
|
||||
"\n",
|
||||
" loss = field.l2_loss(velocities[16] - SOLUTION_T16)*2./N # MSE\n",
|
||||
" print('Optimization step %d, loss: %f' % (optim_step,loss))\n",
|
||||
"\n",
|
||||
" grads.append( math.gradients(loss, velocity.values) )\n",
|
||||
"\n",
|
||||
" velocity = velocity - LR * grads[-1]\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"\n",
|
||||
"Now we can check well the 16th state of the simulation actually matches the target after the 5 update steps. This is what the loss measures, after all. The next graph shows the constraints (i.e. the solution we'd like to obtain) in green, and the reconstructed state after the initial state `velocity` (which we updated five times via the gradient by now) was updated 16 times by the solver. \n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[<matplotlib.lines.Line2D at 0x13fbd43d0>]"
|
||||
"[<matplotlib.lines.Line2D at 0x7fead6b6ae20>]"
|
||||
]
|
||||
},
|
||||
"execution_count": 8,
|
||||
"execution_count": 6,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3xUVfr48c+T3miht9CbSyeEAAoiqAio+911LbsqKqxlpQnSxIrKoihNXRWVddfefruLgAIqWCCUIFW61FBCCKGkZ2bO74+ZxElIaHOTm0me9+uVV2ZumXlyZ3Kfe8655xwxxqCUUqryCrA7AKWUUvbSRKCUUpWcJgKllKrkNBEopVQlp4lAKaUquSC7A7gctWrVMk2bNrU7DKWU8ivr168/YYypXXS5XyaCpk2bkpiYaHcYSinlV0TkQHHLtWpIKaUqOU0ESilVyWkiUEqpSs6SRCAi80XkuIhsLWG9iMhcEdkjIptFpKvXuqEistvzM9SKeJRSSl08q0oE7wIDz7P+BqCV5+d+4HUAEYkGngJ6AHHAUyJSw6KYlFJKXQRLEoEx5gfg5Hk2uRn4t3FbDVQXkfrA9cAyY8xJY0wasIzzJxSllFIWK6s2gobAIa/nSZ5lJS0/h4jcLyKJIpKYkpJSaoEqpVRl4zf9CIwx84B5ALGxsTp2tipTuc5c1h9ez5ZjW9iRsoMcRw43XXET/Vv0JyjAb/6NlCpWWX2DDwONvZ438iw7DFxdZPmKMopJqfPKyM3g21+/5etdX/PT/p/IyMsotP7r3V9TL6oeMwbNoFdML5uiVMp3ZZUIFgAjRORj3A3Dp40xR0VkCTDNq4H4OmByGcWk1DmMMaxLWsfHmz9m6e6lZDmyCta1rtWaHo160LZOWzJzM/lg0wfsT9vPJ5s/0USg/JoliUBEPsJ9ZV9LRJJw3wkUDGCMeQNYDAwC9gCZwL2edSdF5FlgneelphpjztforFSpyHHk8J9t/+Hd9e+yO3V3wfIuDbowpO0Qrm1xLQ2rFW6+almzJfd+cS9pWWllHa5SlrIkERhj7rjAegM8XMK6+cB8K+JQ6lKdzDzJB5s+4L0N75GamQpA7cja/KnDn7i1/a00rt64xH2rh1cH0ESg/J62cqlK6djZY7yx9g0+2/IZ2Y5sAH5X93cMjx3ODa1vIDgw+IKvER0eDWgiUP5PE4GqVFIyUnhjzRt8uOlDcp25APRt1pe/dv8r8Y3jEZGLfq0a4e6mLU0Eyt9pIlCVQmpmKvPWzuP9je8XlAAGtR7EiJ4jaFO7zWW9ZkRwBCGBIWQ7ssnKyyI8ONzKkJUqM5oIVIWW48jhn+v/yetrXic9Nx2AAS0HMKbXGNrVaefTa4sINcJrkJyeTFpWmiYC5bc0EagKyRjDwp0LmfHDDA6fOQzAVU2vYuyVY+lYr6Nl7+OdCBpUbWDZ6ypVljQRqApn/eH1TFsxjY1HNwLQplYbJl89mauaXmX5e2k7gaoINBGoCuPgqYPM+GEGi3ctBqBWRC3GXjmWW9rfQmBAYKm8Z40wdyI4maXdX5T/0kSg/F52XjZvrH2DN9e+Sa4zl9CgUIbHDuf+uPuJCokq1ffOLxGcyjpVqu+jVGnSRKD82vK9y5n67VQOnj4IwM3tbubRqx4ts/p6rRpSFYEmAuWXDp8+zLPLn2XZnmUAtKrZiqkDphLXOK5M49Dexaoi0ESg/EquM5f5ifN5JeEVsh3ZRAZHMqrXKIZ2HXpRvYGtVtC7OFsTgfJfmgiU31h7aC2PL3ucX0/+CsCgNoOYcvUU6lWpZ1tMWjWkKgJNBKrcO5tzlhd+eIGPNn0EQNMaTXmm/zNc2fRKmyPTqiFVMWgiUOXakt1LePqbpzmecZzggGAe7PEgD/V4iNCgULtDA3TgOVUxaCJQ5VJyejLPfPsMS3YvAaBL/S5Mu34arWu1tjmywrRqSFUEmghUueIyLj7d8inTv5/O2ZyzRAZH8uhVj/KXzn8ptU5hvtCB51RFYNUMZQOBOUAg8LYxZnqR9bOAfp6nEUAdY0x1zzonsMWz7qAx5iYrYlL+J+l0EhO/nsjqQ6sBuKb5NTwz4JlyPYaPDjynKgKfE4GIBAKvAdcCScA6EVlgjNmWv40x5hGv7UcCXbxeIssY09nXOJT/Msbw0aaPmP79dDLyMogOj+bJ/k8ypM2QS5ofwC468Jzyd1aUCOKAPcaYvQCeCepvBraVsP0duOc0VorDpw8zeelkVh5YCbjnCHh6wNPUjKhpb2CXQNsJlL+zIhE0BA55PU8CehS3oYg0AZoB33ktDhORRMABTDfG/LeEfe8H7geIiYmxIGxlJ2MMH2/+mOnfTyc9N50a4TV4pv8zDG472O7QLpkOPKf8XVk3Ft8OfG6McXota2KMOSwizYHvRGSLMebXojsaY+YB8wBiY2NN2YSrSsORM0eYvGQyPx34CYCBrQbyzIBnqBVZy+bILo8OPKf8nRWJ4DDQ2Ot5I8+y4twOPOy9wBhz2PN7r4iswN1+cE4iUP7PGMNnWz/j+eXPk56bTvWw6jw94Gm/aQsoiVYNKX9nRSJYB7QSkWa4E8DtwJ+LbiQibYEaQILXshpApjEmR0RqAb2BFy2ISZUzyenJTPp6Ej/s/wGA61pdx7MDnvXbUoA37V2s/J3PicAY4xCREcAS3LePzjfG/CIiU4FEY8wCz6a3Ax8bY7yrddoBb4qICwjA3UZQUiOz8lOLdy7miWVPcCr7FNXCqvF0/6e5se2Nfl0K8KYDzyl/Z0kbgTFmMbC4yLInizx/upj9VgEdrIhBlT9nss/wzHfP8N9t7vb/Pk378MLAF6gTVcfmyKylVUPK32nPYlUqEg4mMP6r8Rw9e5SwoDAeu/ox/tzpzxWmFOBNE4Hyd5oIlKVyHDm89ONLzF8/H4CO9Toyc9BMmkU3szmy0qOJQPk7TQTKMtuOb2PsorHsTt1NoAQyoucIHurxkC0TxpQlTQTK32kiUD5zupy8ve5tZq2cRZ4rj2Y1mjFz0Ew61u9od2hlIiI4gqCAILId2eQ4csrNENlKXSxNBMonR84cYdzicaxNWgvAXZ3vYmLfiZVq8DURITggGIfLgcPlIBRNBMq/aCJQl23xzsVMWTqFMzlnqB1ZmxcGvkDfZn3tDssWQYFB4ACHy2F3KEpdMk0E6pJl5Gbw7HfP8tnWzwDo36I/f7/+7341UJzVggLc/0qaCJQ/0kSgLsnmY5t5ZNEj7E/bT2hQKI/1fYy/dP5Lhbwt9FJoIlD+TBOBuigu4+KtdW8x86eZOFwO2tRqw+whs8vd1JF2KUgETk0Eyv9oIlAXdOzsMcZ/NZ5VB1cBMLTrUCb2mah3x3gpSARGE4HyP5oI1Hkt27OMyUsmk5aVRnR4NC/e8CL9mve78I6VjJYIlD/TRKCKlZWXxbQV0/hw04eAe5ygF294kdqRtW2OrHzSNgLlzzQRqHPsSNnB6C9Hs+fkHkICQ5jQZwJDuw4lQALsDq3c0kSg/JkmAlXAGMOHmz7kueXPkevMpWV0S2YPmU27Ou3sDq3cy08ETpfzAlsqVf5oIlCAe8joyUsm8/XurwG4reNtPNHviUrVQ9gX+Ykgz5VncyRKXTpNBIqNRzcy+svRJJ1JIiokiueve54hbYfYHZZfCQrUqiHlvyyp9BWRgSKyU0T2iMikYtbfIyIpIrLR8zPca91QEdnt+RlqRTzq4riMi3lr53HbR7eRdCaJDnU78OXdX2oSuAyBEghoIlD+yecSgYgEAq8B1wJJwDoRWVDMlJOfGGNGFNk3GngKiAUMsN6zr47nW8pSM1MZ/9V4vt/3PQD3dbuP8X3GExIYYnNk/il/qG1NBMofWVE1FAfsMcbsBRCRj4GbgYuZe/h6YJkx5qRn32XAQOAjC+JSJVh9cDVjF48lOT2Z6mHVmXHDDK5pcY3dYfk1bSxW/syKqqGGwCGv50meZUX9UUQ2i8jnItL4EvdFRO4XkUQRSUxJSbEg7MrH6XIyZ+Uc7vrsLpLTk+neqDsLhy7UJGCBwAB31ZA2Fit/VFY3hn8JNDXGdASWAf+61BcwxswzxsQaY2Jr19ZOTZfq2Nlj3PXpXcxNmIsxhhHxI3j/1vepX6W+3aFVCMEBnqoh7Vms/JAVVUOHgcZezxt5lhUwxqR6PX0beNFr36uL7LvCgpiUl+V7lzPhqwmczDpJ7cjazBw8k14xvewOq0LJLxE4jVYNKf9jRYlgHdBKRJqJSAhwO7DAewMR8b7svAnY7nm8BLhORGqISA3gOs8yZYFcZy7TVkxj+P8bzsmsk1zV9CoW3r1Qk0ApKOhH4NSqIeV/fC4RGGMcIjIC9wk8EJhvjPlFRKYCicaYBcAoEbkJcAAngXs8+54UkWdxJxOAqfkNx8o3h08fZuSXI9l0bBOBEsi4q8bx1+5/1WEiSkl+1ZCWCJQ/sqRDmTFmMbC4yLInvR5PBiaXsO98YL4VcSi37379jke/epTT2adpUKUBc26cQ9cGXe0Oq0LLrxrSNgLlj7RncQXicDmY+dNM3lz7JgDXNL+GGTfMoHp4dZsjq/jySwR615DyR5oIKohjZ48xZtEY1iWtI1ACefSqRxnefbhWBZWRgsZi7Ueg/JAmggpg5YGVjFk4hpNZJ6kbVZc5Q+bQvVF3u8OqVHTQOeXPNBH4MafLyWurX2PuqrkYDL2b9GbmoJnUiqxld2iVTv4QE1oiUP5IE4GfOpFxgnGLx/HTgZ8QhNG9RvNw/MMFVRSqbOmgc8qfaSLwQ4lJiYxaOIrk9GSiw6OZNXgWVza90u6wKrX8Yai1akj5I00EfsQYw1vr3uKlH1/CaZzENoxlzpA51KtSz+7QKj0ddE75M00EfuJU1ikmfD2Bb3/9FoAH4h5g7JVjC05Ayl4FcxZrPwLlh/Qs4gc2H93MiC9HcPjMYaqFVWPGDTPo36K/3WEpLwWJwGgLine truncated
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA4YklEQVR4nO3dd3xUVfr48c+THjqh996LoBEQBSyoCKj723WVXQu74lpWmiACsjZcWYrS1FXR1XXtbXURUEDBiij4pfeA9N4JmZSZeX5/zCQOIQHC3ORmkuf9euU1M7fMPNwZznPPufecI6qKMcaY0ivK7QCMMca4yxKBMcaUcpYIjDGmlLNEYIwxpZwlAmOMKeVi3A7gfFStWlUbNmzodhjGGBNRfv7554OqWi338ohMBA0bNmTp0qVuh2GMMRFFRLbltdyahowxppSzRGCMMaWcJQJjjCnlHEkEIvKqiOwXkdX5rBcRmS4iKSKyUkQuDFnXX0Q2Bf/6OxGPMcaYc+dUjeDfQK8zrL8OaBb8uxt4AUBEkoDHgM5AJ+AxEansUEzGGGPOgSOJQFW/AQ6fYZMbgf9owGKgkojUAq4F5qvqYVU9AsznzAnFGGOMw4rqGkEdYEfI653BZfktP42I3C0iS0Vk6YEDBwotUGOMKW0iph+Bqs4AZgAkJyfb2NmmSGX6Mvl518+s2ruK9QfWk+HN4IbWN3BVk6uIiYqY/0bG5KmofsG7gHohr+sGl+0CLs+1/KsiismYMzqZeZIvN3/J5xs/57ut33Ey6+Qp6z/f9Dk1y9VkUu9JdK3f1aUojQlfUSWCmcBAEXmXwIXhY6q6R0TmAuNCLhBfA4wuopiMOY2qsmTnEt5d+S7zNs3D4/XkrGtetTmd63amZfWWpGWm8daKt9h6ZCvvrXzPEoGJaI4kAhF5h8CZfVUR2UngTqBYAFV9EZgD9AZSgDTgz8F1h0XkSWBJ8K3GquqZLjobUygyvBl8vPZj/v3zv9l0aFPO8o61O9K3ZV+ubnI1dSqeevmqaZWm/PmjP3PEc6SowzXGUY4kAlX9w1nWK3B/PuteBV51Ig5jCupw2mHeWvEWbyx7g0NphwCoVrYav2/3e25uezP1KtXLd99KiZUALBGYiGdXuUyptPfEXl786UU+WPUB6d50ANrUaMNdyXdxXfPriI2OPet7JCUmAZYITOSzRGBKlQMnD/Dijy/y9oq3yfRlAtCjUQ/+cvFf6FKvCyJyzu9VOTFwacsSgYl0lghMqXAo7RAzfprBm8vfzKkB9G7em4GXDKRFtRbn9Z5lYssQFx1HujcdT5aHxNhEJ0M2pshYIjAlWoY3g9d+fo0XfnyB1MxUAHo27cnQrkNpVb1VWO8tIlROrMy+1H0c8RyxRGAiliUCUyKpKrM2zGLSN5PYdXwXAN0admPYZcNoX7O9Y58TmghqV6jt2PsaU5QsEZgS5+ddPzPuq3Es37McgBZVWzD68tF0a9jN8c+y6wSmJLBEYEqM7Ue3M+mbSczZOAeAqmWqMuyyYdzU9iaio6IL5TMrJwQSwWGPdX8xkcsSgYl46VnpvPjTi7z000tk+jKJj4nnruS7uLvT3ZSLK1eon51dIzjqOVqon2NMYbJEYCLawi0LGfvlWLYf2w7Aja1u5MFuDxZZe701DZmSwBKBiUi7ju3iyYVPMj9lPgDNqjRjbM+xdKrXqUjjsN7FpiSwRGAiSqYvk1eXvsqzPzxLujedsrFlGdx1MP0v7H9OvYGdltO7ON0SgYlclghMxPhpx0/8bf7f2Hx4MwC9W/RmzOVjqFm+pmsxWdOQKQksEZhi70TGCSZ8M4F3VrwDQMPKDXniqie4rOFlLkdmTUOmZLBEYIq1uZvm8vgXj7P/5H5io2K5t/O93Nf5PuJj4t0ODbCB50zJYInAFEv7UvfxxJdPMHfTXAA61urIuGvH0bxqc5cjO5U1DZmSwBKBKVb86uf9Ve8z/uvxnMg4QdnYsjzY7UFu7XBroXUKC4cNPGdKAqdmKOsFTAOigVdUdXyu9VOAK4IvywDVVbVScJ0PWBVct11Vb3AiJhN5dh7bycjPR7J4x2IArmx8JU/0fKJYj+FjA8+ZkiDsRCAi0cDzwNXATmCJiMxU1bXZ26jqAyHbDwI6hryFR1U7hBuHiVyqyjsr3mH81+M5mXWSpMQkHr3qUfq26Fug+QHcYgPPmUjnRI2gE5CiqlsAghPU3wiszWf7PxCY09gYdh3bxeh5o/l+2/dAYI6Ax3s+TpUyVVyO7NzZdQIT6ZxIBHWAHSGvdwKd89pQRBoAjYAFIYsTRGQp4AXGq+on+ex7N3A3QP369cOP2rhKVXl35buM/3o8qZmpVE6szBNXPUGfln3cDq3AbOA5E+mK+mJxP+BDVfWFLGugqrtEpDGwQERWqerm3Duq6gxgBkBycrIWTbimMOw+vpvRc0fz3bbvAOjVrBdP9HyCqmWruhzZ+bGB50ykcyIR7ALqhbyuG1yWl37A/aELVHVX8HGLiHxF4PrBaYnARD5V5YPVH/DUwqdIzUylUkIlHu/5eMRcC8iPNQ2ZSOdEIlgCNBORRgQSQD/gj7k3EpGWQGXgh5BllYE0Vc0QkarApcBEB2Iyxcy+1H2M+nwU32z9BoBrml3Dkz2fjNhaQCjrXWwiXdiJQFW9IjIQmEvg9tFXVXWNiIwFlqrqzOCm/YB3VTW0WacV8JKI+IEoAtcI8rvIbCLUnA1zeGT+IxxNP0rFhIo8ftXjXN/y+oiuBYSygedMpHPkGoGqzgHm5Fr2aK7Xj+ex3yKgnRMxmOLnePpxnljwBJ+s/QSA7g27M6HXBKqXq+5uYA6zpiET6axnsSkUP2z/gRGfjWDPiT0kxCTw8OUP88cL/lhiagGhLBGYSGeJwDgqw5vB098+zas/vwpA+5rtmdx7Mo2SGrkcWeGxRGAinSUC45i1+9cybPYwNh3aRLREM/CSgdzX+T5XJowpSpYITKSzRGDC5vP7eGXJK0z5fgpZ/iwaVW7E5N6TaV+rvduhFYkysWWIiYoh3ZtOhjej2AyRbcy5skRgwrL7+G6GzxnOTzt/AuD2DrczssfIUjX4mogQGxWL1+/F6/cSjyUCE1ksEZjzNmfDHMbMG8PxjONUK1uNCb0m0KNRD7fDckVMdAx4wev3uh2KMQVmicAU2MnMkzy54Ek+WP0BAFc1uYp/XPuPiBoozmkxUYH/SpYITCSyRGAKZOXelTww+wG2HtlKfEw8D/d4mFs73FoibwstCEsEJpJZIjDnxK9+Xl7yMpO/m4zX76VF1RZM7Tu12E0d6ZacROCzRGAijyUCc1Z7T+xlxGcjWLR9EQD9L+zPyO4j7e6YEDmJQC0RmMhjicCc0fyU+YyeO5ojniMkJSYx8bqJXNH4irPvWMpYjcBEMksEJk+eLA/jvhrH2yveBgLjBE28biLVylZzObLiya4RmEhmicCcZv2B9Qz5dAgph1OIi47joe4P0f/C/kRJlNuhFVuWCEwks0Rgcqgqb694m78v/DuZvkyaJjVlat+ptKreyu3Qir3sRODz+86ypTHFjyUCAwSGjB49dzSfb/ocgFva38IjVzxSqnoIhyM7EWT5s1yOxJiCs0RgWL5nOUM+HcLO4zspF1eOp655ir4t+7odVkSJibamIRO5HGn0FZFeIrJBRFJEZFQe6/8kIgdEZHnw766Qdf1FZFPwr78T8Zhz41c/M36awS3v3MLO4ztpV6Mdn97xqSWB8xAt0YAlAhOZwq4RiEg08DxwNbATWCIiM/OYcvI9VR2Ya98k4DEgGVDg5+C+Np5vITuUdogRn43g61++BuDOi+5kRPcRxEXHuRxZZMoeatsSgYlETjQNdQJSVHULgIi8C9wInMvcw9cC81X1cHDf+UAv4B0H4jL5WLx9McPmDGNf6j4qJVRi0nWTuLLJlW6HFdHsYrGJZE40DdUBdoS83hlcltvvRGSliHwoIvUKuC8icreILBWRpQcOHHAg7NLH5/cx7ftp3P7B7exL3cfFdS9mVv9ZlgQcEB0VaBqyi8UmEhXVjeGfAg1VtT0wH3i9oG+gqjNUNVlVk6tVs05NBbX3xF5uf/92pv8wHVVlYJeBvHnzm9QqX8vt0EqE2Khg05D1LDYRyImmoV1AvZDXdYPLcqjqoZCXrwATQ/a9PNe+XzkQkwmxcMtCHvrsIQ57DlOtbDUm95lM1/pd3Q6rRMmuEfjUmoZM5HGiRrAEaCYijUQkDugHzAzdQERCTztvANYFn88FrhGRyiJSGbgmuMw4INOXybivxnHXf+/isOcw3Rp2Y9YdsywJFIKcfgQ+axoykSfsGoGqekVkIIECPBp4VVXXiMhYYKmqzgQGi8gNgBc4DPwpuO9hEXmSQDIBGJt94diEZ9exXQz6dBAr9q4gWqIZ3m04f7n4LzZMRCHJbhqyGoGJRI50KFPVOcCcXMseDXk+Ghidz76vAq86EYcJWLB5AQ9+9iDH0o9Ru3xtpl0/jQtrX+h2WCVadtOQXSMwkch6FpcgXr+Xyd9N5qWfXgLgysZXMum6SVRKrORuYKVAdo3A7hoykcgSQQmx98Rehs4eypKdS4iWaB7s9iB3XXyXNQUVkZyLxdaPwEQgSwQlwPfbvmforKEc9hymRrkaTOs7jYvrXux2WKWKDTpnIpklggjm8/t4fvHzTF80HUW5tMGlTO49maplq7odWqmTPcSE1QhMJLJEEKEOnjzI8DnD+W7bdwjCkK5DuL/L/TlNFKZo2aBzJpJZIohAS3cuZfCswexL3UdSYhJT+kzhsoaXuR1WqZY9DLU1DZlIZIkggqgqLy95mae/fRqf+kiuk8y0vtOoWb6m26GVejbonIlklggixFHPUR76/CG+3PwlAPd0uodhlw3LKYCMu3LmLLZ+BCYCWSkSAVbuWcnATwey6/guKiZUZNJ1k7iqyVVuh2VC5CQCtURgIo8lgmJMVXlj2RuM+2oLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 432x288 with 1 Axes>"
|
||||
]
|
||||
@@ -233,32 +265,340 @@
|
||||
"fig = plt.figure().gca()\n",
|
||||
"\n",
|
||||
"# target constraint at t=0.5\n",
|
||||
"fig.plot(pltx, target.velocity.data.flatten(), lw=2, color='forestgreen') \n",
|
||||
"fig.plot(pltx, SOLUTION_T16.values.numpy('x'), lw=2, color='forestgreen') \n",
|
||||
"\n",
|
||||
"# constrained state of simulation\n",
|
||||
"contrained_state = sess.run(states[16], feed_dict={state_in: stateN}).velocity.data\n",
|
||||
"fig.plot(pltx, contrained_state.flatten(), lw=2, color='mediumblue')"
|
||||
"# optimized state of our simulation after 16 steps\n",
|
||||
"fig.plot(pltx, velocities[16].values.numpy('x'), lw=2, color='mediumblue')\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This seems to be going in the right direction! It's not there yet, but we've also only computed 5 steps so far. The two peaks with a positive velocity on the left side of the shock and the negative peak on the right side are starting to show.\n",
|
||||
"This seems to be going in the right direction! It's not there yet, but we've also only computed 5 GD update steps so far. The two peaks with a positive velocity on the left side of the shock and the negative peak on the right side are starting to show.\n",
|
||||
"\n",
|
||||
"This is a good indicator that the backpropagation of gradients through all of our 16 simulated steps is behaving correctly, and driving the solution in the right direction. This hints at how powerful this setup is: the gradient that we obtain from each of the simulation steps (and each operation within them) can easily be chained together into more complex sequences. Here, we're backpropagation through all 16 steps of the simulation, and we could easily enlarge this \"look-ahead\" of the optimization without any additional effort.\n",
|
||||
"This is a good indicator that the backpropagation of gradients through all of our 16 simulated steps is behaving correctly, and that it's driving the solution in the right direction. This hints at how powerful this setup is: the gradient that we obtain from each of the simulation steps (and each operation within them) can easily be chained together into more complex sequences. Here, we're backpropagation through all 16 steps of the simulation, and we could easily enlarge this \"look-ahead\" of the optimization with minor changes to the code.\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## More Optimization Steps\n",
|
||||
"\n",
|
||||
"A simple direct reconstruction problem like this one is always a good initial test for a DP solver, e.g., before moving to more complex setups like coupling it with an NN. If the direct optimization does not converge, there's probably still something fundamentally wrong, and there's no point involving an NN. \n",
|
||||
"\n",
|
||||
"Our Burgers example here looks good - before including an NN, let's address the performance issue so: here we've updated our solutoin via a gradient that is \"manually\" computed with multiple TF session run() calls. TF's own optimizers do this internally, more flexibly and more efficiently. Let's start over, and switch to a `tf.train.Optimizer()` instead."
|
||||
"Before moving on to more complex physics simulations, or involving NNs, let's finish the optimization task at hand, and run more steps to get a better solution.\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Optimization step 0, loss: 0.183535\n",
|
||||
"Optimization step 5, loss: 0.096258\n",
|
||||
"Optimization step 10, loss: 0.054813\n",
|
||||
"Optimization step 15, loss: 0.032847\n",
|
||||
"Optimization step 20, loss: 0.020363\n",
|
||||
"Optimization step 25, loss: 0.012899\n",
|
||||
"Optimization step 30, loss: 0.008259\n",
|
||||
"Optimization step 35, loss: 0.005283\n",
|
||||
"Optimization step 40, loss: 0.003360\n",
|
||||
"Runtime 20.50s\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"import time\n",
|
||||
"start = time.time()\n",
|
||||
"\n",
|
||||
"for optim_step in range(45):\n",
|
||||
" velocities = [velocity]\n",
|
||||
" with math.record_gradients(velocity.values):\n",
|
||||
" for time_step in range(STEPS):\n",
|
||||
" v1 = diffuse.explicit(1.0*velocities[-1], NU, DT)\n",
|
||||
" v2 = advect.semi_lagrangian(v1, v1, DT)\n",
|
||||
" velocities.append(v2)\n",
|
||||
"\n",
|
||||
" loss = field.l2_loss(velocities[16] - SOLUTION_T16)*2./N # MSE\n",
|
||||
" if optim_step%5==0: \n",
|
||||
" print('Optimization step %d, loss: %f' % (optim_step,loss))\n",
|
||||
"\n",
|
||||
" grad = math.gradients(loss, velocity.values)\n",
|
||||
"\n",
|
||||
" velocity = velocity - LR * grad\n",
|
||||
"\n",
|
||||
"end = time.time()\n",
|
||||
"print(\"Runtime {:.2f}s\".format(end-start))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Thinking back to the PINN version, we have a much lower error here after only 50 steps, and the oveall runtime is also much lower (ca. a factor of 4, depending on your hardware).\n",
|
||||
"\n",
|
||||
"Let's plot again how well our solution at $t=0.5$ (blue) matches the constraints (green) now:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[<matplotlib.lines.Line2D at 0x7fead6b6a1c0>]"
|
||||
]
|
||||
},
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABAm0lEQVR4nO3dd3gU1frA8e+bTYUASSB0EnqxgkbwooAUFRHFa0WviheQIlWkIwgoHWkiTUCxF7wqAopeEFARKT8RAZFmCAkhtBAgZZPdPb8/dsNdIKFlk80m7+d58rAzc2b2Zba8e86ZOUeMMSillCq+/LwdgFJKKe/SRKCUUsWcJgKllCrmNBEopVQxp4lAKaWKOX9vB3AtypUrZ6pXr+7tMJRSyqds3br1uDEm8sL1PpkIqlevzpYtW7wdhlJK+RQROZjTem0aUkqpYk4TgVJKFXOaCJRSqpjzSCIQkcUiclREduSyXURklojsE5HtInKL27ZOIrLX9dfJE/EopZS6cp6qEbwDtL3E9vuAOq6/bsBcABGJAF4BmgCNgVdEJNxDMSmllLoCHkkExpj1wMlLFOkAvGucNgJhIlIJuBf43hhz0hiTDHzPpROKUkopDyuoPoIqwCG35XjXutzWX0REuonIFhHZcuzYsXwLVCmlihufuY/AGLMAWAAQExOjY2erApVpz2Rrwlb+OPIHu4/txmqz8uB1D9K6Vmv8/XzmY6RUjgrqHZwAVHNbrupalwDcdcH6tQUUk1KXlJqZyur9q/l2z7f8FPsTqVmp523/du+3VAytyJR2U2ga1dRLUSqVdwWVCJYBvUXkY5wdwynGmEQRWQWMd+sgvgcYVkAxKXURYwyb4zfz8faP+W7vd6Tb0s9tq1uuLk2qNqF++fqkZabxwe8fEJscyyfbP9FEoHyaRxKBiHyE85d9ORGJx3klUACAMWYesBJoB+wD0oB/u7adFJFXgc2uQ401xlyq01mpfGG1Wfli1xe8s/Ud9p7Ye259o8qNaF+/PXfXupsqZc7vvqpdtjb//vzfJKcnF3S4SnmURxKBMebJy2w3QK9cti0GFnsiDqWu1sm0k3zw+we899t7nEg7AUBkyUgeu/ExHr/hcaqFVct137CQMABNBMrnaS+XKpaOnDnCvE3z+OyPz8iwZQBwfYXr6RrTlfvq3keAJQCA2FgrixadYOPGVP7v/9Jo0CCYkSMrcc89pYgIiQA0ESjfp4lAFSvHUo8x79d5fPj7h2TaMwFoUaMFz9/2PLdXux0RASA11c6ECUlMnZqE1fq/i9R+/jmVtm33cdddoXy0tBKgiUD5Pk0Eqlg4kXaCBZsW8P6298/VANrVbUfvf/SmXmS988ouX57CCy/EcehQFgAdO4bz+OPh3HxzCEuXJjN5chJr155lwZtnCSwVSIYtg/SsdEICQgr8/6WUJ2giUEWa1Wbl7a1vM/fXuZzNPAtAm9pt6N+0Pw3KNzivbFJSFn37HuLTT08BcOutJZg1qypNm4aeKzN4cEVuu60krVrtZcaMo9QeWYkTHCQ5PVkTgfJZmghUkWSMYflfy5myfgoJpxMAaFa9GQPuHMBNFW+6qOzbb59g4MAEkpPtlCjhx7hxlenTJxKLRS469l13hdKsWSg//niWlJ/vhdsXkJyeTOXSlQvk/6aUp2kiUEXO1oStjF87nm2J2wCoV64ew+4aRrPqzS4qu2tXOi+8cIh165y1hbZtSzN3bjWqVw/K9fgiwiuvVKRNm33ErrqTyg3f1X4C5dM0EagiI+5UHFPWT2HlnpUAlCtRjgF3DuDRGx7F4mc5r+zZs3ZeffUI06YlYbNBuXL+zJxZlSefDD/XYXwprVqVomnTkmzYAKc33s3JR/T2F+W7NBEon5eRlcG8TfOYv2k+mfZMgvyD6BrTlW6NuxEaGHpeWWMMX3xxiv794zl0KAsR6NGjHOPGVSYi4so/DiJC796RbNiQSsbfDTiVfsrD/yulCo4mAuXTfjjwA2NXjyUuJQ6ADg06MLDZwBzb6/fvt9KnzyG++eY04OwMnju3GrfdVvKanrtu3WAAbMmRJKcfukxppQovTQTKJyWkJPDqD6/y/b7vAahTtg5j24ylcbXGF5XNyHAwaVISEyYcwWo1lCljYfz4ynTvXi7HzuArVb16IABZyeU5mbb9mo+jlLdpIlA+JdOeyeIti3njlzfIsGVQMqAkfZv2pdMtnc7dDezu229T6N07nv37rQA8+2wEkydXoUKFi8terYgICyElHaSnluDIidTL76BUIaWJQPmMTYc28fL3L7P/5H4A2tVrx4i7RlCxVMWLyh46lMmLL8bz+eenALjuumDmzKlGixalPBaPiFCpquHAXxB/0OGx4ypV0DQRqELvjPUMk9ZP4qPfPwKgenh1xrQew53V77yobFaWYcaMo4wZk0hqqoOSJf0YPboS/fqVJyDg2puBclMt2p8DfxmSEiyXL6xUIaWJQBVqq/auYvR/R3M09SgBfgH0aNKDnk16EuR/8XX+69ef4YUXDrFzp3MIiUceCWP69KpUqxaYb/HVrB7MOtI5kZj7fQdKFXaaCFShlHQ2iTGrx7Bq7yoAGlVqxPh7x1O3XN2LyyZlMWhQAu+957yWv1atIGbPrkrbtmXyPc56tUOBdFKOXNuVR0oVBpoIVKHiMA4+/eNTJq6byBnrGUoGlGRgs4H8q+G/LropzG43zJ9/nOHDD5OSYicoSBg2rCJDhlQgONivQOKtUysUOIb1RIQOPKd8lqdmKGsLzAQswEJjzMQLtk8HWroWSwDljTFhrm124A/XtjhjzIOeiEn5nviUeIZ8O4SNhzYC0KpmK8a0GZPjPQGbN6fSs+chtm5NA5xDQ8yeXY1atQq2iaZGDefzZSWX14HnlM/KcyIQEQvwJnA3EA9sFpFlxphd2WWMMS+6le8DNHI7RLoxpmFe41C+yxjDR79/xMR1E0nNSiUiJIJRrUfRvl77i4Z7SE62MWLEYebNO44xULVqADNmVOXhh8OuaGgIT8u+l8CWHMnJtJM68JzySZ6oETQG9hljDgC4JqjvAOzKpfyTOOc0VoqElASGfTeMnw/+DDjnCBjdZjRlS5Q9r5wxhvffP8nAgQkcPWrD3x9efLE8o0ZVIjTUe1fshIf7E1Aig6y0EGIPn+KGi69kVarQ80QiqAK4318fDzTJqaCIRAM1gDVuq4NFZAtgAyYaY77MZd9uQDeAqKiovEetvMoYw8fbP2biuomczTxLeEg4Y1qP4f76919UdudO5wih69c7Rwi9886SzJ0bxQ03FI5mmDIVUjn+dzB/7T9D+1u8HY1SV6+gO4s7AkuNMXa3ddHGmAQRqQmsEZE/jDH7L9zRGLMAWAAQExNjLtyufMfh04cZtmoYPx38CYC2ddoyps0YypUsd1651FTnCKGvv+4cITQy0p8pU6rw7LMRXmkGyk1ExUyO/w2xf1u9HYpS18QTiSABqOa2XNW1LicdgV7uK4wxCa5/D4jIWpz9BxclAuX7jDF8tuMzxv0wjrOZZwkLDmN0m9EX9QUYY/jqqxT69j2UpxFCC0qFqnb2AHEHbd4ORalr4olP1WagjojUwJkAOgJPXVhIROoD4cAvbuvCgTRjjFVEygF3AJM9EJMqZJLOJjH026Gsj10PwD117uHVNq9eVAuIjbXSu/chVqxwjhB6yy0hzJ0bRePGhfc6/apRziR2+FDhqaUodTXynAiMMTYR6Q2swnn56GJjzE4RGQtsMcYscxXtCHxsjHFv1mkAzBcRB+CHs48gt05m5aNW/rWSkd+P5FTGKcoEl2F069E8UP+B82oBdrth1qyjvPxyImlpDkqXdk4X2bNnztNFFibZVw4dT8z7QHZKeYNH6tnGmJXAygvWjbpgeXQO+20AbvREDKrwOZ1xmjFrxvDlri8BaF69OZPaTqJ8aPnzym3fnkbXrnFs3uy8J+CJJ8KZPr0qlSr5xhdr7ZrO2krykcLRea3U1Sp8Da6qSPgl7hcGfTOIxDOJBPsHM/yu4Tx181Pn1QIyMhy89toRJk06gs3mvCdgzpxqPPBAmPcCvwb1apUGbJw5WgpjTKHqyFbqSmgiUB5ltVmZ+uNUFm9dDMBNFW9iWrtp1IiocV65devO0K1bHHv2WBGBXr0iGT++MqVL+94ontEVI8CSgCMzkIwMQ0iIJgLlWzQRKI/ZdXQXA1YMYO+JvVjEQu9/9KZnk57nTRiTkmJn8OAEFiw4DkCDBsG89VYUd9wRmtthC73wkHDEEoexB2Cz6ZXNyvdoIlB5ZnfYWbh5IdN/nk6WI4sa4TWY1m4aN1W66bxyX3xxil69DpGYmEVAgDB8eAWGDatIUFDBDBCXX0oElEAsNgxwNj2DUqUK7xVOSuVEE4HKk8OnD/PSypfYFL8JgGcaPsOQFkPOG3zt8OFM+vSJ5z//OQXA7beXZOHCKK6/vmh0rooI4uecoSwjU+8lUL5HE4G6Ziv/WsmI70Zw2nqayJKRTGo7iRY1Wpzbboxh4cITDBqUQEqKndBQPyZM8I1LQq+WWJyJwKqJQPkgTQTqqqVmpvLqmlf5bMdnALSu1ZoJ9044b6C4PXsy6NYtjnXrnOMD3X9/aebOjcrX2cK8KTsRaI1A+SJNBOqqbD+ynRdXvEhscixB/kEMbzGcfzX817lLJrOyDFOnJjFmTCJWqyEy0p9Zs6ryxBPhRfqySj8/rREo36WJQF0Rh3Hw1ua3mPbTNGwOG/XK1WNG+xnnTR25eXMqXbvGsX17OgCdOkXw+utVKVu26L/NxN85jmKGVROB8j1F/xOq8uzImSMM+mYQG+I2ANDplk4MaT7k3ATyqal2Ro1KZMaMozgcUKNGIPPnR3H33aW9GXaB8vNzXjZqzdJEoHyPJgJ1Sd/v+55hq4aRnJ5MREgEk++bTMuaLc9t/+6703TvHkdsbCZ+fjBwYHlGj65EyZK+d2NYXvj5u/oIrPbLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 432x288 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"fig = plt.figure().gca()\n",
|
||||
"fig.plot(pltx, SOLUTION_T16.values.numpy('x'), lw=2, color='forestgreen') \n",
|
||||
"fig.plot(pltx, velocities[16].values.numpy('x'), lw=2, color='mediumblue')\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Not bad. But how well is the initial state recovered via backpropagation through the 16 simulation steps? This is what we're changing, and because it's only indirectly constrained via the observation later in time there is more room to deviate from a desired or expected solution.\n",
|
||||
"\n",
|
||||
"This is shown in the next plot: "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[<matplotlib.lines.Line2D at 0x7fead71aadf0>]"
|
||||
]
|
||||
},
|
||||
"execution_count": 9,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABPtklEQVR4nO3de3yO9f/A8df7vnc0tpkNGzanOSbUQiQlSUpOnYh0cCahAzpHqZBDkpBzyrHwJRVFqZ/TlDNjmDFmYzabnXd/fn9c92Zsc2j3du/weT4ee7jv6/pc9/2+r831vq/PUZRSaJqmaaWXyd4BaJqmafalE4GmaVoppxOBpmlaKacTgaZpWimnE4GmaVop52DvAP4Lb29vVb16dXuHoWmaVqzs3r37glLK5/rtxTIRVK9eneDgYHuHoWmaVqyIyKnctuuqIU3TtFJOJwJN07RSTicCTdO0Uk4nAk3TtFJOJwJN07RSziaJQETmiUiUiBzIY7+IyBciEioi+0Tkrmz7+ojIMetPH1vEo2mapt06W90RLAA63GD/o0Cg9ac/MBNARLyA94HmQDPgfREpb6OYNE3TtFtgk3EESqk/RaT6DYp0BhYpY87r7SLiKSK+wAPARqVUDICIbMRIKN/bIi7t1iWkJrAvch/HLx7njy2JbP3Rh9YD/8LseoUyTmXwdPGkctnK1PGuQ12fuviW80VE7B22pmk2UFgDyqoAp7M9P2Pdltf2HESkP8bdBP7+/gUTZSmilCLkQgg/H/2Z34//zuHow1iUBYDI+W+SeKg6l3zXU+6eLbkeX82jGq2rt6Zd7XbcF3AfZpO5EKPXNM2Wis3IYqXUbGA2QFBQkF5N5z9KTE1kzeE1LP53MSEXQgBQChzNDtxR6Q7q+dRjSUZ9EoGHKjzPC52fJjE1kR1/pxESmkCZu3/jUNQhTsed5ru93/Hd3u+oVLYSXRt0pXfT3lQuV9m+H1DTtNtWWIkgAqiW7XlV67YIjOqh7Nu3FFJMpUpiaiKL/l3E7F2ziUuOA8DL1Yu6ST353/j7mTW7Ks8+VQmA+bH7gHTMl/15JLAGAMM7HiA8PJWDB3tQt54TB84f4I+Tf7Dm8BrCLoXx9c6vmbd7Hl0bdmVw88FU9ahqr4+qadptKqzuo2uB5629h1oAcUqpc8AvQHsRKW9tJG5v3abZiEVZWLF/BQ9+8yATt04kLjmOJr5NmPLYFP4e+DfVLj7N5Vhh0y/JAKSlKaKi0gE4cSIVgPj4DMLDjcf/939XMJvMNPZtzLCWw9j00iaW9VhGx7odSctIY9m+ZTw872EmbZ1EQmqCfT60pmm3xSZ3BCLyPcY3e28ROYPRE8gRQCn1NfAT0BEIBRKBF637YkRkHLDL+lJjMxuOtfw7duEY72x8h+AIY4K+OyvfyeutX6elf8usht5Tp4wLfOaFPjIyjcxlrE+cSDFe51hK1mtu336Fvn29s56LCEFVggiqEsSJVif4ctuXrDm8hpk7ZrLqwCo+bv8xbWu1LfDPqmnaf2erXkM9brJfAUPy2DcPmGeLODSDRVlY+M9CJvw5gdSMVCqUqcA7D75Dp3qdcvT0yUwAmQkhIiIta9/58+lcuZJBSEhy1rYdO67k+b41vWoy+bHJ9G7am3G/j2Nv5F76/diPrg268v5D71POuZwtP6amaTaiRxaXMBeuXOClVS/x0eaPSM1I5elGT7PxpY08Uf8JRISVKy9x/PjVb/iZiSA8PBWlFBERqde8XlhYKkePXi1/8GAyly9n3DCGpn5NWdFzBW898BbODs78eOhHnlj0BAfPH7ThJ9U0zVZ0IihBDpw/QJdvu7A1bCueLp7M7DyTTx75BA8XDwC2bk3gqadOMmhQOGC0B2TeASQnK6Kj0zl7Nu2a1zxxIoWjR6/eESgFu3YZdwXffHOB1atjc43FbDLzctDLrH9+PQ0rNiQ8Lpwnv3uS5fuX2/pja5qWTzoRlBDrj6znme+f4Vz8OZr6NmV9n/W0D2x/TZnMap09e5IAOHs2FYvl6v7w8NRrqobAaDAOCTHuCIKCygBGO8HWrQn06xfOU0+duOYO43o1vGqwoucKejTuQWpGKmN+GcP4LeOzxixommZ/OhGUAPOC5zFs3TCS05N56o6nWPLMklz78+/daySA6Oh0YmLSs9oFMmVPBA0bugDX3hE8/7wXYCSCjz+OBCA9Hd555+wN43N2cOajhz/i00c+xcHkwNzguQxeM5iktKR8fGpN02xFJ4JiTCnFp398ysdbPgZgTJsxfPLIJzg7OOdaft++qxfekJDkrPaBTKdOpWa1Edx3X1kAtm27Qny8BS8vM48/blQxbdoUzy+/XMbNzYSzs7B06SV27068abxPNXqK+U/Ox93ZnY2hG3lp1UvEp8Tf/gfXNM2mdCIopizKwvub3mfOrjk4mByY3HEyfe/pm+f8P6mpFg4fvlrXHxKSknVH4OxsHBMenpbVRtC6tZEIgoONC3ydOi5Ur+5ExYoOJCcb/UsHD/Zm6FBjHewRI84wadJ5evcO49dfL+cZd0v/lizvsZxKZSux88xOnl/xPLFJsfk4E5qm5ZdOBMWQRVl4d+O7LNm7BCezE7O7zqZzg843PObw4WTS0q7OzGHcERgX/ebN3QA4dSolq2rovvuMbZljCurUcUZEaNHC2O7sLIwcWYkxYyrj4WFm69YE3ngjgm+/jaF37zASEvLuWRToHcjSZ5dSzaMa+yL30XtF76zRzpqmFT6dCIoZpRQfbPqApfuW4uzgzOyus2lTo02OchkZio4dQ+nZ8yRKqaz2ARcX49v/kSNXq4YyL/oHDiSTkGChTBkT/v5OeHldnUiuTh2juqltW2MswIAB3lSu7EiFCg7MmFGNZs3K0LdvBRo1ciEqKp1p06Ju+Dn8Pf1Z+uxSqpevzqGoQ7y48kVdTaRpdqITQTEz+e/JWXcCc7rOoXX11rmW++OPBDZsuMz3319i796krETQsaNRz5+9aiizGihzBLGfnyMiQs2aV9sa6tY1Go8HD/Zh/fpaTJp0dS6h557zYseOesyZE8DUqcaUUhMnRhETk85338XQvv0x9uzJ2YZQuVxlFj+1mKruVdkbuZe+P/TVDciaZgc6ERQjc4Pn8tX2rzCLmemdptMqoFWeZZcvv5T1eNGimKyG4ief9AQgNPRqImjWzA1Hx6ttC1WqOAJQo4ZT1rbMOwJHR6FjR49rymfXtm052rUrR1xcBvfcc4Tnngtj48Z4Ro/OvWeRn7sf3z79LZXLVSY4Ipjh64eTbkm/2anQNM2GdCIoJtYfWc/4LeMB+KzDZ7Sr3S7PsunpilWrYrOef/ddTNbYgZYt3ahWzZG0NEViogV3dxNeXg5Uq+aYVT4zEWS/I6hd2+WWYx0/3g8wxiC4uZlwchJ+/fUyJ0/mPt6gmmc1FnRfgIeLB5tCN/HBbx+glJ5pXNMKi04ExcDuiN28vuF1AEa3GU3Xhl1vWH7LlnguXEinTh1nAgOdOX8+nQsX0vHwMOPv75RVzQPg7+90zb+QMxH4+ztRpsyt/6ncc48bH3zgS7dunuzZU4+nny6PUsZI5LwEegcyu8tsnMxOfL/3e2btnHXL76dpWv7oRFDEhceGM3D1QFIzUunZuCd9g/re9JjMaqFnnilP795eWdvvvNMVEbkmEQQE5EwEfn7G4zvuMMo1bux623G//74vq1bVpHZtFwYMMGYrnTfv4jU9l64XVDWIqY9NBWDS1klsDN142++radrt04mgCLuSeoUBqwcQkxTD/dXv5/2H3s9znEBkZBq//x7PuXNp/PBDLABPP12eXr2uJoLMC3q9elerfG50R3DvvW6sXVuTmTOzryl0+1q1cqN+fRciI9NZsOAi8+ZdYPz4SFJSck4z8UidR3i99esoFCPXj+Rw1OF8vbemaTdXbJaqLG2UUoz6eRRHLxylpldNvuj0BQ6m3H9dyckWWrc+Smjo1Tr4+vVdaNjQBRGhdeuybN2aQNOmRiLI7Y4g81+4mghEhE6dPPP9WUSE/v29GTHiDP37h2dtt1gU77zjm6P8wGYDCb0YyupDqxmwegBreq+hvGv5fMehaVru9B1BEfX1zq/ZcHQDZZ3KMqvLrBvO5T9lShShoSl4eppxczN+pQMGeGfdPcyZ48+771bmueeMu4NbbSOwpeef96JiRQecnIQ2bYzuqp99dp7z59NylBURxrcfz52V7yTicgQj1o8gw3Ljqa81TfvvbJIIRKSDiISISKiIjM5l/xQR2WP9OSoisdn2ZWTbt9YW8RR3u87sYvJfkwGY/NhkanrVzLPs2bOpWRPArVxZg7i4xpw/34hhw3yyytSt68LYsX44Oxu/7qpVHbMaf3NrI/D1tX0i8PJy4Pjxhly8eCdbttShUycPEhIsjB17Ltfyzg7OzHhiBl6uXmwN28r0bdNtHpOmaYZ8JwIRMQMzgEeBBkAPEWmQvYxSaoRSqolSqgkwHfgh2+6kzH1KqSfyG09xdynpEsPXD8eiLAxsNpCHaj10w/KjR5/lyhULXbp48NBD7pjNQsWKjnm2JQCYTEL79uWoUMFMw4ZGdVH16k54ezvQqJELTk4Fc6NYtqyZsmWN0cqffuqHyQSzZl24ZgW07Pzc/Zj6+FRMYmL6tulsDdtaIHFpWmlni//xzYBQpdQJpVQqsBS40cQ3PYDvbfC+JY5SitE/jyYyPpKmfk0Z3mp4nmUtFsVHH51j8eIYnJ2Fzz+vmmfZ3KxaVZPTpxvh4WFcmF1cTBw50oC//66bn49wyxo0cOXllyuQkQE9epzk3LmcVUQArQJa8WrLVwF4Y8MbXEy8WCjxaVppYotEUAU4ne35Geu2HEQkAKgB/J5ts4uIBIvIdhHpktebiEh/a7ng6OhoG4Rd9Cz8dyGbjm/C3dmdqY9NxdGcexVNTEw6jz12nHffPYcITJpU5ZrBX7fCZBJcXa/99Veo4EC5cuY8jrCLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 432x288 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"fig = plt.figure().gca()\n",
|
||||
"pltx = np.linspace(-1,1,N)\n",
|
||||
"\n",
|
||||
"# ground truth state at time=0 , move down\n",
|
||||
"INITIAL_GT = np.asarray( [-np.sin(np.pi * x) * 1. for x in np.linspace(-1,1,N)] ) # 1D array\n",
|
||||
"#t0gt = np.asarray( [ [-math.sin(np.pi * x) * 1.] for x in np.linspace(-1,1,n)] )\n",
|
||||
"fig.plot(pltx, INITIAL_GT.flatten() , lw=2, color='forestgreen') # ground truth initial state of sim\n",
|
||||
"\n",
|
||||
"fig.plot(pltx, velocity.values.numpy('x'), lw=2, color='mediumblue') # manual\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Naturally, this is a tougher task: the optimization receives direct feedback what the state at $t=0.5$ should look like, but due to the non-linear model equation, we typically have a large number of solution that exactly or numerically very closely satisfy the constraints. Hence, our minimizer not necessarily finds the exact state we started from. However, it's still quite close in this Burgers scenario.\n",
|
||||
"\n",
|
||||
"Before measuring the overall error of the reconstruction, let's visualize the full evolution of our system over time as this also yields the solution in the form of a numpy array that we can compare to the other versions:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1IAAAEeCAYAAACAIyQ5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABa/UlEQVR4nO29fZBtZ3Xe+VvvPt33SogvIwfLgjG4TMbGzgRiCnuGjE34SIjjsVw1DMFOHOzg0njKxM43Ip6yPcSk5GTKDlN2OVYBMSQOghA7aBISjPkY11QMQRgmGFy2BQYjRSAEEiDde7v7nHfNH/u01Fzuh/re9+nupX5+qlu3+/TpdV/ts88+77PXWs+KzMQYY4wxxhhjzEOnHfYCjDHGGGOMMaYaFlLGGGOMMcYYs08spIwxxhhjjDFmn1hIGWOMMcYYY8w+sZAyxhhjjDHGmH1iIWWMMcYYY4wx+2Rx2AswxhhjjDHGPHz4Cy/47/Jzd9+379/7wAf+6O2Z+QLBkiRYSBljjDHGGGOG8bm77+N9t/7Dff/eIv7q1YLlyLCQMsYYY4wxxgwjSXpfHfYy5FhIGWOMMcYYYwaSZC4PexFyLKSMMcYYY4wx40jIdEbKGGOMMcYYYx4ySdKdkTLGGGOMMcaY/eDSPmOMMcYYY4zZJxZSxhhjjDHGGLM/MsluIWWMMcYYY4wx+8MZKWOMMcYYY4zZDy7tM8YYY4wxxph9ktB3DnsRciykjDHGGGOMMcPIdEbKGGOMMcYYY/ZJgs0mjDHGGGOMMWYfpIWUMcYYY4wxxuwfl/YZY4wxxhhjzEMnSMIZKWOMMcYYY4zZB8ektK8d9gKMMcYYY4wxphrOSBljjDHGGGMGcjwyUhZSxhhjjDHGmIEkYbMJY4wxxhhjjNkHCfTVYa9CjoWUMcYYY4wxZiB27TPGGGOMMcaYfZLOSBljjDHGGGPMvjgm9ucWUsYYY4wxxpihhDNSxhhjjDHGGLMP0qV9xhhjjDHGGLNvnJEyxhhjjDHGmH3hjJQxxhhjjDHG7IvIPBYZqXbYCzDGGGOMMcY8zOir/f+5CBHxgoj4/Yi4LSJuOMfPfz4iPrT+8wcRce+en632/OyWEf+LzkgZY4wxxhhjxiHISEXEBPwi8HzgduD9EXFLZn70wX82/9ae5/8N4Ol7QpzOzKeNXJMzUsYYY4wxxpixjM9IPRO4LTM/npnbwM3AdRd4/vcBbxz0f3NOLKSMMcYYY4wxA0mi933/uQjXAp/a8/3t68e+goj4OuDJwLv2PHwyIm6NiPdGxPdexv/cA7i0zxhjjDHGGDOO5FJd+66OiFv3fH9TZt50CXFeDLwlM/cu4usy846I+HrgXRHx4cz82KUschcLKWOMMcYYY8xR4O7MfMZ5fnYH8MQ93z9h/di5eDHwo3sfyMw71n9/PCLew9w/dVlCyqV9xhhjjDHGmIGkokfq/cBTIuLJEbHJLJa+wn0vIr4ReCzw23see2xEnFh/fTXwLOCjZ//ufnFGyhhjjDHGGDOUyIv2PO2LzFxGxMuAtwMT8LrM/EhEvBK4NTN3RdWLgZszM/f8+jcBvxwRnTmRdONet79LJb783zDGGGOMMcaYS+dbv/GKfO9rv37fv7f5Zz/6gQuU9h05nJEyxhhjjDHGjOXiLnzlsZAyxhhjjDHGjCPTQsoYY4wxxhhj9ktcmv15KSykjDHGGGOMMQNxRsoYY4wxxhhj9kdiIWWMMcYYY4wx+8MZKWOMMcYYY4zZH8lDGbBbnqYKHBEviIjfj4jbIuIG1b9jjDHGGGOMOToESfS+7z/VkGSkImICfhF4PnA78P6IuGXEBGFjjDHGGGPMEaegMNovqtK+ZwK3ZebHASLiZuA6wELKGGOMMcaYhzM2m7gsrgU+tef724FvO9+TH7NxIr/mxBWipYwmDnsBR4I87AUYY4wxxhwDPr11ii/sbBXbgNpsQkpEXA9cD/A1J67gV771fxwav2et8y2Lrbfa8a1CtfPAaPD7y5gZXxONgR/58LsOewn7J4H+8L/trhJSdwBP3PP9E9aPPUBm3gTcBPDURz46F228s4dqMyK5sIs+K1IUeBLE9ObRmwYllc6vKR7+Hz4Xw+8FA0Cx90Kl64wxcpyRumTeDzwlIp7MLKBeDHz/+Z4ckWxs7AxdQKbGkLD3QuJMGFdBpQ8g2XEtJqgrEYXOr0ro3gu1NtCVqHStrYZvgtTad1TBR/ToIhFSmbmMiJcBb2dOXrwuMz9y3l8IWExjM1I9u+TN3JrmIqESfgpUYrLSRrfSB0XPIAp1tamObaUNTinhKzqu3uzrULwXKl0TjY6eQRS61tah4jF1j9RlkZlvA972UJ4bJK2NPdgNleDRfFj0HH+yqdYaocqe1ckiqj4oFK+ZSkBUy0RU2pgrhG+1ja7svK0kUgvx8N8uPfyQ3Gy2iNJQ8bLlHqmDQ1PaF5KNuSob00SiTyf8FHE1YrIJ9JmuFNOZSQWZUSsjJenDrPP/Dzrh6+yshlLvr5K7UgG+aWXUCJIER42jIaRasthYDo2pEhHZRXFF5YLVhJ8Ci8lamUkFlQQq1BOpCsplUhV4oytBJaZLnVtCqojqSoK6zkr3ks5IHSSLxXghpfiw6CuFXx30lWKj12lNlemqk+1TCd8yfQaRkvdCpfLGiPGuoKDbODVFGlVGlhJ+ivO22ga6zLWrGA//e++Hh+RzoVB2uiQu7Ts4IpJpuJBqhECcNNGmVIUuKzf+zVGp96qS6AOdkUeVEjRZmVikhR8wCe4v1coi1hKTRkOlM7YaFqlFBYmF1MEQkUyb4zNSCiGVq0kzQ2mw2QbosnLZW53yxh4oLsGtqUo86xhuQJ270FIhWaSMRemmZTGpEZMKdAJVVIZYSKBWen9V47hnUeus9EEyj0WL1NEQUkQyje6R6kHE+E+2nPp6cz6WCIHoyybpZeorzYfFSlE2KbOrF23Mu8a2v1K2T0E1IVnOvbFYFlFBlU1ZxKrMWkEnUKtcu2YE769CAlWFqqJAQpFlfgXOSB0MEUlTCKkmuPisGikQPRLRnkn4YlkqK0cTWV8LzoPMkGQNJCJiUm2cNCVd7uMpeAd68HqVJeQ+v6BSsZji2tUE+yMlsvmVRSoKSpJUeptdMkdCSBHQBvdIkUEXXCgiJlKSq9S8FAoTi9YSCgk0hZhMUa9cds1d0hRkpHRZuVqir05/UB3RB6pywTp3oKdCa4VjsV86NKplERXUMuAxD3AMLgxHQkhFS6YTg+dI9SBW41/BLmsu19whDMHFR5aNEZUMpqBkMDMkjcWqskmLvjqiD5CYuUC1Ek9NOVM14aegXOmogkLlqOAsYrWeydGUdRgsuuz9cCSEFJG00UJq2UhR6lqxeVKI9gBJn1hf1bpTmooPIIldvQ5VVk5BpVLMmEQ3VkTlnbIbNse+xBMq9bEc9025impZxGpi8lhT8bCmZr8cES8AXg1MwGsy88azfv6DwD8B7lg/9AuZ+Zr1z14C/O/rx38mM19/ues5EkJK0iPVGn17aEhgtjcdfmJkAwaXNsJc3jg+6nwMCl0su8ZnURBTZ5+r6OvLDMlFUnIMppVMoFXZlKpEX7XSjSqiD2qVM/Vi54HRnLetUM+R0snUrBl8XYg5O/CLwPOB24H3R8QtmfnRs576psx82Vm/+1XATwHPYL7z9YH1795zOWs6EkKKSGJwRkphfQ6QyxQM5e00wUuRXVN+lrLyRk0JngKFMAHNnjQzCMkupxFTDRGRXTRAWTVXTlQ2qaBaVq4WddSJTvTVyfSZWtS5yhZm/HvtmcBtmflxgIi4GbgOOFtInYu/ALwjMz+//t13AC8A3ng5CzoaQqolsTn4LmFfSe6M5M6CWI2/o9lDMUeqkUvNLK3jXt6YqbHBV6ERJyKrdkkZoujGSsaxL5sETb+cKiunoF4P3vGe++USz5kq2XQVlcxnSvZIaUr7rgU+tef724FvO8fz/ueI+A7gD4C/lZmfOs/vXnu5CzoaQiogNgZvoztIyuUABHfeQrARiZ6y+5kKE4tK5Y30kMz+UqH4sOhMmqyBwhAhVU6TzsplBkyCm0uFDFJkiJbrbJ9LPFW4xNM8SFzq5+7VEXHrnu9vysyb9vH7/zfwxszcioj/FXg98JxLWchD4WgIKSAGryT7LCSGs1iBwMSi9fEiIldN4ywHkiFrkv4zUlTeqDEzkaw1Q5Lua9NKszEVhAyynEAbTbUeqUq9cpVKMWWoDoHo/KrksugSTx1lhF+t+xSXy92Z+Yzz/OwO4Il7vn8CD5pKAJCZn9vz7WuAf7znd5991u++53IWCkdFSDXgxNizJFZJCu66BX3ePI1G4WyyCLrgwq6YI7UroEZX4cnukkYSkiyPAFH2LDM0joiqTITqroKiXEx0jZFoyRCNQygm0BQ4K4ezcuWookxmqgm/coy/hr0feEpEPJlZGL0Y+P69T4iIazLzzvW33wP83vrrtwP/KCIeu/7+zwOvuNwFHQ0hFQGbg5fSkxCUiiWpyXQJMlL0EPWJ5fAPjMggl5oPNkXmSDGbCpjtaEdfeCYEBinMG3PFMFrF3C9RLbyy92p4zAiRKcSqXF/b8JiiDXSlrJxKQsjm1RVyxlT1x2isqWv1IVYRfu6RWofMXEbEy5hF0QS8LjM/EhGvBG7NzFuAH4uI72HuF/k88IPr3/18RPxDZjEG8Mpd44nL4QgJqY2xMZcrSc426JrWK4WIaPMmZ3zgS657PXiEd18lH0ItUDQry8obJR9u48/ZEK1VJdAU561qk5c9NFk5TZIHl01qhueqNrq9iQ5CtcycgkqHoIbeMWcjMSLKtwFvO+uxn9zz9Ss4T6YpM18HvG7kei5ZSEXEE4E3AI9n3vndlJmvXvu0vwl4EvAJ4EUX9WiPgMVgTdeapvi1J5JRzZs53m+/p6RDKFjCUpHhELxeyhkRVUoCsqH6FFJsnRRzv5IuMXShae4UyvYLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 1152x360 with 2 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"import pylab\n",
|
||||
"\n",
|
||||
"def show_state(a):\n",
|
||||
" a=np.expand_dims(a, axis=2)\n",
|
||||
" for i in range(4):\n",
|
||||
" a = np.concatenate( [a,a] , axis=2)\n",
|
||||
" a = np.reshape( a, [a.shape[0],a.shape[1]*a.shape[2]] )\n",
|
||||
" fig, axes = pylab.subplots(1, 1, figsize=(16, 5))\n",
|
||||
" im = axes.imshow(a, origin='upper', cmap='inferno')\n",
|
||||
" pylab.colorbar(im) \n",
|
||||
" \n",
|
||||
"# get numpy versions of all states \n",
|
||||
"vels = [ x.values.numpy('x,vector') for x in velocities] \n",
|
||||
"# concatenate along vector/features dimension\n",
|
||||
"vels = np.concatenate(vels, axis=-1) \n",
|
||||
"\n",
|
||||
"# save for comparison with other methods\n",
|
||||
"np.savez_compressed(\"./temp/burgers-diffphys-solution.npz\", np.reshape(vels,[N,STEPS+1])) # remove batch & channel dimension\n",
|
||||
"\n",
|
||||
"show_state(vels)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Physics-Informed vs. Differentiable Physics Reconstruction\n",
|
||||
"\n",
|
||||
"Now we have both versions, the one with the PINN, and the DP version, so let's compare both reconstructions in more detail.\n",
|
||||
"\n",
|
||||
"Let's first look at the solutions side by side. The code below generates an image with 3 versions, from top to bottom: the \"ground truth\" (GT) solution as given by the regular forward simulation, in the middle the PINN reconstruction, and at the bottom the differentiable physics version.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"\n",
|
||||
"Solutions Ground Truth (top), PINN (middle) , DiffPhys (bottom):\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAesAAAEyCAYAAADawvhGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABuyklEQVR4nO29e5Rt11Xe+Ztr71N1r64kS7aMELbBJpjmmdhBwyZNRmIMTgyhUdINRpCAnZgoMOwmD5K2DD2ANnEPQRKIe8BwUGwHmwCyMQErxLTjB24GIxgsg4NfBISxsYRs+aHHfVXV2XvN/mPvK5Wub93zzauzzt67an8aGrfq1DxrrbPPOfub85tzzWXuzowZM2bMmDFjvEhDL2DGjBkzZsyYcXHMZD1jxowZM2aMHDNZz5gxY8aMGSPHTNYzZsyYMWPGyDGT9YwZM2bMmDFyzGQ9Y8aMGTNmjBzFyNrMnmtm/8PM7jSzm0vNM2PGjBkzZhx2WIl91mZWAX8EPAe4C3g38O3u/sG1TzZjxowZM2YccpSKrJ8B3OnuH3b3PeA24IZCc82YMWPGjBmHGnWhcZ8AfGzf73cBz9xvYGY3ATcBHEvVV33+8SsKLWXGjGnAsaGXMOOI4xO7p3lguTupD+LffO5f9E9/6lT4ee95z5++xd2fW2BJRVCKrFfC3W8FbgX40ise46/9qv9ZfJ7+OSpx88uB+SG23nUjutYhMeR1KokpvQczDu/nUMX3vO8dQy8hjE9/6hS/c8ePhp9X29+7psByiqEUWd8NPGnf70/sH7sgqpS54sRpaeCcNeU+u8WIXbR11zMHOa9//qitihKkUurGN3QEOgYCPqykMoZrOyWs+3MwxavvODm3Qy+jOEqR9buBp5rZU+hI+kbgOw4yTlXmsstXk7XnJJO1Z52sc64kO4C21cladwAsdJMaktgjzgrEHBZ9DcM6NhGUIp+hXxeUcZr0b2IMR90JkD8vNsWDnRz3ZuhFFEcRsnb3xsxeDLyF7vv3Gnf/wEH2KWWOX35GGrtttK+zt5VO1m2Sv8y51W8nOUDsqhMCOmGWiOxLEeVUnJWHx11/bWYJxwZKXoPpOEJpYK4e+lqpftUkXRoH9zmyvmS4+5uBNyu2VmW2rxQKBNxod7fE+U0my9xWuHijzI1+ydT5I5G15xQiVvV1lSLAIR2Lkrag3RxChFKIUUoQRQlnBWKfgxIrGJpUh57fmF5k7Th5jqw3A6tatq5aTdaeE9XunjSmtxV5qUXBIbJeapescxZEFSAg2ZNNJ8BAZN9GIns1FRGtG5iIYwE6qVTBNagYWl3Intc+JpRxWKakrkRQRImZaGg9y+AbglWZ6srVMrg3FWmhvSneJvKeSKxthYvEmpfi/DnpZBmw7aJlnazVG1UlE2Ugso84IehRqPr6ISbve9YJaGhSGbrQMUJTU5LhY7lt/fNSYq0plSgwm15kPZP1JlFBeozwwc8ZP6OStZFUsl7WuCqZ79XyjdLVKDybPn9ABVAdkG5ctXCvEFkHFINotC7XLoTGLbDW0PyzYxHB7IQc4qS1O55nst4MKoMrjq22y45Vmgxuy4wvRMl82eBi4ZotGlAj292FNn82mVhTwNbbJJOKRarcC8j7USdAsnODUFQZUTfWH1WVcBYi80dt5TELORZjyMUPvdtBRRLHnWZkDcyR9YZQ1fCYq1bbNQ3UWtU4TYvtBIi90W4ovtvK6pfVWhGSZ5OdBW8rvBFJJRCFp0gULioGEMibFyDrbv4y1fsRxUAfc1hnAQo5DAH+K6dEqLax65rE1zYGx0LCFCPrWQbfIFIiX37larucSeq3o2mgFm/Ue0us0YjV6ix/n11UAdwN3xPX2qYAseuRtTdJ/vJbJW6T8ISJjgWBXHykyj1VohPmRqoi+e0CefMIAYccC902FZDiSzkLKjwn2WEYXoWI5cznnDWAQ14OvYjiGAdZVzX58qtW2+WM6sra3i62pW3zYmcHWtEzq3chax9oS9qYlh1XSaVNWCPm7ZsKVAlWzO9DLBdvC91ZQI0sc2BLWiCyDkX2BRyLEFFUA0fh2QIEWMpZEG3T+kkNyqkQg+4emGBk7T5H1huDp4r2xOMEwyZA1jvYjiaZW73AGtEzS6lzGiTbXc2uabuIXYA1WU7PWKPfpFKVZcnWRXmfnOQI0KtA4Zwnfa2tTlS5kAqgQlVMunFF5y5YkV/kdU3IWYi8LoukOAp0BivhLNhEO5gxF5htCKkmX3b1artA/9dUb2G1FlmnndO4SNaWKkxdhyrZ5ww7AWJvxC9U6yDaRr6kvtQrx9XCNV/oRX4hshTTG+6JSpXBsxVRAajatZNlRAXoxtWdAH3M9RcZUulrsIBcbKDXGRxxZ2E08JmsNwa3Ct++arVhbuTyD69Ok8RiNK8XpL0dbeCU5EhB/ooG5H2aBkQCYm+JLcQvae16bU0S8/s54614k8wWSwVIll3ErqYC5FqAwDYzi5C1658tfcxS+fUAUaipiCoXcwJCY6o+4wichXU7AtOMrJmrwTcFSzW2tTqydm/k8gevt2FvW5y/wtUoPBBZq9xnOWMhslbz65Uu2YuV82THkvguZD2yJ9M5DAKs0W/S8lpBdxayrd0BgK5yXXZCVMmeYN6+gBPgkfdAnT+oFhRxAqzQ9jkLJI5dv7bSqBPMWRuOzZH1hmCJVF8umcpCeHOWbNrLs1RjjShDp0qS4y23cmTtzZ5M1tYssT11rZH8esB2T/xiNC0mEjAZELfPUSFFH4Ye2ZNNPnHIAkVjFogWU4TUIkRRgixLqAAEo/BAfttUf2UMKQM5si/hLEwwsp5l8M3BSFQCWbs3uGuk7ukYnsSmJKnGa1EGh/V/MPIxUtLO86ZZ4rX2uqzeLVQ4J0bhTatL9k2r31Brlyvy2QvcfEKOhbjYnDE5Cg+QqqrEBMb1AKl5W+m2kRx/EceijAqgytVWEXQC1Mg+cK1UGX6CkfVRwSjIGktUabVk7V7LK/bckJNmnFMNzVnNFmTJRVUBrN3tInbFdm8Ha0SyTAnPaipAe/2WWz2/vsiwDBD7uiN7CEj26JF9BqtEAmwDXdQCFfnyewBBxUBUeAKk1jWmEXc7qFX+2fTUEXpFemif84ScgEMdWc/V4JuDYSRFshZlbYBsgTev7iJxBQ7r39PX7OqSfb2F7WnEmkTJPoScdceiWeqkUjdlJPvU6vviA5G1Gq3KlfvQVe+rhQ6BXLhKrJ4NU2/+gag2JO+LuXirAsVV2eSQMeYAGBapc5AMrUyKQ51/kpG1Y3OB2figROAA5pXmAABtruVG8K34oXBv5DG9OoaJKoA1u5hIlmp+HfTKdc9toBiulovhrFnqhXOh/LqYBshZl+wjRXY1AQLWu+OFTpyQFQP9uoby9gEHwJJKfglTUyGBbV4lHADQnYCoAyCvQFQBJlkN7qw/KBkhJkfWKpLV8n2vSsfkSFzOmbu+zYysV7mTKjmyzVaHtjTI+8dV1Au5GM7rALHv7epkHUGE2FWiaAK95FMksg68fpVUPDB/ZP92K6oQ6PvCI4V7uL7WqAOgk7t4DSL7x0uoAFMk67kafHNwnCyQihopR+0zUKnVMvX6ZXAPkLW3Z/FKK4aLVLl3ufh1R+H6ljRrliAWzkUia/mGFtnrnnNABbAYsUvzx3LGqlRqOUjW6rahiGPRVF3uXkGgP3pI3perzBOmXrBALl6vRwjUOKhpAHHmccHnyHpzcNxXX+zWW1kGjyDiBLi1cu7cxd7gboFmL6nW56ffby4gIZaWBDxYa/YCzV6WoBbO0e1NVxCJE+QtcRFiT5FcvEoSrhM7wSI79ZCaxrVDHzIyUbib3E+/azcayNeuO7cLoboB/b0NyPsE9mSrn9cpRtbz1q3NwT3TNKdW2ql5XQAzPWcdgeosZG+oxb3j2RtMVQFyjQfIWpXBI1XucrOXahtLYn/2VGFyZF3hoictR/Y5y8RuuQW1KUlV64fEqEQFumQPoVy8LFkndGIXI3bDQWx4Q+V6ZJ/RSTjrfecjhXuyDD+0vD/N0Hr9KbwRYhRkjWdyFqTdjLxis1rWa0uQeihnLm8w7qCSdSYm2asFcbLklBtcrRyvd+UtaQnWn7POWe8eltsyFfFqFbDYjAOIXafsIBZ4EZDMTXQAIFLdHCEgU7d5dzd9WbLWK9dVJwiIORZqJKxG9lPcuuWzDL4xuLfkvftXGwYi6yRuxYJO2jaRMCPEHnICAo5Fqzg2PVSyLlEM57nB5GYvuryvt3ENfIFzKzcaCeXXm1rPr0dk+Aixy9JysHBORYlcfESyz64Te2Nd1K4gQsCRKndZ2vaAgykW7s2R9WgxCrLGW9i9b7Vd0qNVT41MVF4dk2VoUrlIXEGmq15XIUfLtWibAr53e7ZIlbsK9ya4JU3fEqd2hvMUOH8cyhTORYg9kl8P5FZlqMRuyHloixB7oCJfdyyGjeyBzgFQ1jvFnPVcYLY5WG5IZz+z2jBVZHXvcn2MdnGFZuuXhyNxBcWK4QpE7E2j2Ua2pIWK4VLd2QuwVK89vw6QxGYzHjj/vCuyC7R8VcbMWW83GjkkJlTlHmlMEwjXIkV26vylJHs5w+R6cxyx7303bLB6X8EEydrc58h6Y8gN6awSWQdOD2p35Kgy5wavVLK+Svp+WKplaR3KROuh+eXXr8NToBguLeQ0hwOOWOUe+BLL8nqgap2U9M9tSlKVuxOrAwrJ+5HT39SPV4TYVWQP5HYDxB6AVA3fI0SB6p7s7HqlvRqET1QGnyPrFTCzjwAn6dpgN+5+vZk9Fng98GTgI8Dz3P2iTGy5pTolkrX4pni9hW+Le4zbHZmsW9AjW+3UzWLFcCVy5u6tfK08cP44pkfWEchpgAislgk4Iq9bifawENhmFty/HiHACLEriKgAEWJPAWIPbIcuk4snRuyHFYUiazN7LvAKuk/vq9z9lvP+/pPA1/a/XgZ8jrtf1f+tBd7X/+3Line truncated
|
||||
"text/plain": [
|
||||
"<Figure size 1152x360 with 2 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# note, this requires running the forward sim & PINN notebooks beforehand\n",
|
||||
"sol_gt=npfile=np.load(\"./temp/burgers-groundtruth-solution.npz\")[\"arr_0\"] ; #print(format(sol_gt.shape)) \n",
|
||||
"sol_pi=npfile=np.load(\"./temp/burgers-pinn-solution.npz\")[\"arr_0\"] ; #print(format(sol_pi.shape)) \n",
|
||||
"sol_dp=npfile=np.load(\"./temp/burgers-diffphys-solution.npz\")[\"arr_0\"] ; #print(format(sol_dp.shape)) \n",
|
||||
"\n",
|
||||
"divider = np.ones([10,33])*-1. # we'll sneak in a block of -1s to show a black divider in the image\n",
|
||||
"sbs = np.concatenate( [sol_gt, divider, sol_pi, divider, sol_dp], axis=0)\n",
|
||||
"\n",
|
||||
"print(\"\\nSolutions Ground Truth (top), PINN (middle) , DiffPhys (bottom):\")\n",
|
||||
"show_state(np.reshape(sbs,[N*3+20,33,1]))\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"It's quite clearly visible here that the PINN solution (in the middle) recovers the overall shape of the solution, hence the temporal constraints are at least partially fulfilled. However, it doesn't manage to capture the amplitudes of the GT solution very well.\n",
|
||||
"\n",
|
||||
"The reconstruction from the optimization with a differentiable solver (at the bottom) is much closer to the ground truth thanks to an improved flow of gradients over the whole course of the sequence. In addition, it can leverage the grid-based discretization for both forwards as well as backwards passes, and in this way provide a more accurate signal to the unknown initial state. It is nonetheless visible, that the reconstruction lacks certain \"sharper\" features of the GT version, e.g., visible in the bottom left corner of the solution image.\n",
|
||||
"\n",
|
||||
"Let's quantify these errors over the whole sequence:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"MAE PINN: \t0.21460 \n",
|
||||
"MAE DP: \t0.05909\n",
|
||||
"\n",
|
||||
"Error GT to PINN (top) , DiffPhys (bottom):\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAqMAAAEyCAYAAADOTdkTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABYKElEQVR4nO29fdhtZ13f+fmttfdzTt4JBENIUhPHtIpawzRXsINX5VWj40V0xqaBUUFDo3PBWFvtEGQudKDpdbQtDHNpLadABauGFEUyNp3Ii4zTUTRBmQFCqxGD5DQQEgIkOec8z97r/s0fa51k53Be9/d+9n3vtX8frnCeZz97rX3v9fpd39/Lbe5OEARBEARBEJSgKT2AIAiCIAiCYHMJMRoEQRAEQRAUI8RoEARBEARBUIwQo0EQBEEQBEExQowGQRAEQRAExQgxGgRBEARBEBRj18SomV1jZv/ZzO4xs5t263OCIAiCIAiC9cV2o8+ombXAnwEvBu4D7gRe6u53Z/+wIAiCIAiCYG3ZLWf0auAed/+0u+8AtwDX7tJnBUEQBEEQBGvKZJfWezHw2YXf7wOec7w3n7+15c88Y+8uDeUUEQ1ix/KMQxtEDasQKb8dy2+DDIziS5Q/FlRGsRuCYET82WNfetDdn156HKfKd13zN/2hBx897eU++tG/vMPdr9mFIe0KuyVGT4qZ3QjcCPCMPXv51av+G2l97mVvXOrnpwzjL74NMoiHHNtBIcc2LP0dclD6WIJxbMcaqGFfBuUxi0cjgOf/4Xs/U3oMp8NDDz7KH931xtNebmI/eMEuDGfX2C0xegC4dOH3S4bXHsfd9wP7Ab7xnKf4bF5MF2ehBme09M27hpte6W0AsR1qoYb9EAS1MJbzYdNEteOk1JUexq6zWwrwTuAKM7ucXoReD7zseG92h52urBhtKjjAQ0DUccEcw36oYTuuO7ENg7ExBiG3eeel4z4vPYhdZ1cUoLvPzezVwB1AC7zD3T953PdjKM5ojhNsDCdpoLN5F7pjE9uhHlIFUZd1p4nsXWAc5/XG3asd3MMZXRp3vx24/dTea+x07dKfZaY7m+oBHhe7YEyEAArGRBzPPWO4T41BUJ8OjpPCGV0NyU0K07dNwsSTTBez0uJZ0gQ27okxCE7Apt20guMT18aedRflYxDTp0+E6VeGA7O0fMtT9QRrcPli1VqSlu+wDIK4/IWmhtzbIAiCRcbyYLLponrdxfRyhBhdGWqYvnWXnpjMnLbRTnJ3ky8ULj71mTibVg4h2YkX/Qr0dBAEGYkH1Hysu6jedDG9FO54CjG6EhyYueCMJlWMgrvmbLZNAuFC0Y+hbN6rLIYzhFDGEIYJQR0ETzCG7hQQojoHOcT0RgracEZXg6M5o9MmSQdo45KOzIKZy2NYfgv2aHIcyJBqIDurFYjZMQjqMRAPBUFOahDVIYjX3x0+fSJMvzLctZxR0J6Wmsf/T/h88RrR4LIzSqPJScuQe6uH6dVUA2nxfgyimIybVh2op1PQE6K+HkpfW+K6UgKHNCs9iF2nDjEKzJSTLDXSSWLmIIphldYSSb3qi9ZmK4rZDtO7CojXOt3dzSNoFXK4u6VvWhA3rrFQg6gPQVwH6nUlrgmnj3s4oyvDMeaKGGyS9kXcaEQh1sg3/yaDCBEFtajkzLTcXUAW5KqYBV3QqhfcHFPL1pCuUFoQx41vPJQWxCGG85DjmrB557VDFDCtBnfYTspB2uCiM2oF0wQeX4e4fA1CUBZypb8DLuck5XFnywraGsSsSmkxDJt44xwnOcRwCNo81HBerxQPMboyHK1wZZZA0bKN6c6omisJurvaqJkGSdsGySxD8Y7oUBcWs6C7s2ZePBw2Fne2NKX3Y1APqqANMbvBRJh+NfRN75c/05rGpWa4jesN45vUaLfvJulBerl4SHVGtf0AfbKCQg5ntvRFX3VmaxCzEO5sDmrYj0EQlMNwLJzR1eBuzBUxaia1NXKcTkoTgLmoYDw1tOoMTKqIEZd3jFaUg6VzRhPl3VX183NMwFCDCAoxq5Ncb7cW1EE4qxtKhOlXR19Nv/zylkwM0+utmdoK8ljUMXRiV4JkDsLkBf1KRG9T/Pgx5L1CnlQBhRpEUIjZnhoeLILyuIcgDeqlCjEK0AkKoDMD4cbhbnTiSao6q6qIAl2MquHdxq14r1P5vmnl20tFqkBPaRHk2GgEaRAE60o4oyvDHebCjat1UURZvw6FmZumJ/v4sITccF51dzP0a1ULydQpTRMmtzWIVIFIFThCiNkgCDQciwKm1dBX0y+//CyZlG+pFlABtAadtAZkS03NW1XbW9Ek2ZrUuxI0mgjKoB1yiCBZiFUgZnO4q0EQBBuNA0lWF9VThRgFTYz2TtTyNz53H0L9y9OJzmovA9UxqAJCE1JdauQwvZqtYOaguORkcFcraNyv5ytoJHK0+RLHEM5qFmrYjkGwuUQ1/cpQnVE9LGpymF7NGbUGVFtOHUPbiukOUDxMr7a3ylGEZa5OZaUt3pAhX1M8Fs3Q23xFEVYQBBuPhzO6KlQxamK+paEVUAGys0rSNJABrXjflEPk5rKzqRb/zNXc3QxFWDkmQJCIIqxslHYFo4gqyMUYzseNJFo7rRCHpKT5GcyFu2djWgEVwMz1gJyiAMz0VIEaOgKoeatNk6TcXTOX3VVVOtRQhNWqY3A9TF9DukMNeavrHuoPhzkINCyc0dWQI0yvXK7Vz4feWdWfPJdfgeGyoFad1RwdASYZCqAUAdG4oT6DqgJGbW9Vw0xYWVo7ZRCzkbdaB7EdgmBJPML0K0VxRpP37qayvCpG5fCweNNsMF1Qq1N8iO4uZOgIAFIBk5tLs3mB3qOzQ9P03QhmwsrlSK573moO1t1ZDfIQTe/Xl91yRs3sGuAtQAu8zd33HfX3NwPPH349E/gad3/K8LcO+Pjwt79y95coY6lGjCp4jllzxOut7q6Kjfvlz9dTFfpRlO8IoODoAkLdjhOGQiqBKiYPECeiKD0TVg1FWOo2gHG4iuGsBpvL7jijZtYCvwS8GLgPuNPMbnP3ux//ZPd/uPD+/wl49sIqDrn7lbnGU40YVfLkEiYZiw2aMwu6EOwp3Z5KWvzISLQxiHmrTSOKSbe+RZU0BtGVHMax9OebmnVaR95qlskDlP7D4SQB4ayOiZjffv0w991yRq8G7nH3TwOY2S3AtcDdx3n/S4Gf3Y2BQCVi1NHEYI5UTfVyqRRQQd/aSfkeWdpTyfcM1b3Qx5Ajf1hFLYACTVCrLlL/+eXzVnMUYdWQtyp9foYDUnaYM3QEUMghZsNZzUOI2ULsjhi9GPjswu/3Ac851hvN7GuBy4EPLby818zuAubAPnf/bWUwVYhRFR8aoCxLcr21UysmGuZwJYu3p8Kl3N3O+5msFJIj5ozm6U+pkKOiXxUgepsvxGhF+SIs0LejMjNc//lRhFVDe6scD3ghaHUxCxsoaJd3Ri8YxOIR9rv7/iVHcT3wHndfHMjXuvsBM/s64ENm9nF3/4sl1z8OMQqZwooCqpjMcX7p1fS6hyFGqOkydJdSpwMtnbdKanDhO5ip+1F3d0dRhJVhSlMX77415K3WMHlADe5qUAc5BO3asZwYfdDdrzrB3w8Aly78fsnw2rG4HnjV4gvufmD499Nm9mH6fNLNFqNqNf2RdZRcPocJX7o9laO7UTmmVVVccgfZDmvFHll9DrOyI/TG/aqgbogiLMjgKo6gCCtLm68RiNnS2yBYVxxLOUoZv4o7gSvM7HJ6EXo98LKj32Rm3wCcD/zhwmvnAwfdfdvMLgCeC/yCMphRiFHQxaB8sVMdlAz5kurOVAvBwOUenW2GaVUVIeUZZg7KMQOT8h0a0x29HA9HOTobaMuPowhL+vwKZsJSj8UxiFnQBW2kCmwozq7kjLr73MxeDdxB39rpHe7+STN7A3CXu982vPV64BZ/cpjnG4G3mj0+/+W+xSr8ZRiNGC2N/NwiupKN6UVUar6mnCiILuTM9c4M6vVaTVVoGk0M5niGlvtjFi7iAt1ZzRLZzSBm1WK24kVY4ue7OJEF1OFKjsGdVQkxXBfufjtw+1Gvvf6o33/uGMv9AfAtOccyGjGqRHfNdGdVvnOZSS5K570rqFC6PRXAXAzTt2LOp5nuZqlFK+osUsgTMOhuVjLXBanYYosmyduh+LwnGY7H0v1aVTELdeS9KoxFzKqUFsPrye70Ga2N0YhRhRwJ0eoF3x7/v+Up3Z4qx9z0+ixSYhEXXnwMoFXTN6aLKLWAqSVHzqhG4yZthxxdDSJVYDx5rwo1hNhrSDUIlsN8V3JGq6IKMWro03mWRhW0nXzbgiZDvqVEhrnpVSGnOqs5pjTN0RFAEXIpQ7P2LL1KCzfuz/FsVboIq4Z+rWOYfKC0mIXygrYGdxZC0J42MTf9ZlH6uaNxfQxq8VCWdJ7Ceas1pBqU7ggAefJWVZSbn5yqQIZ2aTnSHcTlk/WCVFpHBWK2BndVIUe/1xyFXAo1uLOQR9BuHLtTTV8VIUYzIbuzeu2PvIIcMwepqGJylkxy2d1ddndLdwRocDrxgt+qB2OGY0mdllXNGc3hRagCok81UHKgXRazqnZIaPnw/To0ckwcUEMhl0oNglZhIwug3EOMrhLNCcpRZSivYq0/H3QHJkfr/hzTqirfoq+mH0FHALWIShTUNHoBU452a2rOqPIdssyklUHMqrsyy/2/glQB6fPN5VSsHIVcCjkmcVBZdzFcil2am74qqhGj646e4aajnqjqdKLJXHcP1GtN0tbRGKh3b7Wa3lx3QdTtqM8ipWlRM5dnFJMLoBDFYIYbp9w7uIIirL5bSflUAYUcvlQNqQYqY3B3149wRleKFFpFu1hkqabPsA7dl1SrLTWHunl8LcszyxDiVkL9remFXHJXArRUgwaYyTmj+tGYhGOhcZMVcY58zzbDOhTGUoSlUkPeq0oNqQYhaNcLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 1152x360 with 2 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"err_pi = np.sum( np.abs(sol_pi-sol_gt)) / (STEPS*N)\n",
|
||||
"err_dp = np.sum( np.abs(sol_dp-sol_gt)) / (STEPS*N)\n",
|
||||
"print(\"MAE PINN: \\t{:7.5f} \\nMAE DP: \\t{:7.5f}\".format(err_pi,err_dp))\n",
|
||||
"\n",
|
||||
"print(\"\\nError GT to PINN (top) , DiffPhys (bottom):\")\n",
|
||||
"show_state(np.reshape( np.concatenate([sol_pi-sol_gt, divider, sol_dp-sol_gt],axis=0) ,[N*2+10,33,1]))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"That's a pretty clear result: the PINN error is ca. 5 times higher than the Differentiable Physics (DP) reconstruction.\n",
|
||||
"\n",
|
||||
"This difference also shows clearly in the jointly visualized image at the bottom: the magnitudes of the errors of the DP reconstruction are much closer to zero, as indicated by the purple color above.\n",
|
||||
"\n",
|
||||
"A simple direct reconstruction problem like this one is always a good initial test for a DP solver, e.g., before moving to more complex setups like coupling it with an NN. If the direct optimization does not converge, there's probably still something fundamentally wrong, and there's no point involving an NN. \n",
|
||||
"\n",
|
||||
"Now we have a first example to show similarities and differences of the two approaches. In the section, we'll present a more in depth discussion of these findings, before moving to more complex cases.\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"## Next Steps\n",
|
||||
"\n",
|
||||
"As with the PINN version, there's variety of things that can be improved and experimented with the code above:\n",
|
||||
"\n",
|
||||
"* You can try to adjust the training parameters to improve the reconstruction.\n",
|
||||
"* As for the PINN case, you can activate a different optimizer, and observe the changing (not necessarily improved) behavior.\n",
|
||||
"* Vary the number of steps, or the resolution of the simulation and reconstruction.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
@@ -277,7 +617,7 @@
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.7.6"
|
||||
"version": "3.8.5"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
|
||||
@@ -336,7 +336,8 @@
|
||||
"\n",
|
||||
"Now we have both versions, so let's compare both reconstructions in more detail.\n",
|
||||
"\n",
|
||||
"Let's first look at the solutions side by side. The code below generates an image with 3 versions, from top to bottom: the \"ground truth\" (GT) solution as given by the regular forward simulation, in the middle the PINN reconstruction, and at the bottom the differentiable physics version."
|
||||
"Let's first look at the solutions side by side. The code below generates an image with 3 versions, from top to bottom: the \"ground truth\" (GT) solution as given by the regular forward simulation, in the middle the PINN reconstruction, and at the bottom the differentiable physics version.\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -481,7 +482,7 @@
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.7.6"
|
||||
"version": "3.8.5"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
|
||||
+13
-13
@@ -25,7 +25,7 @@ TODO, visual overview of DP training
|
||||
|
||||
## Differentiable Operators
|
||||
|
||||
With this direction we build on existing numerical solvers. I.e.,
|
||||
With the DP direction we build on existing numerical solvers. I.e.,
|
||||
the approach is strongly relying on the algorithms developed in the larger field
|
||||
of computational methods for a vast range of physical effects in our world.
|
||||
To start with we need a continuous formulation as model for the physical effect that we'd like
|
||||
@@ -35,7 +35,7 @@ method for discretization of the equation.
|
||||
|
||||
Let's assume we have a continuous formulation $\mathcal P^*(\mathbf{x}, \nu)$ of the physical quantity of
|
||||
interest $\mathbf{u}(\mathbf{x}, t): \mathbb R^d \times \mathbb R^+ \rightarrow \mathbb R^d$,
|
||||
with a model parameter $\nu$ (e.g., a diffusion or viscosity constant).
|
||||
with model parameters $\nu$ (e.g., diffusion, viscosity, or conductivity constants).
|
||||
The component of $\mathbf{u}$ will be denoted by a numbered subscript, i.e.,
|
||||
$\mathbf{u} = (u_1,u_2,\dots,u_d)^T$.
|
||||
%and a corresponding discrete version that describes the evolution of this quantity over time: $\mathbf{u}_t = \mathcal P(\mathbf{x}, \mathbf{u}, t)$.
|
||||
@@ -54,9 +54,11 @@ $\partial \mathcal P_i / \partial \mathbf{u}$.
|
||||
|
||||
Note that we typically don't need derivatives
|
||||
for all parameters of $\mathcal P$, e.g. we omit $\nu$ in the following, assuming that this is a
|
||||
given model parameter, with which the NN should not interact. Naturally, it can vary,
|
||||
by $\nu$ will not be the output of a NN representation. If this is the case, we can omit
|
||||
providing $\partial \mathcal P_i / \partial \nu$ in our solver.
|
||||
given model parameter, with which the NN should not interact.
|
||||
Naturally, it can vary within the solution manifold that we're interested in,
|
||||
but $\nu$ will not be the output of a NN representation. If this is the case, we can omit
|
||||
providing $\partial \mathcal P_i / \partial \nu$ in our solver. However, the following learning process
|
||||
natuarlly transfers to including $\nu$ as a degree of freedom.
|
||||
|
||||
## Jacobians
|
||||
|
||||
@@ -93,7 +95,7 @@ this would cause huge memory overheads and unnecessarily slow down training.
|
||||
Instead, for backpropagation, we can provide faster operations that compute products
|
||||
with the Jacobian transpose because we always have a scalar loss function at the end of the chain.
|
||||
|
||||
[TODO check transpose of Jacobians in equations]
|
||||
**[TODO check transpose of Jacobians in equations]**
|
||||
|
||||
Given the formulation above, we need to resolve the derivatives
|
||||
of the chain of function compositions of the $\mathcal P_i$ at some current state $\mathbf{u}^n$ via the chain rule.
|
||||
@@ -121,7 +123,7 @@ as this [nice survey by Baydin et al.](https://arxiv.org/pdf/1502.05767.pdf).
|
||||
|
||||
## Learning via DP Operators
|
||||
|
||||
Thus, long story short, once the operators of our simulator support computations of the Jacobian-vector
|
||||
Thus, once the operators of our simulator support computations of the Jacobian-vector
|
||||
products, we can integrate them into DL pipelines just like you would include a regular fully-connected layer
|
||||
or a ReLU activation.
|
||||
|
||||
@@ -175,9 +177,6 @@ procedure for a _forward_ solve.
|
||||
Note that to simplify things, we assume that $\mathbf{u}$ is only a function in space,
|
||||
i.e. constant over time. We'll bring back the time evolution of $\mathbf{u}$ later on.
|
||||
%
|
||||
[TODO, write out simple finite diff approx?]
|
||||
[denote discrete d as $\mathbf{d}$ below?]
|
||||
%
|
||||
Let's denote this re-formulation as $\mathcal P$. It maps a state of $d(t)$ into a
|
||||
new state at an evoled time, i.e.:
|
||||
|
||||
@@ -186,7 +185,7 @@ $$
|
||||
$$
|
||||
|
||||
As a simple example of an optimization and learning task, let's consider the problem of
|
||||
finding an motion $\mathbf{u}$ such that starting with a given initial state $d^{~0}$ at $t^0$,
|
||||
finding a motion $\mathbf{u}$ such that starting with a given initial state $d^{~0}$ at $t^0$,
|
||||
the time evolved scalar density at time $t^e$ has a certain shape or configuration $d^{\text{target}}$.
|
||||
Informally, we'd like to find a motion that deforms $d^{~0}$ into a target state.
|
||||
The simplest way to express this goal is via an $L^2$ loss between the two states. So we want
|
||||
@@ -273,8 +272,9 @@ be preferable to actually constructing $A$.
|
||||
|
||||
## A (slightly) more complex example
|
||||
|
||||
[TODO]
|
||||
more complex, matrix inversion, eg Poisson solve
|
||||
**[TODO]**
|
||||
|
||||
a bit more complex, matrix inversion, eg Poisson solve
|
||||
dont backprop through all CG steps (available in phiflow though)
|
||||
rather, re-use linear solver to compute multiplication by inverse matrix
|
||||
|
||||
|
||||
@@ -82,6 +82,12 @@ See also... Test link: {doc}`supervised`
|
||||
---
|
||||
|
||||
|
||||
## TODOs , include
|
||||
|
||||
- general motivation: repeated solves in classical solvers -> potential for ML
|
||||
- PINNs: often need weighting of added loss terms for different parts
|
||||
|
||||
|
||||
## TODOs , Planned content
|
||||
|
||||
Loose collection of notes and TODOs:
|
||||
|
||||
+11
-12
@@ -23,14 +23,14 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"execution_count": 17,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Using phiflow version: 2.0.0\n"
|
||||
"Using phiflow version: 2.0.0rc0\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
@@ -156,7 +156,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[<matplotlib.lines.Line2D at 0x7faed7f5eeb0>]"
|
||||
"[<matplotlib.lines.Line2D at 0x7f7ef47d8160>]"
|
||||
]
|
||||
},
|
||||
"execution_count": 5,
|
||||
@@ -177,9 +177,8 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# we only need \"velocity.data\" from each phiflow state\n",
|
||||
"#vels = [v.values.numpy('x') for v in velocities]\n",
|
||||
"vels = [v.values.numpy('x,vector') for v in velocities] # vel vx\n",
|
||||
"# get \"velocity.values\" from each phiflow state with a channel dimensions, i.e. \"vector\"\n",
|
||||
"vels = [v.values.numpy('x,vector') for v in velocities] # gives a list of 2D arrays \n",
|
||||
"\n",
|
||||
"import pylab\n",
|
||||
"\n",
|
||||
@@ -203,14 +202,14 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"execution_count": 16,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Vels array shape: (128, 33, 1)\n"
|
||||
"resulting image size(128, 528)\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -230,18 +229,18 @@
|
||||
"def show_state(a):\n",
|
||||
" # we only have 33 time steps, blow up by a factor of 2^4 to make it easier to see\n",
|
||||
" # (could also be done with more evaluations of network)\n",
|
||||
" a=np.expand_dims(a, axis=2)\n",
|
||||
" for i in range(4):\n",
|
||||
" a = np.concatenate( [a,a] , axis=2)\n",
|
||||
"\n",
|
||||
" a = np.reshape( a, [a.shape[0],a.shape[1]*a.shape[2]] )\n",
|
||||
" #print(\"resulting image size\" +format(a.shape))\n",
|
||||
" print(\"resulting image size\" +format(a.shape))\n",
|
||||
"\n",
|
||||
" fig, axes = pylab.subplots(1, 1, figsize=(16, 5))\n",
|
||||
" im = axes.imshow(a, origin='upper', cmap='inferno')\n",
|
||||
" pylab.colorbar(im) \n",
|
||||
" \n",
|
||||
"#vels_img = np.asarray( np.stack(vels), dtype=np.float32 ).transpose() # no component channel\n",
|
||||
"vels_img = np.asarray( np.concatenate(vels, axis=-1), dtype=np.float32 ) # vel vx\n",
|
||||
"vels_img = np.reshape(vels_img, list(vels_img.shape)+[1] ) ; print(\"Vels array shape: \"+format(vels_img.shape))\n",
|
||||
"vels_img = np.asarray( np.concatenate(vels, axis=-1), dtype=np.float32 ) \n",
|
||||
"\n",
|
||||
"# save for comparison with reconstructions later on\n",
|
||||
"np.savez_compressed(\"./temp/burgers-groundtruth-solution.npz\", np.reshape(vels_img,[N,STEPS+1])) # remove batch & channel dimension\n",
|
||||
|
||||
Reference in new issue
Block a user