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 . 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.
The MWE is
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-
= BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
= PID 98 RUNNING AT 70a39b6ee647
= EXIT CODE: 11
The MWE for this case is also attached here-
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.