Assemble_scalar returns different value for the same function represented two different ways

General
1

Hello,

In one of my simulations, I noticed some potential inaccuracies after integration and I am trying to understand. Please consider the following simple code that performs integration on the same function represented two different ways. Method A is the interpolation of a function and then integrating, Method B is directly integrating a non-interpolated function. Can someone please help me understand the difference between the two cases? I expect them to be the same!

from dolfinx import mesh, fem, nls, plot, cpp
import ufl

# setup
nx, ny = 20,20
domain = mesh.create_rectangle(comm=MPI.COMM_WORLD, points=((-1., -1.), (1., 1.)),n=(nx, ny), cell_type=mesh.CellType.triangle)

W = fem.VectorFunctionSpace(domain, ("CG",1), dim = 2)
W0,_ = W.sub(0).collapse()

# create some functions
F = fem.Function(W)
f1, f2 = ufl.split(F)
F.sub(0).interpolate(lambda x: x[0]**2)
F.sub(1).interpolate(lambda x: x[1]+2)

# Two ways to represent the same function
A = fem.Function(W0)
B = f1/f2
A.interpolate(fem.Expression(B,W0.element.interpolation_points()))

# Look at the integration:
# Method A:
fem.assemble_scalar(fem.form(A*ufl.dx))
0.7365659808912665
# Method B:
fem.assemble_scalar(fem.form(B*ufl.dx))
0.7360701674392127
2

Dividing a CG-1 function by another CG-1 function is not in CG-1, and thus by interpolating prior to integrating, you introduce an error.
Consider the following:

  • Use a CG2 space to have each of the polynomials being represented exactly, then the integral should be equal to \frac{2\log 3}{3}.
import numpy as np
from dolfinx import mesh, fem
import ufl
from mpi4py import MPI
# setup
nx, ny = 20, 20
domain = mesh.create_rectangle(comm=MPI.COMM_WORLD, points=(
    (-1., -1.), (1., 1.)), n=(nx, ny), cell_type=mesh.CellType.triangle)

W = fem.VectorFunctionSpace(domain, ("Lagrange", 2), dim=2)
W0, _ = W.sub(0).collapse()

F = fem.Function(W)
f1, f2 = ufl.split(F)
F.sub(0).interpolate(lambda x: x[0]**2)
F.sub(1).interpolate(lambda x: x[1]+2)

# Two ways to represent the same function
A = fem.Function(W0)
B = f1/f2
A.interpolate(fem.Expression(B, W0.element.interpolation_points()))

# Look at the integration:
# Method A:
A = fem.assemble_scalar(fem.form(A*ufl.dx))
# Method B:
B = fem.assemble_scalar(fem.form(B*ufl.dx))

exact = 2 * np.log(3)/3
print(f"{exact=} {A=} {B=} {A-exact=:.2e} {B-exact=:.2e}")
3

Thank you. The behaviour is a little different if we consider the case of constants, still with CG1. Printing the error as you did, I would expect to be a couple orders lower for machine precision/exact integration. No?

from dolfinx import mesh, fem, nls, plot, cpp
import ufl
from mpi4py import MPI

# setup
nx, ny = 20,20
domain = mesh.create_rectangle(comm=MPI.COMM_WORLD, points=((-1., -1.), (1., 1.)),n=(nx, ny), cell_type=mesh.CellType.triangle)

W = fem.VectorFunctionSpace(domain, ("CG",1), dim = 2)
W0,_ = W.sub(0).collapse()

# create some functions
F = fem.Function(W)
f1, f2 = ufl.split(F)
F.sub(0).x.set(2)
F.sub(1).x.set(2)

# Two ways to express the same function
A = fem.Function(W0)
B = f1/f2
A.interpolate(fem.Expression(B,W0.element.interpolation_points()))

# Look at the integration:
# Method A:
A = fem.assemble_scalar(fem.form(A*ufl.dx))
# Method B:
B = fem.assemble_scalar(fem.form(B*ufl.dx))

exact = 4.
print(f"{exact=} {A=} {B=} {A-exact=:.2e} {B-exact=:.2e}")

exact=4.0 A=3.9999999999999454 B=3.999999999999812 A-exact=-5.46e-14 B-exact=-1.88e-13
4

This here is quite dangerous, it would not work if you wanted to set different numbers to each sub-space.

Machine precision (1e-16) is rarely achieved when one do sums, divisions etc. The error can accumulate.
Say we have number a and b, which is the same up to 10e-16, i.e. a = c + A \cdot 10^{-16}, b= c + B\cdot 10^{-16}, then

\frac{a}{b} = \frac{c(1+ A/c\cdot 10^{-16})}{c(1+B/c\cdot 10^{-16})}=\mathcal{O}(1)

This error can then accumulate when summing up values.

1 Like