Error reported for solving thermoelastic problem

[RokyC]

I recently referred to weakly coupled implementations of thermoelastic problems on Linear thermoelasticity (weak coupling) — Numerical tours of continuum mechanics using FEniCS master documentation (comet-fenics.readthedocs.io). The temperature field calculations are working fine,but Elastic problem calculation broke down. Ask for the possible reasons.

(Error: error code 10
[0] KSPSolve() at /home/conda/feedstock_root/build_artifacts/petsc_1661243737224/work/src/ksp/ksp/interface/itfunc.c:1078
[0] KSPSolve_Private() at /home/conda/feedstock_root/build_artifacts/petsc_1661243737224/work/src/ksp/ksp/interface/itfunc.c:843
[0] KSPSetUp() at /home/conda/feedstock_root/build_artifacts/petsc_1661243737224/work/src/ksp/ksp/interface/itfunc.c:407
[0] PCSetUp() at /home/conda/feedstock_root/build_artifacts/petsc_1661243737224/work/src/ksp/pc/interface/precon.c:990
[0] PCSetUp_LU() at /home/conda/feedstock_root/build_artifacts/petsc_1661243737224/work/src/ksp/pc/impls/factor/lu/lu.c:93
[0] MatLUFactorSymbolic() at /home/conda/feedstock_root/build_artifacts/petsc_1661243737224/work/src/mat/interface/matrix.c:2966
[0] MatLUFactorSymbolic_SeqAIJ() at /home/conda/feedstock_root/build_artifacts/petsc_1661243737224/work/src/mat/impls/aij/seq/aijfact.c:361
[0] PetscFreeSpaceGet() at /home/conda/feedstock_root/build_artifacts/petsc_1661243737224/work/src/mat/utils/freespace.c:10
[0] PetscMallocA() at /home/conda/feedstock_root/build_artifacts/petsc_1661243737224/work/src/sys/memory/mal.c:401
[0] PetscMallocAlign() at /home/conda/feedstock_root/build_artifacts/petsc_1661243737224/work/src/sys/memory/mal.c:49
[0] PetscFreeSpaceGet
) 

Specifically as follows,The grid used is a hexahedral grid:

Surface,yp = 1, 4  
top = 3   
bottom = 2 
mesh,subdomains,boundaries = read_mesh(root='./Mesh')
T0 = 293.0 
E = 209e3 
nu = 0.3 
lmbda = E*nu/((1+nu)*(1-2*nu))
mu = E/2/(1+nu)
rho = 8030.0e-12  
alpha = 1.84e-5  
kappa = alpha*(2*mu + 3*lmbda)
cp = 502.48e3  
k = 16.27e-3  
Vte = dolfinx.fem.FunctionSpace(mesh, ("CG", 1))
Vue = dolfinx.fem.VectorFunctionSpace(mesh, ("CG", 2)) 
V_von_mises = dolfinx.fem.FunctionSpace(mesh, ("CG", 1))  
dx = ufl.Measure("dx")(subdomain_data=subdomains)
ds = ufl.Measure("ds")(subdomain_data=boundaries)
dt = dolfinx.fem.Constant(mesh,c=0.)
theta_trial = ufl.TrialFunction(Vte)
theta_test = ufl.TestFunction(Vte)
thetaold = dolfinx.fem.Function(Vte)
thetaold.x.array[:] = T0
thetaold.name = "Tempreture"
thetanew = dolfinx.fem.Function(Vte)
h_robin = dolfinx.fem.Function(Vte)
t_wall = dolfinx.fem.Function(Vte)  
mecha_trial = ufl.TrialFunction(Vue)
mecha_test = ufl.TestFunction(Vue)
mechaold = dolfinx.fem.Function(Vue)
mechaold.x.array[:] = T0
mechaold.name = "displacement"
mechanew = dolfinx.fem.Function(Vue)
stresses = dolfinx.fem.Function(V_von_mises)
stresses.name = "VonmiseS"

therm_forml = ufl.inner(cp*rho*theta_trial / dt,theta_test)*dx + k*ufl.dot(ufl.grad(theta_trial), ufl.grad(theta_test))*dx
therm_formr = ufl.inner(cp*rho*(thetaold) / dt,theta_test)*dx - ufl.inner(h_robin*(t_wall-thetaold),theta_test)*ds(Surface)

a_cpp_therm = dolfinx.fem.form(therm_forml)
f_cpp_therm = dolfinx.fem.form(therm_formr)
pCoord = ufl.SpatialCoordinate(mesh)  
sigma_inner = 2.0*mu*ufl.sym(ufl.grad(mecha_trial))+lmbda*ufl.nabla_div(mecha_trial)*ufl.Identity(len(mecha_trial)) # ufl.sym(ufl.grad(mecha_trial)) Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
therm_sigma = kappa*(thetanew-thetaold)*ufl.Identity(len(mecha_trial))
mecha_forml = ufl.inner(sigma_inner, ufl.sym(ufl.grad(mecha_test)))*dx
mecha_formr = ufl.inner(therm_sigma, ufl.grad(mecha_test))*dx + (ufl.inner((pCoord[0]), mecha_test[0])+ufl.inner(pCoord[2], mecha_test[2]))*dx
facets_FixSurface = boundaries.indices[boundaries.values == bottom] 
bdofs_FixSurface = dolfinx.fem.locate_dofs_topological(Vue, mesh.topology.dim - 1, facets_FixSurface) 
bc = dolfinx.fem.dirichletbc(np.zeros(3, dtype=petsc4py.PETSc.ScalarType), bdofs_FixSurface, Vue)

....
for (i, dti) in enumerate(np.diff(t)):
    print("Increment " + str(i+1))
    dt.value = dti
    h_robin.x.array[bdofsx_Surface] =[....]
    t_wall.x.array[bdofsx_Surface] = [....]
    A_therm = dolfinx.fem.petsc.assemble_matrix(a_cpp_therm) 
    A_therm.assemble() 
    F_therm = dolfinx.fem.petsc.assemble_vector(f_cpp_therm)  
    F_therm.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE) 
    solve(A=A_therm,F=F_therm,X=thetanew,mesh=mesh)
    A_mecha = dolfinx.fem.petsc.assemble_matrix(a_cpp_mecha, bcs=[bc])
    A_mecha.assemble()
    F_mecha = dolfinx.fem.petsc.assemble_vector(f_cpp_mecha)
    dolfinx.fem.petsc.apply_lifting(F_mecha, [a_cpp_mecha], [bcs]) 
    F_mecha.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)
    dolfinx.fem.petsc.set_bc(F_mecha, bcs)  
    solve(A=A_mecha,F=F_mecha,X=mechanew,mesh=mesh)

def solve(A,F,X,mesh):
    ksp = petsc4py.PETSc.KSP()
    ksp.create(mesh.comm)
    ksp.setOperators(A)        
    ksp.setType("preonly")  
    ksp.getPC().setType("lu")
    ksp.getPC().setFactorSolverType("mumps")
    ksp.setFromOptions()
    ksp.solve(F, X.vector)
    X.x.scatter_forward()

[dokken]

@RokyC

  A_therm = dolfinx.fem.petsc.assemble_matrix(a_cpp_therm) 
    A_therm.assemble() 
    F_therm = dolfinx.fem.petsc.assemble_vector(f_cpp_therm)  
    F_therm.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE) 
    solve(A=A_therm,F=F_therm,X=thetanew,mesh=mesh)
    A_mecha = dolfinx.fem.petsc.assemble_matrix(a_cpp_mecha, bcs=[bc])
    A_mecha.assemble()
    F_mecha = dolfinx.fem.petsc.assemble_vector(f_cpp_mecha)
    dolfinx.fem.petsc.apply_lifting(F_mecha, [a_cpp_mecha], [bcs]) 
    F_mecha.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)
    dolfinx.fem.petsc.set_bc(F_mecha, bcs)  
    solve(A=A_mecha,F=F_mecha,X=mechanew,mesh=mesh)

def solve(A,F,X,mesh):
    ksp = petsc4py.PETSc.KSP()
    ksp.create(mesh.comm)
    ksp.setOperators(A)        
    ksp.setType("preonly")  
    ksp.getPC().setType("lu")
    ksp.getPC().setFactorSolverType("mumps")
    ksp.setFromOptions()
    ksp.solve(F, X.vector)
    X.x.scatter_forward()

You should refactor this code such that you do not recreate matrices and solvers inside loops, as petsc has changed its garbage collection.

Manipulating vector data of fem.Function

This behavior is not safe, as dolfinx.fem.Function(..).vector does provide a view into data stored in the dolfinx.fem.Function(Vh). This data (and the underlying vector) is deleted when dl.fem.Functon(Vh) goes out of scope, which is right away. (This is due to PETSc’s modern garbage collection (dolfinx/function.py at main · FEniCS/dolfinx · GitHub, garbage collection different between petsc4py 3.17.4 and 3.18.3 (#1309) · Issues · PETSc / petsc · GitLab). Thus, if you want to use the petsc-vect…

[RokyC]

I abandon the loop and try to calculate a time step alone. The code still crashes when running to solve the elastic problem.I wonder if dolfinx’s variational form syntax is different from dolfin’s. Or my grid is a hexahedral grid that causes an error in the calculation of thermal stress by invoking the temperature field.

19:55:15.335 [error] Disposing session as kernel process died ExitCode: undefined, Reason: 
19:55:15.346 [info] Dispose Kernel process 26972.
19:55:15.349 [error] Raw kernel process exited code: undefined
19:55:15.390 [error] Error in waiting for cell to complete [Error: Canceled future for execute_request message before replies were done
	at t.KernelShellFutureHandler.dispose (~/.vscode/extensions/ms-toolsai.jupyter-2023.6.1001652307-linux-x64/out/extension.node.js:2:32375)
	at ~/.vscode/extensions/ms-toolsai.jupyter-2023.6.1001652307-linux-x64/out/extension.node.js:2:51427
	at Map.forEach (<anonymous>)
	at y._clearKernelState (~/.vscode/extensions/ms-toolsai.jupyter-2023.6.1001652307-linux-x64/out/extension.node.js:2:51412)
	at y.dispose (~/.vscode/extensions/ms-toolsai.jupyter-2023.6.1001652307-linux-x64/out/extension.node.js:2:44894)
	at ~/.vscode/extensions/ms-toolsai.jupyter-2023.6.1001652307-linux-x64/out/extension.node.js:24:111857
	at re (~/.vscode/extensions/ms-toolsai.jupyter-2023.6.1001652307-linux-x64/out/extension.node.js:2:1587256)
	at my.dispose (~/.vscode/extensions/ms-toolsai.jupyter-2023.6.1001652307-linux-x64/out/extension.node.js:24:111833)
	at vy.dispose (~/.vscode/extensions/ms-toolsai.jupyter-2023.6.1001652307-linux-x64/out/extension.node.js:24:119116)
	at process.processTicksAndRejections (node:internal/process/task_queues:96:5)]
19:55:15.399 [warn] Cell completed with errors {
  message: 'Canceled future for execute_request message before replies were done'
}
19:55:15.414 [info] End cell 3 execution @ 1687175715409, started @ 1687175665722, elapsed time = 49.687s
19:55:15.417 [warn] Cancel all remaining cells due to cancellation or failure in execution
19:55:15.739 [info] End cell 3 execution @ undefined, started @ undefined, elapsed time = 0s

The code is as follows

Vte = dolfinx.fem.FunctionSpace(mesh, ("CG", 1)) 
Vue = dolfinx.fem.VectorFunctionSpace(mesh, ("CG", 2)) 
V_von_mises = dolfinx.fem.FunctionSpace(mesh, ("CG", 1)) 

dx = ufl.Measure("dx")(subdomain_data=subdomains)
ds = ufl.Measure("ds")(subdomain_data=boundaries)
dt = dolfinx.fem.Constant(mesh,c=0.)

theta_trial = ufl.TrialFunction(Vte)
theta_test = ufl.TestFunction(Vte)
thetaold = dolfinx.fem.Function(Vte)
thetaold.x.array[:] = T0
thetaold.name = "Tempreture"
thetanew = dolfinx.fem.Function(Vte)
h_robin = dolfinx.fem.Function(Vte)  
t_wall = dolfinx.fem.Function(Vte)  
mecha_trial = ufl.TrialFunction(Vue)
mecha_test = ufl.TestFunction(Vue)
mechaold = dolfinx.fem.Function(Vue)
mechaold.x.array[:] = 0
mechaold.name = "displacement"
mechanew = dolfinx.fem.Function(Vue)
stresses = dolfinx.fem.Function(V_von_mises)
stresses.name = "VonmiseS"

therm_forml = ufl.inner(cp*rho*theta_trial / dt,theta_test)*dx + k*ufl.dot(ufl.grad(theta_trial), ufl.grad(theta_test))*dx
therm_formr = ufl.inner(cp*rho*(thetaold) / dt,theta_test)*dx - ufl.inner(h_robin*(t_wall-thetaold),theta_test)*ds(Surface)

a_cpp_therm = dolfinx.fem.form(therm_forml)
f_cpp_therm = dolfinx.fem.form(therm_formr)
....
file_results=dolfinx.io.VTKFile(MPI.COMM_WORLD, "./Temp.pvd", "w")
dt.value = t[1]-t[0]
X_predict = scaler.transform(par_t[1,:].reshape(1, -1))
pre_HCof = Rom1.predict(X_predict)
pre_Tjr = Rom2.predict(X_predict)
heatCof = scipy.interpolate.griddata(inter_coord,pre_HCof,dof_coordinates,fill_value=True,method='nearest')
Twall = scipy.interpolate.griddata(inter_coord,pre_Tjr,dof_coordinates,fill_value=True,method='nearest')

h_robin.x.array[bdofsx_Surface] = -heatCof
t_wall.x.array[bdofsx_Surface] = Twall

A_therm = dolfinx.fem.petsc.assemble_matrix(a_cpp_therm) 
A_therm.assemble() 
F_therm = dolfinx.fem.petsc.assemble_vector(f_cpp_therm)
F_therm.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE) 
solve(A=A_therm,F=F_therm,X=thetanew,mesh=mesh)

sigma_inner = 2.0*mu*ufl.sym(ufl.grad(mecha_trial))+lmbda*ufl.nabla_div(mecha_trial)*ufl.Identity(len(mecha_trial)) # ufl.sym(ufl.grad(mecha_trial)) Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
mecha_forml = ufl.inner(sigma_inner, ufl.sym(ufl.grad(mecha_test)))*dx
mecha_formr = ufl.inner(kappa*(thetanew-293.0)*ufl.Identity(len(mecha_trial)), ufl.grad(mecha_test))*dx

facets_Surface = boundaries.indices[boundaries.values == Surface] 
bdofsx_Surface = dolfinx.fem.locate_dofs_topological(Vte, mesh.topology.dim - 1, facets_Surface) 
dof_coordinates = Vte.tabulate_dof_coordinates()[bdofsx_Surface] 
zero_scalar = petsc4py.PETSc.ScalarType(0)
facets_FixSurface = boundaries.indices[boundaries.values == bottom]  
bdofsx_FixSurface = dolfinx.fem.locate_dofs_topological(Vue.sub(0), mesh.topology.dim - 1, facets_FixSurface)
bcsx = dolfinx.fem.dirichletbc(zero_scalar, bdofsx_FixSurface, Vue.sub(0))
bdofsz_FixSurface = dolfinx.fem.locate_dofs_topological(Vue.sub(2), mesh.topology.dim - 1, facets_FixSurface) 
bcsz = dolfinx.fem.dirichletbc(zero_scalar, bdofsx_FixSurface, Vue.sub(2))
bc = [bcsx, bcsz]

a_cpp_mecha = dolfinx.fem.form(mecha_forml)
f_cpp_mecha = dolfinx.fem.form(mecha_formr)
A_mecha = dolfinx.fem.petsc.assemble_matrix(a_cpp_mecha, bcs=bc) 
A_mecha.assemble()  
F_mecha = dolfinx.fem.petsc.assemble_vector(f_cpp_mecha)
dolfinx.fem.petsc.apply_lifting(F_mecha, [a_cpp_mecha], [bc]) 
F_mecha.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(F_mecha, bc)  
solve(A=A_mecha,F=F_mecha,X=mechanew,mesh=mesh)

file_results.close()

[dokken]

You need to supply a reproducible code, this means that all variables (including the import statements, mesh declaration, and the mesh file) has to be included in the post, for anyone to reproduce your error message.

For instance, can you reproduce the error on a unit cube with hexahedral elements?

I would also try to run the code a pure python script, and not through notebooks, as it can swallow error messages.

[RokyC]

I tried a simple implementation that dolfinx could compute on a grid he generated himself. But using the mesh I imported from meshio(both are hexahedral grids.), dolfinx can only calculate the thermal equation, but not the elastic equation.Even using.py instead of jupyterlab will kill processes without error.
Here is the code for using dolfinx mesh

import numpy as np
import ufl
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import dolfinx
from dolfinx import mesh, fem, plot, io
import petsc4py

L = 1
W = 0.2
T0 = 293.0 
E = 209e3 
nu = 0.3 
lmbda = E*nu/((1+nu)*(1-2*nu))
mu = E/2/(1+nu)
rho = 8030.0e-12  
alpha = 1.84e-5  
kappa = alpha*(2*mu + 3*lmbda)
cp = 502.48e3 
k = 16.27e-3 

def clamped_boundary(x):
    return np.isclose(x[0], 0)
def boundary_D(x):
    return np.isclose(x[0], 1)

domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([L, W, W])],
                  [20,6,6], cell_type=mesh.CellType.hexahedron)

fdim = domain.topology.dim - 1
facets = mesh.locate_entities(domain, fdim, boundary_D)
facet_tag = mesh.meshtags(domain, fdim, facets, 1)
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)

Vte = dolfinx.fem.FunctionSpace(domain, ("CG", 1))
Vue = dolfinx.fem.VectorFunctionSpace(domain, ("CG", 2)) 
V_von_mises = dolfinx.fem.FunctionSpace(domain, ("CG", 1)) 
dt = dolfinx.fem.Constant(domain,c=0.)

theta_trial = ufl.TrialFunction(Vte)
theta_test = ufl.TestFunction(Vte)
thetaold = dolfinx.fem.Function(Vte)
thetaold.x.array[:] = T0
thetanew = dolfinx.fem.Function(Vte)
thetanew.name = "Tempreture"
h_robin = dolfinx.fem.Function(Vte)  
t_wall = dolfinx.fem.Function(Vte)  

mecha_trial = ufl.TrialFunction(Vue)
mecha_test = ufl.TestFunction(Vue)
mechaold = dolfinx.fem.Function(Vue)
mechaold.x.array[:] = 0
mechanew = dolfinx.fem.Function(Vue)
mechanew.name = "displacement"
stresses = dolfinx.fem.Function(V_von_mises)
stresses.name = "VonmiseS"

dx = ufl.dx
ds = ufl.Measure("ds", domain=domain)
# problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
# uh = problem.solve()
dofs_D = fem.locate_dofs_geometrical(Vte, boundary_D)
therm_forml = ufl.inner(cp*rho*theta_trial / dt,theta_test)*dx + k*ufl.dot(ufl.grad(theta_trial), ufl.grad(theta_test))*dx
therm_formr = ufl.inner(cp*rho*(thetaold) / dt,theta_test)*dx - ufl.inner(h_robin*(t_wall-thetaold),theta_test)*ds(subdomain_data=facet_tag)
a_cpp_therm = dolfinx.fem.form(therm_forml)
f_cpp_therm = dolfinx.fem.form(therm_formr)
sigma_inner = 2.0*mu*ufl.sym(ufl.grad(mecha_trial))+lmbda*ufl.nabla_div(mecha_trial)*ufl.Identity(len(mecha_trial)) 
mecha_forml = ufl.inner(sigma_inner, ufl.sym(ufl.grad(mecha_test)))*dx
mecha_formr = ufl.inner(kappa*(thetanew-293.0)*ufl.Identity(len(mecha_trial)), ufl.grad(mecha_test))*dx
a_cpp_mecha = dolfinx.fem.form(mecha_forml)
f_cpp_mecha = dolfinx.fem.form(mecha_formr)

dt.value = 1
h_robin.x.array[dofs_D] = -0.00056841
t_wall.x.array[dofs_D] = 1365

A_therm = dolfinx.fem.petsc.assemble_matrix(a_cpp_therm) 
A_therm.assemble()  
F_therm = dolfinx.fem.petsc.assemble_vector(f_cpp_therm)
F_therm.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE) 
ksp = petsc4py.PETSc.KSP()
ksp.create(domain.comm)
ksp.setOperators(A_therm)        
ksp.setType("preonly")    
ksp.getPC().setType("lu")  
ksp.getPC().setFactorSolverType("mumps")
ksp.setFromOptions()
ksp.solve(F_therm, thetanew.vector)
thetanew.x.scatter_forward()
file_results=dolfinx.io.VTKFile(MPI.COMM_WORLD, "./Temp.pvd", "w")
file_results.write_function(thetanew, 1)
file_results.close()

# fdim = domain.topology.dim - 1
# boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
# u_D = np.array([0,0,0], dtype=ScalarType)
# bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(Vue, fdim, boundary_facets), Vue)

# A_mecha = dolfinx.fem.petsc.assemble_matrix(a_cpp_mecha, bcs=[bc]) 
A_mecha = dolfinx.fem.petsc.assemble_matrix(a_cpp_mecha)
A_mecha.assemble()
F_mecha = dolfinx.fem.petsc.assemble_vector(f_cpp_mecha)
# dolfinx.fem.petsc.apply_lifting(F_mecha, [a_cpp_mecha], [[bc]])
F_mecha.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)
# dolfinx.fem.petsc.set_bc(F_mecha, [bc])
ksp = petsc4py.PETSc.KSP()
ksp.create(domain.comm)
ksp.setOperators(A_mecha)
ksp.setType("preonly")    
ksp.getPC().setType("lu")  
ksp.getPC().setFactorSolverType("mumps")
ksp.setFromOptions()
ksp.solve(F_mecha, mechanew.vector)
mechanew.x.scatter_forward()

file_results=dolfinx.io.VTKFile(MPI.COMM_WORLD, "./Mecha.pvd", "w")
file_results.write_function(mechanew, 1)
file_results.close()

Here is a program that uses the meshio grid

import dolfinx.fem, dolfinx.mesh, dolfinx.io, dolfinx.plot
import numpy as np
import petsc4py.PETSc
import ufl
from mpi4py import MPI
import os

def read_mesh(root=''):
    with dolfinx.io.XDMFFile(MPI.COMM_WORLD, os.path.join(root, 'mesh.xdmf'), 'r') as xdmf:
        mesh = xdmf.read_mesh(name="Grid")  
        subdomains = xdmf.read_meshtags(mesh, name="Grid")
    mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)  
    # mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
    with dolfinx.io.XDMFFile(MPI.COMM_WORLD, os.path.join(root, 'mt.xdmf'), 'r') as xdmf:
        boundaries = xdmf.read_meshtags(mesh, name="Grid")
    
    return mesh,subdomains,boundaries

Surface,yp = 1, 4  
top = 3    
bottom = 2 

mesh,subdomains,boundaries = read_mesh(root='./Mesh')

T0 = 293.0 
E = 209e3 
nu = 0.3 
lmbda = E*nu/((1+nu)*(1-2*nu))
mu = E/2/(1+nu)
rho = 8030.0e-12  
alpha = 1.84e-5  
kappa = alpha*(2*mu + 3*lmbda)
cp = 502.48e3  
k = 16.27e-3  

Vte = dolfinx.fem.FunctionSpace(mesh, ("CG", 1)) 
Vue = dolfinx.fem.VectorFunctionSpace(mesh, ("CG", 2)) 
V_von_mises = dolfinx.fem.FunctionSpace(mesh, ("CG", 1)) 

dx = ufl.Measure("dx")(subdomain_data=subdomains)
ds = ufl.Measure("ds")(subdomain_data=boundaries)

theta_trial = ufl.TrialFunction(Vte)
theta_test = ufl.TestFunction(Vte)
thetaold = dolfinx.fem.Function(Vte)
thetaold.x.array[:] = T0
thetanew = dolfinx.fem.Function(Vte)
thetanew.name = "Tempreture"
h_robin = dolfinx.fem.Function(Vte)  
t_wall = dolfinx.fem.Function(Vte)  
dt = dolfinx.fem.Constant(mesh,c=0.)

mecha_trial = ufl.TrialFunction(Vue)
mecha_test = ufl.TestFunction(Vue)
mechaold = dolfinx.fem.Function(Vue)
mechaold.x.array[:] = 0
mechanew = dolfinx.fem.Function(Vue)
mechanew.name = "displacement"
stresses = dolfinx.fem.Function(V_von_mises)
stresses.name = "VonmiseS"

# problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
# uh = problem.solve()

therm_forml = ufl.inner(cp*rho*theta_trial / dt,theta_test)*dx + k*ufl.dot(ufl.grad(theta_trial), ufl.grad(theta_test))*dx
therm_formr = ufl.inner(cp*rho*(thetaold) / dt,theta_test)*dx - ufl.inner(h_robin*(t_wall-thetaold),theta_test)*ds(Surface)
a_cpp_therm = dolfinx.fem.form(therm_forml)
f_cpp_therm = dolfinx.fem.form(therm_formr)

sigma_inner = 2.0*mu*ufl.sym(ufl.grad(mecha_trial))+lmbda*ufl.nabla_div(mecha_trial)*ufl.Identity(len(mecha_trial)) 
mecha_forml = ufl.inner(sigma_inner, ufl.sym(ufl.grad(mecha_test)))*dx
mecha_formr = ufl.inner(kappa*(thetanew-T0)*ufl.Identity(len(mecha_trial)), ufl.grad(mecha_test))*dx
a_cpp_mecha = dolfinx.fem.form(mecha_forml)
f_cpp_mecha = dolfinx.fem.form(mecha_formr)

T3Surf_index = boundaries.indices[boundaries.values == Surface]  
T3Surf_bdofsT = dolfinx.fem.locate_dofs_topological(Vte, mesh.topology.dim - 1, T3Surf_index) 
dof_coordinates = Vte.tabulate_dof_coordinates()[T3Surf_bdofsT]  

dt.value = 1
h_robin.x.array[T3Surf_bdofsT] = -0.00056841
t_wall.x.array[T3Surf_bdofsT] = 1365

A_therm = dolfinx.fem.petsc.assemble_matrix(a_cpp_therm) 
A_therm.assemble()  
F_therm = dolfinx.fem.petsc.assemble_vector(f_cpp_therm)
F_therm.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE) 
ksp = petsc4py.PETSc.KSP()
ksp.create(mesh.comm)
ksp.setOperators(A_therm)        
ksp.setType("preonly")    
ksp.getPC().setType("lu")  
ksp.getPC().setFactorSolverType("mumps")
ksp.setFromOptions()
ksp.solve(F_therm, thetanew.vector)
thetanew.x.scatter_forward()
file_results=dolfinx.io.VTKFile(MPI.COMM_WORLD, "./Temp.pvd", "w")
file_results.write_function(thetanew, 1)
file_results.close()

# A_mecha = dolfinx.fem.petsc.assemble_matrix(a_cpp_mecha, bcs=[bc]) 
A_mecha = dolfinx.fem.petsc.assemble_matrix(a_cpp_mecha) 
A_mecha.assemble() 
F_mecha = dolfinx.fem.petsc.assemble_vector(f_cpp_mecha)
# dolfinx.fem.petsc.apply_lifting(F_mecha, [a_cpp_mecha], [[bc]]) 
F_mecha.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)
# dolfinx.fem.petsc.set_bc(F_mecha, [bc])  
ksp = petsc4py.PETSc.KSP()
ksp.create(mesh.comm)
ksp.setOperators(A_mecha)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
ksp.setFromOptions()
ksp.solve(F_mecha, mechanew.vector)
mechanew.x.scatter_forward()

file_results=dolfinx.io.VTKFile(MPI.COMM_WORLD, "./Mecha.pvd", "w")
file_results.write_function(mechanew, 1)
file_results.close()

In order to reproduce the problem，what should I do to upload the grid file of ‘mesh.xdmf,mesh.h5,mt.xdmf,mt.h5’?

[dokken]

You Need to either:

1. convert them to plain-text (ascii) xdmf and upload them here
2. share them through Google drive.

@RokyC

Even using.py instead of jupyterlab will kill processes without error.

what line exactly does the code crash in?

[RokyC]

Grid file sharing link below
(Mesh - Google Drive)

[RokyC]

Crashing occurred while running to this line:ksp.solve(F_mecha, mechanew.vector), and the system did not respond after the crash.I tried to change the variational form again (not involving the temperature problem calculation results) and the same problem still occurred

[RokyC]

Oh, I have found the reason. The memory space is insufficient due to the number of grids. I am very sorry for the trouble these days, and I will continue to learn to optimize the memory usage of the code