Add 2D wave solver and source wavelet utilities

Phase 2 Wave Equations progress:
- src/wave/wave2D_devito.py: 2D wave solver using Devito .laplace
- src/wave/sources.py: Ricker wavelet, Gaussian pulse, spectrum analysis
- Updated src/wave/__init__.py with new exports
- Added 16 new tests (29 total wave tests, 88 total)

Also includes typo fixes in legacy files from pre-commit hooks.

Author: ggorman

src/diffu/diffu2D_u0.py

@@ -82,7 +82,7 @@ def solver_dense(
     u_n = np.zeros((Nx+1, Ny+1))      # u at the previous time level
 
     Ix = range(0, Nx+1)
-    Iy = range(0, Ny+1)
+    It = range(0, Ny+1)
     It = range(0, Nt+1)
 
     # Make U_0x, U_0y, U_Lx and U_Ly functions if they are float/int
@@ -101,7 +101,7 @@ def solver_dense(
 
     # Load initial condition into u_n
     for i in Ix:
-        for j in Iy:
+        for j in It:
             u_n[i,j] = I(x[i], y[j])
 
     # Two-dim coordinate arrays for vectorized function evaluations
@@ -130,7 +130,7 @@ def solver_dense(
         p = m(i,j);  A[p, p] = 1
     # Loop over all internal mesh points in y diretion
     # and all mesh points in x direction
-    for j in Iy[1:-1]:
+    for j in It[1:-1]:
         i = 0;  p = m(i,j);  A[p, p] = 1   # boundary
         for i in Ix[1:-1]:                 # interior points
             p = m(i,j)
@@ -152,7 +152,7 @@ def solver_dense(
         j = 0
         for i in Ix:
             p = m(i,j);  b[p] = U_0y(t[n+1])  # boundary
-        for j in Iy[1:-1]:
+        for j in It[1:-1]:
             i = 0;  p = p = m(i,j);  b[p] = U_0x(t[n+1])  # boundary
             for i in Ix[1:-1]:
                 p = m(i,j)                                # interior
@@ -176,7 +176,7 @@ def solver_dense(
 
         # Fill u with vector c
         for i in Ix:
-            for j in Iy:
+            for j in It:
                 u[i,j] = c[m(i,j)]
 
         if user_action is not None:
@@ -228,7 +228,7 @@ def solver_sparse(
     u_n = np.zeros((Nx+1, Ny+1))    # u at the previous time level
 
     Ix = range(0, Nx+1)
-    Iy = range(0, Ny+1)
+    It = range(0, Ny+1)
     It = range(0, Nt+1)
 
     # Make U_0x, U_0y, U_Lx and U_Ly functions if they are float/int
@@ -247,7 +247,7 @@ def solver_sparse(
 
     # Load initial condition into u_n
     for i in Ix:
-        for j in Iy:
+        for j in It:
             u_n[i,j] = I(x[i], y[j])
 
     # Two-dim coordinate arrays for vectorized function evaluations
@@ -271,7 +271,7 @@ def solver_sparse(
 
     m = lambda i, j: j*(Nx+1) + i
     j = 0; main[m(0,j):m(Nx+1,j)] = 1  # j=0 boundary line
-    for j in Iy[1:-1]:             # Interior mesh lines j=1,...,Ny-1
+    for j in It[1:-1]:             # Interior mesh lines j=1,...,Ny-1
         i = 0;   main[m(i,j)] = 1  # Boundary
         i = Nx;  main[m(i,j)] = 1  # Boundary
         # Interior i points: i=1,...,N_x-1
@@ -306,7 +306,7 @@ def solver_sparse(
         j = 0
         for i in Ix:
             p = m(i,j);  b[p] = U_0y(t[n+1])          # Boundary
-        for j in Iy[1:-1]:
+        for j in It[1:-1]:
             i = 0;  p = m(i,j);  b[p] = U_0x(t[n+1])  # Boundary
             for i in Ix[1:-1]:
                 p = m(i,j)                            # Interior
@@ -329,7 +329,7 @@ def solver_sparse(
         f_a_n   = f(xv, yv, t[n])
 
         j = 0; b[m(0,j):m(Nx+1,j)] = U_0y(t[n+1])     # Boundary
-        for j in Iy[1:-1]:
+        for j in It[1:-1]:
             i = 0;   p = m(i,j);  b[p] = U_0x(t[n+1]) # Boundary
             i = Nx;  p = m(i,j);  b[p] = U_Lx(t[n+1]) # Boundary
             imin = Ix[1]
@@ -372,7 +372,7 @@ def CG_callback(c_k):
                       % (CG_iter[-1], CG_tol))
             '''
         # Fill u with vector c
-        #for j in Iy:  # vectorize y lines
+        #for j in It:  # vectorize y lines
         #    u[0:Nx+1,j] = c[m(0,j):m(Nx+1,j)]
         u[:,:] = c.reshape(Ny+1,Nx+1).T
 
@@ -443,7 +443,7 @@ def solver_classic_iterative(
         u_new = np.zeros((Nx+1, Ny+1))  # help array
 
     Ix = range(0, Nx+1)
-    Iy = range(0, Ny+1)
+    It = range(0, Ny+1)
     It = range(0, Nt+1)
 
     # Make U_0x, U_0y, U_Lx and U_Ly functions if they are float/int
@@ -462,7 +462,7 @@ def solver_classic_iterative(
 
     # Load initial condition into u_n
     for i in Ix:
-        for j in Iy:
+        for j in It:
             u_n[i,j] = I(x[i], y[j])
 
     # Two-dim coordinate arrays for vectorized function evaluations
@@ -488,7 +488,7 @@ def solver_classic_iterative(
                 j = 0
                 for i in Ix:
                     u[i,j] = U_0y(t[n+1])  # Boundary
-                for j in Iy[1:-1]:
+                for j in It[1:-1]:
                     i = 0;   u[i,j] = U_0x(t[n+1])  # Boundary
                     i = Nx;  u[i,j] = U_Lx(t[n+1])  # Boundary
                     for i in Ix[1:-1]:

src/diffu/diffu3D_u0.py

@@ -121,7 +121,7 @@ def solver_sparse_CG(
     u_n = np.zeros((Nx + 1, Ny + 1, Nz + 1))
 
     Ix = range(0, Nx + 1)
-    Iy = range(0, Ny + 1)
+    It = range(0, Ny + 1)
     Iz = range(0, Nz + 1)
     It = range(0, Nt + 1)
 
@@ -148,7 +148,7 @@ def solver_sparse_CG(
 
     # Load initial condition into u_n
     for i in Ix:
-        for j in Iy:
+        for j in It:
             for k in Iz:
                 u_n[i, j, k] = I(x[i], y[j], z[k])
 
@@ -181,7 +181,7 @@ def solver_sparse_CG(
     for k in Iz[1:-1]:  # interior mesh layers k=1,...,Nz-1
         j = 0
         main[m(0, j, k) : m(Nx + 1, j, k)] = 1  # j=0 boundary line
-        for j in Iy[1:-1]:  # interior mesh lines j=1,...,Ny-1
+        for j in It[1:-1]:  # interior mesh lines j=1,...,Ny-1
             i = 0
             main[m(i, j, k)] = 1  # boundary node
             i = Nx
@@ -225,7 +225,7 @@ def solver_sparse_CG(
         # Compute b, scalar version
         """
         k = 0                 # k=0 boundary layer
-        for j in Iy:
+        for j in It:
             for i in Ix:
                 p = m(i,j,k);  b[p] = U_0z(t[n+1])
 
@@ -234,7 +234,7 @@ def solver_sparse_CG(
             for i in Ix:
                 p = m(i,j,k);  b[p] = U_0y(t[n+1])
 
-            for j in Iy[1:-1]:   # interior mesh lines j=1,...,Ny-1
+            for j in It[1:-1]:   # interior mesh lines j=1,...,Ny-1
                 i = 0;  p = m(i,j,k);  b[p] = U_0x(t[n+1])  # boundary node
 
                 for i in Ix[1:-1]:           # interior nodes
@@ -253,7 +253,7 @@ def solver_sparse_CG(
                 p = m(i,j,k);  b[p] = U_Ly(t[n+1])
 
         k = Nz                 # k=Nz boundary layer
-        for j in Iy:
+        for j in It:
             for i in Ix:
                 p = m(i,j,k);  b[p] = U_Lz(t[n+1])
 
@@ -270,7 +270,7 @@ def solver_sparse_CG(
         for k in Iz[1:-1]:  # interior mesh layers k=1,...,Nz-1
             j = 0
             b[m(0, j, k) : m(Nx + 1, j, k)] = U_0y(t[n + 1])  # j=0, boundary mesh line
-            for j in Iy[1:-1]:  # interior mesh lines j=1,...,Ny-1
+            for j in It[1:-1]:  # interior mesh lines j=1,...,Ny-1
                 i = 0
                 p = m(i, j, k)
                 b[p] = U_0x(t[n + 1])  # boundary node
@@ -321,7 +321,7 @@ def solver_sparse_CG(
 
         # Fill u with vector c
         # for k in Iz:
-        # for j in Iy:
+        # for j in It:
         #    u[0:Nx+1,j,k] = c[m(0,j,k):m(Nx+1,j,k)]
         u[:, :, :] = c.reshape(Nz + 1, Ny + 1, Nx + 1).T
 

src/nonlin/split_diffu_react.py

@@ -229,7 +229,7 @@ def ordinary_splitting(I, a, b, f, L, dt, dt_Rfactor, F, t, T, user_action=None)
             user_action(u, x, t, n + 1)
 
 
-def Strang_splitting_1stOrder(I, a, b, f, L, dt, dt_Rfactor, F, t, T, user_action=None):
+def Strange_splitting_1stOrder(I, a, b, f, L, dt, dt_Rfactor, F, t, T, user_action=None):
     """Strange splitting while still using FE for the diffusion
     step and for the reaction step. Gives 1st order scheme.
     Introduce an extra time mesh t2 for the diffusion part,
@@ -293,7 +293,7 @@ def Strang_splitting_1stOrder(I, a, b, f, L, dt, dt_Rfactor, F, t, T, user_actio
             user_action(u, x, t, n + 1)
 
 
-def Strang_splitting_2ndOrder(I, a, b, f, L, dt, dt_Rfactor, F, t, T, user_action=None):
+def Strange_splitting_2andOrder(I, a, b, f, L, dt, dt_Rfactor, F, t, T, user_action=None):
     """Strange splitting using Crank-Nicolson for the diffusion
     step (theta-rule) and Adams-Bashforth 2 for the reaction step.
     Gives 2nd order scheme. Introduce an extra time mesh t2 for
@@ -436,9 +436,9 @@ def action(u, x, t, n):
                 T=T,
                 user_action=action,
             )
-        elif scheme == "Strang_splitting_1stOrder":
+        elif scheme == "Strange_splitting_1stOrder":
             print("Running Strange splitting with 1st order schemes...")
-            Strang_splitting_1stOrder(
+            Strange_splitting_1stOrder(
                 I=I,
                 a=a,
                 b=b,
@@ -451,9 +451,9 @@ def action(u, x, t, n):
                 T=T,
                 user_action=action,
             )
-        elif scheme == "Strang_splitting_2ndOrder":
+        elif scheme == "Strange_splitting_2andOrder":
             print("Running Strange splitting with 2nd order schemes...")
-            Strang_splitting_2ndOrder(
+            Strange_splitting_2andOrder(
                 I=I,
                 a=a,
                 b=b,
@@ -488,8 +488,8 @@ def action(u, x, t, n):
     schemes = [
         "diffusion",
         "ordinary_splitting",
-        "Strang_splitting_1stOrder",
-        "Strang_splitting_2ndOrder",
+        "Strange_splitting_1stOrder",
+        "Strange_splitting_2andOrder",
     ]
 
     for scheme in schemes:

src/softeng2/make_wave2D_u0.sh

@@ -60,4 +60,3 @@ else
   echo "Building Cython module $module failed"
   exit 1
 fi
-

src/softeng2/wave2D_u0_adv.py

@@ -105,14 +105,14 @@ def solver(I, V, f, c, Lx, Ly, Nx, Ny, dt, T,
     f_a = zeros((Nx+1,Ny+1), order=order)   # for compiled loops
 
     Ix = range(0, u.shape[0])
-    Iy = range(0, u.shape[1])
+    It = range(0, u.shape[1])
     It = range(0, t.shape[0])
 
     import time; t0 = time.perf_counter()          # for measuring CPU time
     # Load initial condition into u_n
     if version == 'scalar':
         for i in Ix:
-            for j in Iy:
+            for j in It:
                 u_n[i,j] = I(x[i], y[j])
     else: # use vectorized version
         u_n[:,:] = I(xv, yv)
@@ -173,30 +173,30 @@ def solver(I, V, f, c, Lx, Ly, Nx, Ny, dt, T,
 
 def advance_scalar(u, u_n, u_nm1, f, x, y, t, n, Cx2, Cy2, dt2,
                    V=None, step1=False):
-    Ix = range(0, u.shape[0]);  Iy = range(0, u.shape[1])
+    Ix = range(0, u.shape[0]);  It = range(0, u.shape[1])
     if step1:
         dt = sqrt(dt2)  # save
         Cx2 = 0.5*Cx2;  Cy2 = 0.5*Cy2; dt2 = 0.5*dt2  # redefine
         D1 = 1;  D2 = 0
     else:
         D1 = 2;  D2 = 1
     for i in Ix[1:-1]:
-        for j in Iy[1:-1]:
+        for j in It[1:-1]:
             u_xx = u_n[i-1,j] - 2*u_n[i,j] + u_n[i+1,j]
             u_yy = u_n[i,j-1] - 2*u_n[i,j] + u_n[i,j+1]
             u[i,j] = D1*u_n[i,j] - D2*u_nm1[i,j] + \
                      Cx2*u_xx + Cy2*u_yy + dt2*f(x[i], y[j], t[n])
             if step1:
                 u[i,j] += dt*V(x[i], y[j])
     # Boundary condition u=0
-    j = Iy[0]
+    j = It[0]
     for i in Ix: u[i,j] = 0
-    j = Iy[-1]
+    j = It[-1]
     for i in Ix: u[i,j] = 0
     i = Ix[0]
-    for j in Iy: u[i,j] = 0
+    for j in It: u[i,j] = 0
     i = Ix[-1]
-    for j in Iy: u[i,j] = 0
+    for j in It: u[i,j] = 0
     return u
 
 def advance_vectorized(u, u_n, u_nm1, f_a, Cx2, Cy2, dt2,

src/softeng2/wave2D_u0_class.py

@@ -89,14 +89,14 @@ def solver(I, V, f, c, Lx, Ly, Nx, Ny, dt, T,
     f_a = zeros((Nx+1,Ny+1), order=order)   # for compiled loops
 
     Ix = range(0, u.shape[0])
-    Iy = range(0, u.shape[1])
+    It = range(0, u.shape[1])
     It = range(0, t.shape[0])
 
     import time; t0 = time.perf_counter()          # for measuring CPU time
     # Load initial condition into u_1
     if version == 'scalar':
         for i in Ix:
-            for j in Iy:
+            for j in It:
                 u_1[i,j] = I(x[i], y[j])
     else: # use vectorized version
         u_1[:,:] = I(xv, yv)
@@ -158,30 +158,30 @@ def solver(I, V, f, c, Lx, Ly, Nx, Ny, dt, T,
 
 def advance_scalar(u, u_1, u_2, f, x, y, t, n, Cx2, Cy2, dt2,
                    V=None, step1=False):
-    Ix = range(0, u.shape[0]);  Iy = range(0, u.shape[1])
+    Ix = range(0, u.shape[0]);  It = range(0, u.shape[1])
     if step1:
         dt = sqrt(dt2)  # save
         Cx2 = 0.5*Cx2;  Cy2 = 0.5*Cy2; dt2 = 0.5*dt2  # redefine
         D1 = 1;  D2 = 0
     else:
         D1 = 2;  D2 = 1
     for i in Ix[1:-1]:
-        for j in Iy[1:-1]:
+        for j in It[1:-1]:
             u_xx = u_1[i-1,j] - 2*u_1[i,j] + u_1[i+1,j]
             u_yy = u_1[i,j-1] - 2*u_1[i,j] + u_1[i,j+1]
             u[i,j] = D1*u_1[i,j] - D2*u_2[i,j] + \
                      Cx2*u_xx + Cy2*u_yy + dt2*f(x[i], y[j], t[n])
             if step1:
                 u[i,j] += dt*V(x[i], y[j])
     # Boundary condition u=0
-    j = Iy[0]
+    j = It[0]
     for i in Ix: u[i,j] = 0
-    j = Iy[-1]
+    j = It[-1]
     for i in Ix: u[i,j] = 0
     i = Ix[0]
-    for j in Iy: u[i,j] = 0
+    for j in It: u[i,j] = 0
     i = Ix[-1]
-    for j in Iy: u[i,j] = 0
+    for j in It: u[i,j] = 0
     return u
 
 def advance_vectorized(u, u_1, u_2, f_a, Cx2, Cy2, dt2,

src/softeng2/wave2D_u0_loop_c.c

@@ -20,4 +20,3 @@ void advance(double* u, double* u_1, double* u_2, double* f,
   i = 0;  for (j=0; j<=Ny; j++) u[idx(i,j)] = 0;
   i = Nx; for (j=0; j<=Ny; j++) u[idx(i,j)] = 0;
 }
-

src/softeng2/wave2D_u0_loop_c_f2py_signature.f

@@ -7,5 +7,3 @@ subroutine advance(u, u_1, u_2, f, Cx2, Cy2, dt2, Nx, Ny)
 Cf2py intent(c) u, u_1, u_2, f, Cx2, Cy2, dt2, Nx, Ny
       return
       end
-
-