Hi all,
Could I ask the factors causing kernel crash in dolfinx, and the elements number is increased from [4 \times 4 \times 4] to [10 \times 10 \times10] in a 3D domain, then it shows that the kernel crashes in the following information:
I checked the relavent information here, but didn’t fix it. Does anyone could help me with that?
Thanks a lot
Without having the specific code at hand, or any information about how much memory (RAM) that you have available, it is very hard to give you any guidance.
Hi dokken, the code is:
import numpy as np
import matplotlib.pyplot as plt
from dolfinx import mesh, fem, io, plot, default_scalar_type, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
import ufl
import gmsh
import os
from tqdm import tqdm
import csv
import sys
lambda_ = 400889.871e6 ## unit is Pa
mu = 80.194e6
nu = lambda_/(2*(lambda_ + mu)) ## 3D domain
E = mu*(3*lambda_ + 2*mu)/(lambda_ + mu) ## 3D domain
### (1a) define the 3D domain
coords_lower_left = [0.0, 0.0, 0.0]
Lx, Ly, Lz = 1.0e-3, 1.0e-3, 1.0e-3
coords_upper_right = [Lx, Ly, Lz]
nx, ny, nz = 8, 8, 8
num_elements = [nx, ny, nz]
domain = mesh.create_box(MPI.COMM_WORLD, [coords_lower_left, coords_upper_right], num_elements, mesh.CellType.hexahedron)
### (1b) set element types
element_u = ufl.VectorElement("CG", domain.ufl_cell(), degree=2) ## vector: displacement
element_p= ufl.FiniteElement("CG", domain.ufl_cell(), degree=1) ## scalar: pressure
element_mix = ufl.MixedElement([element_u, element_p])
# element_mix = mixed_element([element_u, element_p], gdim=3)
V = fem.FunctionSpace(domain, element_mix)
V_disp = V.sub(0)
V_press = V.sub(1)
num_subs = V.element.num_sub_elements
spaces = []
maps = []
for i in range(num_subs):
space_i, map_i = V.sub(i).collapse()
spaces.append(space_i)
maps.append(map_i)
### displacement component in x, y, z direction
V_disp_x = V_disp.sub(0)
V_disp_y = V_disp.sub(1)
V_disp_z = V_disp.sub(2)
### define surface
def front(x):
return np.isclose(x[0], Lx, atol=1.0e-8)
def right(x):
return np.isclose(x[1], Ly, atol=1.0e-8) ## y=Ly
def bottom(x):
return np.isclose(x[2], 0, atol=1.0e-8) ## z=0
def top(x):
return np.isclose(x[2], Lz, atol=1.0e-8) ## z=Lz
### sub-domain surface where traction force is applied
def Omega_traction(x):
return np.logical_and(np.logical_and(x[0]>=0.5*Lx, x[1]>=0.5*Ly), np.isclose(x[2],Lz, atol=1e-8))
### constraint displacement component on surface
fdim = domain.topology.dim - 1
value_zero= fem.Constant(domain, 0.0)
### front: ux=0
facet_front = mesh.locate_entities_boundary(domain, fdim, front)
dofs_front_x = fem.locate_dofs_topological(V_disp_x, fdim, facet_front)
bc_front_x = fem.dirichletbc(value_zero, dofs_front_x, V_disp_x)
### right: uy=0
facet_right = mesh.locate_entities_boundary(domain, fdim, right)
dofs_right_y = fem.locate_dofs_topological(V_disp_y, fdim, facet_right)
bc_right_y = fem.dirichletbc(value_zero, dofs_right_y, V_disp_y)
### bottom: uz=0
facet_bottom = mesh.locate_entities_boundary(domain, fdim, bottom)
dofs_bottom_z = fem.locate_dofs_topological(V_disp_z, fdim, facet_bottom)
bc_bottom_z = fem.dirichletbc(value_zero, dofs_bottom_z, V_disp_z)
### top: ux=0 & uy=0
facet_top = mesh.locate_entities_boundary(domain, fdim, top)
dofs_top_x = fem.locate_dofs_topological(V_disp_x, fdim, facet_top)
dofs_top_y = fem.locate_dofs_topological(V_disp_y, fdim, facet_top)
bc_top_x = fem.dirichletbc(value_zero, dofs_top_x, V_disp_x)
bc_top_y = fem.dirichletbc(value_zero, dofs_top_y, V_disp_y)
bcs = [bc_front_x, bc_right_y, bc_bottom_z, bc_top_x, bc_top_y]
metadata = {"quadrature_degree": 4}
### 2.2 apply the traction force
traction = fem.Constant(domain, default_scalar_type((0.0, 0.0, 0.0)))
facet_indices_Omega = mesh.locate_entities_boundary(domain, fdim, Omega_traction)
tag_Omega = 5
facet_tag_Omega = mesh.meshtags(domain, fdim, facet_indices_Omega, tag_Omega)
### trail and test function
state = fem.Function(V)
u, p = ufl.split(state)
# u, p = ufl.TrialFunctions(V)
v, q = ufl.TestFunctions(V)
#### Kinematics
d = len(u)
I = ufl.variable(ufl.Identity(d))
F = ufl.variable(I + ufl.grad(u))
C = ufl.variable(F.T * F)
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F)) ##
psi_dev = mu/2*(Ic - 3) - mu*ufl.ln(J)
psi_imp = p * ufl.ln(J) - (1/(2*lambda_))*p**2
psi = psi_dev + psi_imp
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
### test: apply the traction force on the top surface
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag_Omega, metadata=metadata)
PK1 = ufl.diff(psi, F)
a = ufl.inner(ufl.grad(v), PK1)*dx + q * (ufl.ln(J) - (1/lambda_)*p) * dx
L = ufl.dot(v, traction)* ds(tag_Omega)
R = a - L
problem = NonlinearProblem(R, state, bcs)
solver = NewtonSolver(domain.comm, problem)
### Set Newton solver options
solver.atol = 1e-6
solver.rtol = 1e-6
solver.convergence_criterion = "incremental"
# log.set_log_level(log.LogLevel.INFO)
steps = 60
traction_final = -60e6
traction_interval = traction_final/steps
for n in tqdm(range(0, steps+1)):
traction.value[2] = n * traction_interval
num_its, converged = solver.solve(state)
assert (converged)
u_sol, p_sol = state.split()
u_sol.x.scatter_forward() #### update the solution
p_sol.x.scatter_forward()
state.x.scatter_forward()
This is certainly not the code you have been running, as
is not valid.
Secondly, you have not made any specifications of the Krylov Subspace solver, meaning that you are using PETSc defaults.
See for instance: Implementation — FEniCSx tutorial on how to specify solvers.
Secondly, changing your problem from 8x8x8, to 10x10x10, it runs (4.58 s/it).
Hi dokken, I recorrect the above code and thanks for the tutorial about specifying solvers.
When I try to use the parameter for specifying the solver, the residual r(abs) = 0 in the first step, and the solution 0 everywhere, do you know why?
The code is quite long, so I just skimmed through it without reading it carefully.
It seems to me that you have homogeneous Dirichlet BCs and that your traction is zero. I am not sure for what reason you would expect the solution to be nonzero.
I applied the traction force by
steps = 20
traction_final = -320e6
traction_interval = traction_final/steps
for n in tqdm(range(0, steps+1)):
traction.value[2] = n * traction_interval
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(state)
assert (converged)
print(f"Number of interations: {n:d}")
assert (converged)
u_sol, p_sol = state.split()
u_sol.x.scatter_forward() #### update the solution
p_sol.x.scatter_forward()
state.x.scatter_forward()
and the solution is 0 every where.
the above code only works when the element meshing is 4x4x4, but kernel crashes when 8x8x8 or more due to the memory size of my computer, is specifying solver helpful to fix the crashing problem?
You only get a zero solution at the first step, when you apply a zero traction.
I ask again
How much memory do you have available?
16 GB RAM
8-cores M1 CPU
That should be sufficient for such a small problem.
Have you looked at the memory usage while running the code?
Additionally, could you try to run it as a Python script rather than within jupyter?