Instabilities with Summation by Parts thorn 4th order operator

[Miguel Zilhão]

I’ve noticed that instabilities develop when using the 4th order accurate 2nd derivative Summation by Parts operator, from the SummationByParts thorn. I’ve coded a simple wave equation evolution thorn that demonstrates the problem, which can be found here: https://github.com/mzilhao/IST_WaveToy. The thornlist used to checkout and compile this thorn is attached.

2nd derivative SBP operators of order 2, 6 and 8 work fine. With order 4, however, there are instabilities propagating from the boundary that crash the code very early on. The code also seems to run fine with the parameter “SummationByParts::sbp_2nd_deriv = no”.

In the repository I've also added the parameter files tried with the different SBP orders (under ./par).

[Peter Diener]

I’ve checked the optimized diagonal 4th order second derivative operator against the mathematica notebook originally used to derive it and didn’t find any mistakes. However, the 4th order 2nd derivative operator is the only 2nd derivative operator that has both a minimal bandwidth and an optimized operator. The other 2nd derivative operators have only minimal bandwidth operators. In order to explain further, when the SBP operators are derived many have free parameters that can be chosen freely. For the 1st derivative diagonal operators the 2-1 and 4-2 operators are unique, while the 6-3 and 8-4 operators have free parameters. One way to choose the parameters is to set as many coefficients to zero (called minimal bandwidth) , another is to choose the parameters so the sum of the squares of the truncation error coefficients is minimized (called optimal). For the 2nd derivative operators the 4-2, 6-3 and 8-4 operators have free parameters. For some reason an optimized operator was only found for the 4-2 operators, while only minimal bandwidth operators were implemented for the 6-3 and 8-4 cases. Default parameter choices selects the optimized operators where implemented but uses the minimal bandwidth operators where they are not. I tried used the minimal bandwidth diagonal 2nd derivative operator (by setting SummationByParts::operator_type = "Minimal Bandwidth") and that seems to be stable.

It should be noted that stability of any of these operators are not guaranteed (i.e. proven) in multiple dimensions, so it is not impossible that SPB satisfying operators leads to unstable evolutions. However, it is surprising to me that using the free parameter to minimize the sum of squares of the truncation error coefficients leads to such dramatically different stability properties.

Miguel, can you confirm that you also get stable evolution for the minimal bandwidth 4th order diagonal operator by setting SummationByParts::operator_type = "Minimal Bandwidth"?

[Miguel Zilhão]

Thanks for looking into this. Indeed, I can confirm that setting SummationByParts::operator_type = "Minimal Bandwidth" with the 4th order accurate 2nd derivative operator seems to lead to a stable evolution. Should this maybe be the default choice?

With the optimal stencil, the instabilities seem to grow from the boundary. Could this somehow be due to a mismatch between the truncation error in the bulk and boundary regions?

‌

[Peter Diener]

I don’t think I want to change the default operator_type to “Minimal Bandwidth” as this would then also be the default for the 1st derivative operators, where the optimized operators lead to stable (and more accurate) evolutions. What I’m contemplating is just disabling the optimized 4th order second derivative operator and (similarly to the higher order second derivatives) just use the minimal bandwidth, regardless of parameter settings. I believe you might be the first to try to actually use the second derivative operators, so such a change would not affect many people. I will attempt to make a pull request for that.

With SBP operators, there’s always a mismatch between the truncation errors in the interior and the boundary as the order is only half at the boundary compared to the interior. And in fact the optimized operator should have less overall truncation error compared to the minimal bandwidth operator. So I can’t explain why it is less stable. Then again, to have a proper proof of numerical stability in 1 dimension requires that the boundary conditions are implemented using something like a Simultaneous Approximation Term (SAT) methods. It is possible that the optimized operators would be stable if such a boundary method was used.

[Roland Haas]

From minutes:

Roland asked Peter D. about the status of the summation by partos ticket. Peter reported that by using the minimal bandwidth form, there will be minimal change to the code. After a discussion on how to report this change -- as deprecated feature or bug, decision was made to report it as deprecation and include it in the master branch after release.

which will be addressed after the ET_2020_05 release. A deprecation warning has been added to the release notes: http://einsteintoolkit.org/about/releases/ET_2020_05_announcement.html

[Roland Haas]

Since the release has passed, @Peter Diener do you want to update the default parameter in param.ccl or do whatever was required to make the minimal version of the operators being the ones used as you described above?

[Peter Diener]

The pull request is here:
https://bitbucket.org/cactuscode/cactusnumerical/pull-requests/13/dont-call-the-optimized-4th-order-second

‌

[Roland Haas]

Looking at the pull request and based on discussion in today’s ET call, I have the following review comments:

• since now the if and else branches of the 2nd derivatives are indentical, I would remove the if statement and replace it by a comment. Usually I would suggest replacing it by a check that the passed option is either “minimal bandwidth” or “optimal” choice for the operators, however the current code already does “minimal bandwith” or “other” so there seems little point of adding such a check now.
• I would advocate removing the now dead code from the repository, anyone who needs it can get it back from the version control system
• I would check that the documentation is updated if it mentioned the existence of an “optimal” 2nd derivative operators

Fine to apply otherwise.

‌

[Roland Haas]

@Peter Diener any progress on this?

[Roland Haas]

@Peter Diener any progress on this?

[Roland Haas]

@Peter Diener any progress on this?

[Peter Diener]

I have pushed changes that should address the comments to the pull request.