milankl
commented
2 months ago I just printed the time step and time for dynamics and particle advection (PA)
julia> run!(simulation, period=Day(1));
[ Info: (0, DateTime("2000-01-01T00:00:00"))
[ Info: ("PA:", 0, DateTime("2000-01-01T00:00:00"))
[ Info: (0, DateTime("2000-01-01T00:15:00"))
[ Info: (1, DateTime("2000-01-01T00:30:00"))
[ Info: (2, DateTime("2000-01-01T01:00:00"))
[ Info: (3, DateTime("2000-01-01T01:30:00"))
[ Info: (4, DateTime("2000-01-01T02:00:00"))
[ Info: (5, DateTime("2000-01-01T02:30:00"))
[ Info: (6, DateTime("2000-01-01T03:00:00"))
[ Info: (7, DateTime("2000-01-01T03:30:00"))
[ Info: (8, DateTime("2000-01-01T04:00:00"))
[ Info: ("PA:", 8, DateTime("2000-01-01T04:00:00"))
[ Info: (9, DateTime("2000-01-01T04:30:00"))
[ Info: (10, DateTime("2000-01-01T05:00:00"))
[ Info: (11, DateTime("2000-01-01T05:30:00"))
[ Info: (12, DateTime("2000-01-01T06:00:00"))
[ Info: (13, DateTime("2000-01-01T06:30:00"))
[ Info: (14, DateTime("2000-01-01T07:00:00"))
[ Info: (15, DateTime("2000-01-01T07:30:00"))
[ Info: (16, DateTime("2000-01-01T08:00:00"))
[ Info: ("PA:", 16, DateTime("2000-01-01T08:00:00"))
[ Info: (17, DateTime("2000-01-01T08:30:00"))
[ Info: (18, DateTime("2000-01-01T09:00:00"))
[ Info: (19, DateTime("2000-01-01T09:30:00"))
[ Info: (20, DateTime("2000-01-01T10:00:00"))
[ Info: (21, DateTime("2000-01-01T10:30:00"))
[ Info: (22, DateTime("2000-01-01T11:00:00"))
[ Info: (23, DateTime("2000-01-01T11:30:00"))
[ Info: (24, DateTime("2000-01-01T12:00:00"))
[ Info: ("PA:", 24, DateTime("2000-01-01T12:00:00"))
[ Info: (25, DateTime("2000-01-01T12:30:00"))
[ Info: (26, DateTime("2000-01-01T13:00:00"))
[ Info: (27, DateTime("2000-01-01T13:30:00"))
[ Info: (28, DateTime("2000-01-01T14:00:00"))
[ Info: (29, DateTime("2000-01-01T14:30:00"))
[ Info: (30, DateTime("2000-01-01T15:00:00"))
[ Info: (31, DateTime("2000-01-01T15:30:00"))
[ Info: (32, DateTime("2000-01-01T16:00:00"))
[ Info: ("PA:", 32, DateTime("2000-01-01T16:00:00"))
[ Info: (33, DateTime("2000-01-01T16:30:00"))
[ Info: (34, DateTime("2000-01-01T17:00:00"))
[ Info: (35, DateTime("2000-01-01T17:30:00"))
[ Info: (36, DateTime("2000-01-01T18:00:00"))
[ Info: (37, DateTime("2000-01-01T18:30:00"))
[ Info: (38, DateTime("2000-01-01T19:00:00"))
[ Info: (39, DateTime("2000-01-01T19:30:00"))
[ Info: (40, DateTime("2000-01-01T20:00:00"))
[ Info: ("PA:", 40, DateTime("2000-01-01T20:00:00"))
[ Info: (41, DateTime("2000-01-01T20:30:00"))
[ Info: (42, DateTime("2000-01-01T21:00:00"))
[ Info: (43, DateTime("2000-01-01T21:30:00"))
[ Info: (44, DateTime("2000-01-01T22:00:00"))
[ Info: (45, DateTime("2000-01-01T22:30:00"))
[ Info: (46, DateTime("2000-01-01T23:00:00"))
[ Info: (47, DateTime("2000-01-01T23:30:00"))As you can see the first Euler step starts with 0, the "2nd" time step (the first Leapfrog step) isn't counted, so that after time step 1 30min, time step 2 1h, timestep 3 1h30min, etc. The particle advection is then done after each 8 (by default but that's a parameter you can choose) timesteps correctly skipping the not counted first time steps.
The 48th time step that would step from 24h to 24h30 is not executed as the simulation period is set to Day(1). So I am indeed wondering whether we should do the particle advection after time step -1, 7, 15, 23, 31, 39, 47, instead of after 0, 8, 16, 24, 32, 40, 48 (which is right now not executed). I also realised that we're technically using the velocities at n+1 (i.e. 1, 9, 17, 25, ...) because particle_advection! is currently called after gridded! which transforms the spectral state at the end of a time step (i.e. the new time step) to grid-point space. Maybe that should be changed too. The decision when to execute the particle advection is made here
and particle_advection! is called here
Feel free to put some more thought into when to call particle advection and createa pull request for it!
I've notice that the time integration of the model is the leap-frog scheme, which uses:
The scheme for particle tracking is predictor-corrector, which uses Δt/2 for prediction and then Δt for correction.
I am a bit worry about the syncronization of the time for both outputs. For example, after one day, does the output time of particle tracking is exactly the same as the one for the model? Is there a Δt / 2 or Δt offset? I am asking this because I knew that some model may have this problem. As a result, the output of particle position is not at the same time as the output of the dynamic variables. If I want to interpolate the dynamic variable onto particle position (e.g., PV, which is almost invariant along particle trajectory), this would yield a small error.