Enable mass lumping for projection operator

[jwallwork23]

Closes #95.

This PR allows the projection operator (and its adjoint) to be mass lumped. For P1 spaces, this implies that extrema are preserved following the transfer.

While developing these changes, I noticed #113 and #114, so the associated functionality and tests are turned off here.

[jwallwork23]

Still work to do to address the issues referenced above, but hopefully this is sufficient for your purposes @ddundo.

[ddundo]

Thanks very very much @jwallwork23! And sorry for the issues you raised - I tried to warn you about my incompetence in mesh-adaptation/goalie#128 (comment) :)

This actually makes quite a significant difference in my case, which I didn't expect. I will bring it up on the next meeting :)

[jwallwork23]

TODO: Make sure project is the default, not interpolate.

[jwallwork23]

Hopefully this is useful @ddundo! New minimal_diffusion option. A few tweaks still required, including plotting. However, I thought you might find the code of use.

[jwallwork23]

• Project alone means you can conserve mass.
• Lumped project means (additionally) no new extrema.
• This post-processing step then reduces the numerical diffusion as much as possible.

All follows (Farrell et al., 2009).

[ddundo]

Thanks a lot for doing this @jwallwork23! Looks great from what you show here!

I couldn't use it directly in my case because I have a vector field, so I tried to modify it by considering each component as a separate scalar field, but then didn't get very good results. So I assume something more sophisticated is necessary :)

As a side note... I am also confused how tests pass for you here since you didn't merge the commit for #117 here... I had to do it locally.

[jwallwork23]

Thanks a lot for doing this @jwallwork23! Looks great from what you show here!

No worries - I had a long train journey and this was an interesting puzzle to solve during it.

I couldn't use it directly in my case because I have a vector field, so I tried to modify it by considering each component as a separate scalar field, but then didn't get very good results. So I assume something more sophisticated is necessary :)

Ah yes, the vector extension will have to come in a later PR I think. Shouldn't be too hard, though.

As a side note... I am also confused how tests pass for you here since you didn't merge the commit for #117 here... I had to do it locally.

Hm yeah good point. Pass.

[jwallwork23]

I'll probably rename the minimal_diffusion kwarg as bounded (as in the paper) and drop the lumped option. (I don't think there would be a case where you'd want to use lumping on its own without the post-processing.)

[jwallwork23]

@ddundo This is now ready for review. Minimal diffusion correction taken out for now and marked as future work in #124. Note also #123.

Do you know how to get the colourbars formatted correctly?

[ddundo]

Thanks a lot for this @jwallwork23! The problem with the colourbars is that the levels in tricontourf and ticks in colorbar are different. So if you do

fig, axes = plt.subplots(ncols=2, nrows=2, figsize=(10, 10))
levels = [-0.05] + [0.15 * i for i in range(1, 8)]
for ax, f in zip(axes.flatten(), (sensor, f_interp, f_proj, f_bounded)):
    im = tricontourf(f, axes=ax, levels=levels)
    ax.axis(False)
    ax.set_aspect(1)
    fig.colorbar(im, ax=ax, ticks=levels)
# fig.colorbar(im, ax=axes, ticks=levels)  # or replace above with this for a single colourbar
plt.savefig("ping_pong-final.jpg")

you get this

[ddundo]

Thanks very much again for this! I just left minor comments :)

[jwallwork23]

Thanks for the thorough review @ddundo. I think I addressed all points in b775043.

[ddundo]

Looks good! Sorry for nitpicking :)

Oh and I also ran the tests both on my local install and on the docker container, and tests pass for me. I also updated firedrake.

[jwallwork23]

Thanks again @ddundo. Pushed a minor change that fixes the CI, not controversial.