Hello all. I am having a problem solving a large matrix problem using the PETSc CG solver. The matrix assembly craps out at a certain point. It seems that the PETSc solver can only handle matrices below a certain size, since I can reduce the element order and the solution works. The error I get is:
<Mesh of topological dimension 3 (tetrahedra) with 7350 vertices and 37319 cells, ordered>
Solution vector size = 2868756
Traceback (most recent call last):
File "MonopoleLSFEM.py", line 122, in <module>
A, b = assemble_system(a, L, bcs)
File "/usr/lib/python3/dist-packages/dolfin/fem/assembling.py", line 382, in assemble_system
assembler.assemble(A_tensor, b_tensor)
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to successfully call PETSc function 'MatXIJSetPreallocation'.
*** Reason: PETSc error code is: 1 ((null)).
*** Where: This error was encountered inside /build/dolfin-GfMsdO/dolfin-2019.1.0/dolfin/la/PETScMatrix.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------My computer is a Dell laptop (Inspiron 15) with 16GB RAM and 4 processor cores running Ubuntu 18 with Fenics 2019 installed from PPA and working. The python code is
import meshio
msh = meshio.read("Monopole.msh")
meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells["tetra"]}))
from mshr import *
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
import os, sys, traceback
h = 0.15
semicirc_r = 2.0
a = 0.8
b = 0.8
ts = 0.1
lc = 1.0
rc = 0.11
cc = 0.0325
poffset = 0.2
pi = 3.1415926536
tol = 1.0e-12
eta = 377.0
k0 = 1.047
eps0 = 1.0 # WG air filled
eps_c = 2.2 # Teflon coax dielectric
eps_s = 1.00 # PCB substrate
class PEC(SubDomain):
def inside(self, x, on_boundary):
return on_boundary
class InputBC(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], -h-lc, tol)
class OutputBC(SubDomain):
def inside(self, x, on_boundary):
r_b = sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2])
r_c = sqrt(x[0] * x[0] + x[1] * x[1])
return on_boundary and ((near(r_b, 2.0, 5e-2) and x[2] > 0.0) or (near(r_c, 2.0, 5e-2) and between(x[2], (-h, 0))))
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 3)
info(mesh)
#plot(mesh)
#plt.show()
# Mark boundaries
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domains.set_all(3)
pec = PEC()
pec.mark(sub_domains, 0)
in_port = InputBC()
in_port.mark(sub_domains, 1)
out_port = OutputBC()
out_port.mark(sub_domains, 2)
File("BoxSubDomains.pvd").write(sub_domains)
Dk = Expression('x[2] <= 0.0 + tol ? (x[2] <= -h? eps_c : eps_s) : eps0', degree = 0, tol = 1.0e-12, h = h, eps_s = eps_s, eps_c = eps_c, eps0 = eps0) # Dielectric subdomains
# Set up function spaces
# For low order problem
cell = tetrahedron
ele_type = FiniteElement('N1curl', cell, 3) # H(curl) element for EM
V2 = FunctionSpace(mesh, MixedElement([ele_type, ele_type, ele_type, ele_type]))
V = FunctionSpace(mesh, ele_type)
e_r, e_i, h_r, h_i = TrialFunctions(V2)
ev_r, ev_i, hv_r, hv_i = TestFunctions(V2)
#surface integral definitions from boundaries
ds = Measure('ds', domain = mesh, subdomain_data = sub_domains)
# with source and sink terms
u0 = Constant((0.0, 0.0, 0.0)) #PEC definition
e_src = Expression(('x[0] / (2.0 * pi * (pow(x[0], 2.0) + pow(x[1] - poffset,2.0)))', '(x[1] - poffset) / (2.0 * pi *(pow(x[0],2.0) + pow(x[1] - poffset,2.0)))', '0.0'), degree = 2, poffset = poffset)
#Boundary condition dictionary
boundary_conditions = {0: {'PEC' : u0},
1: {'InputBC': (e_src, eps_c)},
2: {'OutputBC': 1.0}}
n = FacetNormal(mesh)
#Build PEC boundary conditions for real and imaginary parts
bcs = []
integral_source = []
integrals_load =[]
for i in boundary_conditions:
if 'PEC' in boundary_conditions[i]:
bb1 = (inner(cross(n, ev_r), cross(n, e_r)) + inner(cross(n, ev_i), cross(n, e_i)) + inner(dot(n, hv_r), dot(n, h_r)) + inner(dot(n, hv_i), dot(n, h_i))) * ds(i)
integrals_load.append(bb1)
# Build input BC source term and loading term
for i in boundary_conditions:
if 'InputBC' in boundary_conditions[i]:
Zm = 1.0/sqrt(Dk)
r, s = boundary_conditions[i]['InputBC']
bb = inner(cross(n, cross(ev_i, n)) - Zm * cross(n, hv_i), r) * ds(i)
integral_source.append(bb)
bb1 = (inner(cross(n, cross(ev_r, n)) - Zm * cross(n, hv_r), cross(n, cross(e_r, n)) - Zm * cross(n, h_r)) + inner(cross(n, cross(ev_i, n)) - Zm * cross(n, hv_i), cross(n, cross(e_i, n)) - Zm * cross(n, h_i))) * ds(i)
integrals_load.append(bb1)
for i in boundary_conditions:
if 'OutputBC' in boundary_conditions[i]:
Zm = 1.0/sqrt(Dk)
bb1 = (inner(cross(n, cross(ev_r, n)) - Zm * cross(n, hv_r), cross(n, cross(e_r, n)) - Zm * cross(n, h_r)) + inner(cross(n, cross(ev_i, n)) - Zm * cross(n, hv_i), cross(n, cross(e_i, n)) - Zm * cross(n, h_i))) * ds(i)
integrals_load.append(bb1)
a = (inner(curl(ev_r) - k0 * hv_i, curl(e_r) - k0 * h_i) + inner(curl(ev_i) + k0 * hv_r, curl(e_i) + k0 * h_r) + inner(curl(hv_r) + k0 * Dk * ev_i, curl(h_r) + k0 * Dk * e_i) + inner(curl(hv_i) - k0 * Dk * ev_r, curl(h_i) - k0 * Dk * e_r)) * dx + sum(integrals_load)
L = sum(integral_source)
u1 = Function(V2)
vdim = u1.vector().size()
print ('Solution vector size = ',vdim)
u1.vector()[:] = np.random.uniform(-1, 1, vdim)
A, b = assemble_system(a, L, bcs)
solver = PETScKrylovSolver("cg","sor")
solver.parameters['absolute_tolerance'] = 1e-5
solver.parameters['relative_tolerance'] = 1e-12
solver.parameters['maximum_iterations'] = 10000
solver.parameters['monitor_convergence'] = True
solver.parameters['nonzero_initial_guess'] = False
solver.parameters['error_on_nonconvergence'] = False
solver.parameters['report'] = True
solver.set_operators(A, A)
solver.solve(u1.vector(),b)
u1_r, u1_i, v1_r, v1_i = u1.split(True)
sys.exit(0)The GMSH mesh script is here:
// Monopole antenna model
h = 0.15; // substrate height
r = 2.0; // Ground plane radius
a = 1.35; // Monopole length
b = 0.8;
t = 0.1;
lc = 1.0; // length coax
rc = 0.11; // Coax shield radius
cc = 0.0325; // Coax center conductor radius
s = 0.2; // Coax probe offset
SetFactory("OpenCASCADE");
Sphere(1) = {0, 0, 0, r};
Box(2) = {-r, -r, 0, 2.0*r, 2.0*r, r};
BooleanIntersection(3) = {Volume{1}; Delete;}{Volume{2}; Delete;};
Cylinder(4) = {0, 0, -h, 0, 0, h, r}; // Substrate
Cylinder(5) = {0, s, -h-lc, 0, 0, lc, rc}; // Coax shield
Cylinder(6) = {0, s, -h-lc, 0, 0, h+lc+t/2, cc}; // coax center
Cylinder(7) = {0, s, 0, 0, 0, a, 2*cc}; // Antenna element
BooleanUnion(8) = {Volume{6}; Delete;}{Volume{7}; Delete;}; //to be subtracted out
BooleanDifference(9) = {Volume{3}; Delete;}{Volume{8};};
BooleanDifference(10) = {Volume{4}; Delete;}{Volume{8};};
BooleanDifference(11) = {Volume{5}; Delete;}{Volume{8}; Delete;};
BooleanFragments{Volume{9, 10, 11}; Delete;}{}
Physical Volume("Patch") = {9, 10, 11};
Mesh.CharacteristicLengthMax = 0.55;
Mesh.MinimumCirclePoints = 12;
Mesh.CharacteristicLengthFromCurvature = 1;
Mesh.Algorithm3D = 4;Is it possible to solve large FE problems using the PETSc solver? The memory in my computer was not from exhausted when the solver died. (Process monitor indicated 14GB memory used when process failed. Swap was not yet heavily used.)
Tasks: 253 total, 1 running, 198 sleeping, 0 stopped, 0 zombie
%Cpu(s): 2,2 us, 0,6 sy, 0,0 ni, 97,2 id, 0,0 wa, 0,0 hi, 0,0 si, 0,0 st
KiB Mem : 16311448 total, 13953404 free, 1460220 used, 897824 buff/cache
KiB Swap: 4986876 total, 4070556 free, 916320 used. 14313228 avail Mem Any help is appreciated.