Firedrake is an automated system for the portable solution of partial differential equations using the finite element method (FEM)
Home Page: https://firedrakeproject.org
License: Other
At the moment, to build an extruded mesh the user must pass in a C kernel describing how the z coordinate of the coordinate field should be built.
For complicated cases, I guess there's no other way round this, however for simple cases (e.g. X equi-spaced layers for a total height of Y) there should be utility methods of doing so.
Related to this, users should probably have to provide an op2.Kernel, rather than just the C code (that way, they can name it whatever they like).
The following fails when run on 3 MPI process:
from firedrake import *
m = UnitSquareMesh(2,2)with an assertion:
Source location: (Python_Interface_f.F90, 364)
Error message: Failed assertion 2*ifacets == facets - efacets
On the process where it fails, the mesh looks like:
+--+--+
| /| /|
|/ |/ |
+--+--+
| /| /
|/ |/
+--+
i.e. almost but not quite the complete serial mesh.
The extruded version of this broke, then I noticed the non-extruded version broke too.
from firedrake import *
mesh = UnitSquareMesh(2, 2)
V1 = FunctionSpace(mesh, "CG", 1)
V2 = VectorFunctionSpace(mesh, "CG", 1)
f = Function(V1)
g = Function(V2)
f.interpolate(Expression("x[0]"))
g.interpolate(Expression(("-x[1]", "x[0]")))
W = V1 * V2
u, p = TrialFunctions(W)
v, q = TestFunctions(W)
a = u*v*dx + dot(p,q)*dx
L = f*v*dx + dot(g,q)*dx
out = Function(W)
solve(a == L, out, solver_parameters={'pc_type': 'jacobi', 'ksp_type': 'cg'}) # BOOM
print out.dat.data
Eventually needed for dynamical core project (but not too soon as it will annoy Jemma!).
"""Tests for mixed Helmholtz convergence on extruded meshes"""
import numpy as np
from firedrake import *
for ii in [1, 2, 3, 4]:
m = UnitSquareMesh(2**ii, 2**ii)
mesh = ExtrudedMesh(m, layers=(2**ii + 1), layer_height = 1.0/2**ii)
horiz_elt = FiniteElement("RT", "triangle", 1)
vert_elt = FiniteElement("DG", "interval", 0)
product_elt = HDiv(OuterProductElement(horiz_elt, vert_elt))
V1 = FunctionSpace(mesh, product_elt)
alt_horiz_elt = FiniteElement("DG", "triangle", 0)
product_elt = OuterProductElement(alt_horiz_elt, vert_elt)
V2 = FunctionSpace(mesh, product_elt)
f = Function(V2)
exact = Function(V2)
f.interpolate(Expression("(1+8*pi*pi)*sin(x[0]*pi*2)*sin(x[1]*pi*2)"))
exact.interpolate(Expression("sin(x[0]*pi*2)*sin(x[1]*pi*2)"))
W = V1 * V2
u, p = TrialFunctions(W)
v, q = TestFunctions(W)
a = (p*q - q*div(u) + dot(v, u) + div(v)*p)*dx
L = f*q*dx
out = Function(W)
solve(a == L, out)
print "L_inf norm: " + str(max(abs(out.dat[1].data - exact.dat.data)))
print "L2 norm: " + str(sqrt(assemble((out[3]-exact)*(out[3]-exact)*dx)))
L_inf norm: 0.131777047356
L2 norm: 0.131777047355
L_inf norm: 0.146867971137
L2 norm: 0.0597095688761
L_inf norm: 0.0389877820493
L2 norm: 0.015987880734
Traceback (most recent call last):
File "helmholtz.py", line 32, in <module>
solve(a == L, out)
File "/home/atm112/local/firedrake/python/firedrake/solving.py", line 650, in solve
_solve_varproblem(*args, **kwargs)
File "/home/atm112/local/firedrake/python/firedrake/solving.py", line 678, in _solve_varproblem
solver.solve()
File "/home/atm112/local/firedrake/python/firedrake/solving.py", line 251, in solve
(self.snes.getIterationNumber(), reason))
RuntimeError: Nonlinear solve failed to converge after 1 nonlinear iterations with reason: unknown reason (petsc4py enum incomplete?)
from firedrake import *
m = UnitTriangleMesh()
V = FunctionSpace(m, 'CG', 1)
u = Function(V)
assemble(sqrt(u)*dx)
=>
RuntimeError: In instant.recompile: The module did not compile with command 'python setup.py build_ext install --install-platlib=.', see '/data/lmitche1/instant-error/instant_module_8be95ee025f7b4457796d4db0755698ebe3eb8c9/compile.log'If I look there I see:
...
for (int ip = 0; ip<6; ip++)
{
double F0 = 0.0;
for (int r = 0; r<3; r++)
{
F0 += (w0[r][0]*FE0[ip][r]);
}
A[0] += (None*det*W6[ip]);
}Note the None in the increment into A[0] which should probably be sqrt.
from firedrake import *
m = UnitSquareMesh(8, 8)
mesh = ExtrudedMesh(m, 3, layer_height=0.5)
V = VectorFunctionSpace(mesh, "CG", 1)
expr = Expression(("x[0]", "x[1]", "x[2]"))
project (expr, V) # BOOM
Longer version:
from firedrake import *
m = UnitSquareMesh(8, 8)
mesh = ExtrudedMesh(m, 3, layer_height=0.5)
V = VectorFunctionSpace(mesh, "CG", 1)
expr = Expression(("x[0]", "x[1]", "x[2]"))
#project (expr, V)
fs = VectorFunctionSpace(mesh, "CG", V.ufl_element().degree()+1, dim=3)
v = Function(fs)
v.interpolate(expr)
p = TestFunction(V)
q = TrialFunction(V)
a = inner(p, q) * dx
L = inner(p, v) * dx
ret = Function(V)
#solve(a == L, ret, bcs=None, solver_parameters={}, form_compiler_parameters=None) # also BOOM
problem = LinearVariationalProblem(a, L, ret, None, form_compiler_parameters=None)
solver = LinearVariationalSolver(problem, parameters={})
solver.solve() # BOOM
Dies inside python/firedrake/solving.py at the following call:
with self._problem.u_ufl.dat.vec as v:
self.snes.solve(None, v)
self.snes is a PETSc thing, which is way outside my knowledge. @wence- ?
Request ability to assemble interior facet integrals for extruded meshes. Required in dynamical core project for temperature advection, pressure gradient term (and possibly vorticity equation).
Mixed Helmholtz attempt.
Using branches hdivcurl of FIAT/UFL/FFC, extrusion-vfs of PyOP2, extrusion-hdivcurl-proper of firedrake
from firedrake import *
m = UnitSquareMesh(4, 4)
mesh = ExtrudedMesh(m, 5, layer_height=0.25)
horiz_elt = FiniteElement("RT", "triangle", 1)
vert_elt = FiniteElement("DG", "interval", 0)
product_elt = HDiv(OuterProductElement(horiz_elt, vert_elt))
V1 = FunctionSpace(mesh, product_elt)
horiz_elt = FiniteElement("DG", "triangle", 0)
vert_elt = FiniteElement("DG", "interval", 0)
product_elt = OuterProductElement(horiz_elt, vert_elt)
V2 = FunctionSpace(mesh, product_elt)
f = Function(V2)
exact = Function(V2)
f.interpolate(Expression("(1+8*pi*pi)*sin(x[0]*pi*2)*sin(x[1]*pi*2)"))
exact.interpolate(Expression("sin(x[0]*pi*2)*sin(x[1]*pi*2)"))
W = V1 * V2
u, p = TrialFunctions(W)
v, q = TestFunctions(W)
a = (p*q - q*div(u) + dot(v, u) + div(v)*p)*dx
L = f*q*dx
out = Function(W)
solve(a == L, out)
print out.dat.data
# print sqrt(assemble((out[2]-exact)*(out[2]-exact)*dx))
Gives
atm112@ubuntu:~$ python mixedbork.py
In instant.recompile: The module did not compile with command 'python setup.py build_ext install --install-platlib=.', see '/home/atm112/.instant/error/instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876/compile.log'
Traceback (most recent call last):
File "mixedbork.py", line 28, in <module>
solve(a == L, out)
File "/home/atm112/local/firedrake/python/firedrake/solving.py", line 644, in solve
_solve_varproblem(*args, **kwargs)
File "/home/atm112/local/firedrake/python/firedrake/solving.py", line 672, in _solve_varproblem
solver.solve()
File "/home/atm112/local/firedrake/python/firedrake/solving.py", line 229, in solve
self.snes.solve(None, v)
File "SNES.pyx", line 404, in petsc4py.PETSc.SNES.solve (src/petsc4py.PETSc.c:124245)
File "petscsnes.pxi", line 223, in petsc4py.PETSc.SNES_Function (src/petsc4py.PETSc.c:27067)
File "/home/atm112/local/firedrake/python/firedrake/solving.py", line 156, in form_function
with self._F_tensor.dat.vec_ro as v:
File "/usr/lib/python2.7/contextlib.py", line 17, in __enter__
return self.gen.next()
File "/home/atm112/local/PyOP2/pyop2/petsc_base.py", line 149, in vecscatter
with acc(d) as v:
File "/usr/lib/python2.7/contextlib.py", line 17, in __enter__
return self.gen.next()
File "/home/atm112/local/PyOP2/pyop2/petsc_base.py", line 84, in vec_context
self._force_evaluation()
File "/home/atm112/local/PyOP2/pyop2/base.py", line 1335, in _force_evaluation
_trace.evaluate(set([self]), set([self]))
File "/home/atm112/local/PyOP2/pyop2/base.py", line 145, in evaluate
comp._run()
File "/home/atm112/local/PyOP2/pyop2/base.py", line 2881, in _run
return self.compute()
File "/home/atm112/local/PyOP2/pyop2/base.py", line 2888, in compute
self._compute_if_not_empty(self.it_space.iterset.core_part)
File "/home/atm112/local/PyOP2/pyop2/base.py", line 2900, in _compute_if_not_empty
self._compute(part)
File "/home/atm112/local/PyOP2/pyop2/sequential.py", line 101, in _compute
fun(*self._jit_args)
File "/home/atm112/local/PyOP2/pyop2/host.py", line 411, in __call__
return self.compile()(*args)
File "/home/atm112/local/PyOP2/pyop2/host.py", line 452, in compile
modulename=self._kernel.name if configuration["debug"] else None)
File "/usr/local/lib/python2.7/dist-packages/instant/inlining.py", line 95, in inline_with_numpy
module = build_module(**kwargs)
File "/usr/local/lib/python2.7/dist-packages/instant/build.py", line 542, in build_module
recompile(modulename, module_path, new_compilation_checksum, build_system)
File "/usr/local/lib/python2.7/dist-packages/instant/build.py", line 123, in recompile
instant_error(msg % (cmd, compile_log_filename_dest))
File "/usr/local/lib/python2.7/dist-packages/instant/output.py", line 57, in instant_error
raise RuntimeError(text)
RuntimeError: In instant.recompile: The module did not compile with command 'python setup.py build_ext install --install-platlib=.', see '/home/atm112/.instant/error/instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876/compile.log'
Content of file:
running build_ext
building '_instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876' extension
creating build
creating build/temp.linux-x86_64-2.7
mpicc -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wstrict-prototypes -fPIC -I/usr/local/lib/python2.7/dist-packages/petsc/include -I/usr/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c mat_utils.cxx -o build/temp.linux-x86_64-2.7/mat_utils.o
cc1plus: warning: command line option ‘-Wstrict-prototypes’ is valid for C/ObjC but not for C++ [enabled by default]
mpicc -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wstrict-prototypes -fPIC -I/usr/local/lib/python2.7/dist-packages/petsc/include -I/usr/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876_wrap.cxx -o build/temp.linux-x86_64-2.7/instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876_wrap.o
cc1plus: warning: command line option ‘-Wstrict-prototypes’ is valid for C/ObjC but not for C++ [enabled by default]
In file included from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarraytypes.h:1728:0,
from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarrayobject.h:17,
from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/arrayobject.h:15,
from instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876_wrap.cxx:3197:
/usr/lib/python2.7/dist-packages/numpy/core/include/numpy/npy_deprecated_api.h:11:2: warning: #warning "Using deprecated NumPy API, disable it by #defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
#warning "Using deprecated NumPy API, disable it by #defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION"
^
instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876_wrap.cxx: In function ‘void wrap_form_cell_integral_0_otherwise__(PyObject*, PyObject*, PyObject*, PyObject*, PyObject*, PyObject*, PyObject*, PyObject*, PyObject*, PyObject*, PyObject*, PyObject*, PyObject*, PyObject*, PyObject*, PyObject*, PyObject*, PyObject*, PyObject*, PyObject*)’:
instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876_wrap.cxx:3254:11: error: ‘xtr_arg0_0_map1_0’ was not declared in this scope
xtr_arg0_0_map1_0[0] = *(arg0_0_map1_0 + i * 1 + 0);;
^
instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876_wrap.cxx:3257:74: warning: left operand of comma operator has no effect [-Wunused-value]
form_cell_integral_0_otherwise(arg0_0 + xtr_arg0_0_map0_0[i_0]*(1, 1), arg1_0_vec, arg2_0_vec, arg3_0_vec, i_0 + 0);
^
instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876_wrap.cxx:3280:40: warning: left operand of comma operator has no effect [-Wunused-value]
arg3_0_vec[0] += _off30[0] * (1, 1);
^
instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876_wrap.cxx:3281:40: warning: left operand of comma operator has no effect [-Wunused-value]
arg3_0_vec[1] += _off30[1] * (1, 1);
^
instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876_wrap.cxx:3282:40: warning: left operand of comma operator has no effect [-Wunused-value]
arg3_0_vec[2] += _off30[2] * (1, 1);
^
instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876_wrap.cxx:3283:40: warning: left operand of comma operator has no effect [-Wunused-value]
arg3_0_vec[3] += _off30[3] * (1, 1);
^
instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876_wrap.cxx:3291:47: error: ‘xtr_arg0_0_map0_1’ was not declared in this scope
form_cell_integral_0_otherwise(arg0_0 + xtr_arg0_0_map0_1[i_0]*(1, 1), arg1_0_vec, arg2_0_vec, arg3_0_vec, i_0 + 3);
^
instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876_wrap.cxx:3291:74: warning: left operand of comma operator has no effect [-Wunused-value]
form_cell_integral_0_otherwise(arg0_0 + xtr_arg0_0_map0_1[i_0]*(1, 1), arg1_0_vec, arg2_0_vec, arg3_0_vec, i_0 + 3);
^
instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876_wrap.cxx:3314:40: warning: left operand of comma operator has no effect [-Wunused-value]
arg3_0_vec[0] += _off30[0] * (1, 1);
^
instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876_wrap.cxx:3315:40: warning: left operand of comma operator has no effect [-Wunused-value]
arg3_0_vec[1] += _off30[1] * (1, 1);
^
instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876_wrap.cxx:3316:40: warning: left operand of comma operator has no effect [-Wunused-value]
arg3_0_vec[2] += _off30[2] * (1, 1);
^
instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876_wrap.cxx:3317:40: warning: left operand of comma operator has no effect [-Wunused-value]
arg3_0_vec[3] += _off30[3] * (1, 1);
^
instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876_wrap.cxx:3205:11: warning: unused variable ‘arg0_1’ [-Wunused-variable]
double *arg0_1 = (double *)(((PyArrayObject *)_arg0_1)->data);
^
instant_module_c6462c135a902b9c4679ed9d63d5d267d43b9876_wrap.cxx:3224:9: warning: unused variable ‘_off31’ [-Wunused-variable]
int * _off31 = (int *)(((PyArrayObject *)off31)->data);
^
error: command 'mpicc' failed with exit status 1
For some reason, the following two notations are not equivalent:
solve(A == b, x)
and
solve(assemble(A), x.vector(), assemble(b))
For the example below, the difference in the norm is about 1e-09.
from firedrake import *
from numpy.linalg import norm
mesh = UnitSquareMesh(50, 50)
V = FunctionSpace(mesh, "CG", 1)
f = Function(V)
f.vector()[:] = 1.
t = TestFunction(V)
q = TrialFunction(V)
a = inner(t, q)*dx
L = inner(f, t)*dx
# Solve the system using forms
sol = Function(V)
solve(a == L, sol)
# And again
sol2 = Function(V)
solve(a == L, sol2)
assert norm(sol.vector()[:] - sol2.vector()[:]) == 0 # Passes
# Solve the system using preassembled objects
sol3 = Function(V)
solve(assemble(a), sol3.vector(), assemble(L).dat)
assert norm(sol.vector()[:] - sol3.vector()[:]) == 0 # Fails
The README is completely out of date and should be replace by a PyOP2 style README in rST format with installation instructions etc.
m = UnitSquareMesh(4, 4)
mesh = ExtrudedMesh(m, 11, layer_height=0.1)
V2_a_horiz = FiniteElement("RT", "triangle", 2)
V2_a_vert = FiniteElement("DG", "interval", 1)
V2_a = HDiv(OuterProductElement(V2_a_horiz, V2_a_vert))
V2_b_horiz = FiniteElement("DG", "triangle", 1)
V2_b_vert = FiniteElement("CG", "interval", 2)
V2_b = HDiv(OuterProductElement(V2_b_horiz, V2_b_vert))
V2_elt = EnrichedElement(V2_a, V2_b)
V2 = FunctionSpace(mesh, V2_elt)
V3 = FunctionSpace(mesh, "DG", 1, vfamily="DG", vdegree=1)
f = Function(V3)
exact = Function(V3)
f.interpolate(Expression("(3*pi*pi)*sin(x[0]*pi)*sin(x[1]*pi)*sin(x[2]*pi)"))
exact.interpolate(Expression("sin(x[0]*pi)*sin(x[1]*pi)*sin(x[2]*pi)"))
W = V2 * V3
u, p = TrialFunctions(W)
v, q = TestFunctions(W)
a = (-q*div(u) + dot(v, u) + div(v)*p)*dx
L = f*q*dx
bcs = [DirichletBC(W[0], Expression(("0.0", "0.0", "pi*sin(x[0]*pi)*sin(x[1]*pi)")), "bottom"),
DirichletBC(W[0], Expression(("0.0", "0.0", "-pi*sin(x[0]*pi)*sin(x[1]*pi)")), "top")]
out = Function(W)
solve(a == L, out, bcs=bcs, solver_parameters={'pc_type': 'fieldsplit',
'pc_fieldsplit_type': 'schur',
'ksp_type': 'cg',
'pc_fieldsplit_schur_fact_type': 'FULL',
'fieldsplit_0_ksp_type': 'cg',
'fieldsplit_1_ksp_type': 'cg'})
u, p = out.split()
print "L_inf norm: " + str(max(abs(p.dat.data - exact.dat.data)))
print "L2 norm: " + str(sqrt(assemble((p-exact)*(p-exact)*dx)))
plotting = project(p, FunctionSpace(mesh, "CG", 1))
file = File("mixed_bcs_poisson.pvd")
file << plotting
Annoyingly this only seems to happen on extruded meshes, which points the finger in our direction.
Anyway, here's what happens: we solve two equations separately and simultaneously. When solved simultaneously, the matrix looks correct, but the answer is different. "WTF?" @wence- @kynan
from firedrake import *
m = UnitSquareMesh(3, 3)
mesh = ExtrudedMesh(m, 11, layer_height=0.1)
V1 = FunctionSpace(mesh, "CG", 1)
V2 = FunctionSpace(mesh, "CG", 2)
f = Function(V1)
g = Function(V2)
f.interpolate(Expression("x[0]"))
g.interpolate(Expression("x[0]"))
d1 = TrialFunction(V1)
d2 = TestFunction(V1)
d3 = assemble(d1*d2*dx)
d4 = assemble(f*d2*dx)
e1 = TrialFunction(V2)
e2 = TestFunction(V2)
e3 = assemble(e1*e2*dx)
e4 = assemble(g*e2*dx)
W = V1 * V2
u, p = TrialFunctions(W)
v, q = TestFunctions(W)
a = u*v*dx + p*q*dx
L = f*v*dx + g*q*dx
out = Function(W)
A = assemble(a)
b = assemble(L)
# Our matrix system looks correct - we compare the blocks and the RHS with the separate components
print max(abs(A.M[0,0].array - d3.M.array))
print max(abs(A.M[1,1].array - e3.M.array))
print max(abs(A.M[0,1].array))
print max(abs(A.M[1,0].array))
print max(abs(b.dat.data[0] - d4.dat.data))
print max(abs(b.dat.data[1] - e4.dat.data))
d5 = Function(V1)
solve(d3,d5,d4)
e5 = Function(V2)
solve(e3,e5,e4)
# if we solve the problems separately, we get an error of 10^-7, the Krylov solver tolerance
print [max(abs(d5.dat.data - f.dat.data)), max(abs(e5.dat.data - g.dat.data))]
# yet if we solve them together....
solve(a == L, out, solver_parameters={'pc_type': 'jacobi', 'ksp_type': 'cg'})
# the error is O(1)
print [max(abs(out.dat.data[0] - f.dat.data)), max(abs(out.dat.data[1] - g.dat.data))]
Required to re-implement our shallow water on the sphere stuff (and first step to extrusion on sphere).
Hi,
this simple script results in a *** FLUIDITY ERROR ***:
from firedrake import *
UnitIntervalMesh(2)
The error message I get is:
*** FLUIDITY ERROR ***
Source location: (Fields_Allocates.F90, 1467)
Error message: Surface element in sndgln is not on the surface.
Backtrace will follow if it is available:
/home/firedrake/src/firedrake/lib/libfluidity.so.1(fprint_backtrace_+0x1f) [0x7fa07dec9a59]
/home/firedrake/src/firedrake/lib/libfluidity.so.1(__fldebug_MOD_flabort_pinpoint+0x277) [0x7fa07dee6ef7]
/home/firedrake/src/firedrake/lib/libfluidity.so.1(__fields_allocates_MOD_add_faces_face_list+0xd55) [0x7fa07df53c15]
/home/firedrake/src/firedrake/lib/libfluidity.so.1(__fields_allocates_MOD_add_faces+0x2056) [0x7fa07df5afd6]
/home/firedrake/src/firedrake/lib/libfluidity.so.1(__read_triangle_MOD_read_triangle_files_to_field+0x1f33) [0x7fa07e02dc43]
/home/firedrake/src/firedrake/lib/libfluidity.so.1(__read_triangle_MOD_read_triangle_simple+0x29b) [0x7fa07e030afb]
/home/firedrake/src/firedrake/lib/libfluidity.so.1(__mesh_files_MOD_read_mesh_simple+0x49e) [0x7fa07dfde81e]
/home/firedrake/src/firedrake/lib/libfluidity.so.1(read_mesh_f+0x447) [0x7fa07e003077]
/home/firedrake/src/firedrake/python/firedrake/core_types.so(+0x24cc5) [0x7fa07e5d3cc5]
python() [0x4ce25e]
python(PyObject_Call+0x36) [0x4abf96]
/home/firedrake/src/firedrake/python/firedrake/core_types.so(+0x3005d) [0x7fa07e5df05d]
python(PyEval_EvalFrameEx+0xa5b) [0x47bf2b]
python(PyEval_EvalCodeEx+0x19a) [0x4e048a]
python() [0x4e0eea]
python() [0x4ce25e]
python() [0x4a0534]
python() [0x4ae24e]
python(PyEval_EvalFrameEx+0xa5b) [0x47bf2b]
python(PyEval_EvalCodeEx+0x19a) [0x4e048a]
python(PyEval_EvalCode+0x32) [0x53ff12]
python() [0x54038b]
python(PyRun_FileExFlags+0x92) [0x4658f4]
python(PyRun_SimpleFileExFlags+0x2d8) [0x465e24]
python(Py_Main+0xb4e) [0x466b99]
/lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf5) [0x7fa07ef14ea5]
python() [0x4e1115]
Use addr2line -e <binary> <address> to decipher.
Error is terminal.
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 16.
NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
I installed the entire PyOP2 + Firedrake system today, so I should be on the most up to date version.
from firedrake import *
m = UnitTriangleMesh()
mesh = ExtrudedMesh(m, layers=3, layer_height=1.0)
V0 = FunctionSpace(mesh, "CG", 1, vfamily="CG", vdegree=1)
V1_a_horiz = FiniteElement("RT", "triangle", 1)
V1_a_vert = FiniteElement("CG", "interval", 1)
V1_a = HCurl(OuterProductElement(V1_a_horiz, V1_a_vert))
V1_b_horiz = FiniteElement("CG", "triangle", 2)
V1_b_vert = FiniteElement("DG", "interval", 0)
V1_b = HCurl(OuterProductElement(V1_b_horiz, V1_b_vert))
V1_elt = EnrichedElement(V1_a, V1_b)
V1 = FunctionSpace(mesh, V1_elt)
u = Function(V0)
u.interpolate(Expression("x[2]"))
v = TrialFunction(V1)
w = TestFunction(V1)
a = dot(v, w)*dx
L = dot(grad(u), w)*dx
v = Function(V1)
parms = {'snes_type': 'ksponly', 'ksp_type': 'preonly', 'pc_type': 'lu'}
solve(a == L, v, solver_parameters=parms)
print np.sort(v.dat.data) # expect 12 1s and 9 zeros
I'm trying to do something similar to the DOLFIN example here so I can set up Python initial conditions and boundary conditions in Firedrake:
from firedrake import *
import libspud
class ScalarExpressionFromOptions(Expression):
def __init__(self, path):
self.path = path
if(libspud.have_option(path + "/constant")):
self.constant = True
self.source_value = libspud.get_option(path + "/constant")
elif(libspud.have_option(path + "/python")):
self.constant = False
self.source_value = libspud.get_option(path + "/python")
else:
print "No value specified"
sys.exit(1)
def eval(self, values, x):
if(self.constant):
values[0] = self.source_value
else:
values[0] = self.get_values_from_python(x)
def get_values_from_python(self, x):
s = self.source_value
exec s
return val(x)
mesh = UnitSquareMesh(10, 10)
FS = FunctionSpace(mesh, "CG", 1)
h_old = Function(FS)
# Load options tree
libspud.load_options("test.swml")
# Set up initial conditions
h_initial = ScalarExpressionFromOptions(path = "/material_phase[0]/scalar_field::FreeSurfacePerturbationHeight/prognostic/initial_condition")
h_old.interpolate(h_initial)
File("test.pvd") << h_old
After the instantiation of a new ScalarExpressionFromOptions object called h_initial, I get the following error message when trying to do h_old.interpolate(h_initial) where h_old is a Function:
christian@elevate ~ $ python expression_subclassing.py
Info : Increasing process stack size (8192 kB < 16 MB)
Traceback (most recent call last):
File "expression_subclassing.py", line 36, in <module>
h_old.interpolate(h_initial)
File "core_types.pyx", line 1532, in firedrake.core_types.Function.interpolate (firedrake/core_types.c:27404)
AttributeError: 'ScalarExpressionFromOptions' object has no attribute 'code'
Within the firedrake docs, cross-linking between firedrake classes that are defined in different files typically doesn't produce an actual link in the sphinx docs.
I think the write way to fix this is to write:
:class:`.FiredrakeClass`instead of
:class:`FiredrakeClass`(Note leading dot in the former case)
Sphinx will then try and resolve the class in the namespace of the current module.
I think this is because FIAT doesn't know about them. In ffc they're built out of Lagrange elements of appropriate degree and then turned into RestrictedElements by removing the edge dofs.
This, to me at least, seems like a horrible conflation of abstractions, but maybe Bubble functions don't fit into FIAT.
In DOLFIN, when one calls Function(f), the result is a copy of f, including the coefficient values.
Currently in Firedrake, the result is a zero Function.
Since the policy of Firedrake is to be DOLFIN compatible, I suggest that this behaviour should be changed.
This is a test script that demonstrates the problem:
from firedrake import *
mesh = UnitIntervalMesh(2)
V = FunctionSpace(mesh, "CG", 1)
f = Function(V)
f.interpolate(Expression("1"))
assert (f.dat.vec.array == [1, 1, 1]).all() # Passes
g = Function(f)
print g.dat.vec.array
assert (g.dat.vec.array == [1, 1, 1]).all() # Fails
I will shortly make a merge request with a fix.
If I run firedrake tests individually, I pass all of them. However, there are a couple of tests (test_projection and test_helmholtz_mixed) that fail when I execute the entire test suite.
Probably this is something to do with caching, as the error I got is a compile-time error due to the buggy code that FFC used to generate (which is now fixed, anyway). So, apparently, firedrake is somewhat picking up old kernel code
typo: "with your exisiting Pyhton environment" -> "with your existing Python environment"
install instruction typo: pip install "Cython=>0.17" decorator "numpy>=1.6" pyyaml
=> should be >=
instructions for FFC are outdated: we use master rather than the pyop2 branch. this affects the following lines:
git clone -b pyop2 https://bitbucket.org/mapdes/ffc.git $FFC_DIR
git+https://bitbucket.org/mapdes/ffc.git@pyop2#egg=ffc
instructions for installing FIAT without pip are missing.
The P1 mass matrix tests fail.
This is because the tolerance is too tight:
assert error[0] < 1.0e-15
actually:
firedrake_identity: Assigning error = [1.2212453270876722e-15, 3.8552338933622061e-07, 3.0702377464247326e-07, 9.549121590369981e-07]
The tolerance changed in the extrusion merge from 1e-14 to 1e-15, why?
The core_types extension module is not compatible with Python 2.6 due to the use of capsules, which were introduced in 2.7.
A possible workaround would be introducing a compatibility header that redefines PyCapsule C API calls in terms of PyCObject.
We don't want to have to set the PETSC_OPTIONS environment variable to pass options to the SNES.
I think the right way to do this is to take an options string (or perhaps dict mapping petsc options to values) then build a unique prefix for this snes solver, set the prefix on the options and the solver and then do setfromoptions.
The following extrusion tests don't test anything:
The problem is that firedrake.adjoint is not overloaded and ufl.adjoint returns a form with new ufl.Arguments. Latter have lost the additional information that firedrake stores in firedrake.Arguments, and a result the output of ufl.adjoint can not be assembled.
Following test script demonstrates this issue:
from firedrake import *
mesh = UnitSquareMesh(2, 2)
V = FunctionSpace(mesh, "CG", 1)
test = TestFunction(V)
trial = TrialFunction(V)
f = Function(V)
form = f*test*trial*dx
print map(type, form.compute_form_data().original_arguments)
print map(type, firedrake.adjoint(form).compute_form_data().original_arguments)
assemble(firedrake.adjoint(form))
Implemented as built-in mesh. Required for 3D dynamical core development.
Request: kernels that assemble local matrices (including mixed function spaces), and then apply direct solver (BLAS) to local matrix/RHS system, and write output to chosen field.
This would be great for reconstructing the mass flux from the density advection to feed to the other advection schemes in the dynamical core project, but I have a possible workaround that is less efficient but will get it working.
The extruded MMS tests take way too long to run:
test_firedrake_extrusion_helmholtztest_firedrake_extrusion_projectFiredrake's implementation of names of Function objects is inconsistent and incompatible with Dolfin:
Firedrake:
f = Function(V)
f.name = "My Function"
str(f ) # "My Function"
Dolfin:
f = Function(V)
f.rename("My Function", "My label")
f.name() # "My Function"
str(f ) # "My Function"
f.label() # "My label"
We're currently only testing the sequential backend and chances are not all the Firedrake functionality works on all the backends. Hence this needs testing.
Prior to that we might want to make the backend selection a bit easier in the UI.
Verify as follows:
mpirun -np 2 python -c 'from firedrake import *; UnitIntervalMesh(20)'
Firedrake Expressions are somewhat defunct at the moment: They derived from ufl.Coefficient but they really are not coefficients i.e. ufl.Coefficient.__init__ is never called and therefore a firedrake.Expression cannot be used as a ufl.Coefficient.
The reason is that information required to initialise a Coefficient (i.e. the element) is not available at the time the Expression is constructed.
Needed for hybridised solvers (to implement timestepping for dynamical core work).
We need Hdiv and Hcurl fields that still apply Piola transformations but are discontinuous. David says this can be done by fiddling with the local DOF diagram (to make all DOFs associated with the element interior) before passing to global numbering, but not before passing to FFC.
Firedrake currently does not implement the Function.vector() functionality that exists in DOLFIN.
An example is:
f = Function(V)
fvec = f.vector()
fvec[2] = 1
print fvec.array()
This feature is currently developed in the vector branch of firedrake..
To populate the initial tracer field this test does the following:
f = firedrake.Function(fs)
populate_p0 = op2.Kernel("""
void populate_tracer(double *x[])
{
x[0][0] = 3;
}""", "populate_tracer")
op2.par_loop(populate_p0, f.cell_set,
f.dat(op2.INC, f.cell_node_map))Let me count the sins:
f.assign(3) ?The following code calculates the integral of z^2 in an extruded unit square in two different ways.
The second one gives the wrong answer
@dham can say more
from firedrake import *
m = UnitSquareMesh(1, 1)
mesh = ExtrudedMesh(m, 3, layer_height=0.5)
expr1 = Expression("x[2]")
exactfs1 = FunctionSpace(mesh, "CG", 1)
exact1 = Function(exactfs1)
exact1.interpolate(expr1)
print assemble(exact1*exact1*dx) # 1/3, correct
expr2 = Expression(("0.0", "0.0", "x[2]"))
exactfs2 = VectorFunctionSpace(mesh, "CG", 1)
exact2 = Function(exactfs2)
exact2.interpolate(expr2)
print assemble(dot(exact2, exact2)*dx) # should be 1/3, but isn't
The tests of pyop2 initialisation are order-dependent. This works fine when running in serial, the first test to run is the one that checks that we haven't initialised pyop2 yet. However, running in parallel, we can't guarantee that this test runs first on the slave process, (in which case it fails).
The following works with a UnitSquareMesh:
from firedrake import *
mesh = UnitSquareMesh(10, 10)
fs = FunctionSpace(mesh, "CG", 1)
h = Function(fs).interpolate(Expression("x[0]"))
File("test.pvd") << h
but when I try with mesh = UnitIntervalMesh(10), Firedrake gives the following error message when trying to write h to a file:
christian@elevate ~ $ python test_file.py
Info : Increasing process stack size (8192 kB < 16 MB)
Traceback (most recent call last):
File "test_file.py", line 8, in <module>
File("test.pvd") << h
File "/home/christian/firedrake/python/firedrake/io.py", line 73, in __lshift__
self._file << data
File "/home/christian/firedrake/python/firedrake/io.py", line 375, in __lshift__
self._update_PVD(data)
File "/home/christian/firedrake/python/firedrake/io.py", line 393, in _update_PVD
new_vtk << function
File "/home/christian/firedrake/python/firedrake/io.py", line 160, in __lshift__
coordinates = self._fd_to_evtk_coord(mesh._coordinates)
File "/home/christian/firedrake/python/firedrake/io.py", line 258, in _fd_to_evtk_coord
if len(fdcoord[0]) == 3:
TypeError: object of type 'numpy.float64' has no len()
Required to do dynamical core on spherical annulus (i.e., the atmosphere).
Apparently this issue has come up before. When using UnitSquareMesh(10, 10) in Firedrake, gmsh successfully creates a .msh file from a .geo file in /tmp:
christian@elevate ~ $ python firedrake_gmsh_test.py
Info : Running 'gmsh /tmp/firedrake-mesh-cache-uid1000/unitsquare_10_10.geo -2 -o /tmp/firedrake-mesh-cache-uid1000/unitsquare_10_10.msh' [1 node(s), max. 4 thread(s)]
Info : Started on Fri Dec 13 15:27:18 2013
Info : Reading '/tmp/firedrake-mesh-cache-uid1000/unitsquare_10_10.geo'...
Info : Done reading '/tmp/firedrake-mesh-cache-uid1000/unitsquare_10_10.geo'
Info : Meshing 1D...
Info : Meshing curve 1 (extruded)
Info : Meshing curve 2 (extruded)
Info : Meshing curve 3 (extruded)
Info : Meshing curve 4 (extruded)
Info : Done meshing 1D (0 s)
Info : Meshing 2D...
Info : Meshing surface 5 (extruded)
Info : Done meshing 2D (0 s)
Info : 121 vertices 244 elements
Info : Writing '/tmp/firedrake-mesh-cache-uid1000/unitsquare_10_10.msh'...
Info : Done writing '/tmp/firedrake-mesh-cache-uid1000/unitsquare_10_10.msh'
Info : Stopped on Fri Dec 13 15:27:18 2013
but Firedrake then seems to stall. I have to kill the Python process (which reports full CPU usage):
PID USER PR NI VIRT RES SHR S %CPU %MEM TIME+ COMMAND
christia 2149 20 0 407m 38m 15m R 101 0.5 4:56.14 python
and then re-run Firedrake to get things working properly.
I've been looking at the Poisson demo recently added by Florian. How come the following works
solve(a == L, w, bcs=[bc0, bc1], solver_parameters={'ksp_type':'gmres', 'pc_type':'jacobi'})
but pre-assembling the system and then solving
A = assemble(a)
b = assemble(L)
solve(A, w, b, bcs=[bc0, bc1], solver_parameters={'ksp_type':'gmres', 'pc_type':'jacobi'})
fails with the following error message?
christian@elevate ~ $ python test_mixed.py
Info : Increasing process stack size (8192 kB < 16 MB)
Traceback (most recent call last):
File "test_mixed.py", line 26, in <module>
solve(A, w, b, bcs=[bc0, bc1], solver_parameters={'ksp_type':'gmres', 'pc_type':'jacobi'}) # Fails
File "/home/christian/firedrake/python/firedrake/solving.py", line 654, in solve
return _la_solve(*args, **_kwargs)
File "/home/christian/firedrake/python/firedrake/solving.py", line 558, in _la_solve
solver.solve(A.M, x.dat, b.dat)
File "/home/christian/PyOP2/PyOP2/pyop2/base.py", line 3244, in solve
self._solve(A, x, b)
File "/home/christian/PyOP2/PyOP2/pyop2/petsc_base.py", line 455, in _solve
raise RuntimeError(msg)
RuntimeError: KSP Solver failed to converge in 1000 iterations: DIVERGED_MAX_IT (Residual norm: 1.417050e-03)
This inconsistency is something I have also seen when running my shallow water model.
The full code is:
from firedrake import *
mesh = UnitSquareMesh(32, 32)
BDM = FunctionSpace(mesh, "BDM", 1)
DG = FunctionSpace(mesh, "DG", 0)
W = BDM * DG
sigma, u = TrialFunctions(W)
tau, v = TestFunctions(W)
f = Function(DG).interpolate(Expression(
"10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)"))
a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
L = - f*v*dx
bc0 = DirichletBC(W.sub(0), Expression(("0.0", "-sin(5*x[0])")), 1)
bc1 = DirichletBC(W.sub(0), Expression(("0.0", "sin(5*x[0])")), 2)
w = Function(W)
A = assemble(a)
b = assemble(L)
#solve(a == L, w, bcs=[bc0, bc1], solver_parameters={'ksp_type':'gmres', 'pc_type':'jacobi'})
solve(A, w, b, bcs=[bc0, bc1], solver_parameters={'ksp_type':'gmres', 'pc_type':'jacobi'}) # Fails
sigma, u = w.split()
File("poisson_mixed.pvd") << u
Is it worth/necessary to add a thin wrapper around PyOP2 types in Firedrake to allow adding additional functionality? In particular for OP2/PyOP2#240 we would probably want to be able to multiply a PyOP2 Mat with a Firedrake Vector (again for DOLFIN compatibility). Thoughts?
Installation of Firedrake + dependencies + running our test suite on empty cache now takes > 50min and therefore times out on Travis.
There's a few potential fixes:
Leads to errors like
ImportError: No module named pyop2, yaml, ...
One solution would be to copy FEniCS's Dorsal by writing out a
add the following line to your .bashrc: source /blah/blah.conf
or
export PYTHONPATH=...
from firedrake import *
m = ExtrudedMesh(UnitTriangleMesh(), layers=2, layer_height=1)
V = FunctionSpace(m, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
f = Function(V)
L = f*v*dx
out = Function(V)
out.dat.data[0] = 1
assembled = np.round(assemble(ufl.action(a, out)).dat.data - assemble(L).dat.data - assemble(ufl.action(a, out) - L).dat.data, decimals=3)
print assembled
assert np.allclose(assembled, 0.0)There is not proper test of the VTK I/O functionality. Only the test_firedrake_extrusion_helmholtz writes a VTK file but doesn't do anything with it.
As a follow-up from #10: Many of the Mesh classes are not tested, in particular UnitIntervalMesh, UnitCircleMesh, UnitTetrahedronMesh, UnitTriangleMesh.
I'm trying to implement a DG version of my shallow water model in Firedrake. The solver always fails after just 1 iteration whenever facet integrals are included, with the following error message:
Residual norms for firedrake_snes_1_ solve.
0 KSP Residual norm -nan
Traceback (most recent call last):
File "test_dS.py", line 24, in <module>
solve(a == L, solution, bc, solver_parameters={'ksp_monitor':True})
File "/home/christian/firedrake/python/firedrake/solving.py", line 645, in solve
_solve_varproblem(*args, **kwargs)
File "/home/christian/firedrake/python/firedrake/solving.py", line 673, in _solve_varproblem
solver.solve()
File "/home/christian/firedrake/python/firedrake/solving.py", line 244, in solve
nonlinear iterations with reason: %s" % (self.snes.getIterationNumber(), reason))
RuntimeError: Nonlinear solve failed to converge after 1 nonlinear iterations with reason: unknown reason (petsc4py enum incomplete?)
When the assembled RHS is printed out, several NaN values are present:
[ 0.05000009 -0.0270154 0.04833402 -0.06815219 0.04449678 0.04382895
0.04272321 0.04119058 0.0392464 0.03691007 0.03420496 0.03115808
0.0836053 nan nan nan nan nan
nan nan nan nan 0.14473835 0.14497394
0.14521259 0.14545125 0.1456899 0.14592855 0.14616721 0.14639499
nan -0.01573746 nan nan nan nan
nan nan nan nan 0.14970218 0.1759264
0.17728686 0.17864733 0.18000779 0.18136826 0.18272872 0.18407939
nan 0.1730736 0.2007407 0.20391732 0.20709394 0.21027056
0.21344718 0.2166238 0.21979181 nan 0.16972407 0.19872415
0.20385172 0.20897929 0.21410685 0.21923442 0.22436199 0.22948222
nan 0.16467871 0.19472201 0.2017493 0.20877658 0.21580386
0.22283115 0.22985843 0.23687972 nan 0.15798794 0.18877427
0.19763106 0.20648785 0.21534463 0.22420142 0.23305821 0.24191039
nan 0.14971861 0.18094037 0.19153816 0.20213596 0.21273375
0.22333155 0.23392934 0.24452398 nan 0.13995333 0.17129856
0.18353148 0.19576439 0.2079973 0.22023022 0.23246313 0.24469437
nan nan nan nan nan nan
nan nan nan nan nan nan
nan nan nan nan nan nan
nan]
I've attached a simplified example which reproduces the issue below - here the code tries to solve for the divergence of a vector field. Can anyone see something wrong with this code?
from firedrake import *
mesh = UnitSquareMesh(10, 10)
VFS = VectorFunctionSpace(mesh, "DG", 1)
SFS = FunctionSpace(mesh, "CG", 1)
n = FacetNormal(mesh.ufl_cell())
w = TestFunction(SFS)
divu = TrialFunction(SFS)
velocity = project(Expression(("cos(x[0])", "x[1]*sin(x[0])")), VFS)
# Integration of w*div(velocity) by parts
M = inner(w, divu)*dx
Ct_facet = inner(jump(w, n), avg(velocity))*dS
Ct_int = -inner(grad(w), velocity)*dx
F = M - Ct_int - Ct_facet
a = lhs(F); L = rhs(F)
print assemble(L).vector().array() # Several NaN values appear in the Vector here when the dS integral is included.
solution = Function(SFS)
bc = DirichletBC(SFS, Expression("0.0"), (1,2,3,4))
solve(a == L, solution, bc, solver_parameters={'ksp_monitor':True})
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
Django
The Web framework for perfectionists with deadlines.
Laravel
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
Microsoft
Open source projects and samples from Microsoft.
Google
Google ❤️ Open Source for everyone.
Alibaba
Alibaba Open Source for everyone
D3
Data-Driven Documents codes.
Tencent
China tencent open source team.