comparison between numerical and analytical solution of the wave equation

mms/acoustic_2D_mms.py

@@ -0,0 +1,203 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from simwave import (
+    SpaceModel, TimeModel, RickerWavelet, Wavelet, Solver, Compiler,
+    Receiver, Source, plot_wavefield, plot_shotrecord, plot_velocity_model
+)
+
+
+def initial_shape(space_model, kx, kz):
+    nbl_pad_width = space_model.nbl_pad_width
+    halo_pad_width = space_model.halo_pad_width
+    bbox = space_model.bounding_box
+    # values
+    xmin, xmax = bbox[0: 2]
+    zmin, zmax = bbox[2: 4]
+    x_coords = np.linspace(xmin, xmax, space_model.shape[0])
+    z_coords = np.linspace(zmin, zmax, space_model.shape[1])
+    X, Z = np.meshgrid(x_coords, z_coords)
+
+    return np.sin(kx * X) * np.sin(kz * Z)
+
+
+class MySource(Source):
+
+    @property
+    def interpolated_points_and_values(self):
+        """ Calculate the amplitude value at all points"""
+
+        points = np.array([], dtype=np.uint)
+        values = np.array([], dtype=self.space_model.dtype)
+
+        nbl_pad_width = self.space_model.nbl_pad_width
+        halo_pad_width = self.space_model.halo_pad_width
+        dimension = self.space_model.dimension
+        bbox = self.space_model.bounding_box
+
+        for dim, dim_length in enumerate(self.space_model.shape):
+
+            # points
+            lpad = halo_pad_width[dim][0] + nbl_pad_width[dim][0]
+            rpad = halo_pad_width[dim][1] + nbl_pad_width[dim][1] + dim_length
+            points = np.append(points, np.array([lpad, rpad], dtype=np.uint))
+            # values
+            coor_min, coor_max = bbox[dim * dimension:2 + dim * dimension]
+            coordinates = np.linspace(coor_min, coor_max, dim_length)
+            amplitudes = np.sin(np.pi / (coor_max - coor_min) * coordinates)
+            values = np.append(values, amplitudes.astype(values.dtype))
+            # offset
+            offset = np.array([0, values.size], dtype=np.uint)
+
+        return points, values, offset
+
+def cosenoid(time_model, amplitude, omega=2*np.pi*1):
+    """Function that generates the signal satisfying s(0) = ds/dt(0) = 0"""
+    return amplitude * np.cos(omega * time_model.time_values)
+
+def create_models(Lx, Lz, vel, tf, h, p):
+    # create the space model
+    space_model = SpaceModel(
+        bounding_box=(0, Lx, 0, Lz),
+        grid_spacing=(h, h),
+        velocity_model=vel,
+        space_order=p,
+        dtype=np.float32
+    )
+
+    # config boundary conditions
+    # (none,  null_dirichlet or null_neumann)
+    space_model.config_boundary(
+        damping_length=0,
+        boundary_condition=(
+            "null_dirichlet", "null_dirichlet",
+            "null_dirichlet", "null_dirichlet"
+        ),
+        damping_polynomial_degree=3,
+        damping_alpha=0.001
+    ) 
+
+    # create the time model
+    time_model = TimeModel(
+        space_model=space_model,
+        tf=tf,
+        saving_stride=1
+    )
+
+    return space_model, time_model
+
+def geometry_acquisition(space_model, sources=None, receivers=None):
+    if sources is None:
+        sources=[(2560, 2560)]
+    if receivers is None:
+        receivers=[(2560, i) for i in range(0, 5120, 10)]
+
+    # create the set of sources
+    # source = Source(
+    #     space_model,
+    #     coordinates=sources,
+    #     window_radius=4
+    # )
+
+    # personalized source
+    my_source = MySource(space_model, coordinates=[])
+
+    # crete the set of receivers
+    receiver = Receiver(
+        space_model=space_model,
+        coordinates=receivers,
+        window_radius=4
+    )
+
+    return my_source, receiver
+
+def build_solver(space_model, time_model, freq, vp, kx, kz, omega):
+
+    # geometry acquisition
+    source, receiver = geometry_acquisition(space_model)
+
+    # create a ricker wavelet with 10hz of peak frequency
+    # ricker = RickerWavelet(freq, time_model)
+
+    # create a cosenoid wavelet
+    amplitude = (kx**2+kz**2) - omega**2/vp**2
+    # amplitude = 0
+    wavelet = Wavelet(cosenoid, time_model=time_model, amplitude=amplitude, omega=omega)
+
+    # set compiler options
+    compiler = Compiler( cc='gcc', language='cpu_openmp', cflags='-O3 -fPIC -ffast-math -Wall -std=c99 -shared')
+
+    # create the solver
+    solver = Solver(
+        space_model=space_model,
+        time_model=time_model,
+        sources=source,
+        receivers=receiver,
+        # wavelet=ricker,
+        wavelet=wavelet,
+        compiler=compiler
+    )
+
+    return solver
+
+
+if __name__ == "__main__":
+
+    # problem params
+    Lx, Lz = 5120, 5120
+    tf = 1.5
+    freq = 5
+    vp = 1500
+
+    # discretization params
+    n = 513    
+    h = 10
+    p = 4
+
+    # harmonic parameters
+    kx = np.pi / Lx
+    kz = np.pi / Lz
+    omega = 2*np.pi*freq
+
+    # velocity model
+    vel = vp * np.ones(shape=(n, n), dtype=np.float32)
+
+    # discretization
+    space_model, time_model = create_models(Lx, Lz, vel, tf, h, p)
+
+    # initial grid
+    # initial_grid = np.zeros(shape=space_model.shape)
+    initial_grid = initial_shape(space_model, kx, kz)
+
+    # get solver
+    solver = build_solver(space_model, time_model, freq, vp, kx, kz, omega)
+
+    # run the forward
+    u_full, recv = solver.forward(
+        [initial_grid * np.cos(-1 * time_model.dt * omega),
+         initial_grid * np.cos(-0 * time_model.dt * omega)
+        ]
+    )
+
+    # plot stuff
+    print("u_full shape:", u_full.shape)
+    plot_wavefield(u_full[-1])
+    plot_shotrecord(recv)
+
+    # save numpy array
+    # print("saving numpy array")
+    # np.save("numerical_harmonic", u_full)    
+
+    # load analytical array
+    u_analytic = np.load("analytical_harmonic.npy")
+    # plot comparison
+    numerical_values = u_full[:, 256, 256]
+    analytic_values = u_analytic[:, 256, 256]
+    time_values = time_model.time_values
+    plt.plot(time_values, analytic_values, time_values, numerical_values)
+    plt.legend(["analytical", "numerical"])
+    plt.title("Comparison between analytical and numerical result (simwave)")
+    plt.xlabel("time")
+    plt.show()
+
+

mms/comparison.sh

@@ -0,0 +1,5 @@
+#!/bin/bash
+
+# create analytical solution, then numerical and compare both
+python harmonic_solution.py
+python acoustic_2D_mms.py

mms/harmonic_solution.py

@@ -0,0 +1,46 @@
+"""Generates an harmonic solution to the wave equation
+
+div(grad(u)) - 1/c^2 * d^2u/dt^2 = f
+
+where u = sin(kx*x) * sin(kz*z) * cos(w*t)
+"""
+
+
+# from simwave import plot_shotrecord
+import numpy as np
+import matplotlib.pyplot as plt
+# fig = plt.figure()
+# ax = fig.add_subplot(projection='3d')
+# ax.set_zlim(-1.01, 1.01)
+
+# Domain info
+dz, dx = (10, 10)
+zmin, zmax, xmin, xmax = (0, 5120, 0, 5120)
+nz = int((zmax - zmin) / dz + 1)
+nx = int((xmax - xmin) / dx + 1)
+zcoords = np.linspace(zmin, zmax, nz)
+xcoords = np.linspace(xmin, xmax, nx)
+Z, X = np.meshgrid(zcoords, xcoords)
+
+# Params
+c = 1500
+freq= 5
+w = 2*np.pi*freq
+kz = np.pi / (zmax - zmin)
+kx = np.pi / (xmax - xmin)
+
+# Time info
+t0 = 0
+tf = 1.5
+nt = 369
+tp = np.linspace(t0, tf, nt)
+
+# Solution
+space_solution = np.sin(kx*X) * np.sin(kz*Z)
+time_solution = np.cos(w*tp)
+solution = (space_solution[:,:,None] * time_solution).T
+
+# Save it
+print("saving analytical solution to .npy file")
+np.save("analytical_harmonic", solution)
+

simwave/kernel/frontend/solver.py

@@ -132,11 +132,11 @@ def forward(self, initial_grid=None):
 
             initial_grid = self.space_model.dtype(initial_grid)
 
-            if initial_grid.shape != self.space_model.shape:
-                raise ValueError(
-                    "initial_grid must have shape = {}").format(
-                    self.space_model.shape
-                )
+            # if initial_grid.shape != self.space_model.shape:
+            #     raise ValueError(
+            #         "initial_grid must have shape = {}").format(
+            #         self.space_model.shape
+            #     )
 
             halo = self.space_model.halo_size[0]
             nbl = self.space_model.nbl
@@ -152,8 +152,8 @@ def forward(self, initial_grid=None):
                 x1 = halo + nbl[2]
                 x2 = halo + nbl[3]
 
-                u[0, z1:-z2, x1:-x2] = initial_grid
-                u[1, z1:-z2, x1:-x2] = initial_grid
+                u[0, z1:-z2, x1:-x2] = initial_grid[0]
+                u[1, z1:-z2, x1:-x2] = initial_grid[1]
 
             else:
 
@@ -164,8 +164,8 @@ def forward(self, initial_grid=None):
                 y1 = halo + nbl[4]
                 y2 = halo + nbl[5]
 
-                u[0, z1:-z2, x1:-x2, y1:-y2] = initial_grid
-                u[1, z1:-z2, x1:-x2, y1:-y2] = initial_grid
+                u[0, z1:-z2, x1:-x2, y1:-y2] = initial_grid[0]
+                u[1, z1:-z2, x1:-x2, y1:-y2] = initial_grid[1]
 
         u_full, recv = self._middleware.exec(
             operator='forward',