SpeedyWeather / SpeedyWeather.jl

Play atmospheric modelling like it's LEGO.
https://speedyweather.github.io/SpeedyWeather.jl/dev
MIT License
400 stars 24 forks source link

Lagrangian sampling of the model state #549

Open miniufo opened 3 weeks ago

miniufo commented 3 weeks ago

When I add Lagrangian particles to the simulation, I also need to sample the model state along each Lagrangian trajectory. For example, I need a time series of vorticity, or maybe a passive tracer (or active tracer like PV) released in the dynamic flow field, along a particle track. Just wondering if there is way to do this.

The tracking scheme requires to transform the fields from spectral space back to physical space, which may allow a straightforward along-track sampling of the flow field. If I want other fields like pressure or eta (in shallow water case), will this increase the computational cost as we need extra transformation of pressure or eta?

milankl commented 2 weeks ago

When I add Lagrangian particles to the simulation, I also need to sample the model state along each Lagrangian trajectory. For example, I need a time series of vorticity, or maybe a passive tracer (or active tracer like PV) released in the dynamic flow field, along a particle track. Just wondering if there is way to do this.

Yes, for every step of the particle advection, the velocities u, v have to be interpolated onto the particle locations and in the same way you can interpolate any other gridded field onto the same locations (otherwise you would have to update the locator inside the interpolator). So I would define a callback that has some preallocated array to store the record of whichever variable you want to track and then you essentially do like in src/dynamics/particle_advection.jl

https://github.com/SpeedyWeather/SpeedyWeather.jl/blob/4be86102c7a48c78b52edf3c6cc052df3fc34337/src/dynamics/particle_advection.jl#L182-L183

with u_new an array of length nparticles and u_grid any gridded field you want to interpolate. So inside your callback callback! you probably would have some lines that look like (using vor, i.e. vorticity as an example)

(; vor_tracked) = callback
(; vor_grid) = diagn.layers[1].grid_variables
(; interpolator) = diagn.particles

interpolate!(vor_tracked, vor_grid, interpolator)

The interpolator will already be set the the particles locations, but you have to make sure that the particle advection is executed on the same time steps as the tracking (or not, depending on what you want).

The tracking scheme requires to transform the fields from spectral space back to physical space, which may allow a straightforward along-track sampling of the flow field.

If you're thinking about hooking in new locations into the spectral transform ... no that's not a good idea. That would allow for some flexibility that'll likely bite back in performance for all other 99% of transforms you do. We don't support this. But that's not a problem: All variables that are available in spectral space are also available in grid space, so you just pick them up there. We have to do the transforms anyway.

If I want other fields like pressure or eta (in shallow water case), will this increase the computational cost as we need extra transformation of pressure or eta?

No for eta you'd use diagn.surface.pres_grid etc. They are all in grid space too.