bench: Add wave3d setup

fast/setup_wave3d.py

@@ -0,0 +1,115 @@
+# Script to save initial data for the Acoustic wave execution benchmark
+# Based on the implementation of the Devito acoustic example implementation
+# Not using Devito's source injection abstraction
+import sys
+import numpy as np
+
+from devito import (TimeFunction, Eq, Operator, solve, norm,
+                    configuration)
+from examples.seismic import RickerSource
+from examples.seismic import Model, TimeAxis
+from fast.bench_utils import plot_3dfunc
+from devito.tools import as_tuple
+
+import argparse
+np.set_printoptions(threshold=np.inf)
+
+
+parser = argparse.ArgumentParser(description='Process arguments.')
+
+parser.add_argument("-d", "--shape", default=(16, 16, 16), type=int, nargs="+",
+                    help="Number of grid points along each axis")
+parser.add_argument("-so", "--space_order", default=4,
+                    type=int, help="Space order of the simulation")
+parser.add_argument("-to", "--time_order", default=2,
+                    type=int, help="Time order of the simulation")
+parser.add_argument("-nt", "--nt", default=20,
+                    type=int, help="Simulation time in millisecond")
+parser.add_argument("-bls", "--blevels", default=1, type=int, nargs="+",
+                    help="Block levels")
+parser.add_argument("-plot", "--plot", default=False, type=bool, help="Plot2D")
+parser.add_argument("-devito", "--devito", default=False, type=bool, help="Devito run")
+parser.add_argument("-xdsl", "--xdsl", default=False, type=bool, help="xDSL run")
+args = parser.parse_args()
+
+
+mpiconf = configuration['mpi']
+
+# Define a physical size
+# nx, ny, nz = args.shape
+nt = args.nt
+
+shape = (args.shape)  # Number of grid point (nx, ny, nz)
+spacing = as_tuple(10.0 for _ in range(len(shape)))  # Grid spacing in m. The domain size is now 1km by 1km
+origin = as_tuple(0.0 for _ in range(len(shape)))  # What is the location of the top left corner.
+# This is necessary to define
+# the absolute location of the source and receivers
+
+# Define a velocity profile. The velocity is in km/s
+v = np.empty(shape, dtype=np.float32)
+v[:, :, :] = 1
+
+# With the velocity and model size defined, we can create the seismic model that
+# encapsulates this properties. We also define the size of the absorbing layer as
+# 10 grid points
+so = args.space_order
+to = args.time_order
+
+model = Model(vp=v, origin=origin, shape=shape, spacing=spacing,
+              space_order=so, nbl=0)
+
+# plot_velocity(model)
+
+t0 = 0.  # Simulation starts a t=0
+tn = nt  # Simulation last 1 second (1000 ms)
+dt = model.critical_dt  # Time step from model grid spacing
+print("dt is:", dt)
+
+time_range = TimeAxis(start=t0, stop=tn, step=dt)
+
+# The source is positioned at a $20m$ depth and at the middle of the
+# $x$ axis ($x_{src}=500m$),
+# with a peak wavelet frequency of $10Hz$.
+f0 = 0.010  # Source peak frequency is 10Hz (0.010 kHz)
+src = RickerSource(name='src', grid=model.grid, f0=f0,
+                   npoint=1, time_range=time_range)
+
+# First, position source centrally in all dimensions, then set depth
+src.coordinates.data[0, :] = np.array(model.domain_size) * .5
+
+# We can plot the time signature to see the wavelet
+# src.show()
+
+# Define the wavefield with the size of the model and the time dimension
+u = TimeFunction(name="u", grid=model.grid, time_order=to, space_order=so)
+# Another one to clone data
+u2 = TimeFunction(name="u", grid=model.grid, time_order=to, space_order=so)
+ub = TimeFunction(name="ub", grid=model.grid, time_order=to, space_order=so)
+
+# We can now write the PDE
+# pde = model.m * u.dt2 - u.laplace + model.damp * u.dt
+pde = u.dt2 - u.laplace
+
+stencil = Eq(u.forward, solve(pde, u.forward))
+# stencil
+
+# Finally we define the source injection and receiver read function to generate
+# the corresponding code
+# print(time_range)
+
+print("Init norm:", np.linalg.norm(u.data[:]))
+src_term = src.inject(field=u.forward, expr=src * dt**2 / model.m)
+op0 = Operator([stencil] + src_term, subs=model.spacing_map, name='SourceDevitoOperator')
+
+# Run with source and plot
+op0.apply(time=time_range.num-1, dt=model.critical_dt)
+
+if len(shape) == 3:
+    if args.plot:
+        plot_3dfunc(u)
+
+# Save Data here
+shape_str = '_'.join(str(item) for item in shape)
+np.save("so%s_critical_dt%s.npy" % (so, shape_str), model.critical_dt, allow_pickle=True)
+np.save("so%s_wave_dat%s.npy" % (so, shape_str), u.data[:], allow_pickle=True)
+np.save("so%s_grid_extent%s.npy" % (so, shape_str), model.grid.extent, allow_pickle=True)

fast/wave2d_b.py

@@ -37,7 +37,7 @@
 # Define a physical size
 # nx, ny, nz = args.shape
 nt = args.nt
-so = args.so
+so = args.space_order
 
 shape = (args.shape)  # Number of grid point (nx, ny, nz)
 shape_str = '_'.join(str(item) for item in shape)

fast/wave3d_b.py

@@ -0,0 +1,129 @@
+# Based on the implementation of the Devito acoustic example implementation
+# Not using Devito's source injection abstraction
+import sys
+import numpy as np
+
+from devito import (TimeFunction, Eq, Operator, solve, norm,
+                    XDSLOperator, configuration, Grid)
+from examples.seismic import RickerSource
+from examples.seismic import Model, TimeAxis, plot_image
+from fast.bench_utils import plot_3dfunc
+from devito.tools import as_tuple
+
+import argparse
+np.set_printoptions(threshold=np.inf)
+
+
+parser = argparse.ArgumentParser(description='Process arguments.')
+
+parser.add_argument("-d", "--shape", default=(16, 16, 16), type=int, nargs="+",
+                    help="Number of grid points along each axis")
+parser.add_argument("-so", "--space_order", default=4,
+                    type=int, help="Space order of the simulation")
+parser.add_argument("-to", "--time_order", default=2,
+                    type=int, help="Time order of the simulation")
+parser.add_argument("-nt", "--nt", default=20,
+                    type=int, help="Simulation time in millisecond")
+parser.add_argument("-bls", "--blevels", default=1, type=int, nargs="+",
+                    help="Block levels")
+parser.add_argument("-plot", "--plot", default=False, type=bool, help="Plot2D")
+parser.add_argument("-devito", "--devito", default=False, type=bool, help="Devito run")
+parser.add_argument("-xdsl", "--xdsl", default=False, type=bool, help="xDSL run")
+args = parser.parse_args()
+
+
+mpiconf = configuration['mpi']
+
+# Define a physical size
+# nx, ny, nz = args.shape
+nt = args.nt
+so = args.space_order
+
+shape = (args.shape)  # Number of grid point (nx, ny, nz)
+shape_str = '_'.join(str(item) for item in shape)
+spacing = as_tuple(10.0 for _ in range(len(shape)))  # Grid spacing in m. The domain size is now 1km by 1km
+origin = as_tuple(0.0 for _ in range(len(shape)))  # What is the location of the top left corner.
+domain_size = tuple((d-1) * s for d, s in zip(shape, spacing))
+extent = np.load("so%s_grid_extent%s.npy" % (so, shape_str), allow_pickle=True)
+
+grid = Grid(shape=shape, extent=as_tuple(extent))
+
+# With the velocity and model size defined, we can create the seismic model that
+# encapsulates this properties. We also define the size of the absorbing layer as
+# 10 grid points
+so = args.space_order
+to = args.time_order
+
+t0 = 0.  # Simulation starts a t=0
+tn = nt  # Simulation last 1 second (1000 ms)
+
+# Define the wavefield with the size of the model and the time dimension
+u = TimeFunction(name="u", grid=grid, time_order=to, space_order=so)
+# Another one to clone data
+u2 = TimeFunction(name="u", grid=grid, time_order=to, space_order=so)
+ub = TimeFunction(name="ub", grid=grid, time_order=to, space_order=so)
+
+# We can now write the PDE
+# pde = model.m * u.dt2 - u.laplace + model.damp * u.dt
+# import pdb;pdb.set_trace()
+pde = u.dt2 - u.laplace
+
+stencil = Eq(u.forward, solve(pde, u.forward))
+
+# print("Init Devito linalg norm 0 :", np.linalg.norm(u.data[0]))
+# print("Init Devito linalg norm 1 :", np.linalg.norm(u.data[1]))
+# print("Init Devito linalg norm 2 :", np.linalg.norm(u.data[2]))
+# print("Norm of initial data:", norm(u))
+
+configuration['mpi'] = 0
+u2.data[:] = u.data[:]
+configuration['mpi'] = mpiconf
+
+u.data[:] = np.load("so%s_wave_dat%s.npy" % (so, shape_str), allow_pickle=True)
+dt = np.load("so%s_critical_dt%s.npy" % (so, shape_str), allow_pickle=True)
+
+# np.save("critical_dt%s.npy" % shape_str, model.critical_dt, allow_pickle=True)
+# np.save("wave_dat%s.npy" % shape_str, u.data[:], allow_pickle=True)
+
+if len(shape) == 3 and args.plot:
+    plot_3dfunc(u)
+
+print("Init norm:", np.linalg.norm(u.data[:]))
+# print("Init linalg norm:", np.linalg.norm(u.data[0]))
+# print("Init linalg norm:", np.linalg.norm(u.data[1]))
+# print("Init linalg norm:", np.linalg.norm(u.data[2]))
+
+
+if args.devito:
+    # Run more with no sources now (Not supported in xdsl)
+    # op1 = Operator([stencil], name='DevitoOperator', subs=grid.spacing_map)
+    op1 = Operator([stencil], name='DevitoOperator')
+    op1.apply(time=nt, dt=dt)
+
+    configuration['mpi'] = 0
+    ub.data[:] = u.data[:]
+    configuration['mpi'] = mpiconf
+
+    if len(shape) == 3 and args.plot:
+        plot_3dfunc(u)
+
+    print("Devito norm:", norm(u))
+    # print("Devito linalg norm 0:", np.linalg.norm(u.data[0]))
+    # print("Devito linalg norm 1:", np.linalg.norm(u.data[1]))
+    # print("Devito linalg norm 2:", np.linalg.norm(u.data[2]))
+
+
+if args.xdsl:
+
+    # Run more with no sources now (Not supported in xdsl)
+    xdslop = XDSLOperator([stencil], name='xDSLOperator')
+    xdslop.apply(time=nt, dt=dt)
+
+    if len(shape) == 3 and args.plot:
+        plot_3dfunc(u)
+
+    print("XDSL norm:", norm(u))
+
+    # print("XDSL output norm 0:", np.linalg.norm(u.data[0]))
+    # print("XDSL output norm 1:", np.linalg.norm(u.data[1]))
+    # print("XDSL output norm 2:", np.linalg.norm(u.data[2]))