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(a)gmail.com> wrote:
> Andrew, where would I do that?
>
>
> On Friday, August 5, 2016, Andrew McRae <A.T.T.McRae(a)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(a)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-alp
>>> ha_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
>>>
>>
>>