Using a precompiled form over different meshes

General
1

In a substructuring problem, I would like to evaluate a set of forms (~ 10) over a very large number of meshes (1.000), with the long-term goal of scaling the process to even larger systems (e.g. 10k meshes).

I’ve found an example from 2015 that suggests using a ‘mask field’ to select particular subdomains, however, this has the side effect of turning the problem from linear to quadratic in complexity, as now one must compute a larger number of larger solutions?

Has any development taken place that would allow to use a precompiled form over multiple domains?

2

Could you clarify what you want to do? In the opening you mention multiple meshes, and at the end multiple domains. Do you want to assemble the same form (kernel) over many meshes, or over many (sub)domains of a common mesh?

3

Thanks for your reply, and apologies for the confusing wording. I have entirely separated meshes for every component. The meshes need not have the same number of elements or topology (else I could just deform one into the other and then run the precompiled form), but they are made of the same types of elements. For every mesh we have then an individual function space. In these function spaces, we would like to compute around 15 integrals (the same integral over every different mesh).

I have taken a look at the UFL C code and it seems to generate a for loop that iterates over building blocks (triangles, tetrahedra and so on) with some switch logic to decide which integral to compute – In this sense, it should be able to evaluate integrals over any function space that is made of the same building blocks, but I have not been able to find a (documented) way to take advantage of this.

For reference, my goal is to model metamaterial systems that can contain up to a million sites (that’s 100 x 100 x 100 unit cells), so I. would like to avoid generating the full geometry at once.

4

Here’s an example. It moves the mesh-independent definitions and code generation outside of the loop over meshes. The two loops below do the same thing, but you’ll that the first is much faster (about 50x faster on my computer).

It’s a bit messy from Python because we don’t expose a number C++ convenience functions in Python since the DOLFINx Python code takes care of some of the data marshalling before going into the underlying C++ layer. It could be made simpler.

import time

import basix
import dolfinx.fem.petsc
import ufl
from dolfinx import cpp as _cpp
from dolfinx import fem, jit, mesh
from mpi4py import MPI
from ufl import dx, grad, inner

# Create meshes
msh0 = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                             points=((0.0, 0.0), (2.0, 1.0)), n=(32, 16),
                             cell_type=mesh.CellType.triangle)
msh1 = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                             points=((0.0, 0.0), (2.0, 1.0)), n=(32, 16),
                             cell_type=mesh.CellType.triangle)

# Create abstract UFL elements, meshes, spaces and forms
e = basix.ufl.element("Lagrange", "triangle", 1)
coord_element = basix.ufl.element("Lagrange", "triangle", 1, rank=1)
mesh_ufl = ufl.Mesh(coord_element)
V_ufl = ufl.FunctionSpace(mesh_ufl, e)
u, v = ufl.TrialFunction(V_ufl), ufl.TestFunction(V_ufl)
a = inner(grad(u), grad(v)) * dx

# Just-in time compile the UFL form
ufcx_form, module_form, code = jit.ffcx_jit(MPI.COMM_WORLD, a)

form_ptr = module_form.ffi.cast("uintptr_t", module_form.ffi.addressof(ufcx_form))

# Just-in time compile the UFL element and dofmap data
(ufcx_element, ufcx_dofmap), module_element, code = jit.ffcx_jit(MPI.COMM_WORLD, e)

# Create C++ element from JIT-ed element
element_cpp = _cpp.fem.FiniteElement_float64(
    module_element.ffi.cast("uintptr_t", module_element.ffi.addressof(ufcx_element)))

# Pointer (integer) to compiled dofmap data
dofmap_ptr = module_element.ffi.cast("uintptr_t", module_element.ffi.addressof(ufcx_dofmap))

t0 = time.time()
meshes = 10*[msh0._cpp_object, msh1._cpp_object]
for msh in meshes:
    # Create C++ dofmap for current mesh
    cpp_dofmap = _cpp.fem.create_dofmap(
        MPI.COMM_WORLD, dofmap_ptr, msh.topology, element_cpp)

    # Create C++ function space for current mesh (holder for a mesh, element and dofmap)
    Vcpp = _cpp.fem.FunctionSpace_float64(msh, element_cpp, cpp_dofmap)

    # Create C++ bilinear form
    a_cpp = _cpp.fem.Form_float64(form_ptr, spaces=[Vcpp, Vcpp],
                                  coefficients=[], constants=[], subdomains={}, mesh=None)

    # Create PETSc sparse matrix
    A = _cpp.fem.petsc.create_matrix(a_cpp)

    # Pack coefficient (empty in this case) and assemble matrix
    coeffs = _cpp.fem.pack_coefficients(a_cpp)
    _cpp.fem.petsc.assemble_matrix(A=A, a=a_cpp, constants=[], coeffs=coeffs, bcs=[])
    A.assemble()
    # print("Norm (new):", A.norm())

t1 = time.time()
print("Time (A):", t1-t0)

# Reference case
t0 = time.time()
P1 = basix.ufl.element("Lagrange", msh0.basix_cell(), 1)
meshes = 10*[msh0, msh1]
for msh in meshes:
    V = fem.FunctionSpace(msh, P1)
    u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
    a = fem.form(inner(grad(u), grad(v)) * dx)

    # Create PETSc sparse matrix
    A = dolfinx.fem.petsc.assemble_matrix(a)
    A.assemble()
    # print("Norm (ref):", A.norm())

t1 = time.time()
print("Time (B):", t1-t0)
4 Likes
5

Thanks for the answer.

This is very, very, very helpful.

6

Thanks again for the very helpful answer. I’m now back to this (the wheels of research grind slowly…) and have a couple follow-up questions:

(1) How would one do the same for scalar and vector assembly --instead of matrices?
(2) Can the C++ function spaces be reused among multiple evaluations, e.g. different forms but same function space? Or are they very light to instantiate?