diff --git a/.github/workflows/verification.yml b/.github/workflows/verification.yml index d7e8b019..920b483f 100644 --- a/.github/workflows/verification.yml +++ b/.github/workflows/verification.yml @@ -75,6 +75,10 @@ jobs: run: | cd fdm-devito-notebooks/02_wave/src-wave/wave1D $RUN_CMD python -m pytest -W ignore::DeprecationWarning -s -v --cov . --cov-config=.coveragerc --cov-report=xml:waves_coverage.xml $SKIP wave1D_u0.py::test_constant + - name: Waves (2.6 to 2.9) + run: | + cd fdm-devito-notebooks/02_wave/src-wave/wave1D + $RUN_CMD python -m pytest -W ignore::DeprecationWarning -s -v --cov . --cov-config=.coveragerc --cov-report=xml:waves_coverage.xml $SKIP wave1D_n0.py - name: Upload coverage to Codecov uses: codecov/codecov-action@v1.0.6 with: diff --git a/fdm-devito-notebooks/02_wave/.f41e97ebc401af3f04047239501937ae936b7ae5_archive.npz b/fdm-devito-notebooks/02_wave/.f41e97ebc401af3f04047239501937ae936b7ae5_archive.npz new file mode 100644 index 00000000..a9a7d8a1 Binary files /dev/null and b/fdm-devito-notebooks/02_wave/.f41e97ebc401af3f04047239501937ae936b7ae5_archive.npz differ diff --git a/fdm-devito-notebooks/02_wave/src-wave/wave1D/.f41e97ebc401af3f04047239501937ae936b7ae5_archive.npz b/fdm-devito-notebooks/02_wave/src-wave/wave1D/.f41e97ebc401af3f04047239501937ae936b7ae5_archive.npz new file mode 100644 index 00000000..d36795c7 Binary files /dev/null and b/fdm-devito-notebooks/02_wave/src-wave/wave1D/.f41e97ebc401af3f04047239501937ae936b7ae5_archive.npz differ diff --git a/fdm-devito-notebooks/02_wave/src-wave/wave1D/wave1D_dn.py b/fdm-devito-notebooks/02_wave/src-wave/wave1D/wave1D_dn.py index 850dfe09..0a219d7e 100644 --- a/fdm-devito-notebooks/02_wave/src-wave/wave1D/wave1D_dn.py +++ b/fdm-devito-notebooks/02_wave/src-wave/wave1D/wave1D_dn.py @@ -36,7 +36,11 @@ solution on the screen (as an animation). """ import numpy as np -import scitools.std as plt +# import scitools.std as plt +import time +from devito import Constant, Grid, TimeFunction, SparseTimeFunction, Function, Eq, solve, Operator, Buffer + +# TODO: Remove scalar vs vectorized version since Devito doesn't require these def solver(I, V, f, c, U_0, U_L, L, dt, C, T, user_action=None, version='scalar'): @@ -49,6 +53,106 @@ def solver(I, V, f, c, U_0, U_L, L, dt, C, T, dx = dt*c/float(C) Nx = int(round(L/dx)) x = np.linspace(0, L, Nx+1) # Mesh points in space + C2 = C**2 # Help variable in the scheme + + # Make sure dx and dt are compatible with x and t + dx = x[1] - x[0] + dt = t[1] - t[0] + + # Wrap user-given f, I, V, U_0, U_L if None or 0 + if f is None or f == 0: + f = (lambda x, t: 0) if version == 'scalar' else \ + lambda x, t: np.zeros(x.shape) + if I is None or I == 0: + I = (lambda x: 0) if version == 'scalar' else \ + lambda x: np.zeros(x.shape) + if V is None or V == 0: + V = (lambda x: 0) if version == 'scalar' else \ + lambda x: np.zeros(x.shape) + if U_0 is not None: + if isinstance(U_0, (float,int)) and U_0 == 0: + U_0 = lambda t: 0 + # else: U_0(t) is a function + if U_L is not None: + if isinstance(U_L, (float,int)) and U_L == 0: + U_L = lambda t: 0 + # else: U_L(t) is a function + + t0 = time.perf_counter() # Measure CPU time + + # Set up grid + grid = Grid(shape=(Nx+1), extent=(L)) + t_s = grid.stepping_dim + + # Create and initialise u + u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2, save=Nt+1) + u.data[0, :] = I(x[:]) + + if user_action is not None: + user_action(u.data[0], x, t, 0) + + x_dim = grid.dimensions[0] + t_dim = grid.time_dim + + # The wave equation we are trying to solve + pde = (1/c**2)*u.dt2-u.dx2 + + # Source term and injection into equation + dt_symbolic = grid.time_dim.spacing + src = SparseTimeFunction(name='f', grid=grid, npoint=Nx+1, nt=Nt+1) + + for i in range(Nt): + src.data[i] = f(x, t[i]) + + src.coordinates.data[:, 0] = x + src_term = src.inject(field=u.forward, expr=src * (dt_symbolic**2)) + stencil = Eq(u.forward, solve(pde, u.forward)) + + # Set up special stencil for initial timestep with substitution for u.backward + v = Function(name='v', grid=grid, npoint=Nx+1, nt=1) + v.data[:] = V(x[:]) + stencil_init = stencil.subs(u.backward, u.forward - dt_symbolic*v) + + # Boundary conditions, depending on arguments + bc = [] + if U_0 is None and U_L is None: + bc += [Eq(u[t_s+1, 0], u[t_s+1, 1])] + else: + if U_0 is not None: + bc += [Eq(u[t_s+1, 0], U_0(t_s+1))] + if U_L is not None: + bc += [Eq(u[t_s+1, L], U_L(t_s+1))] + + op_init = Operator([stencil_init]+src_term+bc) + op = Operator([stencil]+src_term+bc) + + op_init.apply(time_M=1, dt=dt) + + if user_action is not None: + user_action(u.data[1], x, t, 1) + + op.apply(time_m=1,time_M=Nt, dt=dt) + + for n in range(1, Nt+1): + if user_action is not None: + if user_action(u.data[n], x, t, n): + break + + cpu_time = time.perf_counter() - t0 + + return u.data[-1], x, t, cpu_time + +def python_solver(I, V, f, c, U_0, U_L, L, dt, C, T, + user_action=None, version='scalar'): + """ + Solve u_tt=c^2*u_xx + f on (0,L)x(0,T]. + u(0,t)=U_0(t) or du/dn=0 (U_0=None), u(L,t)=U_L(t) or du/dn=0 (u_L=None). + """ + Nt = int(round(T/dt)) + t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time + dx = dt*c/float(C) + Nx = int(round(L/dx)) + x = np.linspace(0, L, Nx+1) # Mesh points in space C2 = C**2; dt2 = dt*dt # Help variables in the scheme # Make sure dx and dt are compatible with x and t dx = x[1] - x[0] @@ -77,10 +181,10 @@ def solver(I, V, f, c, U_0, U_L, L, dt, C, T, u_n = np.zeros(Nx+1) # Solution at 1 time level back u_nm1 = np.zeros(Nx+1) # Solution at 2 time levels back - Ix = range(0, Nx+1) - It = range(0, Nt+1) + Ix = list(range(0, Nx+1)) + It = list(range(0, Nt+1)) - import time; t0 = time.clock() # CPU time measurement + import time; t0 = time.perf_counter() # CPU time measurement # Load initial condition into u_n for i in Ix: @@ -175,7 +279,7 @@ def solver(I, V, f, c, U_0, U_L, L, dt, C, T, # Important to correct the mathematically wrong u=u_nm1 above # before returning u u = u_n - cpu_time = time.clock() - t0 + cpu_time = time.perf_counter() - t0 return u, x, t, cpu_time @@ -228,7 +332,7 @@ def plot_u(u, x, t, n): ext = codec2ext[codec] cmd = '%(movie_program)s -r %(fps)d -i %(filespec)s '\ '-vcodec %(codec)s movie.%(ext)s' % vars() - print cmd + print(cmd) os.system(cmd) return cpu @@ -266,7 +370,7 @@ def assert_no_error(u, x, t, n): solver(I, V, f, c, U_0, U_L, L, dt, C, T, user_action=assert_no_error, version='vectorized') - print U_0, U_L + print(U_0, U_L) def test_quadratic(): """ @@ -361,14 +465,10 @@ def test_plug(): u_s, x, t, cpu = solver( I=I, V=None, f=None, c=0.5, U_0=None, U_L=None, L=L, - dt=dt, C=1, T=4, user_action=None, version='scalar') - u_v, x, t, cpu = solver( - I=I, - V=None, f=None, c=0.5, U_0=None, U_L=None, L=L, - dt=dt, C=1, T=4, user_action=None, version='vectorized') + dt=dt, C=1, T=4, user_action=None) tol = 1E-13 - diff = abs(u_s - u_v).max() - assert diff < tol + # diff = abs(u_s - u_v).max() + # assert diff < tol u_0 = np.array([I(x_) for x_ in x]) diff = np.abs(u_s - u_0).max() assert diff < tol @@ -383,7 +483,7 @@ def guitar(C=1, Nx=50, animate=True, version='scalar', T=2): cpu = viz(I, None, None, c, U_0, U_L, L, dt, C, T, umin=-1.1, umax=1.1, version=version, animate=True) - print 'CPU time: %s version =' % version, cpu + print('CPU time: %s version =' % version, cpu) def moving_end(C=1, Nx=50, reflecting_right_boundary=True, T=2, @@ -412,7 +512,7 @@ def U_0(t): umax = 1.1*0.5 cpu = viz(I, None, None, c, U_0, U_L, L, dt, C, T, umin=-umax, umax=umax, version=version, animate=True) - print 'CPU time: %s version =' % version, cpu + print('CPU time: %s version =' % version, cpu) def sincos(C=1): @@ -451,7 +551,7 @@ def action(u, x, t, n): _dx = L/Nx _dt = C*_dx/c dt.append(_dt) - print dt[-1], E[-1] + print(dt[-1], E[-1]) return dt, E if __name__ == '__main__': diff --git a/fdm-devito-notebooks/02_wave/src-wave/wave1D/wave1D_dn_vc.py b/fdm-devito-notebooks/02_wave/src-wave/wave1D/wave1D_dn_vc.py index 39d33d73..d0377960 100644 --- a/fdm-devito-notebooks/02_wave/src-wave/wave1D/wave1D_dn_vc.py +++ b/fdm-devito-notebooks/02_wave/src-wave/wave1D/wave1D_dn_vc.py @@ -31,8 +31,9 @@ """ import time, glob, shutil, os import numpy as np +from devito import Constant, Grid, TimeFunction, SparseTimeFunction, SparseFunction, Eq, solve, Operator, Buffer, Function -def solver( +def devito_solver( I, V, f, c, U_0, U_L, L, dt, C, T, user_action=None, version='scalar', stability_safety_factor=1.0): @@ -92,7 +93,135 @@ def solver( ('None' if U_L is None else inspect.getsource(U_L)) + \ '_' + str(L) + str(dt) + '_' + str(C) + '_' + str(T) + \ '_' + str(stability_safety_factor) - hashed_input = hashlib.sha1(data).hexdigest() + hashed_input = hashlib.sha1(data.encode('utf-8')).hexdigest() + if os.path.isfile('.' + hashed_input + '_archive.npz'): + # Simulation is already run + return -1, hashed_input + + + grid = Grid(shape=(Nx+1), extent=(L)) + t_s = grid.stepping_dim + + u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2, save=Nt+1) + # Initialise values + u.data[0, :] = I(x[:]) + + if user_action is not None: + user_action(u.data[0], x, t, 0) + + x_dim = grid.dimensions[0] + t_dim = grid.time_dim + + # Initialise variable velocity function + q_f = Function(name='q', grid=grid, npoint=Nx+1, nt=Nt+1) + q_f.data[:] = q[:] + + pde = u.dt2-(q_f*u.dx).dx + stencil = Eq(u.forward, solve(pde, u.forward)) + + # Source term and injection into equation + dt_symbolic = grid.time_dim.spacing + src = SparseTimeFunction(name='f', grid=grid, npoint=Nx+1, nt=Nt+1) + for i in range(Nt): + src.data[i] = f(x, t[i]) + src.coordinates.data[:, 0] = x + src_term = src.inject(field=u.forward, expr=src * (dt_symbolic**2)) + + v = SparseFunction(name='v', grid=grid, npoint=Nx+1, nt=1) + v.data[:] = V(x[:]) + stencil_init = stencil.subs(u.backward, u.forward - 2*dt_symbolic*v) + + # Boundary conditions, depending on arguments + bc = [] + if U_0 is None and U_L is None: + bc += [Eq(u[t_s+1, 0], u[t_s+1, 1])] + else: + if U_0 is not None: + bc += [Eq(u[t_s+1, 0], U_0(t_s+1))] + if U_L is not None: + bc += [Eq(u[t_s+1, L], U_L(t_s+1))] + + op_init = Operator([stencil_init]+bc) + op = Operator([stencil]+src_term+bc) + + cpu = time.perf_counter() + op_init.apply(time_M=1, dt=dt) + + if user_action is not None: + user_action(u.data[1], x, t, 1) + + op.apply(time_m=1,time_M=Nt, dt=dt) + + for n in range(1, Nt+1): + if user_action is not None: + if user_action(u.data[n], x, t, n): + break + + cpu = time.perf_counter() - cpu + + return cpu, hashed_input + +def python_solver( + I, V, f, c, U_0, U_L, L, dt, C, T, + user_action=None, version='scalar', + stability_safety_factor=1.0): + """Solve u_tt=(c^2*u_x)_x + f on (0,L)x(0,T].""" + + # --- Compute time and space mesh --- + Nt = int(round(T/dt)) + t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time + + # Find max(c) using a fake mesh and adapt dx to C and dt + if isinstance(c, (float,int)): + c_max = c + elif callable(c): + c_max = max([c(x_) for x_ in np.linspace(0, L, 101)]) + dx = dt*c_max/(stability_safety_factor*C) + Nx = int(round(L/dx)) + x = np.linspace(0, L, Nx+1) # Mesh points in space + # Make sure dx and dt are compatible with x and t + dx = x[1] - x[0] + dt = t[1] - t[0] + + # Make c(x) available as array + if isinstance(c, (float,int)): + c = np.zeros(x.shape) + c + elif callable(c): + # Call c(x) and fill array c + c_ = np.zeros(x.shape) + for i in range(Nx+1): + c_[i] = c(x[i]) + c = c_ + + q = c**2 + C2 = (dt/dx)**2; dt2 = dt*dt # Help variables in the scheme + + # --- Wrap user-given f, I, V, U_0, U_L if None or 0 --- + if f is None or f == 0: + f = (lambda x, t: 0) if version == 'scalar' else \ + lambda x, t: np.zeros(x.shape) + if I is None or I == 0: + I = (lambda x: 0) if version == 'scalar' else \ + lambda x: np.zeros(x.shape) + if V is None or V == 0: + V = (lambda x: 0) if version == 'scalar' else \ + lambda x: np.zeros(x.shape) + if U_0 is not None: + if isinstance(U_0, (float,int)) and U_0 == 0: + U_0 = lambda t: 0 + if U_L is not None: + if isinstance(U_L, (float,int)) and U_L == 0: + U_L = lambda t: 0 + + # --- Make hash of all input data --- + import hashlib, inspect + data = inspect.getsource(I) + '_' + inspect.getsource(V) + \ + '_' + inspect.getsource(f) + '_' + str(c) + '_' + \ + ('None' if U_0 is None else inspect.getsource(U_0)) + \ + ('None' if U_L is None else inspect.getsource(U_L)) + \ + '_' + str(L) + str(dt) + '_' + str(C) + '_' + str(T) + \ + '_' + str(stability_safety_factor) + hashed_input = hashlib.sha1(data.encode('utf-8')).hexdigest() if os.path.isfile('.' + hashed_input + '_archive.npz'): # Simulation is already run return -1, hashed_input @@ -102,11 +231,11 @@ def solver( u_n = np.zeros(Nx+1) # Solution at 1 time level back u_nm1 = np.zeros(Nx+1) # Solution at 2 time levels back - import time; t0 = time.clock() # CPU time measurement + import time; t0 = time.perf_counter() # CPU time measurement # --- Valid indices for space and time mesh --- - Ix = range(0, Nx+1) - It = range(0, Nt+1) + Ix = list(range(0, Nx+1)) + It = list(range(0, Nt+1)) # --- Load initial condition into u_n --- for i in range(0,Nx+1): @@ -204,7 +333,7 @@ def solver( # Update data structures for next step u_nm1, u_n, u = u_n, u, u_nm1 - cpu_time = time.clock() - t0 + cpu_time = time.perf_counter() - t0 return cpu_time, hashed_input @@ -423,10 +552,10 @@ def make_movie_file(self): os.chdir(directory) # cd directory fps = 24 # frames per second - if self.backend is not None: - from scitools.std import movie - movie('frame_*.png', encoder='html', - output_file='index.html', fps=fps) + # if self.backend is not None: + # from scitools.std import movie + # movie('frame_*.png', encoder='html', + # output_file='index.html', fps=fps) # Make other movie formats: Flash, Webm, Ogg, MP4 codec2ext = dict(flv='flv', libx264='mp4', libvpx='webm', @@ -437,6 +566,7 @@ def make_movie_file(self): ext = codec2ext[codec] cmd = '%(movie_program)s -r %(fps)d -i %(filespec)s '\ '-vcodec %(codec)s movie.%(ext)s' % vars() + print(cmd) os.system(cmd) os.chdir(os.pardir) # move back to parent directory @@ -456,7 +586,7 @@ def close_file(self, hashed_input): archive_name = '.' + hashed_input + '_archive.npz' filenames = glob.glob('.' + self.filename + '*.dat.npz') merge_zip_archives(filenames, archive_name) - print 'Archive name:', archive_name + # print('Archive name:', archive_name) # data = numpy.load(archive); data.files holds names # data[name] extract the array @@ -468,7 +598,7 @@ def demo_BC_plug(C=1, Nx=40, T=4): # Scaled problem: L=1, c=1, max I=1 L = 1. dt = (L/Nx)/C # choose the stability limit with given Nx - cpu, hashed_input = solver( + cpu, hashed_input = python_solver( I=lambda x: 0 if abs(x-L/2.0) > 0.1 else 1, V=0, f=0, c=1, U_0=lambda t: 0, U_L=None, L=L, dt=dt, C=C, T=T, @@ -477,7 +607,7 @@ def demo_BC_plug(C=1, Nx=40, T=4): action.make_movie_file() if cpu > 0: # did we generate new data? action.close_file(hashed_input) - print 'cpu:', cpu + print('cpu:', cpu) def demo_BC_gaussian(C=1, Nx=80, T=4): """Demonstrate u=0 and u_x=0 boundary conditions with a bell function.""" @@ -487,7 +617,7 @@ def demo_BC_gaussian(C=1, Nx=80, T=4): title='u(0,t)=0, du(L,t)/dn=0.', filename='tmpdata') L = 1. dt = (L/Nx)/c # choose the stability limit with given Nx - cpu, hashed_input = solver( + cpu, hashed_input = python_solver( I=lambda x: np.exp(-0.5*((x-0.5)/0.05)**2), V=0, f=0, c=1, U_0=lambda t: 0, U_L=None, L=L, dt=dt, C=C, T=T, @@ -523,7 +653,7 @@ def U_0(t): 'moving_end', -2.3, 2.3, skip_frame=4, title='u(0,t)=0.25*sin(6*pi*t) if t < 1/3 else 0, ' + bc_right, filename='tmpdata') - cpu, hashed_input = solver( + cpu, hashed_input = python_solver( I, V, f, c, U_0, U_L, L, dt, C, T, user_action=action, version=version, stability_safety_factor=1) @@ -676,7 +806,7 @@ def c(x): # Choose the stability limit with given Nx, worst case c # (lower C will then use this dt, but smaller Nx) dt = (L/Nx)/c_0 - cpu, hashed_input = solver( + cpu, hashed_input = devito_solver( I=I, V=None, f=None, c=c, U_0=None, U_L=None, L=L, dt=dt, C=C, T=T, @@ -687,7 +817,7 @@ def c(x): if cpu > 0: # did we generate new data? action.close_file(hashed_input) action.make_movie_file() - print 'cpu (-1 means no new data generated):', cpu + print('cpu (-1 means no new data generated):', cpu) def convergence_rates( u_exact, @@ -720,8 +850,8 @@ def __call__(self, u, x, t, n): E.append(error_calculator.error) h.append(dt) dt /= 2 # halve the time step for next simulation - print 'E:', E - print 'h:', h + print('E:', E) + print('h:', h) r = [np.log(E[i]/E[i-1])/np.log(h[i]/h[i-1]) for i in range(1,num_meshes)] return r @@ -746,9 +876,10 @@ def test_convrate_sincos(): T=1, version='scalar', stability_safety_factor=1.0) - print 'rates sin(x)*cos(t) solution:', \ - [round(r_,2) for r_ in r] + print('rates sin(x)*cos(t) solution:', \ + [round(r_,2) for r_ in r]) assert abs(r[-1] - 2) < 0.002 if __name__ == '__main__': - test_convrate_sincos() + pulse() + # test_convrate_sincos() diff --git a/fdm-devito-notebooks/02_wave/src-wave/wave1D/wave1D_n0.py b/fdm-devito-notebooks/02_wave/src-wave/wave1D/wave1D_n0.py index d8e5875a..ce8224a5 100644 --- a/fdm-devito-notebooks/02_wave/src-wave/wave1D/wave1D_n0.py +++ b/fdm-devito-notebooks/02_wave/src-wave/wave1D/wave1D_n0.py @@ -28,8 +28,88 @@ solution on the screen (as an animation). """ import numpy as np +import time +from devito import Constant, Grid, TimeFunction, SparseTimeFunction, Function, Eq, solve, Operator, Buffer def solver(I, V, f, c, L, dt, C, T, user_action=None): + """Solve u_tt=c^2*u_xx + f on (0,L)x(0,T].""" + Nt = int(round(T/dt)) + t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time + dx = dt*c/float(C) + Nx = int(round(L/dx)) + x = np.linspace(0, L, Nx+1) # Mesh points in space + C2 = C**2 # Help variable in the scheme + + # Make sure dx and dt are compatible with x and t + dx = x[1] - x[0] + dt = t[1] - t[0] + + # Initialising functions f and V if not provided + if f is None or f == 0 : + f = lambda x, t: 0 + if V is None or V == 0: + V = lambda x: 0 + + t0 = time.perf_counter() # Measure CPU time + + # Set up grid + grid = Grid(shape=(Nx+1), extent=(L)) + t_s = grid.stepping_dim + + # Create and initialise u + u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2, save=Nt+1) + u.data[0, :] = I(x[:]) + + if user_action is not None: + user_action(u.data[0], x, t, 0) + + x_dim = grid.dimensions[0] + t_dim = grid.time_dim + + # The wave equation we are trying to solve + pde = (1/c**2)*u.dt2-u.dx2 + + # Source term and injection into equation + dt_symbolic = grid.time_dim.spacing + src = SparseTimeFunction(name='f', grid=grid, npoint=Nx+1, nt=Nt+1) + + for i in range(Nt): + src.data[i] = f(x, t[i]) + + src.coordinates.data[:, 0] = x + src_term = src.inject(field=u.forward, expr=src * (dt_symbolic**2)) + stencil = Eq(u.forward, solve(pde, u.forward)) + + # Set up special stencil for initial timestep with substitution for u.backward + v = Function(name='v', grid=grid, npoint=Nx+1, nt=1) + v.data[:] = V(x[:]) + stencil_init = stencil.subs(u.backward, u.forward - dt_symbolic*v) + + # Boundary condition + bc = [Eq(u[t_s+1, 0], u[t_s+1, 1])] + + # Create and apply operators + op_init = Operator([stencil_init]+src_term+bc) + op = Operator([stencil]+src_term+bc) + + op_init.apply(time_M=1, dt=dt) + + if user_action is not None: + user_action(u.data[1], x, t, 1) + + op.apply(time_m=1, time_M=Nt, dt=dt) + + for n in range(1, Nt+1): + if user_action is not None: + if user_action(u.data[n], x, t, n): + break + + cpu_time = time.perf_counter() - t0 + + return u.data[-1], x, t, cpu_time + + +def python_solver(I, V, f, c, L, dt, C, T, user_action=None): """ Solve u_tt=c^2*u_xx + f on (0,L)x(0,T]. u(0,t)=U_0(t) or du/dn=0 (U_0=None), u(L,t)=U_L(t) or du/dn=0 (u_L=None). @@ -54,7 +134,7 @@ def solver(I, V, f, c, L, dt, C, T, user_action=None): u_n = np.zeros(Nx+1) # Solution at 1 time level back u_nm1 = np.zeros(Nx+1) # Solution at 2 time levels back - import time; t0 = time.clock() # CPU time measurement + import time; t0 = time.perf_counter() # CPU time measurement # Load initial condition into u_n for i in range(0, Nx+1): @@ -96,7 +176,7 @@ def solver(I, V, f, c, L, dt, C, T, user_action=None): # Wrong assignment u = u_nm1 must be corrected before return u = u_n - cpu_time = time.clock() - t0 + cpu_time = time.perf_counter() - t0 return u, x, t, cpu_time from wave1D_u0 import viz diff --git a/fdm-devito-notebooks/02_wave/src-wave/wave1D/wave1D_u0.py b/fdm-devito-notebooks/02_wave/src-wave/wave1D/wave1D_u0.py index afc20030..26f72057 100644 --- a/fdm-devito-notebooks/02_wave/src-wave/wave1D/wave1D_u0.py +++ b/fdm-devito-notebooks/02_wave/src-wave/wave1D/wave1D_u0.py @@ -54,8 +54,11 @@ def solver(I, V, f, c, L, dt, C, T, user_action=None): t_s = grid.stepping_dim # Create and initialise u - u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2) - u.data[:,:] = I(x[:]) + u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2, save=Nt+1) + u.data[0, :] = I(x[:]) + + if user_action is not None: + user_action(u.data[0], x, t, 0) x_dim = grid.dimensions[0] t_dim = grid.time_dim @@ -88,8 +91,17 @@ def solver(I, V, f, c, L, dt, C, T, user_action=None): op = Operator([stencil]+src_term+bc) op_init.apply(time_M=1, dt=dt) + + if user_action is not None: + user_action(u.data[1], x, t, 1) + op.apply(time_m=1, time_M=Nt, dt=dt) + for n in range(1, Nt+1): + if user_action is not None: + if user_action(u.data[n], x, t, n): + break + cpu_time = time.perf_counter() - t0 return u.data[-1], x, t, cpu_time diff --git a/fdm-devito-notebooks/02_wave/wave1D_fd2.ipynb b/fdm-devito-notebooks/02_wave/wave1D_fd2.ipynb index 091094a1..e10da9de 100644 --- a/fdm-devito-notebooks/02_wave/wave1D_fd2.ipynb +++ b/fdm-devito-notebooks/02_wave/wave1D_fd2.ipynb @@ -91,10 +91,7 @@ "\n", "\n", "How can we incorporate the condition ([1](#wave:pde1:Neumann:0))\n", - "in the finite difference scheme? Since we have used central\n", - "differences in all the other approximations to derivatives in the\n", - "scheme, it is tempting to implement ([1](#wave:pde1:Neumann:0)) at\n", - "$x=0$ and $t=t_n$ by the difference" + "in the finite difference scheme? We can see that, on the boundary" ] }, { @@ -106,7 +103,7 @@ "\n", "$$\n", "\\begin{equation}\n", - "[D_{2x} u]^n_0 = \\frac{u_{-1}^n - u_1^n}{2\\Delta x} = 0\n", + "u_x \\approx \\frac{u_1^n - u_0^n}{\\Delta x} = 0\n", "\\thinspace .\n", "\\label{wave:pde1:Neumann:0:cd} \\tag{2}\n", "\\end{equation}\n", @@ -117,190 +114,272 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The problem is that $u_{-1}^n$ is not a $u$ value that is being\n", - "computed since the point is outside the mesh. However, if we combine\n", - "([2](#wave:pde1:Neumann:0:cd)) with the scheme\n", - "" + "Therefore, we can pass the following equation as a boundary condition to the Devito `Operator`:\n", + "```python\n", + "bc = [Eq(u[t_s+1, 0], u[t_s+1, 1])]\n", + "```\n", + "## Solver\n", + "Now, let's take a look at the full Devito solver function. The solver is very similar to the one in the waves [Implementation](wave1D_prog.ipynb) section, except for the boundary condition." ] }, { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "
\n", - "\n", - "$$\n", - "\\begin{equation}\n", - "u^{n+1}_i = -u^{n-1}_i + 2u^n_i + C^2\n", - "\\left(u^{n}_{i+1}-2u^{n}_{i} + u^{n}_{i-1}\\right),\n", - "\\label{wave:pde1:Neumann:0:scheme} \\tag{3}\n", - "\\end{equation}\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "for $i=0$, we can eliminate the fictitious value $u_{-1}^n$. We see that\n", - "$u_{-1}^n=u_1^n$ from ([2](#wave:pde1:Neumann:0:cd)), which\n", - "can be used in ([3](#wave:pde1:Neumann:0:scheme)) to\n", - "arrive at a modified scheme for the boundary point $u_0^{n+1}$:" - ] - }, - { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": 1, "metadata": {}, + "outputs": [], "source": [ - "\n", - "
\n", - "\n", - "$$\n", - "\\begin{equation}\n", - "u^{n+1}_i = -u^{n-1}_i + 2u^n_i + 2C^2\n", - "\\left(u^{n}_{i+1}-u^{n}_{i}\\right),\\quad i=0 \\thinspace . \\label{_auto1} \\tag{4}\n", - "\\end{equation}\n", - "$$" + "# NBVAL_IGNORE_OUTPUT\n", + "import numpy as np\n", + "import time as time\n", + "import matplotlib.pyplot as plt\n", + "from devito import Constant, Grid, TimeFunction, SparseTimeFunction, Function, Eq, solve, Operator, Buffer" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": 2, "metadata": {}, + "outputs": [], "source": [ - "[Figure](#wave:pde1:fig:Neumann:stencil) visualizes this equation\n", - "for computing $u^3_0$ in terms of $u^2_0$, $u^1_0$, and\n", - "$u^2_1$.\n", + "def python_solver(I, V, f, c, L, dt, C, T, user_action=None):\n", + " \"\"\"\n", + " Solve u_tt=c^2*u_xx + f on (0,L)x(0,T].\n", + " u(0,t)=U_0(t) or du/dn=0 (U_0=None), u(L,t)=U_L(t) or du/dn=0 (u_L=None).\n", + " \"\"\"\n", + " Nt = int(round(T/dt))\n", + " t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time\n", + " dx = dt*c/float(C)\n", + " Nx = int(round(L/dx))\n", + " x = np.linspace(0, L, Nx+1) # Mesh points in space\n", + " C2 = C**2; dt2 = dt*dt # Help variables in the scheme\n", + " # Make sure dx and dt are compatible with x and t\n", + " dx = x[1] - x[0]\n", + " dt = t[1] - t[0]\n", "\n", - "\n", - "\n", - "
\n", + " # Wrap user-given f, V\n", + " if f is None or f == 0:\n", + " f = (lambda x, t: 0)\n", + " if V is None or V == 0:\n", + " V = (lambda x: 0)\n", "\n", - "

Modified stencil at a boundary with a Neumann condition.

\n", - "\n", + " u = np.zeros(Nx+1) # Solution array at new time level\n", + " u_n = np.zeros(Nx+1) # Solution at 1 time level back\n", + " u_nm1 = np.zeros(Nx+1) # Solution at 2 time levels back\n", "\n", - "\n", + " import time; t0 = time.perf_counter() # CPU time measurement\n", "\n", + " # Load initial condition into u_n\n", + " for i in range(0, Nx+1):\n", + " u_n[i] = I(x[i])\n", + " print(u_n)\n", "\n", - "Similarly, ([1](#wave:pde1:Neumann:0)) applied at $x=L$\n", - "is discretized by a central difference" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "
\n", + " if user_action is not None:\n", + " user_action(u_n, x, t, 0)\n", "\n", - "$$\n", - "\\begin{equation}\n", - "\\frac{u_{N_x+1}^n - u_{N_x-1}^n}{2\\Delta x} = 0\n", - "\\thinspace .\n", - "\\label{wave:pde1:Neumann:0:cd2} \\tag{5}\n", - "\\end{equation}\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Combined with the scheme for $i=N_x$ we get a modified scheme for\n", - "the boundary value $u_{N_x}^{n+1}$:" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "
\n", + " # Special formula for the first step\n", + " for i in range(0, Nx+1):\n", + " ip1 = i+1 if i < Nx else i-1\n", + " im1 = i-1 if i > 0 else i+1\n", + " u[i] = u_n[i] + dt*V(x[i]) + \\\n", + " 0.5*C2*(u_n[im1] - 2*u_n[i] + u_n[ip1]) + \\\n", + " 0.5*dt2*f(x[i], t[0])\n", + " print(u)\n", "\n", - "$$\n", - "\\begin{equation}\n", - "u^{n+1}_i = -u^{n-1}_i + 2u^n_i + 2C^2\n", - "\\left(u^{n}_{i-1}-u^{n}_{i}\\right),\\quad i=N_x \\thinspace . \\label{_auto2} \\tag{6}\n", - "\\end{equation}\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The modification of the scheme at the boundary is also required for\n", - "the special formula for the first time step. How the stencil moves\n", - "through the mesh and is modified at the boundary can be illustrated by\n", - "an animation in a [web page](${doc_notes}/book/html/mov-wave/N_stencil_gpl/index.html)\n", - "or a [movie file](${docraw}/mov-wave/N_stencil_gpl/movie.ogg).\n", + " if user_action is not None:\n", + " user_action(u, x, t, 1)\n", "\n", + " # Update data structures for next step\n", + " #u_nm1[:] = u_n; u_n[:] = u # safe, but slower\n", + " u_nm1, u_n, u = u_n, u, u_nm1\n", "\n", + " for n in range(1, Nt):\n", + " for i in range(0, Nx+1):\n", + " ip1 = i+1 if i < Nx else i-1\n", + " im1 = i-1 if i > 0 else i+1\n", + " u[i] = - u_nm1[i] + 2*u_n[i] + \\\n", + " C2*(u_n[im1] - 2*u_n[i] + u_n[ip1]) + \\\n", + " dt2*f(x[i], t[n])\n", "\n", + " if user_action is not None:\n", + " if user_action(u, x, t, n+1):\n", + " break\n", "\n", - "## Implementation of Neumann conditions\n", - "
\n", + " # Update data structures for next step\n", + " #u_nm1[:] = u_n; u_n[:] = u # safe, but slower\n", + " u_nm1, u_n, u = u_n, u, u_nm1\n", "\n", - "We have seen in the preceding section\n", - "that the special formulas for the boundary points\n", - "arise from replacing $u_{i-1}^n$ by $u_{i+1}^n$ when computing\n", - "$u_i^{n+1}$ from the stencil formula for $i=0$. Similarly, we\n", - "replace $u_{i+1}^n$ by $u_{i-1}^n$ in the stencil formula\n", - "for $i=N_x$. This observation can conveniently\n", - "be used in the coding: we just work with the general stencil formula,\n", - "but write the code such that it is easy to replace `u[i-1]` by\n", - "`u[i+1]` and vice versa. This is achieved by\n", - "having the indices `i+1` and `i-1` as variables `ip1` (`i` plus 1)\n", - "and `im1` (`i` minus 1), respectively.\n", - "At the boundary we can easily define `im1=i+1` while we use\n", - "`im1=i-1` in the internal parts of the mesh. Here are the details\n", - "of the implementation (note that the updating formula for `u[i]`\n", - "is the general stencil formula):" + " # Wrong assignment u = u_nm1 must be corrected before return\n", + " u = u_n\n", + " cpu_time = time.perf_counter() - t0\n", + " return u, x, t, cpu_time" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": 3, "metadata": {}, + "outputs": [], "source": [ - " i = 0\n", - " ip1 = i+1\n", - " im1 = ip1 # i-1 -> i+1\n", - " u[i] = u_n[i] + C2*(u_n[im1] - 2*u_n[i] + u_n[ip1])\n", + "# %load -s solver, src-wave/wave1D/wave1D_n0.py\n", + "def solver(I, V, f, c, L, dt, C, T, user_action=None):\n", + " \"\"\"Solve u_tt=c^2*u_xx + f on (0,L)x(0,T].\"\"\"\n", + " Nt = int(round(T/dt))\n", + " t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time\n", + " dx = dt*c/float(C)\n", + " Nx = int(round(L/dx))\n", + " x = np.linspace(0, L, Nx+1) # Mesh points in space\n", + " C2 = C**2 # Help variable in the scheme\n", + " \n", + " # Make sure dx and dt are compatible with x and t\n", + " dx = x[1] - x[0]\n", + " dt = t[1] - t[0]\n", + "\n", + " # Initialising functions f and V if not provided\n", + " if f is None or f == 0 :\n", + " f = lambda x, t: 0\n", + " if V is None or V == 0:\n", + " V = lambda x: 0\n", " \n", - " i = Nx\n", - " im1 = i-1\n", - " ip1 = im1 # i+1 -> i-1\n", - " u[i] = u_n[i] + C2*(u_n[im1] - 2*u_n[i] + u_n[ip1])\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We can in fact create one loop over both the internal and boundary\n", - "points and use only one updating formula:" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - " for i in range(0, Nx+1):\n", - " ip1 = i+1 if i < Nx else i-1\n", - " im1 = i-1 if i > 0 else i+1\n", - " u[i] = u_n[i] + C2*(u_n[im1] - 2*u_n[i] + u_n[ip1])\n" + " t0 = time.perf_counter() # Measure CPU time\n", + "\n", + " # Set up grid\n", + " grid = Grid(shape=(Nx+1), extent=(L))\n", + " t_s = grid.stepping_dim\n", + " \n", + " # Create and initialise u\n", + " u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2)\n", + " for i in range(Nx+1):\n", + " u.data[:,i] = I(x[i])\n", + " print(u.data[0])\n", + "\n", + " x_dim = grid.dimensions[0]\n", + " t_dim = grid.time_dim\n", + " \n", + " # The wave equation we are trying to solve\n", + " pde = (1/c**2)*u.dt2-u.dx2\n", + " \n", + " # Source term and injection into equation\n", + " dt_symbolic = grid.time_dim.spacing \n", + " src = SparseTimeFunction(name='f', grid=grid, npoint=Nx+1, nt=Nt+1)\n", + " \n", + " for i in range(Nt):\n", + " src.data[i] = f(x, t[i])\n", + " \n", + " src.coordinates.data[:, 0] = x\n", + " src_term = src.inject(field=u.forward, expr=src * (dt_symbolic**2))\n", + " stencil = Eq(u.forward, solve(pde, u.forward))\n", + "\n", + " # Set up special stencil for initial timestep with substitution for u.backward\n", + " v = Function(name='v', grid=grid, npoint=Nx+1, nt=1)\n", + " v.data[:] = V(x[:])\n", + " stencil_init = stencil.subs(u.backward, u.forward - dt_symbolic*v)\n", + "\n", + " # Boundary condition\n", + " bc = [Eq(u[t_s+1, 0], u[t_s+1, 1])]\n", + " bc += [Eq(u[t_s+1, Nx], u[t_s+1, Nx-1])]\n", + "\n", + " # Create and apply operators\n", + " op_init = Operator([stencil_init]+bc+src_term)\n", + " op = Operator([stencil]+bc+src_term)\n", + " \n", + " op_init.apply(time_M=1, dt=dt)\n", + " \n", + " op.apply(time_m=1, time_M=Nt, dt=dt)\n", + " print(\"u.data\", u.data)\n", + " cpu_time = time.perf_counter() - t0\n", + " \n", + " return u.data[1], x, t, cpu_time\n" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": 4, "metadata": {}, - "source": [ - "The program [`wave1D_n0.py`](${src_wave}/wave1D/wave1D_n0.py)\n", - "contains a complete implementation of the 1D wave equation with\n", - "boundary conditions $u_x = 0$ at $x=0$ and $x=L$.\n", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0. 0. 0. 0. 1. 1. 0. 0. 0. 0. 0.]\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Data type float64 of runtime value `dt` does not match the Constant data type \n", + "Operator `Kernel` run in 0.01 s\n", + "Data type float64 of runtime value `dt` does not match the Constant data type \n", + "Operator `Kernel` run in 0.01 s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "u.data [[1. 1. 0. 0. 0. 0. 0. 0. 1. 0. 0.]\n", + " [0. 0. 0. 1. 0. 0. 1. 0. 0. 0. 0.]\n", + " [0. 0. 1. 0. 0. 0. 0. 1. 0. 0. 0.]]\n", + "[0. 0. 0. 0. 1. 1. 0. 0. 0. 0. 0.]\n", + "[0. 0. 0. 0.5 0.5 0.5 0.5 0. 0. 0. 0. ]\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8KklEQVR4nO2deXxb5ZX3f0e75H1NnNiO7cRxAmQhcQhLadNSILRAmJYXUlq2Fiid0oVCt+lGmW4z0GkhMLRpCx36TllKO5C+wEBbyKe0bAkNCSSxE9txEjuL13jTLj3vH1os3SvZsq3nSro638+HD76LdB/Z0f3d85zz/A4JIcAwDMPkL4ZMD4BhGIbJLCwEDMMweQ4LAcMwTJ7DQsAwDJPnsBAwDMPkOaZMD2CmVFZWioaGhkwPg2EYJqd46623BoQQVYmO5ZwQNDQ0YOfOnZkeBsMwTE5BRIeTHeOpIYZhmDyHhYBhGCbPYSFgGIbJc3IuR8AwDAMAPp8PPT09cLvdmR5KVmGz2VBbWwuz2Zzya1gIGIbJSXp6elBUVISGhgYQUaaHkxUIITA4OIienh40Njam/DppU0NE9DAR9RHRu0mOExHdT0QdRLSHiNbIGgvDMPrD7XajoqKCRSAGIkJFRcWMoySZEcGvATwA4NEkxy8B0Bz+bz2Ah8L/Z5hZcXzEhXtfOICjQ86MXN9qNuBDK2qweV0d35w0gn/PambzO5EmBEKIvxJRwxSnbALwqAj5YL9ORKVEVCOEOC5rTIy++cJjb+PN7qGMjuGVgwNYWGrHe5cmXLfDMFlJJquGFgI4GrPdE96ngohuIaKdRLSzv79fk8ExM+fkb0+i/ZZ2HP3xUQw+OwhnhxNBf1CTa7t9Aew4nFkRiPD3joFMD4HRCCLCHXfcEd2+9957cdddd6Xlve+66y4sXLgQq1evxhlnnIFt27al5X0TkRPJYiHEVgBbAaC1tZU76WQpQy8M4eSjJ+P2kZmw9OdLUXNjjdRrHxqYQLb0WOrsn5B+jdE3RzH0/BDsLXY4ljngWOqA0WGUfl0mHqvVij/84Q/4+te/jsrKyrS//+23344777wT+/fvx/nnn4++vj4YDOl/fs+kEPQCqIvZrg3vY3IUZ7t6bl74BCzzLNKv3dk/Hre9dlEZvnJxi/TrAsDhQSe+8vs90e0uxVhkMPznYXTf1R23z1pvRe3ttaj7Yl3iF+mUhq89K/X9u3/04aTHTCYTbrnlFvzkJz/B97///bhjN9xwAy699FJceeWVAIDCwkKMj49j+/bt+M53voPS0lK88847uOqqq7BixQrcd999cLlcePrpp7F48eK491q+fDlMJhOOHj2KDRs24MCBAzCbzRgdHcWqVaui27Mlk0KwDcBtRPQ4QkniEc4P5C5CCDjbEidpHS0O6dfvUjyFr6otxfqmCunXBYDTFhTHCcGRISd8gSDMRnkzr4lE13PEA2gzE8fE8NnPfhYrV67EV77ylZRfs3v3buzfvx/l5eVoamrCTTfdhDfffBP33XcftmzZgp/+9Kdx57/xxhswGAyor6/Hhg0b8Oyzz+KKK67A448/jo985CNzEgFAbvnoYwBeA9BCRD1E9CkiupWIbg2f8hyALgAdAH4B4J9ljYWRj6/Ph8BIQLWfLARbg0369ZURQVNVgfRrRiiymVFdZI1u+4MChwflVi4lEgIAcCyTL7pMPMXFxbjuuutw//33p/yadevWoaamBlarFYsXL8ZFF10EAFixYgW6u7uj5/3kJz/B6tWrceedd+KJJ54AEeGmm27CI488AgB45JFHcOONN875M8isGvrYNMcFgM/Kuj6jLcmiAXuzHWSUX+KnjAgWVxVKv6byen1jnpjxjGNJtZwxZDr6YtR88YtfxJo1a+JuyiaTCcFgKEQLBoPwer3RY1br5IODwWCIbhsMBvj9/uixSI4glvPOOw/d3d3Yvn07AoEAzjjjjDmPPyeSxUz2U3BGAU7/w+lwtjnhbHfC2eaEq92lyY1JCKGal1+sYUQAhCKQ17oGo9tdA/ISxpmOvrKNqebwtaK8vBxXXXUVfvWrX+GTn/wkgJBl/ltvvYWrrroK27Ztg8/nS9v1rrvuOlxzzTX41re+lZb3YyFg0oK5woyqf4qvnRdCIOiWP2l9ctSDCe/kjbHQakJVzFSNFjQpIpDOPnkJY7ISlmxZAle7Kyq8nqOeaPQlwuVTvNhKW+644w488MAD0e2bb74ZmzZtwqpVq7Bx40YUFKTv4eTjH/84vvnNb+JjH5ty4iVlSGRLzV2KtLa2Cm5Mw8TyascArvnlG9HtVbUleOa292g6hu3tfbjhkR3R7bWLyvD7z5wr/brODicm3pnA+D/GceqVUxAeAWe7EyufX4ni9cXSr59J9u/fj+XLl2d6GBnhqaeewjPPPIPf/OY3CY8n+t0Q0VtCiNZE53NEwKQd/6g/Oj0UO020+pXVMJfOrbohEepEsbb5AUCdk1COSRYdX+zA0LPqhXTONqfuhSBf+dznPofnn38ezz33XNrek4WASTs71+yEu1NteuVqd8G8XoYQKBPF2uYHAGBBqR1WkwGe8ErqU04fhia8KC+Qu4bCscyRWAiSVBUxuc+WLVvS/p7cmIZJO8lKGGXdnLIhIjAaCI2V8QKkRVSQLBmfrKqIYRLBQsCkHa1vTpkuHU12XS1WGGstuow+YSFg0o6WNyeXN4BjI67oNhGwqCIztfTKRWxKgZJBMtH1HPFABHKrEITJHJwjYNKOlhGB0myutswOmzkz5mtKIdBiashcZUb5JeWwLrTCscwRMqFrccDWaNNkIR+jD1gImDkjhIirWXcsc8DWZIOjxRFyxmxxRF0y003XgHIhWWamhRJdW4uIgIiw8rmV0q/DJMZoNGLFihXw+XwwmUy47rrrcPvtt8/KIXTnzp149NFHcf/992P79u2wWCw491z5JcgACwEzR4QQeLXmVViqLHFPpKc9fhqK1hRJfyrt7Iu/2TZVZk4IlMniw0NOeP1BWEw8A6sF22l7wv0bxIZpX+s55sFrC1+b8evsdjvefvttAEBfXx+uueYajI6O4rvf/e60r1XS2tqK1tZQmf/27dtRWFiomRDwv1BmTnhPeuE76cPEuxPof6ofR75/BG3XtWHX+bs0ub4qIqjWvnQ0QpHNjHnFkyuaA0GBI2lum/nuR9/F7o27cfCLB9H7UC+GXx6G55gHubYwVI9UV1dj69ateOCBByCEQCAQwJe//GWsW7cOK1euxM9//nMAwObNm/Hss5PW2TfccAOeeuopbN++HZdeeim6u7vxs5/9LGo498orr6C7uxsf+MAHsHLlSlxwwQU4cuRIWsfOEQEzJ5KanzU7MmI2l8mIIHL9k6NyzOeEEBj+yzACIwEMvzAcd2x9x3rYF9vTch1m9jQ1NSEQCKCvrw/PPPMMSkpKsGPHDng8Hpx33nm46KKLcPXVV+PJJ5/Ehz/8YXi9XvzlL3/BQw89hDfeCK2Ob2howK233orCwsKo4dxll12G66+/Htdffz0efvhhfP7zn8fTTz+dtnFzRMDMCVe7K+F+e4v8m1I2mM0pUSeM05cn8J70stlcDvHiiy/i0UcfxerVq7F+/XoMDg7i4MGDuOSSS/Dyyy/D4/Hg+eefx3vf+17Y7VN/X1577TVcc801AIBrr70Wf/vb39I6Vo4ImDmR1BdfA9dRpdlcUQbM5pTIXEuQVHQ1svpmpqerqwtGoxHV1dUQQmDLli24+OKLVedt2LABL7zwAp544gls3rw5AyONhyMCZk64u9VWEoA2QpCoGU2mHTdllpByM5rspr+/H7feeituu+02EBEuvvhiPPTQQ1H76QMHDmBiIhQhXn311XjkkUfwyiuvYOPGjar3KioqwtjYWHT73HPPxeOPPw4A+O///m+cf/75aR07RwTMnDj996fDe8IbZy7nbHOicJX8uXr1tFBm8wOJxtDZP6Eqr50t3IxmalKp8kmGdYF1Vq93uVxYvXp1tHz02muvxZe+9CUAwE033YTu7m6sWbMGQghUVVVF5/UvuugiXHvttdi0aRMsFrUf1WWXXYYrr7wSzzzzDLZs2YItW7bgxhtvxD333IOqqqpoh7J0wULAzAkigrXGCmuNFWUbyjS9tnL+Xcv2lMlYqDCfG3GFzOcqCuc+ZdXw3QbM+/i8OFdXZ7sTBSvUnzvoC8J9yA0REChYnvnfi14JBNQ5mwgGgwE/+MEP8IMf/EB1zGw2Y2go3ixww4YN2LBhAwBg6dKl2LNnT9zxl156ae4DTgILAaMp/nE/hFfAXD53F1LltEs2RASGsPlc24nJsL5rYCItQmAqMqFobRGK1hapjjnbnTj+8PGoOLg73RB+gcorKnHG/8y9lSGjb1gIGGlMtE1g+MXh0NNr+AnW2+tF/dfq0fTDpjm/v6p0NAuEAAgJUpwQ9I9jXUO51Gt6+7w4+u9HVfvZfI5JBRYCRhojr4yg4wsdqv3puDm5vAH0nsoOszklMktIk5EsT+DqcCHoD8Kg09XN6cq/6InZLC7U578OJitI6kKaBvO5Q4rm8HVljoyZzSnJhB21ucoMU5n6uU74BNyHEld25To2mw2Dg4O8qjoGIQQGBwdhs81sXQlHBIw0ZD6lJiodzRYyEREQERwtDoy+Pqo65mx3wtGcHdFSOqmtrUVPTw/6+/szPZSswmazoba2dkavYSFgpGGuMsNUaoL/lD9uf+QpdS43p2xpRpMIZa7iiEbmc45l8UJgqjDB0eKAwazPwN9sNqOxsTHTw9AFLASMNIhIdXMCAdZ6K3yDPqB59u+tNJvLpoig0GrCvGJr1HMoYj6XLs+hZMy7fh5Kzi+J2n5bKuX2S2b0AwsBMyu8J73w9HpgX2qHqTD5P6P5N8xH+YfKQ70JWhywN9thdMx9Ll81NZRhszklSvO5zjSazyWjbEMZsEHqJRidwkLAzIqBpwdw4NYDAADLQku0CU3FpRWouKQiet6CTy9I+7VDZnOKqaEM2k8nYnF1AV7rGoxuz6VJTcAdQGAkAHO1mStkGCmwEDCzIrYE1NvrhbfXi1MvnYKpzBQnBDI4MeqGU2k2l4YFW+lEGaHMxXNo9PVR7H7/bhhLjHFd34rXF6PsAm1XczP6hIWAmRWZ9L1RLSSrLsy6J+XF1ekrIY24jgZGAhh7cwxjb4YWq1V+tJKFgEkL+iwnYKSTSSdMlbVEZXZNCwFAU6W6hHS29e6ZtPpm8gMWAmbGBNyBpIuUMhERKJ++s4GI+VyEiPncbEgafbH9NJMmpAoBEW0konYi6iCiryU4Xk9ELxPRLiLaQ0QfkjkeJj0ERgIou7AM1kXx8/KWGgtMxfJnG9UVQ9kXEUTM52LpGphdwpgjAkY20r61RGQE8CCACwH0ANhBRNuEEPtiTvsmgCeFEA8R0WkAngPQIGtMTHqwzLNg1QurAAABZwCug6EeBAFXckvedJKtZnNKlOZznX0zN58TAQFjoRFkIQhv/NQSCwGTLmQ+vp0FoEMI0QUARPQ4gE0AYoVAACgO/1wC4JjE8TASMDqMKFxVqEkjGkBtNmfIIrM5Jcr+ybOJCMhIWLd7HURAwH3YHbWZ9vR4YCqZ+uvrG/bB2R5uFnTAiYbvNujWfI6ZGzKFYCGAWF/cHgDrFefcBeBFIvocgAIAH0z0RkR0C4BbAKC+vj7tA2W0QwgBT48HznYnvCe8mP+J+TN6vXJFcW0Wmc0pUUYqnX2zrxwiI8HeZIe9yY6KD01dnvvOpncw+voofH2+uP3zb5ivS88hZu5kunz0YwB+LYT4MRGdA+A3RHSGECIYe5IQYiuArQDQ2trKVoM5RsAdQPsn20NPswecCE6E/rxkJlRvrp7RU6raYyj78gMRVC6ks8wRzBRfv08lAoB+zeeYuSMzTuwFUBezXRveF8unADwJAEKI1wDYAFRKHBOTAQxWA4b+dwjju8ajIgCEzee6Z2aRnCv5AQBoVIhUxHxONslyB+mw/2b0iUwh2AGgmYgaicgCYDOAbYpzjgC4AACIaDlCQsCesjojYpGciJnenLLZflpJxHwuQsh8ToMmNUnKSiML0xhGiTQhEEL4AdwG4AUA+xGqDtpLRHcT0eXh0+4AcDMR7QbwGIAbBHeZ0CXpujkpcwTZZD+dCOX4tOhNYG+xJ9zPEQGTDKk5AiHEcwiVhMbu+3bMz/sAnCdzDEx2kI6bUyKzuWyOCIDQ+F7tnDSfm4vnUKpERdeAUIK5xQ7HMgeK1qib3jMMkPlkMZMnxEYEBocBjqUh87Ti84qneFU8KrM5W/aZzSlRt63UICJYYse6vetgX2yHwcrlosz0sBAwM6L3oV4giOhTpnWhNSXDt5L3lGDliyvhaHHAWmsFGWZuEtfZp04UZ5vZnBJVCakGEYHBZEDBadkdKTHZBQsBMyOO3nM0zmfIUBB6uj/996fD3ph4+gcALJUWlF84s1W1StT5gey/2akWlYXN51IRsOGXhuEf9cOxzAF7kx0GCz/dM3JgIWBSJuAOqMo9gxNBjO8ah7nCLP362dynOBkLSuywmQ1w+0JloxHzuYoUprSO3nsUQ88PhTaMofl+xzIHGu5q4Pl+Jq3wIwaTMq6DrpApiAI2m0uOwUBoqFBbUqdCXCI9EPr9D/5xECLAhXVMemEhYFIm0y6YuWA/nYjZNKlJFH1FYLM5Jt2wEDApk6zmXwtffKfXnzNmc0qUjXNSSRhnOvpi8gv+F8WkTOkHSrHoO4vgbJt0tAw6g0nXCKSTQwqfnrpyB6ym7DSbU6KOCKafGsp09MXkFywETMqUnFOCknNKotsiKODp9cBglx9YqhaS5UB+IIKykX0q5nPWWitqbqqJ2k77+kMmctyVjJEBCwEza8hAsNXZNLmW2mMoN/IDQHLzOcsUrqslZ5eg5OxJ0fUNhXoLTNeDIJZoD4N252RfgjYnqq+pxoKbF8z8gzC6hYWA0ZyAe7KrmbPdCWebEws+vQCl55cmfU0ulo5GKLSaML/YhhOjoeRvxHxuSXXqJaDmcnNcNJYKR398FF1f7VLtty+xAzfP6K0YncNCwGhK55c7cfTHR1WJ0MLVhVMKQS65jiaiqaogKgQA0NE3MyGYDfalbD7HpAZXDTGaYq40J6yGmcqFNBhUm83lUkQAJGpSo6H5nIJkiWgmf2EhYDRlNi6kJ0bdcPnizeYqCy1pH5tMlBGM0jdJBvYmO5CgsMrX74NvSN3BjMlfWAgYTUnaoGaKp9RE0UC2m80pyUREYLAYYF+cRHg5KmBi4BwBoyn2xeGn1ED8/shTqrlc7VmkvGnmWn4AUI95JuZzc6FoXRFMJSY4WkK23/YWe+hnXo/AxMBCwGhK5Ck1MB6I3pwi/zc4EgeonX251ZUsEYnM5wYnvKiU3E/htP97mtT3Z/QBCwEzLRP7J9B2XZvqidLebIfRPvPVvev2rJtRwxTlAqxcsJ9WYjAQGisLsf/4aHRfV/+EdCFgmFRgIWCmZWLvBMZ2jmFs51jc/tL3l2L1S6tn/H4z7ZqljAhyaTFZLE1VBXFC0Nk/jrMa43s0BNwB7P/Y/jjBdSxzJJwyY5h0wULATEuyih4t5pmdXj+OjcQ0wskhszkl6raV6oSxq8OFgacHVPsLVhRg3Z510sbG5DdcNcRMSyZdR3PZbE5Jom5lSpKJrhaNf5j8hYWAmZZkNyctXEeVTVxyyWxOidJ8LpEddSZFl8lfWAiYKRFCJLdE1uDmpJw+ycWKoQjKEtKjwy54/PF1tJkUXSZ/4RwBMy2tu1sn3SvDJnHuw27Y6uU7j6oighwWgoJE5nODTjTPm/QcyqToMvkLCwEzJUQEe6Md9kY7sFH766sjgtydGgKAxdXx5nOd/RNxQtC8pRkT705EBdfZ7oSr05X2xLwICvgGfLBU55ZVByMHFgImqxBBATKEVtsmMpvL5YgACOUJ/t4xGN1W5gmK1xejeH1x3L6gNwgyz34FctATRP/v++PF5YALxgIjzus/b9bvy+gHFgImIwRcAYztGIvrSeBsd8Jaa8WZ288EoDabK85BszklqVQOKTFY5pjKMwBt17dB+ONtX4OuIHyDPq5IYlgImMzgPuzG2+97W7U/MDJ5408UDeSa2ZwSZUSjifmc2QBbkw2uA+qKJGe7EyXnzqzhDaM/uGqIyQhJLZIHJi2Sc70ZTSLUdtTjECJBg4Y0MxvXVyZ/YCFgMkIqFsl6Kh2NEDGfizDq9mNwwiv9ukmb1HC3MgaShYCINhJROxF1ENHXkpxzFRHtI6K9RPRbmeNhsoukT6nhm5OydDTXK4aASfO5WJReSjJI9Ls2FBggfPKjESb7kZYjICIjgAcBXAigB8AOItomhNgXc04zgK8DOE8IMUxE1bLGw2QfjmUODP5xULXfc8QDQJ8RARAStDgX0oEJrG+qkHrNovVFWPj5hXFGdpYFlpzPuTDpQWay+CwAHUKILgAgoscBbAKwL+acmwE8KIQYBgAhRJ/E8TAzxDfkg7HAOGO30FQpPrcYlR+pjOtJYF9qh7nUnNBsrj5HzeaUKBPGWkQEhWcUovm+ZunXYXITmUKwEMDRmO0eAOsV5ywFACL6O0Kpw7uEEP+rfCMiugXALQBQX18vZbCMmoOfP4i+x/pga7TFNZEpv6Qctrq5ryquuqIKVVdUJTymrBjKZbM5JaoS0oGJuPUTDKM1mS4fNQFoBrABQC2AvxLRCiHEqdiThBBbAWwFgNbWVp7U1AhXuwsIAu5ON9ydbgw9NwQAWPH8irQIwVSom9HoY1oISGxHffgHh9H7QK+qa1vR2iJY5uX22gkm+5EpBL0A6mK2a8P7YukB8IYQwgfgEBEdQEgYdkgcF5MCQojkfQg08L1RNaPJYddRJY2Kz3JkyInxwxPwnfRh5OQIRv46Ej225KdLUPuFWq2HyOQZMquGdgBoJqJGIrIA2Axgm+KcpxGKBkBElQhNFXVJHBOTIt7jXgTGA6r9BptBE7M5VURQrZ+IoMBqQk3J5O8wKICRfYlXGLPrKKMF0oRACOEHcBuAFwDsB/CkEGIvEd1NRJeHT3sBwCAR7QPwMoAvCyHUZSSM5iRbaGRfatdkLlvPEQGgWFgmAO/BJH0INOgCxzBScwRCiOcAPKfY9+2YnwWAL4X/Y7II/4gf5nlm+E764vZrcWMKBoWqM5meIgIglCeImM+VjhPIqU59aRV9MUymk8VMlhKp6PGd8k32IWh3omCF/CfzRGZzFQX6SpjGRjiVIwRBACm0wN5sBxm5koiRDwsBMyXmUjPM680qa2SZqD2Gct9sTknsWoKO2iAe+A8THr3ozDg3VmudVeoYAs4AnAdCDYeKzynm6COPYSFgsgohBLoPjKK5x4DD1UF4LfoqHY2gnOo6eGoCjuUOFJwmN+I69otj6P9dqDdBZAU3ALQ83IKaG2ukXpvJXlgImIxz/OHjGPnbSPRpuG7Ij2/Aju993IWO2qAuXEeV1BTbYDMb4PYFAYTM5wbGvagqkhsFuA66MPynYdV+Np/Lb9h9lMk4A08P4MQjJzD62ij8Q/7o/pqh0D9PPUYEBgOhqVK9sEw2SV1I2Y46r2EhYDJOspvT/KFQXkAPrqOJUPUmSKFb2VyZzvGVyU9YCJiMk+zmVDNkgNFAujGbU5LIakI2yUTX3elGMDxNxeQfnCNgMk7SiGDQgLoyu27M5pQoIwLlamoZmCvMMFWY4B8MTcFZa62wt9jhWOZA0BmEoYSfDfMRFgIm4yhtFPwGgZNlAj1VQTRVale2qjXKiEBZNiuL0357GsxVZtib7TAV8i2AYSFgFAghcOLXJ2BvDj0lWirlL+SyVFqw5KdLYGuy4cnBPty77xCC4QfTm3W2ojgWpfnc0SEnPP6A9Aio/KJyqe/P5B4sBEwc3uNetH+yPbptKjfB0eJA4ZmFWPrgUmnXjThs7vvtkagIAOomLnqiwGrC+lM2dMKDwWKBoAE4POjE0nlFmR4ak2ewEDBxKKtH/EN+jL42isCE2olUBsqGNHosHY0ghMCn/ssIi9sBryk0HXZ01wGYzypH/b/Uw2jXZ26EyT44M8TEkayeXIseBInM5vS4mCyC97gXlnA3ToufUNdvgPHFMRy996i09qAMk4iUIgIi+nai/UKIu9M7HCbTJG1Go4Hr6PE8MJuLJdNW3wwTIdWpodjHNBuASxHqMcDojExGBMo6+sXV+jObiyWTosswsaQkBEKIH8duE9G9CDWVYXRG2QVlMNgMcLW74OpwQfhD3sha3JzUzWj0mx8ApogIuCsZozGzTRY7EOpBzOiM+i/XA18O/Rz0B+E+5IazzQnHaRpEBKpmNPrNDwChxVxF64rQt3sUdu9k5BNYpN/pMCY7STVH8A6ASNsMI4AqAJwf0DkGkwGOZgcczdpMVaj6EOg8Iqi/sx71d9bjQz/9K3o7xlEzSKgZMuCWFiOaNbi+EALe495o0yFXuwumMhMavt2gwdWZbCLViODSmJ/9AE6GexIzTNqIlI4agkDVKcLCXT4cef4IbI02VF9ZneHRyWPxvCLsOzGGkUKBtkVBXGjx4RzJ1xzfPY5d5+9CYCy+LNjeYmchyENSzREclj0QJr+Z8PhRttuLz/3JjupTBFOQMPiLTgwCqLi8QtdC0KRYYayF+Zy11qoSAWDSfM5g5vLVfIL/2kxWcGhgAl4TsGDIAFMwvlJI7175ajtq+UIQMZ9TIvwC7i639Osz2QULAZMVdPaP43hFYhtkvVskq+yoNXAhBbhJDTMJCwGTFXT2T2DCDozZheqY3p9SlRFBxHxONtykhonAXkNMVhCZFz9eHkRRr9pjx9nu1O1CK4fFhAUlNhwbCYldUGhjPudoccBgM0SdZh0tDthb7Cg5t0TqdZnsg4WAyQoiFUMnyoNY2msESo0oXl4QvUFpsbI5kzRVFUaFAAgJo2whWPj5hai7s47tLBgWAibEwDMDGPjjwOSNt8UBW5MNBpP82cNgUKBrIBQR/OF8H57c4MX2H3wQlYVW6dfOFpqqCvC3joHothb9i402djdlQrAQMACA4ZeGceJXJ+L2kZmw+J7F0V4Bsjg+6oY7nAw+VSRQYjfr2mxu/J1x9D/VPzkds9SesW5lDAOwEDBhElWKCJ+Aucos/doqj6GqAl2bzY28MoLDd8cvzVk034TLms3447k+ANpEBAwTgauGGACZdcJUuY7quBkNkKQ884QfpphCoa7+cQihrqBiGBmwEDAIuALwHPEkPKaFE6by6VfPzWiA5HX6gzGLp8fcfvSPJ/6bMEy6YSFg4DromrQUjMGy0AJTofzZw0iiOILeIwJXuyvhftMSW9y2sm0nw8hCqhAQ0UYiaieiDiL62hTnfZSIBBG1yhwPkxjbIhtO/5/T0fSjJsy/cT6KzymONq3XAnWfYv1GBAFXAO7DiRfHlSxXrDBmIWA0QtrjHhEZATwI4EIAPQB2ENE2IcQ+xXlFAL4A4A1ZY2GmxlRiQtUVVar9AZf81a0THj+Ox9TPGw2E+nL9CgGCQPMDzXC2hayfnW1OeI54YFlgQX1dIXBg8tRMVg5F8hN6Ttozk8iM+88C0CGE6AIAInocwCYA+xTn/SuAf0O0HQqTLRjt8uvMlc3q68sdsGiwdiFTGAuMWPjPC+P2BSYC8Bz3YGDiVNx+LVxIgVDOYmLvRFSYIr0J1ry+RreruZl4ZArBQgBHY7Z7AKyPPYGI1gCoE0I8S0RJhYCIbgFwCwDU19dLGCqTKdTNaELRgO+UD652V/zNqcOFtTvX6s4i2VhghGOJA029vrj9WpWQtn2qDaN/H1Xtd7bp19aDiSdj6wiIyADgPwDcMN25QoitALYCQGtrK9fU6QjlzW5xdWie/I0lb8A/qO595O5y6/bmpKyW6hkOmc9ZTXIjM8cyR2IhYBfSvEHmo1UvgLqY7drwvghFAM4AsJ2IugGcDWAbJ4zzC+X0RyQiyEeL5Ij5XISI+Zz067ILad4jUwh2AGgmokYisgDYDGBb5KAQYkQIUSmEaBBCNAB4HcDlQoidEsfEZBnqNQShiCBfb05NSquJPvl5gnwUXSYeaUIQ7ml8G4AXAOwH8KQQYi8R3U1El8u6LpM7BIMCh1RrCPI3IgDUpbNaNKlJJrruQ/rtAcHEIzVHIIR4DsBzin3fTnLuBpljYbKPYyOuqNkcAJTYzSgPm81xRBBCi4jA1mhD6ftLQ30JWiZtv62L8sf9Nd9h07k8RwiRsVrxRAvJImNxLHfAttgWZ4vtWKb/vgQqF1INIgKD2YDVL62Wfh0me2EhyHPeXP4myESTN9rwTbdwTaH0Mk1V6WjMTdDR7MDZHWdLvX42oqwc6uobz6hYM/kBC0EeE3AGor43zr3xUy7vGXuPdCFQRwT69hhqu7ENnmOeeNFd5oBlgSV6o59fbIPDYoTTG1rVPeYJmc9VF9mmemuGmRMsBHmM62Bi8zNrrTUjZnN6dx0dfnkYnsMeDL84HLd/7c61KFobaktpMBAaKwuw99hkXX9X/wQLASMVfS3RZGZEssSrFtbTANDZl0dmc87Urb5VCWPuVsZIhoUgj0lWiqnFyt1xjx8nRvPHbG4mVt+qElJ2IWUkw0KQx7i7E9eJayEEh/rzy2wuqegmqILiiIDRGs4R5DEtv2pB0w+bJi2Rw/8VrpGftFU3o9FvNADMrBUoRwSM1rAQ5DFEBMs8CyzzLCh9X6mm105mLaFX6u6oQ8WlFXFuqs42JwpXqT93Y6XafM7tC8Bmlm8LHkvQF4Sr0xUqL16i7/Ub+Q4LAZMRktlPp4J/zA8EQw11cgVjgRFFa4pQtKZo2nMj5nPHwg17IuZzLfOnf+1cGN8zjpO/ORmNDF2dLiAAzL9xPpY9vEzqtZnMkjvfJEZXqNYQVCeOCMbfGcepl07FPUl7j3nR9KMm1H9Vv70pFlcXRoUACLm0yhYC92E3jt57VLVf77YeDAsBkwESmc0liwiG/ncIXV/pUu3Xu/lcU2UBXjk4EN3WImGc1N9J579rhquGmAygNJsrdUyazSlJ6kKq86dUZYSkRcLY1mgDmdRWFv4hP7wDXunXZzIHCwGjOapEcWVBUi+dfH1KbarUvoTUYDbAviTxYkK9C2++w0LAaI6yK9lUHkO2RhvInH9PqYur1SWkQsjv0qpc5WyuNKPkPSUJF8Mx+oFzBIzmKKc5piodNZgNsC+2xz+REmBrsMF30gdLZeIppVwnU+ZzCz69AJWbKqOmeOZys9TrMdkBC0Ee4hvywdXlgqPFAVOR9v8E1PbTU5eO1ny6BoHxQNSt077EDqNd25p6rSFSm8919sk3n6u4pELq+zPZCQtBHjL8p2Hs27wPAGBZYIn2ICi7qAxV/1Ql/foztZ+u+2KdzOFIJegLwjfgg2W+ZcY9BRZXFca7kA6M45zFfKNm0g8LQR4Sm2j1HvPCe8yLUy+fAoyQLgSJzeb0u2p1Yu8E3jrzLRiLjXGd1grXFE779K2MlJRurQyTLlgI8pCZ+N6kG6XZ3CKdm81FGv8ERgMY2zGGsR1jAIDSC0qnFQJlpKT0Z2KYdKHfbyCTlJk4YaabfGtGMxerb1XbSjafYyTBQpBnCCEy2oegsy95n2I9kjT6SkF0leZzR8PmcwyTblgI8ozAeACl55fC1mADYnKXBocB1lqr9Ot3DuRPVzJgbhGBw2LCwtLJun4RNp9jmHTDOYI8w1RkwsrnVwIAAq4AXB0uONuc8A/7QYaZVbXMhnyLCAxWA8hKEJ74FVmpRl9NVQXoPTXZW7pTA/M5Jv9gIchjjHYjClcUonCFNjfjkNnczEpHc501r66BCAi4j7ijDqqugy5Y61KLvhZXFcaZzylXZcvGN+SLur66D7nR8N2GGZfBMtkPCwGjGb2nXPD4UzObmwoRFPD0eOBsD0Uy1VdVp3OYaYeMBHujHfZGOyo2zmwdgKqEVIOEsRACuz+4GxN7JuAb8MUdW/i5hbBU6XM1dz7DQsBoRtccogHfsA8HPnMArnYXnAecCDpDgmIqN2W9EMwFVQmpBhEBEcF73KsSASCU82Ah0B+cLGY0Q3kTm0lXMmOhEQN/GMD42+NREQD0bz6XqIRUC/O5pK6v7EKqS1gIGM1QewylHhFEzOcSoeebU8R8LsKYx4/+MY/06yYrb40skGP0BQsBoxlqj6GZlY7m482JiDKSJ1DaUUfQs+jmM1KFgIg2ElE7EXUQ0dcSHP8SEe0joj1E9BciWiRzPExmmUtEAOTvzSkTTWqiomsE7M12VFxWgbo76zD/hvnSr81oj7RkMREZATwI4EIAPQB2ENE2IcS+mNN2AWgVQjiJ6DMA/h3A1bLGxGSOcY8fJ0cnpzRMBsKiipmtZI6NCAwFhqiJW+FafZegqhPG8iOCwtWFWLd/HexNdhgsPHGgd2RWDZ0FoEMI0QUARPQ4gE0AokIghHg55vzXAXxC4njynuO/Pj7p69/igLXWqskiMkBtNldf7oDZOLMbTPnF5Vj5p5VwLHPAutCaN/XsqoSxBuZzRpsRBcv0veqbmUSmECwEcDRmuwfA+inO/xSA5xMdIKJbANwCAPX19ekaX97Re38vxndN3kQMDgMcSx1Y9ugy6YvKZtqMJhHWGiusNfJtMNLBqb+dgq/fB0dLqJHOXJ6q1TkCdiFl0ktWrCMgok8AaAXwvkTHhRBbAWwFgNbWVu6eOguEEHAeiJ9LDzqDGH97HKZS+f8MZtKnWA8c+89j6HusL7RhRGhBWYsd9V+rR+l7Smf0XsocQc+wC25fADazvru0Mdohc/KvF0Bsa6na8L44iOiDAL4B4HIhhPy6uDzF0+tBcCKo2m9wGGBdqIHZnKpPsb6nHeLM5gKAq8OFoWeHEHSp/wbTYbcYVeZz3YNsSc2kD5lCsANAMxE1EpEFwGYA22JPIKIzAfwcIRHokziWvCepHfJShzZmc3kUEciw+ubeBIxMpAmBEMIP4DYALwDYD+BJIcReIrqbiC4Pn3YPgEIAvyOit4loW5K3Y+ZIslp7LZrRJDKb07Pr6JTR1yytvpXCqXRxZZi5IHVyWAjxHIDnFPu+HfPzB2Ven5mk+OxiNNzVEHLADP8XnAgmrc1PJ0qzubJZms3lCjKiL+XiO6VvE8PMhaxIFjPyKVpbhKK1kz72Qgh4ej0wmOXXiCtvWnqOBgDAXGlGzc01UdtpX1/IvG0u0Zfyd6a1HTWjb1gI8hQigq3Wpsm1VM1oZmA2l4sUrS5Cy9aW6LZv2AdnuxMGWzpLSEPmc7LXUoiAgLvbHe1JEBG3BZ9ZgHmb50m9NqMdLASMdJQLoBZXzz0iCDgDcB10RW9MznYn6u6sQ9GZ2de9y1xmRsnZJXN6j4j5nNMb6lk8Hjafqy6WK+aHvnUIR354RLW/+KxiFgIdwULASKezTzE1NMeIoP2Wdhz/5XFAsaKk9P2lWSkE6SBiPvdu72h0X0f/uHQhsC/NT3+nfINNRBjppDsiMFeYVSIA6P/mlAnPoWR5jWTlsUxuwkLASCWR2Vx9+dxKVpNVOunZjhpQrzDWRAiSrHtwdbkQ9M58cRyTnbAQMFJRVrfMxmxOSdLuWTp/Ss2E55C5zAxzlVl9IAC4OvUtvPkE5wgYqSifWtNROjrdU6pebZNVU0MauJACQPH6Ynj7vVHXWscyR9RMj9EHLASMVNTWEnMvHTWXm2GuNoPMFL0pRf4PHTtTNyqS7FqZz6344wqp789kHhYCnePp9eCdTe+on+aa7TA65LtXqttTpmcx2dmHz4bRll/umxHzud5ToSmZiPncsvnFGR4Zk+uwEOicif0TGH9rHONvxT+ZF55ZiNZ/tEq/fjr6ECQiG0VACIF9V+2DrdEWFV17ix2WyvTZaTRVFUSFAAiV5rIQMHOFhUDnJPW9maUL5kxIZDanZ9dRT68H/U/1q/Zb6604u/vstKwCXlxViFcODkS32WqCSQf6zKoxUTLpOprIbK4sD83mzOXmtFlBsPkcIwMWAp2T7OakheuoelpIv9EAoI3oKn+H3LaSSQcsBDonaYMUDSICdaJY32ZzWohuotXFQnD3VmZucI5A55z5yplxxmzOdidc7S44muULQb5FBOnuSpaIecVWFFiMmIgxn+sb82CeZM8hRt+wEOgc2yIbbItsKL+oXPNryyodzVaaftSE8Y+NT4pumxPuTndao6+Q+Vwh3ukdie7r7B/PiBCIgIBv0AdLtX7zPvkCCwEjDeXKV60a1oug0KQPs5KiNUUoWhPvfhr0BdM+lqaqgjgh6OqfwLmLK9N6DSWBiQAGnh6IEznXQRestVasP7he6rUZ+bAQMFIYc/vSbjanxD/ux9jOMbjaXXE3qIKVBVjxdHashpXRAU5pPqdFwjjoDWL/J/ar9rsOuRD0BGGwcroxl2EhYKSgXD9QXzF3szklE+9OYPf7d6v2k0nHPhMAFlcrSkg1cCE1l4VsPSJtN6OEzecKTtN3IYDeYRlnpKBKFFemPz+QrxbJmYgIAO5NoGdYCBgpqBLF1el/Yow8parQuUVyY2UBYten9Z4Kmc/JJqn9t84bAuUDLASMFFSuoxIiAiA/b052ixELSibXJgihnoqTQaKIwFhoRNCt3+grX+AcASMFLSICIHRzGnllRLXfc8ST4Gz9sLi6MM58rqt/Astr5JrPlbynBAu/sDDOxdZSY0mbfQaTOVgIdIpv2AeD3ZARl85AArM5GTkCACg5vwT+YT/sLfbJ3gQtDpiK9f1Pu6myAH89MGlwp4X5XPFZxSg+i51O9Yi+vy15zKFvHcKxh47B1mCL3hztLXaUX1wOe6Ncn6FjGprNzb92PuZfO1/Ke6dKJtYtKO062HOImQssBDrF1e4CgoC7yw13lxtDzw8BAJY/tly6EKi7kul7RXHvf/bi8N2H47ql2VvsKFpTBOsCq5RrqttWsgspM3tYCHRKJvsQdKr6FOu7xtzZ5oSv34eR/pG4fEXDXQ1o+E6DlGuqXEj7xiGE4Pl6ZlZw1ZAOCUwE4OlJnCx1LNXCdTS/IoJk9tMyrb4j5nMRJrwB9I3pO0HOyIOFQIc4DySOBqx1VhgLtO9TnLeuoxKtviPmc7FwnoCZLXk1NeQf8SdccWoqNqXkleIb9kH41d7vplJTSp4y3gEvkMA63lxhnjbZKISAb8CX8JilKj4R6x/ywzLfAu8Jb9x+LXoQAPL6FKeK75QPwpfg71RigsGSwt95yAcRSPD6MhMMpvjXByYC8BxNEn1JtvpWms91amA+p8Q/5letIzAWGGF0TP/Akenvo2/QBxFUv95cbgYZp59i8/Z7E+43V07fkU4EQ86tKghp7XGdKlKFgIg2ArgPgBHAL4UQP1IctwJ4FMBaAIMArhZCdKd7HKNuH3z+ILqv2o+xF0+pji/63TIUbyyb9n0ObtgN9x7109+Sv62EfdX0N7v9p++EX+nVAmBZx1qY5039xxdBgXerX1cfIGDF6Dnx+9Y70HJwLQIjfng63PAccMFz0AVLgw2D43KnD5yKKQoZZnPTsfeje3HqpVOq/av+vAplF0z/d9513q6EOZZ1+9ahYHn839nd7Q7961Ys7LXWy4++lFNu+4+PSv/7KvEccuPwlfvhOeiO7qv59wZUfqZm2tcevr4do88MqfbX/2YpSq6omPb1nZe8C+cbY6r9TS+ejoJzpi9zbVv7D/gOq39fLe+cCUvD9Lbe79a9DuFRC8np/ethsE0tRP4hH/Yv2qnabywz4rQjZ035WovJgCJbghX1c0CaEBCREcCDAC4E0ANgBxFtE0LsizntUwCGhRBLiGgzgH8DcHW6x/KFx3bh5fZ+3N5hxaoEH/mLT7yN3TunX6L/3RM2LIL6y33NL9/AkXnTr668b8KOkgSzcRf/5BWMFE7dZYoE8AjUYhMUAmu/9+dprw0DgCMAvjf9qelEhtlcNlFwegHe63wvXJ2uSXvmdheMJfKn4JSR1m/fOILfvnFE+nVjsXqBpjUGfPXgZD7knhfa8efBvdO+9rP7rViX4Pv4laf2YOe7038fv9FjQ3OC7+Mn/2snOv4y/ffxnlN2VCX4Pl72wN8xUDp917df+B0wQ/3kf+6PXoJvmjtrgQt4MMH3ecTln/b7fOnKGjxwzZppxzcTZEYEZwHoEEJ0AQARPQ5gE4BYIdgE4K7wz08BeICISHDvPd2g90QxABgsBhQsL1BFC7LJ1O/2gztNWN1hQs0QoWJMvyKfT8j8Ky4EcDRmuye8L+E5Qgg/gBEAqpiQiG4hop1EtLO/v195mMli1i6afiqGmR1NVQUodaR3iiAVKkcNOOOwkUVAR+REslgIsRXAVgBobW2dcbRQaDOjvMCSdIqi0GpCeQrzuUZD4teX2E0oT+FhMFkCqdRhhnGa11OySJeAckmrdueC0UB4b3MVrjtnUaaHolusJiN+evVq/Oj5Nk1LR98+G7hgl4DZn/jfc4HVhPKC6UXCYkryfbSl9n00JSmwKLaZUV4w/W3CMMX3MZhKcJckH1zmsMA/jT47krzWkML3udCa/tu2TCHoBVAXs10b3pfonB4iMgEoQShpnFa2fOxMAMC+ffswPDysPn79SpRfOH1P37dfehsTfvUKzsc+04rCFdOH6W8+/qa6sQeAP91xLiwpJItffeBV9QEC/vGt90977XzCVGKCuVL9TSRzaoutTGVJXp9CJYnWbGipxoaWas2vO/rxUXR/txtj/xhTJcq/esVS3Hfzgmnfo+1QGwb71F/3e685A5WXT1/9tOf1PRhzqZPFv755DYrXT58s3vn0zoTmhP/vC2fDtmj6ZPGrP381YbL4tX85f1qPL9+wD2/+8k3V/opSE/7xLe1bf5Ks6fjwjf0AgAsQuuHvAHCNEGJvzDmfBbBCCHFrOFn8ESHEVVO9b2trq9i5U51tZxiGYZJDRG8JIVoTHZMWEQgh/ER0G4AXECqwe1gIsZeI7gawUwixDcCvAPyGiDoADAHYLGs8DMMwTGKk5giEEM8BeE6x79sxP7sB/B+ZY2AYhmGmhtP+DMMweQ4LAcMwTJ7DQsAwDJPnsBAwDMPkOdLKR2VBRP0ADs/y5ZUABtI4nFyAP3N+wJ85P5jLZ14khKhKdCDnhGAuENHOZHW0eoU/c37Anzk/kPWZeWqIYRgmz2EhYBiGyXPyTQi2ZnoAGYA/c37Anzk/kPKZ8ypHwDAMw6jJt4iAYRiGUcBCwDAMk+foUgiIaCMRtRNRBxF9LcFxKxE9ET7+BhE1ZGCYaSWFz/wlItpHRHuI6C9ElPMdY6b7zDHnfZSIBBHlfKlhKp+ZiK4K/633EtFvtR5juknh33Y9Eb1MRLvC/74/lIlxpgsiepiI+ojo3STHiYjuD/8+9hDR3BsYCyF09R9CltedAJoAWADsBnCa4px/BvCz8M+bATyR6XFr8JnfD8AR/vkz+fCZw+cVAfgrgNcBtGZ63Br8nZsB7AJQFt6uzvS4NfjMWwF8JvzzaQC6Mz3uOX7m9wJYA+DdJMc/BOB5hHqknQ3gjbleU48RwVkAOoQQXUIIL4DHAWxSnLMJwH+Ff34KwAWUrI9kbjDtZxZCvCyEcIY3X0eoY1wuk8rfGQD+FcC/AXBrOThJpPKZbwbwoBBiGACEEH0ajzHdpPKZBYBIS7ISAMc0HF/aEUL8FaH+LMnYBOBREeJ1AKVEVDOXa+pRCBYCOBqz3RPel/AcIYQfwAiACk1GJ4dUPnMsn0LoiSKXmfYzh0PmOiHEs1oOTCKp/J2XAlhKRH8noteJaKNmo5NDKp/5LgCfIKIehPqffE6boWWMmX7fpyUnmtcz6YOIPgGgFcD7Mj0WmRCRAcB/ALghw0PRGhNC00MbEIr6/kpEK4QQpzI5KMl8DMCvhRA/JqJzEOp6eIYQIpjpgeUKeowIegHUxWzXhvclPCfcW7kEgLqLdu6QymcGEX0QwDcAXC6EUHftzi2m+8xFAM4AsJ2IuhGaS92W4wnjVP7OPQC2CSF8QohDCPUNb9ZofDJI5TN/CsCTACCEeA2ADSFzNr2S0vd9JuhRCHYAaCaiRiKyIJQM3qY4ZxuA68M/XwngJRHOwuQo035mIjoTwM8REoFcnzcGpvnMQogRIUSlEKJBCNGAUF7kciHEzswMNy2k8m/7aYSiARBRJUJTRV0ajjHdpPKZjwC4AACIaDlCQtCv6Si1ZRuA68LVQ2cDGBFCHJ/LG+puakgI4Sei2wC8gFDFwcNCiL1EdDeAnUKIbQB+hVD42IFQUmZz5kY8d1L8zPcAKATwu3Be/IgQ4vKMDXqOpPiZdUWKn/kFABcR0T4AAQBfFkLkbLSb4me+A8AviOh2hBLHN+Tygx0RPYaQmFeG8x7fAWAGACHEzxDKg3wIQAcAJ4Ab53zNHP59MQzDMGlAj1NDDMMwzAxgIWAYhslzWAgYhmHyHBYChmGYPIeFgGEYJs9hIWAYhslzWAgYhmHyHBYChpkjRLQu7AtvI6KCcB+AMzI9LoZJFV5QxjBpgIi+h5C1gR1AjxDihxkeEsOkDAsBw6SBsA/ODoT6HpwrhAhkeEgMkzI8NcQw6aECIS+nIoQiA4bJGTgiYJg0QETbEOqe1QigRghxW4aHxDApozv3UYbRGiK6DoBPCPFbIjICeJWIPiCEeCnTY2OYVOCIgGEYJs/hHAHDMEyew0LAMAyT57AQMAzD5DksBAzDMHkOCwHDMEyew0LAMAyT57AQMAzD5Dn/H7px8k79oT3KAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# NBVAL_IGNORE_OUTPUT\n", + "L = 1.0\n", + "I = lambda y: 0 if abs(y-L/2.0) > 0.1 else 1\n", + "\n", + "Nx = 10\n", + "c = 0.5\n", + "C = 1\n", + "dt = C*(L/Nx)/c\n", + "nperiods = 4\n", + "T = L/c*nperiods # One period: c*T = L\n", + "Nt = int(round(T/dt))\n", + "x = np.linspace(0, L, Nx+1) # Mesh points in space\n", + "\n", + "L = 1.0\n", + "c = 0.5\n", + "dt = (L/10)/c # Nx=10\n", + "I = lambda x: 0 if abs(x-L/2.0) > 0.1 else 1\n", + "# I = lambda x: 5*x\n", + "U_0 = lambda t: t\n", + "T = 4\n", + "\n", + "\n", + "u_mine, x_arr, t_arr, cpu_time = solver(I, None, None, c, L, dt, C, T)\n", + "u, x, t, cpu = python_solver(I, None, None, c, L, dt, C, T)\n", + "\n", + "plt.plot(x, u, label='NumPy', linewidth=4)\n", + "plt.plot(x_arr, u_mine, 'm:', label='Devito', linewidth=6)\n", + "plt.xlabel('x')\n", + "plt.ylabel('u')\n", + "plt.legend(loc='best')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The program [`wave1D_n0.py`](src-wave/wave1D/wave1D_n0.py)\n", + "also contains this solver function.\n", "\n", "It would be nice to modify the `test_quadratic` test case from the\n", - "`wave1D_u0.py` with Dirichlet conditions, described in the section [wave:pde1:impl:vec:verify:quadratic](#wave:pde1:impl:vec:verify:quadratic). However, the Neumann\n", + "`wave1D_u0.py` with Dirichlet conditions, described in the section on [verification](wave1D_prog.ipynb#wave:pde1:impl:vec:verify:quadratic). However, the Neumann\n", "conditions require the polynomial variation in the $x$ direction to\n", "be of third degree, which causes challenging problems when\n", "designing a test where the numerical solution is known exactly.\n", @@ -309,204 +388,386 @@ "with a plug wave at rest and see that the initial condition is\n", "reached again perfectly after one period of motion, but such\n", "a test requires $C=1$ (so the numerical solution coincides with\n", - "the exact solution of the PDE, see the section [Numerical dispersion relation](wave_analysis.ipynb#wave:pde1:num:dispersion)).\n", + "the exact solution of the PDE, see the section [Numerical dispersion relation](wave_analysis.ipynb#wave:pde1:num:dispersion)).\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Notice.**\n", "\n", + "The program [`wave1D_dn.py`](src-wave/wave1D/wave1D_dn.py)\n", + "solves the 1D wave equation $u_{tt}=c^2u_{xx}+f(x,t)$ with\n", + "quite general boundary and initial conditions:\n", "\n", - "## Index set notation\n", - "
\n", + " * $x=0$: $u=U_0(t)$ or $u_x=0$\n", "\n", + " * $x=L$: $u=U_L(t)$ or $u_x=0$\n", "\n", - "To improve our mathematical writing and our implementations,\n", - "it is wise to introduce a special notation for index sets. This means\n", - "that we write\n", - "$x_i$, followed by $i\\in\\mathcal{I}_x$, instead of $i=0,\\ldots,N_x$.\n", - "Obviously, $\\mathcal{I}_x$ must be the index set $\\mathcal{I}_x =\\{0,\\ldots,N_x\\}$, but it\n", - "is often advantageous to have a symbol for this set rather than\n", - "specifying all its elements (all the time, as we have done up to\n", - "now). This new notation saves writing and makes\n", - "specifications of algorithms and their implementation as computer code\n", - "simpler.\n", + " * $t=0$: $u=I(x)$\n", "\n", - "The first index in the set will be denoted $\\mathcal{I}_x^0$\n", - "and the last $\\mathcal{I}_x^{-1}$. When we need to skip the first element of\n", - "the set, we use $\\mathcal{I}_x^{+}$ for the remaining subset\n", - "$\\mathcal{I}_x^{+}=\\{1,\\ldots,N_x\\}$. Similarly, if the last element is\n", - "to be dropped, we write $\\mathcal{I}_x^{-}=\\{0,\\ldots,N_x-1\\}$ for the\n", - "remaining indices.\n", - "All the\n", - "indices corresponding to inner grid points are specified by\n", - "$\\mathcal{I}_x^i=\\{1,\\ldots,N_x-1\\}$. For the time domain we find it\n", - "natural to explicitly use 0 as the first index, so we will usually\n", - "write $n=0$ and $t_0$ rather than $n=\\mathcal{I}_t^0$. We also avoid notation\n", - "like $x_{\\mathcal{I}_x^{-1}}$ and will instead use $x_i$, $i={\\mathcal{I}_x^{-1}}$.\n", + " * $t=0$: $u_t=V(x)$\n", "\n", - "The Python code associated with index sets applies the following\n", - "conventions:\n", + "A lot of test examples are also included in the program:\n", "\n", + " * A rectangular plug-shaped initial condition. (For $C=1$ the solution\n", + " will be a rectangle that jumps one cell per time step, making the case\n", + " well suited for verification.)\n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "
Notation Python
$\\mathcal{I}_x$ Ix
$\\mathcal{I}_x^0$ Ix[0]
$\\mathcal{I}_x^{-1}$ Ix[-1]
$\\mathcal{I}_x^{-}$ Ix[:-1]
$\\mathcal{I}_x^{+}$ Ix[1:]
$\\mathcal{I}_x^i$ Ix[1:-1]
\n", - "**Why index sets are useful.**\n", + " * A Gaussian function as initial condition.\n", "\n", - "An important feature of the index set notation is that it\n", - "keeps our formulas and code independent of how\n", - "we count mesh points. For example, the notation $i\\in\\mathcal{I}_x$ or $i=\\mathcal{I}_x^0$\n", - "remains the same whether $\\mathcal{I}_x$ is defined as above or as starting at 1,\n", - "i.e., $\\mathcal{I}_x=\\{1,\\ldots,Q\\}$. Similarly, we can in the code define\n", - "`Ix=range(Nx+1)` or `Ix=range(1,Q)`, and expressions\n", - "like `Ix[0]` and `Ix[1:-1]` remain correct. One application where\n", - "the index set notation is convenient is\n", - "conversion of code from a language where arrays has base index 0 (e.g.,\n", - "Python and C) to languages where the base index is 1 (e.g., MATLAB and\n", - "Fortran). Another important application is implementation of\n", - "Neumann conditions via ghost points (see next section).\n", + " * A triangular profile as initial condition, which resembles the\n", + " typical initial shape of a guitar string.\n", "\n", + " * A sinusoidal variation of $u$ at $x=0$ and either $u=0$ or\n", + " $u_x=0$ at $x=L$.\n", "\n", + " * An analytical solution $u(x,t)=\\cos(m\\pi t/L)\\sin({\\frac{1}{2}}m\\pi x/L)$, which can be used for convergence rate tests.\n", "\n", - "For the current problem setting in the $x,t$ plane, we work with\n", - "the index sets" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "
\n", - "\n", - "$$\n", - "\\begin{equation}\n", - "\\mathcal{I}_x = \\{0,\\ldots,N_x\\},\\quad \\mathcal{I}_t = \\{0,\\ldots,N_t\\},\n", - "\\label{_auto3} \\tag{7}\n", - "\\end{equation}\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "defined in Python as" + "\n" ] }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ - "mathcal{I}_x = range(0, Nx+1)\n", - "mathcal{I}_t = range(0, Nt+1)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "A finite difference scheme can with the index set notation be specified as" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "$$\n", - "\\begin{align*}\n", - "u_i^{n+1} &= u^n_i - \\frac{1}{2}\n", - "C^2\\left(u^{n}_{i+1}-2u^{n}_{i} + u^{n}_{i-1}\\right),\\quad,\n", - "i\\in\\mathcal{I}_x^i,\\ n=0,\\\\ \n", - "u^{n+1}_i &= -u^{n-1}_i + 2u^n_i + C^2\n", - "\\left(u^{n}_{i+1}-2u^{n}_{i}+u^{n}_{i-1}\\right),\n", - "\\quad i\\in\\mathcal{I}_x^i,\\ n\\in\\mathcal{I}_t^i,\\\\ \n", - "u_i^{n+1} &= 0,\n", - "\\quad i=\\mathcal{I}_x^0,\\ n\\in\\mathcal{I}_t^{-},\\\\ \n", - "u_i^{n+1} &= 0,\n", - "\\quad i=\\mathcal{I}_x^{-1},\\ n\\in\\mathcal{I}_t^{-}\\thinspace .\n", - "\\end{align*}\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The corresponding implementation becomes" + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import time\n", + "from devito import Constant, Grid, TimeFunction, SparseTimeFunction, Function, Eq, solve, Operator, Buffer" ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ - "# Initial condition\n", - "for i in mathcal{I}_x[1:-1]:\n", - " u[i] = u_n[i] - 0.5*C2*(u_n[i-1] - 2*u_n[i] + u_n[i+1])\n", + "# %load -s python_solver, src-wave/wave1D/wave1D_dn.py\n", + "def python_solver(I, V, f, c, U_0, U_L, L, dt, C, T,\n", + " user_action=None, version='scalar'):\n", + " \"\"\"\n", + " Solve u_tt=c^2*u_xx + f on (0,L)x(0,T].\n", + " u(0,t)=U_0(t) or du/dn=0 (U_0=None), u(L,t)=U_L(t) or du/dn=0 (u_L=None).\n", + " \"\"\"\n", + " Nt = int(round(T/dt))\n", + " t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time\n", + " dx = dt*c/float(C)\n", + " Nx = int(round(L/dx))\n", + " x = np.linspace(0, L, Nx+1) # Mesh points in space\n", + " C2 = C**2; dt2 = dt*dt # Help variables in the scheme\n", + " # Make sure dx and dt are compatible with x and t\n", + " dx = x[1] - x[0]\n", + " dt = t[1] - t[0]\n", "\n", - "# Time loop\n", - "for n in mathcal{I}_t[1:-1]:\n", - " # Compute internal points\n", - " for i in mathcal{I}_x[1:-1]:\n", - " u[i] = - u_nm1[i] + 2*u_n[i] + \\\n", - " C2*(u_n[i-1] - 2*u_n[i] + u_n[i+1])\n", - " # Compute boundary conditions\n", - " i = mathcal{I}_x[0]; u[i] = 0\n", - " i = mathcal{I}_x[-1]; u[i] = 0" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "**Notice.**\n", + " # Wrap user-given f, I, V, U_0, U_L if None or 0\n", + " if f is None or f == 0:\n", + " f = (lambda x, t: 0) if version == 'scalar' else \\\n", + " lambda x, t: np.zeros(x.shape)\n", + " if I is None or I == 0:\n", + " I = (lambda x: 0) if version == 'scalar' else \\\n", + " lambda x: np.zeros(x.shape)\n", + " if V is None or V == 0:\n", + " V = (lambda x: 0) if version == 'scalar' else \\\n", + " lambda x: np.zeros(x.shape)\n", + " if U_0 is not None:\n", + " if isinstance(U_0, (float,int)) and U_0 == 0:\n", + " U_0 = lambda t: 0\n", + " # else: U_0(t) is a function\n", + " if U_L is not None:\n", + " if isinstance(U_L, (float,int)) and U_L == 0:\n", + " U_L = lambda t: 0\n", + " # else: U_L(t) is a function\n", "\n", - "The program [`wave1D_dn.py`](src-wave/wave1D/python/wave1D_dn.py)\n", - "applies the index set notation and\n", - "solves the 1D wave equation $u_{tt}=c^2u_{xx}+f(x,t)$ with\n", - "quite general boundary and initial conditions:\n", + " u = np.zeros(Nx+1) # Solution array at new time level\n", + " u_n = np.zeros(Nx+1) # Solution at 1 time level back\n", + " u_nm1 = np.zeros(Nx+1) # Solution at 2 time levels back\n", "\n", - " * $x=0$: $u=U_0(t)$ or $u_x=0$\n", + " Ix = list(range(0, Nx+1))\n", + " It = list(range(0, Nt+1))\n", "\n", - " * $x=L$: $u=U_L(t)$ or $u_x=0$\n", + " import time; t0 = time.perf_counter() # CPU time measurement\n", "\n", - " * $t=0$: $u=I(x)$\n", + " # Load initial condition into u_n\n", + " for i in Ix:\n", + " u_n[i] = I(x[i])\n", "\n", - " * $t=0$: $u_t=V(x)$\n", + " if user_action is not None:\n", + " user_action(u_n, x, t, 0)\n", "\n", - "The program combines Dirichlet and Neumann conditions, scalar and vectorized\n", - "implementation of schemes, and the index set notation into one piece of code.\n", - "A lot of test examples are also included in the program:\n", + " # Special formula for the first step\n", + " for i in Ix[1:-1]:\n", + " u[i] = u_n[i] + dt*V(x[i]) + \\\n", + " 0.5*C2*(u_n[i-1] - 2*u_n[i] + u_n[i+1]) + \\\n", + " 0.5*dt2*f(x[i], t[0])\n", "\n", - " * A rectangular plug-shaped initial condition. (For $C=1$ the solution\n", - " will be a rectangle that jumps one cell per time step, making the case\n", - " well suited for verification.)\n", + " i = Ix[0]\n", + " if U_0 is None:\n", + " # Set boundary values du/dn = 0\n", + " # x=0: i-1 -> i+1 since u[i-1]=u[i+1]\n", + " # x=L: i+1 -> i-1 since u[i+1]=u[i-1])\n", + " ip1 = i+1\n", + " im1 = ip1 # i-1 -> i+1\n", + " u[i] = u_n[i] + dt*V(x[i]) + \\\n", + " 0.5*C2*(u_n[im1] - 2*u_n[i] + u_n[ip1]) + \\\n", + " 0.5*dt2*f(x[i], t[0])\n", + " else:\n", + " u[0] = U_0(dt)\n", "\n", - " * A Gaussian function as initial condition.\n", + " i = Ix[-1]\n", + " if U_L is None:\n", + " im1 = i-1\n", + " ip1 = im1 # i+1 -> i-1\n", + " u[i] = u_n[i] + dt*V(x[i]) + \\\n", + " 0.5*C2*(u_n[im1] - 2*u_n[i] + u_n[ip1]) + \\\n", + " 0.5*dt2*f(x[i], t[0])\n", + " else:\n", + " u[i] = U_L(dt)\n", "\n", - " * A triangular profile as initial condition, which resembles the\n", - " typical initial shape of a guitar string.\n", + " if user_action is not None:\n", + " user_action(u, x, t, 1)\n", "\n", - " * A sinusoidal variation of $u$ at $x=0$ and either $u=0$ or\n", - " $u_x=0$ at $x=L$.\n", + " # Update data structures for next step\n", + " #u_nm1[:] = u_n; u_n[:] = u # safe, but slower\n", + " u_nm1, u_n, u = u_n, u, u_nm1\n", "\n", - " * An analytical solution $u(x,t)=\\cos(m\\pi t/L)\\sin({\\frac{1}{2}}m\\pi x/L)$, which can be used for convergence rate tests.\n", + " for n in It[1:-1]:\n", + " # Update all inner points\n", + " if version == 'scalar':\n", + " for i in Ix[1:-1]:\n", + " u[i] = - u_nm1[i] + 2*u_n[i] + \\\n", + " C2*(u_n[i-1] - 2*u_n[i] + u_n[i+1]) + \\\n", + " dt2*f(x[i], t[n])\n", "\n", + " elif version == 'vectorized':\n", + " u[1:-1] = - u_nm1[1:-1] + 2*u_n[1:-1] + \\\n", + " C2*(u_n[0:-2] - 2*u_n[1:-1] + u_n[2:]) + \\\n", + " dt2*f(x[1:-1], t[n])\n", + " else:\n", + " raise ValueError('version=%s' % version)\n", "\n", + " # Insert boundary conditions\n", + " i = Ix[0]\n", + " if U_0 is None:\n", + " # Set boundary values\n", + " # x=0: i-1 -> i+1 since u[i-1]=u[i+1] when du/dn=0\n", + " # x=L: i+1 -> i-1 since u[i+1]=u[i-1] when du/dn=0\n", + " ip1 = i+1\n", + " im1 = ip1\n", + " u[i] = - u_nm1[i] + 2*u_n[i] + \\\n", + " C2*(u_n[im1] - 2*u_n[i] + u_n[ip1]) + \\\n", + " dt2*f(x[i], t[n])\n", + " else:\n", + " u[0] = U_0(t[n+1])\n", + "\n", + " i = Ix[-1]\n", + " if U_L is None:\n", + " im1 = i-1\n", + " ip1 = im1\n", + " u[i] = - u_nm1[i] + 2*u_n[i] + \\\n", + " C2*(u_n[im1] - 2*u_n[i] + u_n[ip1]) + \\\n", + " dt2*f(x[i], t[n])\n", + " else:\n", + " u[i] = U_L(t[n+1])\n", + "\n", + " if user_action is not None:\n", + " if user_action(u, x, t, n+1):\n", + " break\n", + "\n", + " # Update data structures for next step\n", + " #u_nm1[:] = u_n; u_n[:] = u # safe, but slower\n", + " u_nm1, u_n, u = u_n, u, u_nm1\n", + "\n", + " # Important to correct the mathematically wrong u=u_nm1 above\n", + " # before returning u\n", + " u = u_n\n", + " cpu_time = time.perf_counter() - t0\n", + " return u, x, t, cpu_time\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "# %load -s solver, src-wave/wave1D/wave1D_dn.py\n", + "def devito_solver(I, V, f, c, U_0, U_L, L, dt, C, T,\n", + " user_action=None, version='scalar'):\n", + " \"\"\"\n", + " Solve u_tt=c^2*u_xx + f on (0,L)x(0,T].\n", + " u(0,t)=U_0(t) or du/dn=0 (U_0=None), u(L,t)=U_L(t) or du/dn=0 (u_L=None).\n", + " \"\"\"\n", + " Nt = int(round(T/dt))\n", + " t = np.linspace(0, Nt*dt, Nt+1) # Mesh points in time\n", + " dx = dt*c/float(C)\n", + " Nx = int(round(L/dx))\n", + " x = np.linspace(0, L, Nx+1) # Mesh points in space\n", + " C2 = C**2 # Help variable in the scheme\n", + "\n", + " # Make sure dx and dt are compatible with x and t\n", + " dx = x[1] - x[0]\n", + " dt = t[1] - t[0]\n", + "\n", + " # Wrap user-given f, I, V, U_0, U_L if None or 0\n", + " if f is None or f == 0:\n", + " f = (lambda x, t: 0) if version == 'scalar' else \\\n", + " lambda x, t: np.zeros(x.shape)\n", + " if I is None or I == 0:\n", + " I = (lambda x: 0) if version == 'scalar' else \\\n", + " lambda x: np.zeros(x.shape)\n", + " if V is None or V == 0:\n", + " V = (lambda x: 0) if version == 'scalar' else \\\n", + " lambda x: np.zeros(x.shape)\n", + " if U_0 is not None:\n", + " if isinstance(U_0, (float,int)) and U_0 == 0:\n", + " U_0 = lambda t: 0\n", + " # else: U_0(t) is a function\n", + " else:\n", + " U_0 = lambda t: 0\n", + " if U_L is not None:\n", + " if isinstance(U_L, (float,int)) and U_L == 0:\n", + " U_L = lambda t: 0\n", + " else:\n", + " U_L = lambda t: 0\n", + " # else: U_L(t) is a function\n", + "\n", + " t0 = time.perf_counter() # Measure CPU time\n", + " \n", + " # Set up grid\n", + " grid = Grid(shape=(Nx+1), extent=(L))\n", + " t_s = grid.stepping_dim\n", + " \n", + " # Create and initialise u\n", + " u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2)\n", + " print(u.data.shape)\n", + " u.data[:] = I(x)\n", + "# for i in range(Nx+1):\n", + "# u.data[:,i] = I(x[i])\n", + "\n", + " x_dim = grid.dimensions[0]\n", + " t_dim = grid.time_dim\n", + " \n", + " # The wave equation we are trying to solve\n", + " pde = (1/c**2)*u.dt2-u.dx2\n", + " \n", + " # Source term and injection into equation\n", + " dt_symbolic = grid.time_dim.spacing \n", + " src = SparseTimeFunction(name='f', grid=grid, npoint=Nx+1, nt=Nt+1)\n", + "\n", + " for i in range(Nt):\n", + " src.data[i] = f(x, t[i])\n", + "\n", + " src.coordinates.data[:, 0] = x\n", + " src_term = src.inject(field=u.forward, expr=src * (dt_symbolic**2))\n", + " stencil = Eq(u.forward, solve(pde, u.forward))\n", + " \n", + " # Set up special stencil for initial timestep with substitution for u.backward\n", + " v = Function(name='v', grid=grid, npoint=Nx+1, nt=1)\n", + " v.data[:] = V(x[:])\n", + " stencil_init = stencil.subs(u.backward, u.forward - dt_symbolic*v)\n", + " \n", + " import inspect\n", + "# print(inspect.getsource(U_0))\n", + " \n", + " u0 = Function(name='u0', grid=grid, npoint=1)\n", + "# print(u0.data.shape)\n", + " # Boundary conditions, depending on arguments\n", + " bc = []\n", + "# if U_0 is None and U_L is None:\n", + " bc += [Eq(u[t_s+1, 0], u[t_s+1, 1])]\n", + " bc += [Eq(u[t_s+1, Nx], u[t_s+1, Nx-1])]\n", + "# else:\n", + "# if U_0 is not None:\n", + "# bc += [Eq(u[t_s+1, 0], U_0(t_s+1))]\n", + "# if U_L is not None:\n", + "# bc += [Eq(u[t_s+1, L], U_L(t_s+1))]\n", + "\n", + "# bc = [Eq(u[t_s+1, 0], -u[t_s+1, 1])]\n", + "\n", + " \n", + "# print(bc)\n", + " \n", + " op_init = Operator([stencil_init]+src_term+bc)\n", + " op = Operator([stencil]+src_term+bc)\n", + " \n", + " op_init.apply(time_M=1, dt=dt)\n", + " op.apply(time_m=1,time_M=Nt, dt=dt)\n", + "\n", + " cpu_time = time.perf_counter() - t0\n", + " return u.data[1], x, t, cpu_time\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(3, 11)\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Data type float64 of runtime value `dt` does not match the Constant data type \n", + "Operator `Kernel` run in 0.01 s\n", + "Data type float64 of runtime value `dt` does not match the Constant data type \n", + "Operator `Kernel` run in 0.01 s\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbpUlEQVR4nO3dfXRV9Z3v8fcXEhKBABKCIA8GUDQWFZ340LGt1Kepjrd2er3cPsmotaw147jUwdsH72rr7Uzvna52tNbOjMMMY6cuZ/RWWZY71qkg+EAVNCCKJoioKLEoCSAQ5CEP3/vHOcacnRM4kLPPPnvvz2stl+dhc/LdEPhk//Y+n2PujoiIpNeQqAcQEZFoKQhERFJOQSAiknIKAhGRlFMQiIikXEXUAxypcePGeX19fdRjiIjEypo1a9rdvS7fc7ELgvr6epqamqIeQ0QkVszs7YGe09KQiEjKKQhERFJOQSAiknKxO0eQT2dnJ62trezfvz/qUcpOdXU1kydPprKyMupRRKRMJSIIWltbqampob6+HjOLepyy4e5s376d1tZWpk2bFvU4IlKmQlsaMrNqM3vezF4ys1fN7H/l2abKzB40s01mttrM6o/ma+3fv5/a2lqFQICZUVtbqyMlETmkMI8IDgAXunuHmVUCK83sMXdf1WebrwM73f1EM/sS8CPgvx/NF1MI5Kffl2Tbumsff7diExvf64h6FCmRE48byf/+k9OK+pqhBYFn+q0/+u6szP4X7Ly+Erg9e/sh4OdmZq5ubJGC3PLgOla9uSPqMaSEDnb3FP01Q71qyMyGmtk6YBuw1N1XBzaZBGwBcPcuYBdQm+d15ptZk5k1tbW1hTnyUTMzFixY0Hv/Jz/5CbfffntRXvv2229n0qRJzJ49m1mzZrFkyZKivK7E2+b2vQoBKYpQg8Ddu919NjAZOMfMZh3l6yx090Z3b6yry/sO6chVVVWxePFi2tvbQ3n9W265hXXr1vGrX/2K6667jp6e4v9UIPGyrOX9qEeQhCjJVUPu/oGZrQA+B7zS56l3gSlAq5lVAKOB7Uf7deq//eig5jyczX/zxwM+V1FRwfz587nzzjv54Q9/mPPcNddcwxVXXMFVV10FwMiRI+no6ODJJ5/k+9//PmPGjGH9+vXMnTuX0047jbvuuot9+/bxyCOPMGPGjJzXamhooKKigi1btjBnzhw2btxIZWUlu3fv5owzzui9L8n3eHNuEMz/zHQuOmV8RNNIqYyoKv4/26EFgZnVAZ3ZEDgGuITMyeC+lgB/CjwHXAUsj/P5gRtuuIHTTz+db37zmwX/mpdeeomWlhbGjh3L9OnTuf7663n++ee56667uPvuu/npT3+as/3q1asZMmQIU6dOZc6cOTz66KN84Qtf4IEHHuCLX/yiQiAlduw9SNPm3GWhq887gSljh0c0kcRZmEtDE4EVZvYy8AKZcwT/YWY/MLPPZ7dZBNSa2SbgL4FvhzhP6EaNGsW8efP42c9+VvCvOfvss5k4cSJVVVXMmDGDSy+9FIDTTjuNzZs392535513Mnv2bG699VYefPBBzIzrr7+ee++9F4B7772Xa6+9tqj7I+Vr+YZt9PT5kemUCTUKATlqYV419DJwZp7Hv9fn9n7gv4U1QxRuvvlmzjrrrJx/lCsqKnrX9Ht6ejh48GDvc1VVVb23hwwZ0nt/yJAhdHV19T53yy23cOutt+Z8rfPPP5/Nmzfz5JNP0t3dzaxZR3UKRmJoafN7OfcvPfW4iCaRJEjEO4s/cqg1/FIZO3Ysc+fOZdGiRVx33XVApjp7zZo1zJ07lyVLltDZ2Vm0rzdv3jy+8pWv8N3vfrdorynlbX9nN09vzL0o4ZJTJ0Q0jSSBSudCsGDBgpyrh77xjW/w1FNPccYZZ/Dcc88xYsSIon2tr371q+zcuZMvf/nLRXtNKW/PvtHOvs7u3vsTRlUza9KoCCeSuLO4nZttbGz04AfTtLS00NDQENFE0XrooYf49a9/zX333TfgNmn+/Umi7yx+mX9/fkvv/a+dN5W//kJx32kqyWNma9y9Md9ziVoaSpsbb7yRxx57jN/85jdRjyIl0tPjLGvZlvOYloVksBQEMXb33XdHPYKU2LrWD2jbc6D3/siqCs6bPjbCiSQJEnOOIG5LXKWi35dkWRp4E9kFJ9dRVTE0omkkKRIRBNXV1Wzfvl3/6AV89HkE1dXVUY8iRbIsEAS6bFSKIRFLQ5MnT6a1tZVyLaSL0kefUCbxt7l9L69v+7huumKIMWemKiVk8BIRBJWVlfoELkm84LLQOdPGMnq4KkVk8BKxNCSSBsEguETLQlIkCgKRGNix9yBNb+eWzCkIpFgUBCIx8ETL+zklcw0TRzH5WJXMSXEoCERiQMtCEiYFgUiZ29/ZzTOv55bM6bJRKSYFgUiZ+92m3JK5iaOr+cTxKpmT4lEQiJS54LLQxQ3HYWYRTSNJpCAQKWP5S+a0LCTFpSAQKWMvbvmA9o6PS+Zqqio4b3pthBNJEikIRMpYvpK5YRX6ayvFpe8okTK2rEWXjUr4FAQiZeqt9r1sCpbMnaySOSk+BYFImVra/F7O/XOnj2X0MSqZk+JTEIiUqX7vJm7QspCEQ0EgUoa2dxxgzds7cx67WOcHJCQKApEy9MSGbTklc6eqZE5CpCAQKUPBj6TU0YCESUEgUmZUMielpiAQKTMrX88tmTteJXMSMgWBSJnpVzJ3qkrmJFwKApEy0t3jPLFB7yaW0lIQiJSRdVt20t5xsPd+TVUF505TyZyEK7QgMLMpZrbCzJrN7FUzuynPNqPN7P+Z2UvZba4Nax6ROFjanFs5rZI5KYWKEF+7C1jg7mvNrAZYY2ZL3b25zzY3AM3u/l/MrA54zczud/eDeV9RJOGCtRJaFpJSCO1HDXff6u5rs7f3AC3ApOBmQI1lzoSNBHaQCRCR1HmzrYM32vb23lfJnJRKSY45zaweOBNYHXjq50AD8HtgPXCTu/fk+fXzzazJzJra2trCHlckEsGrhc6bXquSOSmJ0IPAzEYCDwM3u/vuwNN/BKwDjgdmAz83s34XTLv7QndvdPfGurq6kCcWiUa/kjktC0mJhBoEZlZJJgTud/fFeTa5FljsGZuAt4BTwpxJpBy1dxxgzTsqmZNohHnVkAGLgBZ3v2OAzd4BLspufxxwMvBmWDOJlKvlG7bhgZK5SWOOiW4gSZUwrxo6H7gaWG9m67KP3QZMBXD3e4C/An5hZusBA77l7u15Xksk0bQsJFEKLQjcfSWZf9wPtc3vgUvDmkEkDvYd7OaZ13MvglAQSCnpnSoiEVu5qZ39nR9fLDdpzDEqmZOSUhCIRCz4JrKLG8arZE5KSkEgEqHuHueJltxaiUtOnRDRNJJWCgKRCK3bspPte3NL5s6ZNjbCiSSNFAQiEXo8cLXQnFPGq2ROSk7fcSIR0mWjUg4UBCIReaOtgzf7lMxVDjXmnKwKFSk9BYFIRPKVzI2qVsmclJ6CQCQiWhaScqEgEIlA254DrA2UzF3UoCCQaCgIRCKwIlAy94njVTIn0VEQiEQgeNmoloUkSgoCkRLbd7CblZtUMiflQ0EgUmLPvN7Wr2Tu1IkqmZPoKAhESizf1UIqmZMoKQhESqi7x1m+Ibdk7mJdLSQRUxCIlNCL7wRK5qorOHe6SuYkWgoCkRIKLgt99uTxVA7VX0OJlr4DRUpI7yaWcqQgECmRTds6eLM9t2TuApXMSRlQEIiUiErmpFwpCERKJPjZxFoWknKhIBApgbY9B3hxywc5j+myUSkXCgKREli+4f2ckrlZk0ZxvErmpEwoCERKoN/VQg0TIppEpD8FgUjIPjzYxTOvt+c8pvMDUk4UBCIhe+b1dg505ZbMNUysiXAikVwKApGQqWROyp2CQCRE+UrmtCwk5UZBIBKite/sZEefkrlR1RWcM00lc1JeQgsCM5tiZivMrNnMXjWzmwbYbo6Zrctu81RY84hEoV/J3CkqmZPyUxHia3cBC9x9rZnVAGvMbKm7N3+0gZmNAf4e+Jy7v2Nm40OcR6Sk3F0lcxILof1o4u5b3X1t9vYeoAWYFNjsK8Bid38nu902RBLijbYO3gqWzM1UyZyUn5Ico5pZPXAmsDrw1EzgWDN70szWmNm8AX79fDNrMrOmtra2fJuIlJ3H85TM1ahkTspQ6EFgZiOBh4Gb3X134OkK4A+APwb+CPiumc0Mvoa7L3T3RndvrKvTT1QSD8sCQXCploWkTIV5jgAzqyQTAve7++I8m7QC2919L7DXzJ4GzgA2hjmXSNjylswpCKRMhXnVkAGLgBZ3v2OAzX4NfMrMKsxsOHAumXMJIrH2REtuydxpk0YzcbRK5qQ8hXlEcD5wNbDezNZlH7sNmArg7ve4e4uZ/SfwMtAD/LO7vxLiTCIloauFJE5CCwJ3Xwkc9n307v5j4MdhzSFSah8e7GLlptySOX32gJQzvbNFpMie3qiSOYkXBYFIkS1rUcmcxIuCQKSI8pXM6bJRKXcKApEiWvN2/5K5s1UyJ2VOQSBSREub38u5f6FK5iQG9B0qUiT5Sub0JjKJAwWBSJFs2tbB5u0f9t5XyZzEhYJApEiCJXOfnDFOJXMSCwoCkSLJd9moSBwoCESKYNue/awLlMxdoncTS0woCESK4ImWbTklc6dPHs2E0dXRDSRyBBQEIkXQ72ohHQ1IjCgIRAZp74H+JXM6PyBxUlD7qJl9L9/j7v6D4o4jEj/PvN7GwT4lc5OPPYZTJqhkTuKj0BrqvX1uVwNXoA+QEQFgaXNut5BK5iRuCgoCd//bvvfN7CfAb0OZSCRGurp7WL5Bl41KvB3tOYLhwORiDiISR2ve3snODzt7748+ppKz61UyJ/FS6DmC9cBHF8cNBeoAnR+Q1AteLfTZk+tUMiexU+g5giv63O4C3nf3rhDmEYkNd2dpv3cTT4hoGpGjV+g5grfDHkQkbl7f1sHbfUrmhg0dwgUnq2RO4kfHsCJHKbgs9MkZtYysKvQgW6R8KAhEjlIwCHS1kMSVgkDkKGzb3b9kTrUSElcKApGjsKwl901kKpmTOFMQiByF4GcTq3Ja4kxBIHKE9h7o4ndvbM957JJPKAgkvhQEIkcoWDI3ZewxnHycSuYkvhQEIkco+NnElzRMUMmcxJqCQOQIZErm+reNisSZgkDkCDS9vZMP+pXMHRvhRCKDpyAQOQLBN5FdeMp4KlQyJzEX2newmU0xsxVm1mxmr5rZTYfY9mwz6zKzq8KaR2Sw3J1l/UrmtCwk8RdmMUoXsMDd15pZDbDGzJa6e3PfjcxsKPAj4PEQZxEZtHwlc5+ZqZI5ib/Qjgjcfau7r83e3kPmoy0n5dn0RuBhYFue50TKRnBZ6A9PVMmcJENJFjfNrB44E1gdeHwS8CfAPxzm1883syYza2prawttTpFD6XfZqJaFJCFCDwIzG0nmJ/6b3X134OmfAt9y955+v7APd1/o7o3u3lhXp0NxKb33d+/nJZXMSUKFelxrZpVkQuB+d1+cZ5NG4IHsm3HGAZebWZe7PxLmXCJHKniS+IzJozlulErmJBlCCwLL/Ou+CGhx9zvybePu0/ps/wvgPxQCUo6WaVlIEizMI4LzgauB9Wa2LvvYbcBUAHe/J8SvLVI0eUvm9NnEkiChBYG7rwQKLmBx92vCmkVkMJ7emFsyN3XscGYeNzLCiUSKS2+JFDmMfB9JqZI5SRIFgcghdHX3sPy13Le46GohSRoFgcghvLA5t2RuzHCVzEnyKAhEDiF42eiFJ6tkTpJH39EiA3D3vOcHRJJGQSAygI3vd/DOjj4lcxUqmZNkUhCIDGBp83s598+fUcsIlcxJAikIRAYQXBa6WMtCklAKApE83t+9n5dad+U8pstGJakUBCJ59CuZmzJGJXOSWAoCkTyCy0KXallIEkxBIBLQcaCLZzcFS+YUBJJcCgKRgKc3tnGw++OSuRNqh3PSeJXMSXIpCEQC+r2JrEElc5JsCgKRPjq7e1i+IVAyp2UhSTgFgUgfTZt3smtfbslc4wkqmZNkUxCI9BFcFrrwFJXMSfLpO1wky91Z2pJbK6HLRiUNFAQiWa+9v4ctO/b13h9WMYRPn6SSOUk+BYFI1tJXc5eFPnXiOJXMSSooCESylgZqJdQtJGmhIBAB3tu1n5f7lcyNj2gakdJSEIjQv2Ru9pQxjFfJnKSEgkCEPO8m1tVCkiIKAkm9jgNdPPdGbsmcLhuVNFEQSOo99VpuyVx97XBOVMmcpIiCQFIveH7gYpXMScooCCTV8pXM6fyApI2CQFLthc07ckrmjh1eyR+oZE5SRkEgqda/ZO44lcxJ6oT2HW9mU8xshZk1m9mrZnZTnm2+amYvm9l6M3vWzM4Iax6RIHfXZaMiQJhFKl3AAndfa2Y1wBozW+ruzX22eQu4wN13mtllwELg3BBnEum14b09tO7MLZn7zMxxEU4kEo3QgsDdtwJbs7f3mFkLMAlo7rPNs31+ySpgcljziAQta+5fMjd8mErmJH1KshhqZvXAmcDqQ2z2deCxAX79fDNrMrOmtra2ECaUNAqWzGlZSNIq9CAws5HAw8DN7r57gG0+SyYIvpXveXdf6O6N7t5YV6d+eBm8YMmcGVykkjlJqVCPg82skkwI3O/uiwfY5nTgn4HL3H17vm1Eii14NDB7yhjG16hkTtIpzKuGDFgEtLj7HQNsMxVYDFzt7hvDmkUkSFcLiXwszCOC84GrgfVmti772G3AVAB3vwf4HlAL/H32Lf1d7t4Y4kwi7NnfyXNvtOc8ppI5SbMwrxpaCRyysMXdrweuD2sGkXye3thOZ7f33q+vHc6MOpXMSXrpLZSSOkub38u5f8mpKpmTdFMQSKrkL5mbENE0IuVBQSCp8sJbO9i9v6v3/tgRw1QyJ6mnIJBUebxfydx4hg7RspCkm4JAUkMlcyL5KQgkNTa8t4d3P/i4ZK6qYgifPkklcyIKAkmN4NGASuZEMhQEkhpaFhLJT0EgqbB11z7WvxssmVMQiICCQFIi+NkDZ04ZQ11NVUTTiJQXBYGkQvCyUb2JTORjCgJJvD37O1n1Zm7Duc4PiHxMQSCJ99TGtpySuWnjRjCjbkSEE4mUFwWBJF6+q4VUMifyMQWBJFpndw8r+pXMaVlIpC8FgSTa84GSudoRwzhrqkrmRPpSEEiiBZeFVDIn0p+CQBJLJXMihVEQSGK1bO1fMvcplcyJ9KMgkMQKHg18+iSVzInkoyCQxFra0v+ziUWkPwWBJNLvP9jHK+/u7r1vBheeoiAQyUdBIIm0rCV3WeisqceqZE5kAAoCSSRdLSRSOAWBJM7uPCVzF+uzB0QGpCCQxHnqtdySuenjRnDi+JERTiRS3hQEkjhaFhI5MgoCSZTO7h5WvKaSOZEjoSCQRHn+rR3sCZTMnamSOZFDUhBIogSXhS5qUMmcyOEoCCQx8pXM6WohkcMLrXjFzKYAvwSOAxxY6O53BbYx4C7gcuBD4Bp3X1vsWToOdHGgs7vYLytlZtO2jpySuerKIXz6pLoIJxKJhzAbuLqABe6+1sxqgDVmttTdm/tscxlwUva/c4F/yP6/qH702AbuW/V2sV9WytynTqzjmGFDox5DpOyFtjTk7ls/+une3fcALcCkwGZXAr/0jFXAGDObGNZMki6X6mohkYKU5ByBmdUDZwKrA09NArb0ud9K/7DAzOabWZOZNbW1tYU2pyTHyKoKLmoYH/UYIrEQejm7mY0EHgZudvfdh9s+H3dfCCwEaGxs9MNs3s/wqqGMHTHsaL60xNDxY6q5+aKZ1I5UyZxIIUINAjOrJBMC97v74jybvAtM6XN/cvaxovrOZQ1857KGYr+siEgihLY0lL0iaBHQ4u53DLDZEmCeZZwH7HL3rWHNJCIi/YV5RHA+cDWw3szWZR+7DZgK4O73AL8hc+noJjKXj14b4jwiIpJHaEHg7iuBQ76l090duCGsGURE5PD0zmIRkZRTEIiIpJyCQEQk5RQEIiIpZ5nztfFhZm3A0RYHjQPaizhOHGif00H7nA6D2ecT3D1vC2PsgmAwzKzJ3RujnqOUtM/poH1Oh7D2WUtDIiIppyAQEUm5tAXBwqgHiID2OR20z+kQyj6n6hyBiIj0l7YjAhERCVAQiIikXCKDwMw+Z2avmdkmM/t2nuerzOzB7POrs5+gFmsF7PNfmlmzmb1sZk+Y2QlRzFlMh9vnPtv9VzNzM4v9pYaF7LOZzc3+Wb9qZv9W6hmLrYDv7almtsLMXsx+f18exZzFYmb/YmbbzOyVAZ43M/tZ9vfjZTM7a9Bf1N0T9R8wFHgDmA4MA14CTg1s8+fAPdnbXwIejHruEuzzZ4Hh2dt/loZ9zm5XAzwNrAIao567BH/OJwEvAsdm74+Peu4S7PNC4M+yt08FNkc99yD3+TPAWcArAzx/OfAYmXbn84DVg/2aSTwiOAfY5O5vuvtB4AHgysA2VwL/mr39EHBR9oN04uqw++zuK9z9w+zdVWQ+DS7OCvlzBvgr4EfA/lIOF5JC9vkbwN+5+04Ad99W4hmLrZB9dmBU9vZo4PclnK/o3P1pYMchNrkS+KVnrALGmNnEwXzNJAbBJGBLn/ut2cfybuPuXcAuoLYk04WjkH3u6+tkfqKIs8Puc/aQeYq7P1rKwUJUyJ/zTGCmmf3OzFaZ2edKNl04Ctnn24GvmVkrmQ+7urE0o0XmSP++H1boH14v5cXMvgY0AhdEPUuYzGwIcAdwTcSjlFoFmeWhOWSO+p42s9Pc/YMohwrZl4FfuPvfmtkngfvMbJa790Q9WFwk8YjgXWBKn/uTs4/l3cbMKsgcTm4vyXThKGSfMbOLgf8JfN7dD5RotrAcbp9rgFnAk2a2mcxa6pKYnzAu5M+5FVji7p3u/hawkUwwxFUh+/x14P8CuPtzQDWZcrakKujv+5FIYhC8AJxkZtPMbBiZk8FLAtssAf40e/sqYLlnz8LE1GH32czOBP6RTAjEfd0YDrPP7r7L3ce5e72715M5L/J5d2+KZtyiKOR7+xEyRwOY2TgyS0VvlnDGYitkn98BLgIwswYyQdBW0ilLawkwL3v10HnALnffOpgXTNzSkLt3mdlfAL8lc8XBv7j7q2b2A6DJ3ZcAi8gcPm4ic1LmS9FNPHgF7vOPgZHAr7Lnxd9x989HNvQgFbjPiVLgPv8WuNTMmoFu4H+4e2yPdgvc5wXAP5nZLWROHF8T5x/szOzfyYT5uOx5j+8DlQDufg+Z8yCXA5uAD4FrB/01Y/z7JSIiRZDEpSERETkCCgIRkZRTEIiIpJyCQEQk5RQEIiIppyAQEUk5BYGISMopCEQGyczOzvbCV5vZiOznAMyKei6RQukNZSJFYGZ/Taba4Big1d3/T8QjiRRMQSBSBNkenBfIfO7BH7p7d8QjiRRMS0MixVFLpsuphsyRgUhs6IhApAjMbAmZT8+aBkx097+IeCSRgiWufVSk1MxsHtDp7v9mZkOBZ83sQndfHvVsIoXQEYGISMrpHIGISMopCEREUk5BICKScgoCEZGUUxCIiKScgkBEJOUUBCIiKff/AQmb+adZlzUhAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "L = 1.0\n", + "c = 0.5\n", + "dt = (L/10)/c # Nx=10\n", + "# I = lambda x: 0 if abs(x-L/2.0) > 0.1 else 1\n", + "I = lambda x: 5*x\n", + "U_0 = lambda t: t\n", + "T = 4\n", + "\n", + "\n", + "u, x, t, cpu = devito_solver(\n", + " I=I,\n", + " V=None, f=None, c=0.5, U_0=None, U_L=None, L=L,\n", + " dt=dt, C=1, T=dt, user_action=None)\n", "\n", - "[hpl 1: Should include some experiments here or make exercises. Qualitative\n", - "behavior of the wave equation can be exemplified.]\n", + "u_corr, x_corr, t_corr, cpu_corr = python_solver(\n", + " I=I,\n", + " V=None, f=None, c=0.5, U_0=None, U_L=None, L=L,\n", + " dt=dt, C=1, T=4*dt, user_action=None)\n", "\n", + "plt.plot(x_corr, u_corr, label='NumPy', linewidth=4)\n", + "# plt.plot(x, u, 'm:', label='Devito', linewidth=6)\n", + "plt.xlabel('x')\n", + "plt.ylabel('u')\n", + "plt.legend(loc='best')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ "## Verifying the implementation of Neumann conditions\n", "
\n", "\n", @@ -514,8 +775,8 @@ "How can we test that the Neumann conditions are correctly implemented?\n", "The `solver` function in the `wave1D_dn.py` program described in the\n", "box above accepts Dirichlet or Neumann conditions at $x=0$ and $x=L$.\n", - "mathcal{I}_t is tempting to apply a quadratic solution as described in\n", - "the sections [wave:pde2:fd](#wave:pde2:fd) and [wave:pde1:impl:verify:quadratic](#wave:pde1:impl:verify:quadratic),\n", + "It is tempting to apply a quadratic solution as described in\n", + "the sections [A slightly generalized model problem](wave1D_fd1.ipynb#wave:pde2:fd) and [Verification: exact quadratic solution](wave1D_prog.ipynb#wave:pde1:impl:verify:quadratic),\n", "but it turns out that this solution is no longer an exact solution\n", "of the discrete equations if a Neumann condition is implemented on\n", "the boundary. A linear solution does not help since we only have\n", @@ -525,10 +786,11 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 11, "metadata": {}, "outputs": [], "source": [ + "# %load -s test_constant, src-wave/wave1D/wave1D_dn.py\n", "def test_constant():\n", " \"\"\"\n", " Check the scalar and vectorized versions for\n", @@ -563,7 +825,7 @@ " solver(I, V, f, c, U_0, U_L, L, dt, C, T,\n", " user_action=assert_no_error,\n", " version='vectorized')\n", - " print U_0, U_L" + " print(U_0, U_L)\n" ] }, { @@ -584,10 +846,11 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 12, "metadata": {}, "outputs": [], "source": [ + "# %load -s test_plug, src-wave/wave1D/wave1D_dn.py\n", "def test_plug():\n", " \"\"\"Check that an initial plug is correct back after one period.\"\"\"\n", " L = 1.0\n", @@ -608,299 +871,26 @@ " assert diff < tol\n", " u_0 = np.array([I(x_) for x_ in x])\n", " diff = np.abs(u_s - u_0).max()\n", - " assert diff < tol" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Other tests must rely on an unknown approximation error, so effectively\n", - "we are left with tests on the convergence rate.\n", - "\n", - "## Alternative implementation via ghost cells\n", - "
\n", - "\n", - "### Idea\n", - "\n", - "Instead of modifying the scheme at the boundary, we can introduce\n", - "extra points outside the domain such that the fictitious values\n", - "$u_{-1}^n$ and $u_{N_x+1}^n$ are defined in the mesh. Adding the\n", - "intervals $[-\\Delta x,0]$ and $[L, L+\\Delta x]$, known as *ghost\n", - "cells*, to the mesh gives us all the needed mesh points, corresponding\n", - "to $i=-1,0,\\ldots,N_x,N_x+1$. The extra points with $i=-1$ and\n", - "$i=N_x+1$ are known as *ghost points*, and values at these points,\n", - "$u_{-1}^n$ and $u_{N_x+1}^n$, are called *ghost values*.\n", - "\n", - "The important idea is\n", - "to ensure that we always have" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "$$\n", - "u_{-1}^n = u_{1}^n\\hbox{ and } u_{N_x+1}^n = u_{N_x-1}^n,\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "because then\n", - "the application of the standard scheme at a boundary point $i=0$ or $i=N_x$\n", - "will be correct and guarantee that the solution is compatible with the\n", - "boundary condition $u_x=0$.\n", - "\n", - "Some readers may find it strange to just extend the domain with ghost\n", - "cells as a general technique, because in some problems there is a\n", - "completely different medium with different physics and equations right\n", - "outside of a boundary. Nevertheless, one should view the ghost cell\n", - "technique as a purely mathematical technique, which is valid in the\n", - "limit $\\Delta x \\rightarrow 0$ and helps us to implement derivatives.\n", - "\n", - "\n", - "### Implementation\n", - "\n", - "The `u` array now needs extra elements corresponding to the ghost\n", - "points. Two new point values are needed:" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "u = zeros(Nx+3)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The arrays `u_n` and `u_nm1` must be defined accordingly.\n", - "\n", - "Unfortunately, a major indexing problem arises with ghost cells.\n", - "The reason is that Python indices *must* start\n", - "at 0 and `u[-1]` will always mean the last element in `u`.\n", - "This fact gives, apparently, a mismatch between the mathematical\n", - "indices $i=-1,0,\\ldots,N_x+1$ and the Python indices running over\n", - "`u`: `0,..,Nx+2`. One remedy is to change the mathematical indexing\n", - "of $i$ in the scheme and write" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "$$\n", - "u^{n+1}_i = \\cdots,\\quad i=1,\\ldots,N_x+1,\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "instead of $i=0,\\ldots,N_x$ as we have previously used. The ghost\n", - "points now correspond to $i=0$ and $i=N_x+1$.\n", - "A better solution is to use the ideas of the section [Index set notation](#wave:indexset):\n", - "we hide the specific index value in an index set and operate with\n", - "inner and boundary points using the index set notation.\n", - "\n", - "To this end, we define `u` with proper length and `mathcal{I}_x` to be the corresponding\n", - "indices for the real physical mesh points ($1,2,\\ldots,N_x+1$):" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - " u = zeros(Nx+3)\n", - " mathcal{I}_x = range(1, u.shape[0]-1)\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "That is, the boundary points have indices `mathcal{I}_x[0]` and `mathcal{I}_x[-1]` (as before).\n", - "We first update the solution at all physical mesh points (i.e., interior\n", - "points in the mesh):" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [], - "source": [ - "for i in mathcal{I}_x:\n", - " u[i] = - u_nm1[i] + 2*u_n[i] + \\\n", - " C2*(u_n[i-1] - 2*u_n[i] + u_n[i+1])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The indexing becomes a bit more complicated when we call functions like\n", - "`V(x)` and `f(x, t)`, as we must remember that the appropriate\n", - "$x$ coordinate is given as `x[i-mathcal{I}_x[0]]`:" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [], - "source": [ - "for i in mathcal{I}_x:\n", - " u[i] = u_n[i] + dt*V(x[i-mathcal{I}_x[0]]) + \\\n", - " 0.5*C2*(u_n[i-1] - 2*u_n[i] + u_n[i+1]) + \\\n", - " 0.5*dt2*f(x[i-mathcal{I}_x[0]], t[0])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "mathcal{I}_t remains to update the solution at ghost points, i.e., `u[0]`\n", - "and `u[-1]` (or `u[Nx+2]`). For a boundary condition $u_x=0$,\n", - "the ghost value must equal the value at the associated inner mesh\n", - "point. Computer code makes this statement precise:" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [], - "source": [ - "i = mathcal{I}_x[0] # x=0 boundary\n", - "u[i-1] = u[i+1]\n", - "i = mathcal{I}_x[-1] # x=L boundary\n", - "u[i+1] = u[i-1]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The physical solution to be plotted is now in `u[1:-1]`, or\n", - "equivalently `u[mathcal{I}_x[0]:mathcal{I}_x[-1]+1]`, so this slice is\n", - "the quantity to be returned from a solver function.\n", - "A complete implementation appears in the program\n", - "[`wave1D_n0_ghost.py`](${src_wave}/wave1D/wave1D_n0_ghost.py).\n", - "\n", - "**Warning.**\n", - "\n", - "We have to be careful with how the spatial and temporal mesh\n", - "points are stored. Say we let `x` be the physical mesh points," - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [], - "source": [ - "x = linspace(0, L, Nx+1)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\"Standard coding\" of the initial condition," - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [], - "source": [ - "for i in mathcal{I}_x:\n", - " u_n[i] = I(x[i])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "becomes wrong, since `u_n` and `x` have different lengths and the index `i`\n", - "corresponds to two different mesh points. In fact, `x[i]` corresponds\n", - "to `u[1+i]`. A correct implementation is" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": {}, - "outputs": [], - "source": [ - "for i in mathcal{I}_x:\n", - " u_n[i] = I(x[i-mathcal{I}_x[0]])" + " assert diff < tol\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Similarly, a source term usually coded as `f(x[i], t[n])` is incorrect\n", - "if `x` is defined to be the physical points, so `x[i]` must be\n", - "replaced by `x[i-mathcal{I}_x[0]]`.\n", - "\n", - "An alternative remedy is to let `x` also cover the ghost points such that\n", - "`u[i]` is the value at `x[i]`.\n", - "\n", - "\n", - "\n", - "The ghost cell is only added to the boundary where we have a Neumann\n", - "condition. Suppose we have a Dirichlet condition at $x=L$ and\n", - "a homogeneous Neumann condition at $x=0$. One ghost cell $[-\\Delta x,0]$\n", - "is added to the mesh, so the index set for the physical points\n", - "becomes $\\{1,\\ldots,N_x+1\\}$. A relevant implementation\n", - "is" - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "metadata": {}, - "outputs": [], - "source": [ - "u = zeros(Nx+2)\n", - "mathcal{I}_x = range(1, u.shape[0])\n", - "...\n", - "for i in mathcal{I}_x[:-1]:\n", - " u[i] = - u_nm1[i] + 2*u_n[i] + \\\n", - " C2*(u_n[i-1] - 2*u_n[i] + u_n[i+1]) + \\\n", - " dt2*f(x[i-mathcal{I}_x[0]], t[n])\n", - "i = mathcal{I}_x[-1]\n", - "u[i] = U_0 # set Dirichlet value\n", - "i = mathcal{I}_x[0]\n", - "u[i-1] = u[i+1] # update ghost value" + "Other tests must rely on an unknown approximation error, so effectively we are left with tests on the convergence rate.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The physical solution to be plotted is now in `u[1:]`\n", - "or (as always) `u[mathcal{I}_x[0]:mathcal{I}_x[-1]+1]`.\n", - "\n", "\n", "# Generalization: variable wave velocity\n", "
\n", "\n", "\n", - "Our next generalization of the 1D wave equation ([wave:pde1](#wave:pde1)) or\n", - "([wave:pde2](#wave:pde2)) is to allow for a variable wave velocity $c$:\n", + "Our next generalization of the 1D wave equation is to allow for a variable wave velocity $c$:\n", "$c=c(x)$, usually motivated by wave motion in a domain composed of\n", "different physical media. When the media differ in physical properties\n", "like density or porosity, the wave velocity $c$ is affected and\n", @@ -956,8 +946,8 @@ "source": [ "This is the most frequent form of a wave\n", "equation with variable wave velocity,\n", - "but other forms also appear, see the section [wave:app:string](#wave:app:string)\n", - "and equation ([wave:app:string:model2](#wave:app:string:model2)).\n", + "but other forms also appear, see the section [Waves on a string](wave_app.ipynb#wave:app:string)\n", + "and [this equation](wave_app.ipynb#wave:app:string:model2).\n", "\n", "As usual, we sample ([8](#wave:pde2:var:c:pde)) at a mesh point," ] @@ -1281,11 +1271,7 @@ "jumps (which is typical for geological media).\n", "The geometric mean is less used, but popular in\n", "discretizations to linearize quadratic\n", - "% if BOOK == \"book\":\n", - "nonlinearities (see the section [vib:ode2:fdm:fquad](#vib:ode2:fdm:fquad) for an example).\n", - "% else:\n", - "nonlinearities.\n", - "% endif\n", + "nonlinearities (see the section [A centred scheme for quadratic damping](../01_vib/vib_gen.ipynb#vib:ode2:fdm:fquad) for an example).\n", "\n", "With the operator notation from ([12](#wave:pde2:var:c:mean:arithmetic))\n", "we can specify the discretization of the complete variable-coefficient\n", @@ -1370,7 +1356,7 @@ "
\n", "\n", "\n", - "The stability criterion derived later (the section [wave:pde1:stability](#wave:pde1:stability))\n", + "The stability criterion derived later (the section [Stability](wave_analysis.ipynb#wave:pde1:stability))\n", "reads $\\Delta t\\leq \\Delta x/c$. If $c=c(x)$, the criterion will depend\n", "on the spatial location. We must therefore choose a $\\Delta t$ that\n", "is small enough such that no mesh cell has $\\Delta t > \\Delta x/c(x)$.\n", @@ -1642,10 +1628,44 @@ "## Implementation of variable coefficients\n", "
\n", "\n", - "The implementation of the scheme with a variable wave velocity $q(x)=c^2(x)$\n", - "may assume that $q$ is available as an array `q[i]` at\n", - "the spatial mesh points. The following loop is a straightforward\n", - "implementation of the scheme ([16](#wave:pde2:var:c:scheme:impl)):" + "The implementation of the scheme with a variable wave velocity $q(x)=c^2(x)$ is very simple using Devito. Previously, we used this simpler form of the wave equation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\begin{equation}\n", + "\\frac{\\partial^2 u}{\\partial t^2} =\n", + "c^2\\frac{\\partial^2 u}{\\partial x^2}\n", + " + f(x,t) \\thinspace .\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We represented this equation in our code as" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "grid = Grid(shape=(Nx+1), extent=(L))\n", + "u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "with our `pde` variable simply being" ] }, { @@ -1654,29 +1674,14 @@ "metadata": {}, "outputs": [], "source": [ - "for i in range(1, Nx):\n", - " u[i] = - u_nm1[i] + 2*u_n[i] + \\\n", - " C2*(0.5*(q[i] + q[i+1])*(u_n[i+1] - u_n[i]) - \\\n", - " 0.5*(q[i] + q[i-1])*(u_n[i] - u_n[i-1])) + \\\n", - " dt2*f(x[i], t[n])" + "pde = (1/c**2)*u.dt2-u.dx2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The coefficient `C2` is now defined as `(dt/dx)**2`, i.e., *not* as the\n", - "squared Courant number, since the wave velocity is variable and appears\n", - "inside the parenthesis.\n", - "\n", - "With Neumann conditions $u_x=0$ at the\n", - "boundary, we need to combine this scheme with the discrete\n", - "version of the boundary condition, as shown in the section [Neumann condition and a variable coefficient](#wave:pde2:var:c:Neumann).\n", - "Nevertheless, it would be convenient to reuse the formula for the\n", - "interior points and just modify the indices `ip1=i+1` and `im1=i-1`\n", - "as we did in the section [Implementation of Neumann conditions](#wave:pde2:Neumann:impl). Assuming\n", - "$dq/dx=0$ at the boundaries, we can implement the scheme at\n", - "the boundary with the following code." + "For incorporating the variable wave velocity term, we first need to create a Devito `Function` to represent our given function `q` and initialise its values." ] }, { @@ -1685,25 +1690,19 @@ "metadata": {}, "outputs": [], "source": [ - "i = 0\n", - "ip1 = i+1\n", - "im1 = ip1\n", - "u[i] = - u_nm1[i] + 2*u_n[i] + \\\n", - " C2*(0.5*(q[i] + q[ip1])*(u_n[ip1] - u_n[i]) - \\\n", - " 0.5*(q[i] + q[im1])*(u_n[i] - u_n[im1])) + \\\n", - " dt2*f(x[i], t[n])" + "# Setting q as identity function for demonstrative purposes\n", + "# In reality, this function q would be passed as an argument to the solver function\n", + "q = lambda x: x\n", + "\n", + "q_f = Function(name='q', grid=grid, npoint=Nx+1, nt=Nt+1)\n", + "q_f.data[:] = q(x[:])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "With ghost cells we can just reuse the formula for the interior\n", - "points also at the boundary, provided that the ghost values of both\n", - "$u$ and $q$ are correctly updated to ensure $u_x=0$ and $q_x=0$.\n", - "\n", - "A vectorized version of the scheme with a variable coefficient\n", - "at internal mesh points becomes" + "We can now set our `pde` variable as" ] }, { @@ -1712,10 +1711,20 @@ "metadata": {}, "outputs": [], "source": [ - "u[1:-1] = - u_nm1[1:-1] + 2*u_n[1:-1] + \\\n", - " C2*(0.5*(q[1:-1] + q[2:])*(u_n[2:] - u_n[1:-1]) -\n", - " 0.5*(q[1:-1] + q[:-2])*(u_n[1:-1] - u_n[:-2])) + \\\n", - " dt2*f(x[1:-1], t[n])" + "pde = u.dt2-(q_f*u.dx).dx" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "and the rest of the steps are now identical to what we did in our `solver` functions before.\n", + "\n", + "**Note**\n", + "\n", + "The coefficient `C2` is now defined as `(dt/dx)**2`, i.e., *not* as the\n", + "squared Courant number, since the wave velocity is variable and appears\n", + "inside the parenthesis." ] }, { @@ -1752,7 +1761,7 @@ "metadata": {}, "source": [ "One example appears when modeling elastic waves in a rod\n", - "with varying density, cf. ([wave:app:string](#wave:app:string)) with $\\varrho (x)$.\n", + "with varying density, cf. ([Waves on a string](wave_app.ipynb#wave:app:string)) with $\\varrho (x)$.\n", "\n", "A natural scheme for ([24](#wave:pde2:var:c:pde2)) is" ] @@ -1891,7 +1900,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "for $i\\in\\mathcal{I}_x^i$ and $n\\geq 1$.\n", + "for $i\\in \\mathcal{I}_x^i$ and $n\\geq 1$.\n", "New equations must be derived for $u^1_i$, and for boundary points in case\n", "of Neumann conditions.\n", "\n", @@ -1904,7 +1913,7 @@ "
\n", "\n", "\n", - "The program [`wave1D_dn_vc.py`](${src_wave}/wave1D/wave1D_dn_vc.py)\n", + "The program [`wave1D_dn_vc.py`](src-wave/wave1D/wave1D_dn_vc.py)\n", "is a fairly general code for 1D wave propagation problems that\n", "targets the following initial-boundary value problem" ] @@ -1984,27 +1993,6 @@ "$$" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The only new feature here is the time-dependent Dirichlet conditions, but\n", - "they are trivial to implement:" - ] - }, - { - "cell_type": "code", - "execution_count": 16, - "metadata": {}, - "outputs": [], - "source": [ - "i = mathcal{I}_x[0] # x=0\n", - "u[i] = U_0(t[n+1])\n", - "\n", - "i = mathcal{I}_x[-1] # x=L\n", - "u[i] = U_L(t[n+1])" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -2016,16 +2004,7 @@ "a variable wave velocity. The different code segments needed\n", "to make these extensions have been shown and commented upon in the\n", "preceding text. We refer to the `solver` function in the\n", - "`wave1D_dn_vc.py` file for all the details. Note in that\n", - " `solver` function, however, that the technique of \"hashing\" is\n", - "used to check whether a certain simulation has been run before, or not.\n", - "% if BOOK == 'book':\n", - "This technique is further explained in the section [softeng2:wave1D:filestorage:hash](#softeng2:wave1D:filestorage:hash).\n", - "% endif\n", - "\n", - "The vectorization is only applied inside the time loop, not for the\n", - "initial condition or the first time steps, since this initial work\n", - "is negligible for long time simulations in 1D problems.\n", + "`wave1D_dn_vc.py` file for all the details.\n", "\n", "The following sections explain various more advanced programming\n", "techniques applied in the general 1D wave equation solver.\n", @@ -2048,7 +2027,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 24, "metadata": {}, "outputs": [], "source": [ @@ -2187,11 +2166,6 @@ "saves the plot to file, and stores the solution in a file for\n", "later retrieval.\n", "\n", - "More details on storing the solution in files appear in\n", - "in\n", - "the document\n", - "[Scientific software engineering; wave equation case](http://tinyurl.com/k3sdbuv/pub/softeng2)\n", - "[[Langtangen_deqbook_softeng2]](#Langtangen_deqbook_softeng2).\n", "\n", "## Pulse propagation in two media\n", "\n", @@ -2212,11 +2186,11 @@ "3. a \"cosine hat\" consisting of one period of the cosine function\n", " (`cosinehat`),\n", "\n", - "4. frac{1}{2} a period of a \"cosine hat\" (`frac{1}{2}-cosinehat`)\n", + "4. half a period of a \"cosine hat\" (`half-cosinehat`)\n", "\n", "These peak-shaped initial conditions can be placed in the middle\n", "(`loc='center'`) or at the left end (`loc='left'`) of the domain.\n", - "With the pulse in the middle, it splits in two parts, each with frac{1}{2}\n", + "With the pulse in the middle, it splits in two parts, each with half\n", "the initial amplitude, traveling in opposite directions. With the\n", "pulse at the left end, centered at $x=0$, and using the symmetry\n", "condition $\\partial u/\\partial x=0$, only a right-going pulse is\n", @@ -2232,10 +2206,11 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 121, "metadata": {}, "outputs": [], "source": [ + "# %load -s pulse, src-wave/wave1D/wave1D_dn_vc.py\n", "def pulse(\n", " C=1, # Maximum Courant number\n", " Nx=200, # spatial resolution\n", @@ -2243,7 +2218,7 @@ " version='vectorized',\n", " T=2, # end time\n", " loc='left', # location of initial condition\n", - " pulse_thinspace .='gaussian', # pulse/init.cond. type\n", + " pulse_tp='gaussian', # pulse/init.cond. type\n", " slowness_factor=2, # inverse of wave vel. in right medium\n", " medium=[0.7, 0.9], # interval for right medium\n", " skip_frame=1, # skip frames in animations\n", @@ -2265,13 +2240,13 @@ " elif loc == 'left':\n", " xc = 0\n", "\n", - " if pulse_thinspace . in ('gaussian','Gaussian'):\n", + " if pulse_tp in ('gaussian','Gaussian'):\n", " def I(x):\n", " return np.exp(-0.5*((x-xc)/sigma)**2)\n", - " elif pulse_thinspace . == 'plug':\n", + " elif pulse_tp == 'plug':\n", " def I(x):\n", " return 0 if abs(x-xc) > sigma else 1\n", - " elif pulse_thinspace . == 'cosinehat':\n", + " elif pulse_tp == 'cosinehat':\n", " def I(x):\n", " # One period of a cosine\n", " w = 2\n", @@ -2279,7 +2254,7 @@ " return 0.5*(1 + np.cos(np.pi*(x-xc)/a)) \\\n", " if xc - a <= x <= xc + a else 0\n", "\n", - " elif pulse_thinspace . == 'frac{1}{2}-cosinehat':\n", + " elif pulse_tp == 'half-cosinehat':\n", " def I(x):\n", " # Half a period of a cosine\n", " w = 4\n", @@ -2287,7 +2262,7 @@ " return np.cos(np.pi*(x-xc)/a) \\\n", " if xc - 0.5*a <= x <= xc + 0.5*a else 0\n", " else:\n", - " raise ValueError('Wrong pulse_thinspace .=\"%s\"' % pulse_thinspace .)\n", + " raise ValueError('Wrong pulse_tp=\"%s\"' % pulse_tp)\n", "\n", " def c(x):\n", " return c_0/slowness_factor \\\n", @@ -2295,7 +2270,7 @@ "\n", " umin=-0.5; umax=1.5*I(xc)\n", " casename = '%s_Nx%s_sf%s' % \\\n", - " (pulse_thinspace ., Nx, slowness_factor)\n", + " (pulse_tp, Nx, slowness_factor)\n", " action = PlotMediumAndSolution(\n", " medium, casename=casename, umin=umin, umax=umax,\n", " skip_frame=skip_frame, screen_movie=animate,\n", @@ -2304,7 +2279,7 @@ " # Choose the stability limit with given Nx, worst case c\n", " # (lower C will then use this dt, but smaller Nx)\n", " dt = (L/Nx)/c_0\n", - " cpu, hashed_input = solver(\n", + " u, x, t, cpu = devito_solver(\n", " I=I, V=None, f=None, c=c,\n", " U_0=None, U_L=None,\n", " L=L, dt=dt, C=C, T=T,\n", @@ -2315,68 +2290,7 @@ " if cpu > 0: # did we generate new data?\n", " action.close_file(hashed_input)\n", " action.make_movie_file()\n", - " print 'cpu (-1 means no new data generated):', cpu\n", - "\n", - "def convergence_rates(\n", - " u_exact,\n", - " I, V, f, c, U_0, U_L, L,\n", - " dt0, num_meshes,\n", - " C, T, version='scalar',\n", - " stability_safety_factor=1.0):\n", - " \"\"\"\n", - " Half the time step and estimate convergence rates for\n", - " for num_meshes simulations.\n", - " \"\"\"\n", - " class ComputeError:\n", - " def __init__(self, norm_type):\n", - " self.error = 0\n", - "\n", - " def __call__(self, u, x, t, n):\n", - " \"\"\"Store norm of the error in self.E.\"\"\"\n", - " error = np.abs(u - u_exact(x, t[n])).max()\n", - " self.error = max(self.error, error)\n", - "\n", - " E = []\n", - " h = [] # dt, solver adjusts dx such that C=dt*c/dx\n", - " dt = dt0\n", - " for i in range(num_meshes):\n", - " error_calculator = ComputeError('Linf')\n", - " solver(I, V, f, c, U_0, U_L, L, dt, C, T,\n", - " user_action=error_calculator,\n", - " version='scalar',\n", - " stability_safety_factor=1.0)\n", - " E.append(error_calculator.error)\n", - " h.append(dt)\n", - " dt /= 2 # halve the time step for next simulation\n", - " print 'E:', E\n", - " print 'h:', h\n", - " r = [np.log(E[i]/E[i-1])/np.log(h[i]/h[i-1])\n", - " for i in range(1,num_meshes)]\n", - " return r\n", - "\n", - "def test_convrate_sincos():\n", - " n = m = 2\n", - " L = 1.0\n", - " u_exact = lambda x, t: np.cos(m*np.pi/L*t)*np.sin(m*np.pi/L*x)\n", - "\n", - " r = convergence_rates(\n", - " u_exact=u_exact,\n", - " I=lambda x: u_exact(x, 0),\n", - " V=lambda x: 0,\n", - " f=0,\n", - " c=1,\n", - " U_0=0,\n", - " U_L=0,\n", - " L=L,\n", - " dt0=0.1,\n", - " num_meshes=6,\n", - " C=0.9,\n", - " T=1,\n", - " version='scalar',\n", - " stability_safety_factor=1.0)\n", - " print 'rates sin(x)*cos(t) solution:', \\\n", - " [round(r_,2) for r_ in r]\n", - " assert abs(r[-1] - 2) < 0.002" + " print('cpu (-1 means no new data generated):', cpu)\n" ] }, { @@ -2440,7 +2354,7 @@ "Consider the wave equation with damping ([27](#wave:pde3)).\n", "The goal is to find an exact solution to a wave problem with damping and zero source term.\n", "A starting point is the standing wave solution from\n", - "[wave:exer:standingwave](#wave:exer:standingwave). mathcal{I}_t becomes necessary to\n", + "[Simulate a standing wave exercise](wave1D_prog.ipynb#wave:exer:standingwave). It becomes necessary to\n", "include a damping term $e^{-\\beta t}$ and also have both a sine and cosine\n", "component in time:" ] @@ -2450,7 +2364,7 @@ "metadata": {}, "source": [ "$$\n", - "\\uex(x,t) = e^{-\\beta t}\n", + "u(x,t) = e^{-\\beta t}\n", "\\sin kx \\left( A\\cos\\omega t\n", "+ B\\sin\\omega t\\right)\n", "\\thinspace .\n", @@ -2715,7 +2629,7 @@ "metadata": {}, "source": [ "$$\n", - "\\uex(x,t) = e^{-\\beta t}\n", + "u(x,t) = e^{-\\beta t}\n", "\\sin kx \\left( A\\cos\\omega t\n", "+ B\\sin\\omega t\\right)\n", "\\thinspace .\n", @@ -2766,7 +2680,7 @@ "$u_t(x,0)=0$ and there is no source term $f$.\n", "The boundary conditions can be set to $u=0$.\n", "The solution to this problem is symmetric around $x=0$.\n", - "This means that we can simulate the wave process in only frac{1}{2}\n", + "This means that we can simulate the wave process in only half\n", "of the domain $[0,L]$.\n", "\n", "\n", @@ -2814,7 +2728,7 @@ "Perform simulations of the complete wave problem\n", "on $[-L,L]$. Thereafter, utilize the\n", "symmetry of the solution and run a simulation\n", - "in frac{1}{2} of the domain $[0,L]$, using a boundary condition\n", + "in half of the domain $[0,L]$, using a boundary condition\n", "at $x=0$. Compare plots from the two solutions and\n", "confirm that they are the same.\n", "\n", @@ -2829,11 +2743,11 @@ "$u(-L,t)=u(L,t)=0$ and $u_x(0,t)=u(L,t)=0$, respectively.\n", "\n", "The original `wave1D_dn.py` code makes a movie by playing all the\n", - "`.png` files in a browser. mathcal{I}_t can then be wise to let the `viz`\n", + "`.png` files in a browser. It can then be wise to let the `viz`\n", "function create a movie directory and place all the frames and HTML\n", "player file in that directory. Alternatively, one can just make\n", "some ordinary movie file (Ogg, WebM, MP4, Flash) with `ffmpeg` or\n", - "`ffmpeg` and give it a name. mathcal{I}_t is a point that the name is\n", + "`ffmpeg` and give it a name. It is a point that the name is\n", "transferred to `viz` so it is easy to call `viz` twice and get two\n", "separate movie files or movie directories.\n", "\n", @@ -3240,9 +3154,21 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 18, "metadata": {}, - "outputs": [], + "outputs": [ + { + "ename": "ModuleNotFoundError", + "evalue": "No module named 'scitools'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;31m#!/usr/bin/env python3\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mscitools\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstd\u001b[0m \u001b[0;32mimport\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 3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0;31m# Add an x0 coordinate for solving the wave equation on [x0, xL]\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[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'scitools'" + ] + } + ], "source": [ "#!/usr/bin/env python\n", "from scitools.std import *\n", @@ -3280,25 +3206,25 @@ " u_1 = zeros(Nx+1) # Solution at 1 time level back\n", " u_2 = zeros(Nx+1) # Solution at 2 time levels back\n", "\n", - " mathcal{I}_x = range(0, Nx+1)\n", - " mathcal{I}_t = range(0, Nt+1)\n", + " Ix = range(0, Nx+1)\n", + " It = range(0, Nt+1)\n", "\n", " import time; t0 = time.clock() # CPU time measurement\n", "\n", " # Load initial condition into u_1\n", - " for i in mathcal{I}_x:\n", + " for i in Ix:\n", " u_1[i] = I(x[i])\n", "\n", " if user_action is not None:\n", " user_action(u_1, x, t, 0)\n", "\n", " # Special formula for the first step\n", - " for i in mathcal{I}_x[1:-1]:\n", + " for i in Ix[1:-1]:\n", " u[i] = u_1[i] + dt*V(x[i]) + \\\n", " 0.5*C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1]) + \\\n", " 0.5*dt2*f(x[i], t[0])\n", "\n", - " i = mathcal{I}_x[0]\n", + " i = Ix[0]\n", " if U_0 is None:\n", " # Set boundary values du/dn = 0\n", " # x=0: i-1 -> i+1 since u[i-1]=u[i+1]\n", @@ -3311,7 +3237,7 @@ " else:\n", " u[0] = U_0(dt)\n", "\n", - " i = mathcal{I}_x[-1]\n", + " i = Ix[-1]\n", " if U_L is None:\n", " im1 = i-1\n", " ip1 = im1 # i+1 -> i-1\n", @@ -3327,10 +3253,10 @@ " # Update data structures for next step\n", " u_2[:], u_1[:] = u_1, u\n", "\n", - " for n in mathcal{I}_t[1:-1]:\n", + " for n in It[1:-1]:\n", " # Update all inner points\n", " if version == 'scalar':\n", - " for i in mathcal{I}_x[1:-1]:\n", + " for i in Ix[1:-1]:\n", " u[i] = - u_2[i] + 2*u_1[i] + \\\n", " C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1]) + \\\n", " dt2*f(x[i], t[n])\n", @@ -3343,7 +3269,7 @@ " raise ValueError('version=%s' % version)\n", "\n", " # Insert boundary conditions\n", - " i = mathcal{I}_x[0]\n", + " i = Ix[0]\n", " if U_0 is None:\n", " # Set boundary values\n", " # x=0: i-1 -> i+1 since u[i-1]=u[i+1] when du/dn=0\n", @@ -3356,7 +3282,7 @@ " else:\n", " u[0] = U_0(t[n+1])\n", "\n", - " i = mathcal{I}_x[-1]\n", + " i = Ix[-1]\n", " if U_L is None:\n", " im1 = i-1\n", " ip1 = im1\n", @@ -3461,14 +3387,14 @@ " return 1\n", "\n", " # Solution on [-L,L]\n", - " cpu = viz(I=I, V=0, f=0, c, U_0=0, U_L=0, -L, L, 2*Nx, C, T,\n", + " cpu = viz(I=I, V=0, f=0, c=c, U_0=0, U_L=0, x0=-L, xL=L, Nx=2*Nx, C=C, T=T,\n", " umin=-1.1, umax=1.1, version=version, animate=animate,\n", " movie_dir='full')\n", "\n", " # Solution on [0,L]\n", - " cpu = viz(I=I, V=0, f=0, c, U_0=None, U_L=0, 0, L, Nx, C, T,\n", + " cpu = viz(I=I, V=0, f=0, c=c, U_0=None, U_L=0, x0=0, xL=L, Nx=Nx, C=C, T=T,\n", " umin=-1.1, umax=1.1, version=version, animate=animate,\n", - " movie_dir='frac{1}{2}')\n", + " movie_dir='half')\n", "\n", "\n", "if __name__ == '__main__':\n", @@ -3499,7 +3425,7 @@ "media with different wave velocities. The (scaled) velocity in\n", "the left medium is 1 while it is $\\frac{1}{s_f}$ in the right medium.\n", "Report what happens with a Gaussian pulse, a \"cosine hat\" pulse,\n", - "frac{1}{2} a \"cosine hat\" pulse, and a plug pulse for resolutions\n", + "half a \"cosine hat\" pulse, and a plug pulse for resolutions\n", "$N_x=40,80,160$, and $s_f=2,4$. Simulate until $T=2$.\n", "\n", "\n", @@ -3519,7 +3445,7 @@ "import wave1D_dn_vc as wave\n", "import os, sys, shutil, glob\n", "\n", - "for pulse_thinspace . in 'gaussian', 'cosinehat', 'frac{1}{2}-cosinehat', 'plug':\n", + "for pulse_thinspace . in 'gaussian', 'cosinehat', 'half-cosinehat', 'plug':\n", " for Nx in 40, 80, 160:\n", " for sf in 2, 4:\n", " if sf == 1 and Nx > 40:\n", @@ -3553,9 +3479,9 @@ "\n", "The experiments performed in [Exercise 3: Send pulse waves through a layered medium](#wave:app:exer:pulse1D) shows\n", "considerable numerical noise in the form of non-physical waves,\n", - "especially for $s_f=4$ and the plug pulse or the frac{1}{2} a \"cosinehat\"\n", + "especially for $s_f=4$ and the plug pulse or the half a \"cosinehat\"\n", "pulse. The noise is much less visible for a Gaussian pulse. Run the\n", - "case with the plug and frac{1}{2} a \"cosinehat\" pulse for $s_f=1$, $C=0.9,\n", + "case with the plug and half a \"cosinehat\" pulse for $s_f=1$, $C=0.9,\n", "0.25$, and $N_x=40,80,160$. Use the numerical dispersion relation to\n", "explain the observations.\n", "Filename: `pulse1D_analysis`.\n", @@ -3743,7 +3669,7 @@ "\n", "**Hint.**\n", "Modify the solver function in\n", - "[`wave1D_dn.py`](${src_wave}/wave1D/wave1D_dn.py).\n", + "[`wave1D_dn.py`](src-wave/wave1D/wave1D_dn.py).\n", "\n", "\n", "\n", @@ -3837,7 +3763,7 @@ "
\n", "\n", "\n", - "mathcal{I}_t is frequently of interest to follow wave motion over large\n", + "It is frequently of interest to follow wave motion over large\n", "distances and long times. A straightforward approach is to\n", "work with a very large domain, but that might lead to a lot of\n", "computations in areas of the domain where the waves cannot\n", @@ -3946,7 +3872,7 @@ "we can just set $[qD_xu]_{i+\\frac{1}{2}}^n=0$.\n", "Derive the complete scheme\n", "using this technique. The implementation of the boundary condition at\n", - "$L-\\Delta x/2$ is $\\Oof{\\Delta x^2}$ accurate, but the interesting question\n", + "$L-\\Delta x/2$ is $\\mathcal{O}({\\Delta x^2})$ accurate, but the interesting question\n", "is what impact the movement of the boundary has on the convergence\n", "rate. Compute the errors as usual over the entire mesh and use $q$ from\n", "a) or b).\n", @@ -3978,12 +3904,12 @@ "this exact solution to machine precision. More precisely, we seek\n", "$u=X(x)T(t)$, with $T(t)$ as a linear function and $X(x)$ as a\n", "parabola that fulfills the boundary conditions. Inserting this $u$ in\n", - "the PDE determines $f$. mathcal{I}_t turns out that $u$ also fulfills the\n", + "the PDE determines $f$. It turns out that $u$ also fulfills the\n", "discrete equations, because the truncation error of the discretized\n", "PDE has derivatives in $x$ and $t$ of order four and higher. These\n", "derivatives all vanish for a quadratic $X(x)$ and linear $T(t)$.\n", "\n", - "mathcal{I}_t would be attractive to use a similar approach in the case of\n", + "It would be attractive to use a similar approach in the case of\n", "Neumann conditions. We set $u=X(x)T(t)$ and seek lower-order\n", "polynomials $X$ and $T$.\n", "To force $u_x$ to vanish at the boundary, we let $X_x$ be\n", @@ -4006,7 +3932,7 @@ "$$\n", "\\begin{equation}\n", "[D_{2x}u]^n_i = u_{x}(x_i,t_n) + \\frac{1}{6}u_{xxx}(x_i,t_n)\\Delta x^2\n", - "+ \\Oof{\\Delta x^4}\\thinspace .\n", + "+ \\mathcal{O}({\\Delta x^4})\\thinspace .\n", "\\label{wave:fd2:exer:verify:cubic:D2x} \\tag{60}\n", "\\end{equation}\n", "$$" @@ -4052,7 +3978,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "as it should. (Note that all the higher-order terms $\\Oof{\\Delta x^4}$\n", + "as it should. (Note that all the higher-order terms $\\mathcal{O}{\\Delta x^4}$\n", "also have higher-order derivatives that vanish for a cubic polynomial.)\n", "So to summarize, the fundamental problem is that $u$ as a product of\n", "a cubic polynomial and a linear or quadratic polynomial in time\n", @@ -4294,7 +4220,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.3" + "version": "3.8.2" } }, "nbformat": 4, diff --git a/fdm-devito-notebooks/02_wave/wave_app.ipynb b/fdm-devito-notebooks/02_wave/wave_app.ipynb index 23183ed0..b6d0dbd4 100644 --- a/fdm-devito-notebooks/02_wave/wave_app.ipynb +++ b/fdm-devito-notebooks/02_wave/wave_app.ipynb @@ -1884,7 +1884,25 @@ ] } ], - "metadata": {}, + "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.3" + } + }, "nbformat": 4, "nbformat_minor": 4 }