Add automated benchmark and validation scripts

experiments/benchmarks.py

@@ -5,7 +5,10 @@
     numpy_residual_function, cython_residual_function, numba_residual_function
 from lid_driven_cavity_problem.staggered_grid import Graph
 from lid_driven_cavity_problem.time_stepper import run_simulation
+import matplotlib.pyplot as plt
+import numpy as np
 
+N_RUNS = 5
 
 def run_solver_and_measure_time(solver_type='petsc', language='c++', grid_size=32, Re=400.0):
     if solver_type == 'petsc':
@@ -40,22 +43,49 @@ def run_solver_and_measure_time(solver_type='petsc', language='c++', grid_size=3
     U_bc = (mi * Re) / (rho * size_x)
 
     graph = Graph(size_x, size_y, nx, ny, dt, rho, mi, U_bc)
-    b = time.time()
-    g = run_simulation(graph, final_time, solver, residual_f)
+    runs = []
+    for i in range(N_RUNS):
+        print("Run %s/%s" % (i + 1, N_RUNS,))
+        start = time.time()
+        g = run_simulation(graph, final_time, solver, residual_f)
+        stop = time.time()
+        runs.append(stop - start)
 
     del g
     del solver
 
-    print(time.time() - b)
+    return np.mean(runs)
 
+available_grid_sizes = [8, 16, 32, 64, 128]
+results = {}
 for language in [
-#     'python',
+    'python',
     'numpy',
     'cython',
     'c++',
     'numba',
 ]:
-    for grid_size in [8, 16, 32, 64, 128]:
-        print("%s@%s" % (language, grid_size,))
-        run_solver_and_measure_time(language=language, grid_size=grid_size)
-        print("")
+    results[language] = []
+    for grid_size in available_grid_sizes:
+        if grid_size == 128 and language == 'python':
+            print("Skipping 128x128 mesh with raw Python, as it'll take too long")
+            continue
+
+        print("Running %s@%sx%s" % (language, grid_size, grid_size,))
+        results[language].append(run_solver_and_measure_time(language=language, grid_size=grid_size))
+
+plt.figure(1)
+plt.title("Time comparison")
+ax = plt.gca()
+ax.set_xlabel('Grid Size')
+ax.set_ylabel('Time (s)')
+for language, times in results.items():
+    plt.plot(np.arange(0, len(times), 1), np.array(times), label=language)
+labels = [item.get_text() for item in ax.get_xticklabels()]
+for i, gs in enumerate(available_grid_sizes):
+    labels[i] = '%sx%s' % (gs, gs,)
+plt.xticks(range(len(available_grid_sizes)))
+ax.set_xticklabels(labels)
+ax.set_yscale("log")
+plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
+plt.show()

experiments/ghia_ghia_shin_results/U.txt

@@ -1,15 +1,17 @@
--3.85275436e-2    -3.85425800e-2    -4.1974991e-2       -9.259926e-2        -2.02330048e-1
--6.9584425e-2     -6.96238561e-2    -7.7125399e-2       -1.78748051e-1      -3.478451e-1
--9.6906717e-2     -9.6983962e-2     -1.09816214e-1      -2.6391720e-1       -3.844094e-1
--1.22595555e-1    -1.22721979e-1    -1.41930064e-1      -3.2122908e-1       -3.189461e-1
--1.47461728e-1    -1.47636199e-1    -1.72712391e-1      -3.2025109e-1       -2.456937e-1
--1.71067124e-1    -1.71260757e-1    -1.98470859e-1      -2.6630635e-1       -1.837321e-1
--1.91535923e-1    -1.91677043e-1    -2.12962392e-1      -1.9073056e-1       -1.2341046e-1
--2.05191715e-1    -2.05164738e-1    -2.091491418e-1     -1.15053628e-1      -6.205613e-2
--2.06089397e-1    -2.05770198e-1    -1.82080595e-1      -4.2568947e-2       5.6180e-4
--1.85581148e-1    -1.84928116e-1    -1.31256301e-1      3.024302e-2         6.5248742e-2
--1.322092275e-1   -1.313892353e-1   -6.0245594e-2       1.0545601e-1        1.3357257e-1
--3.2443684e-2     -3.1879308e-2     2.7874448e-2        1.8130685e-1        2.0791461e-1
-1.27054983e-1     1.26912095e-1     1.40425325e-1       2.5220384e-1        2.884424e-1
-3.55228331e-1     3.54430364e-1     3.1055709e-1        3.1682969e-1        3.625454e-1
-6.51176326e-1     6.50529292e-1     5.97466694e-1       4.69580199e-1       4.229321e-1 
+1.0000  1.00000    1.00000
+0.9766  0.84123    0.75837
+0.9688  0.78871    0.68439
+0.9609  0.73722    0.61756
+0.9531  0.68717    0.55892
+0.8516  0.23151    0.29093
+0.7344  0.00332    0.16256
+0.6172  -0.13641   0.02135
+0.5000  -0.20581   -0.11477
+0.4531  -0.21090   -0.17119
+0.2813  -0.15662   -0.32726 
+0.1719  -0.10150   -0.24299 
+0.1016  -0.06434   -0.14612 
+0.0703  -0.04775   -0.10338 
+0.0625  -0.04192   -0.09266 
+0.0547  -0.03717   -0.08186
+0.0000  0.00000    0.00000

experiments/ghia_ghia_shin_results/V.txt

@@ -1,15 +1,17 @@
-9.4572847e-2     9.2970121e-2     9.4807616e-2     1.85132290e-1    2.807057e-1
-1.55984965e-1    1.52547843e-1    1.4924300e-1     2.6225126e-1     3.650418e-1
-1.82641889e-1    1.78781456e-1    1.74342933e-1    2.9747923e-1     3.678527e-1
-1.78849493e-1    1.76415100e-1    1.79243328e-1    3.0096003e-1     3.0710428e-1
-1.51784706e-1    1.52055820e-1    1.69132064e-1    2.6831096e-1     2.3126839e-1
-1.089092434e-1   1.121477612e-1   1.45730201e-1    2.0657139e-1     1.6056422e-1
-5.66144697e-2    6.21048147e-2    1.087758646e-1   1.30571694e-1    9.296931e-2
-6.3677058e-6     6.3603620e-3     5.7536559e-2     5.2058082e-2     2.579946e-2
--5.66033951e-2   -5.10417285e-2   -7.748504e-3     -2.4714514e-2    -4.184068e-2
--1.089027070e-1  -1.056157259e-1  -8.4066715e-2    -1.00884164e-1   -1.107983e-1
--1.51784274e-1   -1.51622101e-1   -1.63010143e-1   -1.82109238e-1   -1.816797e-1
--1.78854716e-1   -1.81633561e-1   -2.27827313e-1   -2.80990219e-1   -2.533815e-1
--1.82650133e-1   -1.87021651e-1   -2.53768577e-1   -4.0004235e-1    -3.315667e-1
--1.55992321e-1   -1.59898186e-1   -2.18690812e-1   -4.4901185e-1    -4.677756e-1
--9.4576294e-2    -9.6409942e-2    -1.23318170e-1   -2.70354943e-1   -4.5615254e-1 
+1.0000    0.00000    0.00000
+0.9688    -0.05906   -0.12146
+0.9609    -0.07391   -0.15663
+0.9531    -0.08864   -0.19254
+0.9453    -0.10313   -0.22847
+0.9063    -0.16914   -0.23827
+0.8594    -0.22445   -0.44993
+0.8047    -0.24533   -0.38598
+0.5000    0.05454    0.05188
+0.2344    0.17527    0.30174
+0.2266    0.17507    0.30203
+0.1563    0.16077    0.28124
+0.0938    0.12317    0.22965
+0.0781    0.10890    0.20920
+0.0703    0.10091    0.19713
+0.0625    0.09233    0.18360
+0.0000    0.00000    0.00000

experiments/run_solver.py

@@ -16,7 +16,7 @@
 
 PLOT_RESULTS = True
 SOLVER_TYPE = 'petsc'
-LANGUAGE = 'c++'
+LANGUAGE = 'numba'
 
 if SOLVER_TYPE == 'petsc':
     solver = petsc_solver_wrapper.solve
@@ -43,13 +43,13 @@
 
 size_x = 1.0
 size_y = 1.0
-nx = 80
-ny = 80
+nx = 180
+ny = 180
 dt = 1e-2
 rho = 1.0
 final_time = None  # Run until steady state
 mi = 1.0
-Re = 10.0
+Re = 400.0
 U_bc = (mi * Re) / (rho * size_x)
 print("Run Parameters:")
 print("size_x = %s" % (size_x,))
@@ -94,43 +94,31 @@
 U = U.reshape(nx, ny - 1)
 V = V.reshape(nx - 1, ny)
 
-pos_ghia = np.loadtxt('ghia_ghia_shin_results/pos.txt')
 U_ghia = np.loadtxt('ghia_ghia_shin_results/U.txt')
 V_ghia = np.loadtxt('ghia_ghia_shin_results/V.txt')
-if np.isclose(Re, 0.01):
-    U_ghia = U_ghia[:, 0]
-    V_ghia = V_ghia[:, 0]
-elif np.isclose(Re, 10.0):
+pos_U_ghia = U_ghia[:, 0]
+pos_V_ghia = V_ghia[:, 0]
+if np.isclose(Re, 100.0):
     U_ghia = U_ghia[:, 1]
     V_ghia = V_ghia[:, 1]
-elif np.isclose(Re, 100.0):
+elif np.isclose(Re, 400.0):
     U_ghia = U_ghia[:, 2]
     V_ghia = V_ghia[:, 2]
-elif np.isclose(Re, 400.0):
-    U_ghia = U_ghia[:, 3]
-    V_ghia = V_ghia[:, 3]
-elif np.isclose(Re, 1000.0):
-    U_ghia = U_ghia[:, 4]
-    V_ghia = V_ghia[:, 4]
-else:
-    print("Re=%s" % (Re,))
-    U_ghia = None
-    V_ghia = None
 
 plt.figure(3)
 plt.title("U velocity in the mesh center-x")
 U_normalized = U / U_bc
 U_center = U_normalized[:, len(U) // 2]
 plt.plot(U_center, np.linspace(0.0, size_y, len(U_center)))
 if U_ghia is not None:
-    plt.plot(U_ghia, pos_ghia, 'xb')
+    plt.plot(U_ghia, pos_U_ghia, 'xb')
 
 plt.figure(4)
 plt.title("V velocity in the mesh center-y")
 V_normalized = V / U_bc
 V_center = V_normalized[len(V) // 2, :]
 plt.plot(np.linspace(0.0, size_x, len(V_center)), V_center)
 if V_ghia is not None:
-    plt.plot(pos_ghia, V_ghia, 'xb')
+    plt.plot(pos_V_ghia, V_ghia, 'xb')
 
 plt.show()

experiments/validation.py

@@ -0,0 +1,102 @@
+from lid_driven_cavity_problem.nonlinear_solver import petsc_solver_wrapper
+from lid_driven_cavity_problem.residual_function import pure_python_residual_function, \
+    numpy_residual_function, cython_residual_function, numba_residual_function
+from lid_driven_cavity_problem.staggered_grid import Graph
+from lid_driven_cavity_problem.time_stepper import run_simulation
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+grid_size = 129
+Re = 400.0
+size_x = 1.0
+size_y = 1.0
+nx = grid_size
+ny = grid_size
+dt = 1e-2
+rho = 1.0
+final_time = None  # Run until steady state
+mi = 1.0
+U_bc = (mi * Re) / (rho * size_x)
+
+
+def run_solver_and_return_results(language, interpolation_type):
+    solver = petsc_solver_wrapper._PetscSolverWrapper()
+
+    if language == 'python':
+        residual_f = pure_python_residual_function.residual_function
+    elif language == 'numpy':
+        residual_f = numpy_residual_function.residual_function
+    elif language == 'cython':
+        residual_f = cython_residual_function.residual_function
+    elif language == 'c++':
+        from lid_driven_cavity_problem.residual_function import cpp_residual_function
+        residual_f = cpp_residual_function.residual_function
+    elif language == 'numba':
+        residual_f = numba_residual_function.residual_function
+    else:
+        assert False
+
+    graph = Graph(size_x, size_y, nx, ny, dt, rho, mi, U_bc)
+    if interpolation_type == 'cds':
+        graph.use_cds = True
+    g = run_simulation(graph, final_time, solver, residual_f)
+
+    U, V, P = \
+        np.array(g.ns_x_mesh.phi).reshape(nx, ny - 1), \
+        np.array(g.ns_y_mesh.phi).reshape(nx - 1, ny), \
+        np.array(g.pressure_mesh.phi)
+
+    del g
+    del solver
+
+    return U, V, P
+
+
+results = {}
+for language in [
+#     'python',
+#     'numpy',
+#     'cython',
+#     'c++',
+    'numba',
+]:
+    for interpolation_type in ['cds', 'uds']:
+        print("Running %s %s..." % (language, interpolation_type,))
+        results['%s interpolation' % (interpolation_type,)] = run_solver_and_return_results(language, interpolation_type)
+
+print("All done. Preparing plot...")
+
+U_ghia = np.loadtxt('ghia_ghia_shin_results/U.txt')
+V_ghia = np.loadtxt('ghia_ghia_shin_results/V.txt')
+pos_U_ghia = U_ghia[:, 0]
+U_ghia = U_ghia[:, 2]
+pos_V_ghia = V_ghia[:, 0]
+V_ghia = V_ghia[:, 2]
+
+plt.figure(1)
+plt.title("U velocity in the mesh center-x")
+ax = plt.gca()
+ax.set_xlabel('U (m/s)')
+ax.set_ylabel('y (m)')
+for language, (U, V, P) in results.items():
+    U_normalized = U / U_bc
+    U_center = U_normalized[:, len(U) // 2]
+    plt.plot(U_center, np.linspace(0.0, size_y, len(U_center)), label=language)
+plt.plot(U_ghia, pos_U_ghia, 'xb', label='Ghia, Ghia and Shin')
+plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
+
+plt.figure(2)
+plt.title("V velocity in the mesh center-y")
+ax = plt.gca()
+ax.set_xlabel('x (m)')
+ax.set_ylabel('V (m/s)')
+for language, (U, V, P) in results.items():
+    V_normalized = V / U_bc
+    V_center = V_normalized[len(V) // 2, :]
+    plt.plot(np.linspace(0.0, size_x, len(V_center)), V_center, label=language)
+plt.plot(pos_V_ghia, V_ghia, 'xb', label='Ghia, Ghia and Shin')
+plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
+
+print("All done.")
+plt.show()