geodynamics / aspect

A parallel, extensible finite element code to simulate convection in both 2D and 3D models.
https://aspect.geodynamics.org/
Other
217 stars 232 forks source link

Add distance weighted particle interpolator #5815

Closed gassmoeller closed 1 week ago

gassmoeller commented 1 month ago

This is not ready for review.

@anne-glerum: This PR adds a new particle interpolator that interpolates between particles based on a distance weighted average (similar to what radial basis functions would do). The interpolator in its current form produces worse errors (though the same convergence rate) as the 'cell average' interpolator. However, it should be smoother and less affected by small particle movements across cell boundaries.

This PR is currently in a very rough state, I only open this to remind myself to finish it after the hackathon.

gassmoeller commented 3 weeks ago

I cleaned up the PR and made sure I can run the particle benchmarks to show the method is converging. In terms of accuracy this interpolator has a similar convergence rate as 'cell average', but because the particle contribution is weighted by distance and particles of neighbor cells are considered as well, it should be less sensitive to particles crossing cell boundaries.

Ready for review.

gassmoeller commented 3 weeks ago

I think I have addressed all comments, but let me run some more benchmarks with this before merging.

gassmoeller commented 3 weeks ago

Ok, I am done benchmarking. I think this is ready. In the benchmark rigid_shear/transient it reaches the same convergence order as cell average at better accuracy: error_over_resolution_end

If anyone used cell average regularly before, I would now recommend using distance weighted average instead.

Edit: Except that distance weighted average is slow, even slower than bilinear least squares. This is likely just an issue with the implementation though, I can fix that in a follow up PR.

naliboff commented 2 weeks ago

@gassmoeller - thanks for the work to add this feature!

I did a bit of testing related to #4370 and ran into an issue when doing an initial adaptive refinement step on the first time step.

The model runs without issue if set Initial adaptive refinement = 0

The log.txt and a prm that be can quickly run on one processor are attached. Let me know how I can help with additional testing!

test_weighted_particle_average_refinement.prm.txt log.txt

gassmoeller commented 2 weeks ago

@naliboff: Sorry about that I clearly didnt do enough testing for this interpolator, adaptive refinement was broken. I think I fixed it now and added a test. Could you give your model another try? I encountered another problem with your model that might be a bug in deal.II. If you crash still try changing the number of particles per direction to 6 (from 7). I will look into that more.

naliboff commented 2 weeks ago

@gassmoeller - Thanks for the fix and the model runs successfully now with both adaptive refinement and either 6 or 7 initial particles per direction. Is there any other additional testing I can do in the short term? I'll look more carefully at the effect of the weighted average scheme on convergence on the complex models in the coming weeks.

gassmoeller commented 1 week ago

This is ready from my side. Just a quick note that at present this interpolator is quite slow (about twice as slow than the linear least squares interpolator, which does a much more complicated operation). Presumably this is because it has to create the std::set of neighbor cells for every cell. I have some ideas for how to speed this up, but it is not a terrible cost at the moment and the functionality should be fine as is, so I think we can merge this PR and address the performance in a follow up.