Particle advection time stepping corrected

[milankl]

Time stepping for particles is now

• initialized (i.e. $u(t=0), v(t=0)$ interpolated to initial particle locations for the predictor step
• then executed at the end of a time step just before the n-th time step, where u,v are already gridded! at $t=n\Delta t$ (even though the clock time didn't step forward yet)
• ignores the initial euler and unfiltered leapfrog step (as before)

So that (PA Particle advection, otherwise dynamics)

julia> run!(simulation, period=Day(1));
PA: Initial velocities interpolated on initial locations[ Info: (0, DateTime("2000-01-01T00:00:00"))
[ Info: (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"))
PA: [ Info: (7, DateTime("2000-01-01T04:00:00"))
[ Info: (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"))
PA: [ Info: (15, DateTime("2000-01-01T08:00:00"))
[ Info: (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"))
PA: [ Info: (23, DateTime("2000-01-01T12:00:00"))
[ Info: (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"))
PA: [ Info: (31, DateTime("2000-01-01T16:00:00"))
[ Info: (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"))
PA: [ Info: (39, DateTime("2000-01-01T20:00:00"))
[ Info: (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"))
PA: [ Info: (47, DateTime("2000-01-02T00:00:00"))

Note that the notion of time somewhat different

• dynamics: the denoted time is time $t$ where the tendencies are evaluated to leapfrog from $t-\Delta t$ over $t$ to $t+\Delta t$ (hence the integration ends at t = 24h)
• particle advection: the denoted time $t$ is the time where the Heun step ends, so particle advection steps from $t-\Delta$ to $t$ with velocity at $t-\Delta t$ (the predictor step) averaged with the velocity at $t$ (the corrector step)

We may illustrate this as follows. Note two half steps in first_timsteps! to get from 0 to 0 (Euler forward, not counted) and from 0 to 1 (unfiltered leapfrog)

Dynamics (step i): | -- | -- | ------ | ------ | ------ | ------ | ------ | ------ | ------ | 0 0 1 2 3 4 5 6 7 8 Particle advection (step j, every 4th dynamics step) | -- | -- | ------ | ------ | ------ | ------ | ------ | ------ | ------ | 0->P1 C1->1->P2 C2->2

So that to get from particle step j=0 to j=1 it's the average betwee predictor step P1 which uses velocities at i=0 (stored during initialisation) and corrector step C1 which uses velocities at i=4; which is also directly used to interpolate and store velocities for predictor step P2 which is then averaged with corrector step C2 using velocities at i=8 to get from particle step j=1 to j=2 when the dynamics have reached i=8.

[milankl]

I hear you saying that you would use particle_advection = ParticleAdvection2D(spectral_grid, every_n_timesteps=1) i.e. advection on everyt time step, the last commit makes sure that also in that case particle advection is never executed on the 1st step in first_timesteps! (which is a half_step) to have 0, 0p, 1p, 2p, 3p, ... (dynamics 0 0 1 2 3, particle advection p) not 0p, 0p, 1p, 2p, 3p, ... which would cause an initial offset by dt/2.