Accuracy of projection of UFL expression

[afalaize]

Dear Fenics users;

I am facing a behavior I did not expect when projecting UFL expression onto given function space for function assignement.

a (not so) minimal example to reproduce can be found here. Steps are

1. define a symbolic expression with Sympy for scalar function
2. use Sympy codegen capabilities to build associated Dolfin Expression and assign to given Dolfin Function, say “u”, by projecting onto FunctionSpace
3. compute rho = 1 + (dot(grad(u), grad(u))**0.5 by three methods, namely
   3.1) using symbolic evaluation of grad(u) with Sympy and dolfin.project,
   3.2) using numerical evaluation of grad(u) with dolfin.grad and dolfin.project,
   3.3) using numerical evaluation of grad(u) with UFL syntax “as_vector(u.dx(i), i)” and dolfin.project.

Results are quite surprising, as seen from printing of min and max values:

u: min=-2.9994028324045994e-23, max=0.9999999999878484
rho sympy: min=0.9999999999999971, max=63.84829467210787
rho dolfin: min=0.9999999999999982, max=70.89190994721339
rho ufl: min=-2.5106197601269207, max=70.84901819962683

Especially, min value for the method 3.3 is far below 0 while it should be (around) 1.

Could any one explain what I am missing?

Best regards;

Antoine Falaize

[nate]

This isn’t surprising. You’re projecting the gradient to a less rich function space in your dolfin implementation than the exact solution you compute for your sympy expression by symbolic differentiation. Consider enriching the approximation spaces in dolfin by increasing the degree of the approximating polynomial.

Furthermore, your pseudo L_\infty norms aren’t a great method for measuring the projection error. Consider instead measuring in the L_2 and H^1 norms.

[afalaize]

Thank you for reply.

I’am not sure to understand. I detail my question below.

Analytic evaluation

• Target function is \rho(x) = 1 + \sqrt{\nabla u(x) \cdot \nabla u(x) } .

• Implementation of u from Sympy analytic expression:

    order = 2    # usually 2 or 3     
    V_u = FunctionSpace(mesh, "Lagrange", order)
    u = Function(V_u)
    u_c_code = ...  # C code generated from sympy analytic expression
    u_dolfin_expr = Expression(u_c_code, element=u.ufl_element(), ...)
    u.assign(project(u_dolfin_expr, u.function_space()))
  

• Implementation of \rho with Sympy evaluation of derivatives (\rho is one order of approximation lesser than approximation order of function u due to spatial derivative):

    V_rho = FunctionSpace(mesh, "Lagrange", order-1)
    rho = Function(V_rho)
    rho_c_expr = ... # C code generated from analytic Sympy expression
    rho_dolfin_expr = Expression(rho_c_expr, element=rho.ufl_element(), ...)
    rho.assign(project(rho_dolfin_expr, rho.function_space()))
  

Numerical evaluation

• Implementation with UFL expression:

    i = Index()
    rho_ufl_expr = 1 + pow(dot(as_vector(u.dx(i), i), as_vector(u.dx(i), i)), 0.5)
    rho_ufl = Function(V_rho)
    rho_ufl.assign(project(rho_ufl_expr, rho_ufl.function_space()))
  

• Implementation with prior evaluation of \nabla u(x) used in a dolfin expression:

    V_grad = VectorFunctionSpace(mesh, "Lagrange", order-1)
    grad_u = Function(V_grad)
    grad_u.assign(project(as_vector(u.dx(i), i),  V_grad))
    rho_grad_expr = Expression(
        "1+pow(pow(grad_u[0],2) + pow(grad_u[1],2), 0.5)",
        degree=grad_u.ufl_element().degree(),
        grad_u=grad_u)
    rho_grad = Function(V_rho)
    rho_grad.assign(project(rho_grad_expr rho_dolfin.function_space()))
  

Results

The min and max values for rho_sympy, rho_ufl and rho_grad are

rho_sympy: min=0.9999999999999971, max=63.84829467210787
rho_ufl: min=-2.5106197601269207, max=70.84901819962683
rho_grad: min=0.9999999999999982, max=70.89190994721339

Questions

• I expected the implementation of rho_ufl and rho_grad to use the same UFL mechanism and numerical evaluation under the hood. Why are they so different?
• The implementation of rho_ufl with order=2 yields a minimum value of -2.51 for the function 1 \leq \rho(x) = 1 + \sqrt{\nabla u(x) \cdot \nabla u(x) }. How is it possible?

Maybe I should edit my original message…

Best regards;

Antoine