Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Small changes to FTLE method #1295

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Small changes to FTLE method #1295

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 8 additions & 3 deletions opendrift/models/basemodel/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand All @@ -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
Expand Down
27 changes: 25 additions & 2 deletions opendrift/models/physics_methods.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -527,15 +529,36 @@ 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
D = np.dot(np.transpose(J[i,j]), J[i,j])
# 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)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks generally good. But is there a particular reason why you have moved calculation of Jacobian from outside of the loop and into the loop? I fear that this will give much longer computation time, as the overhead of each calculation is now repeated typically up to order of a million times (~if 1000 x 1000).
Did you notice such a slowdown?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see your point and yes, I noticed that it became somewhat slower. There should be a way to take it out of the loop, but I don't see how right now.

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
Expand Down