Fixes for Firedrake 2025.10.1

[cpjordan]

It might be worth starting main/release branches mirroring Firedrake so that it is easier to track upstream changes

Fixes for Firedrake:

• sub_elements re-naming (#4354)
• PointNotInDomainError name due to .at() deprecation (#4543)
• removal of Interpolate.interpolate() (#4531)
• can no longer cast Constants+domain (#4497)
• must pass apply_riesz=True to get the gradient for the adjoint (see e.g. #3941)
• reaching max. iterations in optimisation now causes a SciPyConvergenceError (pyadjoint #218)

[stephankramer]

@cpjordan I hope you don't mind I just pushed an alternative fix for the Interpolator change, now actually reusing an interpolator. Motivation in the commit message, hope it makes sense

[cpjordan]

@cpjordan I hope you don't mind I just pushed an alternative fix for the Interpolator change, now actually reusing an interpolator. Motivation in the commit message, hope it makes sense

Thanks - that is what I should have done in the first place and what I wanted to do anyway but I was just getting rid of error codes instead of using my brain.

[cpjordan]

I think I was mostly confident on the other changes @stephankramer as I was just following the changes Firedrake had made (I've linked all the relevant PRs). The other bit was the change in what was returned by ReducedFunctional.derivative(), where in the inversion tools I have just set apply_riesz=True to keep it consistent with previous implementation (i.e. have the gradient returned):

self.set_initial_state(J, self.reduced_functional.derivative(apply_riesz=True), self.control_coeff_list)

but in the optimisation callbacks you are passed the derivative so you need to get the gradient slightly differently:

        derivatives_fn = [d.riesz_representation('L2') for d in derivatives]
        return [[float(d.dat.data_ro[0]) for d in derivatives_fn]]

I suppose for the inversion tools, you could do the same as the above.

Hopefully that's everything covered and then I'll tidy up my other PRs and do the main/release. I'm happy to have a call to discuss these when appropriate - might make things a bit quicker!

[stephankramer]

Yes so this Riesz representation is a bit funny for constant controls: it depends whether you think about them as just a constant or a function over the domain that is constant everywhere. For the turbines I think the first makes more sense, in which case we should actually be using the little l2 (which is actually do nothing in the sense that you get the same value array but instead of coefficients for the dual basis representing a Cofunction, you then use same coefficients for a linear combination of basis functions that forms the l2-Riesz representative) - which makes sense if you think about the vector space of controls (typically here the turbine coordinates) just being $$R^{2n}$$ (for n turbines) equiped with the standard Euclidean inner product. What we're doing now I think means that we actually multiply the derivative (i.e. the actual partial derivative of the functional w.r.t. to each individual individual coordinate) with the total area of the domain. For the optimisation it probably doesn't really matter, because it's equivalent to just a rescaling of your coordinate space (it might matter a little if you do gradient descent, or for the initial scaling in BFGS).

So come to think it, I actually think that it makes more sense to use l2 in the DerivativeConstantControlOptimisationCallback in which case you don't need to do a riesz map at all you can immediately pick the value from .dat.data[0] of the Cofunction in R space. I think at the moment that's only for diagnostic output anyway.

For the inversion tool: yes, let's keep the current functionality (i.e. apply whatever riesz map is specified when the control is created), but then we probably (in a future PR) want to add the functionality that we can choose which one that is in add_control(). Because the choice may indeed depend on whether it's just a coefficient, or spatial field - but in particular it also depends on whether the optimisation algorithm actually uses the correct inner product. For all scipy algorithms it doesn't do that, so there you might argue that it's actually better to stick with 'l2' (as that's the inner product that scipy will be using internally) even for spatial fields.

[stephankramer]

Anyway, I'll leave that up to you (whether you want to change the behaviour of DerivativeConstantControlOptimisationCallback) - for the rest I agree with everything you've done. So let's get this in first, which at the moment should work with Firedrake release and main, and then afterwards we can worry about setting up a main and release Thetis branch (and change CI accordingly).

[thangel]

Looks good, thanks for implementing these