Skip to content

Commit f518900

Browse files
jakubboberjakubbober
committed
Update notebook
1 parent 5655842 commit f518900

1 file changed

Lines changed: 107 additions & 105 deletions

File tree

‎fdm-devito-notebooks/04_advec/advec.ipynb‎

Lines changed: 107 additions & 105 deletions
Original file line numberDiff line numberDiff line change
@@ -294,7 +294,7 @@
294294
},
295295
{
296296
"cell_type": "code",
297-
"execution_count": 36,
297+
"execution_count": 6,
298298
"metadata": {},
299299
"outputs": [],
300300
"source": [
@@ -307,24 +307,25 @@
307307
},
308308
{
309309
"cell_type": "code",
310-
"execution_count": 37,
310+
"execution_count": 7,
311311
"metadata": {},
312312
"outputs": [],
313313
"source": [
314314
"# %load -s solver_FECS, src-advec/advec1D.py\n",
315315
"def solver_FECS(I, U0, v, L, dt, C, T, user_action=None):\n",
316316
" Nt = int(round(T/float(dt)))\n",
317+
" t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time\n",
317318
" dx = v*dt/C\n",
318319
" Nx = int(round(L/dx))\n",
319-
"\n",
320-
" t = np.linspace(0, Nt*dt, Nt+1)\n",
321-
" x = np.linspace(0, L, Nx+1)\n",
320+
" x = np.linspace(0, L, Nx+1) # Mesh points in space\n",
321+
" \n",
322+
" # Make sure dx and dt are compatible with x and t\n",
322323
" dx = float(x[1] - x[0])\n",
323324
" dt = float(t[1] - t[0])\n",
324325
" C = v*dt/dx\n",
325326
"\n",
326327
" grid = Grid(shape=(Nx+1,), extent=(L,))\n",
327-
" t_s=grid.stepping_dim\n",
328+
" t_s=grid.time_dim\n",
328329
"\n",
329330
" u = TimeFunction(name='u', grid=grid, space_order=2, save=Nt+1)\n",
330331
"\n",
@@ -333,9 +334,12 @@
333334
" stencil = solve(pde, u.forward)\n",
334335
" eq = Eq(u.forward, stencil)\n",
335336
" \n",
336-
" bc = [Eq(u[t_s+1, 0], U0)]\n",
337-
" u.data[1, :] = [I(xs) for xs in x]\n",
338-
"\n",
337+
" # Set initial condition u(x,0) = I(x)\n",
338+
" u.data[1, :] = [I(xi) for xi in x]\n",
339+
" \n",
340+
" # Insert boundary condition\n",
341+
" bc = [Eq(u[t_s, 0], U0)]\n",
342+
" \n",
339343
" op = Operator([eq] + bc)\n",
340344
" op.apply(time_m=1, dt=dt)\n",
341345
" if user_action is not None:\n",
@@ -421,7 +425,7 @@
421425
},
422426
{
423427
"cell_type": "code",
424-
"execution_count": 38,
428+
"execution_count": 8,
425429
"metadata": {},
426430
"outputs": [],
427431
"source": [
@@ -674,45 +678,31 @@
674678
]
675679
},
676680
{
677-
"cell_type": "code",
678-
"execution_count": 141,
681+
"cell_type": "markdown",
679682
"metadata": {},
680-
"outputs": [
681-
{
682-
"ename": "NameError",
683-
"evalue": "name 'scheme' is not defined",
684-
"output_type": "error",
685-
"traceback": [
686-
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
687-
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
688-
"\u001b[0;32m<ipython-input-141-61a51f7851a3>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0;34m...\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mn\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mNt\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----> 9\u001b[0;31m \u001b[0;32mif\u001b[0m \u001b[0mscheme\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'FE'\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 10\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;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mNx\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 11\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mu_1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m0.5\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mC\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mu_1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mu_1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\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",
689-
"\u001b[0;31mNameError\u001b[0m: name 'scheme' is not defined"
690-
]
691-
}
692-
],
693683
"source": [
684+
"```python\n",
694685
"Nt = int(round(T/float(dt)))\n",
695686
"t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time\n",
696687
"...\n",
697-
"u = np.zeros(Nx+1)\n",
698-
"u_1 = np.zeros(Nx+1)\n",
699-
"u_2 = np.zeros(Nx+1)\n",
700-
"...\n",
701-
"for n in range(0, Nt):\n",
702-
" if scheme == 'FE':\n",
703-
" for i in range(1, Nx):\n",
704-
" u[i] = u_1[i] - 0.5*C*(u_1[i+1] - u_1[i-1])\n",
705-
" elif scheme == 'LF':\n",
706-
" if n == 0:\n",
707-
" # Use some scheme for the first step\n",
708-
" for i in range(1, Nx):\n",
709-
" ...\n",
710-
" else:\n",
711-
" for i in range(1, Nx+1):\n",
712-
" u[i] = u_2[i] - C*(u_1[i] - u_1[i-1])\n",
713-
"\n",
714-
" # Switch variables before next step\n",
715-
" u_2, u_1, u = u_1, u, u_2"
688+
"if scheme == 'FE':\n",
689+
" u = TimeFunction(name='u', grid=grid, space_order=2, save=2)\n",
690+
" pde = u.dtr + v*u.dxc\n",
691+
" \n",
692+
" for i in range(1, Nx):\n",
693+
" u[i] = u_1[i] - 0.5*C*(u_1[i+1] - u_1[i-1])\n",
694+
"elif scheme == 'LF':\n",
695+
" # Use some scheme for the first step\n",
696+
" u1 = ...\n",
697+
" pde1 = ...\n",
698+
" ...\n",
699+
" \n",
700+
" # Move to the LF scheme after the next step\n",
701+
" u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2, save=Nt+1)\n",
702+
" u.data[0:2, :] = u1.data\n",
703+
" pde = u.dtc + v*u.dxc\n",
704+
" ...\n",
705+
"```"
716706
]
717707
},
718708
{
@@ -1113,30 +1103,16 @@
11131103
]
11141104
},
11151105
{
1116-
"cell_type": "code",
1117-
"execution_count": 147,
1106+
"cell_type": "markdown",
11181107
"metadata": {},
1119-
"outputs": [
1120-
{
1121-
"ename": "NameError",
1122-
"evalue": "name 'periodic_bc' is not defined",
1123-
"output_type": "error",
1124-
"traceback": [
1125-
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
1126-
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
1127-
"\u001b[0;32m<ipython-input-147-6c4e6ba84a70>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mif\u001b[0m \u001b[0mperiodic_bc\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 2\u001b[0m \u001b[0mi\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mu_2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mC\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mu_1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mu_1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mNx\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\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 4\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;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mNx\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 5\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mu_2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mC\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mu_1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mu_1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\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",
1128-
"\u001b[0;31mNameError\u001b[0m: name 'periodic_bc' is not defined"
1129-
]
1130-
}
1131-
],
11321108
"source": [
1133-
"if periodic_bc:\n",
1134-
" i = 0\n",
1135-
" u[i] = u_2[i] - C*(u_1[i+1] - u_1[Nx-1])\n",
1136-
"for i in range(1, Nx):\n",
1137-
" u[i] = u_2[i] - C*(u_1[i+1] - u_1[i-1])\n",
1138-
"if periodic_bc:\n",
1139-
" u[Nx] = u[0]"
1109+
"```python\n",
1110+
"pbc = [Eq(u[t_s+1, 0], u[t_s-1, 0] - C*(u[t_s, 1] - u[t_s, Nx-1]))]\n",
1111+
"pbc += [Eq(u[t_s, Nx], u[t_s, 0])]\n",
1112+
"...\n",
1113+
"\n",
1114+
"op = Operator(bc + (pbc if periodic_bc else []) + [eq])\n",
1115+
"```"
11401116
]
11411117
},
11421118
{
@@ -1194,23 +1170,12 @@
11941170
]
11951171
},
11961172
{
1197-
"cell_type": "code",
1198-
"execution_count": 143,
1173+
"cell_type": "markdown",
11991174
"metadata": {},
1200-
"outputs": [
1201-
{
1202-
"data": {
1203-
"text/plain": [
1204-
"0.0"
1205-
]
1206-
},
1207-
"execution_count": 143,
1208-
"metadata": {},
1209-
"output_type": "execute_result"
1210-
}
1211-
],
12121175
"source": [
1213-
"dx*(0.5*u[0] + 0.5*u[Nx] + np.sum(u[1:-1]))"
1176+
"```python\n",
1177+
"dx*(0.5*u[0] + 0.5*u[Nx] + np.sum(u[1:-1]))\n",
1178+
"```"
12141179
]
12151180
},
12161181
{
@@ -1227,14 +1192,12 @@
12271192
},
12281193
{
12291194
"cell_type": "code",
1230-
"execution_count": null,
1195+
"execution_count": 12,
12311196
"metadata": {},
12321197
"outputs": [],
12331198
"source": [
1234-
"# %load -s solver, src-advec/advec1D.py\n",
12351199
"def solver(I, U0, v, L, dt, C, T, user_action=None,\n",
1236-
" scheme='FE', periodic_bc=False):\n",
1237-
" print('USING DEVITO')\n",
1200+
" scheme='FE', periodic_bc=True):\n",
12381201
" Nt = int(round(T/float(dt)))\n",
12391202
" t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time\n",
12401203
" dx = v*dt/C\n",
@@ -1249,11 +1212,8 @@
12491212
"\n",
12501213
" integral = np.zeros(Nt+1)\n",
12511214
"\n",
1252-
"# if user_action is not None:\n",
1253-
"# user_action(u_n, x, t, 0)\n",
1254-
"\n",
1255-
" grid = Grid(shape=(Nx+1,), extent=(L,))\n",
1256-
" t_s=grid.stepping_dim\n",
1215+
" grid = Grid(shape=(Nx+1,), extent=(L,), dtype=np.float64)\n",
1216+
" t_s=grid.time_dim\n",
12571217
" u = None\n",
12581218
" pde = None\n",
12591219
" \n",
@@ -1265,30 +1225,74 @@
12651225
" u = u(so=2)\n",
12661226
" pde = u.dtr + v*u.dxc\n",
12671227
" \n",
1228+
" pbc = [Eq(u[t_s+1, 0], u[t_s, 0] - 0.5*C*(u[t_s, 1] - u[t_s, Nx]))]\n",
1229+
" pbc += [Eq(u[t_s, Nx], u[t_s, 0])]\n",
1230+
" \n",
12681231
" elif scheme == 'LF':\n",
1269-
" u = u(to=2, so=2)\n",
1232+
" # Use UP scheme for first timestep\n",
1233+
" u1 = TimeFunction(name='u1', grid=grid, save=2)\n",
1234+
" pde1 = u1.dtr + v*u1.dxl\n",
1235+
" \n",
1236+
" stencil1 = solve(pde1, u1.forward)\n",
1237+
" eq1 = Eq(u1.forward, stencil1)\n",
1238+
" \n",
1239+
" # Set initial condition u(x,0) = I(x)\n",
1240+
" u1.data[0, :] = [I(xi) for xi in x]\n",
1241+
" \n",
1242+
" bc1 = [Eq(u1[t_s+1, 0], U0)]\n",
1243+
" pbc1 = [Eq(u1[t_s, 0], u1[t_s, Nx])]\n",
1244+
" \n",
1245+
" integral[0] = dx*(0.5*u1.data[0][0] + 0.5*u1.data[0][Nx] + np.sum(u1.data[0][1:Nx]))\n",
1246+
" \n",
1247+
" if user_action is not None:\n",
1248+
" user_action(u1.data[0], x, t, 0)\n",
1249+
" \n",
1250+
" op1 = Operator(bc1 + (pbc1 if periodic_bc else []) + [eq1])\n",
1251+
" op1.apply(dt=dt)\n",
1252+
" \n",
1253+
" integral[1] = dx*(0.5*u1.data[1][0] + 0.5*u1.data[1][Nx] + np.sum(u1.data[1][1:Nx]))\n",
1254+
" \n",
1255+
" if user_action is not None:\n",
1256+
" user_action(u1.data[1], x, t, 1)\n",
1257+
" \n",
1258+
" print('I:', integral[1])\n",
1259+
" \n",
1260+
" # Now continue with LF scheme\n",
1261+
" u = u(to=2, so=2)\n",
1262+
" u.data[0:2, :] = u1.data\n",
12701263
" pde = u.dtc + v*u.dxc\n",
1264+
" \n",
1265+
" pbc = [Eq(u[t_s+1, 0], u[t_s-1, 0] - C*(u[t_s, 1] - u[t_s, Nx - 1]))]\n",
1266+
" pbc += [Eq(u[t_s, Nx], u[t_s, 0])]\n",
12711267
" \n",
12721268
" elif scheme == 'UP':\n",
12731269
" u = u()\n",
1274-
" pde = u.dtr - v*u.dxl\n",
1275-
" \n",
1276-
" elif scheme == 'LW':\n",
1277-
" pass\n",
1278-
" \n",
1270+
" pde = u.dtr + v*u.dxl\n",
1271+
" \n",
1272+
" pbc = [Eq(u[t_s, 0], u[t_s, Nx])]\n",
1273+
"\n",
12791274
" else:\n",
12801275
" raise ValueError('scheme=\"%s\" not implemented' % scheme)\n",
12811276
"\n",
12821277
" stencil = solve(pde, u.forward)\n",
12831278
" eq = Eq(u.forward, stencil)\n",
12841279
" \n",
1285-
" bc = [Eq(u[t_s+1, 0], U0)] # non-periodic boundary condition\n",
1286-
" u.data[0:2, :] = [I(xs) for xs in x]\n",
1280+
" if scheme != 'LF':\n",
1281+
" # Set initial condition u(x,0) = I(x)\n",
1282+
" u.data[0, :] = [I(xi) for xi in x]\n",
1283+
" \n",
1284+
" # Compute the integral under the curve\n",
1285+
" integral[0] = dx*(0.5*u.data[0][0] + 0.5*u.data[0][Nx] + np.sum(u.data[0][1:Nx]))\n",
1286+
" \n",
1287+
" if user_action is not None:\n",
1288+
" user_action(u.data[0], x, t, 0)\n",
12871289
"\n",
1288-
" op = Operator([eq] + bc)\n",
1289-
" op.apply(time_m=1, dt=float(dt))\n",
1290-
" print(u.data)\n",
1291-
" for n in range(0, Nt+1):\n",
1290+
" bc = [Eq(u[t_s+1, 0], U0)]\n",
1291+
"\n",
1292+
" op = Operator(bc + (pbc if periodic_bc else []) + [eq])\n",
1293+
" op.apply(time_m=1 if scheme == 'LF' else 0, time_M=Nt-1, dt=float(dt))\n",
1294+
"\n",
1295+
" for n in range(2 if scheme == 'LF' else 1, Nt+1):\n",
12921296
" # Compute the integral under the curve\n",
12931297
" integral[n] = dx*(0.5*u.data[n][0] + 0.5*u.data[n][Nx] + np.sum(u.data[n][1:Nx]))\n",
12941298
"\n",
@@ -1313,7 +1317,7 @@
13131317
},
13141318
{
13151319
"cell_type": "code",
1316-
"execution_count": 151,
1320+
"execution_count": 15,
13171321
"metadata": {},
13181322
"outputs": [],
13191323
"source": [
@@ -1341,7 +1345,6 @@
13411345
" lines = plt.plot(x, u)\n",
13421346
" plt.axis([x[0], x[-1], -0.5, 1.5])\n",
13431347
" plt.xlabel('x'); plt.ylabel('u')\n",
1344-
" plt.axes().set_aspect(0.15)\n",
13451348
" plt.savefig('tmp_%04d.png' % n)\n",
13461349
" plt.savefig('tmp_%04d.pdf' % n)\n",
13471350
" else:\n",
@@ -1401,7 +1404,6 @@
14011404
" plt.figure(2)\n",
14021405
" plt.axis([0, L, -0.5, 1.1])\n",
14031406
" plt.xlabel('$x$'); plt.ylabel('$u$')\n",
1404-
" plt.axes().set_aspect(0.5) # no effect\n",
14051407
" plt.savefig('tmp1.png'); plt.savefig('tmp1.pdf')\n",
14061408
" plt.show()\n",
14071409
" # Make videos from figure(1) animation files\n",

0 commit comments

Comments
 (0)