I am running into a bizarre, unpredictable behaviour I cannot understand.
A py.script runs fine, on the command line or VS Code, but it freezes (9 out of 10 times, totally erratically) when ran on a Jupyter notebook.
It seems the last extract_deflection function is causing all the fuss when ran on a Jupyter notebook. I was wondering if anything obvious sprung to mind.
I have also doubts about how I set the solve_inplace function. I append the intermediate solutions to a list u_hist which the function returns, while the last solution is saved inplace. I wonder if I am doing something funny with the memory.
Second question: My aim is to run plenty of these solve_inplace calls over a wholegrid of input parameters. Do I need to delete the Nonlinear Problem before creating a new one? I noticed for example here, projection.py that the professionals “delete” PETSc object as they become useless. I see a method __del__ exist for NewtonSolver, but I get the error message 'NewtonSolver' object has no attribute '__del__'. Thank you so much
A MWE is attached.The main file is
import numpy as np
import dolfinx
from dolfinx import fem, mesh, plot
from utils_copy import create_mesh
from utils_copy import solve_inplace
from utils_copy import extract_deflection
import matplotlib.pyplot as plt
domain = create_mesh(20,1,5,3)
V = fem.VectorFunctionSpace(domain, ("CG", 2))
u = fem.Function(V)
results = solve_inplace(domain, V, u, 1.0e4,0.3)
V = fem.VectorFunctionSpace(domain, ("CG", 2))
u = fem.Function(V)
results_lin = solve_inplace(domain, V, u, 1.0e4,0.3, LINEAR=True)
linear_deflection = [extract_deflection(domain, results_lin[i])[1] for i in range(9)]
nonlinear_deflection = [extract_deflection(domain, results[i])[1] for i in range(9)]
print("Made it")
and some utilities functions, basically out the Fenicsx tutorials.
import numpy as np
import ufl
import dolfinx
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot
from dolfinx import nls
from dolfinx import log
from dolfinx import geometry
from dolfinx.plot import create_vtk_mesh
import pyvista
pyvista.set_jupyter_backend("panel")
L = 20
def create_mesh(L, H, N_DIV_x, N_DIV_y):
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0,0]), np.array([L, H])],
[N_DIV_x, N_DIV_y], mesh.CellType.quadrilateral)
return domain
def left(x):
return np.isclose(x[0], 0)
def right(x):
return np.isclose(x[0], L)
def solve_inplace(domain, V,u, E, nu, LINEAR = False):
fdim = domain.topology.dim -1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full(len(left_facets), 1, dtype=np.int32), np.full(len(right_facets), 2, dtype=np.int32)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
u_bc = np.array((0,) * domain.geometry.dim, dtype=PETSc.ScalarType)
left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.indices[facet_tag.values==1])
bcs = [fem.dirichletbc(u_bc, left_dofs, V)]
B = fem.Constant(domain, PETSc.ScalarType((0, 0)))
T = fem.Constant(domain, PETSc.ScalarType((0, 0)))
v = ufl.TestFunction(V)
# u = fem.Function(V)
# Spatial dimension
d = len(u)
# Identity tensor
I = ufl.variable(ufl.Identity(d))
# Deformation gradient
F = ufl.variable(I + ufl.grad(u))
# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)
# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))
# Elasticity parameters
E = PETSc.ScalarType(E)
nu = PETSc.ScalarType(nu)
mu = fem.Constant(domain, E/(2*(1 + nu)))
lmbda = fem.Constant(domain, E*nu/((1 + nu)*(1 - 2*nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
# Stress
# Hyper-elasticity
P = ufl.diff(psi, F)
if LINEAR:
P = 2.0 * mu * ufl.sym(ufl.grad(u)) + lmbda * ufl.tr(ufl.sym(ufl.grad(u))) * I
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
# Define form F (we want to find u such that F(u) = 0)
F = ufl.inner(ufl.grad(v), P)*dx - ufl.inner(v, B)*dx - ufl.inner(v, T)*ds(2)
#F = ufl.inner(ufl.grad(v), P)*dx - ufl.inner(v, B)*dx
problem = fem.petsc.NonlinearProblem(F, u, bcs)
solver = nls.petsc.NewtonSolver(domain.comm, problem)
# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"
log.set_log_level(log.LogLevel.INFO)
# log.set_log_level(None)
tval0 = -1
u_hist = []
for n in range(1, 10):
T.value[1] = n * tval0
num_its, converged = solver.solve(u)
assert(converged)
u.x.scatter_forward()
print(f"Time step {n}, Number of iterations {num_its}, Load {T.value}")
# plot_function(n, u)
u_hist.append(u.copy())
# problem.destroy()
return u_hist
### FUNCTION THAT SEEMINGLY CAUSES THE ISSUE
def extract_deflection(domain, uh):
bb_tree = geometry.BoundingBoxTree(domain, domain.topology.dim)
points = np.array([L,0.,0.]).reshape(1,-1)
cells = []
points_on_proc = []
# Find cells whose bounding-box collide with the the points
cell_candidates = geometry.compute_collisions(bb_tree, points.T)
colliding_cells = geometry.compute_colliding_cells(domain, cell_candidates, points.T)
points_on_proc = []
for i, point in enumerate(points):
if len(colliding_cells.links(i))>0:
points_on_proc.append(point)
cells.append(colliding_cells.links(i)[0])
points_on_proc = np.array(points_on_proc, dtype=np.float64)
u_values = uh.eval(points_on_proc, cells)
# displ_max.append(u_values)
return u_values