3D elasticity model with a big mesh

[DiegoRojas]

Hello, nice to meet you. I hope you are well. I am trying to run this code to solve the linear equation of elasticity in 3d and with it I can get the displacements, stress tensor and principal stresses. My .xml file named “Mesh.all3.xml” has the following characteristics:

Expecting 1531431 vertices
Found all vertices
Expecting 7464768 cells
Found all cells
779 megabyte (MB)

Unfortunately when I run the code I get the following error:
Error: Unable to successfully call PETSc function “KSPsolve”.
Reason: PETSc error code is: 83 ((null)).

I attach an image of the error in the code, which tells me it is in row 121. solver.solve(u.vector(), b);

Below I attach my code:´

from IPython.display import clear_output 
try:
    import dolfin
except ImportError:
    !wget "https://fem-on-colab.github.io/releases/fenics-install.sh" -O "/tmp/fenics-install.sh" && bash "/tmp/fenics-install.sh"
    import dolfin

import matplotlib.pyplot as plt
import numpy as np
import pickle
import cmath
import os, time
from ufl import nabla_div

clear_output();

from dolfin import *

mesh = Mesh("Mesh.all3.xml")


#Scaled Variables
E = Constant(20e3) #Módulo de Young
nu = Constant(0.3) #Radio de Poisson
mu = E/(2*(1+nu)) #Módulo de cizalle
lambda_ = E*nu/((1+nu)*(1-2*nu)) #Lamé
rho = Constant(0.027) #densidad

#Strain function
def epsilon(u):
    return sym(nabla_grad(u))

#Stress function
def sigma(u):
    return lambda_*nabla_div(u)*Identity(3) + 2*mu*epsilon(u)

#lambda is a reserved python keyword, naming convention recommends using a single trailing #underscore for such cases

#Load
b_z = -rho
b = Constant((0.0,0.0,b_z))

#Function Spaces
V = VectorFunctionSpace(mesh, "CG", 1)
u_tr = TrialFunction(V)
u_test = TestFunction(V)

a = inner(sigma(u_tr),epsilon(u_test))*dx #bilineal form
l = dot(b,u_test)*dx #lineal form

tol=10e-10
def bottom(x, on_boundary):
    return (on_boundary and near(x[2],0,tol))
bc1 =DirichletBC(V.sub(2), Constant(0.),bottom)

def left(x, on_boundary):
    return (on_boundary and near(x[0], 0,tol))
bc2 =DirichletBC(V.sub(0), Constant(0.),left)

def right(x, on_boundary):
    return (on_boundary and near(x[0], 5760,tol))
bc3 =DirichletBC(V.sub(0), Constant(0.),right)

def north(x, on_boundary):
    return (on_boundary and near(x[1], 5384,tol))
bc4 =DirichletBC(V.sub(1), Constant(0.),north)

def south(x, on_boundary):
    return (on_boundary and near(x[1], 0,tol))
bc5 =DirichletBC(V.sub(1), Constant(0.),south)

bc=[bc1,bc2,bc3,bc4,bc5]

parameters["linear_algebra_backend"] = "PETSc"

def build_nullspace(V, x):
    #Function to build null space for 3D elasticity

    #Create list of vectors for null space
    nullspace_basis = [x.copy() for i in range(6)]

    #Build translational null space basis
    V.sub(0).dofmap().set(nullspace_basis[0], 1.0);
    V.sub(1).dofmap().set(nullspace_basis[1], 1.0);
    V.sub(2).dofmap().set(nullspace_basis[2], 1.0);

    #Build rotational null space basis
    V.sub(0).set_x(nullspace_basis[3], -1.0, 1);
    V.sub(1).set_x(nullspace_basis[3],  1.0, 0);
    V.sub(0).set_x(nullspace_basis[4],  1.0, 2);
    V.sub(2).set_x(nullspace_basis[4], -1.0, 0);
    V.sub(2).set_x(nullspace_basis[5],  1.0, 1);
    V.sub(1).set_x(nullspace_basis[5], -1.0, 2);

    for x in nullspace_basis:
        x.apply("insert")

    #Create vector space basis and orthogonalize
    basis = VectorSpaceBasis(nullspace_basis)
    basis.orthonormalize()

    return basis

#Assemble system, applying boundary conditions and preserving
#symmetry)
A, b = assemble_system(a, l, bc)

#SOLVE
u = Function(V, name="Desplazamiento")

#Create near null space basis (required for smoothed aggregation
#AMG). The solution vector is passed so that it can be copied to
#generate compatible vectors for the nullspace.
null_space = build_nullspace(V, u.vector())

#Attach near nullspace to matrix
as_backend_type(A).set_near_nullspace(null_space)

#Create PETSC smoothed aggregation AMG preconditioner and attach near
#null space
pc = PETScPreconditioner("petsc_amg")

#Use Chebyshev smoothing for multigrid
PETScOptions.set("mg_levels_ksp_type", "chebyshev")
PETScOptions.set("mg_levels_pc_type", "jacobi")

#Improve estimate of eigenvalues for Chebyshev smoothing
PETScOptions.set("mg_levels_esteig_ksp_type", "cg")
PETScOptions.set("mg_levels_ksp_chebyshev_esteig_steps", 50)

#Create CG Krylov solver and turn convergence monitoring on
solver = PETScKrylovSolver("cg", pc)
solver.parameters["monitor_convergence"] = True

#Set matrix operator
solver.set_operator(A);

#Compute solution
solver.solve(u.vector(), b);

#Save solution to VTK format
File("elasticity.pvd", "compressed") << u

Vsig = TensorFunctionSpace(mesh, "DG",0)
sig = Function(Vsig, name="Cauchy_Stress")
sig.assign(project(-sigma(u),V=Vsig, solver_type="cg"));
file_results = XDMFFile("Results.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)
file_results.write(sig, 0.)

What could be the cause of the error?
Could you please recommend me some linear solver for this code?

Thank you very much for your replies. Greetings to the community.

Here is a picture of the error in the code: