Add headland inversion example #371
Conversation
- t passed to cost_fn(t) was out of sync by one timestep so the optimal result was not correct for channel inversion - Now need to pass the cost function through the export_func rather than through update_forcings - Typo fixed in Tohoku Makefile as well
- Details are in the README - inversion of Manning field using different field representations:
- Constant
- Region-based
- Independent Points Scheme
- Free specification at all points in the domain (Hessian regularised)
- Also fix hanging in parallel when using weighting of stations in inversion examples
…hetisproject/thetis into add-headland-inversion-example
|
There are three aspects of this pull request which is intended to provide another example for model calibration using the adjoint-based approach:
Feedback would be appreciated as to whether this is valuable and anything to change. My thoughts:
Please tag anyone else who this might concern/interest! |
| @@ -15,6 +16,32 @@ | |||
| import os | |||
|
|
|||
|
|
|||
| class CostFunctionCallback(DiagnosticCallback): | |||
There was a problem hiding this comment.
Perfect. Maybe we should provide something like this in a more generic way: there's nothing about cost_function that means it specifically has to be a cost function, I mean any other user defined callback function could use exactly the same code - but let's deal with that potentially in a different PR.
There was a problem hiding this comment.
Absolutely. I can update this in a later PR!
| @@ -0,0 +1,827 @@ | |||
| import firedrake as fd | |||
There was a problem hiding this comment.
This needs some more thought. It seems this module is just a copy of the original inversion_tools.py with modifications needed to accommodate vector values instead of scalars. This unfortunately not the way to develop code in a sustainable way, as we have now effectively doubled the amount of code that needs to be maintained. What happens if we further down the line want to further develop/fix bugs/make changes in the scalar inversion tools - or there are interface changes in the internal bits of Thetis that are being used in this module. Do I fix both modules, or do we let things diverge? I also can't easily see where this module is actually different from the original. As a general rule, "Copying code" is almost never the answer - except maybe sometimes for quickly getting something to work, but then if you want to keep it, you need to go back and figure out how to refactor and integrate the two without duplication
There was a problem hiding this comment.
Yes, agree with everything you've said. It is indeed just a copy with modifications for outputting the control and gradients as well as updating from elevation -> velocity.
I intentionally didn't merge this script with the base one because I wasn't sure how would be best to go about certain aspects:
- Do we assume that the calibration parameter is always elevation or velocity (probably) and then automatically detect whether we are dealing with a scalar or vector field? What if we want to use both elevation and velocity data for calibration? Do we limit it to one just now and leave this as a future problem?
- Is the way I am exporting things optimal? It's a bit of a pain because exporting is only set up for when the control parameter is the whole field, but the way I have dealt with it feels horrible. It seems like exporting might be nicer outside of the InversionManager and StationObservationManager classes...
I will work on this on the second batch of commits!
There was a problem hiding this comment.
Made the final updates to deal with this. It is assumed that the field for loss evaluation is either velocity or elevation (default). The exporting is fine in any case, though I'd still be tempted to update it in the future.
One extra thing I did was change the base DiagnosticCallback so that the simulation start time is an attribute when exporting to hdf5. This is easier/more direct to load than taking the units string of the variable data and something that a few people had asked me about.
- Add more introduction - Take the general inversion information out of the full nodal flexibility case - Add some references
- Small README update also - Adds necessary shapefiles for forward run
- So that we don't need geopandas as a dependency - The forward run is tested during the inversion run anyway - Also fix a lint issue in inversion script
- Update examples accordingly
- This is easier to directly load than the units 'time since initial time' - still keep the units
Add an inversion example where the velocity is used in the loss formulation, with different methods for constraining the Manning's field e.g. constant, region-based, independent points scheme