diff --git a/matlab/HD_study.m b/matlab/HD_study.m
deleted file mode 100644
index 9ccde4cd59a3daf79fc99cd7460dc332b583b750..0000000000000000000000000000000000000000
--- a/matlab/HD_study.m
+++ /dev/null
@@ -1,200 +0,0 @@
-red_ = [0.6350 0.0780 0.1840];
-gre_ = [0.4660 0.6740 0.1880];
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% DOUGHERTY
-% nu mu
-% bursts simulations
-b_=[...
- 5e-1, 5e+0;...
- 1e-1, 1e-1;...
- 1e-1, 3e-2;...
- 1e-1, 2e-2;...
- 1e-1, 1e-2;...
- 5e-2, 3e-3;...
- 5e-2, 2e-3;...
- 1e-2, 1e-2;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_DGGK_CLOS_0_mu_1e-02/
- 1e-2, 3e-3;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-02_DGGK_mu_3e-03
- 1e-2, 5e-4;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_mu_5e-04
- 1e-2, 1e-4;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_mu_3e-03
- 1e-3, 2e-2;... % v2.7_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-03_DGGK_CLOS_0_mu_2e-02/
- ];
-% converged turb plateau simulations
-cp_=[...
- 1e+0, 1e-2;...
- 1e+0, 3e-3;...
- 1e+0, 2e-3;...
- 5e-1, 2e-3;...
- 5e-1, 2e-2;...
- 1e-1, 2e-3;...
- 1e-1, 1e-3;...
- 5e-2, 1e-3;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_DGGK_mu_3e-02
- 5e-2, 5e-4;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-02_DGGK_mu_5e-04
- ];
-% moving/no turb plateau
-dp_=[...
- 1e+0, 5e+0;...
- 1e+0, 2e+0;...
- 1e+0, 1e+0;...
- 1e-3, 2e-2;...
- 1e-3, 1e-2;...
- 1e-3, 5e-3;...
- 1e-3, 2.5e-3;...
- 1e-3, 1e-3;... %v2.7_P_6_J_3/200x100_L_60_P_6_J_3_eta_0.6_nu_1e-03_DGGK_CLOS_0_mu_1e-03/
- ];
-% not sure
-rp_=[...
- 0,0;...
- ];
-figure; set(gcf, 'Position', [100, 100, 900, 400])
-title('Hyperdiffusion study, Dougherty GK')
-grid on; xlim([5e-4,5e0]); ylim([5e-5,5e+0]);
-set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-xlabel('$\nu$'); ylabel('$\mu_{HD}$'); hold on;
-
-% Trajectory of simulations
-if 0
-% HD_study/200x100_L_200_P_2_J_1_eta_0.6_nu_1e+00_DGGK_CLOS_0_mu_1e-02/
-mu_ = [0 1 2 5 5];
-nu_ = [1 1 1 1 0.5];
-plot(nu_,mu_,'x--','DisplayName','N=200, L=050, P,J=2,1');
-
-% HD_study/300x150_L_200_P_2_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_1e-02/
-mu_ = [1e-2 1e-3];
-nu_ = [1e-1 1e-1];
-plot(nu_,mu_,'x--','DisplayName','N=300, L=100, P,J=2,1');
-
-% HD_study/300x150_L_200_P_2_J_1_eta_0.6_nu_5e-01_DGGK_CLOS_0_mu_1e-02/
-mu_ = [2e-3 2e-2];
-nu_ = [5e-1 5e-1];
-plot(nu_,mu_,'x--','DisplayName','N=300, L=100, P,J=2,1');
-
-% HD_study/100x50_L_50_P_2_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_1e-02/
-mu_ = [1e-2 5e-3 2e-3];
-nu_ = [1e-1 1e-1 1e-1];
-plot(nu_,mu_,'x-','DisplayName','N=100, L=050, P,J=2,1');
-
-% HD_study/150x75_L_100_P_2_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_3e-02/
-mu_ = [3e-2 1e-3 0];
-nu_ = [1e-1 1e-1 1e-1];
-plot(nu_,mu_,'x--','DisplayName','N=150, L=100, P,J=2,1');
-
-% HD_study/150x75_L_100_P_2_J_1_eta_0.6_nu_1e-02_DGGK_CLOS_0_mu_3e-02/
-mu_ = [3e-3 1e-3 1e-4];
-nu_ = [1e-2 1e-2 1e-2];
-plot(nu_,mu_,'x--','DisplayName','N=150, L=100, P,J=2,1');
-
-% HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_CLOS_0_mu_3e-02/
-mu_ = [3e-3 5e-4 1e-4];
-nu_ = [1e-2 1e-2 1e-2];
-plot(nu_,mu_,'x--','DisplayName','N=150, L=100, P,J=4,2');
-
-% HD_study/150x75_L_100_P_2_J_1_eta_0.6_nu_5e-02_DGGK_CLOS_0_mu_3e-02/
-mu_ = [3e-3 1e-3];
-nu_ = [5e-2 5e-2];
-plot(nu_,mu_,'x--','DisplayName','N=150, L=100, P,J=2,1');
-
-% v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/
-mu_ = [2e-2];
-nu_ = [1e-1];
-% plot(nu_,mu_,'x--','DisplayName','N=200, L=120, P,J=6,3');
-end
-scatter( b_(:,1), b_(:,2),'o',...
- 'MarkerFaceColor',red_,'MarkerEdgeColor',[0 0 0],'SizeData',50,...
- 'DisplayName','Bursts');
-scatter(cp_(:,1),cp_(:,2),'s',...
- 'MarkerFaceColor',gre_,'MarkerEdgeColor',[0 0 0],'SizeData',50,...
- 'DisplayName','Converged Plateau');
-scatter(dp_(:,1),dp_(:,2),'d',...
- 'MarkerFaceColor',gre_,'MarkerEdgeColor',[0 0 0],'SizeData',50,...
- 'DisplayName','Moving Plateau');
-scatter(rp_(:,1),rp_(:,2),'h',...
- 'MarkerFaceColor',[0.9290 0.6940 0.1250],'MarkerEdgeColor',[0 0 0],'SizeData',50,...
- 'DisplayName','not sure');
-plot(0,0,'v','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0], 'DisplayName','$\mu=0$');
-legend('show','Location','NorthWest')
-scatter(1,5e-5,80,'v','MarkerFaceColor',gre_,'MarkerEdgeColor',[0 0 0]);
-% HD_study/150x75_L_100_P_2_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_3e-02/
-scatter(0.1,5e-5,80,'v','MarkerFaceColor',gre_,'MarkerEdgeColor',[0 0 0]);
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% SUGAMA
-% nu mu
-% bursts simulations
-b_=[...
- 1e+0, 1.0e-2;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e+00_SGGK_CLOS_0_mu_1e-02/
- 1e+0, 3.0e-3;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e+00_SGGK_CLOS_0_mu_1e-02/
- 5e-1, 3.0e-2;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_3e-02
- 5e-1, 1.6e-2;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_3e-02
- 5e-1, 0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt
- 1e-1, 2.0e-2;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-01_SGGK_CLOS_0_mu_2e-02/
- 1e-1, 1.6e-2;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_GGK_mu_3e-02/
- 1e-2, 5.0e-3;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-02/
- 1e-2, 3.0e-3;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-02/
- ];
-% converged turb plateau simulations
-cp_=[...
- 1e-1, 0;... % v2.7_P_2_J_1/200x100_L_120_P_2_J_1_eta_0.6_nu_1e-01_SGGK_CLOS_0_mu_0e+00/
- ];
-% moving/no turb plateau
-dp_=[...
- 1e-1, 0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt
- 1e-2, 1.0e-2;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-02/
- 1e-2, 5.0e-3;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-02/
- 1e-2, 3.0e-3;... % v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-02/
- 1e-3, 1.0e-3;... % v2.7_P_6_J_3/200x100_L_60_P_6_J_3_eta_0.6_nu_1e-03_SGGK_CLOS_0_mu_1e-03/
- ];
-% not sure
-rp_=[...
- 1e-1, 3.0e-2;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_3e-02
- 5e-2, 1.0e-3;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_DGGK_mu_3e-02
- ];
-figure; set(gcf, 'Position', [100, 100, 900, 400])
-title('Hyperdiffusion study, Sugama GK')
-grid on; xlim([5e-4,5e0]); ylim([5e-5,5e+0]);
-set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-xlabel('$\nu$'); ylabel('$\mu_{HD}$'); hold on;
-
-% Trajectory of simulations
-if 0
-% v2.7_P_2_J_1/200x100_L_120_P_2_J_1_eta_0.6_nu_1e-01_SGGK_CLOS_0_mu_0e+00/
-mu_ = [0.0];
-nu_ = [0.1];
-plot(nu_,mu_,'x--','DisplayName','N=200, L=050, P,J=2,1');
-
-% Trajectory of simulations
-% v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-01_SGGK_CLOS_0_mu_2e-02/
-mu_ = [0.02];
-nu_ = [0.1];
-plot(nu_,mu_,'x--','DisplayName','N=200, L=050, P,J=2,1');
-
-mu_ = [0.0];
-nu_ = [0.1];
-plot(nu_,mu_,'x--','DisplayName','N=200, L=050, P,J=2,1');
-end
-
-scatter( b_(:,1), b_(:,2),'o',...
- 'MarkerFaceColor',red_,'MarkerEdgeColor',[0 0 0],'SizeData',80,...
- 'DisplayName','Bursts');
-% scatter(cp_(:,1),cp_(:,2),'s',...
-% 'MarkerFaceColor',gre_,'MarkerEdgeColor',[0 0 0],'SizeData',60,...
-% 'DisplayName','Converged Plateau');
-scatter(dp_(:,1),dp_(:,2),'d',...
- 'MarkerFaceColor',gre_,'MarkerEdgeColor',[0 0 0],'SizeData',60,...
- 'DisplayName','Plateau');
-% scatter(rp_(:,1),rp_(:,2),'h',...
-% 'MarkerFaceColor',[0.9290 0.6940 0.1250],'MarkerEdgeColor',[0 0 0],'SizeData',60,...
-% 'DisplayName','not sure');
-scatter(0,0,'v','MarkerFaceColor',[0 0 0],'MarkerEdgeColor',[0 0 0],...
- 'DisplayName','$\mu=0$');
-legend('show','Location','NorthWest')
-
-scatter(0.075,5e-5,80,'v','MarkerFaceColor',gre_,'MarkerEdgeColor',[0 0 0],...
- 'DisplayName','$\mu=0$');
-scatter(0.1,5e-5,80,'v','MarkerFaceColor',gre_,'MarkerEdgeColor',[0 0 0],...
- 'DisplayName','$\mu=0$');
-scatter(0.25,5e-5,80,'v','MarkerFaceColor',red_,'MarkerEdgeColor',[0 0 0],...
- 'DisplayName','$\mu=0$');
-scatter(0.5,5e-5,80,'v','MarkerFaceColor',red_,'MarkerEdgeColor',[0 0 0],...
- 'DisplayName','$\mu=0$');
-scatter(1.0,5e-5,80,'v','MarkerFaceColor',red_,'MarkerEdgeColor',[0 0 0],...
- 'DisplayName','$\mu=0$');
diff --git a/matlab/assemble_cosolver_matrices.m b/matlab/assemble_cosolver_matrices.m
deleted file mode 100644
index 2dad7e0bc559ab23543810b2df56a5a964907c67..0000000000000000000000000000000000000000
--- a/matlab/assemble_cosolver_matrices.m
+++ /dev/null
@@ -1,71 +0,0 @@
-codir = '/home/ahoffman/cosolver/';
-matname= 'gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0';
-matdir = [codir,matname,'/'];
-kperp = load([matdir,'kperp.log']);
-Nkp = numel(kperp);
-outdir = '/home/ahoffman/HeLaZ/iCa/';
-outfile= [outdir,matname,'.h5'];
-
-if(exist(outfile)>0)
- system(['rm ',outdir,matname,'.h5']);
-end
-
-Nmax = numel(kperp);
-for n_=1:Nmax
- n_string = sprintf('%5.5d',n_-1); disp(n_string);
- scandir = ['scanfiles_',n_string,'/'];
- if(exist([matdir,scandir, 'ei.h5'])>0)
- % Load matrices and write them in one h5 file
-
- olddname = '/Ceipj/CeipjF'; infile = [matdir,scandir, 'ei.h5'];
- newdname = ['/',n_string,olddname];
- load_write_mat(infile,olddname,outfile,newdname);
-
- olddname = '/Ceipj/CeipjT'; infile = [matdir,scandir, 'ei.h5'];
- newdname = ['/',n_string,olddname];
- load_write_mat(infile,olddname,outfile,newdname);
-
- olddname = '/Ciepj/CiepjT'; infile = [matdir,scandir, 'ie.h5'];
- newdname = ['/',n_string,olddname];
- load_write_mat(infile,olddname,outfile,newdname);
-
- olddname = '/Ciepj/CiepjF'; infile = [matdir,scandir, 'ie.h5'];
- newdname = ['/',n_string,olddname];
- load_write_mat(infile,olddname,outfile,newdname);
-
- olddname = '/Caapj/Ceepj'; infile = [matdir,scandir,'self.h5'];
- newdname = ['/',n_string,olddname];
- mat_ee = load_write_mat(infile,olddname,outfile,newdname);
-
- olddname = '/Caapj/Ciipj'; infile = [matdir,scandir,'self.h5'];
- newdname = ['/',n_string,olddname];
- mat_ii = load_write_mat(infile,olddname,outfile,newdname);
-
- % to verify symmetry
- verif = max(abs(imag(eig(mat_ee)))) + max(abs(imag(eig(mat_ii))));
- if(verif>0) disp(['Warning, non sym matrix at ', n_string]); end;
- % to verify negative EV
- verif = max(real(eig(mat_ee))) + max(real(eig(mat_ii)));
- if(verif>0) disp(['Warning, positive EV at ', n_string]); end;
- else
- break
- end
-end
-olddname = '/Ceipj/CeipjF'; infile = [matdir,scandir, 'ei.h5'];
-Pmaxe = h5readatt(infile,olddname,'Pmaxe');
-Jmaxe = h5readatt(infile,olddname,'Jmaxe');
-Pmaxi = h5readatt(infile,olddname,'Pmaxi');
-Jmaxi = h5readatt(infile,olddname,'Jmaxi');
-dims_e = [0 0];
-dims_e(1)= Pmaxe; dims_e(2) = Jmaxe;
-dims_i = [Pmaxi; Jmaxi];
-h5create(outfile,'/dims_e',numel(dims_e));
-h5write (outfile,'/dims_e',dims_e);
-h5create(outfile,'/dims_i',numel(dims_i));
-h5write (outfile,'/dims_i',dims_i);
-%
-h5create(outfile,'/coordkperp',numel(kperp(1:Nmax)));
-h5write (outfile,'/coordkperp',kperp(1:Nmax)');
-fid = H5F.open(outfile);
-
-H5F.close(fid);
\ No newline at end of file
diff --git a/matlab/continue_run.m b/matlab/continue_run.m
deleted file mode 100644
index bf3c6a5303a1540bea9e22e792fe10e4e2e63dd5..0000000000000000000000000000000000000000
--- a/matlab/continue_run.m
+++ /dev/null
@@ -1,92 +0,0 @@
-%% Functions to modify preexisting fort.90 input file and launch on marconi
-function [] = continue_run(outfilename)
- EXECNAME = 'helaz_3.8';
- %% CLUSTER PARAMETERS
- CLUSTER.PART = 'prod'; % dbg or prod
- CLUSTER.TIME = '24:00:00'; % allocation time hh:mm:ss
- if(strcmp(CLUSTER.PART,'dbg')); CLUSTER.TIME = '00:30:00'; end;
- CLUSTER.MEM = '64GB'; % Memory
- CLUSTER.JNAME = 'HeLaZ';% Job name
- NP_P = 2; % MPI processes along p
- NP_KX = 24; % MPI processes along kx
- % Compute processes distribution
- Ntot = NP_P * NP_KX;
- Nnodes = ceil(Ntot/48);
- Nppn = Ntot/Nnodes;
- CLUSTER.NODES = num2str(Nnodes); % MPI process along p
- CLUSTER.NTPN = num2str(Nppn); % MPI process along kx
- CLUSTER.CPUPT = '1'; % CPU per task
- %%
- RESDIR = ['../',outfilename(46:end-8),'/'];
- BASIC.RESDIR = RESDIR;
- FORT90 = [RESDIR,'fort.90'];
-
- % Read txt into cell A
- fid = fopen(FORT90,'r');
- i = 1;
- tline = fgetl(fid);
- A{i} = tline;
- while ischar(tline)
- i = i+1;
- tline = fgetl(fid);
- A{i} = tline;
- end
- fclose(fid);
-
- % Find previous job2load
- if( numel(A{5}) == numel(' RESTART = .false.') )
- A{5} = sprintf(' RESTART = .true.');
- J2L = 0;
- else
- line = A{39};
- line = line(end-2:end);
- if(line(1) == '='); line = line(end); end;
- J2L = str2num(line)+1;
- end
- % Change job 2 load in fort.90
- A{39} = [' job2load = ',num2str(J2L)];
- disp(A{39})
- % Change time step
- line_= A{3};
- dt_old = str2num(line_(13:end));
- A{3} = [' dt = ',num2str(dt_old)];
- % Increase endtime
- A{4} = [' tmax = 20000'];
- % Change collision operator
- line_= A{43};
- CO_old = str2num(line_(13:end));
- A{43} = [' CO = ',num2str(2)];
- % Put non linear term back
- A{45} = [' NL_CLOS = -1'];
- % change HD
- line_= A{47};
- mu_old = str2num(line_(13:end));
- A{47} = [' mu = ',num2str(0)];
- % change L
- line_= A{14};
- L_old = str2num(line_(12:end));
- A{14} = [' Lx = ',num2str(L_old)];
- A{16} = [' Ly = ',num2str(L_old)];
- % change eta N
- line_= A{57};
- etan_old = str2num(line_(13:end));
- A{57} = [' eta_n = ',num2str(etan_old)];
- % change eta B
- line_= A{59};
- etab_old = str2num(line_(13:end));
- A{59} = [' eta_B = ',num2str(etab_old)];
- % Rewrite fort.90
- fid = fopen('fort.90', 'w');
- for i = 1:numel(A)
- if A{i+1} == -1
- fprintf(fid,'%s', A{i});
- break
- else
- fprintf(fid,'%s\n', A{i});
- end
- end
- % Copy fort.90 into marconi
- write_sbash_marconi
- % Launch the job
- system('ssh ahoffman@login.marconi.cineca.it sh HeLaZ/wk/setup_and_run.sh');
-end
diff --git a/matlab/dbg_zBC_map.m b/matlab/dbg_zBC_map.m
deleted file mode 100644
index 63de86cbfff105b1e2fd9cdcaf260ccca4451dec..0000000000000000000000000000000000000000
--- a/matlab/dbg_zBC_map.m
+++ /dev/null
@@ -1,84 +0,0 @@
-Nkx = 6;
-Nky = 3;
-
-my = 0:(Nky-1);
-mx = zeros(1,Nkx);
-
-PERIODIC = 1;
-Npol = 1;
-Nexc = 0;
-
-shear = 0.8;
-Ly = 120;
-Lx = 120;
-dky = 2*pi/Ly;
-
-
-for ix = 1:Nkx
- if(mod(Nkx,2) == 0)%even
- mx_max = (Nkx/2);
- if(ix<=Nkx/2+1)
- mx(ix) = (ix-1);
- else
- mx(ix) = ix-Nkx-1;
- end
- else %odd
- mx_max = (Nkx-1)/2;
- if(ix<=(Nkx-1)/2+1)
- mx(ix) = (ix-1);
- else
- mx(ix) = ix-Nkx-1;
- end
- end
-end
-disp(mx)
-
-
-if Nexc == 0 %% Adapt Nexc
- dkx = 2*pi/Lx;
- Nexc = ceil(0.9*2*pi*shear*dky/dkx);
-else
- dkx = 2*pi*shear*dky/Nexc;
-end
-
-kx = mx*dkx;
-ky = my*dky;
-
-kx_max = mx_max*dkx;
-ikx_zBC_R = zeros(Nky,Nkx);
-for iy = 1:Nky
- shift = 2*pi*shear*ky(iy)*Npol;
- for ix = 1:Nkx
- kx_shift = kx(ix) + shift;
- if ((kx_shift > kx_max) && (~PERIODIC))
- ikx_zBC_R(iy,ix) = nan;
- else
- ikx_zBC_R(iy,ix) = ix+(iy-1)*Nexc;
- if(ikx_zBC_R(iy,ix) > Nkx)
-% ikx_zBC_R(iy,ix) = ikx_zBC_R(iy,ix) - Nkx;
- ikx_zBC_R(iy,ix) = mod(ikx_zBC_R(iy,ix)-1,Nkx)+1;
- end
- end
- end
-end
-disp(ikx_zBC_R)
-
-kx_min = (-mx_max+(1-mod(Nkx,2)))*dkx;
-ikx_zBC_L = zeros(Nky,Nkx);
-for iy = 1:Nky
- shift = 2*pi*shear*ky(iy)*Npol;
- for ix = 1:Nkx
- kx_shift = kx(ix) - shift;
- if ((kx_shift < kx_min) && (~PERIODIC))
- ikx_zBC_L(iy,ix) = nan;
- else
- ikx_zBC_L(iy,ix) = ix-(iy-1)*Nexc;
- if(ikx_zBC_L(iy,ix) < 1)
-% ikx_zBC_L(iy,ix) = ikx_zBC_L(iy,ix) + Nkx;
- ikx_zBC_L(iy,ix) = mod(ikx_zBC_L(iy,ix)-1,Nkx)+1;
- end
- end
- end
-end
-disp(ikx_zBC_L)
-
diff --git a/matlab/evolve_tracers.m b/matlab/evolve_tracers.m
index 5b7a4b402187cd864cca2be93f793b9ed8533649..bfa78656322403cb2adbbc551dad4d945f1480bf 100644
--- a/matlab/evolve_tracers.m
+++ b/matlab/evolve_tracers.m
@@ -1,27 +1,36 @@
% Options
SHOW_FILM = 1;
-field2plot ='phi';
-INIT = 'lin'; % lin (for a line)/ round (for a small round)/ gauss for random
-U_TIME = 500; % >0 for frozen velocity at a given time, -1 for evolving field
+field2plot ='phi';
+% INIT = 'lin'; % lin (for a line)/ round (for a small round)/ gauss for random
+INIT = 'gauss'; % lin (for a line)/ round (for a small round)/ gauss for random
+U_TIME = 50; % >0 for frozen velocity at a given time, -1 for evolving field
Evolve_U = 1; % 0 for frozen velocity at a given time, 1 for evolving field
-Tfin = 2000;
+Tfin = 300;
dt_ = 0.1;
Nstep = ceil(Tfin/dt_);
% Init tracers
-Np = 3; %number of tracers
+Np = 10; %number of tracers
% color = tcolors;
color = jet(Np);
tcolors = distinguishable_colors(Np); %Their colors
dimmed = 0; % to dimm the colormap in the background (infty = white, 0 normal color)
-Na = 1000/dt_; %length of trace
+Na = 10/dt_; %length of trace
+% load data
+[data.PHI, data.Ts3D] = compile_results_3D (DATADIR,J0,J1,'phi');
+[data.DENS_ , ~] = compile_results_3Da(DATADIR,J0,J1,'dens');
+[data.Na00, ~] = compile_results_3Da(DATADIR,J0,J1,'Na00');
+
+data.Ni00 = reshape( data.Na00(1,:,:,:,:),size(data.PHI));
+data.DENS_I = reshape(data.DENS_(1,:,:,:,:),size(data.PHI));
+% setup
Traj_x = zeros(Np,Nstep);
Traj_y = zeros(Np,Nstep);
Disp_x = zeros(Np,Nstep);
Disp_y = zeros(Np,Nstep);
-xmax = max(data.x); xmin = min(data.x);
-ymax = max(data.y); ymin = min(data.y);
+xmax = max(data.grids.x); xmin = min(data.grids.x);
+ymax = max(data.grids.y); ymin = min(data.grids.y);
switch INIT
case 'lin'
@@ -48,8 +57,9 @@ switch INIT
end
% position grid and velocity field
-[YY_, XX_ ,ZZ_] = meshgrid(data.y,data.x,data.z);
-[KX,KY] = meshgrid(data.kx,data.ky);
+% [YY_, XX_ ,ZZ_] = meshgrid(data.grids.y,data.grids.x,data.grids.z);
+[YY_, XX_] = meshgrid(data.grids.y,data.grids.x);
+[KX,KY] = meshgrid(data.grids.kx,data.grids.ky);
Ux = zeros(size(XX_));
Uy = zeros(size(XX_));
Uz = zeros(size(XX_));
@@ -57,10 +67,7 @@ ni = zeros(size(XX_));
[~,itu_] = min(abs(U_TIME-data.Ts3D));
% computing the frozen velocity field
-for iz = 1:data.Nz
-% Ux(:,:,iz) = real(ifft2( 1i*KY.*(data.PHI(:,:,iz,itu_)),data.Nx,data.Ny));
-% Uy(:,:,iz) = real(ifft2(-1i*KX.*(data.PHI(:,:,iz,itu_)),data.Nx,data.Ny));
-% ni(:,:,iz) = real(ifft2(data.DENS_I(:,:,iz,itu_),data.Nx,data.Ny));
+for iz = 1:data.grids.Nz
Ux(:,:,iz) = ifourier_GENE( 1i*KY.*(data.PHI(:,:,iz,itu_)))';
Uy(:,:,iz) = ifourier_GENE(-1i*KX.*(data.PHI(:,:,iz,itu_)))';
ni(:,:,iz) = ifourier_GENE(data.DENS_I(:,:,iz,itu_))';
@@ -92,11 +99,10 @@ while(t_<Tfin && it <= Nstep)
end
if Evolve_U && (itu_old ~= itu_)
% updating the velocity field and density field
- for iz = 1:data.Nz
+ for iz = 1:data.grids.Nz
Ux(:,:,iz) = ifourier_GENE( 1i*KY.*(data.PHI(:,:,iz,itu_)))';
Uy(:,:,iz) = ifourier_GENE(-1i*KX.*(data.PHI(:,:,iz,itu_)))';
-% ni(:,:,iz) = ifourier_GENE(data.DENS_I(:,:,iz,itu_))';
- ni(:,:,iz) = ifourier_GENE(data.Ni00(:,:,iz,itu_))';
+ ni(:,:,iz) = ifourier_GENE(data.DENS_I(:,:,iz,itu_))';
end
end
% evolve each tracer
@@ -104,11 +110,11 @@ while(t_<Tfin && it <= Nstep)
% locate the tracer
% find corners of the cell
x_ = Xp(ip);
- [e_x,ixC] = min(abs(x_-data.x));
+ [e_x,ixC] = min(abs(x_-data.grids.x));
if e_x == 0 % on the face
ix0 = ixC;
ix1 = ixC;
- elseif x_ > data.x(ixC) % right from grid point
+ elseif x_ > data.grids.x(ixC) % right from grid point
ix0 = ixC;
ix1 = ixC+1;
else % left
@@ -116,11 +122,11 @@ while(t_<Tfin && it <= Nstep)
ix1 = ixC;
end
y_ = Yp(ip);
- [e_y,iyC] = min(abs(y_-data.y));
+ [e_y,iyC] = min(abs(y_-data.grids.y));
if e_y == 0 % on the face
iy0 = iyC;
iy1 = iyC;
- elseif y_ > data.y(iyC) % above
+ elseif y_ > data.grids.y(iyC) % above
iy0 = iyC;
iy1 = iyC+1;
else % under
@@ -128,20 +134,20 @@ while(t_<Tfin && it <= Nstep)
iy1 = iyC;
end
z_ = Zp(ip,1);
- [e_z,izC] = min(abs(z_-data.z));
+ [e_z,izC] = min(abs(z_-data.grids.z));
if e_z == 0 % on the face
iz0 = izC;
iz1 = izC;
- elseif z_ > data.z(izC) % before
+ elseif z_ > data.grids.z(izC) % before
iz0 = izC;
iz1 = izC+1;
else % behind
iz0 = izC-1;
iz1 = izC;
end
- x0 = data.x(ix0); x1 = data.x(ix1); %left right
- y0 = data.y(iy0); y1 = data.y(iy1); %down top
- z0 = data.z(iz0); z1 = data.z(iz1); %back front
+ x0 = data.grids.x(ix0); x1 = data.grids.x(ix1); %left right
+ y0 = data.grids.y(iy0); y1 = data.grids.y(iy1); %down top
+ z0 = data.grids.z(iz0); z1 = data.grids.z(iz1); %back front
if(e_x > 0)
ai__ = (x_ - x0)/(x1-x0); % interp coeff x
else
@@ -182,7 +188,6 @@ while(t_<Tfin && it <= Nstep)
% push the particle
q = sign(-u___(3));
-% q =1;
x_ = x_ + dt_*u___(1)*q;
y_ = y_ + dt_*u___(2)*q;
diff --git a/matlab/fig_post_processing.m b/matlab/fig_post_processing.m
deleted file mode 100644
index e51be5273ac9e66874f84050ffbb6185da7d4ea1..0000000000000000000000000000000000000000
--- a/matlab/fig_post_processing.m
+++ /dev/null
@@ -1,78 +0,0 @@
-%% Load figure
-figpath = 'C:\Users\antoi\Desktop\gamma_eta_05_nu_1e-01_trunc';
-fig = openfig(figpath);
-
-%% Load data
-axObjs = fig.Children;
-dataObjs = findobj(fig,'-property','YData');
-Nlines = numel(dataObjs);
-
-%% Post processing
-sigma = zeros(Nlines,1);
-mu = sigma;
-tmin = 350;
-for i = 1:Nlines
- x = dataObjs(i).XData;
-
- [~, itmin] = min(abs(tmin-x));
-
- y = dataObjs(i).YData;
-
- sigma(i) = std(y(itmin:end));
- mu(i) = mean(y(itmin:end));
-end
-mu = flip(mu')
-sigma = flip(sigma')
-
-%% Plot mean with error bar
-
-
-
-
-if 0
-%% Handwritten results for nu = 1.0, 150x75, L=100, DGGK
-Results_150x75.Gamma = [0.3794 0.3194 0.3226, 0.0098 0.0221 0.0321 0.0400 0.2897 0.2886 0.2569 0.0104 0.0086 0.0276 0.0320 0.1375 0.1633 0.0848];
-Results_150x75.error = [0.1909 0.1207 0.1336, 0.0028 0.0038 0.0058 0.0086 0.0832 0.0624 0.0557 0.0021 0.0023 0.0068 0.0088 0.0821 0.0278 0.0083];
-Results_150x75.P = [2, 4, 6, 2, 4, 6 8 2, 4, 6 2, 4, 6 8 2, 4, 10];
-Results_150x75.J = [1, 2, 3 1 2 3 4 1, 2, 3 1, 2, 3 4 1, 2 5];
-Results_150x75.etaB = [0.49, 0.49 0.49 0.59 0.59 0.59 0.59 0.50, 0.50 0.50 0.60, 0.60 0.60 0.60 0.51 0.51 0.51];
-Results_150x75.nu = [1.0, 1.0 1.0 1.0 1.0 1.0 1.0 0.5, 0.5 0.5 0.5, 0.5 0.5 0.5 0.1 0.1 0.1];
-Results_150x75.mrkx = [ '*', '*', '*', '*', '*', '*', '*', 'o', 'o', 'o', 'o', 'o', 'o', 'o' 's' 's' 's'];
-Results_150x75.iclr = [ 1, 2, 3, 1 2 3 4 1, 2, 3 1, 2, 3 4 1 2 5];
-
-% Ricci_Rogers.Gamma = [2 1e-1];
-% Ricci_Rogers.etaB = [0.5 1.0];
-Ricci_Rogers.Gamma = [10 1e-2];
-Ricci_Rogers.etaB = [0.5 1.25];
-
-if 1
-% Fig 3 of Ricci Rogers 2006
-SCALING = 2*sqrt(2);
-fig = figure;
-semilogy(Ricci_Rogers.etaB,Ricci_Rogers.Gamma,'--','color',[0,0,0]+0.6);
-hold on;
-plot(10,10,'color',line_colors(1,:)); plot(10,10,'color',line_colors(2,:));
-plot(10,10,'color',line_colors(3,:)); plot(10,10,'color',line_colors(4,:));
-plot(10,10,'color',line_colors(5,:));
-plot(10,10,'*k','MarkerSize',10, 'LineWidth',1.0);
-plot(10,10,'ok','MarkerSize',10, 'LineWidth',1.0);
-plot(10,10,'sk','MarkerSize',10, 'LineWidth',1.0);
-
-res = Results_150x75;
-for i = 1:numel(res.Gamma)
- errorbar(res.etaB(i),res.Gamma(i)*SCALING,res.error(i)*SCALING,...
- res.mrkx(i),'DisplayName','256x128', 'color', line_colors(res.iclr(i),:),...
- 'MarkerSize',12, 'LineWidth',2.0);
- hold on;
-end
-
- xlabel('$L_n/L_B$'); ylabel('$2\sqrt(2)\Gamma^\infty_{part}$')
-end
-grid on; title('$L=100$, $150\times75$, $\nu_{hyp}=0.1$');
-xlim([0,1.75]); ylim([1e-6,100])
-legend('Mix. Length, Ricci 2006','$P=2$, $J=1$','$P=4$, $J=2$','$P=6$, $J=3$','$P=8$, $J=4$','$P=10$, $J=5$',...
- '$\nu_{DGGK}=1.0$', '$\nu_{DGGK}=0.5$','$\nu_{DGGK}=0.1$');
-plot([0.3 0.3],[1e-6,1e2],'r')
-plot([1.6 1.6],[1e-6,1e2],'r')
-plot([0.5],[0.3965],'--','color',[0,0,0]+0.6)
-end
\ No newline at end of file
diff --git a/matlab/flux_results.m b/matlab/flux_results.m
deleted file mode 100644
index 594b3b62b5a69226a163cfca70774bd147354d6c..0000000000000000000000000000000000000000
--- a/matlab/flux_results.m
+++ /dev/null
@@ -1,171 +0,0 @@
-default_plots_options
-if 1
-%% Compute time average and std of the mean flow
-t0 = 1000; t1 = 3000; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D));
-range = it0:it1;
-avg = mean(GFlux_ri(range))
-stdev = std(GFlux_ri(range))^(.5)
-figure
-hist(GFlux_ri(range),20); xlabel('$\Gamma$')
-end
-if 0
-%% Handwritten results for nu = 0.01
-% High definition results (256x128)
-SCALING = 2*sqrt(2);
-Results_256x128.Gamma = [0.02, 0.03, 0.20, 0.037, 2.7, 2.25, 4, 5e-4, 2e-3, 0.03];
-Results_256x128.L = [ 66, 66, 66, 50, 66, 66, 66, 66, 66, 66];
-Results_256x128.P = [ 2, 3, 4, 5, 2, 3, 4, 2, 3, 4];
-Results_256x128.J = [ 1, 2, 2, 3, 1, 2, 2, 1, 2, 2];
-Results_256x128.etaB = [ 0.5, 0.5, 0.5, 0.5, 0.4, 0.4, 0.4, 0.6, 0.6, 0.6];
-Results_256x128.mrkx = [ 'v', '>', '^', 'o', 'v', '>', '^', 'v', '>', '^'];
-Results_256x128.clr = [ 'k', 'k', 'k', 'r', 'r', 'r', 'r', 'k', 'k', 'k'];
-% Low definition results (128x64)
-% Results_128x64.Gamma = [0.29, 0.05, 7e-4, 0.31, 3.7, 2e-3];
-% Results_128x64.L = [ 25, 25, 25, 33 33, 33];
-% Results_128x64.P = [ 2, 2, 2, 2, 2, 2];
-% Results_128x64.J = [ 1, 1, 1, 1, 1, 1];
-% Results_128x64.NU = [0.01, 0.1, 0.01, 0.01, 0.01, 0.01];
-% Results_128x64.etaB = [ 0.5, 0.5, 0.67, 0.5 0.4, 0.6];
-% Results_128x64.mrkx = [ 's', 's', 's', 's', 's', 's'];
-% Results_128x64.clr = [ 'b', 'b', 'b', 'r', 'r', 'r'];
-
-% Ricci_Rogers.Gamma = [2.5 1 1e-2];
-% Ricci_Rogers.etaB = [0.4 0.5 1.0];
-Ricci_Rogers.Gamma = [10 1e-1];
-Ricci_Rogers.etaB = [0.5 1.0];
-if 1
-% Fig 3 of Ricci Rogers 2006
-fig = figure;
-semilogy(Ricci_Rogers.etaB,Ricci_Rogers.Gamma,'--','color',[0,0,0]+0.6);
-hold on;
-res = Results_256x128;
-for i = 1:numel(res.Gamma)
- semilogy(res.etaB(i),res.Gamma(i)*SCALING,...
- res.mrkx(i),'DisplayName','256x128', 'color', res.clr(i));
- hold on;
-end
-% res = Results_128x64;
-% for i = 1:numel(res.Gamma)
-% if res.NU(i) == 0.01
-% semilogy(res.etaB(i),res.Gamma(i),...
-% res.mrkx(i),'DisplayName','128x64', 'color', res.clr(i));
-% end
-% hold on;
-% end
- xlabel('$\eta_B$'); ylabel('$\Gamma^\infty_{part}$')
-end
-grid on; title('$\nu = 0.01$')
-legend('Mix. Length, Ricci 2006','$P=2$, $J=1$','$P=3$, $J=2$','$P=4$, $J=2$','$P=5$, $J=3$')
-FIGNAME = [SIMDIR,'flux_study_nu_1e-2.png']; ylim([1e-4, 10])
-saveas(fig,FIGNAME);
-disp(['Figure saved @ : ',FIGNAME])
-end
-
-if 0
-%% Handwritten results for nu = 0.1
-Results_256x128.Gamma = [0.026,0.026, 1e-2, 1, 1, 1, 2e-2, 1, 0.15, 3e-3];
-Results_256x128.P = [ 2, 3, 4, 2, 3, 4, 2, 3, 4, 4];
-Results_256x128.J = [ 1, 2, 2, 1, 2, 2, 1, 2, 2, 2];
-Results_256x128.etaB = [ 0.5, 0.5, 0.5, 0.4, 0.4, 0.4, 0.6, 0.6, 0.6, 0.7];
-Results_256x128.mrkx = [ 'v', '>', '^', 'v', '>', '^', 'v', '>', '^', '^'];
-Results_256x128.clr = [ 'k', 'k', 'k', 'b', 'b', 'b', 'r', 'b', 'r', 'r'];
-
-% Ricci_Rogers.Gamma = [2 1e-1];
-% Ricci_Rogers.etaB = [0.5 1.0];
-Ricci_Rogers.Gamma = [10 1e-1];
-Ricci_Rogers.etaB = [0.5 1.0];
-
-if 1
-SCALING = 2*sqrt(2);
-% Fig 3 of Ricci Rogers 2006
-fig = figure;
-semilogy(Ricci_Rogers.etaB,Ricci_Rogers.Gamma,'--','color',[0,0,0]+0.6);
-hold on;
-res = Results_256x128;
-for i = 1:numel(res.Gamma)
- semilogy(res.etaB(i),res.Gamma(i)*SCALING,...
- res.mrkx(i),'DisplayName','256x128', 'color', res.clr(i));
- hold on;
-end
-
- xlabel('$\eta_B$'); ylabel('$\Gamma^\infty_{part}$')
-end
-grid on; title('$\nu = 0.1$')
-legend('Mix. Length, Ricci 2006','$P=2$, $J=1$','$P=3$, $J=2$','$P=4$, $J=2$')
-FIGNAME = [SIMDIR,'flux_study_nu_1e-1.png'];
-saveas(fig,FIGNAME);
-disp(['Figure saved @ : ',FIGNAME])
-end
-
-
-if 0
-%% Handwritten results for nu = 1.0, eta = 0.5, 200x100, L=100
-Results_256x128.Gamma = [0.026,0.026, 1e-2, 1, 1, 1, 2e-2, 1, 0.15, 3e-3];
-Results_256x128.P = [ 2, 3, 4, 2, 3, 4, 2, 3, 4, 4];
-Results_256x128.J = [ 1, 2, 2, 1, 2, 2, 1, 2, 2, 2];
-Results_256x128.etaB = [ 0.5, 0.5, 0.5, 0.4, 0.4, 0.4, 0.6, 0.6, 0.6, 0.7];
-Results_256x128.mrkx = [ 'v', '>', '^', 'v', '>', '^', 'v', '>', '^', '^'];
-Results_256x128.clr = [ 'k', 'k', 'k', 'b', 'b', 'b', 'r', 'b', 'r', 'r'];
-
-% Ricci_Rogers.Gamma = [2 1e-1];
-% Ricci_Rogers.etaB = [0.5 1.0];
-Ricci_Rogers.Gamma = [10 1e-1];
-Ricci_Rogers.etaB = [0.5 1.0];
-
-if 1
-% Fig 3 of Ricci Rogers 2006
-SCALING = 2*sqrt(2);
-fig = figure;
-semilogy(Ricci_Rogers.etaB,Ricci_Rogers.Gamma,'--','color',[0,0,0]+0.6);
-hold on;
-res = Results_256x128;
-for i = 1:numel(res.Gamma)
- semilogy(res.etaB(i),res.Gamma(i)*SCALING,...
- res.mrkx(i),'DisplayName','256x128', 'color', res.clr(i));
- hold on;
-end
-
- xlabel('$\eta_B$'); ylabel('$\Gamma^\infty_{part}$')
-end
-grid on; title('$\nu = 0.1$')
-legend('Mix. Length, Ricci 2006','$P=2$, $J=1$','$P=3$, $J=2$','$P=4$, $J=2$')
-FIGNAME = [SIMDIR,'flux_study_nu_1e-1.png'];
-saveas(fig,FIGNAME);
-disp(['Figure saved @ : ',FIGNAME])
-end
-
-if 0
-%% Handwritten results for nu = 1.0, eta = 0.5, 150x75, L=100, DGGK
-Results_150x75.Gamma = [0.026,0.026, 1e-2, 1, 1, 1, 2e-2, 1, 0.15, 3e-3];
-Results_150x75.P = [ 2, 3, 4, 2, 3, 4, 2, 3, 4, 4];
-Results_150x75.J = [ 1, 2, 2, 1, 2, 2, 1, 2, 2, 2];
-Results_150x75.etaB = [ 0.5, 0.5, 0.5, 0.4, 0.4, 0.4, 0.6, 0.6, 0.6, 0.7];
-Results_150x75.mrkx = [ 'v', '>', '^', 'v', '>', '^', 'v', '>', '^', '^'];
-Results_150x75.clr = [ 'k', 'k', 'k', 'b', 'b', 'b', 'r', 'b', 'r', 'r'];
-
-% Ricci_Rogers.Gamma = [2 1e-1];
-% Ricci_Rogers.etaB = [0.5 1.0];
-Ricci_Rogers.Gamma = [10 1e-1];
-Ricci_Rogers.etaB = [0.5 1.0];
-
-if 1
-% Fig 3 of Ricci Rogers 2006
-SCALING = 2*sqrt(2);
-fig = figure;
-semilogy(Ricci_Rogers.etaB,Ricci_Rogers.Gamma,'--','color',[0,0,0]+0.6);
-hold on;
-res = Results_150x75;
-for i = 1:numel(res.Gamma)
- semilogy(res.etaB(i),res.Gamma(i)*SCALING,...
- res.mrkx(i),'DisplayName','256x128', 'color', res.clr(i));
- hold on;
-end
-
- xlabel('$\eta_B$'); ylabel('$\Gamma^\infty_{part}$')
-end
-grid on; title('$\nu = 0.1$')
-legend('Mix. Length, Ricci 2006','$P=2$, $J=1$','$P=3$, $J=2$','$P=4$, $J=2$')
-FIGNAME = [SIMDIR,'flux_study_nu_1e-1.png'];
-saveas(fig,FIGNAME);
-disp(['Figure saved @ : ',FIGNAME])
-end
\ No newline at end of file
diff --git a/matlab/load/compile_results_3Da.m b/matlab/load/compile_results_3Da.m
index 3fffa329aa3ff394b07d8f232ab1a9a723d595ce..7c2e4b89f4fedfe45db15fcee651b125a807e5b7 100644
--- a/matlab/load/compile_results_3Da.m
+++ b/matlab/load/compile_results_3Da.m
@@ -29,6 +29,8 @@ while(CONTINUE)
JOBNUM = JOBNUM + 1;
end
+
+
if(JOBFOUND == 0)
disp('no results found, please verify the paths');
return;
diff --git a/matlab/load/compile_results_low_mem.m b/matlab/load/compile_results_low_mem.m
index 0c80bd79b489899c6d3289c367cbf1eb921e2066..dcc9ad2b55b280d0466c265e265797ea097f1f2c 100644
--- a/matlab/load/compile_results_low_mem.m
+++ b/matlab/load/compile_results_low_mem.m
@@ -145,7 +145,11 @@ else
DATA.grids.ikx0 = ikx0; DATA.grids.iky0 = iky0;
DATA.grids.Nx = Nx; DATA.grids.Ny = Ny; DATA.grids.Nz = Nz; DATA.grids.Nkx = Nkx; DATA.grids.Nky = Nky;
DATA.grids.Np = numel(Parray);DATA.grids.Nj = numel(Jarray);
- DATA.grids.deltap = Parray(2)-Parray(1);
+ if(numel(Parray)>1)
+ DATA.grids.deltap = Parray(2)-Parray(1);
+ else
+ DATA.grids.deltap = 1
+ end
DATA.dir = DIRECTORY;
DATA.localdir = DIRECTORY;
DATA.param_title=['$\nu_{',DATA.inputs.CONAME,'}=$', num2str(DATA.inputs.NU), ...
diff --git a/matlab/load_write_mat.m b/matlab/load_write_mat.m
deleted file mode 100644
index 3315b172c18a9f52f2fea33ade6b55cbce4d4f5b..0000000000000000000000000000000000000000
--- a/matlab/load_write_mat.m
+++ /dev/null
@@ -1,15 +0,0 @@
-function mat_ = load_write_mat(infile,olddname,outfile,newdname)
- mat_ = h5read (infile,olddname);
-
- h5create (outfile,newdname,size(mat_));
- h5write (outfile,newdname,mat_);
-
- if(~strcmp(olddname(end-3:end-2),'ii'))
- neFLR = h5readatt(infile,olddname,'neFLR');
- h5writeatt(outfile,newdname,'neFLR',neFLR);
- end
- if(~strcmp(olddname(end-3:end-2),'ee'))
- niFLR = h5readatt(infile,olddname,'niFLR');
- h5writeatt(outfile,newdname,'niFLR',niFLR);
- end
-end
diff --git a/matlab/marconi_scaling.m b/matlab/marconi_scaling.m
deleted file mode 100644
index b33ec6ed6ae25774ca04eafee3cffb093e6007cd..0000000000000000000000000000000000000000
--- a/matlab/marconi_scaling.m
+++ /dev/null
@@ -1,60 +0,0 @@
-%% Load results np = 24
-np_p_24 = [1 2 3 4];
-np_kx_24= [24 12 8 6];
-el_ti_np_24= zeros(numel(np_p_24),1);
-for i_ = 1:numel(np_p_24)
- npp_ = np_p_24(i_); npkx = np_kx_24(i_);
-
- %% Load from Marconi
- outfile =['/marconi_scratch/userexternal/ahoffman/HeLaZ/results/Marconi_parallel_scaling_2D/',...
- sprintf('%d_%d',npp_,npkx),...
- '_200x100_L_120_P_12_J_5_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-03/out.txt'];
- BASIC.RESDIR = load_marconi(outfile);
- compile_results
- load_params
- el_ti_np_24(i_) = CPUTIME;
-end
-%% Load results np = 48
-np_p_48 = [1 2 3 4 6];
-np_kx_48= [48 24 16 12 8];
-el_ti_np_48= zeros(numel(np_p_48),1);
-for i_ = 1:numel(np_p_48)
- npp_ = np_p_48(i_); npkx = np_kx_48(i_);
-
- %% Load from Marconi
- outfile =['/marconi_scratch/userexternal/ahoffman/HeLaZ/results/Marconi_parallel_scaling_2D/',...
- sprintf('%d_%d',npp_,npkx),...
- '_200x100_L_120_P_12_J_5_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-03/out.txt'];
- BASIC.RESDIR = load_marconi(outfile);
- compile_results
- load_params
- el_ti_np_48(i_) = CPUTIME;
-end
-
-%% Load results np = 72
-np_p_72 = [ 2 3];
-np_kx_72= [36 24];
-el_ti_np_72= zeros(numel(np_p_72),1);
-for i_ = 1:numel(np_p_72)
- npp_ = np_p_72(i_); npkx = np_kx_72(i_);
-
- %% Load from Marconi
- outfile =['/marconi_scratch/userexternal/ahoffman/HeLaZ/results/Marconi_parallel_scaling_2D/',...
- sprintf('%d_%d',npp_,npkx),...
- '_200x100_L_120_P_12_J_5_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-03/out.txt'];
- BASIC.RESDIR = load_marconi(outfile);
- compile_results
- load_params
- el_ti_np_72(i_) = CPUTIME;
-end
-
-%% Plot
-figure
-plt = @(x) (x/el_ti_np_24(1)-1)*100;
-plot(np_p_24,plt(el_ti_np_24),'o--','DisplayName','Ncpu = 24'); hold on;
-plot(np_p_48,plt(el_ti_np_48),'o--','DisplayName','Ncpu = 48 ')
-plot(np_p_72,plt(el_ti_np_72),'o--','DisplayName','Ncpu = 72 ')
-xlabel('Num. proc. p')
-ylabel('Variation from 1 24 [$\%$]')
-title('CPU time change from 1D paralel')
-legend('show')
\ No newline at end of file
diff --git a/matlab/molix_plot_phi.m b/matlab/molix_plot_phi.m
deleted file mode 100644
index b94bd844fa836632e2bae91b3bfd310649f7c59b..0000000000000000000000000000000000000000
--- a/matlab/molix_plot_phi.m
+++ /dev/null
@@ -1,127 +0,0 @@
-function [phib, b_angle] = molix_plot_phi(resfile,varargin)
-% plot ballooning representation at the specified iput by time2plot for
-% each ky modes
-
-
-set(0,'defaulttextInterpreter','latex')
-set(0, 'defaultAxesTickLabelInterpreter','latex');
-set(0, 'defaultLegendInterpreter','latex');
-set(0,'defaultAxesFontSize',16);
-set(0, 'DefaultLineLineWidth', 1.5);
-
-
-% read data and attributes
-coordkx = h5read(resfile,'/data/var2d/phi/coordkx');
-dimskx = size(coordkx);
-Nkx = (dimskx(1) - 1)/2;
-% Nkx = (length(coordkx)-1)/2
-coordz = h5read(resfile,'/data/var2d/phi/coordz');
-coordky = h5read(resfile,'/data/var2d/phi/coordky');
-Nky = length(coordky);
-Nz = length(coordz);
-coordtime =h5read(resfile,'/data/var2d/time');
-Nt = length(coordtime);
-pmax = double(h5readatt(resfile,'/data/input','pmaxi'));
-jmax = double(h5readatt(resfile,'/data/input','jmaxi'));
-epsilon = double(h5readatt(resfile,'/data/input','inv. aspect ratio'));
-safety_fac = double(h5readatt(resfile,'/data/input','safety factor'));
-nu = double(h5readatt(resfile,'/data/input','nu'));
-shear = double(h5readatt(resfile,'/data/input','magnetic shear'));
-
-RTi= double(h5readatt(resfile,'/data/input','RTi'));
-RTe= double(h5readatt(resfile,'/data/input','RTe'));
-Rn= double(h5readatt(resfile,'/data/input','Rni'));
-
-% read electrostatic pot.
-if(nargin == 1 )
- % plot end of the simulation
- if( Nt >1)
- [~,idxte] = min(abs(coordtime -coordtime(end-1)));
- else
- idxte =0;
- end
-elseif(nargin ==2)
- time2plot = varargin{1};
- [~,idxte] = min(abs(coordtime -time2plot));
-end
-
-iframe = idxte-1;
-
-dataset = ['/data/var2d/phi/',num2str(iframe,'%06d')];
-phi = h5read(resfile,dataset);
-
-% read eigenvalu
-% dataset = ['/data/var2d/eig/growth_rate'];
-% growth_rate = h5read(resfile,dataset);
-% dataset = ['/data/var2d/eig/freq'];
-% frequency = h5read(resfile,dataset);
-
-dataset = ['/data/var2d/time'];
-time2d= h5read(resfile,dataset);
-
-% Apply baollooning tranform
-for iky=1:Nky
- dims = size(phi.real);
- phib_real = zeros( dims(1)*Nz ,1);
- phib_imag= phib_real;
- b_angle = phib_real;
-
- midpoint = floor((dims(1)*Nz )/2)+1;
- for ip =1: dims(1)
- start_ = (ip -1)*Nz +1;
- end_ = ip*Nz;
- phib_real(start_:end_) = phi.real(ip,iky,:);
- phib_imag(start_:end_) = phi.imaginary(ip,iky, :);
- end
-
- % Define ballooning angle
- idx = -Nkx:1:Nkx;
- for ip =1: dims(1)
- for iz=1:Nz
- ii = Nz*(ip -1) + iz;
- b_angle(ii) =coordz(iz) + 2*pi*idx(ip);
- end
- end
-
- % normalize real and imaginary parts at chi =0
- [~,idxLFS] = min(abs(b_angle -0));
- phib = phib_real(:) + 1i * phib_imag(:);
- % Normalize to the outermid plane
- phib_norm = phib(:);% / phib( idxLFS) ;
- phib_real_norm(:) = real( phib_norm);%phib_real(:)/phib_real(idxLFS);
- phib_imag_norm(:) = imag( phib_norm);%phib_imag(:)/ phib_imag(idxLFS);
-%
-%
-% figure; hold all;
-% plot(b_angle/pi, phib_real_norm,'-b');
-% plot(b_angle/pi, phib_imag_norm ,'-r');
-% plot(b_angle/pi, sqrt(phib_real_norm .^2 + phib_imag_norm.^2),'--k');
-% legend('real','imag','norm')
-% xlabel('$\chi / \pi$')
-% ylabel('$\phi_B (\chi)$');
-% % title(['$(P,J) =(',num2str(pmax),', ', num2str(jmax),'$), $\nu =',num2str(nu),...
-% % '$, $\epsilon = ',num2str(epsilon),'$, $k_y = ', num2str(coordky( iky)),'$, $q =',num2str(safety_fac),'$, $s = ', num2str(shear),'$, $R_N = ', ...
-% % num2str(Rn),'$, $R_{T_i} = ', num2str(RTi),'$, $N_z =',num2str(Nz),'$']);
-% title(['molix, t=',num2str(coordtime(idxte))]);
-% %set(gca,'Yscale','log')
- %
-
-% %Check symmetry of the mode at the outter mid plane
-% figure; hold all;
-% right = phib_real(midpoint+1:end);
-% left = fliplr(phib_real(1:midpoint-1)');
-% up_down_symm = right(1:end) - left(1:end-1)';
-% %plot(b_angle(midpoint+1:end)/pi,phib_real(midpoint+1:end),'-xb');
-% plot(b_angle(midpoint+1:end)/pi,up_down_symm ,'-xb');
- %plot(abs(b_angle(1:midpoint-1))/pi,phib_real(1:midpoint-1),'-xb');
- %
- %
- % figure; hold all
- % plot(b_angle/pi, phib_imag.^2 + phib_real.^2 ,'xk');
- % %set(gca,'Yscale','log')
- % xlabel('$\chi / \pi$')
- % ylabel('$\phi_B (\chi)$');
- % title(['$(P,J) =(',num2str(pmax),', ', num2str(jmax),'$), $\nu =',num2str(nu),'$, $\epsilon = ',num2str(epsilon),'$, $q =',num2str(safety_fac),'$, $s = ', num2str(shear),'$, $k_y =',num2str(ky),'$']);
-end
-
-end
\ No newline at end of file
diff --git a/matlab/new_flux_results.m b/matlab/new_flux_results.m
deleted file mode 100644
index c52bec38b39e551769904dd2dbb83826ca9b78f3..0000000000000000000000000000000000000000
--- a/matlab/new_flux_results.m
+++ /dev/null
@@ -1,31 +0,0 @@
-%% eta = 0.6
-%nu Gammainf mu
-SGGK_transport = [...
- 1.0e+0, 5.6e-1, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e+00_SGGK_mu_0e+00/
- 5.0e-1, 4.0e-1, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_0e+00/ before wiping ZF
- 5.0e-1, 5.2e-1, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_5e-01_SGGK_mu_0e+00/ after wiping ZF
- 2.5e-1, 3.6e-1, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_2e-01_SGGK_mu_0e+00/
- 1.0e-1, 2.2e-1, 0.0e+0;...
- 1.0e-1, 1.9e-1, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_0e+00
- 7.5e-2, 1.5e-1, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_7e-02_SGGK_mu_0e+00/
- 5.0e-2, 1.1e-1, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_7e-02_SGGK_mu_0e+00/
- 1.0e-2, 4.6e-2, 3.2e-2;...
- ];
-
-DGGK_transport = [...
- 1.0e+0, 5.7e+00, 0.0e+0;... % HeLaZ 2.8 P,J=10,5
- 5.0e-1, 1.1e+01, 0.0e+0;... % simulation_B/cw_DGGK
- 2.5e-1, 6.4e+00, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_2e-01_SGGK_mu_0e+00
- 1.0e-1, 3.0e+00, 0.0e+0;... % HeLaZ 2.8 P,J=2,1
- 1.0e-1, 2.5e+00, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_0e+00
- 7.5e-2, 1.8e+00, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_7e-02_SGGK_mu_0e+00
- 1.0e-2, 3.3e-01, 0.0e+0;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_mu_5e-04
- 1.0e-2, 2.2e-01, 1.0e-4;... % HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_mu_3e-03
- 1.0e-2, 1.4e-01, 3.0e-3;... % HeLaZ 2.8 P,J=6,3
- ];
-
-figure;
-semilogy(SGGK_transport(:,1),SGGK_transport(:,2),'.','MarkerSize',30); hold on;
-semilogy(DGGK_transport(:,1),DGGK_transport(:,2),'.','MarkerSize',30);
-grid on; xlabel('$\nu$'); ylabel('$\Gamma_x$');
-legend(['SGGK';'DGGK'])
\ No newline at end of file
diff --git a/matlab/plot/create_film.m b/matlab/plot/create_film.m
index 5f7aa15abd0e5e20dd5a0b701b48f9ce86f6aa92..58b1baef595e674b081bce2c28b24c75f693ab3f 100644
--- a/matlab/plot/create_film.m
+++ b/matlab/plot/create_film.m
@@ -12,7 +12,7 @@ switch OPTIONS.PLAN
case 'RZ'
toplot = poloidal_plot(DATA,OPTIONS);
otherwise
- toplot = process_field_v2(DATA,OPTIONS);
+ toplot = process_field(DATA,OPTIONS);
end
%%
FILENAME = [DATA.localdir,toplot.FILENAME,format];
diff --git a/matlab/plot/photomaton.m b/matlab/plot/photomaton.m
index a844fafe64824c2347266edf3bc03254ca8997d0..fdcad8f2065ddc7bf6f0c381cee1de77d5584394 100644
--- a/matlab/plot/photomaton.m
+++ b/matlab/plot/photomaton.m
@@ -5,7 +5,7 @@ switch OPTIONS.PLAN
case 'RZ'
toplot = poloidal_plot(DATA,OPTIONS);
otherwise
- toplot = process_field_v2(DATA,OPTIONS);
+ toplot = process_field(DATA,OPTIONS);
end
FNAME = toplot.FILENAME;
FRAMES = toplot.FRAMES;
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index 58305dae1f943add87b63787eae14dbfcd4ca686..84dacb445fbe4bd7e70b278d31302ae1c9c24991 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -64,7 +64,7 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
OPTIONS.COMP = 'avg';
OPTIONS.TIME = DATA.Ts3D;
OPTIONS.POLARPLOT = 0;
- toplot = process_field_v2(DATA,OPTIONS);
+ toplot = process_field(DATA,OPTIONS);
f2plot = toplot.FIELD;
dframe = ite3D - its3D;
clim = max(max(max(abs(plt(f2plot(:,:,:))))));
diff --git a/matlab/plot/show_moments_spectrum.m b/matlab/plot/show_moments_spectrum.m
index d4d71fccd12d133015d88c0d4da74ccce9fb0964..ba4f4008f15091ebb69968937134aabf524ea155 100644
--- a/matlab/plot/show_moments_spectrum.m
+++ b/matlab/plot/show_moments_spectrum.m
@@ -11,16 +11,16 @@ if OPTIONS.ST == 0
end
if OPTIONS.LOGSCALE
logname = 'log';
- compress = @(x,ia) squeeze(log(sum(abs(x(ia,:,:,:,:)),4)));
+ compress = @(x,ia) (log(sum(abs(x(ia,:,:,:,:)),4)));
% compress = @(x,ia) log(sum(abs(squeeze(x(:,:,:,:))),3));
else
logname = '';
- compress = @(x,ia) squeeze(sum(abs(x(ia,:,:,:,:)),4));
+ compress = @(x,ia) (sum(abs(x(ia,:,:,:,:)),4));
% compress = @(x,ia) sum(abs(squeeze(x(:,:,:,:))),3);
end
for ia = 1:DATA.inputs.Na
-% Napjz = sum(abs(squeeze(DATA.Napjz(ia,:,:,:,:))),3);
- Napjz =compress(DATA.Napjz,ia);
+ Napjz = compress(DATA.Napjz,ia);
+ Napjz = reshape(Napjz,DATA.grids.Np,DATA.grids.Nj,numel(DATA.Ts3D));
subplot(double(DATA.inputs.Na),1,double(ia))
plotname = [logname,'$\langle\sum_k |N_',species_name{ia},'^{pj}|\rangle_{p+2j=const}$'];
%We order the moments (0,0) (1,0) (2,0) (0,1) (3,0) (1,1) (4,0) (2,1) (0,2) etc.
@@ -60,12 +60,22 @@ for ia = 1:DATA.inputs.Na
for i_ = 1:Nmoments
Na_ST(i_,:) = Napjz((HL_deg(i_,1)/DATA.grids.deltap)+1,HL_deg(i_,2)+1,:);
end
+
+ if OPTIONS.FILTER
+ % Experimental filtering
+ nt_fourth = ceil(numel(Time_)/4);
+ nt_half = ceil(numel(Time_)/2);
+ Na_ST_avg = mean(Na_ST(:,nt_fourth:nt_half),2);
+ Na_ST = (Na_ST - Na_ST_avg)./Na_ST;
+ OPTIONS.NORMALIZED = 0;
+ %
+ end
% plots
% scaling
if OPTIONS.NORMALIZED
- plt = @(x,i) squeeze(x(i,:)./max(x(i,:)));
+ plt = @(x,i) reshape(x(i,:)./max(x(i,:)),numel(p2j),numel(Time_));
else
- plt = @(x,i) squeeze(x(i,:));
+ plt = @(x,i) reshape(x(i,:),numel(p2j),numel(Time_));
end
if OPTIONS.ST
imagesc(Time_,p2j,plt(Na_ST,1:numel(p2j)));
@@ -75,6 +85,11 @@ for ia = 1:DATA.inputs.Na
yticks(p2j);
yticklabels(ticks_labels)
end
+ if OPTIONS.FILTER
+ caxis([0 0.2]);
+ title('Relative fluctuation from the average');
+ end
+ colorbar
else
colors_ = jet(numel(p2j));
for i = 1:numel(p2j)
@@ -83,7 +98,7 @@ for ia = 1:DATA.inputs.Na
'color',colors_(i,:)); hold on;
end
title(plotname)
-end
+ end
suptitle(DATA.paramshort)
end
diff --git a/matlab/plot/statistical_transport_averaging.m b/matlab/plot/statistical_transport_averaging.m
index cf851f0960c4b635eca357fc4cc1f09475cb7c2c..7207ac6002ec25ba3e12fe5a2e41cdd4e97c0fe2 100644
--- a/matlab/plot/statistical_transport_averaging.m
+++ b/matlab/plot/statistical_transport_averaging.m
@@ -50,6 +50,6 @@ res.Gx_std = std (gamma);
disp(['G_x=',sprintf('%2.2e',res.Gx_avg),'+-',sprintf('%2.2e',res.Gx_std)]);
res.Qx_avg = mean(Qx);
res.Qx_std = std (Qx);
-disp(['G_x=',sprintf('%2.2e',res.Qx_avg),'+-',sprintf('%2.2e',res.Qx_std)]);
+disp(['Q_x=',sprintf('%2.2e',res.Qx_avg),'+-',sprintf('%2.2e',res.Qx_std)]);
end
diff --git a/matlab/post_processing.m b/matlab/post_processing.m
deleted file mode 100644
index 1d051ecea14eb3888dd6e8ce1291623a62bb08ea..0000000000000000000000000000000000000000
--- a/matlab/post_processing.m
+++ /dev/null
@@ -1,205 +0,0 @@
-%% Retrieving max polynomial degree and sampling info
-Npe = numel(Pe); Nje = numel(Je); [JE,PE] = meshgrid(Je,Pe);
-Npi = numel(Pi); Nji = numel(Ji); [JI,PI] = meshgrid(Ji,Pi);
-Ns5D = numel(Ts5D);
-Ns3D = numel(Ts3D);
-% renaming and reshaping quantity of interest
-Ts5D = Ts5D';
-Ts3D = Ts3D';
-
-%% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-disp('Analysis :')
-disp('- iFFT')
-% IFFT (Lower case = real space, upper case = frequency space)
-ne00 = zeros(Nx,Ny,Nz,Ns3D); % Gyrocenter density
-ni00 = zeros(Nx,Ny,Nz,Ns3D); % "
-phi = zeros(Nx,Ny,Nz,Ns3D); % Electrostatic potential
-Z_phi = zeros(Nx,Ny,Nz,Ns3D); % Zonal "
-dens_e = zeros(Nx,Ny,Nz,Ns3D); % Particle density
-dens_i = zeros(Nx,Ny,Nz,Ns3D); %"
-Z_n_e = zeros(Nx,Ny,Nz,Ns3D); % Zonal "
-Z_n_i = zeros(Nx,Ny,Nz,Ns3D); %"
-temp_e = zeros(Nx,Ny,Nz,Ns3D); % Temperature
-temp_i = zeros(Nx,Ny,Nz,Ns3D); % "
-Z_T_e = zeros(Nx,Ny,Nz,Ns3D); % Zonal "
-Z_T_i = zeros(Nx,Ny,Nz,Ns3D); %"
-dyTe = zeros(Nx,Ny,Nz,Ns3D); % Various derivatives
-dyTi = zeros(Nx,Ny,Nz,Ns3D); % "
-dyni = zeros(Nx,Ny,Nz,Ns3D); % "
-dxphi = zeros(Nx,Ny,Nz,Ns3D); % "
-dyphi = zeros(Nx,Ny,Nz,Ns3D); % "
-dx2phi = zeros(Nx,Ny,Nz,Ns3D); % "
-
-for it = 1:numel(Ts3D)
- for iz = 1:numel(z)
- if KIN_E
- NE_ = Ne00(:,:,iz,it);
- ne00 (:,:,iz,it) = real(fftshift(ifft2((NE_),Nx,Ny)));
- if(W_DENS)
- DENS_E_ = DENS_E(:,:,iz,it);
- dens_e (:,:,iz,it) = real(fftshift(ifft2((DENS_E_),Nx,Ny)));
- Z_n_e (:,:,iz,it) = real(fftshift(ifft2((DENS_E_.*(KY==0)),Nx,Ny)));
- end
- if(W_TEMP)
- TEMP_E_ = TEMP_E(:,:,iz,it);
- dyTe(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(TEMP_E_),Nx,Ny)));
- temp_e (:,:,iz,it) = real(fftshift(ifft2((TEMP_E_),Nx,Ny)));
- Z_T_e (:,:,iz,it) = real(fftshift(ifft2((TEMP_E_.*(KY==0)),Nx,Ny)));
- end
- end
- NI_ = Ni00(:,:,iz,it); PH_ = PHI(:,:,iz,it);
- ni00 (:,:,iz,it) = real(fftshift(ifft2((NI_),Nx,Ny)));
- phi (:,:,iz,it) = real(fftshift(ifft2((PH_),Nx,Ny)));
- Z_phi (:,:,iz,it) = real(fftshift(ifft2((PH_.*(KY==0)),Nx,Ny)));
- dxphi (:,:,iz,it) = real(fftshift(ifft2(1i*KX.*(PH_),Nx,Ny)));
- dx2phi(:,:,iz,it) = real(fftshift(ifft2(-KX.^2.*(PH_),Nx,Ny)));
- dyphi (:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(PH_),Nx,Ny)));
- if(W_DENS)
- DENS_I_ = DENS_I(:,:,iz,it);
- dyni (:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(DENS_I_),Nx,Ny)));
- dens_i (:,:,iz,it) = real(fftshift(ifft2((DENS_I_),Nx,Ny)));
- Z_n_i (:,:,iz,it) = real(fftshift(ifft2((DENS_I_.*(KY==0)),Nx,Ny)));
- end
- if(W_TEMP)
- TEMP_I_ = TEMP_I(:,:,iz,it);
- dyTi(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(TEMP_I_),Nx,Ny)));
- temp_i (:,:,iz,it) = real(fftshift(ifft2((TEMP_I_),Nx,Ny)));
- Z_T_i (:,:,iz,it) = real(fftshift(ifft2((TEMP_I_.*(KY==0)),Nx,Ny)));
- end
- end
-end
-
-
-% Post processing
-disp('- post processing')
-% particle flux
-Gamma_x= zeros(Nx,Ny,Nz,Ns3D); % Radial particle transport
-Gamma_y= zeros(Nx,Ny,Nz,Ns3D); % Azymuthal particle transport
-
-% heat flux
-Q_x = zeros(Nx,Ny,Nz,Ns3D);
-Q_y = zeros(Nx,Ny,Nz,Ns3D);
-
-% Epot averages
-phi_maxx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of phi
-phi_avgx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of phi
-phi_maxx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of phi
-phi_avgx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of phi
-
-shear_maxx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of shear
-shear_avgx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of shear
-shear_maxx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of shear
-shear_avgx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of shear
-
-Ne_norm = zeros(Pe_max,Je_max,Ns5D); % Time evol. of the norm of Napj
-Ni_norm = zeros(Pi_max,Ji_max,Ns5D); % .
-
-% Kperp spectrum interpolation
-%full kperp points
-kperp = reshape(sqrt(KX.^2+KY.^2),[numel(KX),1]);
-% interpolated kperps
-nk_noAA = floor(2/3*numel(kx));
-kp_ip = kx;
-[thg, rg] = meshgrid(linspace(0,pi,2*nk_noAA),kp_ip);
-[xn,yn] = pol2cart(thg,rg);
-[ky_s, sortIdx] = sort(ky);
-[xc,yc] = meshgrid(ky_s,kx);
-phi_kp_t = zeros(numel(kp_ip),Nz,Ns3D);
-%
-for it = 1:numel(Ts3D) % Loop over 2D aX_XYays
- for iz = 1:numel(z)
- phi_maxx_maxy(iz,it) = max( max(squeeze(phi(:,:,iz,it))));
- phi_avgx_maxy(iz,it) = max(mean(squeeze(phi(:,:,iz,it))));
- phi_maxx_avgy(iz,it) = mean( max(squeeze(phi(:,:,iz,it))));
- phi_avgx_avgy(iz,it) = mean(mean(squeeze(phi(:,:,iz,it))));
-
- if(W_DENS)
- Gamma_x(:,:,iz,it) = -dens_i(:,:,iz,it).*dyphi(:,:,iz,it);
- Gamma_y(:,:,iz,it) = dens_i(:,:,iz,it).*dxphi(:,:,iz,it);
- end
-
- if(W_TEMP)
- Q_x(:,:,iz,it) = -temp_e(:,:,iz,it).*dyphi(:,:,iz,it);
- Q_y(:,:,iz,it) = temp_i(:,:,iz,it).*dxphi(:,:,iz,it);
- end
-
- shear_maxx_maxy(iz,it) = max( max(squeeze(-(dx2phi(:,:,iz,it)))));
- shear_avgx_maxy(iz,it) = max(mean(squeeze(-(dx2phi(:,:,iz,it)))));
- shear_maxx_avgy(iz,it) = mean( max(squeeze(-(dx2phi(:,:,iz,it)))));
- shear_avgx_avgy(iz,it) = mean(mean(squeeze(-(dx2phi(:,:,iz,it)))));
-
- Z_rth = interp2(xc,yc,squeeze(mean((abs(PHI(:,sortIdx,iz,it))).^2,3)),xn,yn);
- phi_kp_t(:,iz,it) = mean(Z_rth,2);
- end
-end
-%
-for it = 1:numel(Ts5D) % Loop over 5D aX_XYays
- [~, it2D] = min(abs(Ts3D-Ts5D(it)));
- if KIN_E
- Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4)/Nkx/Nky;
- end
- Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4)/Nkx/Nky;
-end
-
-%% Compute primary instability growth rate
-% disp('- growth rate')
-% % Find max value of transport (end of linear mode)
-% [tmp,tmax] = max(GGAMMA_RI*(2*pi/Nx/Ny)^2);
-% [~,itmax] = min(abs(Ts3D-tmax));
-% tstart = 0.1 * Ts3D(itmax); tend = 0.5 * Ts3D(itmax);
-% [~,its3D_lin] = min(abs(Ts3D-tstart));
-% [~,ite3D_lin] = min(abs(Ts3D-tend));
-%
-% g_I = zeros(Nkx,Nky,Nz);
-% for ikx = 1:Nkx
-% for iky = 1:Nky
-% for iz = 1:Nz
-% [g_I(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin)',squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin))));
-% end
-% end
-% end
-% [gmax_I,ikmax_I] = max(max(g_I(1,:,:),[],2),[],3);
-% kmax_I = abs(ky(ikmax_I));
-% Bohm_transport = K_N*gmax_I/kmax_I^2;
-
-%% Compute secondary instability growth rate
-% disp('- growth rate')
-% % Find max value of transport (end of linear mode)
-% % [tmp,tmax] = max(GGAMMA_RI*(2*pi/Nx/Ny)^2);
-% % [~,itmax] = min(abs(Ts2D-tmax));
-% % tstart = Ts2D(itmax); tend = 1.5*Ts2D(itmax);
-% [~,its3D_lin] = min(abs(Ts3D-tstart));
-% [~,ite3D_lin] = min(abs(Ts3D-tend));
-%
-% g_II = zeros(Nkx,Nky);
-% for ikx = 1:Nkx
-% for iky = 1
-% for iz = 1:Nz
-% [g_II(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin)',squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin))));
-% end
-% end
-% end
-% [gmax_II,ikmax_II] = max(max(g_II(1,:,:),[],2),[],3);
-% kmax_II = abs(kx(ikmax_II));
-
-%% zonal vs nonzonal energies for phi(t)
-Ephi_Z = zeros(1,Ns3D);
-Ephi_NZ_kgt0 = zeros(1,Ns3D);
-Ephi_NZ_kgt1 = zeros(1,Ns3D);
-Ephi_NZ_kgt2 = zeros(1,Ns3D);
-high_k_phi = zeros(1,Ns3D);
-for it = 1:numel(Ts3D)
-% Ephi_NZ(it) = sum(sum(((KY~=0).*abs(PHI(:,:,1,it)).^2)));
-% Ephi_Z(it) = sum(sum(((KY==0).*abs(PHI(:,:,1,it)).^2)));
- [amp,ikzf] = max(abs((kx~=0).*PHI(:,1,1,it)));
-% Ephi_NZ(it) = sum(sum(((KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)));
- Ephi_NZ_kgt0(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>0.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2))));
- Ephi_NZ_kgt1(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>1.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2))));
- Ephi_NZ_kgt2(it) = squeeze(sum(sum(((sqrt(KX.^2+KY.^2)>2.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2))));
-% Ephi_Z(it) = kx(ikzf)^2*abs(PHI(ikzf,1,1,it)).^2;
- Ephi_Z(it) = squeeze(sum(sum(((KX~=0).*(KY==0).*(KX.^2).*abs(PHI(:,:,1,it)).^2))));
-% Ephi_NZ(it) = sum(sum(((KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)))-Ephi_Z(it);
-% high_k_phi(it) = squeeze(abs(PHI(18,18,1,it)).^2); % kperp sqrt(2)
-% high_k_phi(it) = abs(PHI(40,40,1,it)).^2;% kperp 3.5
-
-end
diff --git a/matlab/process_field.m b/matlab/process_field.m
index baf1ab809c2be250c9d65234ddc048764ade15b3..5b2d103e3e38e158a342115b942ac53e546671ee 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -9,35 +9,35 @@ for i = 1:numel(OPTIONS.TIME)
end
FRAMES = unique(FRAMES);
%% Setup the plot geometry
-[~,KY] = meshgrid(DATA.kx,DATA.ky);
-directions = {'x','y','z'};
-Nx = DATA.Nx; Ny = DATA.Ny; Nz = DATA.Nz; Nt = numel(FRAMES);
+[KX, KY] = meshgrid(DATA.grids.kx, DATA.grids.ky);
+directions = {'y','x','z'};
+Nx = DATA.grids.Nx; Ny = DATA.grids.Ny; Nz = DATA.grids.Nz; Nt = numel(FRAMES);
POLARPLOT = OPTIONS.POLARPLOT;
LTXNAME = OPTIONS.NAME;
switch OPTIONS.PLAN
case 'xy'
XNAME = '$x$'; YNAME = '$y$';
- [X,Y] = meshgrid(DATA.x,DATA.y);
+ [X,Y] = meshgrid(DATA.grids.x,DATA.grids.y);
REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 1;
case 'xz'
XNAME = '$x$'; YNAME = '$z$';
- [Y,X] = meshgrid(DATA.z,DATA.x);
+ [Y,X] = meshgrid(DATA.grids.z,DATA.grids.x);
REALP = 1; COMPDIM = 1; SCALE = 0;
case 'yz'
XNAME = '$y$'; YNAME = '$z$';
- [Y,X] = meshgrid(DATA.z,DATA.y);
+ [Y,X] = meshgrid(DATA.grids.z,DATA.grids.y);
REALP = 1; COMPDIM = 2; SCALE = 0;
case 'kxky'
XNAME = '$k_x$'; YNAME = '$k_y$';
- [X,Y] = meshgrid(DATA.kx,DATA.ky);
+ [X,Y] = meshgrid(DATA.grids.kx,DATA.grids.ky);
REALP = 0; COMPDIM = 3; POLARPLOT = 0; SCALE = 1;
case 'kxz'
XNAME = '$k_x$'; YNAME = '$z$';
- [Y,X] = meshgrid(DATA.z,DATA.kx);
+ [Y,X] = meshgrid(DATA.grids.z,DATA.grids.kx);
REALP = 0; COMPDIM = 1; POLARPLOT = 0; SCALE = 0;
case 'kyz'
XNAME = '$k_y$'; YNAME = '$z$';
- [Y,X] = meshgrid(DATA.z,DATA.ky);
+ [Y,X] = meshgrid(DATA.grids.z,DATA.grids.ky);
REALP = 0; COMPDIM = 2; POLARPLOT = 0; SCALE = 0;
case 'sx'
XNAME = '$v_\parallel$'; YNAME = '$\mu$';
@@ -46,7 +46,7 @@ switch OPTIONS.PLAN
for i = 1:numel(OPTIONS.TIME)
[~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts5D));
end
- FRAMES = unique(FRAMES);
+ FRAMES = unique(FRAMES); Nt = numel(FRAMES);
end
Xmax_ = max(max(abs(X))); Ymax_ = max(max(abs(Y)));
Lmin_ = min([Xmax_,Ymax_]);
@@ -73,7 +73,30 @@ else
FIELD = zeros(sz(1),sz(2),Nt);
end
%% Process the field to plot
-
+% short term writing --
+% b_e = DATA.sigma_e*sqrt(2*DATA.tau_e)*sqrt(KX.^2+KY.^2);
+% adiab_e = kernel(0,b_e);
+% pol_e = kernel(0,b_e).^2;
+% for n = 1:DATA.Jmaxe
+% adiab_e = adiab_e + kernel(n,b_e).^2;
+% pol_e = pol_e + kernel(n,b_e).^2;
+% end
+% adiab_e = DATA.q_e/DATA.tau_e .* adiab_e;
+% pol_e = DATA.q_e^2/DATA.tau_e * (1 - pol_e);
+%
+% b_i = DATA.sigma_i*sqrt(2*DATA.tau_i)*sqrt(KX.^2+KY.^2);
+% adiab_i = kernel(0,b_i);
+% pol_i = kernel(0,b_i).^2;
+% for n = 1:DATA.Jmaxi
+% adiab_i = adiab_i + kernel(n,b_i).^2;
+% pol_i = pol_i + kernel(n,b_i).^2;
+% end
+% pol_i = DATA.q_i^2/DATA.tau_i * (1 - pol_i);
+% adiab_i = DATA.q_i/DATA.tau_i .* adiab_i;
+% poisson_op = (pol_i + pol_e);
+adiab_e =0; adiab_i =0; poisson_op=1;
+SKIP_COMP = 0; % Turned on only for kin. distr. func. plot
+% --
switch OPTIONS.COMP
case 'avg'
compr = @(x) mean(x,COMPDIM);
@@ -89,24 +112,24 @@ switch OPTIONS.COMP
i = OPTIONS.COMP;
compr = @(x) x(i,:,:);
if REALP
- COMPNAME = sprintf(['y=','%2.1f'],DATA.x(i));
+ COMPNAME = sprintf(['y=','%2.1f'],DATA.grids.x(i));
else
- COMPNAME = sprintf(['k_y=','%2.1f'],DATA.kx(i));
+ COMPNAME = sprintf(['k_y=','%2.1f'],DATA.grids.kx(i));
end
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
case 2
i = OPTIONS.COMP;
compr = @(x) x(:,i,:);
if REALP
- COMPNAME = sprintf(['x=','%2.1f'],DATA.y(i));
+ COMPNAME = sprintf(['x=','%2.1f'],DATA.grids.y(i));
else
- COMPNAME = sprintf(['k_x=','%2.1f'],DATA.ky(i));
+ COMPNAME = sprintf(['k_x=','%2.1f'],DATA.grids.ky(i));
end
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
case 3
i = OPTIONS.COMP;
compr = @(x) x(:,:,i);
- COMPNAME = sprintf(['z=','%2.1f'],DATA.z(i));
+ COMPNAME = sprintf(['z=','%2.1f'],DATA.grids.z(i));
FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
end
end
@@ -136,384 +159,113 @@ switch REALP
end
switch OPTIONS.NAME
- case '\phi'
+ case '\phi' %ES pot
NAME = 'phi';
- if COMPDIM == 3
- for it = 1:numel(FRAMES)
- tmp = squeeze(compr(DATA.PHI(:,:,:,FRAMES(it))));
- FIELD(:,:,it) = squeeze(process(tmp));
- end
- else
- if REALP
- tmp = zeros(Ny,Nx,Nz);
- else
- tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
- end
- for it = 1:numel(FRAMES)
- for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(DATA.PHI(:,:,iz,FRAMES(it))));
- end
- FIELD(:,:,it) = squeeze(compr(tmp));
- end
- end
- case '\psi'
+ FLD_ = DATA.PHI(:,:,:,FRAMES); % data to plot
+ OPE_ = 1; % Operation on data
+ case '\psi' %EM pot
NAME = 'psi';
- if COMPDIM == 3
- for it = 1:numel(FRAMES)
- tmp = squeeze(compr(DATA.PSI(:,:,:,FRAMES(it))));
- FIELD(:,:,it) = squeeze(process(tmp));
- end
- else
- if REALP
- tmp = zeros(Ny,Nx,Nz);
- else
- tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
- end
- for it = 1:numel(FRAMES)
- for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(DATA.PSI(:,:,iz,FRAMES(it))));
- end
- FIELD(:,:,it) = squeeze(compr(tmp));
- end
- end
- case 'N_i^{00}'
+ FLD_ = DATA.PSI(:,:,:,FRAMES);
+ OPE_ = 1;
+ case '\phi^{NZ}' % non-zonal ES pot
+ NAME = 'phiNZ';
+ FLD_ = DATA.PHI(:,:,:,FRAMES);
+ OPE_ = (KY~=0);
+ case 'v_{Ey}' % y-comp of ExB velocity
+ NAME = 'vy';
+ FLD_ = DATA.PHI(:,:,:,FRAMES);
+ OPE_ = -1i*KX;
+ case 'v_{Ex}' % x-comp of ExB velocity
+ NAME = 'vx';
+ FLD_ = DATA.PHI(:,:,:,FRAMES);
+ OPE_ = -1i*KY;
+ case 's_{Ey}' % y-comp of ExB shear
+ NAME = 'sy';
+ FLD_ = DATA.PHI(:,:,:,FRAMES);
+ OPE_ = KX.^2;
+ case 's_{Ex}' % x-comp of ExB shear
+ NAME = 'sx';
+ FLD_ = DATA.PHI(:,:,:,FRAMES);
+ OPE_ = KY.^2;
+ case '\omega_z' % ES pot vorticity
+ NAME = 'vorticity';
+ FLD_ = DATA.PHI(:,:,:,FRAMES);
+ OPE_ = -(KX.^2+KY.^2);
+ case 'N_i^{00}' % first ion gyromoment
NAME = 'Ni00';
- if COMPDIM == 3
- for it = 1:numel(FRAMES)
- tmp = squeeze(compr(DATA.Ni00(:,:,:,FRAMES(it))));
- FIELD(:,:,it) = squeeze(process(tmp));
- end
- else
- if REALP
- tmp = zeros(Ny,Nx,Nz);
- else
- tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
- end
- for it = 1:numel(FRAMES)
- for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(DATA.Ni00(:,:,iz,FRAMES(it))));
- end
- FIELD(:,:,it) = squeeze(compr(tmp));
- end
- end
- case 'N_e^{00}'
+ FLD_ = DATA.Ni00(:,:,:,FRAMES);
+ OPE_ = 1;
+ case 'N_e^{00}' % first electron gyromoment
NAME = 'Ne00';
- if COMPDIM == 3
- for it = 1:numel(FRAMES)
- tmp = squeeze(compr(DATA.Ne00(:,:,:,FRAMES(it))));
- FIELD(:,:,it) = squeeze(process(tmp));
- end
- else
- if REALP
- tmp = zeros(Ny,Nx,Nz);
- else
- tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
- end
- for it = 1:numel(FRAMES)
- for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(DATA.Ne00(:,:,iz,FRAMES(it))));
- end
- FIELD(:,:,it) = squeeze(compr(tmp));
- end
- end
- case 'n_e'
+ FLD_ = DATA.Ne00(:,:,:,FRAMES);
+ OPE_ = 1;
+ case 'N_i^{00}-N_e^{00}' % first electron gyromoment
+ NAME = 'Ni00-Ne00';
+ FLD_ = (DATA.Ni00(:,:,:,FRAMES)-DATA.Ne00(:,:,:,FRAMES))./(poisson_op+(poisson_op==0));
+ OPE_ = 1;
+ case 'n_i' % ion density
+ NAME = 'ni';
+ FLD_ = DATA.DENS_I(:,:,:,FRAMES) - adiab_i.* DATA.PHI(:,:,:,FRAMES);
+ OPE_ = 1;
+ case 'n_e' % electron density
NAME = 'ne';
- if COMPDIM == 3
- for it = 1:numel(FRAMES)
- tmp = squeeze(compr(DATA.DENS_E(:,:,:,FRAMES(it))));
- FIELD(:,:,it) = squeeze(process(tmp));
- end
- else
- if REALP
- tmp = zeros(Ny,Nx,Nz);
- else
- tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
- end
- for it = 1:numel(FRAMES)
- for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(DATA.DENS_E(:,:,iz,FRAMES(it))));
- end
- FIELD(:,:,it) = squeeze(compr(tmp));
- end
- end
- case 'k^2n_e'
+ FLD_ = DATA.DENS_E(:,:,:,FRAMES) - adiab_e.* DATA.PHI(:,:,:,FRAMES);
+ OPE_ = 1;
+ case 'k^2n_e' % electron vorticity
NAME = 'k2ne';
- [KX, KY] = meshgrid(DATA.kx, DATA.ky);
- if COMPDIM == 3
- for it = 1:numel(FRAMES)
- for iz = 1:DATA.Nz
- tmp = squeeze(compr(-(KX.^2+KY.^2).*DATA.DENS_E(:,:,iz,FRAMES(it))));
- end
- FIELD(:,:,it) = squeeze(process(tmp));
- end
- else
- if REALP
- tmp = zeros(Ny,Nx,Nz);
- else
- tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
- end
- for it = 1:numel(FRAMES)
- for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(-(KX.^2+KY.^2).*DATA.DENS_E(:,:,iz,FRAMES(it))));
- end
- FIELD(:,:,it) = squeeze(compr(tmp));
- end
- end
- case 'n_e^{NZ}'
- NAME = 'neNZ';
- if COMPDIM == 3
- for it = 1:numel(FRAMES)
- tmp = squeeze(compr(DATA.DENS_E(:,:,:,FRAMES(it)).*(KY~=0)));
- FIELD(:,:,it) = squeeze(process(tmp));
- end
- else
- if REALP
- tmp = zeros(Ny,Nx,Nz);
- else
- tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
- end
- for it = 1:numel(FRAMES)
- for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(DATA.DENS_E(:,:,iz,FRAMES(it)).*(KY~=0)));
- end
- FIELD(:,:,it) = squeeze(compr(tmp));
- end
- end
- case 'n_i'
- NAME = 'ni';
- if COMPDIM == 3
- for it = 1:numel(FRAMES)
- tmp = squeeze(compr(DATA.DENS_I(:,:,:,FRAMES(it))));
- FIELD(:,:,it) = squeeze(process(tmp));
- end
- else
- if REALP
- tmp = zeros(Ny,Nx,Nz);
- else
- tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
- end
- for it = 1:numel(FRAMES)
- for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(DATA.DENS_I(:,:,iz,FRAMES(it))));
- end
- FIELD(:,:,it) = squeeze(compr(tmp));
- end
- end
- case 'T_i'
+ FLD_ = DATA.DENS_E(:,:,:,FRAMES);
+ OPE_ = -(KX.^2+KY.^2);
+ case 'n_i-n_e' % polarisation
+ NAME = 'pol';
+ OPE_ = 1;
+ FLD_ = ((DATA.DENS_I(:,:,:,FRAMES)- adiab_i.* DATA.PHI(:,:,:,FRAMES))...
+ -(DATA.DENS_E(:,:,:,FRAMES)- adiab_e.* DATA.PHI(:,:,:,FRAMES)));
+ case 'T_i' % ion temperature
NAME = 'Ti';
- if COMPDIM == 3
- for it = 1:numel(FRAMES)
- tmp = squeeze(compr(DATA.TEMP_I(:,:,:,FRAMES(it))));
- FIELD(:,:,it) = squeeze(process(tmp));
- end
- else
- if REALP
- tmp = zeros(Ny,Nx,Nz);
- else
- tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
- end
- for it = 1:numel(FRAMES)
- for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(DATA.TEMP_I(:,:,iz,FRAMES(it))));
- end
- FIELD(:,:,it) = squeeze(compr(tmp));
- end
- end
- case 'n_i^{NZ}'
- NAME = 'niNZ';
- if COMPDIM == 3
- for it = 1:numel(FRAMES)
- tmp = squeeze(compr(DATA.DENS_I(:,:,:,FRAMES(it)).*(KY~=0)));
- FIELD(:,:,it) = squeeze(process(tmp));
- end
- else
- if REALP
- tmp = zeros(Ny,Nx,Nz);
- else
- tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
- end
- for it = 1:numel(FRAMES)
- for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(DATA.DENS_I(:,:,iz,FRAMES(it)).*(KY~=0)));
- end
- FIELD(:,:,it) = squeeze(compr(tmp));
- end
- end
- case '\phi^{NZ}'
- NAME = 'phiNZ';
- if COMPDIM == 3
- for it = 1:numel(FRAMES)
- tmp = squeeze(compr(DATA.PHI(:,:,:,FRAMES(it)).*(KY~=0)));
- FIELD(:,:,it) = squeeze(process(tmp));
- end
- else
- if REALP
- tmp = zeros(Ny,Nx,Nz);
- else
- tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
- end
- for it = 1:numel(FRAMES)
- for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(DATA.PHI(:,:,iz,FRAMES(it)).*(KY~=0)));
- end
- FIELD(:,:,it) = squeeze(compr(tmp));
- end
- end
- case 'v_y'
- NAME = 'vy';
- [KX, ~] = meshgrid(DATA.kx, DATA.ky);
- vE = zeros(DATA.Ny,DATA.Nx,DATA.Nz,numel(FRAMES));
- for it = 1:numel(FRAMES)
- for iz = 1:DATA.Nz
- vE(:,:,iz,it) = real(ifft2(-1i*KX.*(DATA.PHI(:,:,iz,FRAMES(it))),DATA.Ny,DATA.Nx));
- end
- end
- if REALP % plot in real space
- for it = 1:numel(FRAMES)
- FIELD(:,:,it) = squeeze(compr(vE(:,:,:,it)));
- end
- else % Plot spectrum
- process = @(x) abs(fftshift(x,2));
- tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
- for it = 1:numel(FRAMES)
- for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(-1i*KX.*(DATA.PHI(:,:,iz,FRAMES(it)))));
- end
- FIELD(:,:,it) = squeeze(compr(tmp));
- end
- end
- case 'v_x'
- NAME = 'vx';
- [~, KY] = meshgrid(DATA.kx, DATA.ky);
- vE = zeros(DATA.Ny,DATA.Nx,DATA.Nz,numel(FRAMES));
- for it = 1:numel(FRAMES)
- for iz = 1:DATA.Nz
- vE(:,:,iz,it) = real(ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,FRAMES(it)))));
- end
- end
- if REALP % plot in real space
- for it = 1:numel(FRAMES)
- FIELD(:,:,it) = squeeze(compr(vE(:,:,:,it)));
- end
- else % Plot spectrum
- process = @(x) abs(fftshift(x,2));
- tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
- for it = 1:numel(FRAMES)
- for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(-1i*KY.*(DATA.PHI(:,:,iz,FRAMES(it)))));
- end
- FIELD(:,:,it) = squeeze(compr(tmp));
- end
- end
- case 's_y'
- NAME = 'sy';
- [KX, ~] = meshgrid(DATA.kx, DATA.ky);
- vE = zeros(DATA.Ny,DATA.Nx,DATA.Nz,numel(FRAMES));
- for it = 1:numel(FRAMES)
- for iz = 1:DATA.Nz
- vE(:,:,iz,it) = real(ifourier_GENE(KX.^2.*(DATA.PHI(:,:,iz,FRAMES(it)))));
- end
- end
- if REALP % plot in real space
- for it = 1:numel(FRAMES)
- FIELD(:,:,it) = squeeze(compr(vE(:,:,:,it)));
- end
- else % Plot spectrum
- process = @(x) abs(fftshift(x,2));
- tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
- for it = 1:FRAMES
- for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(KX.^2.*(DATA.PHI(:,:,iz,FRAMES(it)))));
- end
- FIELD(:,:,it) = squeeze(compr(tmp));
- end
- end
- case '\omega_z'
- NAME = 'vorticity';
- [KX, KY] = meshgrid(DATA.kx, DATA.ky);
- wE = zeros(DATA.Ny,DATA.Nx,DATA.Nz,numel(FRAMES));
+ FLD_ = DATA.TEMP_I(:,:,:,FRAMES);
+ OPE_ = 1;
+ case 'G_x' % ion particle flux
+ NAME = 'Gx';
+ FLD_ = 0.*DATA.PHI(:,:,:,FRAMES);
+ OPE_ = 1;
for it = 1:numel(FRAMES)
+ tmp = zeros(DATA.Ny,DATA.Nx,Nz);
for iz = 1:DATA.Nz
- wE(:,:,iz,it) = real(ifourier_GENE(-(KX.^2+KY.^2).*(DATA.PHI(:,:,iz,FRAMES(it)))));
- end
- end
- if REALP % plot in real space
- for it = 1:numel(FRAMES)
- FIELD(:,:,it) = squeeze(compr(wE(:,:,:,it)));
- end
- else % Plot spectrum
- process = @(x) abs(fftshift(x,2));
- tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
- for it = 1:FRAMES
- for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = squeeze(process(-(KX.^2+KY.^2).*(DATA.PHI(:,:,iz,FRAMES(it)))));
- end
- FIELD(:,:,it) = squeeze(compr(tmp));
- end
- end
- case '\Gamma_x'
- NAME = 'Gx';
- [~, KY] = meshgrid(DATA.kx, DATA.ky);
- Gx = zeros(DATA.Ny,DATA.Nx,DATA.Nz,numel(FRAMES));
+ vx_ = real((ifourier_GENE(-1i*KY.*(DATA.PHI (:,:,iz,FRAMES(it))))));
+ ni_ = real((ifourier_GENE( DATA.DENS_I(:,:,iz,FRAMES(it)))));
+ gx_ = vx_.*ni_;
+% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.Ny,DATA.Nx))),2));
+ tmp(:,:,iz) = abs((squeeze(fft2(gx_,DATA.Ny,DATA.Nx))));
+ end
+ FLD_(:,:,it)= squeeze(compr(tmp(1:DATA.Nky,1:DATA.Nkx,:)));
+ end
+ case 'Q_x' % ion heat flux
+ NAME = 'Qx';
+ FLD_ = 0.*DATA.PHI(:,:,:,FRAMES);
+ OPE_ = 1;
for it = 1:numel(FRAMES)
- for iz = 1:DATA.Nz
- Gx(:,:,iz,it) = real((ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,FRAMES(it))))))...
- .*real((ifourier_GENE(DATA.DENS_I(:,:,iz,FRAMES(it)))));
- end
- end
- if REALP % plot in real space
- for it = 1:numel(FRAMES)
- FIELD(:,:,it) = squeeze(compr(Gx(:,:,:,it)));
- end
- else % Plot spectrum
- process = @(x) abs(fftshift(x,2));
- shift_x = @(x) fftshift(x,2);
- shift_y = @(x) fftshift(x,2);
tmp = zeros(DATA.Ny,DATA.Nx,Nz);
- for it = 1:numel(FRAMES)
- for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = process(squeeze(fft2(Gx(:,:,iz,it),DATA.Ny,DATA.Nx)));
- end
- FIELD(:,:,it) = squeeze(compr(tmp(1:DATA.Nky,1:DATA.Nkx,:)));
- end
- end
- case 'Q_x'
- NAME = 'Qx';
- [~, KY] = meshgrid(DATA.kx, DATA.ky);
- Qx = zeros(DATA.Ny,DATA.Nx,DATA.Nz,numel(FRAMES));
- for it = 1:numel(FRAMES)
for iz = 1:DATA.Nz
- Qx(:,:,iz,it) = real(ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,FRAMES(it)))))...
- .*real(ifourier_GENE(DATA.DENS_I(:,:,iz,FRAMES(it))))...
- .*real(ifourier_GENE(DATA.TEMP_I(:,:,iz,FRAMES(it))));
- end
- end
- if REALP % plot in real space
- for it = 1:numel(FRAMES)
- FIELD(:,:,it) = squeeze(compr(Qx(:,:,:,it)));
- end
- else % Plot spectrum
- process = @(x) abs(fftshift(x,2));
- shift_x = @(x) fftshift(x,2);
- shift_y = @(x) fftshift(x,2);
- tmp = zeros(DATA.Ny,DATA.Nx,Nz);
- for it = 1:numel(FRAMES)
- for iz = 1:numel(DATA.z)
- tmp(:,:,iz) = process(squeeze(fft2(Qx(:,:,iz,it),DATA.Ny,DATA.Nx)));
- end
- FIELD(:,:,it) = squeeze(compr(tmp(1:DATA.Nky,1:DATA.Nkx,:)));
- end
- end
+ vx_ = real((ifourier_GENE(-1i*KY.*(DATA.PHI (:,:,iz,FRAMES(it))))));
+ ni_ = real((ifourier_GENE( DATA.DENS_I(:,:,iz,FRAMES(it)))));
+ Ti_ = real((ifourier_GENE( DATA.TEMP_I(:,:,iz,FRAMES(it)))));
+ qx_ = vx_.*ni_.*Ti_;
+% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.Ny,DATA.Nx))),2));
+ tmp(:,:,iz) = abs((squeeze(fft2(qx_,DATA.Ny,DATA.Nx))));
+ end
+ FLD_(:,:,it)= squeeze(compr(tmp(1:DATA.Nky,1:DATA.Nkx,:)));
+ end
case 'f_i'
+ SKIP_COMP = 1;
shift_x = @(x) x; shift_y = @(y) y;
NAME = 'fi'; OPTIONS.SPECIE = 'i';
- for it = 1:numel(OPTIONS.TIME)
- [~,it0_] =min(abs(OPTIONS.TIME(it)-DATA.Ts5D));
- OPTIONS.T = DATA.Ts5D(it0_);
+ for it = 1:numel(FRAMES)
+ OPTIONS.T = DATA.Ts5D(FRAMES(it));
OPTIONS.Z = OPTIONS.COMP;
[~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS);
end
case 'f_e'
+ SKIP_COMP = 1;
shift_x = @(x) x; shift_y = @(y) y;
NAME = 'fe'; OPTIONS.SPECIE = 'e';
[~,it0_] =min(abs(OPTIONS.TIME(1)-DATA.Ts5D));
@@ -525,6 +277,30 @@ switch OPTIONS.NAME
OPTIONS.T = DATA.Ts5D(FRAMES(it));
[~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS);
end
+ otherwise
+ disp('Fieldname not recognized');
+ return
+end
+% Process the field according to the 2D plane and the space (real/cpx)
+if ~SKIP_COMP
+if COMPDIM == 3
+ for it = 1:numel(FRAMES)
+ tmp = squeeze(compr(OPE_.*FLD_(:,:,:,it)));
+ FIELD(:,:,it) = squeeze(process(tmp));
+ end
+else
+ if REALP
+ tmp = zeros(Ny,Nx,Nz);
+ else
+ tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
+ end
+ for it = 1:numel(FRAMES)
+ for iz = 1:numel(DATA.grids.z)
+ tmp(:,:,iz) = squeeze(process(OPE_.*FLD_(:,:,iz,it)));
+ end
+ FIELD(:,:,it) = squeeze(compr(tmp));
+ end
+end
end
TOPLOT.FIELD = FIELD;
TOPLOT.FRAMES = FRAMES;
diff --git a/matlab/process_field_v2.m b/matlab/process_field_v2.m
deleted file mode 100644
index 5b2d103e3e38e158a342115b942ac53e546671ee..0000000000000000000000000000000000000000
--- a/matlab/process_field_v2.m
+++ /dev/null
@@ -1,318 +0,0 @@
-function [ TOPLOT ] = process_field( DATA, OPTIONS )
-%UNTITLED4 Summary of this function goes here
-% Detailed explanation goes here
-%% Setup time
-%%
-FRAMES = zeros(size(OPTIONS.TIME));
-for i = 1:numel(OPTIONS.TIME)
- [~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts3D));
-end
-FRAMES = unique(FRAMES);
-%% Setup the plot geometry
-[KX, KY] = meshgrid(DATA.grids.kx, DATA.grids.ky);
-directions = {'y','x','z'};
-Nx = DATA.grids.Nx; Ny = DATA.grids.Ny; Nz = DATA.grids.Nz; Nt = numel(FRAMES);
-POLARPLOT = OPTIONS.POLARPLOT;
-LTXNAME = OPTIONS.NAME;
-switch OPTIONS.PLAN
- case 'xy'
- XNAME = '$x$'; YNAME = '$y$';
- [X,Y] = meshgrid(DATA.grids.x,DATA.grids.y);
- REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 1;
- case 'xz'
- XNAME = '$x$'; YNAME = '$z$';
- [Y,X] = meshgrid(DATA.grids.z,DATA.grids.x);
- REALP = 1; COMPDIM = 1; SCALE = 0;
- case 'yz'
- XNAME = '$y$'; YNAME = '$z$';
- [Y,X] = meshgrid(DATA.grids.z,DATA.grids.y);
- REALP = 1; COMPDIM = 2; SCALE = 0;
- case 'kxky'
- XNAME = '$k_x$'; YNAME = '$k_y$';
- [X,Y] = meshgrid(DATA.grids.kx,DATA.grids.ky);
- REALP = 0; COMPDIM = 3; POLARPLOT = 0; SCALE = 1;
- case 'kxz'
- XNAME = '$k_x$'; YNAME = '$z$';
- [Y,X] = meshgrid(DATA.grids.z,DATA.grids.kx);
- REALP = 0; COMPDIM = 1; POLARPLOT = 0; SCALE = 0;
- case 'kyz'
- XNAME = '$k_y$'; YNAME = '$z$';
- [Y,X] = meshgrid(DATA.grids.z,DATA.grids.ky);
- REALP = 0; COMPDIM = 2; POLARPLOT = 0; SCALE = 0;
- case 'sx'
- XNAME = '$v_\parallel$'; YNAME = '$\mu$';
- [Y,X] = meshgrid(OPTIONS.XPERP,OPTIONS.SPAR);
- REALP = 1; COMPDIM = 3; POLARPLOT = 0; SCALE = 0;
- for i = 1:numel(OPTIONS.TIME)
- [~,FRAMES(i)] =min(abs(OPTIONS.TIME(i)-DATA.Ts5D));
- end
- FRAMES = unique(FRAMES); Nt = numel(FRAMES);
-end
-Xmax_ = max(max(abs(X))); Ymax_ = max(max(abs(Y)));
-Lmin_ = min([Xmax_,Ymax_]);
-Rx = Xmax_/Lmin_ * SCALE + (1-SCALE)*1.2;
-Ry = Ymax_/Lmin_ * SCALE + (1-SCALE)*1.2;
-ASPECT = [Rx, Ry, 1];
-DIMENSIONS = [500, 1000, OPTIONS.RESOLUTION*Rx, OPTIONS.RESOLUTION*Ry];
-% Polar grid
-POLARNAME = [];
-if POLARPLOT
- POLARNAME = 'polar';
- X__ = (X+DATA.a).*cos(Y);
- Y__ = (X+DATA.a).*sin(Y);
- X = X__;
- Y = Y__;
- XNAME='X';
- YNAME='Z';
- DIMENSIONS = [100, 100, OPTIONS.RESOLUTION, OPTIONS.RESOLUTION];
- ASPECT = [1,1,1];
- sz = size(X);
- FIELD = zeros(sz(1),sz(2),Nt);
-else
- sz = size(X);
- FIELD = zeros(sz(1),sz(2),Nt);
-end
-%% Process the field to plot
-% short term writing --
-% b_e = DATA.sigma_e*sqrt(2*DATA.tau_e)*sqrt(KX.^2+KY.^2);
-% adiab_e = kernel(0,b_e);
-% pol_e = kernel(0,b_e).^2;
-% for n = 1:DATA.Jmaxe
-% adiab_e = adiab_e + kernel(n,b_e).^2;
-% pol_e = pol_e + kernel(n,b_e).^2;
-% end
-% adiab_e = DATA.q_e/DATA.tau_e .* adiab_e;
-% pol_e = DATA.q_e^2/DATA.tau_e * (1 - pol_e);
-%
-% b_i = DATA.sigma_i*sqrt(2*DATA.tau_i)*sqrt(KX.^2+KY.^2);
-% adiab_i = kernel(0,b_i);
-% pol_i = kernel(0,b_i).^2;
-% for n = 1:DATA.Jmaxi
-% adiab_i = adiab_i + kernel(n,b_i).^2;
-% pol_i = pol_i + kernel(n,b_i).^2;
-% end
-% pol_i = DATA.q_i^2/DATA.tau_i * (1 - pol_i);
-% adiab_i = DATA.q_i/DATA.tau_i .* adiab_i;
-% poisson_op = (pol_i + pol_e);
-adiab_e =0; adiab_i =0; poisson_op=1;
-SKIP_COMP = 0; % Turned on only for kin. distr. func. plot
-% --
-switch OPTIONS.COMP
- case 'avg'
- compr = @(x) mean(x,COMPDIM);
- COMPNAME = ['avg',directions{COMPDIM}];
- FIELDNAME = ['\langle ',LTXNAME,'\rangle_',directions{COMPDIM}];
- case 'max'
- compr = @(x) max(x,[],COMPDIM);
- COMPNAME = ['max',directions{COMPDIM}];
- FIELDNAME = ['\max_',directions{COMPDIM},' ',LTXNAME];
- otherwise
- switch COMPDIM
- case 1
- i = OPTIONS.COMP;
- compr = @(x) x(i,:,:);
- if REALP
- COMPNAME = sprintf(['y=','%2.1f'],DATA.grids.x(i));
- else
- COMPNAME = sprintf(['k_y=','%2.1f'],DATA.grids.kx(i));
- end
- FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
- case 2
- i = OPTIONS.COMP;
- compr = @(x) x(:,i,:);
- if REALP
- COMPNAME = sprintf(['x=','%2.1f'],DATA.grids.y(i));
- else
- COMPNAME = sprintf(['k_x=','%2.1f'],DATA.grids.ky(i));
- end
- FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
- case 3
- i = OPTIONS.COMP;
- compr = @(x) x(:,:,i);
- COMPNAME = sprintf(['z=','%2.1f'],DATA.grids.z(i));
- FIELDNAME = [LTXNAME,'(',COMPNAME,')'];
- end
-end
-
-switch REALP
- case 1 % Real space plot
- INTERP = OPTIONS.INTERP;
- process = @(x) real(fftshift(ifourier_GENE(x)));
- shift_x = @(x) x;
- shift_y = @(x) x;
- case 0 % Frequencies plot
- INTERP = 0;
- switch COMPDIM
- case 1
- process = @(x) abs(fftshift(x,2));
- shift_x = @(x) fftshift(x,1);
- shift_y = @(x) fftshift(x,1);
- case 2
- process = @(x) abs((x));
- shift_x = @(x) (x);
- shift_y = @(x) (x);
- case 3
- process = @(x) abs(fftshift(x,2));
- shift_x = @(x) fftshift(x,2);
- shift_y = @(x) fftshift(x,2);
- end
-end
-
-switch OPTIONS.NAME
- case '\phi' %ES pot
- NAME = 'phi';
- FLD_ = DATA.PHI(:,:,:,FRAMES); % data to plot
- OPE_ = 1; % Operation on data
- case '\psi' %EM pot
- NAME = 'psi';
- FLD_ = DATA.PSI(:,:,:,FRAMES);
- OPE_ = 1;
- case '\phi^{NZ}' % non-zonal ES pot
- NAME = 'phiNZ';
- FLD_ = DATA.PHI(:,:,:,FRAMES);
- OPE_ = (KY~=0);
- case 'v_{Ey}' % y-comp of ExB velocity
- NAME = 'vy';
- FLD_ = DATA.PHI(:,:,:,FRAMES);
- OPE_ = -1i*KX;
- case 'v_{Ex}' % x-comp of ExB velocity
- NAME = 'vx';
- FLD_ = DATA.PHI(:,:,:,FRAMES);
- OPE_ = -1i*KY;
- case 's_{Ey}' % y-comp of ExB shear
- NAME = 'sy';
- FLD_ = DATA.PHI(:,:,:,FRAMES);
- OPE_ = KX.^2;
- case 's_{Ex}' % x-comp of ExB shear
- NAME = 'sx';
- FLD_ = DATA.PHI(:,:,:,FRAMES);
- OPE_ = KY.^2;
- case '\omega_z' % ES pot vorticity
- NAME = 'vorticity';
- FLD_ = DATA.PHI(:,:,:,FRAMES);
- OPE_ = -(KX.^2+KY.^2);
- case 'N_i^{00}' % first ion gyromoment
- NAME = 'Ni00';
- FLD_ = DATA.Ni00(:,:,:,FRAMES);
- OPE_ = 1;
- case 'N_e^{00}' % first electron gyromoment
- NAME = 'Ne00';
- FLD_ = DATA.Ne00(:,:,:,FRAMES);
- OPE_ = 1;
- case 'N_i^{00}-N_e^{00}' % first electron gyromoment
- NAME = 'Ni00-Ne00';
- FLD_ = (DATA.Ni00(:,:,:,FRAMES)-DATA.Ne00(:,:,:,FRAMES))./(poisson_op+(poisson_op==0));
- OPE_ = 1;
- case 'n_i' % ion density
- NAME = 'ni';
- FLD_ = DATA.DENS_I(:,:,:,FRAMES) - adiab_i.* DATA.PHI(:,:,:,FRAMES);
- OPE_ = 1;
- case 'n_e' % electron density
- NAME = 'ne';
- FLD_ = DATA.DENS_E(:,:,:,FRAMES) - adiab_e.* DATA.PHI(:,:,:,FRAMES);
- OPE_ = 1;
- case 'k^2n_e' % electron vorticity
- NAME = 'k2ne';
- FLD_ = DATA.DENS_E(:,:,:,FRAMES);
- OPE_ = -(KX.^2+KY.^2);
- case 'n_i-n_e' % polarisation
- NAME = 'pol';
- OPE_ = 1;
- FLD_ = ((DATA.DENS_I(:,:,:,FRAMES)- adiab_i.* DATA.PHI(:,:,:,FRAMES))...
- -(DATA.DENS_E(:,:,:,FRAMES)- adiab_e.* DATA.PHI(:,:,:,FRAMES)));
- case 'T_i' % ion temperature
- NAME = 'Ti';
- FLD_ = DATA.TEMP_I(:,:,:,FRAMES);
- OPE_ = 1;
- case 'G_x' % ion particle flux
- NAME = 'Gx';
- FLD_ = 0.*DATA.PHI(:,:,:,FRAMES);
- OPE_ = 1;
- for it = 1:numel(FRAMES)
- tmp = zeros(DATA.Ny,DATA.Nx,Nz);
- for iz = 1:DATA.Nz
- vx_ = real((ifourier_GENE(-1i*KY.*(DATA.PHI (:,:,iz,FRAMES(it))))));
- ni_ = real((ifourier_GENE( DATA.DENS_I(:,:,iz,FRAMES(it)))));
- gx_ = vx_.*ni_;
-% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.Ny,DATA.Nx))),2));
- tmp(:,:,iz) = abs((squeeze(fft2(gx_,DATA.Ny,DATA.Nx))));
- end
- FLD_(:,:,it)= squeeze(compr(tmp(1:DATA.Nky,1:DATA.Nkx,:)));
- end
- case 'Q_x' % ion heat flux
- NAME = 'Qx';
- FLD_ = 0.*DATA.PHI(:,:,:,FRAMES);
- OPE_ = 1;
- for it = 1:numel(FRAMES)
- tmp = zeros(DATA.Ny,DATA.Nx,Nz);
- for iz = 1:DATA.Nz
- vx_ = real((ifourier_GENE(-1i*KY.*(DATA.PHI (:,:,iz,FRAMES(it))))));
- ni_ = real((ifourier_GENE( DATA.DENS_I(:,:,iz,FRAMES(it)))));
- Ti_ = real((ifourier_GENE( DATA.TEMP_I(:,:,iz,FRAMES(it)))));
- qx_ = vx_.*ni_.*Ti_;
-% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.Ny,DATA.Nx))),2));
- tmp(:,:,iz) = abs((squeeze(fft2(qx_,DATA.Ny,DATA.Nx))));
- end
- FLD_(:,:,it)= squeeze(compr(tmp(1:DATA.Nky,1:DATA.Nkx,:)));
- end
- case 'f_i'
- SKIP_COMP = 1;
- shift_x = @(x) x; shift_y = @(y) y;
- NAME = 'fi'; OPTIONS.SPECIE = 'i';
- for it = 1:numel(FRAMES)
- OPTIONS.T = DATA.Ts5D(FRAMES(it));
- OPTIONS.Z = OPTIONS.COMP;
- [~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS);
- end
- case 'f_e'
- SKIP_COMP = 1;
- shift_x = @(x) x; shift_y = @(y) y;
- NAME = 'fe'; OPTIONS.SPECIE = 'e';
- [~,it0_] =min(abs(OPTIONS.TIME(1)-DATA.Ts5D));
- [~,it1_]=min(abs(OPTIONS.TIME(end)-DATA.Ts5D));
- dit_ = 1;%ceil((it1_-it0_+1)/10);
- FRAMES = it0_:dit_:it1_;
- iz = 1;
- for it = 1:numel(FRAMES)
- OPTIONS.T = DATA.Ts5D(FRAMES(it));
- [~,~,FIELD(:,:,it)] = compute_fa(DATA,OPTIONS);
- end
- otherwise
- disp('Fieldname not recognized');
- return
-end
-% Process the field according to the 2D plane and the space (real/cpx)
-if ~SKIP_COMP
-if COMPDIM == 3
- for it = 1:numel(FRAMES)
- tmp = squeeze(compr(OPE_.*FLD_(:,:,:,it)));
- FIELD(:,:,it) = squeeze(process(tmp));
- end
-else
- if REALP
- tmp = zeros(Ny,Nx,Nz);
- else
- tmp = zeros(DATA.Nky,DATA.Nkx,Nz);
- end
- for it = 1:numel(FRAMES)
- for iz = 1:numel(DATA.grids.z)
- tmp(:,:,iz) = squeeze(process(OPE_.*FLD_(:,:,iz,it)));
- end
- FIELD(:,:,it) = squeeze(compr(tmp));
- end
-end
-end
-TOPLOT.FIELD = FIELD;
-TOPLOT.FRAMES = FRAMES;
-TOPLOT.X = shift_x(X);
-TOPLOT.Y = shift_y(Y);
-TOPLOT.FIELDNAME = FIELDNAME;
-TOPLOT.XNAME = XNAME;
-TOPLOT.YNAME = YNAME;
-TOPLOT.FILENAME = [NAME,'_',OPTIONS.PLAN,'_',COMPNAME,'_',POLARNAME];
-TOPLOT.DIMENSIONS= DIMENSIONS;
-TOPLOT.ASPECT = ASPECT;
-TOPLOT.FRAMES = FRAMES;
-TOPLOT.INTERP = INTERP;
-end
-
diff --git a/matlab/setup.m b/matlab/setup.m
index aa1c410ad82fe19fbfbeee48503839deda3cb3dc..6256ca693e0a805317d6550bd85439c1b2770eec 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -1,17 +1,14 @@
%% _______________________________________________________________________
SIMDIR = ['../results/',SIMID,'/'];
% Grid parameters
-GRID.pmaxe = PMAXE; % Electron Hermite moments
-GRID.jmaxe = JMAXE; % Electron Laguerre moments
-GRID.pmaxi = PMAXI; % Ion Hermite moments
-GRID.jmaxi = JMAXI; % Ion Laguerre moments
-GRID.Nx = NX; % x grid resolution
-GRID.Lx = LX; % x length
-GRID.Nexc = NEXC; % to extend Lx when s>0
-GRID.Ny = NY; % y ''
-GRID.Ly = LY; % y ''
-GRID.Nz = NZ; % z resolution
-GRID.Npol = NPOL; % z resolution
+GRID.pmax = PMAX; % Hermite moments
+GRID.jmax = JMAX; % Laguerre moments
+GRID.Nx = NX; % x grid resolution
+GRID.Lx = LX; % x length
+GRID.Nexc = NEXC; % to extend Lx when s>0
+GRID.Ny = NY; % y ''
+GRID.Ly = LY; % y ''
+GRID.Nz = NZ; % z resolution
if SG; GRID.SG = '.true.'; else; GRID.SG = '.false.';end;
% Geometry
GEOM.geom = ['''',GEOMETRY,''''];
@@ -23,12 +20,11 @@ GEOM.delta = DELTA; % triangularity
GEOM.zeta = ZETA; % squareness
GEOM.parallel_bc = ['''',PARALLEL_BC,''''];
GEOM.shift_y = SHIFT_Y;
+GEOM.Npol = NPOL;
% Model parameters
-MODEL.CLOS = CLOS;
-MODEL.NL_CLOS = NL_CLOS;
MODEL.LINEARITY = ['''',LINEARITY,''''];
-MODEL.KIN_E = KIN_E;
-if KIN_E; MODEL.KIN_E = '.true.'; else; MODEL.KIN_E = '.false.';end;
+MODEL.Na = NA;
+if ADIAB_E; MODEL.ADIAB_E = '.true.'; else; MODEL.ADIAB_E = '.false.';end;
MODEL.beta = BETA;
MODEL.mu_x = MU_X;
MODEL.mu_y = MU_Y;
@@ -56,6 +52,11 @@ MODEL.K_Te = K_Te;
MODEL.k_gB = k_gB; % Magnetic gradient
MODEL.k_cB = k_cB; % Magnetic curvature
MODEL.lambdaD = LAMBDAD;
+% CLOSURE parameters
+CLOSURE.hierarchy_closure = ['''',HRCY_CLOS,''''];
+CLOSURE.nonlinear_closure = ['''',NLIN_CLOS,''''];
+CLOSURE.dmax = DMAX;
+CLOSURE.nmax = NMAX;
% Collision parameters
COLL.collision_model = ['''',CO,''''];
if (GKCO); COLL.GK_CO = '.true.'; else; COLL.GK_CO = '.false.';end;
@@ -93,17 +94,10 @@ end
if ~ABCO
CONAME = [CONAME,'aa'];
end
-if (CLOS == 0); CLOSNAME = 'Trunc.';
-elseif(CLOS == 1); CLOSNAME = 'Clos. 1';
-elseif(CLOS == 2); CLOSNAME = 'Clos. 2';
-end
+CLOSNAME = [HRCY_CLOS,' dmax=',num2str(DMAX)];
+NLCLOSNAME = [NLIN_CLOS,' nmax=',num2str(NMAX)];
% Hermite-Laguerre degrees naming
-if (PMAXE == PMAXI) && (JMAXE == JMAXI)
- HLdeg_ = ['_',num2str(PMAXE+1),'x',num2str(JMAXE+1)];
-else
- HLdeg_ = ['_Pe_',num2str(PMAXE+1),'_Je_',num2str(JMAXE+1),...
- '_Pi_',num2str(PMAXI+1),'_Ji_',num2str(JMAXI+1)];
-end
+HLdeg_ = ['_',num2str(PMAX+1),'x',num2str(JMAX+1)];
% temp. dens. drives
drives_ = [];
if abs(K_Ni) > 0; drives_ = [drives_,'_kN_',num2str(K_Ni)]; end;
@@ -115,7 +109,7 @@ coll_ = sprintf(coll_,NU);
lin_ = [];
if ~LINEARITY; lin_ = '_lin'; end
adiabe_ = [];
-if ~KIN_E; adiabe_ = '_adiabe'; end
+if ADIAB_E; adiabe_ = '_adiabe'; end
% resolution and boxsize
res_ = [num2str(GRID.Nx),'x',num2str(GRID.Ny)];
if (LX ~= LY)
@@ -158,11 +152,11 @@ BASIC.maxruntime = str2num(CLUSTER.TIME(1:2))*3600 ...
+ str2num(CLUSTER.TIME(4:5))*60 ...
+ str2num(CLUSTER.TIME(7:8));
% Outputs parameters
-OUTPUTS.nsave_0d = floor(1.0/SPS0D/DT);
-OUTPUTS.nsave_1d = -1;
-OUTPUTS.nsave_2d = floor(1.0/SPS2D/DT);
-OUTPUTS.nsave_3d = floor(1.0/SPS3D/DT);
-OUTPUTS.nsave_5d = floor(1.0/SPS5D/DT);
+OUTPUTS.dtsave_0d = DTSAVE0D;
+OUTPUTS.dtsave_1d = -1;
+OUTPUTS.dtsave_2d = DTSAVE2D;
+OUTPUTS.dtsave_3d = DTSAVE3D;
+OUTPUTS.dtsave_5d = DTSAVE5D;
if W_DOUBLE; OUTPUTS.write_doubleprecision = '.true.'; else; OUTPUTS.write_doubleprecision = '.false.';end;
if W_GAMMA; OUTPUTS.write_gamma = '.true.'; else; OUTPUTS.write_gamma = '.false.';end;
if W_HF; OUTPUTS.write_hf = '.true.'; else; OUTPUTS.write_hf = '.false.';end;
@@ -184,7 +178,7 @@ end
% mkdir(BASIC.MISCDIR)
% end
%% Compile and WRITE input file
-INPUT = write_fort90(OUTPUTS,GRID,GEOM,MODEL,COLL,INITIAL,TIME_INTEGRATION,BASIC);
+INPUT = write_fort90(OUTPUTS,GRID,GEOM,MODEL,CLOSURE,COLL,INITIAL,TIME_INTEGRATION,BASIC);
nproc = 1;
MAKE = 'cd ..; make; cd wk';
% system(MAKE);
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 23bbd924d1d4ce8c384719e89753396c321791e0..badb02a29c0858294f58e5a26fd4be6180b8fe0d 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -1,26 +1,24 @@
-function [INPUT] = write_fort90(OUTPUTS,GRID,GEOM,MODEL,COLL,INITIAL,TIME_INTEGRATION,BASIC)
+function [INPUT] = write_fort90(OUTPUTS,GRID,GEOM,MODEL,CLOSURE,COLL,INITIAL,TIME_INTEGRATION,BASIC)
% Write the input script "fort.90" with desired parameters
INPUT = ['fort_',sprintf('%2.2d',OUTPUTS.job2load+1),'.90'];
fid = fopen(INPUT,'wt');
fprintf(fid,'&BASIC\n');
-fprintf(fid,[' nrun = ', num2str(BASIC.nrun),'\n']);
-fprintf(fid,[' dt = ', num2str(BASIC.dt),'\n']);
-fprintf(fid,[' tmax = ', num2str(BASIC.tmax),'\n']);
+fprintf(fid,[' nrun = ', num2str(BASIC.nrun),'\n']);
+fprintf(fid,[' dt = ', num2str(BASIC.dt),'\n']);
+fprintf(fid,[' tmax = ', num2str(BASIC.tmax),'\n']);
fprintf(fid,[' maxruntime = ', num2str(BASIC.maxruntime),'\n']);
+fprintf(fid,[' job2load = ', num2str(OUTPUTS.job2load),'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&GRID\n');
-fprintf(fid,[' pmaxe = ', num2str(GRID.pmaxe),'\n']);
-fprintf(fid,[' jmaxe = ', num2str(GRID.jmaxe),'\n']);
-fprintf(fid,[' pmaxi = ', num2str(GRID.pmaxi),'\n']);
-fprintf(fid,[' jmaxi = ', num2str(GRID.jmaxi),'\n']);
+fprintf(fid,[' pmax = ', num2str(GRID.pmax),'\n']);
+fprintf(fid,[' jmax = ', num2str(GRID.jmax),'\n']);
fprintf(fid,[' Nx = ', num2str(GRID.Nx),'\n']);
fprintf(fid,[' Lx = ', num2str(GRID.Lx),'\n']);
fprintf(fid,[' Ny = ', num2str(GRID.Ny),'\n']);
fprintf(fid,[' Ly = ', num2str(GRID.Ly),'\n']);
fprintf(fid,[' Nz = ', num2str(GRID.Nz),'\n']);
-fprintf(fid,[' Npol = ', num2str(GRID.Npol),'\n']);
fprintf(fid,[' SG = ', GRID.SG,'\n']);
fprintf(fid,[' Nexc = ', num2str(GRID.Nexc),'\n']);
fprintf(fid,'/\n');
@@ -34,32 +32,29 @@ fprintf(fid,[' kappa = ', num2str(GEOM.kappa),'\n']);
fprintf(fid,[' delta = ', num2str(GEOM.delta),'\n']);
fprintf(fid,[' zeta = ', num2str(GEOM.zeta),'\n']);
fprintf(fid,[' parallel_bc = ', GEOM.parallel_bc,'\n']);
+fprintf(fid,[' shift_y = ', num2str(GEOM.shift_y),'\n']);
+fprintf(fid,[' Npol = ', num2str(GEOM.Npol),'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&OUTPUT_PAR\n');
-fprintf(fid,[' nsave_0d = ', num2str(OUTPUTS.nsave_0d),'\n']);
-fprintf(fid,[' nsave_1d = ', num2str(OUTPUTS.nsave_1d),'\n']);
-fprintf(fid,[' nsave_2d = ', num2str(OUTPUTS.nsave_2d),'\n']);
-fprintf(fid,[' nsave_3d = ', num2str(OUTPUTS.nsave_3d),'\n']);
-fprintf(fid,[' nsave_5d = ', num2str(OUTPUTS.nsave_5d),'\n']);
+fprintf(fid,[' dtsave_0d = ', num2str(OUTPUTS.dtsave_0d),'\n']);
+fprintf(fid,[' dtsave_1d = ', num2str(OUTPUTS.dtsave_1d),'\n']);
+fprintf(fid,[' dtsave_2d = ', num2str(OUTPUTS.dtsave_2d),'\n']);
+fprintf(fid,[' dtsave_3d = ', num2str(OUTPUTS.dtsave_3d),'\n']);
+fprintf(fid,[' dtsave_5d = ', num2str(OUTPUTS.dtsave_5d),'\n']);
fprintf(fid,[' write_doubleprecision = ', OUTPUTS.write_doubleprecision,'\n']);
fprintf(fid,[' write_gamma = ', OUTPUTS.write_gamma,'\n']);
fprintf(fid,[' write_hf = ', OUTPUTS.write_hf,'\n']);
fprintf(fid,[' write_phi = ', OUTPUTS.write_phi,'\n']);
fprintf(fid,[' write_Na00 = ', OUTPUTS.write_Na00,'\n']);
fprintf(fid,[' write_Napj = ', OUTPUTS.write_Napj,'\n']);
-fprintf(fid,[' write_Sapj = ', OUTPUTS.write_Sapj,'\n']);
fprintf(fid,[' write_dens = ', OUTPUTS.write_dens,'\n']);
fprintf(fid,[' write_temp = ', OUTPUTS.write_temp,'\n']);
-fprintf(fid,[' job2load = ', num2str(OUTPUTS.job2load),'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&MODEL_PAR\n');
-fprintf(fid,' ! Collisionality\n');
-fprintf(fid,[' CLOS = ', num2str(MODEL.CLOS),'\n']);
-fprintf(fid,[' NL_CLOS = ', num2str(MODEL.NL_CLOS),'\n']);
fprintf(fid,[' LINEARITY = ', MODEL.LINEARITY,'\n']);
-fprintf(fid,[' KIN_E = ', MODEL.KIN_E,'\n']);
+fprintf(fid,[' Na = ', num2str(MODEL.Na),'\n']);
fprintf(fid,[' mu_x = ', num2str(MODEL.mu_x),'\n']);
fprintf(fid,[' mu_y = ', num2str(MODEL.mu_y),'\n']);
fprintf(fid,[' N_HD = ', num2str(MODEL.N_HD),'\n']);
@@ -68,22 +63,40 @@ fprintf(fid,[' HYP_V = ', MODEL.HYP_V,'\n']);
fprintf(fid,[' mu_p = ', num2str(MODEL.mu_p),'\n']);
fprintf(fid,[' mu_j = ', num2str(MODEL.mu_j),'\n']);
fprintf(fid,[' nu = ', num2str(MODEL.nu),'\n']);
-fprintf(fid,[' tau_e = ', num2str(MODEL.tau_e),'\n']);
-fprintf(fid,[' tau_i = ', num2str(MODEL.tau_i),'\n']);
-fprintf(fid,[' sigma_e = ', num2str(MODEL.sigma_e),'\n']);
-fprintf(fid,[' sigma_i = ', num2str(MODEL.sigma_i),'\n']);
-fprintf(fid,[' q_e = ', num2str(MODEL.q_e),'\n']);
-fprintf(fid,[' q_i = ', num2str(MODEL.q_i),'\n']);
-fprintf(fid,[' K_Ne = ', num2str(MODEL.K_Ne),'\n']);
-fprintf(fid,[' K_Te = ', num2str(MODEL.K_Te),'\n']);
-fprintf(fid,[' K_Ni = ', num2str(MODEL.K_Ni),'\n']);
-fprintf(fid,[' K_Ti = ', num2str(MODEL.K_Ti),'\n']);
-fprintf(fid,[' k_gB = ', num2str(MODEL.k_gB),'\n']);
-fprintf(fid,[' k_cB = ', num2str(MODEL.k_cB),'\n']);
+fprintf(fid,[' k_gB = ', num2str(MODEL.k_gB),'\n']);
+fprintf(fid,[' k_cB = ', num2str(MODEL.k_cB),'\n']);
fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),'\n']);
fprintf(fid,[' beta = ', num2str(MODEL.beta),'\n']);
+fprintf(fid,[' ADIAB_E = ', MODEL.ADIAB_E,'\n']);
fprintf(fid,'/\n');
+fprintf(fid,'&CLOSURE_PAR\n');
+fprintf(fid,[' hierarchy_closure=',CLOSURE.hierarchy_closure,'\n']);
+fprintf(fid,[' dmax =',num2str(CLOSURE.dmax),'\n']);
+fprintf(fid,[' nonlinear_closure=',CLOSURE.nonlinear_closure,'\n']);
+fprintf(fid,[' nmax =',num2str(CLOSURE.nmax),'\n']);
+fprintf(fid,'/\n');
+
+fprintf(fid,'&SPECIES\n');
+fprintf(fid, ' name_ = ions \n');
+fprintf(fid,[' tau_ = ', num2str(MODEL.tau_i),'\n']);
+fprintf(fid,[' sigma_ = ', num2str(MODEL.sigma_i),'\n']);
+fprintf(fid,[' q_ = ', num2str(MODEL.q_i),'\n']);
+fprintf(fid,[' K_N_ = ', num2str(MODEL.K_Ni),'\n']);
+fprintf(fid,[' K_T_ = ', num2str(MODEL.K_Ti),'\n']);
+fprintf(fid,'/\n');
+
+if(MODEL.Na > 1)
+ fprintf(fid,'&SPECIES\n');
+ fprintf(fid, ' name_ = electrons');
+ fprintf(fid,[' tau_ = ', num2str(MODEL.tau_e),'\n']);
+ fprintf(fid,[' sigma_ = ', num2str(MODEL.sigma_e),'\n']);
+ fprintf(fid,[' q_ = ', num2str(MODEL.q_e),'\n']);
+ fprintf(fid,[' K_N_ = ', num2str(MODEL.K_Ne),'\n']);
+ fprintf(fid,[' K_T_ = ', num2str(MODEL.K_Te),'\n']);
+ fprintf(fid,'/\n');
+end
+
fprintf(fid,'&COLLISION_PAR\n');
fprintf(fid,[' collision_model = ', COLL.collision_model,'\n']);
fprintf(fid,[' GK_CO = ', COLL.GK_CO,'\n']);
diff --git a/matlab/write_fort90_COSOlver.m b/matlab/write_fort90_COSOlver.m
deleted file mode 100644
index 066f825df26a15accef92e6210f4222fda192bbf..0000000000000000000000000000000000000000
--- a/matlab/write_fort90_COSOlver.m
+++ /dev/null
@@ -1,84 +0,0 @@
-INPUT = [COSOlver_path,'fort.90'];
-fid = fopen(INPUT,'wt');
-
-fprintf(fid,'&BASIC\n');
-fprintf(fid,['pmaxe = ', num2str(COSOLVER.pmaxe),'\n']);
-fprintf(fid,['jmaxe = ', num2str(COSOLVER.jmaxe),'\n']);
-fprintf(fid,['pmaxi = ', num2str(COSOLVER.pmaxi),'\n']);
-fprintf(fid,['jmaxi = ', num2str(COSOLVER.jmaxi),'\n']);
-fprintf(fid,'JEmaxx=12\n');
-fprintf(fid,'PMmaxx=16\n');
-fprintf(fid,'\n');
-
-fprintf(fid,['neFLR=',num2str(COSOLVER.neFLR),'\n']);
-fprintf(fid,['neFLRs =',num2str(COSOLVER.neFLRs),'\n']);
-fprintf(fid,['npeFLR=',num2str(COSOLVER.npeFLR),'\n']);
-fprintf(fid,['niFLR=',num2str(COSOLVER.niFLR),'\n']);
-fprintf(fid,['niFLRs=',num2str(COSOLVER.niFLRs),'\n']);
-fprintf(fid,['npiFLR=',num2str(COSOLVER.npiFLR),'\n']);
-fprintf(fid,'\n');
-fprintf(fid,['eecolls=.true.','\n']);
-fprintf(fid,['iicolls=.true.','\n']);
-fprintf(fid,['eicolls=.true.','\n']);
-fprintf(fid,['iecolls=.true.','\n']);
-fprintf(fid,'/\n');
-fprintf(fid,'\n');
-
-fprintf(fid,'&BASIS_TRANSFORMATION_PAR\n');
-fprintf(fid,['T5dir = ','''/misc/coeffs_backup/T5src/''','\n']);
-fprintf(fid,['T4dir = ','''/misc/T4/NNT4_L000x200_K000x200_P000x200_J000x200/''','\n']);
-fprintf(fid,['idxT4max = ',num2str(COSOLVER.idxT4max),'\n']);
-fprintf(fid,'idxT5max = 0\n');
-fprintf(fid,'IFT4 = .true.\n');
-fprintf(fid,'IFT5 = .false.\n');
-fprintf(fid,'/\n');
-fprintf(fid,'\n');
-
-fprintf(fid,'&MODEL_PAR \n');
-fprintf(fid,'nu=1\n');
-fprintf(fid,'mu=0.023338\n');
-fprintf(fid,'tau=1\n');
-fprintf(fid,['kperp=',num2str(COSOLVER.kperp,16),'\n']);
-fprintf(fid,'/\n');
-fprintf(fid,'\n');
-
-fprintf(fid,'&OPERATOR_MODEL \n');
-fprintf(fid,['ETEST=', num2str(COSOLVER.co),'\n']);
-fprintf(fid,['EBACK=', num2str(COSOLVER.co),'\n']);
-fprintf(fid,['ITEST=', num2str(COSOLVER.co),'\n']);
-fprintf(fid,['IBACK=', num2str(COSOLVER.co),'\n']);
-fprintf(fid,'\n');
-fprintf(fid,['ESELF=', num2str(COSOLVER.co),'\n']);
-fprintf(fid,['ISELF=', num2str(COSOLVER.co),'\n']);
-fprintf(fid,'\n');
-fprintf(fid,['GKE = ', num2str(COSOLVER.gk),'\n']);
-fprintf(fid,['GKI = ', num2str(COSOLVER.gk),'\n']);
-fprintf(fid,'DKTEST = .F.\n');
-fprintf(fid,'DKBACK = .F.\n');
-fprintf(fid,'ADDTEST = .T.\n');
-fprintf(fid,'ADDBACK = .T.\n');
-fprintf(fid,'only_gk_part = .false.\n');
-fprintf(fid,'only_symmetric = .false.\n');
-fprintf(fid,'/\n');
-fprintf(fid,'\n');
-
-fprintf(fid,'&MPI\n');
-fprintf(fid,'Nprocs_j_ii = 1\n');
-fprintf(fid,'MPI_balanced = .false.\n');
-fprintf(fid,'/\n');
-fprintf(fid,'\n');
-
-fprintf(fid,'&OUTPUT_PAR\n');
-fprintf(fid,'OVERWRITE=.true.\n');
-fprintf(fid,'nsave_ei=0\n');
-fprintf(fid,'nsave_ie=0\n');
-fprintf(fid,'nsave_ee=0\n');
-fprintf(fid,'nsave_ii=2\n');
-fprintf(fid,'ifrestart=.false.\n');
-fprintf(fid,['suffix_resfile=','''''','\n']);
-fprintf(fid,['outdir = ','''../../../HeLaZ/iCa''','\n']);
-fprintf(fid,'/\n');
-fprintf(fid,'\n');
-
-
-fclose(fid);
\ No newline at end of file
diff --git a/matlab/write_sbash_daint.m b/matlab/write_sbash_daint.m
deleted file mode 100644
index 536ff2f725889c9bd7ab3dbc24b4fb3d920719f0..0000000000000000000000000000000000000000
--- a/matlab/write_sbash_daint.m
+++ /dev/null
@@ -1,57 +0,0 @@
-% Write the input script "fort.90" with desired parameters
-INPUT = 'setup_and_run.sh';
-fid = fopen(INPUT,'wt');
-
-fprintf(fid,[...
-'#!/bin/bash\n',...
-'mkdir -p $SCRATCH/HeLaZ/wk\n',...
-...
-'cd $SCRATCH/HeLaZ/wk/\n',...
-...
-'mkdir -p ', BASIC.RESDIR,'\n',...
-'cd ',BASIC.RESDIR,'\n',...
-'cp $HOME/HeLaZ/wk/fort*.90 .\n',...
-'cp $HOME/HeLaZ/wk/batch_script.sh .\n',...
-...
-'jid=$(sbatch batch_script.sh)\n',...
-'echo $jid\n',...
-'echo to check output log :\n',...
-'echo tail -f $SCRATCH/HeLaZ/results/',BASIC.SIMID,'/',BASIC.PARAMS,'/out.txt']);
-
-fclose(fid);
-system(['cp setup_and_run.sh ',BASIC.RESDIR,'/.']);
-
-% Write the sbatch script
-INPUT = 'batch_script.sh';
-fid = fopen(INPUT,'wt');
-
-fprintf(fid,[...
-'#!/bin/bash\n',...
-'#SBATCH --job-name="',CLUSTER.JNAME,'"\n',...
-'#SBATCH --time=', CLUSTER.TIME,'\n',...
-'#SBATCH --nodes=', CLUSTER.NODES,'\n',...
-'#SBATCH --cpus-per-task=', CLUSTER.CPUPT,'\n',...
-'#SBATCH --ntasks-per-node=', CLUSTER.NTPN,'\n',...
-'#SBATCH --ntasks-per-core=', CLUSTER.NTPC,'\n',...
-'#SBATCH --mem=', CLUSTER.MEM,'\n',...
-'#SBATCH --error=err.txt\n',...
-'#SBATCH --output=out.txt\n',...
-'#SBATCH --account="s882"\n',...
-'#SBATCH --constraint=mc\n',...
-'#SBATCH --hint=nomultithread\n',...
-'#SBATCH --partition=',CLUSTER.PART,'\n',...
-'export OMP_NUM_THREADS=$SLURM_CPUS_PER_TASK\n',...
-...% '#SBATCH --job-name=',PARAMS,'\n\n',...
-'module purge\n',...
-'module load PrgEnv-intel\n',...
-'module load cray-hdf5-parallel\n',...
-'module load cray-mpich\n',...
-'module load craype-x86-skylake\n',...
-'module load cray-fftw\n',...
-'srun --cpu-bind=cores ./../../../bin/helaz ',num2str(NP_P),' ',num2str(NP_KX)]);
-%'srun ./../../../bin/helaz']);
-
-fclose(fid);
-system(['cp batch_script.sh ',BASIC.RESDIR,'/.']);
-
-system('scp {fort*.90,setup_and_run.sh,batch_script.sh} ahoffman@ela.cscs.ch:HeLaZ/wk');
diff --git a/matlab/write_sbash_izar.m b/matlab/write_sbash_izar.m
deleted file mode 100644
index 02224f474568df05d54ca53d02805abffcfadb77..0000000000000000000000000000000000000000
--- a/matlab/write_sbash_izar.m
+++ /dev/null
@@ -1,46 +0,0 @@
-% Write the input script "fort.90" with desired parameters
-INPUT = 'setup_and_run.sh';
-fid = fopen(INPUT,'wt');
-
-fprintf(fid,[...
-'#!/bin/bash\n',...
-'mkdir -p $SCRATCH/ahoffman/HeLaZ/wk\n',...
-...
-'cd $SCRATCH/ahoffman/HeLaZ/wk/\n',...
-...
-'mkdir -p ', BASIC.RESDIR,'\n',...
-'cd ',BASIC.RESDIR,'\n',...
-'cp $HOME/HeLaZ/wk/fort.90 .\n',...
-'cp $HOME/HeLaZ/wk/batch_script.sh .\n',...
-...
-'sbatch batch_script.sh\n',...
-'echo tail -f $SCRATCH/ahoffman/HeLaZ/results/',BASIC.SIMID,'/',BASIC.PARAMS,'/out.txt']);
-
-fclose(fid);
-system(['cp setup_and_run.sh ',BASIC.RESDIR,'/.']);
-
-% Write the sbatch script
-INPUT = 'batch_script.sh';
-fid = fopen(INPUT,'wt');
-
-fprintf(fid,[...
-'#!/bin/bash\n',...
-'#SBATCH --chdir $SCRATCH/ahoffman/HeLaZ/results/',BASIC.SIMID,'/',BASIC.PARAMS,'\n',...
-'#SBATCH --job-name ',CLUSTER.JNAME,'\n',...
-'#SBATCH --time ', CLUSTER.TIME,'\n',...
-'#SBATCH --nodes ', CLUSTER.NODES,'\n',...
-'#SBATCH --cpus-per-task ', CLUSTER.CPUPT,'\n',...
-'#SBATCH --ntasks-per-node ', CLUSTER.NTPN,'\n',...
-'#SBATCH --mem ', CLUSTER.MEM,'\n',...
-'#SBATCH --partition ',CLUSTER.PART,'\n',...
-'module purge\n',...
-'module load intel/19.0.5\n',...
-'module load intel-mpi/2019.5.281\n',...
-'module load fftw/3.3.8-mpi-openmp\n',...
-'module load hdf5/1.10.6-mpi\n',...
-'srun --cpu-bind=cores ./../../../bin/helaz ',num2str(NP_P),' ',num2str(NP_KX)]);
-
-fclose(fid);
-system(['cp batch_script.sh ',BASIC.RESDIR,'/.']);
-
-system('scp {fort.90,setup_and_run.sh,batch_script.sh} ahoffman@izar.epfl.ch:/home/ahoffman/HeLaZ/wk');
\ No newline at end of file
diff --git a/matlab/write_sbash_marconi.m b/matlab/write_sbash_marconi.m
deleted file mode 100644
index 28a55760a6c3c80fcf0165fa208a89c21ab8bdf2..0000000000000000000000000000000000000000
--- a/matlab/write_sbash_marconi.m
+++ /dev/null
@@ -1,46 +0,0 @@
-% Write the input script "fort.90" with desired parameters
-INPUT = 'setup_and_run.sh';
-fid = fopen(INPUT,'wt');
-
-fprintf(fid,[...
-'#!/bin/bash\n',...
-'mkdir -p $CINECA_SCRATCH/HeLaZ/wk\n',...
-...
-'cd $CINECA_SCRATCH/HeLaZ/wk/\n',...
-...
-'mkdir -p ', BASIC.RESDIR,'\n',...
-'cd ',BASIC.RESDIR,'\n',...
-'cp $HOME/HeLaZ/wk/fort*.90 .\n',...
-'cp $HOME/HeLaZ/wk/batch_script.sh .\n',...
-'rm $HOME/HeLaZ/wk/fort*.90\n',...
-...
-SBATCH_CMD,...
-'echo tail -f $CINECA_SCRATCH/HeLaZ',BASIC.RESDIR(3:end)]);
-
-fclose(fid);
-system(['cp setup_and_run.sh ',BASIC.RESDIR,'/.']);
-
-% Write the sbatch script
-INPUT = 'batch_script.sh';
-fid = fopen(INPUT,'wt');
-
-fprintf(fid,[...
-'#!/bin/bash\n',...
-'#SBATCH --job-name=',CLUSTER.JNAME,'\n',...
-'#SBATCH --time=', CLUSTER.TIME,'\n',...
-'#SBATCH --nodes=', CLUSTER.NODES,'\n',...
-'#SBATCH --cpus-per-task=', CLUSTER.CPUPT,'\n',...
-'#SBATCH --ntasks-per-node=', CLUSTER.NTPN,'\n',...
-'#SBATCH --mem=', CLUSTER.MEM,'\n',...
-'#SBATCH --error=err_',sprintf('%2.2d',JOB2LOAD+1),'.txt\n',...
-'#SBATCH --output=out_',sprintf('%2.2d',JOB2LOAD+1),'.txt\n',...
-'#SBATCH --account=FUA35_TSVVT421\n',...
-'#SBATCH --partition=skl_fua_',CLUSTER.PART,'\n',...
-'module load autoload hdf5 fftw\n',...
-'srun --cpu-bind=cores ./../../../bin/',EXECNAME,' ',num2str(NP_P),' ',num2str(NP_KX),' ',num2str(JOB2LOAD+1)]);
-
-fclose(fid);
-system(['cp batch_script.sh ',BASIC.RESDIR,'/.']);
-
-system('scp {fort*.90,setup_and_run.sh,batch_script.sh} ahoffman@login.marconi.cineca.it:/marconi/home/userexternal/ahoffman/HeLaZ/wk > trash.txt');
-system('rm trash.txt');
\ No newline at end of file
diff --git a/testcases/ITG_zpinch/fort.90 b/testcases/ITG_zpinch/fort.90
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/testcases/ITG_zpinch/fort_00.90 b/testcases/ITG_zpinch/fort_00.90
new file mode 100644
index 0000000000000000000000000000000000000000..536bb949c3e8f41656a66fc285698f82dad53429
--- /dev/null
+++ b/testcases/ITG_zpinch/fort_00.90
@@ -0,0 +1,104 @@
+&BASIC
+ nrun = 99999999
+ dt = 0.01
+ tmax = 250
+ maxruntime = 72000
+ job2load = -1
+/
+&GRID
+ pmax = 2
+ jmax = 1
+ Nx = 64
+ Lx = 120
+ Ny = 48
+ Ly = 120
+ Nz = 16
+ SG = .f.
+ Nexc = 1
+/
+&GEOMETRY
+ geom = 'Z-pinch'
+ q0 = 0.0
+ shear = 0.0
+ eps = 0.0
+ kappa = 1.0
+ s_kappa= 0.0
+ delta = 0.0
+ s_delta= 0.0
+ zeta = 0.0
+ s_zeta = 0.0
+ parallel_bc = 'shearless'
+ shift_y= 0.0
+/
+&OUTPUT_PAR
+ dtsave_0d = 1
+ dtsave_1d = -1
+ dtsave_2d = -1
+ dtsave_3d = 1
+ dtsave_5d = 10
+ write_doubleprecision = .f.
+ write_gamma = .t.
+ write_hf = .t.
+ write_phi = .t.
+ write_Na00 = .t.
+ write_Napj = .t.
+ write_dens = .t.
+ write_fvel = .t.
+ write_temp = .t.
+/
+&MODEL_PAR
+ LINEARITY = 'nonlinear'
+ Na = 1 ! number of species
+ mu_x = 0.0
+ mu_y = 0.0
+ N_HD = 4
+ mu_z = 0.0
+ HYP_V = 'hypcoll'
+ mu_p = 0.0
+ mu_j = 0.0
+ nu = 0.05
+ beta = 0.0
+ ADIAB_E = .t.
+ tau_e = 1.0
+/
+&CLOSURE_PAR
+ !hierarchy_closure='truncation'
+ hierarchy_closure='max_degree'
+ dmax = 2
+ nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing)
+ nmax = -1
+/
+&SPECIES
+ ! ions
+ name_ = 'ions'
+ tau_ = 1.0
+ sigma_= 1.0
+ q_ = 1.0
+ k_N_ = 0.0
+ k_T_ = 4.0
+/
+&SPECIES
+ ! electrons
+ name_ = 'electrons'
+ tau_ = 1.0
+ sigma_= 0.023338
+ q_ =-1.0
+ k_N_ = 1.6
+ k_T_ = 0.4
+/
+&COLLISION_PAR
+ collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+ GK_CO = .t.
+ INTERSPECIES = .true.
+ mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+/
+&INITIAL_CON
+ INIT_OPT = 'blob' !(phi,blob)
+ ACT_ON_MODES = 'donothing'
+ init_background = 0.0
+ init_noiselvl = 0.005
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
diff --git a/testcases/ITG_zpinch/gyacomo23_sp b/testcases/ITG_zpinch/gyacomo23_sp
new file mode 120000
index 0000000000000000000000000000000000000000..d3982f650fe02a339058b25297c18d43a3f4c36d
--- /dev/null
+++ b/testcases/ITG_zpinch/gyacomo23_sp
@@ -0,0 +1 @@
+../../bin/gyacomo23_sp
\ No newline at end of file
diff --git a/testcases/cyclone_example/fort_00.90 b/testcases/cyclone_example/fort_00.90
index c449b1150b46cd02fe5c8531682113f35b3278fb..07b169f32d8c555582743c52014e868c4de4581b 100644
--- a/testcases/cyclone_example/fort_00.90
+++ b/testcases/cyclone_example/fort_00.90
@@ -18,6 +18,7 @@
/
&GEOMETRY
geom = 's-alpha'
+ !geom = 'miller'
q0 = 1.4
shear = 0.8
eps = 0.18
@@ -52,11 +53,11 @@
mu_x = 0.0
mu_y = 0.0
N_HD = 4
- mu_z = 2.0
+ mu_z = 1.0
HYP_V = 'hypcoll'
mu_p = 0.0
mu_j = 0.0
- nu = 0.1
+ nu = 0.05
beta = 0.0
ADIAB_E = .t.
tau_e = 1.0
@@ -79,7 +80,7 @@
&COLLISION_PAR
collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
- GK_CO = .t.
+ GK_CO = .f.
INTERSPECIES = .true.
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
@@ -91,5 +92,6 @@
iseed = 42
/
&TIME_INTEGRATION_PAR
- numerical_scheme = 'RK4'
+ !numerical_scheme = 'RK4'
+ numerical_scheme = 'SSP_RK3'
/
diff --git a/testcases/cyclone_example/fort_01.90 b/testcases/cyclone_example/fort_01.90
new file mode 100644
index 0000000000000000000000000000000000000000..07b169f32d8c555582743c52014e868c4de4581b
--- /dev/null
+++ b/testcases/cyclone_example/fort_01.90
@@ -0,0 +1,97 @@
+&BASIC
+ nrun = 99999999
+ dt = 0.01
+ tmax = 50
+ maxruntime = 72000
+ job2load = -1
+/
+&GRID
+ pmax = 2
+ jmax = 1
+ Nx = 64
+ Lx = 120
+ Ny = 48
+ Ly = 120
+ Nz = 16
+ SG = .f.
+ Nexc = 0
+/
+&GEOMETRY
+ geom = 's-alpha'
+ !geom = 'miller'
+ q0 = 1.4
+ shear = 0.8
+ eps = 0.18
+ kappa = 1.0
+ s_kappa= 0.0
+ delta = 0.0
+ s_delta= 0.0
+ zeta = 0.0
+ s_zeta = 0.0
+ parallel_bc = 'dirichlet'
+ shift_y= 0.0
+/
+&OUTPUT_PAR
+ dtsave_0d = 1
+ dtsave_1d = -1
+ dtsave_2d = -1
+ dtsave_3d = 1
+ dtsave_5d = 20
+ write_doubleprecision = .f.
+ write_gamma = .t.
+ write_hf = .t.
+ write_phi = .t.
+ write_Na00 = .t.
+ write_Napj = .t.
+ write_dens = .t.
+ write_fvel = .t.
+ write_temp = .t.
+/
+&MODEL_PAR
+ LINEARITY = 'nonlinear'
+ Na = 1 ! number of species
+ mu_x = 0.0
+ mu_y = 0.0
+ N_HD = 4
+ mu_z = 1.0
+ HYP_V = 'hypcoll'
+ mu_p = 0.0
+ mu_j = 0.0
+ nu = 0.05
+ beta = 0.0
+ ADIAB_E = .t.
+ tau_e = 1.0
+/
+&CLOSURE_PAR
+ hierarchy_closure='truncation'
+ dmax = -1
+ nonlinear_closure='truncation'
+ nmax = 0
+/
+&SPECIES
+ ! ions
+ name_ = 'ions'
+ tau_ = 1.0
+ sigma_= 1.0
+ q_ = 1.0
+ k_N_ = 2.22
+ k_T_ = 6.96
+/
+
+&COLLISION_PAR
+ collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+ GK_CO = .f.
+ INTERSPECIES = .true.
+ mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+/
+&INITIAL_CON
+ INIT_OPT = 'blob'
+ ACT_ON_MODES = 'donothing'
+ init_background = 0.0
+ init_noiselvl = 0.005
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ !numerical_scheme = 'RK4'
+ numerical_scheme = 'SSP_RK3'
+/
diff --git a/testcases/zpinch_example/fort_00.90 b/testcases/zpinch_example/fort_00.90
index 2f9ce8c395d79989e6f7f3e9dd818b53657d8119..d7e89e580e73382299e1eabf2e4db3b79807a940 100644
--- a/testcases/zpinch_example/fort_00.90
+++ b/testcases/zpinch_example/fort_00.90
@@ -1,7 +1,7 @@
&BASIC
nrun = 99999999
dt = 0.01
- tmax = 25
+ tmax = 250
maxruntime = 72000
job2load = -1
/
@@ -11,7 +11,7 @@
Nx = 64
Lx = 200
Ny = 48
- Ly = 60
+ Ly = 120
Nz = 1
SG = .f.
Nexc = 1
@@ -62,11 +62,11 @@
tau_e = 1.0
/
&CLOSURE_PAR
- hierarchy_closure='truncation'
- !hierarchy_closure='max_degree'
+ !hierarchy_closure='truncation'
+ hierarchy_closure='max_degree'
dmax = 2
nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing)
- nmax = 0
+ nmax = -1
/
&SPECIES
! ions
@@ -74,7 +74,7 @@
tau_ = 1.0
sigma_= 1.0
q_ = 1.0
- k_N_ = 2.0
+ k_N_ = 1.6
k_T_ = 0.4
/
&SPECIES
@@ -83,7 +83,7 @@
tau_ = 1.0
sigma_= 0.023338
q_ =-1.0
- k_N_ = 2.0
+ k_N_ = 1.6
k_T_ = 0.4
/
&COLLISION_PAR
diff --git a/testcases/zpinch_example/fort_01.90 b/testcases/zpinch_example/fort_01.90
new file mode 100644
index 0000000000000000000000000000000000000000..657cb60135ee94e5fb8cac4481b775bfb9e64da8
--- /dev/null
+++ b/testcases/zpinch_example/fort_01.90
@@ -0,0 +1,104 @@
+&BASIC
+ nrun = 99999999
+ dt = 0.01
+ tmax = 500
+ maxruntime = 72000
+ job2load = 0
+/
+&GRID
+ pmax = 2
+ jmax = 1
+ Nx = 64
+ Lx = 200
+ Ny = 48
+ Ly = 120
+ Nz = 1
+ SG = .f.
+ Nexc = 1
+/
+&GEOMETRY
+ geom = 'Z-pinch'
+ q0 = 0.0
+ shear = 0.0
+ eps = 0.0
+ kappa = 1.0
+ s_kappa= 0.0
+ delta = 0.0
+ s_delta= 0.0
+ zeta = 0.0
+ s_zeta = 0.0
+ parallel_bc = 'shearless'
+ shift_y= 0.0
+/
+&OUTPUT_PAR
+ dtsave_0d = 1
+ dtsave_1d = -1
+ dtsave_2d = -1
+ dtsave_3d = 1
+ dtsave_5d = 10
+ write_doubleprecision = .f.
+ write_gamma = .t.
+ write_hf = .t.
+ write_phi = .t.
+ write_Na00 = .t.
+ write_Napj = .t.
+ write_dens = .t.
+ write_fvel = .t.
+ write_temp = .t.
+/
+&MODEL_PAR
+ LINEARITY = 'nonlinear'
+ Na = 2 ! number of species
+ mu_x = 1.0
+ mu_y = 1.0
+ N_HD = 4
+ mu_z = 0.0
+ HYP_V = 'hypcoll'
+ mu_p = 0.0
+ mu_j = 0.0
+ nu = 0.1
+ beta = 0.0
+ ADIAB_E = .f.
+ tau_e = 1.0
+/
+&CLOSURE_PAR
+ !hierarchy_closure='truncation'
+ hierarchy_closure='max_degree'
+ dmax = 2
+ nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing)
+ nmax = -1
+/
+&SPECIES
+ ! ions
+ name_ = 'ions'
+ tau_ = 1.0
+ sigma_= 1.0
+ q_ = 1.0
+ k_N_ = 1.6
+ k_T_ = 0.4
+/
+&SPECIES
+ ! electrons
+ name_ = 'electrons'
+ tau_ = 1.0
+ sigma_= 0.023338
+ q_ =-1.0
+ k_N_ = 1.6
+ k_T_ = 0.4
+/
+&COLLISION_PAR
+ collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+ GK_CO = .t.
+ INTERSPECIES = .true.
+ mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+/
+&INITIAL_CON
+ INIT_OPT = 'blob' !(phi,blob)
+ ACT_ON_MODES = 'donothing'
+ init_background = 0.0
+ init_noiselvl = 0.005
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
diff --git a/testcases/zpinch_example/fort_02.90 b/testcases/zpinch_example/fort_02.90
new file mode 100644
index 0000000000000000000000000000000000000000..649d9c5763d93f46345daa55adc3b902659925b4
--- /dev/null
+++ b/testcases/zpinch_example/fort_02.90
@@ -0,0 +1,104 @@
+&BASIC
+ nrun = 99999999
+ dt = 0.01
+ tmax = 700
+ maxruntime = 72000
+ job2load = 1
+/
+&GRID
+ pmax = 2
+ jmax = 1
+ Nx = 64
+ Lx = 200
+ Ny = 48
+ Ly = 120
+ Nz = 1
+ SG = .f.
+ Nexc = 1
+/
+&GEOMETRY
+ geom = 'Z-pinch'
+ q0 = 0.0
+ shear = 0.0
+ eps = 0.0
+ kappa = 1.0
+ s_kappa= 0.0
+ delta = 0.0
+ s_delta= 0.0
+ zeta = 0.0
+ s_zeta = 0.0
+ parallel_bc = 'shearless'
+ shift_y= 0.0
+/
+&OUTPUT_PAR
+ dtsave_0d = 1
+ dtsave_1d = -1
+ dtsave_2d = -1
+ dtsave_3d = 1
+ dtsave_5d = 10
+ write_doubleprecision = .f.
+ write_gamma = .t.
+ write_hf = .t.
+ write_phi = .t.
+ write_Na00 = .t.
+ write_Napj = .t.
+ write_dens = .t.
+ write_fvel = .t.
+ write_temp = .t.
+/
+&MODEL_PAR
+ LINEARITY = 'nonlinear'
+ Na = 2 ! number of species
+ mu_x = 1.0
+ mu_y = 1.0
+ N_HD = 4
+ mu_z = 0.0
+ HYP_V = 'hypcoll'
+ mu_p = 0.0
+ mu_j = 0.0
+ nu = 0.1
+ beta = 0.0
+ ADIAB_E = .f.
+ tau_e = 1.0
+/
+&CLOSURE_PAR
+ hierarchy_closure='truncation'
+ !hierarchy_closure='max_degree'
+ dmax = 2
+ nonlinear_closure='anti_laguerre_aliasing' !(truncation,full_sum,anti_laguerre_aliasing)
+ nmax = -1
+/
+&SPECIES
+ ! ions
+ name_ = 'ions'
+ tau_ = 1.0
+ sigma_= 1.0
+ q_ = 1.0
+ k_N_ = 1.6
+ k_T_ = 0.4
+/
+&SPECIES
+ ! electrons
+ name_ = 'electrons'
+ tau_ = 1.0
+ sigma_= 0.023338
+ q_ =-1.0
+ k_N_ = 1.6
+ k_T_ = 0.4
+/
+&COLLISION_PAR
+ collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+ GK_CO = .t.
+ INTERSPECIES = .true.
+ mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+/
+&INITIAL_CON
+ INIT_OPT = 'blob' !(phi,blob)
+ ACT_ON_MODES = 'donothing'
+ init_background = 0.0
+ init_noiselvl = 0.005
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
diff --git a/wk/fast_analysis.m b/wk/fast_analysis.m
index af8a6e41f1a863aead95f9067732c1288ca0119e..4789ce4411052c19afcd4d70f4ccd7ffccfd8d60 100644
--- a/wk/fast_analysis.m
+++ b/wk/fast_analysis.m
@@ -18,7 +18,7 @@ PARTITION = '/home/ahoffman/gyacomo/';
% resdir = 'paper_2_GYAC23/precision_study/test_3x2x128x64x24_sp_muz_2.0';
% resdir = 'paper_2_GYAC23/precision_study/3x2x128x64x24_sp_clos_1';
-%%
+%% Marconi results
% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/5x3x128x64x24_dp_muz_2.0_muxy_0';
% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/5x3x128x64x24_dp_SG';
% resdir = 'paper_2_GYAC23/collisionless/kT_5.3/5x3x128x64x24_dp_muz_2.0_full_NL';
@@ -31,15 +31,40 @@ PARTITION = '/home/ahoffman/gyacomo/';
% resdir = 'paper_2_GYAC23/collisionless/CBC/9x5x128x64x24_dp';
% resdir = 'paper_2_GYAC23/collisionless/CBC/11x6x128x64x24_dp';
% resdir = 'paper_2_GYAC23/collisionless/CBC/9x5x192x96x32_dp';
+
+%% low precision 3D ITG
+% resdir = 'results/paper_2_GYAC23/3x2x64x48x16/CBC_3x2x64x48x16_CLOS_1';
+% resdir = 'results/paper_2_GYAC23/3x2x64x48x16/kT_0.0';
+% resdir = 'results/paper_2_GYAC23/3x2x64x48x16/kT_3.0';
+% resdir = 'results/paper_2_GYAC23/3x2x64x48x16/kT_3.5';
+resdir = 'results/paper_2_GYAC23/3x2x64x48x16/kT_4.0';
+% resdir = 'results/paper_2_GYAC23/3x2x64x48x16/kT_4.5';
+% resdir = 'results/paper_2_GYAC23/3x2x64x48x16/kT_5.3';
+% resdir = 'results/paper_2_GYAC23/3x2x64x48x16/CBC';
+
+% resdir = 'results/paper_2_GYAC23/5x2x64x48x16/kT_3.5';
+% resdir = 'results/paper_2_GYAC23/5x2x64x48x16/kT_4.0';
+% resdir = 'results/paper_2_GYAC23/5x2x64x48x16/kT_4.5';
+% resdir = 'results/paper_2_GYAC23/5x2x64x48x16/kT_5.3';
+% resdir = 'results/paper_2_GYAC23/5x2x64x48x16/CBC';
+
+% resdir = 'results/paper_2_GYAC23/9x2x64x48x16/kT_3.5';
+% resdir = 'results/paper_2_GYAC23/9x2x64x48x16/kT_4.0';
+% resdir = 'results/paper_2_GYAC23/9x2x64x48x16/kT_4.5';
+% resdir = 'results/paper_2_GYAC23/9x2x64x48x16/kT_5.3';
+% resdir = 'results/paper_2_GYAC23/9x2x64x48x16/CBC';
+
%% testcases
-% resdir = 'testcases/zpinch_example';
-% resdir = 'testcases/cyclone_example';
-resdir = 'testcases/DLRA_zpinch/base_case';
+% resdir = 'testcases/ITG_zpinch';
+% resdir = 'testcases/zpinch_example/results_trunc';
+% resdir = 'testcases/zpinch_example/results_maxd=2';
+% resdir = 'testcases/DLRA_zpinch/base_case';
% resdir = 'testcases/DLRA_zpinch/nsv_filter_2';
% resdir = 'testcases/DLRA_zpinch/nsv_filter_6';
+% resdir = 'testcases/cyclone_example';
%%
-J0 = 00; J1 = 10;
+J0 = 00; J1 = 01;
% Load basic info (grids and time traces)
DATADIR = [PARTITION,resdir,'/'];
@@ -90,9 +115,9 @@ end
if 0
%% fields snapshots
% Options
-[data.Na00, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'Na00');
+[data.Na00, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Na00');
+data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
-data.Ni00 = reshape(squeeze(data.Na00(1,:,:,:,:)),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
options.INTERP = 0;
options.POLARPLOT = 0;
options.AXISEQUAL = 0;
@@ -115,12 +140,11 @@ if 0
%% Hermite-Laguerre spectrum
[data.Napjz, data.Ts3D] = compile_results_3Da(DATADIR,J0,J1,'Napjz');
% [data.Napjz, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'Nipjz');
-options.ST = 0;
+options.ST = 1;
options.NORMALIZED = 0;
options.LOGSCALE = 1;
+options.FILTER = 1; %filter the 50% time-average of the spectrum from
fig = show_moments_spectrum(data,options);
-% fig = show_napjz(data,options);
-% save_figure(data,fig,'.png');
end
if 0
@@ -129,8 +153,8 @@ if 0
options.NORMALIZED = 0;
options.TIME = [000:9000];
-options.KX_TW = [1 20]; %kx Growth rate time window
-options.KY_TW = [0 20]; %ky Growth rate time window
+options.KX_TW = [1 80]; %kx Growth rate time window
+options.KY_TW = [0 80]; %ky Growth rate time window
options.NMA = 1;
options.NMODES = 800;
options.iz = 'avg'; % avg or index
diff --git a/wk/lin_Ivanov.m b/wk/lin_Ivanov.m
new file mode 100644
index 0000000000000000000000000000000000000000..a6050c3fbd7ee37caeec43ce754150d219d75e37
--- /dev/null
+++ b/wk/lin_Ivanov.m
@@ -0,0 +1,173 @@
+%% QUICK RUN SCRIPT
+% This script creates a directory in /results and runs a simulation directly
+% from the Matlab framework. It is meant to run only small problems in linear
+% for benchmarking and debugging purposes since it makes Matlab "busy".
+
+%% Set up the paths for the necessary Matlab modules
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
+addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
+addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
+addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
+
+%% Set simulation parameters
+SIMID = 'dbg'; % Name of the simulation
+RUN = 1; % To run or just to load
+default_plots_options
+% EXECNAME = 'gyacomo23_dp'; % double precision
+% To compile single precision gyacomo do 'make clean; make sp' in the /gyacomo folder
+EXECNAME = 'gyacomo23_sp'; % single precision
+%% Set up physical parameters
+CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
+TAU = 1e-3; % e/i temperature ratio
+NU = 0.01/TAU; % Collision frequency
+K_Ne = 0*2.22; % ele Density
+K_Te = 0*6.96; % ele Temperature
+K_Ni = 0*2.22; % ion Density gradient drive
+K_Ti = 1.0*2/TAU; % ion Temperature
+SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+NA = 1; % number of kinetic species
+ADIAB_E = (NA==1); % adiabatic electron model
+BETA = 0.0; % electron plasma beta
+%% Set up grid parameters
+P = 2;
+J = 1;%P/2;
+PMAX = P; % Hermite basis size
+JMAX = J; % Laguerre basis size
+NX = 2; % real space x-gridpoints
+NY = 40; % real space y-gridpoints
+LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
+LY = 2*pi/0.15; % Size of the squared frequency domain in y direction
+NZ = 1; % number of perpendicular planes (parallel grid)
+SG = 0; % Staggered z grids option
+NEXC = 0; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+
+%% GEOMETRY
+GEOMETRY= 'Z-pinch';
+% GEOMETRY= 'miller';
+EPS = 0.0; % inverse aspect ratio
+Q0 = 0.0; % safety factor
+SHEAR = 0.0; % magnetic shear
+KAPPA = 1.0; % elongation
+DELTA = 0.0; % triangularity
+ZETA = 0.0; % squareness
+PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
+SHIFT_Y = 0.0; % Shift in the periodic BC in z
+NPOL = 1; % Number of poloidal turns
+
+%% TIME PARAMETERS
+TMAX = 50; % Maximal time unit
+DT = 1e-2; % Time step
+DTSAVE0D = 1; % Sampling per time unit for 0D arrays
+DTSAVE2D = -1; % Sampling per time unit for 2D arrays
+DTSAVE3D = 1; % Sampling per time unit for 3D arrays
+DTSAVE5D = 20; % Sampling per time unit for 5D arrays
+JOB2LOAD = -1; % Start a new simulation serie
+
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+GKCO = 0; % Gyrokinetic operator
+ABCO = 1; % INTERSPECIES collisions
+INIT_ZF = 0; % Initialize zero-field quantities
+% HRCY_CLOS = 'max_degree'; % Closure model for higher order moments
+HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+DMAX = -1;
+NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+NMAX = 0;
+KERN = 0; % Kernel model (0 : GK)
+INIT_OPT = 'phi'; % Start simulation with a noisy mom00/phi/allmom
+NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5)
+
+%% OUTPUTS
+W_DOUBLE = 1; % Output flag for double moments
+W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy)
+W_HF = 1; % Output flag for high-frequency potential energy
+W_PHI = 1; % Output flag for potential
+W_NA00 = 1; % Output flag for nalpha00 (density of species alpha)
+W_DENS = 1; % Output flag for total density
+W_TEMP = 1; % Output flag for temperature
+W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha)
+W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% UNUSED PARAMETERS
+% These parameters are usually not to play with in linear runs
+MU = 0.0; % Hyperdiffusivity coefficient
+MU_X = MU; % Hyperdiffusivity coefficient in x direction
+MU_Y = MU; % Hyperdiffusivity coefficient in y direction
+N_HD = 4; % Degree of spatial-hyperdiffusivity
+MU_Z = 2.0; % Hyperdiffusivity coefficient in z direction
+HYP_V = 'dvpar4'; % Kinetic-hyperdiffusivity model
+MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
+MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
+LAMBDAD = 0.0; % Lambda Debye
+NOISE0 = 1.0e-5; % Initial noise amplitude
+BCKGD0 = 0.0; % Initial background
+k_gB = 1.0; % Magnetic gradient strength
+k_cB = 1.0; % Magnetic curvature strength
+COLL_KCUT = 1000; % Cutoff for collision operator
+
+%%-------------------------------------------------------------------------
+%% RUN
+setup
+% system(['rm fort*.90']);
+% Run linear simulation
+if RUN
+ MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
+% RUN =['time mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
+ RUN =['time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 4 1 0;'];
+% RUN =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
+% RUN =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
+ MVOUT='cd ../../../wk;';
+ system([MVIN,RUN,MVOUT]);
+end
+
+%% Analysis
+% load
+filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
+LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
+FIGDIR = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
+% Load outputs from jobnummin up to jobnummax
+J0 = 0; J1 = 0;
+data = {}; % Initialize data as an empty cell array
+% load grids, inputs, and time traces
+data = compile_results_low_mem(data,LOCALDIR,J0,J1);
+
+if 0
+%% Plot heat flux evolution
+figure
+semilogy(data.Ts0D,data.HFLUX_X);
+xlabel('$tc_s/R$'); ylabel('$Q_x$');
+end
+if 1 % Activate or not
+%% plot mode evolution and growth rates
+% Load phi
+[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
+options.NORMALIZED = 0;
+options.TIME = data.Ts3D;
+ % Time window to measure the growth of kx/ky modes
+options.KX_TW = [0.5 1]*data.Ts3D(end);
+options.KY_TW = [0.5 1]*data.Ts3D(end);
+options.NMA = 1; % Set NMA option to 1
+options.NMODES = 999; % Set how much modes we study
+options.iz = 'avg'; % Compressing z
+options.ik = 1; %
+options.fftz.flag = 0; % Set fftz.flag option to 0
+fig = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig
+end
+
+if 0
+%% Hermite-Laguerre spectrum
+[data.Napjz, data.Ts3D] = compile_results_3Da(LOCALDIR,J0,J1,'Napjz');
+options.ST = 1;
+options.NORMALIZED = 0;
+options.LOGSCALE = 1;
+options.FILTER = 0; %filter the 50% time-average of the spectrum from
+fig = show_moments_spectrum(data,options);
+end
+
+
+
diff --git a/wk/local_run.m b/wk/local_run.m
new file mode 100644
index 0000000000000000000000000000000000000000..765e0d30ac7ec14ae476d588011a915f227f5e82
--- /dev/null
+++ b/wk/local_run.m
@@ -0,0 +1,162 @@
+%% QUICK RUN SCRIPT
+% This script creates a directory in /results and runs a simulation directly
+% from the Matlab framework. It is meant to run only small problems in linear
+% for benchmarking and debugging purposes since it makes Matlab "busy".
+
+%% Set up the paths for the necessary Matlab modules
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
+addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
+addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
+addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
+
+%% Set simulation parameters
+SIMID = 'dbg'; % Name of the simulation
+RUN = 1; % To run or just to load
+default_plots_options
+% EXECNAME = 'gyacomo23_sp'; % single precision
+EXECNAME = 'gyacomo23_dp'; % double precision
+
+%% Set up physical parameters
+CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
+NU = 0.05; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+K_Ne = 0*2.22; % ele Density
+K_Te = 0*6.96; % ele Temperature
+K_Ni = 1*2.22; % ion Density gradient drive
+K_Ti = 3.6; % ion Temperature
+SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+NA = 1; % number of kinetic species
+ADIAB_E = (NA==1); % adiabatic electron model
+BETA = 0.0; % electron plasma beta
+%% Set up grid parameters
+P = 10;
+J = 1;%P/2;
+PMAX = P; % Hermite basis size
+JMAX = J; % Laguerre basis size
+NX = 32; % real space x-gridpoints
+NY = 10; % real space y-gridpoints
+LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
+LY = 2*pi/0.1; % Size of the squared frequency domain in y direction
+NZ = 16; % number of perpendicular planes (parallel grid)
+SG = 0; % Staggered z grids option
+NEXC = 0; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+
+%% GEOMETRY
+GEOMETRY= 's-alpha';
+% GEOMETRY= 'miller';
+EPS = 0.18; % inverse aspect ratio
+Q0 = 1.4; % safety factor
+SHEAR = 0.8; % magnetic shear
+KAPPA = 1.0; % elongation
+DELTA = 0.0; % triangularity
+ZETA = 0.0; % squareness
+PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
+SHIFT_Y = 0.0; % Shift in the periodic BC in z
+NPOL = 1; % Number of poloidal turns
+
+%% TIME PARAMETERS
+TMAX = 200; % Maximal time unit
+DT = 2e-2; % Time step
+DTSAVE0D = 1; % Sampling per time unit for 0D arrays
+DTSAVE2D = -1; % Sampling per time unit for 2D arrays
+DTSAVE3D = 1; % Sampling per time unit for 3D arrays
+DTSAVE5D = 20; % Sampling per time unit for 5D arrays
+JOB2LOAD = -1; % Start a new simulation serie
+
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+GKCO = 0; % Gyrokinetic operator
+ABCO = 1; % INTERSPECIES collisions
+INIT_ZF = 0; % Initialize zero-field quantities
+HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+DMAX = -1;
+NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+NMAX = 0;
+KERN = 0; % Kernel model (0 : GK)
+INIT_OPT = 'phi'; % Start simulation with a noisy mom00/phi/allmom
+NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5)
+
+%% OUTPUTS
+W_DOUBLE = 1; % Output flag for double moments
+W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy)
+W_HF = 1; % Output flag for high-frequency potential energy
+W_PHI = 1; % Output flag for potential
+W_NA00 = 1; % Output flag for nalpha00 (density of species alpha)
+W_DENS = 1; % Output flag for total density
+W_TEMP = 1; % Output flag for temperature
+W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha)
+W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% UNUSED PARAMETERS
+% These parameters are usually not to play with in linear runs
+MU = 0.0; % Hyperdiffusivity coefficient
+MU_X = MU; % Hyperdiffusivity coefficient in x direction
+MU_Y = MU; % Hyperdiffusivity coefficient in y direction
+N_HD = 4; % Degree of spatial-hyperdiffusivity
+MU_Z = 2.0; % Hyperdiffusivity coefficient in z direction
+HYP_V = 'dvpar4'; % Kinetic-hyperdiffusivity model
+MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
+MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
+LAMBDAD = 0.0; % Lambda Debye
+NOISE0 = 1.0e-5; % Initial noise amplitude
+BCKGD0 = 0.0; % Initial background
+k_gB = 1.0; % Magnetic gradient strength
+k_cB = 1.0; % Magnetic curvature strength
+COLL_KCUT = 1000; % Cutoff for collision operator
+
+%%-------------------------------------------------------------------------
+%% RUN
+setup
+% system(['rm fort*.90']);
+% Run linear simulation
+if RUN
+ MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
+% RUN =['time mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
+ RUN =['time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 4 1 0;'];
+% RUN =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
+% RUN =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
+ MVOUT='cd ../../../wk;';
+ system([MVIN,RUN,MVOUT]);
+end
+
+%% Analysis
+% load
+filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
+LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
+FIGDIR = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
+% Load outputs from jobnummin up to jobnummax
+J0 = 0; J1 = 0;
+data = {}; % Initialize data as an empty cell array
+% load grids, inputs, and time traces
+data = compile_results_low_mem(data,LOCALDIR,J0,J1);
+
+if 1
+%% Plot heat flux evolution
+figure
+semilogy(data.Ts0D,data.HFLUX_X);
+xlabel('$tc_s/R$'); ylabel('$Q_x$');
+end
+if 0 % Activate or not
+%% plot mode evolution and growth rates
+% Load phi
+[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
+options.NORMALIZED = 0;
+options.TIME = data.Ts3D;
+ % Time window to measure the growth of kx/ky modes
+options.KX_TW = [0.5 1]*data.Ts3D(end);
+options.KY_TW = [0.5 1]*data.Ts3D(end);
+options.NMA = 1; % Set NMA option to 1
+options.NMODES = 999; % Set how much modes we study
+options.iz = 'avg'; % Compressing z
+options.ik = 1; %
+options.fftz.flag = 0; % Set fftz.flag option to 0
+fig = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig
+end
+
+
+
diff --git a/wk/multiple_sim_analysis.m b/wk/multiple_sim_analysis.m
new file mode 100644
index 0000000000000000000000000000000000000000..f2f82b19917ed4376bf112674e13f1f2eba1bb0d
--- /dev/null
+++ b/wk/multiple_sim_analysis.m
@@ -0,0 +1,102 @@
+clrs_ = lines(10);
+kN=2.22;
+figure
+for i = 1:3
+if i==1
+prefix = '/home/ahoffman/gyacomo/results/paper_2_GYAC23/3x2x64x48x16/';
+resdirs = {...
+ 'CBC';...
+ 'kT_5.3';...
+ 'kT_4.5';...
+ 'kT_4.0';...
+ 'kT_3.5';...
+ 'kT_3.0';...
+ };
+pname = '3x2x64x48x16'; clr_ = clrs_(1,:);
+ITG_threshold = 3.2;
+end
+
+if i==2
+prefix = '/home/ahoffman/gyacomo/results/paper_2_GYAC23/5x2x64x48x16/';
+resdirs = {...
+ 'CBC';...
+ 'kT_5.3';...
+ 'kT_4.5';...
+ 'kT_4.0';...
+ 'kT_3.5';...
+ 'kT_3.0';...
+ };
+pname = '5x2x64x48x16'; clr_ = clrs_(2,:);
+ITG_threshold = 3.7;
+end
+
+if i==3
+prefix = '/home/ahoffman/gyacomo/results/paper_2_GYAC23/9x2x64x48x16/';
+resdirs = {...
+ 'CBC';...
+ 'kT_5.3';...
+ 'kT_4.5';...
+ 'kT_4.0';...
+ 'kT_3.5';...
+ 'kT_3.0';...
+ };
+pname = '9x2x64x48x16'; clr_ = clrs_(3,:);
+ITG_threshold = 4.2;
+end
+
+if i==4
+
+ITG_threshold > 3.6;
+end
+
+x = [...
+ 6.96,...
+ 5.3,...
+ 4.5,...
+ 4.0,...
+ 3.5,...
+ 3.0...
+ ];
+
+J0 = 00; J1 = 10;
+
+Nseg = 5;
+
+Qx_avg = 0*(1:numel(resdirs));
+Qx_std = 0*(1:numel(resdirs));
+
+for i = 1:numel(resdirs)
+ data = compile_results_low_mem(data,[prefix,resdirs{i},'/'],J0,J1);
+ Trange = data.Ts0D(end)*[0.3 1.0];
+ %
+ [~,it0] = min(abs(Trange(1) -data.Ts0D));
+ [~,it1] = min(abs(Trange(end)-data.Ts0D));
+ %
+ if 0
+ Qx = data.HFLUX_X(it0:it1);
+ Qxa_ = 0*(1:Nseg);
+ for n = 1:Nseg
+ ntseg = floor((it1-it0)/n);
+ for m = 1:n
+ Qxa_(n) = Qxa_(n) + mean(Qx((1:ntseg)+(m-1)*ntseg));
+ end
+ Qxa_(n) = Qxa_(n)/n;
+ end
+ Qx_avg(i) = mean(Qxa_);
+ Qx_std(i) = std(Qxa_);
+ else
+ Qx_avg(i) = mean(data.HFLUX_X(it0:it1));
+ Qx_std(i) = std(data.HFLUX_X(it0:it1));
+ end
+end
+Chi_avg = Qx_avg./x/kN;
+Chi_std = Qx_std./x/kN;
+% plot;
+errorbar(x,Chi_avg,Chi_std,'DisplayName',pname,'color',clr_); hold on;
+plot(ITG_threshold*[1 1],[0 20],'-.','DisplayName','$\kappa_T^{crit}$',...
+ 'color',clr_);
+end
+ylabel('$\chi$');
+xlabel('$\kappa_T (\kappa_N=2.22)$');
+ylim([0,10]);
+legend('show');
diff --git a/wk/heat_flux_postproc.m b/wk/old scripts/heat_flux_postproc.m
similarity index 100%
rename from wk/heat_flux_postproc.m
rename to wk/old scripts/heat_flux_postproc.m
diff --git a/wk/old scripts/local_run.m b/wk/old scripts/local_run.m
deleted file mode 100644
index 21aca7320d10d1013afc7ef8aa1d1c802e241bbd..0000000000000000000000000000000000000000
--- a/wk/old scripts/local_run.m
+++ /dev/null
@@ -1,81 +0,0 @@
-addpath(genpath('../matlab')) % ... add
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Set Up parameters
-CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% PHYSICAL PARAMETERS
-NU = 0.005; % Collision frequency
-K_N = 1.4; % Density gradient drive
-K_T = 0.0; % Temperature '''
-K_E = 0.0; % Electrostat gradient
-SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
-NU_HYP = 0.0;
-KIN_E = 1; % Kinetic (1) or adiabatic (2) electron model
-%% GRID PARAMETERS
-NX = 150; % Spatial radial resolution ( = 2x radial modes)
-LX = 100; % Radial window size
-NY = 150; % Spatial azimuthal resolution (= azim modes)
-LY = 100; % Azimuthal window size
-NZ = 1; % number of perpendicular planes (parallel grid)
-P = 4;
-J = 2;
-%% GEOMETRY PARAMETERS
-Q0 = 1.0; % safety factor
-SHEAR = 0.0; % magnetic shear
-EPS = 0.0; % inverse aspect ratio
-k_gB = 1.0; % Magnetic gradient
-k_cB = 1.0; % Magnetic curvature
-SG = 0; % Staggered z grids option
-%% TIME PARAMETERS
-TMAX = 120; % Maximal time unit
-DT = 2e-2; % Time step
-SPS0D = 1; % Sampling per time unit for profiler
-SPS2D = 1; % Sampling per time unit for 2D arrays
-SPS3D = 1; % Sampling per time unit for 3D arrays
-SPS5D = 1/200; % Sampling per time unit for 5D arrays
-SPSCP = 0; % Sampling per time unit for checkpoints/10
-JOB2LOAD= -1;
-%% OPTIONS AND NAMING
-% Collision operator
-% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-CO = 'SG';
-GKCO = 1; % gyrokinetic operator
-ABCO = 1; % INTERSPECIES collisions
-NL_CLOS = 0; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
-SIMID = 'kobayashi_2015_fig1'; % Name of the simulation
-% SIMID = 'debug'; % Name of the simulation
-LINEARITY = 'linear'; % (nonlinear, semilinear, linear)
-% INIT options
-INIT_PHI = 1; % Start simulation with a noisy phi (0= noisy moments 00)
-INIT_ZF = 0; ZF_AMP = 0.0;
-INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 'donothing';
-%% OUTPUTS
-W_DOUBLE = 1;
-W_GAMMA = 1; W_HF = 1;
-W_PHI = 1; W_NA00 = 1;
-W_DENS = 1; W_TEMP = 1;
-W_NAPJ = 1; W_SAPJ = 0;
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% unused
-PMAXE = P; % Highest electron Hermite polynomial degree
-JMAXE = J; % Highest '' Laguerre ''
-PMAXI = P; % Highest ion Hermite polynomial degree
-JMAXI = J; % Highest '' Laguerre ''
-KERN = 0; % Kernel model (0 : GK)
-KX0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
-KPAR = 0.0; % Parellel wave vector component
-LAMBDAD = 0.0;
-kmax = NX*pi/LX;% Highest fourier mode
-HD_CO = 0.5; % Hyper diffusivity cutoff ratio
-MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
-NOISE0 = 0;%1.0e-4;
-BCKGD0 = 1e-4; % Init background
-TAU = 1.0; % e/i temperature ratio
-MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
-MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
-%% Setup and file management
-setup
-system('rm fort*.90');
-outfile = [BASIC.RESDIR,'out.txt'];
-disp(outfile);