search

firedrakeproject / firedrake

Firedrake is an automated system for the portable solution of partial differential equations using the finite element method (FEM)
https://firedrakeproject.org
Other
496 stars 157 forks source link

Various manual-quadrature-degree weirdness #581

Closed dorugeber closed 8 years ago

dorugeber commented 9 years ago

Toy code (heavily stripped down 3D atmospheric stuff), running from a clean cache

import numpy as np
from firedrake import *

# standard things
parameters["assembly_cache"]["enabled"] = False
parameters["pyop2_options"]["lazy_evaluation"] = False
# parameters["pyop2_options"]["debug"] = True

# Mesh size
refinements = 2  # number of horizontal cells = 20*(4^refinements)
nlayers = 10  # number of vertical levels
dt = 4.0
T = 3600.0

newtonits = 5
alpha = Constant(0.5)        # off-centering; 1 = backward Euler, 0.5 = midpoint

# All test case constants
a_ref = 6.37122e6  # Radius of the Earth (m)
g = 9.80616  # gravity (m/s^2)
p_0 = 1000.0 * 100.0  # Reference pressure (Pa, not hPa)
c_p = 1004.5  # SHC of dry air at constant pressure (J/kg/K)
kappa = 2.0/7.0  # R_d/c_p

# Individual test case constants
z_top = 1.0e4  # Height position of the model top
X = 125.0  # Reduced-size Earth reduction factor
a = a_ref/X  # Scaled radius of the Earth

# Make mesh
m = IcosahedralSphereMesh(radius=a, refinement_level=refinements)

mesh = ExtrudedMesh(m, layers=nlayers, layer_height=z_top/nlayers, extrusion_type="radial")

S1 = FiniteElement("BDFM", triangle, 2)
S2 = FiniteElement("DG", triangle, 1)

T0 = FiniteElement("CG", interval, 2)
T1 = FiniteElement("DG", interval, 1)

W2h_elt = HDiv(OuterProductElement(S1, T1))
W2t_elt = OuterProductElement(S2, T0)
W2v_elt = HDiv(W2t_elt)
W2_elt = W2h_elt + W2v_elt
W3_elt = OuterProductElement(S2, T1)

W2h = FunctionSpace(mesh, W2h_elt)
# W2v = FunctionSpace(mesh, W2v_elt)
W2t = FunctionSpace(mesh, W2t_elt)
W2 = FunctionSpace(mesh, W2_elt)
W3 = FunctionSpace(mesh, W3_elt)
W23 = W2 * W3  # mixed u-p space

W_VectorDG0 = VectorFunctionSpace(mesh, "DG", 0)
W_VectorCG1 = VectorFunctionSpace(mesh, "CG", 1)

zhat = Function(W_VectorDG0)

# Initial/current conditions
u0, theta0, rho0 = Function(W2), Function(W2t), Function(W3)
# Increments
du, dtheta, drho = Function(W2), Function(W2t), Function(W3)

theta_b = Function(W2t)

rho_b = Function(W3)

theta_prime = Function(W2t)

theta0.assign(theta_b + theta_prime)
rho0.assign(rho_b)

# Set up equations
dtc = Constant(dt)

uimp = u0 + alpha*du; thetaimp = theta0 + alpha*dtheta; rhoimp = rho0 + alpha*drho
ustar = TrialFunction(W2)
w = TestFunction(W2)
n = FacetNormal(mesh)

# u_star solver setup

a_ustar = dot(w, ustar)*dx \
          + alpha*dtc*dot(ustar, curl(cross(uimp, w)))*dx \
          - 2*alpha*dtc*dot(avg(cross(n, cross(uimp, w))), avg((sign(dot(uimp, n))+1)*ustar))*(dS_h + dS_v)  # upwinding; "sign(...) + 1" is 0 or 2; avg makes this 0 or 1

exnerimp = (p_0/(kappa*c_p))**(kappa/(kappa-1)) * pow(rhoimp*thetaimp, kappa/(1-kappa))

L_ustar = dot(w, u0)*dx \
          - (1-alpha)*dtc*dot(u0, curl(cross(uimp, w)))*dx \
          + 2*(1-alpha)*dtc*dot(avg(cross(n, cross(uimp, w))), avg((sign(dot(uimp, n))+1)*u0))*(dS_h + dS_v) \
          + 0.5*dtc*div(w)*dot(uimp, uimp)*dx \
          + dtc*c_p*div(thetaimp*w)*exnerimp*dx \
          - dtc*c_p*jump(w*thetaimp, n=n)*avg(exnerimp)*dS_v \
          - dtc*dot(w, g*zhat)*dx

ustar_bcs = [DirichletBC(W2, 0.0, "bottom"),
             DirichletBC(W2, 0.0, "top")]

u_star = Function(W2)
ustar_prob = LinearVariationalProblem(a_ustar, L_ustar, u_star, bcs=ustar_bcs)
ustar_solv = LinearVariationalSolver(ustar_prob)

ustar_solv.solve()  # bork?

This works.

Cleaning cache and replacing all instances of dx with dx(degree=(4, 4)), I get the following error trace:

Traceback (most recent call last):
  File "horiz_nodeg.py", line 103, in <module>
    ustar_solv.solve()  # bork?
  File "<string>", line 2, in solve
  File "/home/atm112/local/PyOP2/pyop2/profiling.py", line 203, in wrapper
    return f(*args, **kwargs)
  File "/home/atm112/local/firedrake/firedrake/variational_solver.py", line 172, in solve
    self.snes.solve(None, v)
  File "PETSc/SNES.pyx", line 520, in petsc4py.PETSc.SNES.solve (src/petsc4py.PETSc.c:165186)
  File "PETSc/petscsnes.pxi", line 253, in petsc4py.PETSc.SNES_Function (src/petsc4py.PETSc.c:31721)
  File "/home/atm112/local/firedrake/firedrake/solving_utils.py", line 184, in form_function
    nest=problem._nest)
  File "/home/atm112/local/firedrake/firedrake/assemble.py", line 67, in assemble
    inverse=inverse, nest=nest)
  File "/home/atm112/local/firedrake/firedrake/assemble.py", line 449, in _assemble
    return thunk(bcs)
  File "/home/atm112/local/firedrake/firedrake/assembly_cache.py", line 324, in inner
    return thunk(bcs)
  File "/home/atm112/local/firedrake/firedrake/assemble.py", line 279, in thunk
    op2.par_loop(*args)
  File "<string>", line 2, in par_loop
  File "/home/atm112/local/PyOP2/pyop2/versioning.py", line 154, in modifies_arguments
    retval = func(*args, **kwargs)
  File "/home/atm112/local/PyOP2/pyop2/op2.py", line 278, in par_loop
    return backends._BackendSelector._backend.par_loop(kernel, iterset, *args, **kwargs)
  File "/home/atm112/local/PyOP2/pyop2/base.py", line 4351, in par_loop
    return _make_object('ParLoop', kernel, it_space, *args, **kwargs).enqueue()
  File "/home/atm112/local/PyOP2/pyop2/base.py", line 79, in enqueue
    _trace.append(self)
  File "/home/atm112/local/PyOP2/pyop2/base.py", line 96, in append
    computation._run()
  File "/home/atm112/local/PyOP2/pyop2/base.py", line 4014, in _run
    return self.compute()
  File "/home/atm112/local/PyOP2/pyop2/base.py", line 4038, in compute
    fun = self._jitmodule
  File "/home/atm112/local/PyOP2/pyop2/utils.py", line 64, in __get__
    obj.__dict__[self.__name__] = result = self.fget(obj)
  File "/home/atm112/local/PyOP2/pyop2/sequential.py", line 158, in _jitmodule
    direct=self.is_direct, iterate=self.iteration_region)
  File "/home/atm112/local/PyOP2/pyop2/caching.py", line 203, in __new__
    obj = make_obj()
  File "/home/atm112/local/PyOP2/pyop2/caching.py", line 193, in make_obj
    obj.__init__(*args, **kwargs)
  File "/home/atm112/local/PyOP2/pyop2/host.py", line 704, in __init__
    self.compile()
  File "/home/atm112/local/PyOP2/pyop2/host.py", line 802, in compile
    compiler=compiler.get('name'))
  File "/home/atm112/local/PyOP2/pyop2/compilation.py", line 269, in load
    dll = compiler.get_so(src, extension)
  File "/home/atm112/local/PyOP2/pyop2/compilation.py", line 186, in get_so
    return ctypes.CDLL(soname)
  File "/usr/lib/python2.7/ctypes/__init__.py", line 365, in __init__
    self._handle = _dlopen(self._name, mode)
OSError: /tmp/pyop2-cache-uid1000/774ed18f78b19aed133013130dc427a8.so: undefined symbol: form_cell_integral_0_otherwise

Indeed, 774ed18f78b19aed133013130dc427a8.c contains a form_interior_facet_horiz_integral_0_otherwise.

Cleaning cache, going back to the original code and replacing all instances of dS_v with dS_v(degree=(4, 4)), I get the previous This function is not capable of handling multiple integrals. error.

If I set the degree to (6, 6) for all different types of integrals (dx, dS_v, dS_h), I sporadically get the undefined symbol error, but the C file corresponding to the .so doesn't exist.

miklos1 commented 9 years ago

I think you should use the same measure object throughout the form. Try:

dx = dx(degree=(4, 4))

before the form instead of textually replacing everywhere.

wence- commented 8 years ago

Now possible after #691