Error while doing modal analysis through parallelization

Hi all,
I am performing modal analysis of a 3D box with a distributed mass on its interior faces of cylindrical holes. The geometry and mesh files are here- geometry and mesh . It has DoF of 1.1 million. When I am running it through parallelization , it is giving an error-

===================================================================================
=   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
=   PID 456 RUNNING AT ec89dea7792c
=   EXIT CODE: 11
=   CLEANING UP REMAINING PROCESSES
=   YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
===================================================================================
YOUR APPLICATION TERMINATED WITH THE EXIT STRING: Segmentation fault (signal 11)
This typically refers to a problem with your application.
Please see the FAQ page for debugging suggestions

The MWE is-

from dolfin import *
import numpy as np
import json



# -------------------------------------------------------------------
# MPI functions -----------------------------------------------------
# -------------------------------------------------------------------
comm = MPI.comm_world
rank = MPI.rank(comm)


def mprint(*argv):
    if rank == 0:
        out = ""
        for arg in argv:
            out = out + str(arg)
        # this forces program to output when run in parallel
        print(out, flush=True)

# -------------------------------------------------------------------
# Set some common FEniCS flags --------------------------------------
# -------------------------------------------------------------------
set_log_level(50)
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True


mesh = Mesh()
with XDMFFile("mesh/tetra.xdmf") as infile:
    infile.read(mesh)

elem_P = VectorElement('Lagrange', mesh.ufl_cell(), 2)
elem_R = VectorElement('Real', mesh.ufl_cell(), 0)
elem = MixedElement([elem_P, elem_R, elem_R])
V = FunctionSpace(mesh, elem)

u_, U_, C_ = TrialFunctions(V)
du, dU, dC = TestFunctions(V)


E, nu = 200e9, 0.3  
rho = Constant(7850)

mu = Constant(E/2./(1+nu))
lmbda = Constant(E*nu/(1+nu)/(1-2*nu))

def eps(v):
    return sym(grad(v))
def sigma(v):
    dim = v.geometric_dimension()
    return 2.0*mu*eps(v) + lmbda*tr(eps(v))*Identity(dim)


support = CompiledSubDomain("x[2] == 1")
mf = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
mf.set_all(0)
support.mark(mf, 1)

bc = DirichletBC(V.sub(0), Constant((0.,0.,0.)), support)

force_face1 = CompiledSubDomain("(x[0] - 5e-2)*(x[0] - 5e-2) + (x[1] - 0.7)*(x[1] - 0.7) + (x[2] - 0.25043)*(x[2] - 0.25043)  < 0.15*0.15")
line_mf = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
line_mf.set_all(0)
force_face1.mark(line_mf, 2)


force_face2 = CompiledSubDomain("(x[0] - 5e-2)*(x[0] - 5e-2) + (x[1] - 0.5)*(x[1] - 0.5) + (x[2] - 0.51831)*(x[2] - 0.51831)  < 0.15*0.15")
force_face2.mark(line_mf, 3)


force_face3 = CompiledSubDomain("(x[0] - 5e-2)*(x[0] - 5e-2) + (x[1] - 0.3)*(x[1] - 0.3) + (x[2] - 0.83302)*(x[2] - 0.83302)  < 0.15*0.15")
force_face3.mark(line_mf, 4)


hole_radius1 = 0.1
plate_thk1 = 0.1
mass1 = 1   # kg
area1 = 2*np.pi*hole_radius1*plate_thk1
M1 = Constant(mass1)
A1 = Constant(area1)

hole_radius2 = 0.1
plate_thk2 = 0.1
mass2 = 1   # kg
area2 = 2*np.pi*hole_radius2*plate_thk2
M2 = Constant(mass2)
A2 = Constant(area2)

hole_radius3 = 0.1
plate_thk3 = 0.1
mass3 = 10   # kg
area3 = 2*np.pi*hole_radius3*plate_thk3
M3 = Constant(mass3)
A3 = Constant(area3)

ds = Measure('ds', mesh, subdomain_data=line_mf)
k_form = inner(sigma(du),eps(u_))*dx - (inner(C_, du-dU) + inner(u_-U_, dC))*ds(2) - (inner(C_, du-dU) + inner(u_-U_, dC))*ds(3) - (inner(C_, du-dU) + inner(u_-U_, dC))*ds(4) 
m_form = rho*inner(u_,du)*dx + M1/A1*inner(U_,du)*ds(2) + M1/A1*inner(U_,du)*ds(2) + M2/A2*inner(U_,du)*ds(3) + M3/A3*inner(U_,du)*ds(4)

K = PETScMatrix()
assemble(k_form, tensor=K)

M = PETScMatrix()
assemble(m_form, tensor=M)


bc.zero(M)
bc.apply(K)

# Solution
eigensolver = SLEPcEigenSolver(K, M)
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 0.

N_eig = 40   # number of eigenvalues
mprint("Computing {} first eigenvalues...".format(N_eig))
eigensolver.solve(N_eig)

file_results = XDMFFile("modal_analysis.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True

# Extraction
for i in range(N_eig):
    # Extract eigenpair
    r, c, rx, cx = eigensolver.get_eigenpair(i)
    
    # 3D eigenfrequency
    freq_3D = np.sign(r)*np.sqrt(np.abs(r))/2/pi
    mprint("Solid FE ({0:d}): {1:8.5f} [Hz]".format(i+1, freq_3D))
    
    # Initialize function and assign eigenvector
    eigenmode = Function(V,name="Eigenvector "+str(i))
    eigenmode.vector()[:] = rx
    file_results.write(eigenmode.split()[0], 0)

How much memory do you have on your system?

Does the code run in serial?

Does the code run on a smaller problem?

Dear Sir,
I am running this code on a 192 GB RAM system. The same code is running correctly for coarser mesh of the same geometry (degrees of freedom=9189) in parallel. Also, the code runs correctly at DoF of 9189 in serial too.

With increasing fineness in the same geometry (at degrees of freedom=245613), the code runs in serial correctly
serial

But when I am running it through parallelization , it is giving the same error as giving in 1.1 million DoF.