# Large system solver issue

**URL:** <https://fenicsproject.discourse.group/t/large-system-solver-issue/7676>\
**Category:** Linear Algebra\
**Created:** [February 12, 2022, 7:32am UTC](https://fenicsproject.discourse.group/t/large-system-solver-issue/7676 "2022-02-12T07:32:58Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![Shubham\_Saurabh](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/shubham_saurabh/32/3460_2.png) [@Shubham\_Saurabh](https://fenicsproject.discourse.group/u/Shubham_Saurabh)\
**Post date:** [February 12, 2022, 7:32am UTC](https://fenicsproject.discourse.group/t/large-system-solver-issue/7676/1 "2022-02-12T07:32:59Z")

</div>

I am solving a eigen value problem of 3D box with 3 point mass of some other geometries distributed in the internal cylindrical faces. The geometry and mesh files are here- [geometry and mesh](https://drive.google.com/drive/folders/1SLXbInv-xVl8VYADQbNLz5CBskD5rVmy?usp=sharing). The code is working well for coarser mesh of DoF of 2.5 lakhs. But when I am making the geometry finer of DoF 1.1 millions, the same code is giving an error- `unable to extract eigen pairs from SLEPc eigen value solver`. Sometimes kernel gets dead too.

 ![Screenshot 2022-01-08 232635](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/5/5190912bad4f7063076f64ff85ba4dd63b8cf107.png)  
The MWE is

```auto
from dolfin import *
import numpy as np
import json

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")

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)

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
print("Computing {} first eigenvalues...".format(N_eig))
eigensolver.solve(N_eig)

# Set up file for exporting results
file_results = XDMFFile("modal_analysis.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True

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
    print("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)

```

When I am running it through **parallelization** , it is giving an error-

```auto
= BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
= PID 98 RUNNING AT 70a39b6ee647
= EXIT CODE: 11

```

 ![mpi_error](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/c/ca2401abd2fd74fb795f27463f72fa92bb645ced.png)

The MWE for this case is also attached here-

```auto

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)

```

Can you please help in rectifying both the modes of solving this problem, with or without parallelization.