A simple question about a boundary condition

Hi,

I have an influx boundary condition as follows:

\frac{\partial A}{\partial {\bf n}}= \kappa (A-A_0)

where \kappa=1 and A_0=1 are constants and {\bf n} is the unit outward normal vector. There are several ways this can be incorporated in the weak form. For a function space V

#Case 1
A = Function(V)
v = Testfunction(V)
A0 = 1
kappa = 1 
weak_form = kappa*(A-A0)*v*ds

#Case 2
A = Function(V)
v = Testfunction(V)
A0 = Constant(1)
kappa = 1
weak_form = kappa*(A-A0)*v*ds

#Case 3
A = Function(V)
v = Testfunction(V)
A0 = project(Constant(1),V)
kappa = 1
weak_form = kappa*(A-A0)*v*ds

#Case 4
A = Function(V)
v = Testfunction(V)
A0 = interpolate(Constant(1),V)
kappa = 1
weak_form = kappa*(A-A0)*v*ds

Now, all of these cases work for any problem, and I receive no errors. However, the results can sometimes be slightly different. I was wondering which one of these is better and why.
Thank you in advance for your help.

It looks like only the definition of A0 varies between cases. I would favor Case 2, because Cases 3 and 4 introduce unnecessary complexity to represent a constant value, while Case 1 would need to recompile any time the numerical value of A0 is changed.

In terms of the results, these should all be equivalent. Are the differences small enough to attribute to floating point round-off error? If not, the only other source of differences I can think of would be if UFL is estimating different quadrature degrees for the different cases. To debug, you might try all the cases with a fixed quadrature degree by setting

ds = ds(metadata={"quadrature_degree":2*k})

where k is the degree of the function space V. However, the integrals here are simple enough that they should always be integrated exactly (up to floating point round-off error) by the default quadrature rule.

2 Likes

Thank you so much. The experiment gave similar results. But to be safe, I will also stick with case 2.