That did the trick, thanks Andrew

On Fri, Aug 5, 2016 at 5:01 PM, Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:

If you want your a quadrature rule which is exact for *form* polynomial degree 5 in the horizontal and 3 in the vertical, use dx(degree=(5, 3)).

UFL symbolically analyses the form to obtain the polynomial degree, and by default uses quadrature which is exact for this.  When non-polynomial operations are present, it uses a heuristic, and when lots are present, it tends to be far too conservative.

E.g. deg(a/b) is estimated as deg(a) + deg(b), and deg(sqrt(f)) is estimated as 2 + deg(f).  These are reasonable on their own, but now your classic unit vector v/sqrt(dot(v, v)) has degree
deg(v) + deg(sqrt(dot(v, v)))

= deg(v) + 2 + deg(dot(v, v))

= deg(v) + 2 + deg(v) + deg(v)

If deg(v) = 1, this is 5(!)

In your case, I guess it's gone for a quadrature rule with 12 points in each direction, so it's estimated the total polynomial degree as 22 or 23... probably far more than you want.

On 5 August 2016 at 22:52, Justin Chang <jychang48@gmail.com> wrote:

Andrew, where would I do that?

On Friday, August 5, 2016, Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:

Set a manual quadrature degree, you probably don't want 1728 points per cell.

On 5 August 2016 at 22:47, Justin Chang <jychang48@gmail.com> wrote:

Hi all,

I have the following 3D advection diffusion problem with SUPG formulation:

# Mesh

meshbase = UnitSquareMesh(seed,seed,quadrilateral=True)
mesh = ExtrudedMesh(meshbase,seed)
V = FunctionSpace(mesh, 'CG', 1)
Q = VectorFunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
sol = Function(V, name="solution")

# Velocity

velocity = interpolate(Expression((  ... )),Q)

# Dispersion tensor

aT = 0.0001
aL = 1
aD = 0.0001
alpha_t = Constant(aT)
alpha_L = Constant(aL)
alpha_D = Constant(aD)
normv = sqrt(dot(velocity,velocity))
Id = Identity(mesh.geometric_dimension())
D = (alpha_D + alpha_t*normv)*Id+(alpha_L-alpha_t)*outer(velocity,velocity)/normv

# Forcing function

f = Constant(0.0)

# Weak form

h = CellSize(V.mesh())
Pe = h/(2*normv)*dot(velocity,grad(v))
a_r = Pe*(dot(velocity,grad(u)) - div(D*grad(u)))*dx
a = a_r + v*dot(velocity,grad(u))*dx + dot(grad(v),D*grad(u))*dx
L = (v + Pe)*f*dx

...

But when I try to run this code i get the following error:

Traceback (most recent call last):
  File "3D_advection_diffusion_ex1.py", line 122, in <module>
    execfile('driver/AD_solver.py')
  File "driver/AD_solver.py", line 36, in <module>
    solver_parameters=solver_parameters, options_prefix="linear_")
  File "/home/justin/Software/firedrake/src/firedrake/firedrake/variational_solver.py", line 261, in __init__
    super(LinearVariationalSolver, self).__init__(*args, **kwargs)
  File "/home/justin/Software/firedrake/src/firedrake/firedrake/variational_solver.py", line 127, in __init__
    ctx = solving_utils._SNESContext(problem)
  File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving_utils.py", line 107, in __init__
    for problem in problems)
  File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving_utils.py", line 107, in <genexpr>
    for problem in problems)
  File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", line 67, in assemble
    inverse=inverse, nest=nest)
  File "<decorator-gen-295>", line 2, in _assemble
  File "/home/justin/Software/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper
    return f(*args, **kwargs)
  File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", line 101, in _assemble
    inverse=inverse)
  File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 288, in compile_form
    number_map).kernels
  File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", line 198, in __new__
    obj = make_obj()
  File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", line 188, in make_obj
    obj.__init__(*args, **kwargs)
  File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 210, in __init__
    tree = tsfc_compile_form(form, prefix=name, parameters=parameters)
  File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", line 54, in compile_form
    kernels.append(compile_integral(integral_data, fd, prefix, parameters))
  File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", line 198, in compile_integral
    quadrature_degree))
  File "/home/justin/Software/firedrake/src/tsfc/tsfc/fiatinterface.py", line 108, in create_quadrature
    raise RuntimeError("Requested a quadrature rule with %d points" % len(rule.get_points()))
RuntimeError: Requested a quadrature rule with 1728 points

If I had used a UnitCubeMesh I had no problems. What is going on here and how do I fix/circumvent this issue?

Thanks,

Justin

_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake