### Mixed Function Spaces in Fenics 2.0

765
views
1
8 months ago by
I am trying to solve a modified elasticity problem. The code was written back in 2016, and since then, "MixedFunctionSpaces" was deprecated. I followed the solution in:

https://fenicsproject.org/qa/11983/mixedfunctionspace-in-2016-2-0

V = VectorFunctionSpace(mesh, 'CG', 1)
F = FunctionSpace(mesh, 'CG', 1)
MS = MixedFunctionSpace([V, V, F, F])​

I used:

p = 1
V_ele = VectorElement("CG", mesh.ufl_cell(), p)
F_ele = FiniteElement("CG", mesh.ufl_cell(), p)
MS = FunctionSpace(mesh, MixedElement([V_ele, V_ele, F_ele, F_ele]))

, but the system kills the process a few minutes after showing "Solving linear variation problem.". There are no errors other than the message "Killed", but I get no results either. This code used to run fine and fast on an older version of fenics.

Since I am using more than one function in my mixed function space, I cannot use the "*" operator.

I would appreciate any help about what is best to use instead of MixedFunctionSpaces. If indeed the above solution is the best way, what is causing the system to kill the run at the solution step?

Community: FEniCS Project
1
Here is the full code body:

from dolfin import *

# Material properties:
mu = 1.0				# shear modulus
rom2 = 36*pi*pi				# rho * omega^2 = 36pi^2 gives wavelength = 1/3.
ks = sqrt(rom2/mu)

# Create functions for boundary conditions
top = 0.5
moden = 1.0
sn = moden*pi/(ks*top)
cn = sqrt(1.0-sn*sn)
zero = Constant(0)

# Exact Solution:
uexactr = Expression(("(-1.0*s/c*cos(m*pi*x[1]/height)*sin(c*k*x[0]))", "(sin(m*pi*x[1]/height)*cos(c*k*x[0]))", "0.0"),
m = moden, s = sn, c = cn, k = ks, height = top, degree=1)

uexacti = Expression( ( "(s/c*cos(m*pi*x[1]/height)*cos(c*k*x[0]))", "(sin(m*pi*x[1]/height)*sin(c*k*x[0]))", "0.0"),
m = moden, s = sn, c = cn, k = ks, height = top, degree=1)

# Geometry and boundaries:
top = 0.50
length = 2.0
width  = 0.5
resolution = 8			# Number of elements along x axis per (real or total) wavelength = 2pi/ks.

# Mesh and mesh parameters:
# Rectangle geometry
nel = int(round(resolution*length*ks/(2*pi)))	# Number of elements along x axis.
h = length/nel					# Element size.
nelx = nel					# Number of elements along x axis.
nely = int(round(top/h))			# Number of elements along y axis.

mesh  = BoxMesh(Point(0, 0, 0), Point(length, top, width), nelx, nely, nely)

# Right boundary for boundary condition:
class Right(SubDomain):
def inside(self, x, on_boundary):
return x[0] > length - DOLFIN_EPS

right = Right()

# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
right.mark(boundaries, 1)

ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

# Computational Parameters:
c = cos(ks*h)
alpha = ks*ks*h*h/6.0
tau1 = h*h/mu/2.0
tau2 = ((1+tau1*rom2)*alpha*(2.0+c)/(1.0-c) - 1.0)/(rom2*ks*ks)

# ========================================================================================================
# Define function spaces

# Old version (before fenics 2.0) - used to work:
#Vp = FunctionSpace(mesh, "CG", 1)
#Vu = VectorFunctionSpace(mesh, "CG", 1)
#V = MixedFunctionSpace([Vu, Vu, Vp, Vp])

# New version (suggested in the forums) - does not work apparently:
p = 1
V_ele = VectorElement("CG", mesh.ufl_cell(), p)
F_ele = FiniteElement("CG", mesh.ufl_cell(), p)
V = FunctionSpace(mesh, MixedElement([V_ele, V_ele, F_ele, F_ele]))

# Define Dirichlet boundaries: top, bottom, right, left:
def lboundary(x):
return (x[0] < DOLFIN_EPS)

def rboundary(x):
return x[0] > length - DOLFIN_EPS

def tbboundary(x):
return x[1] < DOLFIN_EPS or x[1] > top - DOLFIN_EPS

def frontbackboundary(x):
return x[2] > width - DOLFIN_EPS or x[2] < DOLFIN_EPS

def pboundary(x):
return bool(x[1] < DOLFIN_EPS or x[1] > top - DOLFIN_EPS or x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS)

# Enforce Dirichlet condition on y=x1 component of the displacement on top/bottom boundaries:
bctbr = DirichletBC(V.sub(0).sub(1), zero, tbboundary)     # x1 or y component
bctbi = DirichletBC(V.sub(1).sub(1), zero, tbboundary)     # x1 or y component

# Enforce exact solution as Dirichlet condition on left boundary:
bclr = DirichletBC(V.sub(0), uexactr, lboundary)
bcli = DirichletBC(V.sub(1), uexacti, lboundary)

# Enforce zero normal displacement as Dirichlet condition on front+back boundaries:
bcfbr = DirichletBC(V.sub(0).sub(2), zero, frontbackboundary)
bcfbi = DirichletBC(V.sub(1).sub(2), zero, frontbackboundary)

# Collect boundary conditions
# bcs = [bctbr, bctbi, bclr, bcli]
bcs = [bctbr, bctbi, bclr, bcli, bcfbr, bcfbi]

# Define variational problem
(ur,ui,pr,pi) = TrialFunctions(V)
(vr,vi,qr,qi) =  TestFunctions(V)
f = Constant((0, 0, 0))

ar = ( \
- 2.0/3.0*mu*div(vr)*div(ur)					\
- inner(vr, rom2*ur) - div(vr)*pr 					\
- qr*div(ur) 							\
+ ks*mu*(2.0*cn*vr[0]*ui[0] + (cn*cn-sn*sn)/cn*vr[1]*ui[1])*ds(1)

ai = ( \
- 2.0/3.0*mu*div(vi)*div(ui)					\
- inner(vi, rom2*ui) - div(vi)*pi 					\
- qi*div(ui) 							\
- ks*mu*(2.0*cn*vi[0]*ur[0] + (cn*cn-sn*sn)/cn*vi[1]*ur[1])*ds(1)

a = ar + ai

L = inner(vr + tau1*grad(qr), f)*dx

# Compute solution
w = Function(V)
solve(a == L, w, bcs)​
written 8 months ago by Navid Nazari
1
Code finishes fine on my machine. It is possible that it gets killed by OOM killer because it eats your free memory. The code took up to 1.8 GB on my machine. (Tested with FEniCS 2017.2.0+.)
written 8 months ago by Jan Blechta
1
You can try also
solve(a == L, w, bcs, solver_parameters={"linear_solver": "mumps"})​
to use solver with more modest memory requirements than default UMFPACK.
written 8 months ago by Jan Blechta
Thank you very much for your quick response! It worked like a charm. Would you please share how you logged the memory usage? Also, if you reply it under the post, that would be great. Thanks!
written 8 months ago by Navid Nazari
https://en.wikipedia.org/wiki/Htop
written 8 months ago by Jan Blechta