Thanks for the suggestion, it helped me make a simple MWE which surprisingly reproduces the bug quite easily

```
import dolfinx, mpi4py, ufl, gmsh, petsc4py, pyvista
import pyvista
from dolfinx.io import gmshio
import numpy as np
import matplotlib.pyplot as plt
import utils
from dolfinx.fem.petsc import LinearProblem
from tqdm.auto import tqdm
def create_geom(ro=50, fov=None, disp=False, write=False, rot_angle=0):
'''a simple geometry where the inner circle rotates by the specified angle'''
if fov is None:
fov = (400, 300)
fov_x, fov_y = fov
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0) # disable output messages
gmsh.clear()
gdim = 2
model_rank = 0
t_gap = 5 # thickness of airgap
occ = gmsh.model.occ
# first we create the outer domain
outer_box = occ.add_rectangle(-fov_x/2, -fov_y/2, 0, fov_x, fov_y)
air_shell = occ.addDisk(0, 0, 0, ro+t_gap, ro+t_gap)
occ.synchronize()
occ.rotate([(gdim, air_shell)], 0, 0, 0, 0, 0, 1, rot_angle)
occ.synchronize()
occ.fragment([(gdim, outer_box)], [(gdim, air_shell)])
occ.synchronize()
# number all domains
all_doms = gmsh.model.getEntities(gdim)
for j, dom in enumerate(all_doms):
gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1) # create the main group/node
# number all boundaries
all_edges = gmsh.model.getEntities(gdim - 1)
for j, edge in enumerate(all_edges):
gmsh.model.addPhysicalGroup(edge[0], [edge[1]], j + 1) # create the main group/node
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.refine()
gmsh.model.mesh.refine()
if write:
gmsh.write('geom/mwe_multiphenicsx.msh')
model_rank = 0
mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)
if disp:
gmsh.fltk.run() # visible inspection, for debugging
gmsh.finalize()
return mesh, ct, ft
from scipy.interpolate import griddata
def interp_dofs(uh_data, uh_target):
'''return dofs on target function based on the data function'''
points_data = uh_data.function_space.tabulate_dof_coordinates()[:,:2]
points_target = uh_target.function_space.tabulate_dof_coordinates()[:,:2]
return griddata(points_data, uh_data.x.array, points_target)
t_vec = np.linspace(0,1000, 10000)*1e-3
w0 = 2*np.pi*30 # angular frequency
dt = t_vec[1]-t_vec[0]
# boundary tags
bnd_idx_l = (2, ) # left boundary
bnd_idx_r = (3, ) # right boundary
uh_t = t_vec.size*[[]]
for idx_t in tqdm(range(t_vec.size), total=t_vec.size):
rot_angle = w0*t_vec[idx_t]
mesh, ct, ft = create_geom(rot_angle=rot_angle)
Omega = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))
u = ufl.TrialFunction(Omega)
v = ufl.TestFunction(Omega)
uh0 = dolfinx.fem.Function(Omega) # recreate with new mesh to store the previous solution
if idx_t == 0:
uh0.x.array[:] = 0.0 # first time step, initialize to zero
else:
uh0.x.array[:] = interp_dofs(uh_t[idx_t-1], uh0) # interpolate the old solution to this mesh
uh0.x.scatter_forward() # is it necessary?
# some physics
F = ufl.dot(ufl.grad(u), ufl.grad(v))*ufl.dx - ufl.dot(uh0, v)*ufl.dx
a = ufl.lhs(F)
L = ufl.rhs(F)
# boundary conditions
dofs_l = dolfinx.fem.locate_dofs_topological(Omega, ft.dim, ft.find(bnd_idx_l[0]))
dbc_l = dolfinx.fem.dirichletbc(1.0, dofs_l, Omega)
dofs_r = dolfinx.fem.locate_dofs_topological(Omega, ft.dim, ft.find(bnd_idx_r[0]))
dbc_r = dolfinx.fem.dirichletbc(0., dofs_r, Omega)
bcs = [dbc_l, dbc_r]
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu",
"pc_factor_mat_solver_type": "mumps"})
uh = problem.solve()
uh.x.scatter_forward() # why do we need this?
uh_t[idx_t] = uh.copy() # for next step and postprocessing
```

I get the crash at 224th iteration. Thanks in advance for the help!