I am running an elasticity simulation in the Jupyter Notebook kernel provided in the FEniCSx tutorial. The solution is found using a custom Newton solver that has been validated with other scripts.
Early in the day, the simulation was running fine, but recently, when I start to run the solver loop, my Jupyter kernel dies immediately and restarts. I am still able to run those other validation scripts, so there doesn’t seem to be a bug in the solver. The mesh I am working with is relatively small (and smaller that previous meshes I have worked with), so I do not think it is a memory overflow issue.
I went ahead and directly installed dolfinx through anacoda so that I would be using a local Jupyter kernel, but the same problem is persisting.
Any suggestions for next steps to take would be greatly appreciated! The custom test script that is killing the kernel is included below
import gmsh
import numpy as np
from mpi4py import MPI
import matplotlib.pyplot as plt
import pyvista
from dolfinx.io.gmshio import model_to_mesh
from dolfinx.plot import create_vtk_mesh
import numpy as np
import ufl
import dolfinx
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import petsc4py
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot, io, cpp, nls, la
import matplotlib.pyplot as plt
x = np.linspace(-1,1,50)
y = 1.1 + 0.3*x**2
x = np.append(x,[1,-1,-1])
y = np.append(y,[2,2,1.4])
rank = MPI.COMM_WORLD.rank
gmsh.initialize()
############## Mesh Generation ##############
gdim = 2 # Geometric dimension of the mesh
model_rank = 0
mesh_comm = MPI.COMM_WORLD
if mesh_comm.rank == model_rank:
for i in range(len(x)-1):
gmsh.model.occ.addPoint(x[i],y[i],0,i+1)
for i in range(len(x)-2):
gmsh.model.occ.addLine(i+1,i+2,i+1)
gmsh.model.occ.addLine(len(x)-1,1,len(x)-1)
indenter_loop = gmsh.model.occ.addCurveLoop(np.arange(1,len(x)),1)
gmsh.model.occ.synchronize()
indenter = gmsh.model.occ.addPlaneSurface([indenter_loop],1)
floor = gmsh.model.occ.addRectangle(-2,0,0,4,1,2)
background = gmsh.model.occ.addRectangle(-1,1,0,2,1,3)
gmsh.model.occ.synchronize()
bodies = [(2,indenter),(2,floor)]
gmsh.model.occ.synchronize()
whole_domain = gmsh.model.occ.fragment([(2,background)],bodies)
gmsh.model.occ.synchronize()
background_surfaces = []
other_surfaces = []
for i,domain in enumerate(whole_domain[0]):
com = gmsh.model.occ.getCenterOfMass(domain[0],domain[1])
gmsh.model.addPhysicalGroup(domain[0], [domain[1]], tag=i)
if (np.abs(com[1]-1.12)<0.001):
background_surfaces.append(domain)
else:
other_surfaces.append(domain)
gmsh.model.mesh.field.add("Distance", 1)
r = 0.1
edges = gmsh.model.getBoundary(other_surfaces, oriented=False)
gmsh.model.mesh.field.setNumbers(1, "EdgesList", [e[1] for e in edges])
gmsh.model.mesh.field.add("Threshold", 2)
gmsh.model.mesh.field.setNumber(2, "IField", 1)
gmsh.model.mesh.field.setNumber(2, "LcMin", r * 0.75 )
gmsh.model.mesh.field.setNumber(2, "LcMax", 2 * r)
gmsh.model.mesh.field.setNumber(2, "DistMin", 2 * r)
gmsh.model.mesh.field.setNumber(2, "DistMax", 4 * r)
gmsh.model.mesh.field.setAsBackgroundMesh(2)
gmsh.option.setNumber("Mesh.Algorithm", 7)
gmsh.model.mesh.generate(2)
domain, ct, _ = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
gmsh.finalize()
pyvista.set_jupyter_backend("ipygany")
plotter = pyvista.Plotter()
topology, cell_types, geometry = create_vtk_mesh(domain, domain.topology.dim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry )
num_local_cells = domain.topology.index_map(domain.topology.dim).size_local
grid.cell_data["Marker"] = ct.values[ct.indices<num_local_cells]
grid.set_active_scalars("Marker")
actor = plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
pyvista.start_xvfb()
cell_tag_fig = plotter.screenshot("cell_tags.png")
### Defining the material properties ###
QDoF = 2
def safeln(x):
return ufl.ln(x + 1e-8)
def safesqrt(x):
return ufl.sqrt(x + 1e-8)
B = fem.Constant(domain,PETSc.ScalarType((0,-1)))
W_rigid = lambda F: 0.5*((210e9)/(2*(1+0.25)))*(ufl.tr(F.T * F) -2 - 2*ufl.ln(ufl.det(F)))+ (((210e9)*0.25/((1+0.25)*(1-2*0.25)))/2)*((ufl.det(F)-1)**2)
def W_air(F):
C = F.T * F
I1 = ufl.tr(C)
a1 = 1e6
W = a1*(I1 - 2 - 2*ufl.ln(ufl.det(F)))
return W
def top(x):
return np.isclose(x[1],2)
def bottom(x):
return np.isclose(x[1],0)
fdim = domain.topology.dim -1
top_facets = mesh.locate_entities_boundary(domain,fdim,top)
bottom_facets = mesh.locate_entities_boundary(domain,fdim,bottom)
marked_facets = np.hstack([top_facets,bottom_facets])
marked_values = np.hstack([np.full(len(top_facets),4,dtype=np.int32),\
np.full(len(bottom_facets),5,dtype=np.int32)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain,fdim,marked_facets[sorted_facets],marked_values[sorted_facets])
### Defining the function spaces and functions for elastodynamics ###
V_vec = fem.VectorFunctionSpace(domain,("CG",QDoF)) # function space for vector functions
metadata = {"quadrature_degree":QDoF**2}
ds = ufl.Measure('ds',domain=domain,metadata=metadata,subdomain_data=ct)
dx = ufl.Measure("dx",domain=domain,metadata=metadata,subdomain_data=ct)
vals = ct.values[ct.indices<num_local_cells]
rigid_cells = np.where((vals==0)|(vals==1))[0]
air_cells = np.where(vals>=2)[0]
u_bc_bottom = fem.Constant(domain, PETSc.ScalarType((0,0)))
u_bc_top = fem.Constant(domain, PETSc.ScalarType((0)))
top_dofs = fem.locate_dofs_topological(V_vec.sub(1),facet_tag.dim,facet_tag.indices[facet_tag.values==4])
bottom_dofs = fem.locate_dofs_topological(V_vec.sub(1),facet_tag.dim,facet_tag.indices[facet_tag.values==5])
bcs = [fem.dirichletbc(u_bc_top,top_dofs,V_vec.sub(1)),\
fem.dirichletbc(u_bc_bottom,bottom_dofs,V_vec)]
u_ = ufl.TestFunction(V_vec) # test function for the displacement
u = fem.Function(V_vec) # displacement at current time step
def sigma(r,W):
I = ufl.variable(ufl.Identity(2))
F = ufl.variable(I + ufl.grad(r))
psi = W(F)
return ufl.diff(psi,F)
def k(u,u_):
kk_ = ufl.inner(sigma(u,W_air),ufl.grad(u_))*dx(0) + \
ufl.inner(sigma(u,W_air),ufl.grad(u_))*dx(1) + \
ufl.inner(sigma(u,W_air),ufl.grad(u_))*dx(2)
return kk_
res = k(u,u_)
### Custom Solver ###
F_form = fem.form(res) # form of the residual function
J = ufl.derivative(res,u) # jacobian of the residual function
J_form = fem.form(J) # form of the jacobian
du = fem.Function(V_vec) # increment of the solution
u_start = fem.Function(V_vec) # holds previous solution for reverting in case of convergence failure
A = fem.petsc.create_matrix(J_form) # preallocating sparse jacobian matrix
L = fem.petsc.create_vector(F_form) # preallocating residual vector
solver = PETSc.KSP().create(domain.comm)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
opts = PETSc.Options()
prefix = f"solver_{id(solver)}"
solver.setOptionsPrefix(prefix)
option_prefix = solver.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
solver.setFromOptions()
solver.setOperators(A)
# iteration parameters
max_it = int(1000)
tol = 1e-4
min_alpha = 1e-2
def NewtonSolve(print_steps=True,print_solution=True,print_alphas=False):
i = 0 # number of iterations of the Newton solver
converged = False
residuals = []
alphas = []
### Initializing ###
with L.localForm() as loc_L:
loc_L.set(0)
A.zeroEntries()
fem.petsc.assemble_matrix(A, J_form, bcs=bcs)
A.assemble()
fem.petsc.assemble_vector(L, F_form)
L.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
L.scale(-1)
# Compute b - J(u_D-u_(i-1))
fem.petsc.apply_lifting(L, [J_form], [bcs], x0=[u.vector], scale=1)
# Set dx|_bc = u_{i-1}-u_D
fem.petsc.set_bc(L, bcs, u.vector, 1.0)
L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)
err_prev = 1e2
while i<max_it:
# Solve linear problem
solver.solve(L, du.vector)
du.x.scatter_forward()
u_start.x.array[:] = u.x.array[:]
err = np.Inf
alpha = 2
# looking for decreasing step size to set relaxation parameter
while ((err>err_prev) and (alpha > min_alpha)):
# search for the largest alpha that decreases the magnitude of the residual
alpha = max(min_alpha,0.5*alpha)
if print_alphas:
print(f" Droppping alpha to {alpha}")
u.x.array[:] = u_start.x.array[:]
u.x.array[:] += alpha * du.x.array
with L.localForm() as loc_L:
loc_L.set(0)
fem.petsc.assemble_vector(L, F_form)
L.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
L.scale(-1)
# Compute b - J(u_D-u_(i-1))
fem.petsc.apply_lifting(L, [J_form], [bcs], x0=[u.vector], scale=1)
# Set dx|_bc = u_{i-1}-u_D
fem.petsc.set_bc(L, bcs, u.vector, 1.0)
L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)
err = np.abs(np.dot(L,du.x.array))
# Compute norm of update
correction_norm = du.vector.norm(0)
error_norm = L.norm(0)
residuals.append(error_norm)
alphas.append(alpha)
# building A for next iteration
A.zeroEntries()
fem.petsc.assemble_matrix(A, J_form, bcs=bcs)
A.assemble()
i += 1
if print_steps:
print(f" Iteration {i}: Correction norm {correction_norm}, Residual {error_norm}, Relaxation Parameters {alpha}")
if correction_norm < tol:
converged = True
break
if print_solution:
if converged:
print(f" Solution reached in {i} iterations.")# (Residual norm {error_norm})")
if not converged:
print(f"No solution found after {i} iterations. Revert to previous solution and adjust solver parameters.")
plt.plot(np.arange(i),residuals)
plt.xlabel('Iterations')
plt.ylabel('Residual')
plt.show()
plt.plot(np.arange(i),alphas)
plt.xlabel('Iterations')
plt.ylabel('Alpha')
plt.show()
return converged, u_start, i, residuals, alphas
T = 1
t = 0
dt = 1e-4
count = 0
N=0
# iteration
while t<T:
top_ = - 0.1*t #1*np.sin(2*np.pi*t/3)**2
u_bc_top.value = top_
converged,u_prev,n,residuals, alphas = NewtonSolve(print_steps=False,print_solution=False)#solver.solve(u)
if converged:
N += n
u.x.scatter_forward()
count +=1
if (count%np.floor(1*(1/dt)/100))==0:
print('t = ',t,' (',t/T*100,'%) | top = ',top_)
print("Average number of interations: ",N/np.floor(1*(1/dt)/100))
N=0
if (count%np.floor(10*(1/dt)/100))==0:
plot_function(t,u,grid=False,ref_grid=True)
else:
print('Failed to converge')
break
t += dt