Merge pull request #34 from xdslproject/emilien/correctness

Put back correctness checks in every bench file.

fast/diffusion_2D_wBCs.py

@@ -5,7 +5,7 @@
 import argparse
 import numpy as np
 
-from devito import Grid, TimeFunction, Eq, solve, Operator, Constant, norm, XDSLOperator
+from devito import Grid, TimeFunction, Eq, solve, Operator, Constant, norm, XDSLOperator, configuration
 from examples.cfd import init_hat
 from fast.bench_utils import plot_2dfunc
 
@@ -26,6 +26,8 @@
 parser.add_argument("-xdsl", "--xdsl", default=False, type=bool, help="xDSL run")
 args = parser.parse_args()
 
+mpiconf = configuration['mpi']
+
 # Some variable declarations
 nx, ny = args.shape
 nt = args.nt
@@ -42,6 +44,8 @@
 
 grid = Grid(shape=(nx, ny), extent=(2., 2.))
 u = TimeFunction(name='u', grid=grid, space_order=so)
+u2 = TimeFunction(name='u', grid=grid, space_order=so)
+devito_out = TimeFunction(name='u', grid=grid, space_order=so)
 
 a = Constant(name='a')
 # Create an equation with second-order derivatives
@@ -53,6 +57,11 @@
 # Reset our data field and ICs
 init_hat(field=u.data[0], dx=dx, dy=dy, value=1.)
 
+if args.xdsl and args.devito:
+    configuration['mpi'] = 0
+    u2.data[:] = u.data[:]
+    configuration['mpi'] = mpiconf
+
 if args.devito:
     op = Operator([eq_stencil], name='DevitoOperator')
     op.apply(time=nt, dt=dt, a=nu)
@@ -61,10 +70,19 @@
     if args.plot:
         plot_2dfunc(u)
 
+    if args.xdsl:
+        configuration['mpi'] = 0
+        devito_out.data[:] = u.data[:]
+        u.data[:] = u2.data[:]
+        configuration['mpi'] = mpiconf
+
 # Reset data
 init_hat(field=u.data[0], dx=dx, dy=dy, value=1.)
 
 if args.xdsl:
     xdslop = XDSLOperator([eq_stencil], name='XDSLOperator')
     xdslop.apply(time=nt, dt=dt, a=nu)
     print("XDSL Field norm is:", norm(u))
+
+if args.xdsl and args.devito:
+    print("Max error: ", np.max(np.abs(u.data - devito_out.data)))

fast/diffusion_3D_wBCs.py

@@ -6,7 +6,7 @@
 import numpy as np
 
 from devito import (Grid, TimeFunction, Eq, solve, Operator, Constant,
-                    norm, XDSLOperator)
+                    norm, XDSLOperator, configuration)
 from fast.bench_utils import plot_3dfunc
 
 parser = argparse.ArgumentParser(description='Process arguments.')
@@ -26,6 +26,8 @@
 parser.add_argument("-xdsl", "--xdsl", default=False, type=bool, help="xDSL run")
 args = parser.parse_args()
 
+mpiconf = configuration['mpi']
+
 # Some variable declarations
 nx, ny, nz = args.shape
 nt = args.nt
@@ -49,6 +51,7 @@
 
 grid = Grid(shape=(nx, ny, nz), extent=(2., 2., 2.), topology=topology)
 u = TimeFunction(name='u', grid=grid, space_order=so)
+devito_out = TimeFunction(name='u', grid=grid, space_order=so)
 
 a = Constant(name='a')
 # Create an equation with second-order derivatives
@@ -68,6 +71,11 @@
     if args.plot:
         plot_3dfunc(u)
 
+    if args.xdsl:
+        configuration['mpi'] = 0
+        devito_out.data[:] = u.data[:]
+        configuration['mpi'] = mpiconf
+
 if args.xdsl:
     # Reset field
     u.data[:, :, :, :] = 0
@@ -78,3 +86,6 @@
     print("XDSL Field norm is:", norm(u))
     if args.plot:
         plot_3dfunc(u)
+
+if args.xdsl and args.devito:
+    print("Max error: ", np.max(np.abs(u.data - devito_out.data)))

fast/wave2d_b.py

@@ -59,7 +59,7 @@
 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)
+devito_out = TimeFunction(name="u", 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
@@ -73,13 +73,14 @@
 # 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)
 
+if args.xdsl and args.devito:
+    configuration['mpi'] = 0
+    u2.data[:] = u.data[:]
+    configuration['mpi'] = mpiconf
+
 # 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)
 
@@ -94,14 +95,17 @@
     op1 = Operator([stencil], name='DevitoOperator')
     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))
+
+    if args.xdsl:
+        configuration['mpi'] = 0
+        devito_out.data[:] = u.data[:]
+        u.data[:] = u2.data[:]
+        configuration['mpi'] = mpiconf
+
     # 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]))
@@ -125,3 +129,6 @@
     # 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]))
+
+if args.xdsl and args.devito:
+    print("Max error: ", np.max(np.abs(u.data - devito_out.data)))

fast/wave3d_b.py

@@ -62,11 +62,12 @@
 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
 # 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)
+devito_out = TimeFunction(name="u", 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
@@ -80,12 +81,12 @@
 # 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.npz" % (so, shape_str), allow_pickle=True)['arr_0']
 dt = np.load("so%s_critical_dt%s.npy" % (so, shape_str), allow_pickle=True)
+if args.xdsl and args.devito:
+    configuration['mpi'] = 0
+    u2.data[:] = u.data[:]
+    configuration['mpi'] = mpiconf
 
 # 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)
@@ -105,14 +106,16 @@
     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))
+
+    if args.xdsl:
+        configuration['mpi'] = 0
+        devito_out.data[:] = u.data[:]
+        u.data[:] = u2.data[:]
+        configuration['mpi'] = mpiconf
     # 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]))
@@ -132,3 +135,6 @@
     # 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 and args.devito:
+    print("Max error: ", np.max(np.abs(u.data - devito_out.data)))