diff --git a/src/parcels/kernels/_advection.py b/src/parcels/kernels/_advection.py index 10e7ee676..05df14c39 100644 --- a/src/parcels/kernels/_advection.py +++ b/src/parcels/kernels/_advection.py @@ -17,76 +17,62 @@ ] -def _constrain_dt_to_within_time_interval(time_interval, time, dt): - """Helper function to make sure dt does not go outside time_interval. - - This is especially relevant for higher-order RK methods (RK2, RK4, RK45), - which require interpolations at time + dt. If time is at the edges of the - time_interval (typically the last integration step), such an operation would - lead to an OutofTimeError. - """ - if time_interval: - dt = np.where(time + dt <= time_interval.time_length_as_flt, dt, time_interval.time_length_as_flt - time) - dt = np.where(time + dt >= 0, dt, time) - return dt - - def AdvectionRK2(particles, fieldset): # pragma: no cover """Advection of particles using second-order Runge-Kutta integration.""" - dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.t, particles.dt) (u1, v1) = fieldset.UV[particles] - x1, y1 = (particles.x + u1 * 0.5 * dt, particles.y + v1 * 0.5 * dt) - (u2, v2) = fieldset.UV[particles.t + 0.5 * dt, particles.z, y1, x1, particles] - particles.dx += u2 * dt - particles.dy += v2 * dt + x1 = particles.x + u1 * 0.5 * particles.dt + y1 = particles.y + v1 * 0.5 * particles.dt + (u2, v2) = fieldset.UV[particles.t + 0.5 * particles.dt, particles.z, y1, x1, particles] + particles.dx += u2 * particles.dt + particles.dy += v2 * particles.dt def AdvectionRK2_3D(particles, fieldset): # pragma: no cover """Advection of particles using second-order Runge-Kutta integration including vertical velocity.""" - dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.t, particles.dt) (u1, v1, w1) = fieldset.UVW[particles] - x1 = particles.x + u1 * 0.5 * dt - y1 = particles.y + v1 * 0.5 * dt - z1 = particles.z + w1 * 0.5 * dt - (u2, v2, w2) = fieldset.UVW[particles.t + 0.5 * dt, z1, y1, x1, particles] - particles.dx += u2 * dt - particles.dy += v2 * dt - particles.dz += w2 * dt + x1 = particles.x + u1 * 0.5 * particles.dt + y1 = particles.y + v1 * 0.5 * particles.dt + z1 = particles.z + w1 * 0.5 * particles.dt + (u2, v2, w2) = fieldset.UVW[particles.t + 0.5 * particles.dt, z1, y1, x1, particles] + particles.dx += u2 * particles.dt + particles.dy += v2 * particles.dt + particles.dz += w2 * particles.dt def AdvectionRK4(particles, fieldset): # pragma: no cover """Advection of particles using fourth-order Runge-Kutta integration.""" - dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.t, particles.dt) (u1, v1) = fieldset.UV[particles] - x1, y1 = (particles.x + u1 * 0.5 * dt, particles.y + v1 * 0.5 * dt) - (u2, v2) = fieldset.UV[particles.t + 0.5 * dt, particles.z, y1, x1, particles] - x2, y2 = (particles.x + u2 * 0.5 * dt, particles.y + v2 * 0.5 * dt) - (u3, v3) = fieldset.UV[particles.t + 0.5 * dt, particles.z, y2, x2, particles] - x3, y3 = (particles.x + u3 * dt, particles.y + v3 * dt) - (u4, v4) = fieldset.UV[particles.t + dt, particles.z, y3, x3, particles] - particles.dx += (u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * dt - particles.dy += (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * dt + x1 = particles.x + u1 * 0.5 * particles.dt + y1 = particles.y + v1 * 0.5 * particles.dt + (u2, v2) = fieldset.UV[particles.t + 0.5 * particles.dt, particles.z, y1, x1, particles] + x2 = particles.x + u2 * 0.5 * particles.dt + y2 = particles.y + v2 * 0.5 * particles.dt + (u3, v3) = fieldset.UV[particles.t + 0.5 * particles.dt, particles.z, y2, x2, particles] + x3 = particles.x + u3 * particles.dt + y3 = particles.y + v3 * particles.dt + (u4, v4) = fieldset.UV[particles.t + particles.dt, particles.z, y3, x3, particles] + particles.dx += (u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * particles.dt + particles.dy += (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * particles.dt def AdvectionRK4_3D(particles, fieldset): # pragma: no cover """Advection of particles using fourth-order Runge-Kutta integration including vertical velocity.""" - dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.t, particles.dt) (u1, v1, w1) = fieldset.UVW[particles] - x1 = particles.x + u1 * 0.5 * dt - y1 = particles.y + v1 * 0.5 * dt - z1 = particles.z + w1 * 0.5 * dt - (u2, v2, w2) = fieldset.UVW[particles.t + 0.5 * dt, z1, y1, x1, particles] - x2 = particles.x + u2 * 0.5 * dt - y2 = particles.y + v2 * 0.5 * dt - z2 = particles.z + w2 * 0.5 * dt - (u3, v3, w3) = fieldset.UVW[particles.t + 0.5 * dt, z2, y2, x2, particles] - x3 = particles.x + u3 * dt - y3 = particles.y + v3 * dt - z3 = particles.z + w3 * dt - (u4, v4, w4) = fieldset.UVW[particles.t + dt, z3, y3, x3, particles] - particles.dx += (u1 + 2 * u2 + 2 * u3 + u4) / 6 * dt - particles.dy += (v1 + 2 * v2 + 2 * v3 + v4) / 6 * dt - particles.dz += (w1 + 2 * w2 + 2 * w3 + w4) / 6 * dt + x1 = particles.x + u1 * 0.5 * particles.dt + y1 = particles.y + v1 * 0.5 * particles.dt + z1 = particles.z + w1 * 0.5 * particles.dt + (u2, v2, w2) = fieldset.UVW[particles.t + 0.5 * particles.dt, z1, y1, x1, particles] + x2 = particles.x + u2 * 0.5 * particles.dt + y2 = particles.y + v2 * 0.5 * particles.dt + z2 = particles.z + w2 * 0.5 * particles.dt + (u3, v3, w3) = fieldset.UVW[particles.t + 0.5 * particles.dt, z2, y2, x2, particles] + x3 = particles.x + u3 * particles.dt + y3 = particles.y + v3 * particles.dt + z3 = particles.z + w3 * particles.dt + (u4, v4, w4) = fieldset.UVW[particles.t + particles.dt, z3, y3, x3, particles] + particles.dx += (u1 + 2 * u2 + 2 * u3 + u4) / 6 * particles.dt + particles.dy += (v1 + 2 * v2 + 2 * v3 + v4) / 6 * particles.dt + particles.dz += (w1 + 2 * w2 + 2 * w3 + w4) / 6 * particles.dt def AdvectionEE(particles, fieldset): # pragma: no cover @@ -105,8 +91,7 @@ def AdvectionRK45(particles, fieldset): # pragma: no cover Time-step dt is halved if error is larger than fieldset.RK45_tol, and doubled if error is smaller than 1/10th of tolerance. """ - dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.t, particles.dt) - sign_dt = np.sign(dt) + sign_dt = np.sign(particles.dt) c = [1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0] A = [ @@ -120,42 +105,37 @@ def AdvectionRK45(particles, fieldset): # pragma: no cover b5 = [16.0 / 135.0, 0.0, 6656.0 / 12825.0, 28561.0 / 56430.0, -9.0 / 50.0, 2.0 / 55.0] (u1, v1) = fieldset.UV[particles] - x1, y1 = (particles.x + u1 * A[0][0] * dt, particles.y + v1 * A[0][0] * dt) - (u2, v2) = fieldset.UV[particles.t + c[0] * dt, particles.z, y1, x1, particles] - x2, y2 = ( - particles.x + (u1 * A[1][0] + u2 * A[1][1]) * dt, - particles.y + (v1 * A[1][0] + v2 * A[1][1]) * dt, - ) - (u3, v3) = fieldset.UV[particles.t + c[1] * dt, particles.z, y2, x2, particles] - x3, y3 = ( - particles.x + (u1 * A[2][0] + u2 * A[2][1] + u3 * A[2][2]) * dt, - particles.y + (v1 * A[2][0] + v2 * A[2][1] + v3 * A[2][2]) * dt, - ) - (u4, v4) = fieldset.UV[particles.t + c[2] * dt, particles.z, y3, x3, particles] - x4, y4 = ( - particles.x + (u1 * A[3][0] + u2 * A[3][1] + u3 * A[3][2] + u4 * A[3][3]) * dt, - particles.y + (v1 * A[3][0] + v2 * A[3][1] + v3 * A[3][2] + v4 * A[3][3]) * dt, - ) - (u5, v5) = fieldset.UV[particles.t + c[3] * dt, particles.z, y4, x4, particles] - x5, y5 = ( - particles.x + (u1 * A[4][0] + u2 * A[4][1] + u3 * A[4][2] + u4 * A[4][3] + u5 * A[4][4]) * dt, - particles.y + (v1 * A[4][0] + v2 * A[4][1] + v3 * A[4][2] + v4 * A[4][3] + v5 * A[4][4]) * dt, - ) - (u6, v6) = fieldset.UV[particles.t + c[4] * dt, particles.z, y5, x5, particles] - - x_4th = (u1 * b4[0] + u2 * b4[1] + u3 * b4[2] + u4 * b4[3] + u5 * b4[4]) * dt - y_4th = (v1 * b4[0] + v2 * b4[1] + v3 * b4[2] + v4 * b4[3] + v5 * b4[4]) * dt - x_5th = (u1 * b5[0] + u2 * b5[1] + u3 * b5[2] + u4 * b5[3] + u5 * b5[4] + u6 * b5[5]) * dt - y_5th = (v1 * b5[0] + v2 * b5[1] + v3 * b5[2] + v4 * b5[3] + v5 * b5[4] + v6 * b5[5]) * dt + x1 = particles.x + u1 * A[0][0] * particles.dt + y1 = particles.y + v1 * A[0][0] * particles.dt + (u2, v2) = fieldset.UV[particles.t + c[0] * particles.dt, particles.z, y1, x1, particles] + x2 = particles.x + (u1 * A[1][0] + u2 * A[1][1]) * particles.dt + y2 = particles.y + (v1 * A[1][0] + v2 * A[1][1]) * particles.dt + (u3, v3) = fieldset.UV[particles.t + c[1] * particles.dt, particles.z, y2, x2, particles] + x3 = particles.x + (u1 * A[2][0] + u2 * A[2][1] + u3 * A[2][2]) * particles.dt + y3 = particles.y + (v1 * A[2][0] + v2 * A[2][1] + v3 * A[2][2]) * particles.dt + (u4, v4) = fieldset.UV[particles.t + c[2] * particles.dt, particles.z, y3, x3, particles] + x4 = particles.x + (u1 * A[3][0] + u2 * A[3][1] + u3 * A[3][2] + u4 * A[3][3]) * particles.dt + y4 = particles.y + (v1 * A[3][0] + v2 * A[3][1] + v3 * A[3][2] + v4 * A[3][3]) * particles.dt + (u5, v5) = fieldset.UV[particles.t + c[3] * particles.dt, particles.z, y4, x4, particles] + x5 = particles.x + (u1 * A[4][0] + u2 * A[4][1] + u3 * A[4][2] + u4 * A[4][3] + u5 * A[4][4]) * particles.dt + y5 = particles.y + (v1 * A[4][0] + v2 * A[4][1] + v3 * A[4][2] + v4 * A[4][3] + v5 * A[4][4]) * particles.dt + (u6, v6) = fieldset.UV[particles.t + c[4] * particles.dt, particles.z, y5, x5, particles] + + x_4th = (u1 * b4[0] + u2 * b4[1] + u3 * b4[2] + u4 * b4[3] + u5 * b4[4]) * particles.dt + y_4th = (v1 * b4[0] + v2 * b4[1] + v3 * b4[2] + v4 * b4[3] + v5 * b4[4]) * particles.dt + x_5th = (u1 * b5[0] + u2 * b5[1] + u3 * b5[2] + u4 * b5[3] + u5 * b5[4] + u6 * b5[5]) * particles.dt + y_5th = (v1 * b5[0] + v2 * b5[1] + v3 * b5[2] + v4 * b5[3] + v5 * b5[4] + v6 * b5[5]) * particles.dt kappa = np.sqrt(np.pow(x_5th - x_4th, 2) + np.pow(y_5th - y_4th, 2)) - good_particles = (kappa <= fieldset.RK45_tol) | (np.fabs(dt) <= np.fabs(fieldset.RK45_min_dt)) + good_particles = (kappa <= fieldset.RK45_tol) | (np.fabs(particles.dt) <= np.fabs(fieldset.RK45_min_dt)) particles.dx += np.where(good_particles, x_5th, 0) particles.dy += np.where(good_particles, y_5th, 0) increase_dt_particles = ( - good_particles & (kappa <= fieldset.RK45_tol / 10) & (np.fabs(dt * 2) <= np.fabs(fieldset.RK45_max_dt)) + good_particles + & (kappa <= fieldset.RK45_tol / 10) + & (np.fabs(particles.dt * 2) <= np.fabs(fieldset.RK45_max_dt)) ) particles.next_dt = np.where(increase_dt_particles, particles.dt * 2, particles.dt) particles.next_dt = np.where( diff --git a/src/parcels/kernels/_sigmagrids.py b/src/parcels/kernels/_sigmagrids.py index 4f321b71a..72ea991d7 100644 --- a/src/parcels/kernels/_sigmagrids.py +++ b/src/parcels/kernels/_sigmagrids.py @@ -2,8 +2,6 @@ import numpy as np -from parcels.kernels._advection import _constrain_dt_to_within_time_interval - def convert_z_to_sigma_croco(fieldset, t, z, y, x, particle): """Calculate local sigma level of the particles, by linearly interpolating the @@ -49,27 +47,26 @@ def AdvectionRK2_3D_CROCO(particles, fieldset): # pragma: no cover category=RuntimeWarning, ) # Needed because of linear sampling of W with sigma conversion - dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.t, particles.dt) sigma = particles.z / fieldset.h[particles.t, np.zeros_like(particles.z), particles.y, particles.x] sig = convert_z_to_sigma_croco(fieldset, particles.t, particles.z, particles.y, particles.x, particles) (u1, v1) = fieldset.UV[particles.t, sig, particles.y, particles.x, particles] w1 = fieldset.W[particles.t, sig, particles.y, particles.x, particles] w1 *= sigma / fieldset.h[particles.t, np.zeros_like(particles.z), particles.y, particles.x] - x1 = particles.x + u1 * 0.5 * dt - y1 = particles.y + v1 * 0.5 * dt - sig_dep1 = sigma + w1 * 0.5 * dt + x1 = particles.x + u1 * 0.5 * particles.dt + y1 = particles.y + v1 * 0.5 * particles.dt + sig_dep1 = sigma + w1 * 0.5 * particles.dt dep1 = sig_dep1 * fieldset.h[particles.t, np.zeros_like(particles.z), y1, x1] - sig1 = convert_z_to_sigma_croco(fieldset, particles.t + 0.5 * dt, dep1, y1, x1, particles) - (u2, v2) = fieldset.UV[particles.t + 0.5 * dt, sig1, y1, x1, particles] - w2 = fieldset.W[particles.t + 0.5 * dt, sig1, y1, x1, particles] - w2 *= sig_dep1 / fieldset.h[particles.t + 0.5 * dt, np.zeros_like(particles.z), y1, x1] - x2 = particles.x + u2 * 0.5 * dt - y2 = particles.y + v2 * 0.5 * dt - sig_dep2 = sigma + w2 * 0.5 * dt - dep2 = sig_dep2 * fieldset.h[particles.t + 0.5 * dt, np.zeros_like(particles.z), y2, x2] + sig1 = convert_z_to_sigma_croco(fieldset, particles.t + 0.5 * particles.dt, dep1, y1, x1, particles) + (u2, v2) = fieldset.UV[particles.t + 0.5 * particles.dt, sig1, y1, x1, particles] + w2 = fieldset.W[particles.t + 0.5 * particles.dt, sig1, y1, x1, particles] + w2 *= sig_dep1 / fieldset.h[particles.t + 0.5 * particles.dt, np.zeros_like(particles.z), y1, x1] + x2 = particles.x + u2 * 0.5 * particles.dt + y2 = particles.y + v2 * 0.5 * particles.dt + sig_dep2 = sigma + w2 * 0.5 * particles.dt + dep2 = sig_dep2 * fieldset.h[particles.t + 0.5 * particles.dt, np.zeros_like(particles.z), y2, x2] - particles.dx += u2 * dt - particles.dy += v2 * dt + particles.dx += u2 * particles.dt + particles.dy += v2 * particles.dt particles.dz += (dep1 - particles.z) + (dep2 - particles.z)