update physical loss chapter
This commit is contained in:
NT committed
1 parent
48cf5ea80b
commit
2374ecac67
3 files changed
+903
-1361
No files matched your search
@@ -1,600 +0,0 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Burgers Optimization with a Physics-Informed NN\n",
|
||||
"\n",
|
||||
"To illustrate how the physics-informed losses work, let's consider a \n",
|
||||
"reconstruction example with a simple yet non-linear equation in 1D:\n",
|
||||
"Burgers equation $\\frac{\\partial u}{\\partial{t}} + u \\nabla u = \\nu \\nabla \\cdot \\nabla u$\n",
|
||||
"for which we have a series of _observations_ at time $t=0.5$ which should be fulfilled in the solution. In addition, let's impose Dirichlet boundary conditions $u=0$\n",
|
||||
"at the sides of our computational domain, and define the solution in\n",
|
||||
"the time interval $t \\in [0,1]$.\n",
|
||||
"\n",
|
||||
"Note that similar to the previous forward simulation example, \n",
|
||||
"we will still be sampling the solution with 128 points ($n=128), but now we have a discretization via the NN. So we could also sample points inbetween without having to explicitly choose a basis function for interpolation. The discretization via the NN now internally determines how to use its degrees of freedom to construct the basis functions. So we have no direct control over the reconstruction.\n",
|
||||
"\n",
|
||||
"## Preliminaries\n",
|
||||
"\n",
|
||||
"Let's just load TF and phiflow for now, and initialize the random sampling."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"#from phi.tf.flow import *\n",
|
||||
"import numpy as np\n",
|
||||
"import numpy.random as rnd\n",
|
||||
"import tensorflow as tf\n",
|
||||
"\n",
|
||||
"#rnd = TF_BACKEND # for phiflow: sample different points in the domain each iteration\n",
|
||||
"#rnd = math.choose_backend(1) # use same random points for all iterations"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We're importing phiflow here, but we won't use it to compute a solution to the PDE. Below we'll instead use the derivatives of an NN (as explained in the previous section) to set up a loss formulation to train an NN.\n",
|
||||
"\n",
|
||||
"Next, we set up a simple NN with 8 fully connected layers and `tanh` activations with 20 units each. \n",
|
||||
"\n",
|
||||
"We'll also define the `boundary_tx` function which gives an array of constraints for the solution (all for $=0.5$ in this example), and the `open_boundary` function which contains constraints for $x= \\pm1$ being 0."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def network(x, t):\n",
|
||||
" \"\"\" Dense neural network with 8 hidden layers and 3021 parameters in total.\n",
|
||||
" Parameters will only be allocated once (auto reuse).\n",
|
||||
" \"\"\"\n",
|
||||
" y = tf.stack([x, t], axis=-1)\n",
|
||||
" for i in range(8):\n",
|
||||
" #y = tf.keras.layers.dense(y, 20, activation=tf.math.tanh, name='layer%d' % i, reuse=tf.AUTO_REUSE)\n",
|
||||
" y = tf.compat.v1.layers.dense(y, 20, activation=tf.math.tanh, name='layer%d' % i, reuse=True)\n",
|
||||
" return tf.compat.v1.layers.dense(y, 1, activation=None, name='layer_out', reuse=True)\n",
|
||||
"\n",
|
||||
"def boundary_tx(N):\n",
|
||||
" x = np.linspace(-1,1,128)\n",
|
||||
" # precomputed solution from forward simulation:\n",
|
||||
" u = 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",
|
||||
" t = np.asarray( np.ones(x.shape)) * 0.5\n",
|
||||
" perm = np.random.permutation(128) \n",
|
||||
" return (x[perm])[0:N], (t[perm])[0:N], (u[perm])[0:N]\n",
|
||||
"\n",
|
||||
"def _ALT_t0(N): # alternative, impose original initial state at t=0\n",
|
||||
" x = rnd.uniform(-1, 1, [N])\n",
|
||||
" t = np.zeros(x.shape)\n",
|
||||
" u = -np.sin(np.pi * x)\n",
|
||||
" return x, t, u\n",
|
||||
"\n",
|
||||
"def open_boundary(N):\n",
|
||||
" t = rnd.uniform(0, 1, [N])\n",
|
||||
" x = np.concatenate([np.zeros([N//2]) + 1, np.zeros([N//2]) - 1], axis=0)\n",
|
||||
" u = np.zeros([N])\n",
|
||||
" return x, t, u\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Most importantly, we can now also construct the residual loss function `f` that we'd like to minimize in order to guide the NN to retrieve a solution for our model equation. As can be seen in the equation at the top, we need derivatives w.r.t. $t$, $x$ and a second derivative for $x$. The first three lines of `f` below do just that.\n",
|
||||
"\n",
|
||||
"Afterwards, we simply combine the derivates according to Burgers equation:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def f(u, x, t):\n",
|
||||
" \"\"\" Physics-based loss function with Burgers equation \"\"\"\n",
|
||||
" #with tf.GradientTape() as tape: \n",
|
||||
" u_t = tf.compat.v1.gradients(u, t)\n",
|
||||
" u_x = tf.compat.v1.gradients(u, x)\n",
|
||||
" u_xx = tf.compat.v1.gradients(u_x, x)\n",
|
||||
" return u_t + u*u_x - (0.01 / np.pi) * u_xx\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Next, let's set up the sampling points in the inner domain, such that we can compare the solution with the previous forward simulation in phiflow. The last `math.expand_dims()` call simply adds another `batch` dimension, so that the resulting tensor is compatible with the following examples.\n",
|
||||
"\n",
|
||||
"A bit more interesting: grid_u will afterwards contain a full graph to evaluate our NN at $128 \\times 33$ positions, and will retrieve the results in a $[1,128,33,1]$ array once we run it through `session.run`."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"WARNING:tensorflow:From <ipython-input-2-5de5bf78dcb3>:8: dense (from tensorflow.python.keras.legacy_tf_layers.core) is deprecated and will be removed in a future version.\n",
|
||||
"Instructions for updating:\n",
|
||||
"Use keras.layers.Dense instead.\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/keras/legacy_tf_layers/core.py:187: Layer.apply (from tensorflow.python.keras.engine.base_layer) is deprecated and will be removed in a future version.\n",
|
||||
"Instructions for updating:\n",
|
||||
"Please use `layer.__call__` method instead.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# generate array with positions: \n",
|
||||
"# -1 to 1 spatial with 128 cells\n",
|
||||
"# 0 to 1 in time with 1+32 timesteps\n",
|
||||
"N=128\n",
|
||||
"grids_xt = np.meshgrid(np.linspace(-1, 1, N), np.linspace(0, 1, 33), indexing='ij')\n",
|
||||
"grid_x, grid_t = [tf.convert_to_tensor(t, tf.float32) for t in grids_xt]\n",
|
||||
"\n",
|
||||
"# create 4D tensor with batch and channel dimensions in addition to space and time\n",
|
||||
"# in this case gives shape=(1, n, 33, 1)\n",
|
||||
"grid_u = tf.reshape(network(grid_x, grid_t), (1,N,33,1) )"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Let's give this a try: we can initialize a TF session, evaluate `grid_u` and show it in an image, just like the phiflow solution we computed previously. \n",
|
||||
"\n",
|
||||
"(Note, as before the x axis does not show actual simulation time, but is showing 32 steps \"blown\" up by a factor of 16 to make the changes over time easier to see in the image.)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {
|
||||
"scrolled": true
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Size of grid_u: (1, 128, 33, 1)\n",
|
||||
"WARNING:tensorflow:From /Users/thuerey/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/util/tf_should_use.py:247: initialize_all_variables (from tensorflow.python.ops.variables) is deprecated and will be removed after 2017-03-02.\n",
|
||||
"Instructions for updating:\n",
|
||||
"Use `tf.global_variables_initializer` instead.\n",
|
||||
"Randomly initialized network state:\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAB0CAYAAAB+I3LkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbfklEQVR4nO2dbahlV3nH/89MfGsUbaJOpyZtlAZFStUqvqAfokWJIuaLBG2haQnMlwoWCnWk0NJ+Sr/UWlrEgQYVWl9KGwwi1TQq/dJqfDe+REeJmCE6KNEKOvecvdfTD3vfeE3mzvk9c5991777rD8Mc++5a9be68w9v/U///Wsdczd1dTU1NS0LB2rfQNNTU1NTflqcG9qampaoBrcm5qamhaoBvempqamBarBvampqWmBanBvampqWqAmgbuZ3Whm95nZWTM7PcU1mpqampr2l2XXuZvZcUnflPRqSQ9IukfSm939a6kXampqamraV1M49xdLOuvu33H3laQPSLppgus0NTU1Ne2jKeD+DEnf2/P9A+NjTU1NTU2HpCtqXdjMTkk6JUlXXvn4Fz77Ob+x+R/BCMlo0uQuCTYO9ZnYTg6vHbgu7a/WmKdoN/cxl+zr0u4C/aU/N7nd8ftLvq4kubHLJuPmyw/96Ifu/rSL/WwKuJ+TdO2e768ZH/slufsZSWck6YUvfLZ/+n/etbnnvkM3YKVIBbxaSpGctHOpY9dG7bzIup71t17B6/ZwzM6ex+IyMhbanyStYLv1Orc/2q7vpY793jhp5z78v2zsz6UV+33wNWy3grNFV9DE4l1BxPHO5R1Dk69YuwJ/HWg774y97HtjgC+m0gO4F9ZfcdafJP36Hbd/d7+fTQH3eyRdb2bP1AD1N0n6/Uv/E2eQ1QjuTaJw9zK8sLL62+1zY3/O+yP3J433CNr2XWDM2e2yxwyfR/r/3LF79K4MEwG5PzIJSPLQmGm7zX16B9v1ketubuZFcuhvHM7NKgFow7nZu83JtfsAbtKu71l/BbTbpHS4u3tnZm+R9DFJxyXd7u5f3fgPwYtqEkeePFlgR45BRx0YvS4bs9F3Al74tQkQJd4fbsfg7l0PJ0jYLgB39XASgK5YnbMJo4ftOmf+q7B7dOrnyhiPAHiWTnIQj5Se9Yeddm8qJJYppgKv25cZwl2S3P2jkj5K25tXdOQoekieLCQOOuqM4ZitCzh35E7phAsjCmmSyYKBLjIW1s4htIWhTScBHqOwiY/C2JEj98Ic+eCec2MU746xdhDufQ/7c4PO3dSDSWCTqi2o/rIgFGm2iyHSpU4W+J2FByKKbCDSCajvYV4ciCio28XQDoCOqBQMbdLOI7EMbgfzbBijDHBn7Vj+DKFdpLJmAMMxSm8oHqF5dl9gf7SdM+c+9DlT5x6WizsmogniG7a4CEEnpbtTy87cOwh3GstQ5x4CYuA5xJk7cbG7WcEGFb64iJ17cozinWDmzl4q3g2QJaKQjU0Cm9tSaNPsu7ihOMhdKG5ZnnNHkA04sMyMPNuRhxYXAw6fPIc0V470h3Pl7IkvUhHC4I7jG5IrZzt351UwBRZZcYc/ZtUbL2wq8A026k98EujA4udufwjGMPt2F8vcIbSHdkty7uQFTcsRoTsNLdDS/DkzypD4Am0PF0BDi4akPzhm55EHjlGyF1R7OgnAe6TOvbh8DeMW+BLwnpYjQkc+Zt8b2xUGY1qJIgUcPm53DMG4FLhQGnDa2+ncs2MUBCa4SkTr3AM13yjmkWLrDDRLz4xlaEZeXB6pNyfKjlG6gmFMnbuvk7N03E7IkdPKEVqJMpQ4wsoROAmQRUiJRSiS1MF2fTnGqm8CcQubLKTeF+XcwQsfu124UErdKXXkkbLAKdplxzKZtdylpMcoDjf+aN3D6IFOFi4BZ+x4ElC1DT1lRSOKetDuocMfYEzbMWiTdn1h7wQi0F6WcwdgCrldXAqZWRESqHOP5M9zj1GII59gt6ZgjELL/XDc4tBBFw9k37Bdx96E1CwfRLnyBDXf3JGzCS0C7aLk+GZJzh1BMTuioG43eyt+cWkFrdU6sMU+NZYprBolEMtU22K/ogulgUoUWO6Hsu/CF0Bp5UhPK0zW2wlt2l8NaJdlxTIw045EFImxjEVimfRcObu6hY2F59T0PBbH0MZlgRDuvoZwD5QPUmhn13xHoI3aQefe98fSoc0rR2C1TH8cX5dAm8Y8xU0O+otAe1mxDIEdjjxgtQzdrdl3ubXcEoc2PmyLjSW0xR7HMombjgrP0rN3YZaVeM03du65lR6RjTqZuytphUms3K8etGl/taDdged6k+YBdxpTTAHt1Gze8080hLGMr+iYAxt66C7MzBilOI5bfB2ANuwvc6MOLQuUmHPfXaysAe3sssBILfcawzgX2l0A7miBVnwS6Jbl3DdDzGhOHdqAk5kr00XIwrN0ei4KzuYZ3IeTD5NhTOFOa75XEMZwo05Z8cXK7N2VPYxROriwuY3Q7nC7+UO7X4xzd2fgnvsuzMg5K9C5+xquH0SgnVzzPQW0icoObBeANq4wqVAWuLsdPhPapRxj59JBaBc3DCa6+5OCk08WPBqhkwCaIMWhvZxYhu5QXa/zoV2p5hsvLl5g71Z81WNo45pvUqPdFb64iGu5udNm/eWXBS4F2rzmm0G7ZuXIkLlvFl0Ajbhs8ptdnN3fbp8H1UzgDrPqFdxRGoF26uJiANrUuU+xUScb2sjhj9Am7RYG7Y3354bL/bJ3V2aX+4WcO3bkuf0NscxmlWRou1vgfLglwZ3ALuLcs6FNP1aNQvsCXVCF0KYVJiteYVIT2kTdKlDLnXhaYCnsXG6Jw5i3Y5UjkTwbLdA6eycQgztrh8cSuC569yOeuaNoS4afQzgHXFLzgTt17qkVIX06tMlGnVCFyYU60K5Vy+1F6sEBVVI9aGeX+02xCIkdPoQ2bTe40/nDnahzpWfu9NTy5Tj34ixbXlJZYDa0V6XKBpzSC58WSKHd7zB3So947bv8PJvmymvotGmlxxpDm42lh2eRD+6UAMxSK0JcfHExHe44vmGTgBQ5tn8pcHdH4J49tCMbcEiFicN2ksoFulgptliZDe1ioXI/olrQLtBp00qPEliEjLhYAu25V4RMkVNntyts/R+fQTNMFgfXPOBeXJ7o3HFZYE/r0iG0i/gRrwDaodMCd1Cz0CIk/XAGAu3hhAkKd+h2e75oSM8cQXBXLtwlDu3I4mJmvjtFTh1x0ETwCPsA3Fl/BWbpNJYpYgUAmzQPuEPnjqG9pjXftHRRqZ+UUxPa/Tq/ciQT2h6sHEHRQ2jREE4WidvSh3b5zp3GMrg/0I5GFJExRzLyzP6ogy7QkTt0+BKfWC6l+cAdVI9gaK/yz+VGMKaVI5IcVI5EDpSaonKEwd0gtAOVHhG4Jx4AlV0RQhcXIzk1zZ8jOXV2pUctB83HzPqjDnp4t8LabeGCquQ7JHOH0D4C5X79DvxF7Ezkk3J6eHRr5MwREj30xdAntZcI3JMrQtbwXO7I4mJNR05EKz0ikQfO5lGP84c7ddqRjDzyruGgOhDczex+ST+V1Evq3P1FZnaVpA9Kuk7S/ZJudveHLtWPu6PqEd/heXZm5Yh3sNzPpQIddL+CAKPtOpg/w6NbI5tg6Lb07IoQPAmUY+hFmu7clV/pwTP3XOfO3Wm9/JlPAqzdcI+5JY6Zz80mZTj3V7r7D/d8f1rS3e5+m5mdHr9/2yV7KK7ycwB3mFP7ylmuTD8cuGef6O5u+HQ/Wha4XtOywOPw+PXjeFt6dkUIhntgcZE6cuzcwXUjG3UISGLlg9CRT5A/14tbWLsp3HM2jPk9woaX0BSxzE2Sbhi/fq+kT2kj3CUHpXyRM0cyT/cbsu/kipAVrwghYOKVIzyioHDnMM5vl13pQV742YuGscVF1CwUUWQuBk6xaJgdZeDMPbCgSpTt8DfpoHB3SR83M5f0bnc/I+mEuz84/vz7kk5s7MSFTvjraeUI/NDfocIEtOv5R4yRD/OlZ4lEKkfWPXPu68La9X4MuUkXPHNEE8A9OXqgk4Arv4KDxy2sv/xyv4iLzd6tmdsuG8ZTwD2B7QeG+yvc/ZyZPV3SXWb2jb0/dHcfwf8omdkpSack6donXqn+52TRENZow8XF0tPt5vllfKtALTdxVrgsMBBRMNfJ21E3mQ13utNwivx5CqedWZlRE8bZkUctuO+2Jdel6wf03c+ldCC4u/u58e/zZnaHpBdL+oGZnXT3B83spKTz+/zbM5LOSNILrr7a1xeI084t46MVIbSML1IRsgqUBZIIYIqcOhXuyl80jCwusgXV/MXAbNeZ7dzJ9g2Jj7kmjHEVDGuWDuNIVU1J2MV02XA3syslHXP3n45fv0bS30i6U9Itkm4b//7wpr7cDZ0nssbb14/z405huV9mRUhkcXEnUDnCPzpss3paow23h9NFw93Ig8Uom9tIu1vsYbtU514PxhHQZUO72iIkfG6oK54Cxpm/X5t0EOd+QtIdZrbbz7+6+3+a2T2SPmRmt0r6rqSbN3XkblqBBcZaFSHZi4shuCcvLnb4zBFeJ01BQmOUSP6MNjEFnF+NKCNSiZIN4+yIIgRtPDnD/lizEIjpRIAd/lHI3N39O5Ked5HHfyTp9yJ9lSLtrDbfCq0IoYuLvbOzRCKLixRMO/S0wEBGTqOHiIPeJAqmSEQRgztrx/qjuzADGWtiO6micw/EE7k5NStr3r02um4AxNkTxpGAe6bcTTtrBvfsxUW605A48gjckcNXJH+OlAWSdswVF48sbKJm2KnhrfgBIGZDu1ZFSGR3JepP9aKHbCBGQBzJ54kc9jiHUsgUuUw73eZbiSwu0tK3XLjnV3rsgIXh3f5QZIU36kDnLujcYX+SBA/WDAExE8ZTRBT5tdd1ooej4HYjIM6OmPhzfXC6zwLuxQ1Vj0QWFwlIpqgIIU47VhbIc2Vcyw2uG8mLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 432x288 with 1 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {
|
||||
"needs_background": "light"
|
||||
},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# test image output\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"print(\"Size of grid_u: \"+format(grid_u.shape))\n",
|
||||
"\n",
|
||||
"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",
|
||||
" for i in range(4):\n",
|
||||
" a = np.concatenate( [a,a] , axis=3)\n",
|
||||
" a = np.reshape( a, [a.shape[1],a.shape[2]*a.shape[3]] )\n",
|
||||
" #print(a.shape)\n",
|
||||
" plt.imshow(a, origin='upper', cmap='magma')\n",
|
||||
" \n",
|
||||
"print(\"Randomly initialized network state:\")\n",
|
||||
"#show_state(session.run(grid_u))\n",
|
||||
"show_state(grid_u)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This visualization already shows a smooth transition over space and time. So far, this is purely the random initialization of the NN that we're sampling here. So it has nothing to do with a solution of our PDE-based model up to now.\n",
|
||||
"\n",
|
||||
"The next steps will actually evaluate the constraints in terms of data (from the `boundary` functions), and the model constraints from `f` to retrieve an actual solution to the PDE."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## The loss functions & training\n",
|
||||
"\n",
|
||||
"As objective for the learning process we can now combine the _direct_ constraints, i.e., the solution at $t=0.5$ and the Dirchlet $u=0$ boundary conditions with the loss from the PDE residuals. For both boundary constraints we'll use 100 points below, and then sample the solution in the inner region with an additional 1000 points.\n",
|
||||
"\n",
|
||||
"The direct constraints are evaluated via `network(x, t)[:, 0] - u`, where `x` and `t` are the \"position\" where we'd like to sample out solution, and `u` provides the corresponding ground truth value.\n",
|
||||
"\n",
|
||||
"For the physical loss points, we have no ground truth solutions, but we'll only evaluate the PDE residual via the NN derivatives, to see whether the solution satisfies PDE model. If not, this directly gives us an error to be reduced via a update step in the optimization. The corresponding expression is of the form `f(network(x, t)[:, 0], x, t)` below. Note that for both data and physics terms the `network()[:, 0]` simply discard the last size-1 dimension of the $(n,1)$ tensor returned by the network."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"ERROR:tensorflow:==================================\n",
|
||||
"Object was never used (type <class 'tensorflow.python.framework.ops.Operation'>):\n",
|
||||
"<tf.Operation 'init' type=NoOp>\n",
|
||||
"If you want to mark it as used call its \"mark_used()\" method.\n",
|
||||
"It was originally created here:\n",
|
||||
" File \"/Users/thuerey/miniconda3/envs/tf/lib/python3.8/site-packages/IPython/core/interactiveshell.py\", line 3155, in run_cell_async\n",
|
||||
" has_raised = await self.run_ast_nodes(code_ast.body, cell_name, File \"/Users/thuerey/miniconda3/envs/tf/lib/python3.8/site-packages/IPython/core/interactiveshell.py\", line 3347, in run_ast_nodes\n",
|
||||
" if (await self.run_code(code, result, async_=asy)): File \"/Users/thuerey/miniconda3/envs/tf/lib/python3.8/site-packages/IPython/core/interactiveshell.py\", line 3427, in run_code\n",
|
||||
" exec(code_obj, self.user_global_ns, self.user_ns) File \"<ipython-input-7-351d69eca4ac>\", line 9, in <module>\n",
|
||||
" tf.compat.v1.initialize_all_variables() File \"/Users/thuerey/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/util/tf_should_use.py\", line 247, in wrapped\n",
|
||||
" return _add_should_use_warning(fn(*args, **kwargs),\n",
|
||||
"==================================\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"ename": "ValueError",
|
||||
"evalue": "Variable layer0/kernel does not exist, or was not created with tf.get_variable(). Did you mean to set reuse=tf.AUTO_REUSE in VarScope?",
|
||||
"output_type": "error",
|
||||
"traceback": [
|
||||
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
|
||||
"\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)",
|
||||
"\u001b[0;32m<ipython-input-7-351d69eca4ac>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 11\u001b[0m \u001b[0;31m#with app.model_scope():\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m \u001b[0mloss_u\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0.5\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mnetwork\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx_bc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt_bc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mu_bc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m \u001b[0;31m# normalizes by first dimension, N_bc\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 13\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 14\u001b[0m \u001b[0;31m# Physics loss inside of domain\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m<ipython-input-2-5de5bf78dcb3>\u001b[0m in \u001b[0;36mnetwork\u001b[0;34m(x, t)\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m8\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0;31m#y = tf.keras.layers.dense(y, 20, activation=tf.math.tanh, name='layer%d' % i, reuse=tf.AUTO_REUSE)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 8\u001b[0;31m \u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompat\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mv1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlayers\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdense\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m20\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mactivation\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmath\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtanh\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'layer%d'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mreuse\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 9\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mtf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompat\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mv1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlayers\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdense\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mactivation\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'layer_out'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mreuse\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/util/deprecation.py\u001b[0m in \u001b[0;36mnew_func\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m 322\u001b[0m \u001b[0;34m'in a future version'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mdate\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m'after %s'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mdate\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 323\u001b[0m instructions)\n\u001b[0;32m--> 324\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 325\u001b[0m return tf_decorator.make_decorator(\n\u001b[1;32m 326\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnew_func\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'deprecated'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/keras/legacy_tf_layers/core.py\u001b[0m in \u001b[0;36mdense\u001b[0;34m(inputs, units, activation, use_bias, kernel_initializer, bias_initializer, kernel_regularizer, bias_regularizer, activity_regularizer, kernel_constraint, bias_constraint, trainable, name, reuse)\u001b[0m\n\u001b[1;32m 185\u001b[0m \u001b[0m_scope\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 186\u001b[0m _reuse=reuse)\n\u001b[0;32m--> 187\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mlayer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mapply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 188\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 189\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/util/deprecation.py\u001b[0m in \u001b[0;36mnew_func\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m 322\u001b[0m \u001b[0;34m'in a future version'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mdate\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m'after %s'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mdate\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 323\u001b[0m instructions)\n\u001b[0;32m--> 324\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 325\u001b[0m return tf_decorator.make_decorator(\n\u001b[1;32m 326\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnew_func\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'deprecated'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/keras/engine/base_layer_v1.py\u001b[0m in \u001b[0;36mapply\u001b[0;34m(self, inputs, *args, **kwargs)\u001b[0m\n\u001b[1;32m 1699\u001b[0m \u001b[0mOutput\u001b[0m \u001b[0mtensor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1700\u001b[0m \"\"\"\n\u001b[0;32m-> 1701\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1702\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1703\u001b[0m @deprecation.deprecated(\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/keras/legacy_tf_layers/base.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, inputs, *args, **kwargs)\u001b[0m\n\u001b[1;32m 545\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 546\u001b[0m \u001b[0;31m# Actually call layer\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 547\u001b[0;31m \u001b[0moutputs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msuper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mLayer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 548\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 549\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mcontext\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexecuting_eagerly\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/keras/engine/base_layer_v1.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 755\u001b[0m \u001b[0;31m# Build layer if applicable (if the `build` method has been\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 756\u001b[0m \u001b[0;31m# overridden).\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 757\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_maybe_build\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 758\u001b[0m \u001b[0mcast_inputs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_maybe_cast_inputs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 759\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/keras/engine/base_layer_v1.py\u001b[0m in \u001b[0;36m_maybe_build\u001b[0;34m(self, inputs)\u001b[0m\n\u001b[1;32m 2096\u001b[0m \u001b[0;31m# operations.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2097\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mtf_utils\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmaybe_init_scope\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2098\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuild\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minput_shapes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2099\u001b[0m \u001b[0;31m# We must set self.built since user defined build functions are not\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2100\u001b[0m \u001b[0;31m# constrained to set self.built.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/keras/layers/core.py\u001b[0m in \u001b[0;36mbuild\u001b[0;34m(self, input_shape)\u001b[0m\n\u001b[1;32m 1169\u001b[0m 'should be defined. Found `None`.')\n\u001b[1;32m 1170\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minput_spec\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mInputSpec\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmin_ndim\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxes\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m{\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mlast_dim\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1171\u001b[0;31m self.kernel = self.add_weight(\n\u001b[0m\u001b[1;32m 1172\u001b[0m \u001b[0;34m'kernel'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1173\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mlast_dim\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munits\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/keras/legacy_tf_layers/base.py\u001b[0m in \u001b[0;36madd_weight\u001b[0;34m(self, name, shape, dtype, initializer, regularizer, trainable, constraint, use_resource, synchronization, aggregation, partitioner, **kwargs)\u001b[0m\n\u001b[1;32m 446\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0minitializer\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 447\u001b[0m \u001b[0minitializer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscope\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minitializer\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 448\u001b[0;31m variable = super(Layer, self).add_weight(\n\u001b[0m\u001b[1;32m 449\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 450\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/keras/engine/base_layer_v1.py\u001b[0m in \u001b[0;36madd_weight\u001b[0;34m(self, name, shape, dtype, initializer, regularizer, trainable, constraint, partitioner, use_resource, synchronization, aggregation, **kwargs)\u001b[0m\n\u001b[1;32m 429\u001b[0m \u001b[0mcaching_device\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 430\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 431\u001b[0;31m variable = self._add_variable_with_custom_getter(\n\u001b[0m\u001b[1;32m 432\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 433\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/training/tracking/base.py\u001b[0m in \u001b[0;36m_add_variable_with_custom_getter\u001b[0;34m(self, name, shape, dtype, initializer, getter, overwrite, **kwargs_for_getter)\u001b[0m\n\u001b[1;32m 743\u001b[0m \u001b[0minitializer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcheckpoint_initializer\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 744\u001b[0m \u001b[0mshape\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 745\u001b[0;31m new_variable = getter(\n\u001b[0m\u001b[1;32m 746\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 747\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/ops/variable_scope.py\u001b[0m in \u001b[0;36mget_variable\u001b[0;34m(name, shape, dtype, initializer, regularizer, trainable, collections, caching_device, partitioner, validate_shape, use_resource, custom_getter, constraint, synchronization, aggregation)\u001b[0m\n\u001b[1;32m 1554\u001b[0m \u001b[0msynchronization\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mVariableSynchronization\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mAUTO\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1555\u001b[0m aggregation=VariableAggregation.NONE):\n\u001b[0;32m-> 1556\u001b[0;31m return get_variable_scope().get_variable(\n\u001b[0m\u001b[1;32m 1557\u001b[0m \u001b[0m_get_default_variable_store\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1558\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/ops/variable_scope.py\u001b[0m in \u001b[0;36mget_variable\u001b[0;34m(self, var_store, name, shape, dtype, initializer, regularizer, reuse, trainable, collections, caching_device, partitioner, validate_shape, use_resource, custom_getter, constraint, synchronization, aggregation)\u001b[0m\n\u001b[1;32m 1297\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mdtype\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1298\u001b[0m \u001b[0mdtype\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_dtype\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1299\u001b[0;31m return var_store.get_variable(\n\u001b[0m\u001b[1;32m 1300\u001b[0m \u001b[0mfull_name\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1301\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/ops/variable_scope.py\u001b[0m in \u001b[0;36mget_variable\u001b[0;34m(self, name, shape, dtype, initializer, regularizer, reuse, trainable, collections, caching_device, partitioner, validate_shape, use_resource, custom_getter, constraint, synchronization, aggregation)\u001b[0m\n\u001b[1;32m 552\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mcustom_getter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0mcustom_getter_kwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 553\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 554\u001b[0;31m return _true_getter(\n\u001b[0m\u001b[1;32m 555\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 556\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/ops/variable_scope.py\u001b[0m in \u001b[0;36m_true_getter\u001b[0;34m(name, shape, dtype, initializer, regularizer, reuse, trainable, collections, caching_device, partitioner, validate_shape, use_resource, constraint, synchronization, aggregation)\u001b[0m\n\u001b[1;32m 505\u001b[0m \"name was already created with partitioning?\" % name)\n\u001b[1;32m 506\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 507\u001b[0;31m return self._get_single_variable(\n\u001b[0m\u001b[1;32m 508\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 509\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0;32m~/miniconda3/envs/tf/lib/python3.8/site-packages/tensorflow/python/ops/variable_scope.py\u001b[0m in \u001b[0;36m_get_single_variable\u001b[0;34m(self, name, shape, dtype, initializer, regularizer, partition_info, reuse, trainable, collections, caching_device, validate_shape, use_resource, constraint, synchronization, aggregation)\u001b[0m\n\u001b[1;32m 885\u001b[0m \u001b[0;31m# The code below handles only the case of creating a new variable.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 886\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mreuse\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 887\u001b[0;31m raise ValueError(\"Variable %s does not exist, or was not created with \"\n\u001b[0m\u001b[1;32m 888\u001b[0m \u001b[0;34m\"tf.get_variable(). Did you mean to set \"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 889\u001b[0m \"reuse=tf.AUTO_REUSE in VarScope?\" % name)\n",
|
||||
"\u001b[0;31mValueError\u001b[0m: Variable layer0/kernel does not exist, or was not created with tf.get_variable(). Did you mean to set reuse=tf.AUTO_REUSE in VarScope?"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# Boundary loss\n",
|
||||
"num_sample_points_bnd = 100\n",
|
||||
"x_bc, t_bc, u_bc = [np.concatenate([v_t0, v_x], axis=0) for v_t0, v_x in zip(boundary_tx(num_sample_points_bnd), open_boundary(num_sample_points_bnd))]\n",
|
||||
"x_bc, t_bc, u_bc = np.asarray(x_bc,dtype=np.float32), np.asarray(t_bc,dtype=np.float32) ,np.asarray(u_bc,dtype=np.float32)\n",
|
||||
"\n",
|
||||
"tf.compat.v1.disable_eager_execution()\n",
|
||||
"\n",
|
||||
"session = tf.compat.v1.Session(None)\n",
|
||||
"tf.compat.v1.initialize_all_variables()\n",
|
||||
"\n",
|
||||
"#with app.model_scope():\n",
|
||||
"loss_u = 0.5 * (network(x_bc, t_bc)[:, 0] - u_bc)**2 # normalizes by first dimension, N_bc\n",
|
||||
"\n",
|
||||
"# Physics loss inside of domain\n",
|
||||
"num_sample_points_inner = 1000\n",
|
||||
"x_ph, t_ph = tf.convert_to_tensor(rnd.uniform(-1, 1, [num_sample_points_inner])), tf.convert_to_tensor(rnd.uniform(0, 1, [num_sample_points_inner]))\n",
|
||||
"loss_ph = 0.5 * (f(network(x_ph, t_ph)[:, 0], x_ph, t_ph))**2 # normalizes by first dimension, N_ph\n",
|
||||
"\n",
|
||||
"# Combine\n",
|
||||
"ph_factor = 1.\n",
|
||||
"loss = loss_u + ph_factor * loss_ph # allows us to control the relative influence of loss_ph \n",
|
||||
"\n",
|
||||
"optim = tf.train.GradientDescentOptimizer(learning_rate=0.05).minimize(loss)\n",
|
||||
"#optim = tf.train.AdamOptimizer(learning_rate=0.001).minimize(loss) # alternative, but not much benefit here\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The code above just initializes the evaluation of the loss, we still didn't do any optimization steps, but we're finally in a good position to get started with this.\n",
|
||||
"\n",
|
||||
"Note: despite the simple equation, the convergence is typicaly very slow. It needs a _lot_ of iterations. To keep the runtime in a reasonable range, we only do 10k iterations by default below (`iters`). You can increase this value to get better results."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"session.initialize_variables()\n",
|
||||
"\n",
|
||||
"import time\n",
|
||||
"start = time.time()\n",
|
||||
"\n",
|
||||
"iters = 10000\n",
|
||||
"for optim_step in range(iters+1):\n",
|
||||
" _, loss_value = session.run([optim, loss])\n",
|
||||
" if optim_step<3 or optim_step%1000==0: \n",
|
||||
" print('Step %d, loss: %f' % (optim_step,loss_value))\n",
|
||||
" #show_state(grid_u)\n",
|
||||
" \n",
|
||||
"end = time.time()\n",
|
||||
"print(\"Runtime {:.2f}s\".format(end-start))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This training can take a significant amount of time, around 1 minute on a typical notebook, but at least the error goes down significantly (roughly from around 0.2 to ca. 0.03), and the network seems to successfully converge to a solution.\n",
|
||||
"\n",
|
||||
"Let's show the reconstruction of the network, by evaluating the network at the centers of a regular grid, so that we can show the solution as an image. Note that this is actually fairly expensive, we have to run through the whole network with a few thousand weights for all of the $128 \\times 32$ points in the grid.\n",
|
||||
"\n",
|
||||
"It looks pretty good on first sight, though. There's been a very noticeable change compared to the random initialization shown above:\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"show_state(session.run(grid_u)) "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"---\n",
|
||||
"\n",
|
||||
"## Evaluation\n",
|
||||
"\n",
|
||||
"Let's compare solution in a bit more detail. Here are the actual sample points used for constraining the solution (at time step 16, $t=1/2$) shown in gray, versus the reconstructed solution in blue:\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"u = session.run(grid_u)\n",
|
||||
"\n",
|
||||
"# solution is imposed at t=1/2 , which is 16 in the array\n",
|
||||
"bc_tx = 16 \n",
|
||||
"uT = u[0,:,bc_tx,0]\n",
|
||||
"\n",
|
||||
"fig = plt.figure().gca()\n",
|
||||
"fig.plot(np.linspace(-1,1,len(uT)), uT, lw=2, color='blue')\n",
|
||||
"fig.scatter(x_bc[0:100], u_bc[0:100], color='gray')\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Not too bad at the sides of the domain (the Dirichlet boundary conditions $u=0$ are fulfilled), but the shock in the center (at $x=0$) is not well represented.\n",
|
||||
"\n",
|
||||
"Let's check how well the initial state at $t=0$ was reconstructed. That's the most interesting, and toughest part of the problem (the rest basically follows from the model equation and boundary conditions, given the first state).\n",
|
||||
"\n",
|
||||
"It turns out, the accuracy of the initial state is actually not that great: the blue curve from the PINN is quite far away from the constraits (shown in gray)... The solution will get better with larger number of iterations, but it requires a surprisingly large number of them for a fairly simple case. \n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# ground truth solution at t0\n",
|
||||
"t0gt = np.asarray( [ [-np.sin(np.pi * x) * 1.] for x in np.linspace(-1,1,n)] )\n",
|
||||
"\n",
|
||||
"velP0 = u[0,:,0,0]\n",
|
||||
"\n",
|
||||
"fig = plt.figure().gca()\n",
|
||||
"fig.plot(np.linspace(-1,1,len(velP0)), velP0, lw=2, color='blue')\n",
|
||||
"fig.plot(np.linspace(-1,1,len(t0gt)), t0gt, lw=2, color='gray') # optionally show GT, compare to blue\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Especially the maximum / minimum at $x=\\pm 1/2$ are far off, and the boudaries at $x=\\pm 1$ are not fulfilled: the solution is not at zero.\n",
|
||||
"\n",
|
||||
"We have the forward simulator for this simulation, so we can use the $t=0$ solution of the network to \n",
|
||||
"evaluate how well the temporal evoluation was reconstructed by the PINN. This measures how well the temporal evolution of the model equation was captured via the soft constraints of the PINN loss.\n",
|
||||
"\n",
|
||||
"The graph below shows the initial state in blue, and two evolved states at $t=8/32$ and $t=15/32$. Note that this is all from the simulated version, we'll show the PINN version next. \n",
|
||||
"\n",
|
||||
"(Note: The code segments below also have some optional code to show the states at `[steps//4]`. It's commented out by default, you can uncomment or add additional ones to visualize more of the time evolution if you like.)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# re-simulate with phiflow from solution at t=0\n",
|
||||
"dt = 1./32.\n",
|
||||
"steps = 32-bc_tx # depends on where BCs were imposed\n",
|
||||
"initial = u[...,bc_tx:(bc_tx+1),0] # np.reshape(u0, [1,len(u0),1]) \n",
|
||||
"print(initial.shape)\n",
|
||||
"\n",
|
||||
"domain = Domain([n], boundaries=PERIODIC, box=box[-1:1])\n",
|
||||
"state = [BurgersVelocity(domain, velocity=initial, viscosity=0.01/np.pi)]\n",
|
||||
"physics = Burgers()\n",
|
||||
"\n",
|
||||
"for i in range(steps):\n",
|
||||
" state.append( physics.step(state[-1],dt=dt) )\n",
|
||||
"\n",
|
||||
"# we only need \"velocity.data\" from each phiflow state\n",
|
||||
"vel_resim = [x.velocity.data for x in state]\n",
|
||||
"\n",
|
||||
"fig = plt.figure().gca()\n",
|
||||
"pltx = np.linspace(-1,1,len(vels[0].flatten()))\n",
|
||||
"fig.plot(pltx, vel_resim[ 0].flatten(), lw=2, color='blue')\n",
|
||||
"#fig.plot(pltx, vel_resim[steps//4].flatten(), lw=2, color='green')\n",
|
||||
"fig.plot(pltx, vel_resim[steps//2].flatten(), lw=2, color='cyan')\n",
|
||||
"fig.plot(pltx, vel_resim[steps-1].flatten(), lw=2, color='purple')\n",
|
||||
"#fig.plot(pltx, t0gt, lw=2, color='gray') # optionally show GT, compare to blue\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"And here is the PINN output from `u` at the same time steps:\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"velP = [u[0,:,x,0] for x in range(33)]\n",
|
||||
"print(velP[0].shape)\n",
|
||||
"\n",
|
||||
"fig = plt.figure().gca()\n",
|
||||
"fig.plot(pltx, velP[bc_tx+ 0].flatten(), lw=2, color='blue')\n",
|
||||
"#fig.plot(pltx, velP[bc_tx+steps//4].flatten(), lw=2, color='green')\n",
|
||||
"fig.plot(pltx, velP[bc_tx+steps//2].flatten(), lw=2, color='cyan')\n",
|
||||
"fig.plot(pltx, velP[bc_tx+steps-1].flatten(), lw=2, color='purple')\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Not surprisingly, the errors grow over time. Especially the steepening of the solution near the shock at $x=0$ is not \"captured\" well. It's a bit difficult to see in these two graphs, though, let's quantify the error and show the actual difference:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"error = np.sum( np.abs( np.asarray(vel_resim[0:16]).flatten() - np.asarray(velP[bc_tx:bc_tx+steps]).flatten() )) / (steps*n)\n",
|
||||
"print(\"Mean absolute error for re-simulation across {} steps: {:7.5f}\".format(steps,error))\n",
|
||||
"\n",
|
||||
"fig = plt.figure().gca()\n",
|
||||
"fig.plot(pltx, (vel_resim[0 ].flatten()-velP[bc_tx ].flatten()), lw=2, color='blue')\n",
|
||||
"#fig.plot(pltx, (vel_resim[steps//4].flatten()-velP[bc_tx+steps//4].flatten()), lw=2, color='green')\n",
|
||||
"fig.plot(pltx, (vel_resim[steps//2].flatten()-velP[bc_tx+steps//2].flatten()), lw=2, color='cyan')\n",
|
||||
"fig.plot(pltx, (vel_resim[steps-1 ].flatten()-velP[bc_tx+steps-1 ].flatten()), lw=2, color='purple')\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This should show a mean absolute error of ca. $1.5 \\cdot 10^{-2}$ between ground truth re-simulation and the PINN evolution, which is significant for the value range of the simulation.\n",
|
||||
"\n",
|
||||
"And for comparison with the forward simulation and following cases, here are also all steps over time with a color map."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"scrolled": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# show re-simulated solution again as full image over time\n",
|
||||
"sn = np.concatenate(vel_resim, axis=-1)\n",
|
||||
"sn = np.reshape(sn, list(sn.shape)+[1] ) # print(sn.shape)\n",
|
||||
"show_state(sn)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Next, we'll store the full solution over the course of the $t=0 \\dots 1$ time interval, so that we can compare it later on to the full solution from a regular forward solve and compare it to the differential physics solution.\n",
|
||||
"\n",
|
||||
"Thus, stay tuned for the full evaluation and the comparison. It will follow in the file `diffphys-code-tf.ipynb`, after we've discussed the details of how to run the differential physics optimization."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"vels = session.run(grid_u) # special for showing NN results, run through TF \n",
|
||||
"vels = np.reshape( vels, [vels.shape[1],vels.shape[2]] )\n",
|
||||
"\n",
|
||||
"# save for comparison with other methods\n",
|
||||
"np.savez_compressed(\"./temp/burgers-pinn-solution.npz\",vels) ; print(\"Vels array shape: \"+format(vels.shape))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"---\n",
|
||||
"\n",
|
||||
"## Next steps\n",
|
||||
"\n",
|
||||
"This setup is just a starting point for PINNs and physical soft-constraints, of course. The parameters of the setup were chosen to compare to run relatively quickly. As we'll show in the next sections, the behavior of such an inverse solve can be improved substantially by a tighter integration of solver and learning. \n",
|
||||
"\n",
|
||||
"The solution of the PINN setup above can also directly be improved, however. E.g., try to:\n",
|
||||
"\n",
|
||||
"* Adjust parameters of the training to further decrease the error without making the solution diverge\n",
|
||||
"* Adapt the NN architecture for further improvements (keep track of the weight count, though!)\n",
|
||||
"* Activate a different optimizer, and observe the changing behavior (this typically requires adjusting the learning rate). Note that the more complex optimizers don't necessarily do better in this relatively simple example.\n",
|
||||
"* Or modify the setup to make the test case more interesting: e.g., move the boundary conditions further back in time, to let the reconstruction go \"further backward\" in time."
|
||||
]
|
||||
},
|
||||
{
|
||||
"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
|
||||
}
|
||||
+897
-755
File diff suppressed because it is too large.
Load diff
@@ -9,12 +9,12 @@ starting point.
|
||||
On the positive side, we can leverage DL frameworks with backpropagation to compute
|
||||
the derivatives of the model. At the same time, this puts us at the mercy of the learned
|
||||
representation regarding the reliability of these derivatives. Also, each derivative
|
||||
requires backpropagation through the full network, which can be very slow. Especially so
|
||||
requires backpropagation through the full network, which can be very expensive. Especially so
|
||||
for higher-order derivatives.
|
||||
|
||||
And while the setup is relatively simple, it is generally difficult to control. The NN
|
||||
has flexibility to refine the solution by itself, but at the same time, tricks are necessary
|
||||
when it doesn't pick the right regions of the solution.
|
||||
when it doesn't focus on the right regions of the solution.
|
||||
|
||||
## Is it "Machine Learning"?
|
||||
|
||||
@@ -36,7 +36,7 @@ about how well our model will generalize to "real-world" cases that we will enco
|
||||
we deploy it into an application.
|
||||
|
||||
In contrast, for the PINN training as described here, we reconstruct a single solution in a known
|
||||
and given space-time time. As such, any samples from this domain follow the same distribution
|
||||
and given space-time region. As such, any samples from this domain follow the same distribution
|
||||
and hence don't really represent test or OOD sampes. As the NN directly encodes the solution,
|
||||
there is also little hope that it will yield different solutions, or perform well outside
|
||||
of the training distribution. If we're interested in a different solution, we most likely
|
||||
@@ -47,7 +47,7 @@ have to start training the NN from scratch.
|
||||
## Summary
|
||||
|
||||
Thus, the physical soft constraints allow us to encode solutions to
|
||||
PDEs with the tools of ANNs.
|
||||
PDEs with the tools of NNs.
|
||||
An inherent drawback of this approach is that it yields single solutions,
|
||||
and that it does not combine with traditional numerical techniques well.
|
||||
E.g., learned representation is not suitable to be refined with
|
||||
@@ -60,12 +60,12 @@ goals of the next sections.
|
||||
|
||||
✅ Pro:
|
||||
- Uses physical model.
|
||||
- Derivatives can be conveniently compute via backpropagation.
|
||||
- Derivatives can be conveniently computed via backpropagation.
|
||||
|
||||
❌ Con:
|
||||
- Quite slow ...
|
||||
- Physical constraints are enforced only as soft constraints.
|
||||
- Largely incompatible _classical_ numerical methods.
|
||||
- Largely incompatible with _classical_ numerical methods.
|
||||
- Accuracy of derivatives relies on learned representation.
|
||||
|
||||
Next, let's look at how we can leverage numerical methods to improve the DL accuracy and efficiency
|
||||
|
||||
Reference in new issue
Block a user