Add files via upload

Advect_and_Calculate_2DSteady.m

@@ -152,9 +152,9 @@
 time_note = cat(2,time_note{:});
 TSE=cat(2,TSE_spmd_NM{:})/(tspan(end)-tspan(1));
 TSE_Bar=cat(2,TSE_spmd{:})/(tspan(end)-tspan(1));
-TRA_Bar=cat(2,TRA_spmd{:});
+TRA_Bar=cat(2,TRA_spmd{:})/(tspan(end)-tspan(1));
 
-TRA=real(acos(V1(:,1).*V2(:,1)+V1(:,2).*V2(:,2)));
+TRA=real(acos(V1(:,1).*V2(:,1)+V1(:,2).*V2(:,2)))/(tspan(end)-tspan(1));
 
 clear x_spmd y_spmd
 

Advect_and_Calculate_2DUnsteady.m

@@ -7,8 +7,8 @@
 % Output arguments:
 %   xt,yt    : xx-component, yy-component of trajectory final position
 %   time_note    : tspan time if a trajectory left interpolant domain
-%   TSE_Bar_EPS,TRA_Bar_EPS,TRA_EPS    : Single trajectory metrics from
-%   section IV in [1]
+%   TSE_Bar_EPS,TRA_Bar_EPS,TSE_EPS    : Single trajectory metrics from
+%   section IV in [1] (Theorem 3)
 
 %--------------------------------------------------------------------------
 % Author: Nikolas Aksamit  [email protected]
@@ -25,7 +25,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
-function     [xt,yt,time_note,TSE_Bar_EPS,TRA_Bar_EPS,TRA_EPS] = Advect_and_Calculate_2DUnsteady(tspan,xx,yy,U_Interp,V_Interp,NCores)
+function     [xt,yt,time_note,TSE_Bar_EPS,TRA_Bar_EPS,TSE_EPS] = Advect_and_Calculate_2DUnsteady(tspan,xx,yy,U_Interp,V_Interp,NCores)
 
 %%% Calculate v_0 as spatio-temporal mean of the flow in your domain of
 %%% interest
@@ -131,19 +131,14 @@
 
 V1 = cat(1,V1_spmd{:});
 V2 = cat(1,V2_spmd{:});
-V1_EPS=[V1/v_0, ones(size(V1,1),1)];
-V2_EPS=[V2/v_0, ones(size(V1,1),1)];
-
-V1_EPS=V1_EPS./(sqrt(sum(V1_EPS.*V1_EPS,2)));
-V2_EPS=V2_EPS./(sqrt(sum(V2_EPS.*V2_EPS,2)));
-TRA_EPS=real(acos(V1_EPS(:,1).*V2_EPS(:,1)+V1_EPS(:,2).*V2_EPS(:,2)));
+TSE_EPS=1/(tspan(end)-tspan(1))*log(sqrt((sum(V2.*V2,2)+v_0.^2)./(sum(V1.*V1,2).^2+v_0.^2)));
 
 xt = cat(2,x_spmd{:});
 yt = cat(2,y_spmd{:});
 time_note = cat(2,time_note{:});
 
 TSE_Bar_EPS=cat(2,TSE_EPS_spmd{:})/(tspan(end)-tspan(1));
-TRA_Bar_EPS=cat(2,TRA_EPS_spmd{:});
+TRA_Bar_EPS=cat(2,TRA_EPS_spmd{:})/(tspan(end)-tspan(1));
 
 clear x_spmd y_spmd
 

Advect_and_Calculate_3DSteady.m

@@ -54,9 +54,7 @@
     x_spmd=xt(1,Range);
     y_spmd=yt(1,Range);
     z_spmd=zt(1,Range);
-
-
-    Acc_spmd=0*z_spmd;
+
     time_note=zeros(size(x_spmd));
     TSE_spmd=zeros(size(x_spmd));
     TRA_spmd=zeros(size(x_spmd));
@@ -96,8 +94,7 @@
         y_spmd(2,:) = y_spmd(1,:) + deltay(2,:);
         z_spmd(2,:) = z_spmd(1,:) + deltaz(2,:);
 
-
-
+
         %%%% If particle leaves domain, keep record of last position.
         %%%% Record this time of leaving domain in time_note.
         x_spmd(1,~isnan(x_spmd(2,:)) & ~isnan(y_spmd(2,:)))=x_spmd(2,~isnan(x_spmd(2,:)) & ~isnan(y_spmd(2,:)));
@@ -150,8 +147,7 @@
         end
     end
 
-
-
+
     TSE_spmd=TSE_spmd(1,:);
     TRA_spmd=TRA_spmd(1,:);
     TSE_spmd_NM=TSE_spmd_NM(1,:);
@@ -168,9 +164,9 @@
 time_note = cat(2,time_note{:});
 TSE=cat(2,TSE_spmd_NM{:})/(tspan(end)-tspan(1));
 TSE_Bar=cat(2,TSE_spmd{:})/(tspan(end)-tspan(1));
-TRA_Bar=cat(2,TRA_spmd{:});
+TRA_Bar=cat(2,TRA_spmd{:})/(tspan(end)-tspan(1));
 
-TRA=real(acos(V1(:,1).*V2(:,1)+V1(:,2).*V2(:,2)));
+TRA=real(acos(V1(:,1).*V2(:,1)+V1(:,2).*V2(:,2)))/(tspan(end)-tspan(1));
 
 clear x_spmd y_spmd
 

Advect_and_Calculate_3DUnsteady.m

@@ -7,8 +7,8 @@
 % Output arguments:
 %   xt,yt,zt    : xx-component, yy-component, zz-component of trajectory final position
 %   time_note    : tspan time if a trajectory left interpolant domain
-%   TSE_Bar_EPS,TRA_Bar_EPS,TRA_EPS    : Single trajectory metrics from
-%   section IV in [1]
+%   TSE_Bar_EPS,TRA_Bar_EPS,TSE_EPS    : Single trajectory metrics from
+%   section IV in [1] (Theorem 3)
 
 %--------------------------------------------------------------------------
 % Author: Nikolas Aksamit  [email protected]
@@ -25,7 +25,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
-function     [xt,yt,zt,time_note,TSE_Bar_EPS,TRA_Bar_EPS,TRA_EPS] = Advect_and_Calculate_3DUnsteady(tspan,xx,yy,zz,U_Interp,V_Interp,W_Interp,NCores)
+function     [xt,yt,zt,time_note,TSE_Bar_EPS,TRA_Bar_EPS,TSE_EPS] = Advect_and_Calculate_3DUnsteady(tspan,xx,yy,zz,U_Interp,V_Interp,W_Interp,NCores)
 
 %%% Calculate v_0 as spatio-temporal mean of the flow in your domain of
 %%% interest
@@ -147,12 +147,7 @@
 
 V1 = cat(1,V1_spmd{:});
 V2 = cat(1,V2_spmd{:});
-V1_EPS=[V1/v_0, ones(size(V1,1),1)];
-V2_EPS=[V2/v_0, ones(size(V1,1),1)];
-
-V1_EPS=V1_EPS./(sqrt(sum(V1_EPS.*V1_EPS,2)));
-V2_EPS=V2_EPS./(sqrt(sum(V2_EPS.*V2_EPS,2)));
-TRA_EPS=real(acos(V1_EPS(:,1).*V2_EPS(:,1)+V1_EPS(:,2).*V2_EPS(:,2)+V1_EPS(:,3).*V2_EPS(:,3)));
+TSE_EPS=1/(tspan(end)-tspan(1))*log(sqrt((sum(V2.*V2,2)+v_0.^2)./(sum(V1.*V1,2).^2+v_0.^2)));
 
 xt = cat(2,x_spmd{:});
 yt = cat(2,y_spmd{:});

Demo_Steady3D.m

@@ -22,8 +22,8 @@
 
 %%
 disp('Running Advection')
-x_sparse = linspace(2.26,2.36,199);
-y_sparse = linspace(0.33,0.47,200);
+x_sparse = linspace(2.26,2.36,299);
+y_sparse = linspace(0.33,0.47,300);
 z_sparse = 2.55;
 
 [x_gr,y_gr,z_gr]=ndgrid(x_sparse,y_sparse,z_sparse);
@@ -46,7 +46,7 @@
 surf(x_gr,y_gr,TRA_Bar)
 view(0,90)
 shading interp
-title('$\overline{\mathrm{TRA}}_0^{10^{-3}}$','Interpreter','latex')
+title('$\overline{\mathrm{TRA}}_0^{10^{-3}}(\mathbf{x}_0)$','Interpreter','latex')
 axis tight
 colormap(ax(1),[gray_col(:,1)+col, gray_col(:,2), gray_col(:,3)])
 daspect([1 1 1])
@@ -60,7 +60,7 @@
 surf(x_gr,y_gr,TRA)
 view(0,90)
 shading interp
-title('$\mathrm{TRA}_0^{10^{-3}}$','Interpreter','latex')
+title('$\mathrm{TRA}_0^{10^{-3}}(\mathbf{x}_0)$','Interpreter','latex')
 axis tight
 colormap(ax(2),[gray_col(:,1)+col, gray_col(:,2), gray_col(:,3)])
 daspect([1 1 1])
@@ -75,7 +75,7 @@
 surf(x_gr,y_gr,TSE_Bar)
 view(0,90)
 shading interp
-title('$\overline{\mathrm{TSE}}_0^{10^{-3}}$','Interpreter','latex')
+title('$\overline{\mathrm{TSE}}_0^{10^{-3}}(\mathbf{x}_0)$','Interpreter','latex')
 axis tight
 colormap(ax(3),[gray_col(:,1), gray_col(:,2), gray_col(:,3)+col])
 daspect([1 1 1])
@@ -89,7 +89,7 @@
 surf(x_gr,y_gr,TSE)
 view(0,90)
 shading interp
-title('$\mathrm{TSE}_0^{10^{-3}}$','Interpreter','latex')
+title('$\mathrm{TSE}_0^{10^{-3}}(\mathbf{x}_0)$','Interpreter','latex')
 axis tight
 colormap(ax(4),[gray_col(:,1), gray_col(:,2), gray_col(:,3)+col])
 daspect([1 1 1])

Demo_Steady3D_Geometry.m

@@ -33,16 +33,16 @@
 
 %%
 disp('Running Advection')
-x_sparse = linspace(2.26,2.36,199);
-y_sparse = linspace(0.33,0.47,200);
+x_sparse = linspace(2.26,2.36,299);
+y_sparse = linspace(0.33,0.47,300);
 z_sparse = 2.55;
 
 [x_gr,y_gr,z_gr]=ndgrid(x_sparse,y_sparse,z_sparse);
 [x_mesh,y_mesh,z_mesh]=meshgrid(x_sparse,y_sparse,z_sparse);
 
 
 tfin=0.75;
-sVec=linspace(0,tfin,10000);
+sVec=linspace(0,tfin,15000);
 [xt,yt,zt,time_note,TSE_Bar,~,TRA_Bar] = Advect_and_Calculate_3DSteady(sVec,x_gr(:),y_gr(:),z_gr(:),u_interp,v_interp,w_interp,NCores);
 
 disp(' ')
@@ -57,7 +57,7 @@
 surf(x_gr,y_gr,TRA_Bar)
 view(0,90)
 shading interp
-title('$\overline{\mathrm{TRA}}_0^{0.75}$','Interpreter','latex')
+title('$\overline{\mathrm{NTRA}}_0^{0.75}(\mathbf{x}_0)$','Interpreter','latex')
 axis tight
 colormap(ax(1),[gray_col(:,1)+col, gray_col(:,2), gray_col(:,3)])
 daspect([1 1 1])
@@ -72,7 +72,7 @@
 surf(x_gr,y_gr,TSE_Bar)
 view(0,90)
 shading interp
-title('$\overline{\mathrm{TSE}}_0^{0.75}$','Interpreter','latex')
+title('$\overline{\mathrm{NTSE}}_0^{0.75}(\mathbf{x}_0)$','Interpreter','latex')
 axis tight
 colormap(ax(2),[gray_col(:,1), gray_col(:,2), gray_col(:,3)+col])
 daspect([1 1 1])

Demo_Unsteady2D.m

@@ -38,21 +38,21 @@
 tf = t0+days_transport;
 tspan=linspace(t0,tf,tsteps);
 disp(' ')
-[xt,yt,time_note,TSE_Bar_EPS,TRA_Bar_EPS,TRA_EPS] = Advect_and_Calculate_2DUnsteady(tspan,x_grid(:),y_grid(:),u_interp,v_interp,NCores);
+[xt,yt,time_note,TSE_Bar_EPS,TRA_Bar_EPS,TSE_EPS] = Advect_and_Calculate_2DUnsteady(tspan,x_grid(:),y_grid(:),u_interp,v_interp,NCores);
 
 TRA_Bar_EPS=real(reshape(TRA_Bar_EPS,size(x_grid)));
-TRA_EPS=real(reshape(TRA_EPS,size(x_grid)));
+TSE_EPS=real(reshape(TSE_EPS,size(x_grid)));
 TSE_Bar_EPS=reshape(TSE_Bar_EPS,size(x_grid));
 
 %% Plotting
 col=linspace(0.25,0,256)';
 gray_col=colormap(gray(256));
 figure
 ax(1)=subplot(1,3,1);
-surf(x_grid,y_grid,TRA_EPS)
+surf(x_grid,y_grid,TSE_EPS)
 view(0,90)
 shading interp
-title('$\mathrm{TRA}_0^{30}$','Interpreter','latex')
+title('$\mathrm{TSE}_0^{30}(\mathbf{x}_0,\mathbf{v}_0)$','Interpreter','latex')
 axis tight
 colormap(ax(1),[gray_col(:,1), gray_col(:,2)+col, gray_col(:,3)])
 daspect([1 1 1])
@@ -61,11 +61,11 @@
 ylabel('Lat','Interpreter','latex')
 colorbar
 
-ax(2)=subplot(1,3,2);
+ax(2)=subplot(1,3,3);
 surf(x_grid,y_grid,TRA_Bar_EPS)
 view(0,90)
 shading interp
-title('$\overline{\mathrm{TRA}}_0^{30}$','Interpreter','latex')
+title('$\overline{\mathrm{TRA}}_0^{30}(\mathbf{x}_0,\mathbf{v}_0)$','Interpreter','latex')
 axis tight
 colormap(ax(2),[gray_col(:,1)+col, gray_col(:,2), gray_col(:,3)])
 daspect([1 1 1])
@@ -75,12 +75,12 @@
 hold on
 colorbar
 
-ax(3)=subplot(1,3,3);
+ax(3)=subplot(1,3,2);
 TSE_Bar_EPS=reshape(real(TSE_Bar_EPS),size(x_grid));
 surf(x_grid,y_grid,TSE_Bar_EPS)
 view(0,90)
 shading interp
-title('$\overline{\mathrm{TSE}}_0^{30}$','Interpreter','latex')
+title('$\overline{\mathrm{TSE}}_0^{30}(\mathbf{x}_0,\mathbf{v}_0)$','Interpreter','latex')
 axis tight
 colormap(ax(3),[gray_col(:,1), gray_col(:,2), gray_col(:,3)+col])
 daspect([1 1 1])

ReadMe

@@ -0,0 +1,54 @@
+ --------------------------------------------------------------------------
+ Author: Nikolas Aksamit  [email protected]
+--------------------------------------------------------------------------
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%% References:
+ [1] Haller, G., Aksamit, N. O., & Bartos, A. P. E. (2021). Quasi-Objective Coherent Structure Diagnostics from Single Trajectories. 
+  Chaos, 31, 043131-1–17. https://doi.org/10.1063/5.0044151
+ [2] Aksamit, N.O., Haller, G. (2021). Objective Momentum Barriers in Wall Turbulence. 
+  In Review,  http://arxiv.org/abs/2106.07372
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%% Calculate TRA and TSE diagnostics for 2D and 3D, steady and unsteady flows. 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+The four “Advect_and_Calculate” scripts provide the tools to calculate TRA, and TSE values as described in [1] for both 2D and 3D, steady and unsteady flows. For the unsteady cases, the extended phase space versions are calculated which require a choice of v_0. The default choice is the spatio-temporal mean of the flow velocity magnitude in the domain of interest.
+
+Each Advect_and_Calculate script uses similar input formats:
+
+Advect_and_Calculate_2DSteady.m
+% Input arguments:
+     tspan                        : Discrete advection timesteps used for RK4 advection scheme. This must include all intermediate steps between initial and final times(e.g. tspan=linspace(t_initial,t_final,100). For the steady case, these are dummy times as the velocity field is autonomous.
+     xx,yy                        : Initial positions for advection. These must be formatted as Nx1 vectors
+     U_Interp,V_Interp            : Velocity field interpolants with inputs (xx,yy)
+     NCores                       : Number of Cores for parpool
+
+    % Output arguments:
+       xt,yt    : xx-component, yy-component of trajectory final position. 
+       time_note    : tspan time if a trajectory left interpolant domain
+       TSE_Bar,TSE,TRA_Bar,TRA    : Single trajectory metrics from   section II and III in [1]
+
+Advect_and_Calculate_2DUnsteady.m
+% Input arguments:
+     tspan                        : Discrete advection timesteps used for RK4 advection scheme. This must include all intermediate steps between initial and final times(e.g. tspan=linspace(t_initial,t_final,100). For the unsteady case, these are times used in the velocity field interplant.
+     xx,yy                        : Initial positions for advection. These must be formatted as Nx1 vectors
+     U_Interp,V_Interp            : Velocity field interpolants with inputs (time,xx,yy)
+     NCores                       : Number of Cores for parpool
+
+    % Output arguments:
+       xt,yt    : xx-component, yy-component of trajectory final position. 
+       time_note    : tspan time if a trajectory left interpolant domain
+       TSE_Bar,TSE,TRA_Bar    : Extended phase-space single trajectory metrics from section IV [1]
+
+Note: For unsteady cases you can modify v_0, or leave as the mean of values used to define flow field interpolants. Advect_and_Calculate_3DSteady.m and Advect_and_Calculate_3DUnsteady.m are of the same format, with an extra z-component.
+
+Demo_Unsteady2D.m
+Demonstration of unsteady TRA and TSE calculations for geostrophic ocean current data as in [1].
+
+Demo_Steady3D.m
+Demonstration of steady TRA and TSE calculations for Johns Hopkins Turbulence Database Re=1000 Channel flow as in [2].
+
+Demo_Steady3D_Geometry.m
+Demonstration of steady NTRA and NTSE (normalized) calculations for Johns Hopkins Turbulence Database Re=1000 Channel flow as in [2].
+