How does the symbolic simplification in UFL works?

variational formulation
4

There are different quadrature rules in play, as illustrated by:



from ufl.algorithms.compute_form_data import estimate_total_polynomial_degree

import dolfinx, ufl
from mpi4py.MPI import COMM_WORLD
import numpy

mesh=dolfinx.UnitSquareMesh(COMM_WORLD,10,10)
FE=ufl.FiniteElement("Lagrange",mesh.ufl_cell(),2)
V=dolfinx.FunctionSpace(mesh,FE)
u=ufl.TrialFunction(V)
v=ufl.TestFunction(V)

x = ufl.SpatialCoordinate(mesh)

f=dolfinx.Constant(mesh,-6)
a1=ufl.inner(ufl.grad(u)/x[0],ufl.grad(v))*x[0]*ufl.dx
a2=ufl.inner(ufl.grad(u),     ufl.grad(v))*ufl.dx
a3=ufl.inner(ufl.grad(u)/x[0],ufl.grad(v))*ufl.dx
a4=ufl.inner(ufl.grad(u)/x[0]**2,     ufl.grad(v))*x[0]*ufl.dx

import basix
for a in [a1,a2,a3,a4]:
    degree = estimate_total_polynomial_degree(a)
    quadrature = basix.make_quadrature(
            basix.quadrature.string_to_type("default"), basix.cell.string_to_type("triangle"), degree)
    print(degree, quadrature[0].shape)

yielding:

4 (6, 2)
2 (3, 2)
3 (6, 2)
5 (7, 2)

UFL does no symbolic simplification. The differences that arises in your problem is likely due to the singularity at x[0]=0.
Inside the integration kernel, a1=ufl.inner(ufl.grad(u)/x[0],ufl.grad(v))*x[0]*ufl.dx yields a temporary variable on the form
double a = x0/x0;
which is always equal to one, while
a1=ufl.inner(ufl.grad(u)/x[0]**2,ufl.grad(v))*x[0]*ufl.dx
yields a variable on the form
double a = x0/(x0*x0);.
As the quadrature points is getting close to the singularity, you might hit floating precision issues, as x0*x0 can become smaller than machine precision.
Different quadrature rules get within different proximity of the singularity, and can hence yield slightly different results.

1 Like