implement WENO limiter for DG method

[YiminJin]

This PR implements WENO (Weighted Essential Non-Oscillatory) limiter for DG method. The basic idea of WENO limiter is to replace the solution in troubled cells by a polynomial reconstruction that takes the neighbor cells into account. My implementation is based on the simple WENO scheme proposed by Wang and Shu, 2013 (https://doi.org/10.1016/j.jcp.2012.08.028). It is almost the same as the one in elASPECT.

The WENO limiter is able to smooth compositional fields without upper/lower limits, such as plastic strain and viscoelastic stress. I created a new file source/simulator/limiters.cc and moved both the BP limiter and the WENO limiter into this file. I also changed the parameter file so that users can choose whether to smooth a field by BP limiter or WENO limiter, or not to smooth it at all (See the sample prm file below).
kaus_2010.prm.txt

This PR was meant to solve the problem in #5734 . Unfortunately, the WENO limiter cannot improve the convergence rate in that case. Anyway, I think the WENO limiter offers more options to the users and is worth to be implemented.

[gassmoeller]

Hi @YiminJin: This is very interesting and I am happy to consider to add WENO as one of the available limiters to ASPECT. Before I do a full review of this PR, could you resolve the following problems wit this PR:

• the tester no-unity-build is failing, likely because you are missing necessary headers in limiters.cc (see here for the error message). You can reproduce this problem locally if you configure ASPECT with the cmake option -DASPECT_UNITY_BUILD=OFF -DASPECT_PRECOMPILE_HEADERS=OFF. Please make sure you include all necessary headers and you can build without unity build and precompiling headers
• The standard tester fails, because you removed two of the input parameters (the old parameters controlling the limiter). Usually we try to keep parameter files compatible (i.e. working without changes) for new ASPECT versions as much as possible, but we are currently reorganizing a lot for the upcoming release anyway, so we may as well include this. However, please implemented one of these two options:
  • Better: add a new python update script similar to the move_particles.py script that can automatically convert an old parameter file into the new format (removing the old parameter and - if necessary - adding the new parameter with the value of bound preserving
  • Less good: Keep the old input parameter for now, and make sure that if the old input parameter is set to true the new parameter is set to bound preserving.
• In order to merge the WENO limiter please add some tests and at least one benchmark that proves the limiter is working as intended. In the benchmark please also explain how WENO works and what the tuning parameters mean that you have added. I realize we do not have such a benchmark for the bound preserving limiter, which was an oversight at the time we implemented it. As a template you could use the benchmark for operator splitting in benchmarks/operator_splitting (choose a different setup that can illustrate the WENO limiter, maybe the viscoelastic_bending_beam benchmark? or some other setup that shows how WENO prevents oscillations)

I realize this will take some time and is a bit of effort, feel free to ask questions during the process. I quite like the reorganization of the files you did and think it is worthwhile to add different limiters, we just need to make sure the WENO results are included in the test suite and well documented.

Thanks for the contributions!

[gassmoeller]

I think it would need to be bound preserving. boundary preserving refers to a geometrical boundary while bound preserving mean preserving the upper and lower bounds (values) of a solution.

[YiminJin]

@gassmoeller Thank you for the detailed instructions. I have changed "boundary_preserving" to "bound_preserving" and added a python script to reformat the .prm file. I may need some time to design a suitable benchmark for the WENO limiter.

[YiminJin]

Hi @gassmoeller , I have added three benchmarks that show the performance of WENO limiter:

• The viscoelastic_bending_beam benchmark. With WENO limiter, the flow field can keep stable until 19000 yr, while the original viscoelastic_bending_beam.prm can only keep stable before 11000 yr. The results are shown in bending_beam.png (left).
• The viscoelastic_beam_modified benchmark. The models are still running, but from the velocity field we can see that the WENO limiter has greatly improved the stability. The results are shown in bending_beam.png (right).
• the rotate_shape benchmark. I added two benchmarks with BP limiter and WENO limiter, respectively. The results are shown in rotate_shape.png. As seen, the result of BP limiter is the best, the results of WENO limiter and EV method are comparable, while the result of SUPG method is the worst.

I think the original .prm files of bending beam benchmarks can be replaced by the new ones with WENO limiter, what do you think?

I have two questions:

• I don't quite understand the operator_splitting benchmarks you mentioned in the previous comment. They show the improvement in convergence rate when using operator splitting method, but cannot show the performance of stabilization methods (one of them is even without stabilization method).
• I haven't added test problems for the WENO limiter, because all the tests with limiters are for chemical fields, which is not suit for WENO limiter. Do you think it is OK to use those test problems, even if in those situations one should use the BP limiter?

I have also modified the file doc/manual/citing_aspect.bib to be able to cite the papers of WENO limiter in source/simulator/parameters.cc. Please check if I did it correctly.

[YiminJin]

I put a wrong picture for the rotate_shape benchmark with EV method in the previous comment. The one below is correct:

It should also be noticed that in this benchmark the EV method and the DG+WENO method both use their default stabilization parameters. Their performances may be improved if the stabilization parameters are optimized.

[tjhei]

This is great! Let me know if there is anything I can help with.

[gassmoeller]

Some more comments. I have not looked at limiters.cc itself yet.

To your questions:

I have added three benchmarks that show the performance of WENO limiter

Great, thank you! Please add the figures you showed to the documentation of those benchmarks.

I think the original .prm files of bending beam benchmarks can be replaced by the new ones with WENO limiter, what do you think?

If the new limiter performs better like you say then yes, definitely.

I don't quite understand the operator_splitting benchmarks you mentioned in the previous comment. They show the improvement in convergence rate when using operator splitting method, but cannot show the performance of stabilization methods (one of them is even without stabilization method).

I only meant if you want to see the structure of how to write such a benchmark you can use this benchmark as a template. Scientifically it is on a completely different topic of course.

I haven't added test problems for the WENO limiter, because all the tests with limiters are for chemical fields, which is not suit for WENO limiter. Do you think it is OK to use those test problems, even if in those situations one should use the BP limiter?

Can you explain why you say they are not suited for the WENO limiter? Is that because you expect to know the global bounds for the fields, because they are essentially indicator functions (0 or 1)? But that would only mean that BP is maybe better suited for those models, but WENO would still work to reduce oscillations, wouldnt it? (maybe a bit worse than BP). Anyway, if WENO works (independent of the question if it is the best limiter) for these models then you can certainly adapt a few of the existing BP tests to the WENO limiter.

I have also modified the file doc/manual/citing_aspect.bib to be able to cite the papers of WENO limiter in source/simulator/parameters.cc. Please check if I did it correctly.

Looks correct to me. You could make sure by writing a small manual page that explains what the limiting does to discontinuous compositional fields, similar to this page for the stabilization of the continuous fields.

• We also need a changelog entry in doc/modules/changes that announces that WENO is now available and explains where to find more information about WENO

[gassmoeller]

Please add an explanation what the KXRCF indicator is. What does it stand for and what does it mean to have a high KXRCF in a cell?

[gassmoeller]

As above, please extend this documentation by what the KXRCF indicator means and under what circumstances it makes sense to look at it.

[gassmoeller]

Have you tested this script? You will likely need to apply it anyway to all files that currently use the BP limiter in benchmarks/ cookbooks/ and tests/. An easy way to apply the script is to run the file contrib/utilities/update_scripts/update_all_files.sh.

[gassmoeller]

This is likely not true at the moment. You specify that this entry has to follow the pattern Double, i.e. just a single number. Either update the description, or update the input pattern (and add functionality to support different indicator thresholds for different fields).

Later edit: Oh, I saw below you do support different values for the indicator threshold for different compositions. Then you need to fix the Pattern in line 1218 to say Patterns::List(Patterns::Double(0.0)), otherwise no one will be able to actually put in different numbers.

[gassmoeller]

This is not quite how we handle distinctions between temperature and composition in SimulatorAccess at the moment. Ideally we would be able to hand over AdvectionField arguments, but they are currently part of the Simulator class, which is a problem (see #5887). Can you model this function after get_artificial_viscosity and get_artificial_viscosity_composition, i.e. create two functions for temperature and composition? It is not ideal, but until #5887 is fixed I would rather have all functions in SimulatorAccess solve this in the same way.

[anne-glerum]

Hi @YiminJin thank you for adding a limiter for fields without actual limits! I wanted to look at the implementation, but it seems hundreds of back up files were added to this PR (.bak extension) that should not be there.

[gassmoeller]

Hi @YiminJin, thanks for addressing my comments. As Anne wrote you accidentally added all the backup files that the update script created to this PR. This makes reviewing the PR very hard, because the large number of changed files freezes my browser window if I go into review mode. Could you remove all the files that end with .bak from the PR? Happy to give a review on the WENO implementation afterwards.

[YiminJin]

@anne-glerum @gassmoeller Sorry for the mistake. Now the .bak files are removed. In the previous commit I missed the case that there might be comments behind the parameter value, so some of the modified .prm files are incorrect. Now the mistakes are corrected, and the script contrib/utilities/update_scripts/prm_files/reformat_limiter_settings.py has been modified accordingly.

There are still more than 80 changed files due to conflicts. Is there some recommended methods to deal with the conflicts @gassmoeller ?

[gassmoeller]

Thanks for the update. I see a lot of changes that have already been applied to the main branch in #5941. I think the best way to get rid of the conflicts will be to rebase your branch. Your branch has accumulated a lot of commits and rebasing all of them will be quite painful. What I would personally do is start by squashing all commits on this branch into one. Then rebase to the latest version of the main branch. Something like:

git checkout main
git pull upstream main
git checkout weno
git rebase -i `git merge-base main HEAD`
git rebase main

(in the second to last step you will have to squash all of your commits)

[YiminJin]

@gassmoeller I am sorry to say that I have ruined this branch in the previous commit reformat the .prm file for limiter settings. I committed the changes without testing. Now I find that something went wrong after excuting the script contrib/utilities/update_scripts/update_all_files.sh. For example, the 27-th line of include/aspect/volume_of_fluid/handler.h changes from #include <aspect/volume_of_fluid/assembly.h> to #include <aspect/volume_of_fluid/simulator/assemblers/interface.h>, which is clearly incorrect.

I will reset this branch to one of the previous versions and update the .prm files manually. I suggest that the script contrib/utilities/update_scripts/update_all_files.sh is out of date and should not be used for now.

[gassmoeller]

Hi Yimin,
Dont worry, if the branch got messed up I would suggest you keep a copy of the parameter update script, and then just remove the commits that are related to the parameter file updates. It looks like this commit is the last one before the updates to the parameter files: f5cbc1c? If that is correct you can do the following while on your branch (careful): git reset --hard f5cbc1cb05f73165dc1c93a05c82fb43b9deaf67 this should get you back to the last state that was working. I would be ok with this PR if you only included the changes to the update script, but do not apply the script to all parameter files (we can do that separately later).

[YiminJin]

@gassmoeller Thank you for the instruction. I have reset the branch to f5cbc1c
and updated the .prm files manually. Besides, I have make 3 other changes:

1. There was a bug in contrib/utilities/update_scripts/prm_files/reformat_limiter_settings.py (reuse of loop indicator i) and it is fixed now;
2. I mistakenly assumed that the citation in parameter descriptions is in doc/citation_list.tex, but now I find that it is actually in doc/sphinx/references.bib. The references for WENO limiter have been moved there.
3. I added a test for WENO limiter.

Now the branch works fine. Please check it @gassmoeller @anne-glerum . If there is no other problem, I will do the rebase and deal with the conflicts.