From 64b0549b8f09f7132279e3f8209f8f314bfe8b31 Mon Sep 17 00:00:00 2001 From: Mateusz Matuszak Date: Fri, 10 May 2024 11:52:40 +0000 Subject: [PATCH] Small changes to FTLE method --- opendrift/models/basemodel/__init__.py | 11 ++++++++--- opendrift/models/physics_methods.py | 27 ++++++++++++++++++++++++-- 2 files changed, 33 insertions(+), 5 deletions(-) diff --git a/opendrift/models/basemodel/__init__.py b/opendrift/models/basemodel/__init__.py index c526e7809..171e494d1 100644 --- a/opendrift/models/basemodel/__init__.py +++ b/opendrift/models/basemodel/__init__.py @@ -4661,7 +4661,10 @@ def calculate_ftle(self, o.run(duration=duration, time_step=time_step) f_x1, f_y1 = proj(o.history['lon'].T[-1].reshape(X.shape), o.history['lat'].T[-1].reshape(X.shape)) - lcs['RLCS'][i, :, :] = ftle(f_x1 - X, f_y1 - Y, delta, T) + #changes by Mateusz M + #Only final particle position, not difference between new and initial. + #lcs['RLCS'][i, :, :] = ftle(f_x1 - X, f_y1 - Y, delta, T) + lcs['RLCS'][i, :, :] = ftle(f_x1, f_y1, delta, T) # Backwards if ALCS is True: o = self.clone() @@ -4672,8 +4675,10 @@ def calculate_ftle(self, o.run(duration=duration, time_step=-time_step) b_x1, b_y1 = proj(o.history['lon'].T[-1][::-1].reshape(X.shape), o.history['lat'].T[-1][::-1].reshape(X.shape)) - lcs['ALCS'][i, :, :] = ftle(b_x1 - X, b_y1 - Y, delta, T) - + #changes by Mateusz M + #Only final particle position, not difference between new and initial. + #lcs['ALCS'][i, :, :] = ftle(b_x1 - X, b_y1 - Y, delta, T) + lcs['ALCS'][i, :, :] = ftle(b_x1, b_y1, delta, T) lcs['RLCS'] = np.ma.masked_invalid(lcs['RLCS']) lcs['ALCS'] = np.ma.masked_invalid(lcs['ALCS']) # Flipping ALCS left-right. Not sure why this is needed diff --git a/opendrift/models/physics_methods.py b/opendrift/models/physics_methods.py index 804debd49..bd831da51 100644 --- a/opendrift/models/physics_methods.py +++ b/opendrift/models/physics_methods.py @@ -515,7 +515,9 @@ def ftle(X, Y, delta, duration): # From Johannes Rohrs nx = X.shape[0] ny = X.shape[1] - J = np.empty([nx,ny,2,2], np.float32) + + #Changes by Mateusz M + """J = np.empty([nx,ny,2,2], np.float32) FTLE = np.empty([nx,ny], np.float32) # gradient @@ -527,7 +529,7 @@ def ftle(X, Y, delta, duration): J[:,:,1,0] = dy[0] / (2*delta) J[:,:,0,1] = dx[1] / (2*delta) J[:,:,1,1] = dy[1] / (2*delta) - + for i in range(0,nx): for j in range(0,ny): # Green-Cauchy tensor @@ -535,7 +537,28 @@ def ftle(X, Y, delta, duration): # its largest eigenvalue lamda = np.linalg.eigvals(D) FTLE[i,j] = np.log(np.sqrt(max(lamda)))/np.abs(duration) + """ + #from np.empty to np.zeros because new method doesn't work on boundaries. Empry left strange values. + FTLE = np.zeros([nx,ny]) + for i in range(1,nx-1): + for j in range(1,ny-1): + #Centering grid on particle and calculating Jacobian of flow map + J = np.zeros([2,2]) + J[0,0] = (X[i+1, j]-X[i-1, j])/(2*delta) + J[0,1] = (X[i, j+1]-X[i, j-1])/(2*delta) + J[1,0] = (Y[i+1, j]-Y[i-1, j])/(2*delta) + J[1,1] = (Y[i, j+1]-Y[i, j-1])/(2*delta) + # Green-Cauchy tensor + D = np.dot(np.transpose(J), J) + # Needed this, but don't remember the details + try: + lamda, v = np.linalg.eig(D) + except: + lamda = 0 + # its largest eigenvalue + FTLE[i,j] = np.log(np.sqrt(np.max(lamda)))/np.abs(duration) + return FTLE __stokes_coefficients__ = None