How does the symbolic simplification in UFL works?

variational formulation
1

I have a axisymmetric problem with 1/r^2 singularities at the centerline. I have perused some of this forum’s thread’s on the subject :

I have resolved to do a transformation \phi\rightarrow r\phi for my test functions in order to take care of the singularity (this seems to amount to multiplying the variational formulation by r). However, this does not produce the expected improvement ; I don’t understand how UFL simplifies some expressions but not others. Here’s a first example :

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)

uD = dolfinx.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
uD.x.scatter_forward()

fdim = mesh.topology.dim - 1
# Create facet to cell connectivity required to determine boundary facets
mesh.topology.create_connectivity(fdim, mesh.topology.dim)
boundary_facets = numpy.where(numpy.array(dolfinx.cpp.mesh.compute_boundary_facets(mesh.topology)) == 1)[0]
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = dolfinx.DirichletBC(uD, boundary_dofs)

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
L=ufl.inner(f,v)*ufl.dx

problem1 = dolfinx.fem.LinearProblem(a1, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh1 = problem1.solve()

problem2 = dolfinx.fem.LinearProblem(a2, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh2 = problem2.solve()

L2_error = ufl.inner(uh1- uh2, uh1- uh2) * ufl.dx
print(numpy.sqrt(dolfinx.fem.assemble_scalar(L2_error)))

In Docker dolfinx with real mode this gives an error of about 10^{-15}, so machine precision. The same code with :

a1=ufl.inner(ufl.grad(u)/x[0]**2,ufl.grad(v))*x[0]*ufl.dx
a2=ufl.inner(ufl.grad(u)/x[0],     ufl.grad(v))*ufl.dx

Leads to much higher error ; about 10^{-4}. In my application, I think this imperfect simplification is my dominant source of error. Anyone has an idea on how to “force” UFL to simplify a form ?

2

The reason you get different results with the forms

a1=ufl.inner(ufl.grad(u)/x[0]**2,ufl.grad(v))*x[0]*ufl.dx
a2=ufl.inner(ufl.grad(u)/x[0],     ufl.grad(v))*ufl.dx

is because they cannot be integrated exactly using Gaussian quadrature on each element (because the integrand is no longer a polynomial), and UFL’s heuristic for estimating an appropriate quadrature degree is picking different degrees for a1 and a2. If you force the same quadrature degree for both,

dx_fixed_deg = ufl.dx(metadata={"quadrature_degree":4})
a1=ufl.inner(ufl.grad(u)/x[0]**2,ufl.grad(v))*x[0]*dx_fixed_deg
a2=ufl.inner(ufl.grad(u)/x[0],     ufl.grad(v))*dx_fixed_deg

you will see that the results are once again equivalent up to machine precision.

1 Like
3

If that was the only case I would be happy to take that answer. However, the equality in the first case :

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

Makes me think there has to be symbolic simplification somewhere. If it was simply a quadrature problem it would arise in this case too.

In my mind this is linked to Dealing with removable singularities in cylindrical coordinates - #7 by Xuefeng_LIU. In the 1/r singularity case it’s all good, but in the 1/r^2\cdot r it doesn’t work ?

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