Inaccurate linear solves problem

Hi team,

I am trying to solve the periodic Helmholtz problem: u:[0,1]^2\to\mathbb R

u_{xx}+u_{yy}-\sigma u = f,

where f is a nice periodic function, and sigma is rather small, e.g. sigma = 0.0001 or smaller.

The problem I am finding is that the accuracy of the scheme (measured via a reliable error estimator that is equivalent to the piecewise H^2 error) falls as low as 1e-07, but then increases to 1e-05 as the mesh size is reduced.

It also looks at though the residual r = |Ax-b|_{2} from solving the linear system is quite large, varying from 2e-08 to 1e-04, and it also increases at the mesh size decreases.

In the code snippet I also print out the first 10 values of: the right-hand side vector; the values of Ax; and the residual. There appears to be some discrepancy (so the increasing value of r is not just linked to the notion of summing more and more small values as the number of degrees of freedom increases).

In contrast, when sigma is say equal to 1, the scheme and the corresponding results are much better.

It is not clear to me if the Inverse() function is solving the linear system exactly, or just to some tolerance?

Any help would be greatly appreciated, and I have attached a mimimum working code example.

All the best,

Ellya

Hi Ellya,

it sounds to me as if there is a problem with your formulation. I'm not sure I'm able to solve your problem completely, but I have a few pointers which might help:

- You can try to solve the system more accurately by "re-iterating the system". You can do this by using the preconditioned Richards iteration and setting the direct solver as the "preconditioner"
- It looks like your trying a Nitsche formulation:

• Are you not missing the consistency (and symmetry) terms? (I've only worked with 2nd order problems personally, so don't take my word for it).
• From my experience, the penalty parameter you choose as 10, seems a but small
• As you have a fourth order problem, and the problem appears not to be solvable for small h, maybe it's worth checking out the h scaling in the derivation of your penalty term? Should it not be 1/h**2 to be consistent? (Then also k**4, I guess.)

Using the re-iteration, increasing the penalty term to 100 * k**4 /(h**2), the error estimator decreased in each step down to 10^-13.

As a sidenote: You can compute the 2-norm of the residual more easily with
I hope this helps you to solve the problem!

Best wishes,
Henry

Hi Henry,

Thank you for your help, I'll have a go with this reiteration, this is something I was trying to do so glad to hear it works!

I have a couple of questions/ responses:

1. What tolerance does the direct solve .Inverse() solve the linear system to, is it exact?
2. Does multiple iterations of the Richardson solve lead to even better accuracy, or does it only "work" once.
3. I didn't think to go that high with the penalty so I'll try that.
4. The scheme treats the PDE in strong form, so it is not missing any symmetry/ consistency terms, as no integration by parts is used to derive the method. This also means that the 1/h penalty is of the right order, since then the facet integral of (1/h)*jump(du/dn)^2 is of the same order as the element integral of |D^2u|^2. So perhaps I'll try boosting the penalisation coefficient.
5. Is k here the polynomial degree?

Thanks again for your help, it is really appreciated.

All the best,

Ellya

Hi Ellya,

I'll try and answer as best I can:

1. I would say that the exactness will depend on the sparse direct solver (factorisation) behind the .Inverse() and on the condition number of the system you're trying to solve. This solver will probably be either UMFPACK or Pardiso, depending on either how your NGSolve was built.
2. with the keyword argument "printing=True", you will see the behaviour of the residual after each iteration, which will tell you if multiple iterations help.
5. Sorry, yes k was the polynomial degree.

Best wishes
Henry

Hi Henry,

Please accept my apologies, as I had accidentally missed out the element error in the error estimator. The full estimator is still behaving poorly, and only the jumps part is small (which makes sense due to the large penalty). I'm not sure if the Richardson iteration is improving the linear solve (the linear system residual is still of the same order).

Is there any documentation on the linear solvers past the website examples?

All the best,

Ellya

Hi Ellya,

when I ran your example, the residual in the Richardson iteration got smaller by some orders of magnitude after one step, but then did not decrease after further iterations (using pardiso for the direct solver). My best guess is that the system is just too ill conditioned.

Sorry I can't be of more help.

Bets wishes,
Henry