Fix gradient in test_anisotropic

[jwallwork23]

While #381 re-enables the turbines in test_anisotropic, I found there's still something not right because when I use a similar setup for PDE-constrained optimisation (mesh-adaptation/opt_adapt#38), the optimisation progress is non-monotonic and the control turbine area occasionally goes to zero.

I added a Taylor test and found that the gradient of the QoI wasn't being computed correctly with the existing setup. I tracked the problem down to the turbine_density expression, which includes several terms. By projecting this expression into $\mathbb{P}1_{DG}$ space, the Taylor test passed.

[stephankramer]

I am a little confused by this test. It implements its own version of the bump function - which means it's really testing something that's not part of Thetis which has its own implementation of the bump function. Note there's a difference in this version, which cuts off the bump function at a radius r which is a little odd because the bump-function expression is not actually axisymmetric. Another potential difference is the way that quadrature degree is used: in Thetis we use a fixed quadrature degree (to avoid excessive estimated degrees) in the equations and also in the turbine diagnostics. For discrete tidal turbines there's an option, options.discrete_tidal_turbines.degree, which defaults to 10, which changes the quadrature degree in the diagnostics and SWE drag term. The way the test sets things up, because it uses tidal_turbines, which is really intended for a continuous density representation of a turbine farm, you probably have a quadrature degree that's too low where the drag is applied in the Thetis equations, as it will just use the same as for other terms, but then are using an automatic quadrature degree in the functional.

Finally, note that the added Taylor tests are not actually run: test_gradient() has taylor_test=True but also solve_adjoint=False so it returns before getting to the if taylor_test block. Fixing that the taylor test won't actually run because the model hasn't actually been taped (the adjoint is hand-rolled and doesn't use pyadjoint).

[cpjordan]

I think Stephan is probably going to be more suited to review for this one!

[jwallwork23]

Thanks for the comments @stephankramer. This test predates the discrete turbine implementation in Thetis. I'd misunderstood that the discrete turbines implementation only worked with meshed subdomains for the turbine footprints. Reading through the discrete turbine examples, I realised this isn't the case. I figured out how to rework this example to use that approach the other day, so will close this PR and reopen one that does it properly.

Thanks for pointing out the issue with the Taylor test. To be honest, the gradient computation is a bit out of place here so I think I'll move that stuff into a new example in the discrete turbines example subdir. The error estimation is also a bit out of place so I'll have a think about what to do with that.