diff --git a/matlab/load_params.m b/matlab/load_params.m
index d1eee972ddda5a92ec3e73fa1f2c45a59c2ed955..41b51696a9c91bc7434ad864d1dcb6d6a8db043b 100644
--- a/matlab/load_params.m
+++ b/matlab/load_params.m
@@ -9,8 +9,12 @@ end
ETAB = h5readatt(filename,'/data/input','eta_B');
ETAN = h5readatt(filename,'/data/input','eta_n');
ETAT = h5readatt(filename,'/data/input','eta_T');
-
+PMAXI = h5readatt(filename,'/data/input','pmaxi');
+JMAXI = h5readatt(filename,'/data/input','jmaxi');
+PMAXE = h5readatt(filename,'/data/input','pmaxe');
+JMAXE = h5readatt(filename,'/data/input','jmaxe');
NON_LIN = h5readatt(filename,'/data/input','NON_LIN');
+NU = h5readatt(filename,'/data/input','nu')/0.532;
if NON_LIN == 'y'
NON_LIN = 1;
else
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 920e07aeaa92765acfab67156748b96cf727e0ad..b73460d9672a18138776f8d08b20aed0f6e6da34 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -7,7 +7,6 @@ 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,[' RESTART = ', num2str(BASIC.RESTART),'\n']);
fprintf(fid,[' maxruntime = ', num2str(BASIC.maxruntime),'\n']);
fprintf(fid,'/\n');
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index 7d72fb78872d56c30639fd8ebb3c71957b7ba087..6c6ab63fe803a3eb60d4f60ba76b9ff73db19149 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -1,7 +1,10 @@
%% Load results
-if 1
+if 0
%%
- outfile = '/marconi_scratch/userexternal/ahoffman/HeLaZ/results/Marconi/200x100_L_100_Pe_2_Je_1_Pi_2_Ji_1_nB_0.66_nN_1_nu_1e-01_FC_mu_1e-03/out.txt';
+ outfile ='';
+ outfile ='';
+ outfile ='';
+ outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/Marconi/200x100_L_100_Pe_2_Je_1_Pi_2_Ji_1_nB_0.4_nN_1_nu_1e-01_FC_mu_1e-03/out.txt';
BASIC.RESDIR = load_marconi(outfile);
end
%%
@@ -25,6 +28,7 @@ Nkr = numel(kr); Nkz = numel(kz);
[KZ,KR] = meshgrid(kz,kr);
Lkr = max(kr)-min(kr); Lkz = max(kz)-min(kz);
dkr = Lkr/(Nkr-1); dkz = Lkz/(Nkz-1);
+KPERP2 = KZ.^2+KR.^2;
Lk = max(Lkr,Lkz);
dr = 2*pi/Lk; dz = 2*pi/Lk;
@@ -37,11 +41,13 @@ z = dz*(-Nz/2:(Nz/2-1)); Lz = max(z)-min(z);
disp('Analysis :')
disp('- iFFT')
% IFFT (Lower case = real space, upper case = frequency space)
-ne00 = zeros(Nr,Nz,Ns2D);
+ne00 = zeros(Nr,Nz,Ns2D); % Gyrocenter density
ni00 = zeros(Nr,Nz,Ns2D);
+np_i = zeros(Nr,Nz,Ns5D); % Ion particle density
si00 = zeros(Nr,Nz,Ns5D);
phi = zeros(Nr,Nz,Ns2D);
drphi = zeros(Nr,Nz,Ns2D);
+dr2phi = zeros(Nr,Nz,Ns2D);
dzphi = zeros(Nr,Nz,Ns2D);
for it = 1:numel(Ts2D)
@@ -50,11 +56,20 @@ for it = 1:numel(Ts2D)
ni00(:,:,it) = real(fftshift(ifft2((NI_),Nr,Nz)));
phi (:,:,it) = real(fftshift(ifft2((PH_),Nr,Nz)));
drphi(:,:,it) = real(fftshift(ifft2(1i*KR.*(PH_),Nr,Nz)));
+ dr2phi(:,:,it)= real(fftshift(ifft2(-KR.^2.*(PH_),Nr,Nz)));
dzphi(:,:,it) = real(fftshift(ifft2(1i*KZ.*(PH_),Nr,Nz)));
end
for it = 1:numel(Ts5D)
+ [~, it2D] = min(abs(Ts2D-Ts5D(it)));
si00(:,:,it) = real(fftshift(ifft2(squeeze(Si00(:,:,it)),Nr,Nz)));
+
+ Np_i = zeros(Nkr,Nkz); % Ion particle density in Fourier space
+ for ij = 1:Nji
+ Kn = (KPERP2/2.).^(ij-1) .* exp(-KPERP2/2)/(factorial(ij-1));
+ Np_i = Np_i + Kn.*squeeze(Nipj(1,ij,:,:,it));
+ end
+ np_i(:,:,it) = real(fftshift(ifft2(squeeze(Np_i(:,:)),Nr,Nz)));
end
% Post processing
@@ -62,10 +77,11 @@ disp('- post processing')
E_pot = zeros(1,Ns2D); % Potential energy n^2
E_kin = zeros(1,Ns2D); % Kinetic energy grad(phi)^2
ExB = zeros(1,Ns2D); % ExB drift intensity \propto |\grad \phi|
-Flux_ri = zeros(1,Ns2D); % Particle flux Gamma = <ni drphi>
-Flux_zi = zeros(1,Ns2D); % Particle flux Gamma = <ni dzphi>
-Flux_re = zeros(1,Ns2D); % Particle flux Gamma = <ne drphi>
-Flux_ze = zeros(1,Ns2D); % Particle flux Gamma = <ne dzphi>
+GFlux_ri = zeros(1,Ns2D); % Gyrocenter flux Gamma = <ni drphi>
+GFlux_zi = zeros(1,Ns2D); % Gyrocenter flux Gamma = <ni dzphi>
+GFlux_re = zeros(1,Ns2D); % Gyrocenter flux Gamma = <ne drphi>
+GFlux_ze = zeros(1,Ns2D); % Gyrocenter flux Gamma = <ne dzphi>
+PFlux_ri = zeros(1,Ns5D); % Particle flux
Ne_norm = zeros(Npe,Nje,Ns5D);% Time evol. of the norm of Napj
Ni_norm = zeros(Npi,Nji,Ns5D);% .
Se_norm = zeros(Npe,Nje,Ns5D);% Time evol. of the norm of Sapj
@@ -81,10 +97,10 @@ for it = 1:numel(Ts2D) % Loop over 2D arrays
E_pot(it) = pi/Lr/Lz*sum(sum(abs(NI_).^2))/Nkr/Nkr; % integrate through Parseval id
E_kin(it) = pi/Lr/Lz*sum(sum(abs(Ddr.*PH_).^2+abs(Ddz.*PH_).^2))/Nkr/Nkr;
ExB(it) = max(max(max(abs(phi(3:end,:,it)-phi(1:end-2,:,it))/(2*dr))),max(max(abs(phi(:,3:end,it)-phi(:,1:end-2,it))'/(2*dz))));
- Flux_ri(it) = sum(sum(ni00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz;
- Flux_zi(it) = sum(sum(-ni00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz;
- Flux_re(it) = sum(sum(ne00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz;
- Flux_ze(it) = sum(sum(-ne00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz;
+ GFlux_ri(it) = sum(sum(ni00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz;
+ GFlux_zi(it) = sum(sum(-ni00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz;
+ GFlux_re(it) = sum(sum(ne00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz;
+ GFlux_ze(it) = sum(sum(-ne00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz;
end
E_kin_KZ = mean(mean(abs(Ddr.*PHI(:,:,it)).^2+abs(Ddz.*PHI(:,:,it)).^2,3),2);
@@ -92,12 +108,15 @@ E_kin_KR = mean(mean(abs(Ddr.*PHI(:,:,it)).^2+abs(Ddz.*PHI(:,:,it)).^2,3),2);
dEdt = diff(E_pot+E_kin)./dt2D;
for it = 1:numel(Ts5D) % Loop over 5D arrays
+ [~, it2D] = min(abs(Ts2D-Ts5D(it)));
Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4)/Nkr/Nkz;
Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4)/Nkr/Nkz;
Se_norm(:,:,it)= sum(sum(abs(Sepj(:,:,:,:,it)),3),4)/Nkr/Nkz;
Si_norm(:,:,it)= sum(sum(abs(Sipj(:,:,:,:,it)),3),4)/Nkr/Nkz;
Sne00_norm(it) = sum(sum(abs(Se00(:,:,it))))/Nkr/Nkz;
Sni00_norm(it) = sum(sum(abs(Si00(:,:,it))))/Nkr/Nkz;
+ % Particle flux
+ PFlux_ri(it) = sum(sum(np_i(:,:,it).*dzphi(:,:,it2D)))*dr*dz/Lr/Lz;
end
%% Compute growth rate
@@ -147,7 +166,6 @@ set(gcf, 'Position', [100, 100, 900, 800])
plot(Ts2D,Flux_re,'-','DisplayName','$\Gamma_{re}$')
plot(Ts2D,Flux_ze,'-','DisplayName','$\Gamma_{ze}$')
grid on; xlabel('$t c_s/R$'); ylabel('$\Gamma$'); %legend('show');
-if strcmp(OUTPUTS.write_non_lin,'.true.')
subplot(224)
for ip = 1:Npi
for ij = 1:Nji
@@ -157,15 +175,6 @@ if strcmp(OUTPUTS.write_non_lin,'.true.')
end
end
grid on; xlabel('$t c_s/R$'); ylabel('$\sum_{k_r,k_z}|S_i^{pj}|$'); %legend('show');
-else
-%% Growth rate
- subplot(224)
- [~,ikr0KH] = min(abs(kr-KR0KH));
- plot(kz(1:Nz/2),gkr0kz_Ni00(1:Nz/2),...
- 'DisplayName',['J = ',num2str(JMAXI)]);
- xlabel('$k_z$'); ylabel('$\gamma_{Ni}$');
- xlim([0. 1.0]); %ylim([0. 0.04]);
-end
suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB)]);
save_figure
end
@@ -236,7 +245,7 @@ save_figure
end
%%
-t0 = 40;
+t0 = 0;
skip_ = 1;
DELAY = 0.01*skip_;
FRAMES = floor(t0/(Ts2D(2)-Ts2D(1)))+1:skip_:numel(Ts2D);
@@ -256,13 +265,6 @@ FIELDNAME = '$n_e$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
create_gif
end
if 0
-%% Density ion - electron
-GIFNAME = ['ni-ne',sprintf('_%.2d',JOBNUM)]; INTERP = 1;
-FIELD = real(ni00+ne00); X = RR; Y = ZZ; T = Ts2D;
-FIELDNAME = '$n_i-n_e$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
-create_gif
-end
-if 0
%% Phi
GIFNAME = ['phi',sprintf('_%.2d',JOBNUM)];INTERP = 1;
FIELD = real(phi); X = RR; Y = ZZ; T = Ts2D;
@@ -270,6 +272,14 @@ FIELDNAME = '$\phi$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
create_gif
end
if 0
+%% phi @ z = 0
+GIFNAME = ['phi_r0',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
+FIELD =(squeeze(real(phi(:,1,:)))); linestyle = '-.';
+X = (r); T = Ts2D; YMIN = -1.1; YMAX = 1.1; XMIN = min(r); XMAX = max(r);
+FIELDNAME = '$\phi(r=0)$'; XNAME = '$r/\rho_s$';
+create_gif_1D
+end
+if 0
%% Density ion frequency
GIFNAME = ['Ni00',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
FIELD =ifftshift((abs(Ni00)),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts2D;
@@ -288,18 +298,22 @@ end
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-t0 = 1; t1 = 400; [~,it0] = min(abs(t0-Ts2D)); [~,it1] = min(abs(t1-Ts2D));
fig = figure; FIGNAME = 'space_time_drphi';set(gcf, 'Position', [100, 100, 1200, 400])
subplot(211)
- plot(Ts2D(it0:it1),Flux_ri(it0:it1));
+ yyaxis left
+ plot(Ts2D,GFlux_ri); hold on
+ plot(Ts5D,PFlux_ri,'o');
ylabel('$\Gamma_r$'); grid on
-% title(['$\eta=',num2str(ETAB),'\quad',...
-% '\nu_{',CONAME,'}=',num2str(NU),'$'])
-% legend(['$P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$'])
- ylim([0,1.1*max(Flux_ri(it0:it1))]);
+ ylim([0,1.1*max(GFlux_ri)]);
+ yyaxis right
+ plot(Ts2D,squeeze(mean(max(dr2phi))))
+ ylabel('$s\sim\langle\partial_r^2\phi\rangle_z$'); grid on
+ title(['$\eta=',num2str(ETAB),'\quad',...
+ '\nu_{',CONAME,'}=',num2str(NU),'$'])
+ legend(['$P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$'])
subplot(212)
- [TY,TX] = meshgrid(r,Ts2D(it0:it1));
- pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,it0:it1),2))'); set(pclr, 'edgecolor','none'); %colorbar;
+ [TY,TX] = meshgrid(r,Ts2D);
+ pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,:),2))'); set(pclr, 'edgecolor','none'); %colorbar;
xlabel('$t c_s/R$'), ylabel('$r/\rho_s$')
legend('$\langle\partial_r \phi\rangle_z$')
save_figure
diff --git a/wk/debug_script.m b/wk/debug_script.m
index da2f256042484c0b0eeea6c89c7e13ad14722b3d..5aa5d036a7b0db39db3724f8273e24c298db211d 100644
--- a/wk/debug_script.m
+++ b/wk/debug_script.m
@@ -5,9 +5,9 @@ addpath(genpath('../matlab')) % ... add
CLUSTER.TIME = '24:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 1e2; % Collision frequency
+NU = 0e2; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-ETAB = 0.5; % Magnetic gradient
+ETAB = 0.4; % Magnetic gradient
ETAN = 1.0; % Density gradient
ETAT = 0.0; % Temperature gradient
HD_CO = 0.5; % Hyper diffusivity cutoff ratio
@@ -20,7 +20,7 @@ JMAXE = 1; % Highest '' Laguerre ''
PMAXI = 2; % Highest ion Hermite polynomial degree
JMAXI = 1; % Highest '' Laguerre ''
%% TIME PARAMETERS
-TMAX = 200; % Maximal time unit
+TMAX = 100; % Maximal time unit
DT = 1e-2; % Time step
SPS0D = 1/DT; % Sampling per time unit for profiler
SPS2D = 2; % Sampling per time unit for 2D arrays
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 88d87fd4e76d65c923d704322963965a53635d19..b2b28af60d5b7426ab43a3520266b0f45af36dd7 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -11,21 +11,18 @@ CLUSTER.NTPN = '24'; % N tasks per node
CLUSTER.PART = 'prod'; % dbg or prod
CLUSTER.MEM = '16GB'; % Memory
CLUSTER.JNAME = 'gamma_inf'; % Job name
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
NU = 1e-1; % Collision frequency
-TAU = 1.0; % e/i temperature ratio
-ETAB = 0.66; % Magnetic gradient
-ETAN = 1.0; % Density gradient
-ETAT = 0.0; % Temperature gradient
-HD_CO = 0.5; % Hyper diffusivity cutoff ratio
-NOISE0 = 1.0e-5;
+ETAB = 0.66; % Magnetic gradient
+NU_HYP = 0.1; % Hyperdiffusivity coefficient
%% GRID PARAMETERS
N = 200; % Frequency gridpoints (Nkr = N/2)
L = 100; % Size of the squared frequency domain
-P = 4; % Electron and Ion highest Hermite polynomial degree
-J = 2; % Electron and Ion highest Laguerre polynomial degree
+P = 8; % Electron and Ion highest Hermite polynomial degree
+J = 4; % Electron and Ion highest Laguerre polynomial degree
%% TIME PARAMETERS
-TMAX = 500; % Maximal time unit
+TMAX = 400; % Maximal time unit
DT = 1e-2; % Time step
SPS0D = 10; % Sampling per time unit for profiler
SPS2D = 1; % Sampling per time unit for 2D arrays
@@ -34,11 +31,13 @@ SPSCP = 1/10; % Sampling per time unit for checkpoints
RESTART = 0; % To restart from last checkpoint
JOB2LOAD= 0;
%% OPTIONS
-SIMID = 'Marconi'; % Name of the simulation
+SIMID = 'Marconi_new_AA'; % Name of the simulation
CO = -1; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% unused
+%% fixed parameters (for current study)
KR0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
NO_E = 0; % Remove electrons dynamic
% DK = 0; % Drift kinetic model (put every kernel_n to 0 except n=0 to 1)
@@ -51,7 +50,12 @@ JMAXE = J; % Highest '' Laguerre ''
PMAXI = P; % Highest ion Hermite polynomial degree
JMAXI = J; % Highest '' Laguerre ''
kmax = N*pi/L;% Highest fourier mode
-MU = 0.1/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient
+HD_CO = 0.5; % Hyper diffusivity cutoff ratio
+MU = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient
+NOISE0 = 1.0e-5;
+ETAT = 0.0; % Temperature gradient
+ETAN = 1.0; % Density gradient
+TAU = 1.0; % e/i temperature ratio
%% Run following scripts
setup
diff --git a/wk/setup.m b/wk/setup.m
index f11e2da7aa8ac95d35543b485189111ab1cbb698..7e4ec8eb3c1917754a53709a3dc79a94ef1ad029 100644
--- a/wk/setup.m
+++ b/wk/setup.m
@@ -39,18 +39,22 @@ INITIAL.only_Na00 = '.false.';
INITIAL.initback_moments = 0.0e-5;
INITIAL.initnoise_moments = NOISE0;
INITIAL.iseed = 42;
+fcmat_pmaxe = 25;
+fcmat_jmaxe = 18;
+fcmat_pmaxi = 25;
+fcmat_jmaxi = 18;
INITIAL.selfmat_file = ...
- ['''../../../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_',num2str(GRID.pmaxe),...
- '_Jmaxe_',num2str(GRID.jmaxe),'_Pmaxi_',num2str(GRID.pmaxi),'_Jmaxi_',...
- num2str(GRID.jmaxi),'_pamaxx_10.h5'''];
+ ['''../../../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_',num2str(fcmat_pmaxe),...
+ '_Jmaxe_',num2str(fcmat_jmaxe),'_Pmaxi_',num2str(fcmat_pmaxi),'_Jmaxi_',...
+ num2str(fcmat_jmaxi),'_pamaxx_10.h5'''];
INITIAL.eimat_file = ...
- ['''../../../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_',num2str(GRID.pmaxe),...
- '_Jmaxe_',num2str(GRID.jmaxe),'_Pmaxi_',num2str(GRID.pmaxi),'_Jmaxi_',...
- num2str(GRID.jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5'''];
+ ['''../../../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_',num2str(fcmat_pmaxe),...
+ '_Jmaxe_',num2str(fcmat_jmaxe),'_Pmaxi_',num2str(fcmat_pmaxi),'_Jmaxi_',...
+ num2str(fcmat_jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5'''];
INITIAL.iemat_file = ...
- ['''../../../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_',num2str(GRID.pmaxe),...
- '_Jmaxe_',num2str(GRID.jmaxe),'_Pmaxi_',num2str(GRID.pmaxi),'_Jmaxi_',...
- num2str(GRID.jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5'''];
+ ['''../../../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_',num2str(fcmat_pmaxe),...
+ '_Jmaxe_',num2str(fcmat_jmaxe),'_Pmaxi_',num2str(fcmat_pmaxi),'_Jmaxi_',...
+ num2str(fcmat_jmaxi),'_pamaxx_10_tau_1.0000_mu_0.0233.h5'''];
% Naming and creating input file
if (CO == 0); CONAME = 'LB';
elseif(CO == -1); CONAME = 'FC';