add performance test with scikitfem

Author: lukeolson

performance/demo.py

@@ -0,0 +1,67 @@
+
+"""A simple performance test adopted from sciket-fem.
+"""
+from timeit import timeit
+import numpy as np
+import skfem as skf
+from skfem.models.poisson import laplace, unit_load
+import pyamg
+
+def pre(N=3):
+    m = skf.MeshTet.init_tensor(*(3 * (np.linspace(0., 1., N),)))
+    return m
+
+def assembler(m):
+    basis = skf.Basis(m, skf.ElementTetP1())
+    return (
+        laplace.assemble(basis),
+        unit_load.assemble(basis),
+    )
+
+times  = []
+fw = 10
+kmin =  6
+kmax = 20
+Nlist = [int(2 ** (k / 3)) for k in range(kmin, kmax)]
+# print(Nlist)
+
+print('| N | DoFs | Assembly | Solve prep | Solve setup | Solve |')
+print('| ---- | ---- | ---- | ---- | ---- | ---- |')
+for N in Nlist:
+    m = pre(N)
+
+    assemble_time = timeit(lambda: assembler(m), number=1)
+    A, b = assembler(m)
+    D = m.boundary_nodes()
+
+    condense_time = timeit(lambda: skf.condense(A, b, D=D), number=1)
+    A, b, _, _ = skf.condense(A, b, D=D)
+
+    setup_time = timeit(lambda: pyamg.smoothed_aggregation_solver(A).aspreconditioner(), number=1)
+    ml = pyamg.smoothed_aggregation_solver(A,
+                                           max_coarse=10).aspreconditioner()
+
+    mlsolver = skf.solver_iter_pcg(verbose=False, M=ml, rtol=1e-8)
+    solve_time = timeit(lambda: skf.solve(A, b, solver=mlsolver), number=1)
+
+    times.append([N, len(b), assemble_time, condense_time, setup_time, solve_time])
+    print(f'| {N:10d} | {len(b):10d} | {assemble_time:>{fw}.5f} | {condense_time:{fw}.5f} | {setup_time:{fw}.5f} | {solve_time:{fw}.5f} |')
+
+import matplotlib.pyplot as plt
+fig, ax = plt.subplots()
+n = [t[1] for t in times]
+ax.loglog(n, [t[2] for t in times], label='Assembly')
+ax.loglog(n, [t[3] for t in times], label='Solve prep')
+ax.loglog(n, [t[4] for t in times], label='Solve setup')
+ax.loglog(n, [t[5] for t in times], label='Solve')
+ax.set_xlabel('# DoFs')
+ax.set_ylabel('time (s)')
+ax.grid(True)
+plt.legend()
+
+figname = f'./output/performance.png'
+import sys
+if '--savefig' in sys.argv:
+    plt.savefig(figname, bbox_inches='tight', dpi=150)
+else:
+    plt.show()

performance/output/performance.png

@@ -0,0 +1,9 @@
+This demo shows the performance of finite element assembly (by way of scikit-fem)
+and multigrid setup/solve (within PCG).  The solver is run to a tolerance of `1e-8`.
+The figure shows:
+
+- `DoFs`: the total number of degrees of freedom in the system
+- `Assembly`: the total time to assemble the FE matrix (scikit-fem)
+- `Solve prep`: the total time to condense the system to non-Dirichlet nodes (scikit-fem)
+- `Solve setup`: the total time for the AMG setup phase (pyamg)
+- `Solve`: the total time for the AMG solve phase (pyamg) withing PCG (scikit-fem)

performance/readme.md

@@ -3,3 +3,4 @@ matplotlib
 vedo
 pyinstrument
 pyyaml
+scikit-fem

readme.md

@@ -1,9 +1,10 @@
 import subprocess
 import yaml
 import os
+import sys
 from glob import glob
 
-def exectute_demo(exampledir, name='demo.py'):
+def execute_demo(exampledir, name='demo.py'):
     """Exectue a demo in a particular directory."""
     nameall = name.split()
     output = subprocess.run(['python3'] + nameall + ['--savefig'],
@@ -75,6 +76,8 @@ def exectute_demo(exampledir, name='demo.py'):
 Other:
   - dir: profile_pyamg
     title: Profiling Performance
+  - dir: performance
+    title: Scaling performance of AMG and FE assembly
 """)
 
 # generate table and header
@@ -89,6 +92,20 @@ def exectute_demo(exampledir, name='demo.py'):
 
 main = '\n'
 
+dirs = None
+if len(sys.argv) > 1:
+    dirs = sys.argv[1:]
+    for d in dirs:
+        found = False
+        for section in toc:
+            for demo in toc[section]:
+                if d == demo['dir']:
+                    found = True
+        if not found or not os.path.isdir(d):
+            print('usage: runner.py')
+            print('usage: runner.py dir1 [dir2] ...')
+            exit()
+
 for section in toc:
 
     # add to the TOC
@@ -100,7 +117,7 @@ def exectute_demo(exampledir, name='demo.py'):
 
     if toc[section] is not None:
         for demo in toc[section]:
-            print(f'Processing {demo["dir"]}')
+            print(f'Processing {demo["dir"]}.', end=' ')
             title = demo.get('title', None)
             if title:
                 hrefid = title.replace(' ', '').lower()
@@ -117,10 +134,18 @@ def exectute_demo(exampledir, name='demo.py'):
             main += readmeoutput
 
             # get the demo output
+            runit = True
+            if dirs is not None:
+                if demo['dir'] not in dirs:
+                    runit = False
+                    print('')
+            if runit:
+                print('[--->rerunning]')
             for demoname in demonames:
-                output = exectute_demo(demo['dir'], name=demoname)
-                if len(output) > 0:
-                    main += '\n```\n' + output + '```\n'
+                if runit:
+                    output = execute_demo(demo['dir'], name=demoname)
+                    if len(output) > 0:
+                        main += '\n```\n' + output + '```\n'
 
             # get the output figs
             figs = glob(os.path.join(f'{demo["dir"]}', 'output') +'/*.png')