Hi,
We have installed Dolfin, and for a simple poisson problem the code
works, however, for the attached example, and run in parallel, we get
the error below. Below the error is the cmake output which shows which
modules are installed.
Is there any module that should be installed additionally to make this work?
Thanks,
Bart
<<<<ERROR MESSAGE>>>>
Traceback (most recent call last):
File "code_steady.py", line 125, in <module>
Traceback (most recent call last):
File "code_steady.py", line 125, in <module>
pi = Function(DG)
File
"/projects/uoa00295/dolfin-1.5.0/install/dolfin/functions/function.py",
line 274, in __init__
pi = Function(DG)
File
"/projects/uoa00295/dolfin-1.5.0/install/dolfin/functions/function.py",
line 274, in __init__
self.__init_from_function_space(V)
File
"/projects/uoa00295/dolfin-1.5.0/install/dolfin/functions/function.py",
line 319, in __init_from_function_space
self.__init_from_function_space(V)
File
"/projects/uoa00295/dolfin-1.5.0/install/dolfin/functions/function.py",
line 319, in __init_from_function_space
cpp.Function.__init__(self, V)
File
"/projects/uoa00295/dolfin-1.5.0/install/dolfin/cpp/function.py", line
626, in __init__
cpp.Function.__init__(self, V)
File
"/projects/uoa00295/dolfin-1.5.0/install/dolfin/cpp/function.py", line
626, in __init__
_function.Function_swiginit(self,_function.new_Function(*args))
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
***
*** [email protected]
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
***
-------------------------------------------------------------------------
*** Error: Unable to creating uBLASVector.
*** Reason: Distributed uBLASVector is not supported.
*** Where: This error was encountered inside uBLASVector.cpp.
*** Process: unknown
***
*** DOLFIN version: 1.5.0
*** Git changeset:
***
-------------------------------------------------------------------------
******************************************************************************************
******************************************************************************************
******************************************************************************************
******************************************************************************************
******************************************************************************************
<<<<MODULES>>>>>
-- The following optional packages were found:
-- -------------------------------------------
-- (OK) OPENMP
-- (OK) MPI
-- (OK) PETSC
-- (OK) UMFPACK
-- (OK) CHOLMOD
-- (OK) SCOTCH
-- (OK) PARMETIS
-- (OK) ZLIB
-- (OK) PYTHON
-- (OK) VTK
--
-- The following optional packages were not enabled:
-- -------------------------------------------------
-- (--) PETSC4PY
-- (--) SLEPC
-- (--) SLEPC4PY
-- (--) TRILINOS
-- (--) PASTIX
-- (--) SPHINX
-- (--) HDF5
-- (--) QT
--
Dr. Bart Verleye
Centre for eResearch
Level G, Room 409-G21
24 Symonds St.
Auckland 1010
New Zealand
+64 (0) 9 923 9740 ext 89740
"""This demo program solves the mixed formulation of a transient nonlinear Poisson's
equation:
sigma - grad(u) = 0
DpDt - div(u*sigma) = f
The corresponding weak (variational problem)
<sigma, tau> + <div(tau), u> = 0 for all tau
<DpDt, v> - <u*div(sigma), v> - <sigma, sigma*v> = <f, v> for all v
is solved using BDM (Brezzi-Douglas-Marini) elements of degree k for
(sigma, tau) and DG (discontinuous Galerkin) elements of degree k - 1
for (u, v).
"""
from dolfin import *
import matplotlib.pyplot as plt
# Constants
P_FIN = 1.0e1
P_INI = 1.0e0
perm = 1.0e0
visc = 1.0e0
E = 10.0
nu = 0.3
mu = E / (2.0*(1.0 + nu))
lmbda = E*nu / ((1.0 + nu)*(1.0 - 2.0*nu))
# Time variables
dt = 1.0e-2; t = dt; T = 3.0
# Create mesh
mesh = UnitCubeMesh(5, 5, 5)
# Boundary info
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0)
class Front(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)
class Rear(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 1.0)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[2], 0.0)
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[2], 1.0)
# class Well(SubDomain):
# def inside(self, x, on_boundary):
# return (between(x[1], (0.5, 0.8)) and between(x[0], (0.2, 0.5)))
class Well(SubDomain):
def inside(self, x, on_boundary):
return True if pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2) <= 0.04 else False
left = Left()
top = Top()
right = Right()
bottom = Bottom()
front = Front()
rear = Rear()
well = Well()
# Initialize mesh function for interior domains
domains = CellFunction("size_t", mesh)
domains.set_all(0)
well.mark(domains, 1)
# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
front.mark(boundaries, 5)
rear.mark(boundaries, 6)
# plot(boundaries, interactive = True, title = 'boundaries')
# plot(domains, interactive = True, title = 'domains')
# Define function spaces and mixed (product) space
BDM = FunctionSpace(mesh, "BDM", 2)
DG = FunctionSpace(mesh, "DG", 1)
CG= VectorFunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace([BDM, DG, CG])
#
#--------------------------------------------
class Initial_p(Expression):
def eval(self, values, x):
# values[0] = pow(99.0*x[0] + 1.0, 0.5)
values[0] = P_INI
class Initial_u(Expression):
def eval(self, values, x):
# values[0] = 49.5/pow(99.0*x[0] + 1.0, 0.5)
values[0] = 0.0
values[1] = 0.0
values[2] = 0.0
def value_shape(self):
return (3,)
class Initial_us(Expression):
def eval(self, values, x):
values[0] = 0.0
values[1] = 0.0
values[2] = 0.0
def value_shape(self):
return (3,)
pi = Function(DG)
pi.interpolate(Initial_p())
ui = Function(BDM)
ui.interpolate(Initial_u())
usi = Function(CG)
usi.interpolate(Initial_us())
w = Function(W)
w0= Function(W)
# assign initial condition
assign( w.sub(2),usi)
assign( w.sub(1),pi)
assign( w.sub(0),ui)
assign(w0.sub(2),usi)
assign(w0.sub(1),pi)
assign(w0.sub(0),ui)
# ------------------------------------------- end of initilisation
# Define trial and test functions
# w = Function(W)
# (sigma, u) = split(w)
(u, p, us) = split(w)
(u0, p0, us0) = split(w0)
# (tau, v) = TestFunctions(W)
(v, q, vs) = TestFunctions(W)
# # Plot the initial values
# plot(u, axes=True, title = "Initial Velocity")
# plot(p, axes=True, title = "Initial Pressure")
# plot(us, axes=True, mode = 'displacement', title = "Initial Displacement")
# interactive()
# Define source function
# f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
# f = Constant(100000.0)
ff = Constant(0.0)
fs = Constant((0.0, 0.0, 0.0))
class p_top_expression(Expression):
def __init__(self, t):
self.t = t
def eval(self, values, x):
if self.t <=1.0:
values[0] = P_INI + (P_FIN - P_INI)*self.t
else:
values[0] = P_FIN
# p_top = p_top_expression(t)
p_top = Constant(P_FIN)
p_bottom = Constant(P_INI)
# External load on top face
class Traction(Expression):
def eval(self, values, x):
values[0] = 0.0
values[1] = 0.0
values[2] = 0.0
def value_shape(self):
return (3,)
force = Traction()
# Facet normal
n = FacetNormal(mesh)
# Define measures
dx = Measure("dx")[domains]
ds = Measure("ds")[boundaries]
# Define sigma
def sigma(us, p):
return 2.0*mu*sym(grad(us)) + lmbda*tr(sym(grad(us)))*Identity(us.geometric_dimension()) - p*Identity(us.geometric_dimension())
def sigma_e(us):
return 2.0*mu*sym(grad(us)) + lmbda*tr(sym(grad(us)))*Identity(us.geometric_dimension())
def mobility():
return perm/visc
# Define variational form
# linear
# F = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx + f*v*dx- u_left*dot(tau, n)*ds(1) - u_right*dot(tau, n)*ds(3)
# nonlinear
# # define nonlinear part
# def coe(u):
# # return Constant(1.0)
# return Constant(1.0)*u
# F = (dot(sigma, tau) + div(tau)*u + div(coe(u)*sigma)*v)*dx + f*v*dx- u_left*dot(tau, n)*ds(1) - u_right*dot(tau, n)*ds(3)
# F = (dot(sigma, tau) + div(tau)*u + u*nabla_div(sigma)*v + inner(sigma*v, sigma))*dx + f*v*dx- u_left*dot(tau, n)*ds(1) - u_right*dot(tau, n)*ds(3)
# Transient for compressible flow
# F = (dot(u, v) + div(v)*p + (p-p0)*q/dt - p*nabla_div(u)*q - inner(u*q, u))*dx(0) +(dot(u, v) + div(v)*p + (p-p0)*q/dt - p*nabla_div(u)*q - inner(u*q, u))*dx(1) - f*q*dx(1)- p_left*dot(v, n)*ds(1) - p_right*dot(v, n)*ds(3)
# Steady for coupled
F = (dot(u, v) + mobility()*div(v)*p - p*nabla_div(u)*q - (1.0/mobility())*inner(u*q, u) + inner(sigma(us, p), grad(vs)))*dx(0) + (dot(u, v) + mobility()*div(v)*p - p*nabla_div(u)*q - mobility()*inner(u*q, u) + inner(sigma(us, p), grad(vs)))*dx(1) - ff*q*dx(1)- mobility()*p_top*dot(v, n)*ds(2) - mobility()*p_bottom*dot(v, n)*ds(4) - inner(fs, vs)*dx(0) - inner(fs, vs)*dx(1)
# Be careful, volume integral for dx(0) and dx(1)
# u*nabla_div(sigma)*v*dx + inner(sigma, sigma)*v*dx
# Define function G such that G \cdot n = g
# This is to constrain the normal component of the flux vector
class BoundarySource(Expression):
def __init__(self, mesh):
self.mesh = mesh
def eval_cell(self, values, x, ufc_cell):
cell = Cell(self.mesh, ufc_cell.index)
n = cell.normal(ufc_cell.local_facet)
# g = sin(5*x[0])
g = 0.0
values[0] = g*n[0]
values[1] = g*n[1]
values[2] = g*n[2]
def value_shape(self):
return (3,)
G = BoundarySource(mesh)
class FixedFace(Expression):
def eval(self, values, x):
values[0] = 0.0
values[1] = 0.0
values[2] = 0.0
def value_shape(self):
return (3,)
fixed = FixedFace()
class Zero_disp(Expression):
def eval(self, values, x):
values[0] = .0
zero_d = Zero_disp()
# bc = DirichletBC(W.sub(0), G, boundary)
# Boundary condition for fluid
bc1 = [DirichletBC(W.sub(0), G, left), DirichletBC(W.sub(0), G, right), DirichletBC(W.sub(0), G, front), DirichletBC(W.sub(0), G, rear)]
# Boundary condition for solid
# bc.append(DirichletBC(W.sub(2).sub(2), Constant(0), bottom))
bc2 = [DirichletBC(W.sub(2), fixed, bottom), DirichletBC(W.sub(2).sub(0), Constant(0), left), DirichletBC(W.sub(2).sub(0), Constant(0), right), DirichletBC(W.sub(2).sub(1), Constant(0), front), DirichletBC(W.sub(2).sub(1), Constant(0), rear)]
bc = bc1 + bc2
# Compute solution
# solve(F == 0, w, bc)
dw = TrialFunction(W)
J = derivative(F, w, dw)
problem = NonlinearVariationalProblem(F, w, bc, J)
solver = NonlinearVariationalSolver(problem)
#solver.parameters["newton_solver"]["linear_solver"] = "umfpack"
#solver.parameters["newton_solver"]["maximum_iterations"] = 100
# solver.solve()
fig_p = plot(w.sub(1), axes=True, title = 'Pressure')
fig_us= plot(w.sub(2), axes=True, mode = 'displacement', title = 'Displacement')
# Create files for storing results
file_p = File("results/pressure.pvd")
file_p << (w.sub(1), 0.0)
file_us = File("results/displacement.pvd")
file_us << (w.sub(2), 0.0)
# file << (w.sub(1), 0.0)
# solve
solver.solve()
w0.vector()[:] = w.vector()
# (u, p) = w.split()
## If split is not used, p is intouchable to 'fig.plot()'' or 'file <<'
fig_p.plot()
fig_us.plot()
# Save to file
file_p << (w.sub(1), t)
file_us << (w.sub(2), t)
# Project and write stress field to post-processing file
ST = TensorFunctionSpace(mesh, "Discontinuous Lagrange", 0)
total_stress = project(sigma(w.sub(2), w.sub(1)), V=ST)
File("total_streess.pvd") << total_stress
effective_stress = project(sigma_e(w.sub(2)), V=ST)
File("effective_streess.pvd") << effective_stress
# Stress cannot be plotted directely because is is a tensor of order 3
'''
# While loop
while t <= T:
print 'Time : %.8f' %plot
p_left.t = t
solver.solve()
w0.vector()[:] = w.vector()
# (u, p) = w.split()
## If split is not used, p is intouchable to 'fig.plot()'' or 'file <<'
fig_p.plot()
fig_us.plot()
# Save to file
file_p << (w.sub(1), t)
file_us << (w.sub(2), t)
# print
for x in range(101):
print '%.3f %.4f ' %(x/100.0, p([x/100.0, 0.5, 0.5]))
# Update t
t += dt
'''
(u, p, us) = w.split()
# # print
# for x in range(101):
# print '%.3f %.4f ' %(x/100.0, p([0.5, 0.5, x/100.0]))
# # Plot sigma and u
# plot(u, axes=True, title = "Calculated Velocity")
# plot(p, axes=True, title = "Calculated Pressure")
# plot(us, axes=True, mode = 'displacement', title = "Calculated Displacement")
# interactive()
# plot in matplot
data_x = []
data_p = []
# for x in range(101):
# data_x.append(x/100.0)
# data_p.append(p([0.5, 0.5, x/100.0]))
# plt.figure()
# plt.plot(data_x, data_p)
# plt.show()
# Save solution
File('saved_solution.xml') << w
_______________________________________________
fenics-support mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics-support
