H1 Error differences between pde and py interfaces

[Dow Drake]

Using Netgen 6.1-dev Ngsolve 6.1

We're seeing some significant differences (up to a factor of 10 or so) in the H1 error calculated by ngsolve when switching between the pde file and python file interfaces.

The tests I did involve solving the Poisson equation on an H1 space. I varied the order, the mesh size and the type of element (square vs. triangular).

Using a square element and an H1 order 3 space, the computed errors are:
python 4.235E-04
pde 3.376E-03

Please see my detailed post here: https://sourceforge.net/p/netgen-mesher/discussion/905307/thread/a335ae92/ (I accidentally posted to the wrong forum)

The differences are most noticeable for square elements with a space of order 3, but there are also significant differences for triangular elements when the order is 4.

I would like to know which of these methods is recommended for best accuracy, and if there something we can fix in our implementation to make the errors consistent across the interfaces?

Thanks!

[Christoph Wintersteiger]

Hi,

there is one point you have to be aware of. The python Integrate(...) hat a default order for the used integration rule, which is 5. Thus it works for spaces of order 2, since you integrate a squared quadratic function. For a space with order=3 you need exact integration of 6th order polynomials, which is not provided by the standard implementation of the python Integrate(...) function. But you can increase the order of exact integration in the following way:

print ("L2-error:", sqrt (Integrate ( (u-exact)*(u-exact), mesh, order=6)))

print ("H1-error:",sqrt(Integrate((grad(Uh)- gradU)*(grad(Uh)-gradU), mesh, order=6)))

That's the best way to calculate the errors in python. The results are:
L2-error: 7.04e-05
H1-error: 3.375e-03

A few comments on your pde-version:
You can use the numproc difference to calculate the L2-error as well. It has to be done the following way:

bilinearform m -fespace=v -symmetric -nonassemble
mass 1

fespace verr -l2ho -order=0
gridfunction errl2 -fespace=verr
numproc difference npdiffl2 -bilinearform1=m -solution=u -function=uex -diff=errl2

Using the numproc difference I get the results:
Call numproc Calc Difference npdiffl2
difference = 8.81417e-05
Call numproc Calc Difference npdiffh1
difference = 0.00337643

If you use the numproc integrate you have the calculate the squared norm of the error and take the sqrt afterwards. The way you did it is not the norm you wanted!

coefficient errl2 ((u-uex)*(u-uex))
numproc integrate  calcL2err -coefficient=errl2  # squared norm!!
coefficient errh1 ( (grad_u-graduex)*(grad_u-graduex) )  
numproc integrate  calcH1err -coefficient=errh1  # squared norm!!

If I do so, I get these results:
L2-error: 8.731e-05
H1-error: 3.378e-03

I did the tests on the "squareunif.vol.gz" mesh and with a order 3 space and the results match very good. So I guess it was just the order of the Integrate function and the sqrt inside the numproc integrate in your pde-file.

Christoph

[Christoph Wintersteiger]

A way to find such default arguments is to call help of the function.
help(Integrate) gives the following output:

Help on built-in function Integrate in module ngsolve.comp:

Integrate(...)
    Integrate( (CoefficientFunction)cf, (Mesh)mesh [, (VorB)VOL_or_BND=ngsolve.comp.VorB.VOL [, (int)order=5 [, (bool)region_wise=False [, (bool)element_wise=False]]]]) -> object
(END)

The arguments in the square brackets are optional and they have default values. By default the function calculates a volume integral with exact integration up to order=5 and so on.

Christoph

[Dow Drake]

Christoph,

Thank you very much for taking the time to explain all this and pointing out how I can correct my code!

Best,
Dow