From 38bf444b45a75eee9c2b74c78586d490886bc81b Mon Sep 17 00:00:00 2001 From: Emilien Bauer Date: Fri, 11 Aug 2023 19:22:37 +0100 Subject: [PATCH] Get scripts WGPU+MacOS ready --- fast/diffusion_2D_wBCs.py | 6 +- fast/wave2d.py | 158 -------------------------------------- fast/wave2d_b.py | 20 +++++ fast/wave3d.py | 148 ----------------------------------- fast/wave3d_b.py | 24 +++++- 5 files changed, 44 insertions(+), 312 deletions(-) delete mode 100644 fast/wave2d.py delete mode 100644 fast/wave3d.py diff --git a/fast/diffusion_2D_wBCs.py b/fast/diffusion_2D_wBCs.py index d3c0816cb4..665ba20720 100644 --- a/fast/diffusion_2D_wBCs.py +++ b/fast/diffusion_2D_wBCs.py @@ -58,7 +58,7 @@ if args.devito: op = Operator([eq_stencil], name='DevitoOperator', opt=('advanced', {'par-tile': (32,4,8)})) op.apply(time=nt, dt=dt, a=nu) - print("Devito Field norm is:", norm(u)) + print("Devito Field norm is:", np.linalg.norm(u.data)) if args.plot: plot_2dfunc(u) @@ -66,7 +66,7 @@ if args.wgpu: op = WGPUOperator([eq_stencil], name='WGPUOperator', opt=('advanced', {'par-tile': (32,4,8)})) op.apply(time=nt, dt=dt, a=nu) - print("WGPU Field norm is:", norm(u)) + print("WGPU Field norm is:", np.linalg.norm(u.data)) if args.plot: plot_2dfunc(u) @@ -77,4 +77,4 @@ if args.xdsl: xdslop = XDSLOperator([eq_stencil], name='XDSLOperator') xdslop.apply(time=nt, dt=dt, a=nu) - print("XDSL Field norm is:", norm(u)) + print("XDSL Field norm is:", np.linalg.norm(u.data)) diff --git a/fast/wave2d.py b/fast/wave2d.py deleted file mode 100644 index 3c60dc7c3b..0000000000 --- a/fast/wave2d.py +++ /dev/null @@ -1,158 +0,0 @@ -# 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) -from examples.seismic import RickerSource -from examples.seismic import Model, TimeAxis, plot_image -from fast.bench_utils import plot_2dfunc -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), 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 -# import pdb;pdb.set_trace() -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) == 2: - if args.plot: - plot_2dfunc(u) - -# 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 - -if args.devito: - # Run more with no sources now (Not supported in xdsl) - op1 = Operator([stencil], name='DevitoOperator', opt=('advanced', {'par-tile': (32,4,8)})) - op1.apply(time=time_range.num-1, dt=model.critical_dt) - - configuration['mpi'] = 0 - ub.data[:] = u.data[:] - configuration['mpi'] = mpiconf - - if len(shape) == 2 and args.plot: - plot_2dfunc(u) - - # print("After Operator 1: 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: - # Reset initial data - configuration['mpi'] = 0 - u.data[:] = u2.data[:] - configuration['mpi'] = mpiconf - # v[:, ..., :] = 1 - # print("Reinitialise data: Devito norm:", norm(u)) - # print("XDSL init linalg norm:", np.linalg.norm(u.data[0])) - # print("XDSL init linalg norm:", np.linalg.norm(u.data[1])) - # print("XDSL init linalg norm:", np.linalg.norm(u.data[2])) - - # Run more with no sources now (Not supported in xdsl) - xdslop = XDSLOperator([stencil], name='xDSLOperator') - xdslop.apply(time=time_range.num-1, dt=model.critical_dt) - - if len(shape) == 2 and args.plot: - plot_2dfunc(u) - - print("XDSL output norm 0:", np.linalg.norm(u.data[0]), "vs:", np.linalg.norm(ub.data[0])) - print("XDSL output norm 1:", np.linalg.norm(u.data[1]), "vs:", np.linalg.norm(ub.data[1])) - print("XDSL output norm 2:", np.linalg.norm(u.data[2]), "vs:", np.linalg.norm(ub.data[2])) diff --git a/fast/wave2d_b.py b/fast/wave2d_b.py index 411d61ef14..f5fbdbf468 100644 --- a/fast/wave2d_b.py +++ b/fast/wave2d_b.py @@ -5,6 +5,7 @@ from devito import (TimeFunction, Eq, Operator, solve, norm, XDSLOperator, configuration, Grid) +from devito.operator.wgpu_operator import WGPUOperator from examples.seismic import RickerSource from examples.seismic import Model, TimeAxis, plot_image from fast.bench_utils import plot_2dfunc @@ -29,6 +30,7 @@ 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") +parser.add_argument("-wgpu", "--wgpu", default=False, type=bool, help="xDSL run") args = parser.parse_args() @@ -108,6 +110,24 @@ # print("Devito linalg norm 1:", np.linalg.norm(u.data[1])) # print("Devito linalg norm 2:", np.linalg.norm(u.data[2])) +if args.wgpu: + # Run more with no sources now (Not supported in xdsl) + # op1 = Operator([stencil], name='DevitoOperator', subs=grid.spacing_map) + op1 = WGPUOperator([stencil], name='WGPUOperator') + op1.apply(time=nt, dt=dt) + + configuration['mpi'] = 0 + ub.data[:] = u.data[:] + configuration['mpi'] = mpiconf + + if len(shape) == 2 and args.plot: + plot_2dfunc(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: # print("Reinitialise data: Devito norm:", norm(u)) diff --git a/fast/wave3d.py b/fast/wave3d.py deleted file mode 100644 index 71602ca30f..0000000000 --- a/fast/wave3d.py +++ /dev/null @@ -1,148 +0,0 @@ -# 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) -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=(11, 11, 11), 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=200, - type=int, help="Simulation time in millisecond") -parser.add_argument("-bls", "--blevels", default=2, type=int, nargs="+", - help="Block levels") -parser.add_argument("-plot", "--plot", default=False, type=bool, help="Plot3D") -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 - -# 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 -# import pdb;pdb.set_trace() -pde = u.dt2 - u.laplace - -stencil = Eq(u.forward, solve(pde, u.forward)) - -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) - -# devito_norm = norm(u) -# print("Init linalg norm 0 (inlined) :", norm(u)) -# print("Init linalg norm 0 :", np.linalg.norm(u.data[0])) -# print("Init linalg norm 1 :", np.linalg.norm(u.data[1])) -# print("Init linalg norm 2 :", np.linalg.norm(u.data[2])) -# print("Norm of initial data:", np.linalg.norm(u.data[:])) - -configuration['mpi'] = 0 -u2.data[:] = u.data[:] -configuration['mpi'] = mpiconf - -# Run more with no sources now (Not supported in xdsl) -op1 = Operator([stencil], name='DevitoOperator', opt=('advanced', {'par-tile': (32,4,8)})) -op1.apply(time=time_range.num-1, dt=model.critical_dt) - -configuration['mpi'] = 0 -ub.data[:] = u.data[:] -configuration['mpi'] = mpiconf - -if len(shape) == 3: - if args.plot: - plot_3dfunc(u) - -# print("After Operator 1: Devito norm:", np.linalg.norm(u.data[:])) -#print("Devito norm 0:", np.linalg.norm(u.data[0])) -#print("Devito norm 1:", np.linalg.norm(u.data[1])) -#print("Devito norm 2:", np.linalg.norm(u.data[2])) - -# Reset initial data -configuration['mpi'] = 0 -u.data[:] = u2.data[:] -configuration['mpi'] = mpiconf - -# print("Reinitialise data for XDSL:", np.linalg.norm(u.data[:])) -# print("Init XDSL linalg norm 0:", np.linalg.norm(u.data[0])) -# print("Init XDSL linalg norm 1:", np.linalg.norm(u.data[1])) -# print("Init XDSL linalg norm 2:", np.linalg.norm(u.data[2])) - -# Run more with no sources now (Not supported in xdsl) -xdslop = XDSLOperator([stencil], name='XDSLOperator') -xdslop.apply(time=time_range.num-1, dt=model.critical_dt) - -print("XDSL output norm 0:", np.linalg.norm(u.data[0]), "vs:", np.linalg.norm(ub.data[0])) -print("XDSL output norm 1:", np.linalg.norm(u.data[1]), "vs:", np.linalg.norm(ub.data[1])) -print("XDSL output norm 2:", np.linalg.norm(u.data[2]), "vs:", np.linalg.norm(ub.data[2])) diff --git a/fast/wave3d_b.py b/fast/wave3d_b.py index 3e3081f597..89f726f34d 100644 --- a/fast/wave3d_b.py +++ b/fast/wave3d_b.py @@ -5,6 +5,7 @@ from devito import (TimeFunction, Eq, Operator, solve, norm, XDSLOperator, configuration, Grid) +from devito.operator.wgpu_operator import WGPUOperator from examples.seismic import RickerSource from examples.seismic import Model, TimeAxis, plot_image from fast.bench_utils import plot_3dfunc @@ -29,6 +30,7 @@ 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") +parser.add_argument("-wgpu", "--wgpu", default=False, type=bool, help="xDSL run") args = parser.parse_args() @@ -72,7 +74,7 @@ # 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)) +# print("Norm of initial data:", np.linalg.norm(u.data)) configuration['mpi'] = 0 u2.data[:] = u.data[:] @@ -106,12 +108,28 @@ if len(shape) == 3 and args.plot: plot_3dfunc(u) - print("Devito norm:", norm(u)) + print("Devito norm:", np.linalg.norm(u.data)) # 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.wgpu: + + # Run more with no sources now (Not supported in xdsl) + xdslop = WGPUOperator([stencil], name='WGPUOperator') + xdslop.apply(time=nt, dt=dt) + + if len(shape) == 3 and args.plot: + plot_3dfunc(u) + + print("XDSL norm:", np.linalg.norm(u.data)) + + # 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])) + + if args.xdsl: # Run more with no sources now (Not supported in xdsl) @@ -121,7 +139,7 @@ if len(shape) == 3 and args.plot: plot_3dfunc(u) - print("XDSL norm:", norm(u)) + print("XDSL norm:", np.linalg.norm(u.data)) # print("XDSL output norm 0:", np.linalg.norm(u.data[0])) # print("XDSL output norm 1:", np.linalg.norm(u.data[1]))