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)
