If you're planning to use the tet code (UnitCubeMesh) then you should do the same. I guess you snuck just under the 1000-point error check on that, due to small differences in degree estimation. Without a check like that, users would be none-the-wiser and just complain "Firedrake is slow" :) On 5 August 2016 at 23:09, Justin Chang <jychang48@gmail.com> wrote:
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-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
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake