evalf
Slow (hierarchical) basis calculation on highly refined meshes
When a mesh or refined in a small region around some feature, the calculation of the hierarchical basis slows down considerably with the increase of the refinement. After about ~10 refinement iterations the increase in time is exponential, while the increase in ndof is typically more linear or quadratic. This issue seems to stem from the way a hierarchical basis is constructed that involved the calculation of all basis functions on each level of the mesh, most of which are not required. For higher refinement levels it might therefore be beneficial to have an option to only calculate the basis functions that are actually needed, or to calculate them on-the-fly. Such an approach would most likely slightly degrade the performance at low refinement levels, thus having it as an option would be great.
To illustrate the problem I am facing I have made a simple script to benchmark the performance at different levels. This code starts with a 2x2 mesh, which is subsequently refined continuously in the corner. An overview of the ndof and time to calculate the basis is plotted below, clearly showing the exponential increase in time.
from nutils import mesh, function, solver, util, export, cli, testing
from nutils.expression_v2 import Namespace
import numpy
import treelog
import time
import matplotlib.pyplot as plt
def main():
domain, geom = mesh.unitsquare(2, 'square')
x, y = geom
plt.figure()
with treelog.iter.fraction('level', range(24)) as lrange:
for irefine in lrange:
if irefine:
refdom = domain.refined
ns.refbasis = refdom.basis('th-std', degree=3)
indicator = refdom.integral('refbasis_n (x_0 x_0 + x_1 x_1)' @ ns, degree=3*2).eval()
supp = ns.refbasis.get_support(indicator == numpy.min(indicator))
domain = domain.refined_by(refdom.transforms[supp])
ns = Namespace()
ns.x = geom
t = time.time()
nrun=1
for i in range(nrun):
ns.basis = domain.basis('th-std', degree=3)
btime = (time.time() - t)/nrun
treelog.user('Level', irefine, 'ndof', len(ns.basis) , 'time', btime)
plt.semilogy(irefine,btime,'r.')
plt.semilogy(irefine,len(ns.basis),'b.')
plt.xlabel('Refinement step')
plt.ylabel('(blue) ndof; (red) basis time [s]')
plt.tight_layout(pad=.1)
plt.savefig("Refinement_benchmark.png", dpi=300)
plt.savefig("Refinement_benchmark.pdf", dpi=300)
plt.show()
if __name__ == '__main__':
cli.run(main)
Though a basis is indeed formed at every level, we tried to make this a lightweight operation via #380 precisely to avoid exponential slowdown. While this did improve things a lot, it is certainly possible that we still do more work than strictly necessary, as your results suggest, and we will use your benchmarking code to try to identify remaining bottlenecks. Keeping this issue open until we have a clearer picture on this.
