Skip to content

Commit 1e3859d

Browse files
jakubboberjakubbober
committed
Update notebook
1 parent 63bb632 commit 1e3859d

1 file changed

Lines changed: 65 additions & 78 deletions

File tree

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

Lines changed: 65 additions & 78 deletions
Original file line numberDiff line numberDiff line change
@@ -294,7 +294,7 @@
294294
},
295295
{
296296
"cell_type": "code",
297-
"execution_count": 6,
297+
"execution_count": 2,
298298
"metadata": {},
299299
"outputs": [],
300300
"source": [
@@ -307,11 +307,11 @@
307307
},
308308
{
309309
"cell_type": "code",
310-
"execution_count": 7,
310+
"execution_count": 3,
311311
"metadata": {},
312312
"outputs": [],
313313
"source": [
314-
"# %load -s solver_FECS, src-advec/advec1D.py\n",
314+
"# %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",
317317
" t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time\n",
@@ -338,10 +338,11 @@
338338
" u.data[1, :] = [I(xi) for xi in x]\n",
339339
" \n",
340340
" # Insert boundary condition\n",
341-
" bc = [Eq(u[t_s, 0], U0)]\n",
341+
" bc = [Eq(u[t_s+1, 0], U0)]\n",
342342
" \n",
343343
" op = Operator([eq] + bc)\n",
344-
" op.apply(time_m=1, dt=dt)\n",
344+
" op.apply(dt=dt, x_m=1, x_M=Nx-1)\n",
345+
" \n",
345346
" if user_action is not None:\n",
346347
" for n in range(0, Nt + 1):\n",
347348
" user_action(u.data[n], x, t, n)\n"
@@ -425,11 +426,10 @@
425426
},
426427
{
427428
"cell_type": "code",
428-
"execution_count": 8,
429+
"execution_count": 4,
429430
"metadata": {},
430431
"outputs": [],
431432
"source": [
432-
"# %load -s run_FECS, src-advec/advec1D.py\n",
433433
"def run_FECS(case):\n",
434434
" \"\"\"Special function for the FECS case.\"\"\"\n",
435435
" if case == 'gaussian':\n",
@@ -686,20 +686,17 @@
686686
"t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time\n",
687687
"...\n",
688688
"if scheme == 'FE':\n",
689-
" u = TimeFunction(name='u', grid=grid, space_order=2, save=2)\n",
689+
" u = TimeFunction(name='u', grid=grid, space_order=2, save=Nt+1)\n",
690690
" 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",
691+
" ...\n",
694692
"elif scheme == 'LF':\n",
693+
" u = ...\n",
694+
" \n",
695695
" # Use some scheme for the first step\n",
696-
" u1 = ...\n",
697-
" pde1 = ...\n",
696+
" pde0 = ...\n",
698697
" ...\n",
699698
" \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",
699+
" # Move to the LF scheme after the first step\n",
703700
" pde = u.dtc + v*u.dxc\n",
704701
" ...\n",
705702
"```"
@@ -1088,18 +1085,18 @@
10881085
"cell_type": "markdown",
10891086
"metadata": {},
10901087
"source": [
1091-
"mathcal{I}_t means that we in the first equation, involving $u_0^n$, insert $u_{N_x}^n$,\n",
1088+
"It means that we in the first equation, involving $u_0^n$, insert $u_{N_x}^n$,\n",
10921089
"and that we in the last equation, involving $u^{n+1}_{N_x}$ insert $u^{n+1}_0$.\n",
1093-
"Normally, we can do this in the simple way that `u_1[0]` is updated as\n",
1094-
"`u_1[Nx]` at the beginning of a new time level.\n",
1090+
"Normally, we can do this in the simple way that `u[t_s, 0]` is updated as\n",
1091+
"`u[t_s, Nx]` at the beginning of a new time level.\n",
10951092
"\n",
10961093
"In some schemes we may need $u^{n}_{N_x+1}$ and $u^{n}_{-1}$. Periodicity\n",
10971094
"then means that these values are equal to $u^n_1$ and $u^n_{N_x-1}$,\n",
10981095
"respectively. For the upwind scheme, it is sufficient to set\n",
1099-
"`u_1[0]=u_1[Nx]` at a new time level before computing `u[1]`. This ensures\n",
1100-
"that `u[1]` becomes right and at the next time level `u[0]` at the current\n",
1096+
"`u[t_s, 0]=u[t_s, Nx]` at a new time level before computing `u[t_s+1, 1]`. This ensures\n",
1097+
"that `u[t_s+1, 1]` becomes right and at the next time level `u[t_s+1, 0]` at the current\n",
11011098
"time level is correctly updated.\n",
1102-
"For the Leapfrog scheme we must update `u[0]` and `u[Nx]` using the scheme:"
1099+
"For the Leapfrog scheme we must update `u[t_s+1, 0]` and `u[t_s+1, Nx]` using the scheme:"
11031100
]
11041101
},
11051102
{
@@ -1108,10 +1105,8 @@
11081105
"source": [
11091106
"```python\n",
11101107
"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",
1108+
"pbc += [Eq(u[t_s+1, Nx], u[t_s+1, 0])]\n",
11121109
"...\n",
1113-
"\n",
1114-
"op = Operator(bc + (pbc if periodic_bc else []) + [eq])\n",
11151110
"```"
11161111
]
11171112
},
@@ -1174,15 +1169,15 @@
11741169
"metadata": {},
11751170
"source": [
11761171
"```python\n",
1177-
"dx*(0.5*u[0] + 0.5*u[Nx] + np.sum(u[1:-1]))\n",
1172+
"dx*(0.5*u.data[n][0] + 0.5*u.data[n][Nx] + np.sum(u.data[n][1:Nx]))\n",
11781173
"```"
11791174
]
11801175
},
11811176
{
11821177
"cell_type": "markdown",
11831178
"metadata": {},
11841179
"source": [
1185-
"if `u` is an array holding the solution.\n",
1180+
"if `u` is a `TimeFunction` with the `save` parameter set to $Nx+1$ and `n` indicates a current timestep.\n",
11861181
"\n",
11871182
"### The code\n",
11881183
"\n",
@@ -1192,13 +1187,14 @@
11921187
},
11931188
{
11941189
"cell_type": "code",
1195-
"execution_count": 12,
1190+
"execution_count": 7,
11961191
"metadata": {},
11971192
"outputs": [],
11981193
"source": [
1194+
"# %load -s solver src-advec/advec1D.py\n",
11991195
"def solver(I, U0, v, L, dt, C, T, user_action=None,\n",
12001196
" scheme='FE', periodic_bc=True):\n",
1201-
" Nt = int(round(T/float(dt)))\n",
1197+
" Nt = int(round(T/np.float64(dt)))\n",
12021198
" t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time\n",
12031199
" dx = v*dt/C\n",
12041200
" Nx = int(round(L/dx))\n",
@@ -1211,11 +1207,10 @@
12111207
" print('dt=%g, dx=%g, Nx=%d, C=%g' % (dt, dx, Nx, C))\n",
12121208
"\n",
12131209
" integral = np.zeros(Nt+1)\n",
1214-
"\n",
1210+
" \n",
12151211
" grid = Grid(shape=(Nx+1,), extent=(L,), dtype=np.float64)\n",
1212+
"\n",
12161213
" t_s=grid.time_dim\n",
1217-
" u = None\n",
1218-
" pde = None\n",
12191214
" \n",
12201215
" def u(to=1, so=1):\n",
12211216
" u = TimeFunction(name='u', grid=grid, time_order=to, space_order=so, save=Nt+1)\n",
@@ -1226,73 +1221,66 @@
12261221
" pde = u.dtr + v*u.dxc\n",
12271222
" \n",
12281223
" 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",
1224+
" pbc += [Eq(u[t_s+1, Nx], u[t_s+1, 0])]\n",
12301225
" \n",
12311226
" elif scheme == 'LF':\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",
1227+
" # Use UP scheme for the first timestep\n",
1228+
" u = u(to=2, so=2)\n",
1229+
" pde0 = u.dtr(fd_order=1) + v*u.dxl(fd_order=1)\n",
12541230
" \n",
1255-
" if user_action is not None:\n",
1256-
" user_action(u1.data[1], x, t, 1)\n",
1231+
" stencil0 = solve(pde0, u.forward)\n",
1232+
" eq0 = Eq(u.forward, stencil0).subs(t_s, 0)\n",
1233+
"\n",
1234+
" pbc0 = [Eq(u[t_s, 0], u[t_s, Nx]).subs(t_s, 0)]\n",
12571235
" \n",
1258-
" print('I:', integral[1])\n",
1259-
" \n",
12601236
" # Now continue with LF scheme\n",
1261-
" u = u(to=2, so=2)\n",
1262-
" u.data[0:2, :] = u1.data\n",
12631237
" pde = u.dtc + v*u.dxc\n",
12641238
" \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",
1239+
" pbc = [Eq(u[t_s+1, 0], u[t_s-1, 0] - C*(u[t_s, 1] - u[t_s, Nx-1]))]\n",
1240+
" pbc += [Eq(u[t_s+1, Nx], u[t_s+1, 0])]\n",
12671241
" \n",
12681242
" elif scheme == 'UP':\n",
12691243
" u = u()\n",
12701244
" pde = u.dtr + v*u.dxl\n",
12711245
" \n",
12721246
" pbc = [Eq(u[t_s, 0], u[t_s, Nx])]\n",
1273-
"\n",
1247+
" \n",
1248+
" elif scheme == 'LW':\n",
1249+
" u = u(so=2)\n",
1250+
" pde = u.dtr + v*u.dxc - 0.5*dt*v**2*u.dx2\n",
1251+
" \n",
1252+
" pbc = [Eq(u[t_s+1, 0], u[t_s, 0] - 0.5*C*(u[t_s, 1] - u[t_s, Nx-1]) + \\\n",
1253+
" 0.5*C**2*(u[t_s, 1] - 2*u[t_s, 0] + u[t_s, Nx-1]))]\n",
1254+
" pbc += [Eq(u[t_s+1, Nx], u[t_s+1, 0])]\n",
1255+
" \n",
12741256
" else:\n",
12751257
" raise ValueError('scheme=\"%s\" not implemented' % scheme)\n",
12761258
"\n",
12771259
" stencil = solve(pde, u.forward)\n",
12781260
" eq = Eq(u.forward, stencil)\n",
12791261
" \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",
1262+
" bc_init = [Eq(u[t_s+1, 0], U0).subs(t_s, 0)]\n",
1263+
" \n",
1264+
" # Set initial condition u(x,0) = I(x)\n",
1265+
" u.data[0, :] = [I(xi) for xi in x]\n",
12831266
" \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",
1267+
" # Compute the integral under the curve\n",
1268+
" integral[0] = dx*(0.5*u.data[0][0] + 0.5*u.data[0][Nx] + np.sum(u.data[0][1:Nx]))\n",
12861269
" \n",
1287-
" if user_action is not None:\n",
1288-
" user_action(u.data[0], x, t, 0)\n",
1270+
" if user_action is not None:\n",
1271+
" user_action(u.data[0], x, t, 0)\n",
12891272
"\n",
12901273
" bc = [Eq(u[t_s+1, 0], U0)]\n",
1274+
" \n",
1275+
" if scheme == 'LF':\n",
1276+
" op = Operator((pbc0 if periodic_bc else []) + [eq0] + (bc_init if not periodic_bc else []) \\\n",
1277+
" + (pbc if periodic_bc else []) + [eq] + (bc if not periodic_bc else []))\n",
1278+
" else: \n",
1279+
" op = Operator(bc_init + (pbc if periodic_bc else []) + [eq] + (bc if not periodic_bc else []))\n",
1280+
" \n",
1281+
" op.apply(dt=dt, x_m=1, x_M=Nx if scheme == 'UP' else Nx-1)\n",
12911282
"\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",
1283+
" for n in range(1, Nt+1):\n",
12961284
" # Compute the integral under the curve\n",
12971285
" integral[n] = dx*(0.5*u.data[n][0] + 0.5*u.data[n][Nx] + np.sum(u.data[n][1:Nx]))\n",
12981286
"\n",
@@ -1317,11 +1305,10 @@
13171305
},
13181306
{
13191307
"cell_type": "code",
1320-
"execution_count": 15,
1308+
"execution_count": 8,
13211309
"metadata": {},
13221310
"outputs": [],
13231311
"source": [
1324-
"# %load -s run, src-advec/advec1D.py\n",
13251312
"def run(scheme='UP', case='gaussian', C=1, dt=0.01):\n",
13261313
" \"\"\"General admin routine for explicit and implicit solvers.\"\"\"\n",
13271314
"\n",
@@ -2397,7 +2384,7 @@
23972384
},
23982385
{
23992386
"cell_type": "code",
2400-
"execution_count": 152,
2387+
"execution_count": 17,
24012388
"metadata": {},
24022389
"outputs": [],
24032390
"source": [

0 commit comments

Comments
 (0)