From c85255c5af0957e612d8930a838bf2fa999f2953 Mon Sep 17 00:00:00 2001 From: Michele Mastropietro Date: Fri, 22 Feb 2013 10:57:51 +0100 Subject: [PATCH] Added matlab code samples. All of these code samples currently are mis-identified in my repositories. I'm donating them to the cause. --- samples/Matlab/Check_plot.m | 12 ++ samples/Matlab/FTLEH.m | 149 ++++++++++++++++++++++++ samples/Matlab/FTLE_reg.m | 178 +++++++++++++++++++++++++++++ samples/Matlab/Integrate1.m | 40 +++++++ samples/Matlab/Integrate2.m | 60 ++++++++++ samples/Matlab/Lagr.m | 28 +++++ samples/Matlab/Lagrangian_points.m | 16 +++ samples/Matlab/Poincare.m | 18 +++ samples/Matlab/RK4.m | 24 ++++ samples/Matlab/distance.m | 13 +++ samples/Matlab/double_gyre.m | 49 ++++++++ samples/Matlab/example.m | 41 +++++++ samples/Matlab/gpu_RKF45_FILE.m | 156 +++++++++++++++++++++++++ samples/Matlab/test_rk_par.m | 19 +++ 14 files changed, 803 insertions(+) create mode 100644 samples/Matlab/Check_plot.m create mode 100644 samples/Matlab/FTLEH.m create mode 100644 samples/Matlab/FTLE_reg.m create mode 100644 samples/Matlab/Integrate1.m create mode 100644 samples/Matlab/Integrate2.m create mode 100644 samples/Matlab/Lagr.m create mode 100644 samples/Matlab/Lagrangian_points.m create mode 100644 samples/Matlab/Poincare.m create mode 100644 samples/Matlab/RK4.m create mode 100644 samples/Matlab/distance.m create mode 100644 samples/Matlab/double_gyre.m create mode 100644 samples/Matlab/example.m create mode 100644 samples/Matlab/gpu_RKF45_FILE.m create mode 100644 samples/Matlab/test_rk_par.m diff --git a/samples/Matlab/Check_plot.m b/samples/Matlab/Check_plot.m new file mode 100644 index 00000000..1f79e01b --- /dev/null +++ b/samples/Matlab/Check_plot.m @@ -0,0 +1,12 @@ +x_0=linspace(0,100,101); +vx_0=linspace(0,100,101); +z=zeros(101,101); +for i=1:101 + for j=1:101 + z(i,j)=x_0(i)*vx_0(j); + end +end + +figure +pcolor(x_0,vx_0,z) +shading flat \ No newline at end of file diff --git a/samples/Matlab/FTLEH.m b/samples/Matlab/FTLEH.m new file mode 100644 index 00000000..12bab16c --- /dev/null +++ b/samples/Matlab/FTLEH.m @@ -0,0 +1,149 @@ +tic +clear all +%% Choice of the mass parameter +mu=0.1; + +%% Computation of Lagrangian Points +[xl1,yl1,xl2,yl2,xl3,yl3,xl4,yl4,xl5,yl5] = Lagr(mu); + +%% Computation of initial total energy +E_L1=-Omega(xl1,yl1,mu); +E=E_L1+0.03715; % Offset as in figure 2.2 "LCS in the ER3BP" + +%% Initial conditions range +x_0_min=-0.8; +x_0_max=-0.2; + +vx_0_min=-2; +vx_0_max=2; + +y_0=0; + +% Elements for grid definition +n=200; + +% Dimensionless integrating time +T=2; + +% Grid initializing +[x_0,vx_0]=ndgrid(linspace(x_0_min,x_0_max,n),linspace(vx_0_min,vx_0_max,n)); +vy_0=sqrt(2*E+2*Omega(x_0,y_0,mu)-vx_0.^2); + +% Kinetic energy computation +E_cin=E+Omega(x_0,y_0,mu); + +%% Transforming into Hamiltonian variables +px_0=vx_0-y_0; +py_0=vy_0+x_0; + +% Inizializing +x_T=zeros(n,n); +y_T=zeros(n,n); +px_T=zeros(n,n); +py_T=zeros(n,n); +filtro=ones(n,n); +E_T=zeros(n,n); +a=zeros(n,n); % matrix of numbers of integration steps for each integration +np=0; % number of integrated points + +fprintf(' con n = %i\n',n) + +%% Energy tolerance setting +energy_tol=inf; + +%% Computation of the Jacobian of the system +options=odeset('Jacobian',@cr3bp_jac); + +%% Parallel integration of equations of motion +parfor i=1:n + for j=1:n + if E_cin(i,j)>0 && isreal(vy_0(i,j)) % Check for real velocity and positive Kinetic energy + [t,Y]=ode45(@fH,[0 T],[x_0(i,j); y_0; px_0(i,j); py_0(i,j)],options); + % Try to obtain the name of the solver for a following use +% sol=ode45(@f,[0 T],[x_0(i,j); y_0; vx_0(i,j); vy_0(i,j)],options); +% Y=sol.y'; +% solver=sol.solver; + a(i,j)=length(Y); + %Saving solutions + x_T(i,j)=Y(a(i,j),1); + px_T(i,j)=Y(a(i,j),3); + y_T(i,j)=Y(a(i,j),2); + py_T(i,j)=Y(a(i,j),4); + %Computation of final total energy and difference with + %initial one + E_T(i,j)=EnergyH(x_T(i,j),y_T(i,j),px_T(i,j),py_T(i,j),mu); + delta_E=abs(E_T(i,j)-E); + if delta_E > energy_tol; %Check of total energy conservation + fprintf(' Ouch! Wrong Integration: i,j=(%i,%i)\n E_T=%.2f \n delta_E=%.2f\n\n',i,j,E_T(i,j),delta_E); + filtro(i,j)=2; %Saving position of the point + end + np=np+1; + else + filtro(i,j)=0; % 1=interesting point; 0=non-sense point; 2= bad integration point + end + end +end + +t_integrazione=toc; +fprintf(' n = %i\n',n) +fprintf(' energy_tol = %.2f\n',energy_tol) +fprintf('total \t%i\n',n^2) +fprintf('nunber \t%i\n',np) +fprintf('time to integrate \t%.2f s\n',t_integr) + +%% Back to Lagrangian variables +vx_T=px_T+y_T; +vy_T=py_T-x_T; +%% FTLE Computation +fprintf('adesso calcolo ftle\n') +tic +dphi=zeros(2,2); +ftle=zeros(n-2,n-2); + +for i=2:n-1 + for j=2:n-1 + if filtro(i,j) && ... % Check for interesting point + filtro(i,j-1) && ... + filtro(i,j+1) && ... + filtro(i-1,j) && ... + filtro(i+1,j) + + dphi(1,1)=(x_T(i-1,j)-x_T(i+1,j))/(x_0(i-1,j)-x_0(i+1,j)); + + dphi(1,2)=(x_T(i,j-1)-x_T(i,j+1))/(vx_0(i,j-1)-vx_0(i,j+1)); + + dphi(2,1)=(vx_T(i-1,j)-vx_T(i+1,j))/(x_0(i-1,j)-x_0(i+1,j)); + + dphi(2,2)=(vx_T(i,j-1)-vx_T(i,j+1))/(vx_0(i,j-1)-vx_0(i,j+1)); + + if filtro(i,j)==2 % Manual setting to visualize bad integrated points + ftle(i-1,j-1)=-Inf; + else + ftle(i-1,j-1)=1/(2*T)*log(max(abs(eig(dphi'*dphi)))); + end + end + end +end + +%% Plotting results +% figure +% plot(t,Y) +% figure +% plot(Y(:,1),Y(:,2)) +% figure + +xx=linspace(x_0_min,x_0_max,n); +vvx=linspace(vx_0_min,vx_0_max,n); +[x,vx]=ndgrid(xx(2:n-1),vvx(2:n-1)); +figure +pcolor(x,vx,ftle) +shading flat + +t_ftle=toc; +fprintf('tempo per integrare \t%.2f s\n',t_integrazione) +fprintf('tempo per calcolare ftle \t%.2f s\n',t_ftle) + +% save(['var_' num2str(n) '_' num2str(clock(4)]) + +nome=['var_xvx_', 'ode00', '_n',num2str(n),'_e',num2str(energy_tol),'_H']; +save(nome) \ No newline at end of file diff --git a/samples/Matlab/FTLE_reg.m b/samples/Matlab/FTLE_reg.m new file mode 100644 index 00000000..37376be2 --- /dev/null +++ b/samples/Matlab/FTLE_reg.m @@ -0,0 +1,178 @@ +tic +clear all +%% Elements for grid definition +n=100; + +%% Dimensionless integrating time +T=2; + +%% Choice of the mass parameter +mu=0.1; + +%% Computation of Lagrangian Points +[xl1,yl1,xl2,yl2,xl3,yl3,xl4,yl4,xl5,yl5] = Lagr(mu); + +%% Computation of initial total energy +E_L1=-Omega(xl1,yl1,mu); +C_L1=-2*E_L1; % C_L1 = 3.6869532299 from Szebehely +E=E_L1+0.03715; % Offset as in figure 2.2 "LCS in the ER3BP" + +%% Initial conditions range +x_0_min=-0.8; +x_0_max=-0.2; + +vx_0_min=-2; +vx_0_max=2; + +y_0=0; + +% Grid initializing +[x_0,vx_0]=ndgrid(linspace(x_0_min,x_0_max,n),linspace(vx_0_min,vx_0_max,n)); +vy_0=sqrt(2*E+2.*Omega(x_0,y_0,mu)-vx_0.^2); +% Kinetic energy computation +E_cin=E+Omega(x_0,y_0,mu); + +% Inizializing +x_T=zeros(n,n); +y_T=zeros(n,n); +vx_T=zeros(n,n); +vy_T=zeros(n,n); +filtro=ones(n,n); +E_T=zeros(n,n); +delta_E=zeros(n,n); +a=zeros(n,n); % matrix of numbers of integration steps for each integration +np=0; % number of integrated points + +fprintf('integro con n = %i\n',n) + +%% Energy tolerance setting +energy_tol=0.1; + +%% Setting the options for the integrator +RelTol=1e-12;AbsTol=1e-12; % From Short +% RelTol=1e-13;AbsTol=1e-22; % From JD James Mireles +% RelTol=3e-14;AbsTol=1e-16; % HIGH accuracy from Ross +options=odeset('AbsTol',AbsTol,'RelTol',RelTol); +%% Parallel integration of equations of motion +h=waitbar(0,'','Name','Integration in progress, please wait!'); +S=zeros(n,n); +r1=zeros(n,n); +r2=zeros(n,n); +g=zeros(n,n); +for i=1:n + waitbar(i/n,h,sprintf('Computing i=%i',i)); + parfor j=1:n + r1(i,j)=sqrt((x_0(i,j)+mu).^2+y_0.^2); + r2(i,j)=sqrt((x_0(i,j)-1+mu).^2+y_0.^2); + g(i,j)=((1-mu)./(r1(i,j).^3)+mu./(r2(i,j).^3)); + if E_cin(i,j)>0 && isreal(vy_0(i,j)) % Check for real velocity and positive Kinetic energy + S(i,j)=g(i,j)*T; + [s,Y]=ode45(@f_reg,[0 S(i,j)],[x_0(i,j); y_0; vx_0(i,j); vy_0(i,j)],options,mu); + a(i,j)=length(Y); +% if s(a(i,j)) < 2 +% filtro(i,j)=3; +% end + % Saving solutions + x_T(i,j)=Y(a(i,j),1); + vx_T(i,j)=Y(a(i,j),3); + y_T(i,j)=Y(a(i,j),2); + vy_T(i,j)=Y(a(i,j),4); + + % Computation of final total energy and difference with + % initial one + E_T(i,j)=Energy(x_T(i,j),y_T(i,j),vx_T(i,j),vy_T(i,j),mu); + delta_E(i,j)=abs(E_T(i,j)-E); + if delta_E(i,j) > energy_tol; % Check of total energy conservation + fprintf(' Ouch! Wrong Integration: i,j=(%i,%i)\n E_T=%.2f \n delta_E=%f\n\n',i,j,E_T(i,j),delta_E(i,j)); + filtro(i,j)=2; % Saving position of the point + end + np=np+1; + else + filtro(i,j)=0; % 1 = interesting point; 0 = non-sense point; 2 = bad integration point + end + end +end +close(h); +t_integrazione=toc; +%% +filtro_1=filtro; +for i=2:n-1 + for j=2:n-1 + if filtro(i,j)==2 || filtro (i,j)==3 + filtro_1(i,j)=2; + filtro_1(i+1,j)=2; + filtro_1(i-1,j)=2; + filtro_1(i,j+1)=2; + filtro_1(i,j-1)=2; + end + end +end + +fprintf('integato con n = %i\n',n) +fprintf('integato con energy_tol = %f\n',energy_tol) +fprintf('numero punti totali \t%i\n',n^2) +fprintf('numero punti integrati \t%i\n',np) +fprintf('tempo per integrare \t%.2f s\n',t_integrazione) + +%% FTLE Computation +fprintf('adesso calcolo ftle\n') +tic +dphi=zeros(2,2); +ftle=zeros(n-2,n-2); +ftle_norm=zeros(n-2,n-2); + +ds_x=(x_0_max-x_0_min)/(n-1); +ds_vx=(vx_0_max-vx_0_min)/(n-1); + +for i=2:n-1 + for j=2:n-1 + if filtro_1(i,j) && ... % Check for interesting point + filtro_1(i,j-1) && ... + filtro_1(i,j+1) && ... + filtro_1(i-1,j) && ... + filtro_1(i+1,j) + % La direzione dello spostamento la decide il denominatore + + % TODO spiegarsi teoricamente come mai la matrice pu� + % essere ridotta a 2x2 + dphi(1,1)=(x_T(i+1,j)-x_T(i-1,j))/(2*ds_x); %(x_0(i-1,j)-x_0(i+1,j)); + + dphi(1,2)=(x_T(i,j+1)-x_T(i,j-1))/(2*ds_vx); %(vx_0(i,j-1)-vx_0(i,j+1)); + + dphi(2,1)=(vx_T(i+1,j)-vx_T(i-1,j))/(2*ds_x); %(x_0(i-1,j)-x_0(i+1,j)); + + dphi(2,2)=(vx_T(i,j+1)-vx_T(i,j-1))/(2*ds_vx); %(vx_0(i,j-1)-vx_0(i,j+1)); + + if filtro_1(i,j)==2 % Manual setting to visualize bad integrated points + ftle(i-1,j-1)=0; + else + ftle(i-1,j-1)=(1/abs(T))*log(max(sqrt(abs(eig(dphi*dphi'))))); + ftle_norm(i-1,j-1)=(1/abs(T))*log(norm(dphi)); + end + end + end +end + +%% Plotting results +% figure +% plot(t,Y) +% figure +% plot(Y(:,1),Y(:,2)) +% figure + +xx=linspace(x_0_min,x_0_max,n); +vvx=linspace(vx_0_min,vx_0_max,n); +[x,vx]=ndgrid(xx(2:n-1),vvx(2:n-1)); +figure +pcolor(x,vx,ftle) +shading flat + +t_ftle=toc; +fprintf('tempo per integrare \t%.2f s\n',t_integrazione) +fprintf('tempo per calcolare ftle \t%.2f s\n',t_ftle) + +% ora=fstringf %TODO +% save(['var_' num2str(n) '_' num2str(clock(4)]) + +nome=['var_xvx_', 'ode00', '_n',num2str(n)]; +save(nome) \ No newline at end of file diff --git a/samples/Matlab/Integrate1.m b/samples/Matlab/Integrate1.m new file mode 100644 index 00000000..8f894cda --- /dev/null +++ b/samples/Matlab/Integrate1.m @@ -0,0 +1,40 @@ +function [ x_T, y_T, vx_T, e_T, filter ] = Integrate_FILE( x_0, y_0, vx_0, e_0, T, N, mu, options) +%Integrate +% This function performs Runge-Kutta-Fehlberg integration for given +% initial conditions to compute FILE +nx=length(x_0); +ny=length(y_0); +nvx=length(vx_0); +ne=length(e_0); +vy_0=zeros(nx,ny,nvx,ne); +x_T=zeros(nx,ny,nvx,ne); +y_T=zeros(nx,ny,nvx,ne); +vx_T=zeros(nx,ny,nvx,ne); +vy_T=zeros(nx,ny,nvx,ne); +e_T=zeros(nx,ny,nvx,ne); +%% Look for phisically meaningful points +filter=zeros(nx,ny,nvx,ne); %0=meaningless point 1=meaningful point + +%% Integrate only meaningful points +h=waitbar(0,'','Name','Integration in progress, please wait!'); +for i=1:nx + waitbar(i/nx,h,sprintf('Computing i=%i',i)); + for j=1:ny + parfor k=1:nvx + for l=1:ne + vy_0(i,j,k,l)=sqrt(2*Potential(x_0(i),y_0(j),mu)+2*e_0(l)-vx_0(k)^2); + if isreal(vy_0(i,j,k,l)) + filter(i,j,k,l)=1; + ci=[x_0(i), y_0(j), vx_0(k), vy_0(i,j,k,l)]; + [t,Y,te,ye,ie]=ode45(@f,[0 T], ci, options, mu); + x_T(i,j,k,l)=ye(N+1,1); + y_T(i,j,k,l)=ye(N+1,2); + vx_T(i,j,k,l)=ye(N+1,3); + vy_T(i,j,k,l)=ye(N+1,4); + e_T(i,j,k,l)=0.5*(vx_T(i,j,k,l)^2+vy_T(i,j,k,l)^2)-Potential(x_T(i,j,k,l),y_T(i,j,k,l),mu); + end + end + end + end +end +close(h); diff --git a/samples/Matlab/Integrate2.m b/samples/Matlab/Integrate2.m new file mode 100644 index 00000000..10e8bea3 --- /dev/null +++ b/samples/Matlab/Integrate2.m @@ -0,0 +1,60 @@ +function [ x_T, y_T, vx_T, e_T, filter, delta_e ] = Integrate_FTLE_Gawlick_ell( x_0, y_0, vx_0, e_0, T, mu, ecc, nu, options) +%Integrate +% This function performs Runge-Kutta-Fehlberg integration for given +% initial conditions to compute FTLE to obtain the image in the Gawlick's +% article "Lagrangian Coherent Structures in the Elliptic Restricted +% Three-Body Problem". +nx=length(x_0); +ny=length(y_0); +nvx=length(vx_0); +ne=length(e_0); +vy_0=zeros(nx,ny,nvx,ne); +x_T=zeros(nx,ny,nvx,ne); +y_T=zeros(nx,ny,nvx,ne); +vx_T=zeros(nx,ny,nvx,ne); +vy_T=zeros(nx,ny,nvx,ne); +e_T=zeros(nx,ny,nvx,ne); +delta_e=zeros(nx,ny,nvx,ne); +%% Look for phisically meaningful points +filter=zeros(nx,ny,nvx,ne); %0=meaningless point 1=meaningful point +useful=ones(nx,ny,nvx,ne); +%% Integrate only useful points +useful(:,1,:,1)=0; +useful(:,1,:,3)=0; +useful(:,3,:,1)=0; +useful(:,3,:,3)=0; + +%% Integrate only meaningful points +h=waitbar(0,'','Name','Integration in progress, please wait!'); +for i=1:nx + waitbar(i/nx,h,sprintf('Computing i=%i',i)); + for j=1:ny + parfor k=1:nvx + for l=1:ne + if useful(i,j,k,l) + vy_0(i,j,k,l)=-sqrt(2*(Omega(x_0(i),y_0(j),mu)/(1+ecc*cos(nu)))+2*e_0(l)-vx_0(k)^2); + if isreal(vy_0(i,j,k,l)) + filter(i,j,k,l)=1; + + ci=[x_0(i), y_0(j), vx_0(k), vy_0(i,j,k,l)]; + [t,Y]=ode45(@f_ell,[0 T], ci, options, mu, ecc); + + if abs(t(end)) < abs(T) % Consider also negative time + filter(i,j,k,l)=3 + end + + x_T(i,j,k,l)=Y(end,1); + y_T(i,j,k,l)=Y(end,2); + vx_T(i,j,k,l)=Y(end,3); + vy_T(i,j,k,l)=Y(end,4); + e_T(i,j,k,l)=0.5*(vx_T(i,j,k,l)^2+vy_T(i,j,k,l)^2)-Omega(x_T(i,j,k,l),y_T(i,j,k,l),mu); + + % Compute the goodness of the integration + delta_e(i,j,k,l)=abs(e_T(i,j,k,l)-e_0(l)); + end + end + end + end + end +end +close(h); \ No newline at end of file diff --git a/samples/Matlab/Lagr.m b/samples/Matlab/Lagr.m new file mode 100644 index 00000000..5c4b2387 --- /dev/null +++ b/samples/Matlab/Lagr.m @@ -0,0 +1,28 @@ +function [xl1,yl1,xl2,yl2,xl3,yl3,xl4,yl4,xl5,yl5] = Lagr(mu) +% [xl1,yl1,xl2,yl2,xl3,yl3,xl4,yl4,xl5,yl5] = Lagr(mu) +% Lagr This function computes the coordinates of the Lagrangian points, +% given the mass parameter +yl1=0; +yl2=0; +yl3=0; +yl4=sqrt(3)/2; +yl5=-sqrt(3)/2; +c1=roots([1 mu-3 3-2*mu -mu 2*mu -mu]); +c2=roots([1 3-mu 3-2*mu -mu -2*mu -mu]); +c3=roots([1 2+mu 1+2*mu mu-1 2*mu-2 mu-1]); +xl1=0; +xl2=0; +for i=1:5 + if isreal(c1(i)) + xl1=1-mu-c1(i); + end + if isreal(c2(i)) + xl2=1-mu+c2(i); + end + if isreal(c3(i)) + xl3=-mu-c3(i); + end +end +xl4=0.5-mu; +xl5=xl4; +end \ No newline at end of file diff --git a/samples/Matlab/Lagrangian_points.m b/samples/Matlab/Lagrangian_points.m new file mode 100644 index 00000000..0e71e71b --- /dev/null +++ b/samples/Matlab/Lagrangian_points.m @@ -0,0 +1,16 @@ +% Plot dei Lagrangian points +n=5; +mu=linspace(0,0.5,n); +for i=1:n + [xl1,yl1,xl2,yl2,xl3,yl3,xl4,yl4,xl5,yl5] = Lagr(mu(i)); + figure (1) + hold all + plot(xl1, yl1, 's') + plot(xl2, yl2, 's') + plot(xl3, yl3, 's') + plot(xl4, yl4, 's') + plot(xl5, yl5, 's') + plot(-mu,0,'o') + plot(1-mu,0, 'o') + plot([-mu(i) xl4],[0 yl4]) +end \ No newline at end of file diff --git a/samples/Matlab/Poincare.m b/samples/Matlab/Poincare.m new file mode 100644 index 00000000..a9fb41f8 --- /dev/null +++ b/samples/Matlab/Poincare.m @@ -0,0 +1,18 @@ +clear +%% Initial Conditions +mu=0.012277471; +T=10; +N=5; +C=3.17; +x_0=0.30910452642073; +y_0=0.07738174525518; +vx_0=-0.72560796964234; +vy_0=sqrt(-C-vx_0^2+2*Potential(x_0,y_0,mu)); +k=0; +%% Integration +options=odeset('AbsTol',1e-22,'RelTol',1e-13,'Events',@cross_y); +[t,y,te,ye,ie]=ode113(@f,[0 T],[x_0; y_0; vx_0; vy_0],options,mu); + +figure +%plot(ye(:,1),ye(:,3),'rs') +plot(ye(:,1),0,'rs') \ No newline at end of file diff --git a/samples/Matlab/RK4.m b/samples/Matlab/RK4.m new file mode 100644 index 00000000..d21ea755 --- /dev/null +++ b/samples/Matlab/RK4.m @@ -0,0 +1,24 @@ +function x = RK4( fun, tspan, ci, mu ) +%RK4 4th-order Runge Kutta integrator +% Detailed explanation goes here +h=1e-5; +t=tspan(1); +T=tspan(length(tspan)); +dim=length(ci); +%x=zeros(l,dim); +x(:,1)=ci; +i=1; +while t value=0 +D=sqrt((y(1)+mu).^2+y(2).^2); % distance from the largest primary + +value=d-D; +isterminal=1; +direction=0; +end \ No newline at end of file diff --git a/samples/Matlab/double_gyre.m b/samples/Matlab/double_gyre.m new file mode 100644 index 00000000..b0fc82e0 --- /dev/null +++ b/samples/Matlab/double_gyre.m @@ -0,0 +1,49 @@ +clear all +tic +% initialize integration time T, f(x,t), discretization size n ---------------- +T = 8; +x_min=0; +x_max=2; +y_min=0; +y_max=1; +n=50; % how many points per one measure unit (both in x and in y) +ds=1/(n-1); +x_res=(x_max-x_min)*n; +y_res=(y_max-y_min)*n; +grid_x=linspace(x_min,x_max,x_res); +grid_y=linspace(y_min,y_max,y_res); + +advected_x=zeros(x_res,y_res); +advected_y=zeros(x_res,y_res); +% integrate all initial points for t in [0,T] -------------------------------- +parfor i = 1:x_res + for j = 1:y_res + [t,X] = ode45(@dg,[0,T],[grid_x(i),grid_y(j)]); + % store advected positions as they would appear in (x,y) coords ------ + advected_x(i,j) = X(length(X(:,1)),1); + advected_y(i,j) = X(length(X(:,2)),2); + end +end +%% Compute FTLE +sigma=zeros(x_res,y_res); +% at each point in interior of grid, store FTLE ------------------------------ +for i = 2:x_res-1 + for j = 2:y_res-1 + % compute Jacobian phi ----------------------------------------------- + phi(1,1) = (advected_x(i+1,j)-advected_x(i-1,j))/(2*ds); + phi(1,2) = (advected_x(i,j-1)-advected_x(i,j+1))/(2*ds); + phi(2,1) = (advected_y(i+1,j)-advected_y(i-1,j))/(2*ds); + phi(2,2) = (advected_y(i,j-1)-advected_y(i,j+1))/(2*ds); + % find max eigenvalue of phi'*phi ------------------------------------ + lambda_max = max(abs(eig(phi'*phi))); + % store FTLE --------------------------------------------------------- + sigma(i,j) = log(lambda_max)/abs(2*T); + end +end +toc +%% plot FTLE field ------------------------------------------------------------ +figure +contourf(grid_x,grid_y,sigma'); +colorbar('location','EastOutside'); +axis equal +shading flat diff --git a/samples/Matlab/example.m b/samples/Matlab/example.m new file mode 100644 index 00000000..5547bdeb --- /dev/null +++ b/samples/Matlab/example.m @@ -0,0 +1,41 @@ +clear all +tic +% initialize integration time T, f(x,t), discretization size n ---------------- +T = 20; +f_x_t = inline('[v(2);-sin(v(1))]','t','v'); +grid_min = -3.4; +grid_max = 3.4; +grid_width = grid_max-grid_min; +n = 35; +grid_spacing = grid_min:(grid_width/(n-1)):grid_max; +advected_x=zeros(n,n); +advected_y=zeros(n,n); +% integrate all initial points for t in [0,T] -------------------------------- +for i = 1:n + for j = 1:n + [t,x] = ode45(f_x_t,[0,T],[grid_spacing(i),grid_spacing(j)]); + % store advected positions as they would appear in (x,y) coords ------ + advected_x(n-j+1,i) = x(length(x(:,1)),1); + advected_y(n-j+1,i) = x(length(x(:,2)),2); + end +end +sigma=zeros(n,n); +% at each point in interior of grid, store FTLE ------------------------------ +for i = 2:n-1 + for j = 2:n-1 + % compute Jacobian phi ----------------------------------------------- + phi(1,1) = (advected_x(i,j+1)-advected_x(i,j-1))/(2*grid_width/(n-1)); + phi(1,2) = (advected_x(i-1,j)-advected_x(i+1,j))/(2*grid_width/(n-1)); + phi(2,1) = (advected_y(i,j+1)-advected_y(i,j-1))/(2*grid_width/(n-1)); + phi(2,2) = (advected_y(i-1,j)-advected_y(i+1,j))/(2*grid_width/(n-1)); + % find max eigenvalue of phi'*phi ------------------------------------ + lambda_max = max(abs(eig(phi'*phi))); + % store FTLE --------------------------------------------------------- + sigma(i,j) = log(lambda_max)/abs(T); + end +end +toc +%% plot FTLE field ------------------------------------------------------------ +figure +contourf(grid_spacing,grid_spacing,sigma); +colorbar('location','EastOutside'); \ No newline at end of file diff --git a/samples/Matlab/gpu_RKF45_FILE.m b/samples/Matlab/gpu_RKF45_FILE.m new file mode 100644 index 00000000..89280d1a --- /dev/null +++ b/samples/Matlab/gpu_RKF45_FILE.m @@ -0,0 +1,156 @@ +tic +clear +%% Range definition +n=200; + +mu=0.1; +[xl1,yl1,xl2,yl2,xl3,yl3,xl4,yl4,xl5,yl5]=Lagr(mu); +C_L1=2*Omega(xl1,yl1,mu); +E_0=-C_L1/2+0.03715; +Y_0=0; + +nx=n; +x_0_min=-0.8; +x_0_max=-0.15; +x_0=linspace(x_0_min, x_0_max, nx); +dx=(x_0_max-x_0_min)/(nx-1); + +nvx=n; +vx_0_min=-2; +vx_0_max=2; +vx_0=linspace(vx_0_min, vx_0_max, nvx); +dvx=(vx_0_max-vx_0_min)/(nvx-1); + +ny=3; +dy=(dx+dvx)/2; +y_0=[Y_0-dy Y_0 Y_0+dy]; + + + +ne=3; +de=dy; +e_0=[E_0-de E_0 E_0+de]; + +%% Definition of arrays of initial conditions + +%In this approach, only useful pints are stored and integrated + +m=1; +% x=zeros(1,nx*ny*nvx*ne); +% y=zeros(1,nx*ny*nvx*ne); +% vx=zeros(1,nx*ny*nvx*ne); +% e=zeros(1,nx*ny*nvx*ne); +% vy=zeros(1,nx*ny*nvx*ne); +filter=zeros(nx,3,nvx,3); + +for i=1:nx + for j=1:ny + for k=1:nvx + for l=1:ne + v_y=-sqrt(2*Omega(x_0(i),y_0(j),mu)+2*e_0(l)-vx_0(k)^2); + if ~((j~=2) && (l~=2)) && isreal(v_y) + x(m)=x_0(i); + y(m)=y_0(j); + vx(m)=vx_0(k); + e(m)=e_0(l); + vy(m)=v_y; + filter(i,j,k,l)=1; + m=m+1; + end + end + end + end +end + +%% Selection of useful points + +%% Data transfer to GPU +x_gpu=gpuArray(x); +y_gpu=gpuArray(y); +vx_gpu=gpuArray(vx); +vy_gpu=gpuArray(vy); + +%% Integration on GPU +N=1; +t0=0; + +[x_f,y_f,vx_f,vy_f]=arrayfun(@RKF45_FILE_gpu,t0,N,x_gpu,y_gpu,vx_gpu,vy_gpu,mu); + +%% Data back to CPU and GPU memory cleaning +clear x_gpu y_gpu vx_gpu vy_gpu +x_T=gather(x_f); +clear x_f +y_T=gather(y_f); +clear y_f +vx_T=gather(vx_f); +clear vx_f +vy_T=gather(vy_f); +clear vy_f + +%% Construction of matrix for FTLE computation + +X_T=zeros(nx,ny,nvx,ne); +Y_T=zeros(nx,ny,nvx,ne); +VX_T=zeros(nx,ny,nvx,ne); +VY_T=zeros(nx,ny,nvx,ne); +E_T=zeros(nx,ny,nvx,ne); +m=1; +for i=1:nx + for j=1:ny + for k=1:nvx + for l=1:ne + if filter(i,j,k,l)==1 + X_T(i,j,k,l)=x_T(m); + Y_T(i,j,k,l)=y_T(m); + VX_T(i,j,k,l)=vx_T(m); + VY_T(i,j,k,l)=vy_T(m); + E_T(i,j,k,l)=0.5*(VX_T(i,j,k,l)^2+VY_T(i,j,k,l)^2)-Omega(X_T(i,j,k,l),Y_T(i,j,k,l),mu); + m=m+1; + end + end + end + end +end + +%% Compute filter for FTLE +filter_ftle=filter; +for i=2:(nx-1) + for j=2:(ny-1) + for k=2:(nvx-1) + for l=2:(ne-1) + if filter(i,j,k,l)==0 || filter (i,j,k,l)==3 + filter_ftle(i,j,k,l)=0; + + filter_ftle(i+1,j,k,l)=0; + filter_ftle(i-1,j,k,l)=0; + + filter_ftle(i,j+1,k,l)=0; + filter_ftle(i,j-1,k,l)=0; + + filter_ftle(i,j,k+1,l)=0; + filter_ftle(i,j,k-1,l)=0; + + filter_ftle(i,j,k,l+1)=0; + filter_ftle(i,j,k,l-1)=0; + end + end + + end + end +end +%% FTLE computation + +[ftle, dphi]=Compute_FILE_gpu( X_T, Y_T, VX_T, E_T, dx, dy, dvx, de, N, filter_ftle); + +%% Plot results +figure +FTLE=squeeze(ftle(:,2,:,2)); +FTLE(1,:)=[]; +% FTLE(2,:)=[]; +FTLE(:,1)=[]; +% FTLE(:,2)=[]; +x_0(1)=[]; +vx_0(1)=[]; +pcolor(x_0, vx_0, FTLE') +shading flat +toc \ No newline at end of file diff --git a/samples/Matlab/test_rk_par.m b/samples/Matlab/test_rk_par.m new file mode 100644 index 00000000..31564189 --- /dev/null +++ b/samples/Matlab/test_rk_par.m @@ -0,0 +1,19 @@ +clear +mu=0.1; +x_0=linspace(-0.8, -0.15, 2) +y_0=zeros(1,2) +vx_0=linspace(-2, 2, 2) +vy_0=zeros(1,2) +ci=[1-mu-0.05 0 0.005 0.5290] +t0=[0;0] +T=[2;2] +tspan=2 +arg1={@f;@f} +%tspan={[0 2],[0 2]}; +arg=[mu;mu] +[X]=arrayfun(RK4_par,t0,T,x_0',y_0',vx_0',vy_0',arg) +% [X]=arrayfun(@f,[0;1],[ci;ci],[mu;mu]); +%Y=RK4(@f,tspan,ci,mu); +% figure +% plot(Y(:,1),Y(:,2)) +% Y(end,1) \ No newline at end of file