I'm trying to solve an elasticity equation on a 3D domain, but even with 94GB of RAM I run out of memory when trying to use a mesh that is 140x140x140. This seems completely unreasonable. Below is my code. I'm running it using mpirun with 8 processors. If anyone has any ideas of where all of the memory is going and how I can make things more efficent, please let me know.
from dolfin import *
import numpy, copy, time
parameters['form_compiler']['optimize'] = True
parameters['form_compiler']['cpp_optimize'] = True
# define periodic boundary
class PeriodicBoundary(SubDomain):
def inside(self, x, on_boundary):
# return True if on left or bottom boundary AND NOT on one of the two slave edges
return bool((near(x[0], 0.) or near(x[1], 0.) or near(x[2], 0.)) and
(not ((near(x[1], 1.) and near(x[2], 0.)) or
(near(x[1], 0.) and near(x[2], 1.)) or
(near(x[0], 0.) and near(x[2], 1.)) or
(near(x[0], 1.) and near(x[2], 0.)) or
(near(x[0], 0.) and near(x[1], 1.)) or
(near(x[0], 1.) and near(x[1], 0.)))) and on_boundary)
def map(self, x, y):
if near(x[0], 1.) and near(x[1], 1.) and near(x[2], 1.):
y[0] = x[0] - 1.
y[1] = x[1] - 1.
y[2] = x[2] - 1.
elif near(x[1], 1.) and near(x[2], 1.):
y[0] = x[0]
y[1] = x[1] - 1.
y[2] = x[2] - 1.
elif near(x[0], 1.) and near(x[1], 1.):
y[0] = x[0] - 1.
y[1] = x[1] - 1.
y[2] = x[2]
elif near(x[0], 1.) and near(x[2], 1.):
y[0] = x[0] - 1.
y[1] = x[1]
y[2] = x[2] - 1.
elif near(x[0], 1.):
y[0] = x[0] - 1.
y[1] = x[1]
y[2] = x[2]
elif near(x[1], 1.):
y[0] = x[0]
y[1] = x[1] - 1.
y[2] = x[2]
elif near(x[2], 1):
y[0] = x[0]
y[1] = x[1]
y[2] = x[2] - 1.
else:
y[0] = -1000.
y[1] = -1000.
y[2] = -1000.
# Create mesh
nn = 140
mesh = UnitCubeMesh(nn, nn, nn)
V = VectorFunctionSpace(mesh, "CG", 1, constrained_domain=PeriodicBoundary())
class Fiber(SubDomain):
def inside(self, x, on_boundary):
return True if ((x[0]-0.5)**2+(x[1]-0.5)**2 <= 0.25**2) else False
fib = Fiber()
domains = CellFunction("size_t", mesh)
domains.set_all(0)
fib.mark(domains, 1)
dy = Measure("dx")[domains]
# define tensors
bbEm = numpy.zeros(81)
bbEm = bbEm.reshape(3,3,3,3)
bbEf = numpy.zeros(81)
bbEf = bbEf.reshape(3,3,3,3)
delta = [[1,0,0],[0,1,0],[0,0,1]]
# Elastic parameters
Em, Ef, nu = 10.5, 150000.0, 0.3
mum, lmbdam = Em/(2*(1 + nu)), Em*nu/((1 + nu)*(1 - 2*nu))
muf, lmbdaf = Ef/(2*(1 + nu)), Ef*nu/((1 + nu)*(1 - 2*nu))
# define bbEm
for ii in range(3):
for jj in range(3):
for kk in range(3):
for ll in range(3):
bbEm[ii,jj,kk,ll] = mum*(delta[ii][kk]*delta[jj][ll]+delta[ii][ll]*delta[kk][jj])+lmbdam*delta[ii][jj]*delta[kk][ll]
# define bbEf
for ii in range(3):
for jj in range(3):
for kk in range(3):
for ll in range(3):
bbEf[ii,jj,kk,ll] = muf*(delta[ii][kk]*delta[jj][ll]+delta[ii][ll]*delta[kk][jj])+lmbdaf*delta[ii][jj]*delta[kk][ll]
bbEmC = Constant(bbEm)
bbEfC = Constant(bbEf)
# define basis for V
e0 = Constant([1,0,0])
e1 = Constant([0,1,0])
e2 = Constant([0,0,1])
e = [e0, e1, e2]
# Define functions
u = TrialFunction(V)
v = TestFunction(V)
# Set up linear problem
a = bbEmC[i,j,k,l]*grad(u[k])[l]*grad(v[i])[j]*dy(0) \
+ bbEfC[i,j,k,l]*grad(u[k])[l]*grad(v[i])[j]*dy(1)
A = assemble(a)
z = Function(V)
b = None
# solve the problem
ii, jj = 0, 0
L = - bbEmC[i,j,k,l]*outer(e[ii],e[jj])[k,l]*grad(v[i])[j]*dy(0) \
- bbEfC[i,j,k,l]*outer(e[ii],e[jj])[k,l]*grad(v[i])[j]*dy(1)
b = assemble(L, tensor=b)
# Create Krylov solver
solver = KrylovSolver(A, "gmres", "amg")
# Create vector that spans the null space
null_vec = Vector(z.vector())
V.dofmap().set(null_vec, 1.0)
null_vec *= 1.0/null_vec.norm("l2")
# Create null space basis object and attach to Krylov solver
null_space = VectorSpaceBasis([null_vec])
solver.set_nullspace(null_space)
# Orthogonalize b with respect to the null space (this gurantees that
# a solution exists)
null_space.orthogonalize(b);
solver.solve(z.vector(), b)
