From 5b5157b79b2fb9c3e81f8ce91a41193a0e8ee79d Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Mon, 5 Oct 2020 14:52:06 +0200
Subject: [PATCH] HeLaZ can now run in 1D over kz for growth rate study as MOLI
---
matlab/create_gif_1D.m | 51 +++++++++
src/grid_mod.F90 | 44 ++++---
wk/analysis_1D.m | 218 +++++++++++++++++++++++++++++++++++
wk/analysis_2D.m | 252 +++++++++++++++++++++++++++++++++++++++++
wk/setup.m | 82 +++++++-------
5 files changed, 591 insertions(+), 56 deletions(-)
create mode 100644 matlab/create_gif_1D.m
create mode 100644 wk/analysis_1D.m
create mode 100644 wk/analysis_2D.m
diff --git a/matlab/create_gif_1D.m b/matlab/create_gif_1D.m
new file mode 100644
index 00000000..d303e286
--- /dev/null
+++ b/matlab/create_gif_1D.m
@@ -0,0 +1,51 @@
+title1 = GIFNAME;
+FIGDIR = ['../results/',BASIC.SIMID,'/'];
+if ~exist(FIGDIR, 'dir')
+ mkdir(FIGDIR)
+end
+
+GIFNAME = [FIGDIR, GIFNAME,'.gif'];
+
+% Set colormap boundaries
+hmax = max(max(max(FIELD)));
+hmin = min(min(min(FIELD)));
+flag = 0;
+if hmax == hmin
+ disp('Warning : h = hmin = hmax = const')
+else
+% Setup figure frame
+fig = figure;
+ plot(X,FIELD(:,1)); % to set up
+ axis tight manual % this ensures that getframe() returns a consistent size
+ xlabel(XNAME); ylabel(FIELDNAME);
+
+ in = 1;
+ nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:)));
+ for n = FRAMES % loop over selected frames
+ scale = max(FIELD(:,n));
+ plot(X,FIELD(:,n)/scale);
+ ylim([YMIN,YMAX]); xlim([XMIN,XMAX]);
+ title(['$t \approx$', sprintf('%.3d',ceil(T(n))), ', scaling = ',sprintf('%.1e',scale)]);
+ drawnow
+ % Capture the plot as an image
+ frame = getframe(fig);
+ im = frame2im(frame);
+ [imind,cm] = rgb2ind(im,64);
+ % Write to the GIF File
+ if in == 1
+ imwrite(imind,cm,GIFNAME,'gif', 'Loopcount',inf);
+ else
+ imwrite(imind,cm,GIFNAME,'gif','WriteMode','append', 'DelayTime',DELAY);
+ end
+ % terminal info
+ while nbytes > 0
+ fprintf('\b')
+ nbytes = nbytes - 1;
+ end
+ nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,:)));
+ in = in + 1;
+ end
+ disp(' ')
+ disp(['Gif saved @ : ',GIFNAME])
+end
+
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 6fee2179..7845b73b 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -94,10 +94,14 @@ CONTAINS
ENDIF
! Discretized r positions ordered as dx*(0 1 2 -3 -2 -1)
ALLOCATE(rarray(irs:ire))
- DO ir = irs,ire
- rarray(ir) = deltar*(MODULO(ir-1,Nr/2)-Nr/2*FLOOR(2.*real(ir-1)/real(Nr)))
- END DO
- rarray(Nr/2+1) = -rarray(Nr/2+1)
+ IF (NR .GT. 1) THEN
+ DO ir = irs,ire
+ rarray(ir) = deltar*(MODULO(ir-1,Nr/2)-Nr/2*FLOOR(2.*real(ir-1)/real(Nr)))
+ END DO
+ rarray(Nr/2+1) = -rarray(Nr/2+1)
+ ELSE
+ rarray(1) = 0._dp
+ ENDIF
END SUBROUTINE set_rgrid
@@ -134,21 +138,29 @@ CONTAINS
ikrs = 1
ikre = Nkr
! Grid spacings
- deltakr = 2._dp*PI/Nr/deltar
-
+ IF ( Nr .GT. 1 ) THEN ! To avoid case with 0 intervals
+ deltakr = 2._dp*PI/Nr/deltar
+ ELSE
+ deltakr = 0;
+ ENDIF
! Discretized kr positions ordered as dk*(0 1 2)
ALLOCATE(krarray(ikrs:ikre))
- DO ikr = ikrs,ikre
- ! krarray(ikr) = REAL(ikr-1,dp) * deltakr
- krarray(ikr) = deltakr*(MODULO(ikr-1,Nkr/2)-Nkr/2*FLOOR(2.*real(ikr-1)/real(Nkr)))
- if (krarray(ikr) .EQ. 0) THEN
- ikr_0 = ikr
- ENDIF
- END DO
- krarray(Nr/2+1) = -krarray(Nr/2+1)
+ IF ( Nr .GT. 1 ) THEN
+ DO ikr = ikrs,ikre
+ ! krarray(ikr) = REAL(ikr-1,dp) * deltakr
+ krarray(ikr) = deltakr*(MODULO(ikr-1,Nkr/2)-Nkr/2*FLOOR(2.*real(ikr-1)/real(Nkr)))
+ if (krarray(ikr) .EQ. 0) THEN
+ ikr_0 = ikr
+ ENDIF
+ END DO
+ krarray(Nr/2+1) = -krarray(Nr/2+1)
+ ELSE
+ krarray(1) = 0._dp
+ ikr_0 = 1
+ ENDIF
! Orszag 2/3 filter
- two_third_krmax = 2._dp/3._dp*deltakr*Nkr
+ two_third_krmax = 2._dp/3._dp*deltakr*(Nkr/2)
ALLOCATE(AA_r(ikrs:ikre))
DO ikr = ikrs,ikre
IF ( (krarray(ikr) .GT. -two_third_krmax) .AND. (krarray(ikr) .LT. two_third_krmax) ) THEN
@@ -183,7 +195,7 @@ CONTAINS
! kzarray(Nz/2+1) = -kzarray(Nz/2+1)
! Orszag 2/3 filter
- two_third_kzmax = 2._dp/3._dp*deltakz*(Nkz/2);
+ two_third_kzmax = 2._dp/3._dp*deltakz*Nkz;
ALLOCATE(AA_z(ikzs:ikze))
DO ikz = ikzs,ikze
IF ( (kzarray(ikz) .GT. -two_third_kzmax) .AND. (kzarray(ikz) .LT. two_third_kzmax) ) THEN
diff --git a/wk/analysis_1D.m b/wk/analysis_1D.m
new file mode 100644
index 00000000..bcc714c8
--- /dev/null
+++ b/wk/analysis_1D.m
@@ -0,0 +1,218 @@
+default_plots_options
+%% load results
+JOBNUM = 00;
+filename = [BASIC.SIMID,'_','%.2d.h5'];
+filename = sprintf(filename,JOBNUM); disp(['Analysing ',filename])
+[Nipj, p_, j_, kr, kz, Ts, dt] = load_5D_data(filename, 'moments_i');
+Nepj = load_5D_data(filename, 'moments_e');
+Ni = squeeze(Nipj(1,1,:,:,:));
+Ne = squeeze(Nepj(1,1,:,:,:));
+PH = squeeze(load_2D_data(filename, 'phi'));
+Ts = Ts';
+Ns = numel(Ts);
+dt_samp = mean(diff(Ts));
+% Build grids
+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);
+Lk = max(Lkr,Lkz);
+
+dr = 2*pi/Lk; dz = 2*pi/Lk;
+Nr = max(Nkr,Nkz) * (Nkr > 1) + (Nkr == 1); Nz = Nkz;
+r = dr*(-Nr/2:(Nr/2-1)); Lr = max(r)-min(r);
+z = dz*(-Nz/2:(Nz/2-1)); Lz = max(z)-min(z);
+[YY,XX] = meshgrid(z,r);
+% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% IFFT
+ne = zeros(Nz, Ns);
+ni = zeros(Nz, Ns);
+phi = zeros(Nz, Ns);
+
+for it = 1:numel(PH(1,:))
+ NE_ = Ne(:,it); NN_ = Ni(:,it); PH_ = PH(:,it);
+ F_ = (ifft((NE_),Nz));
+ ne(:,it)= real(fftshift(F_));
+ F_ = (ifft((NN_),Nz));
+ ni(:,it) = real(fftshift(F_));
+ F_ = (ifft((PH_),Nz));
+ phi(:,it) = real(fftshift(F_));
+end
+
+% Post processing
+ne_00 = zeros(1,Ns); % Time evolution of ne(r,z) at origin
+Ne_gm = zeros(1,Ns); % Time evolution of Ne(k) max gamma (max real)
+ni_00 = zeros(1,Ns); % .
+Ni_gm = zeros(1,Ns); % .
+phi_00 = zeros(1,Ns); % .
+E_pot = zeros(1,Ns); % Potential energy n^2
+E_kin = zeros(1,Ns); % Kinetic energy grad(phi)^2
+ExB = zeros(1,Ns); % ExB drift intensity \propto |\grad \phi|
+CFL = zeros(1,Ns); % CFL time step
+Ddz = 1i*kz;
+[~,iz0] = min(abs(z)); % index of z==0
+[~,ikz1] = min(abs(kz-round(1/dkz)*dkz)); % index of kz==1
+for it = 1:numel(PH(1,:))
+ NE_ = squeeze(Ne(:,it)); NN_ = squeeze(Ni(:,it)); PH_ = squeeze(PH(:,it));
+
+ ne_00(it) = ne(iz0,it);
+ Ne_gm(it) = max(real(Ne(:,it)));
+
+ ni_00(it) = ni(iz0,it);
+ Ni_gm(it) = max(real(Ni(:,it)));
+
+ phi_00(it) = phi(iz0,it);
+
+ E_pot(it) = pi/Lz*sum(abs(NN_).^2)/Nkz; % integrate through Parseval id
+ E_kin(it) = pi/Lz*sum(abs(Ddz.*PH_).^2)/Nkz;
+ ExB(it) = max(abs(phi(3:end,it)-phi(1:end-2,it))'/(2*dz));
+end
+E_kin_KZ = abs(Ddz.*PH(:,it)).^2;
+dEdt = diff(E_pot+E_kin)./diff(Ts);
+
+% Growth rate
+gamma_Ne = zeros(1,Nkz);
+gamma_Ni = zeros(1,Nkz);
+gamma_PH = zeros(1,Nkz);
+for ikz = 1:Nkz
+ gamma_Ne(ikz) = real(LinearFit_s(Ts,Ne(ikz,:)));
+ gamma_Ni(ikz) = real(LinearFit_s(Ts,Ni(ikz,:)));
+ gamma_PH(ikz) = real(LinearFit_s(Ts,PH(ikz,:)));
+end
+
+%% PLOTS
+%% Time evolutions
+fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM)];
+ subplot(221)
+ semilogy(Ts,abs(ne_00),'-','DisplayName','$n_e^{00}$'); hold on;
+ semilogy(Ts,abs(ni_00),'-','DisplayName','$n_i^{00}$');
+ grid on; xlabel('$t$'); ylabel('$|n_a(x=0,y=0)|$');
+ subplot(222)
+ semilogy(Ts,abs(Ni_gm),'-','DisplayName','$\phi$')
+ grid on; xlabel('$t$'); ylabel('$|\tilde n(k_r\approx 1,k_z\approx 1)|$');
+ subplot(223)
+ semilogy(Ts,E_kin+E_pot,'-','DisplayName','$\sum|ik\tilde\phi_i|^2+\sum|\tilde n_i|^2$')
+ hold on;
+ grid on; xlabel('$t$'); ylabel('$E$'); legend('show');
+FMT = '.fig'; save_figure
+
+%% Growth rate
+fig = figure; FIGNAME = ['growth_rate',sprintf('_%.2d',JOBNUM)];
+ subplot(221)
+ plot(kz,gamma_Ne,'-'); hold on;
+ grid on; xlabel('$k_z$'); ylabel('$\gamma(N_e^{00})$');
+ subplot(222)
+ plot(kz,gamma_Ni,'-'); hold on;
+ grid on; xlabel('$k_z$'); ylabel('$\gamma(N_i^{00})$');
+ subplot(223)
+ plot(kz,gamma_PH,'-'); hold on;
+ grid on; xlabel('$k_z$'); ylabel('$\gamma(\tilde\phi)$'); legend('show');
+FMT = '.fig'; save_figure
+
+%% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+DELAY = 0.07; skip_ = 10;
+FRAMES = 1:skip_:numel(Ts);
+if 0
+%% Density electron
+GIFNAME = ['ne',sprintf('_%.2d',JOBNUM)];
+FIELD = real(ne); X = z; T = Ts;
+FIELDNAME = '$n_e^{00}$'; XNAME = '$z$';
+XMIN = -20; XMAX = 20; YMIN = -0.6; YMAX = 1.1;
+create_gif_1D
+%% Density ion
+GIFNAME = ['ni',sprintf('_%.2d',JOBNUM)];
+FIELD = real(ni); X = z; T = Ts;
+FIELDNAME = '$n_i^{00}$'; XNAME = '$z$';
+XMIN = -20; XMAX = 20; YMIN = -0.6; YMAX = 1.1;
+create_gif_1D
+%% Phi
+GIFNAME = ['phi',sprintf('_%.2d',JOBNUM)];
+FIELD = real(phi); X = z; T = Ts;
+FIELDNAME = '$\phi$'; XNAME = '$z$';
+XMIN = -20; XMAX = 20; YMIN = -1.1; YMAX = 1.1;
+create_gif_1D
+%% Density electron frequency
+GIFNAME = ['Ni',sprintf('_%.2d',JOBNUM)];
+FIELD = fftshift(real(Ni)); X = kr\z; T = Ts;
+FIELDNAME = '$N_i^{00}$'; XNAME = '$k_z$';
+create_gif_1D
+end
+
+if 0
+%% Space-Time diagrams at r = 0
+plt = @(x) real(x);
+fig = figure; FIGNAME = ['z_space_time_diag',sprintf('_%.2d',JOBNUM)];
+ [TY,TX] = meshgrid(Ts,z);
+subplot(221)% density
+ pclr = pcolor(TX,TY,(plt(ne))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$z\,(r=0)$'); ylabel('$t$'); title('$n_e^{00}$');
+subplot(222)% density
+ pclr = pcolor(TX,TY,(plt(ni))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$z\,(r=0)$'); ylabel('$t$'); title('$n_i^{00}$');
+subplot(223)% density
+ pclr = pcolor(TX,TY,(plt(phi))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$z\,(r=0)$'); ylabel('$t$'); title('$\phi$');
+FMT = '.fig'; save_figure
+
+%% Space-Time diagrams at kr=0
+plt = @(x) abs(x);
+fig = figure; FIGNAME = ['kz_space_time_diag',sprintf('_%.2d',JOBNUM)];
+ [TY,TX] = meshgrid(Ts,kz);
+subplot(221)% density
+ pclr = pcolor(TX,TY,(plt(Ne))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$kz$'); ylabel('$t$'); title('$\textrm{Re}(N_e^{00})|_{k_r=0}$');
+subplot(222)% density
+ pclr = pcolor(TX,TY,(plt(Ni))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$kz$'); ylabel('$t$'); title('$\textrm{Re}(N_i^{00})|_{k_r=0}$');
+subplot(223)% density
+ pclr = pcolor(TX,TY,(plt(PH))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$kz$'); ylabel('$t$'); title('$\textrm{Re}(\tilde\phi)|_{k_r=0}$');
+FMT = '.fig'; save_figure
+end
+
+if 0
+%% Mode time evolution
+[~,ik05] = min(abs(kz-0.50));
+[~,ik10] = min(abs(kz-1.00));
+[~,ik15] = min(abs(kz-1.50));
+[~,ik20] = min(abs(kz-2.00));
+
+fig = figure; FIGNAME = ['frame',sprintf('_%.2d',JOBNUM)];
+ subplot(221); plt = @(x) (abs(x));
+ semilogy(Ts,plt(PH(ik05,:)),'DisplayName', '$k_z = 0.5$'); hold on
+ semilogy(Ts,plt(PH(ik10,:)),'DisplayName', '$k_z = 1.0$')
+ semilogy(Ts,plt(PH(ik15,:)),'DisplayName', '$k_z = 1.5$')
+ semilogy(Ts,plt(PH(ik20,:)),'DisplayName', '$k_z = 2.0$')
+ xlabel('$t$'); ylabel('$\hat\phi$'); legend('show');
+ semilogy(Ts,abs(PH(ik05,end)).*exp(gamma_PH(ik05).*(Ts-Ts(end))),'--k')
+ semilogy(Ts,abs(PH(ik10,end)).*exp(gamma_PH(ik10).*(Ts-Ts(end))),'--k')
+ subplot(222); plt = @(x) (abs(x));
+ semilogy(Ts,plt(Ni(ik05,:)),'DisplayName', '$k_z = 0.5$'); hold on
+ semilogy(Ts,plt(Ni(ik10,:)),'DisplayName', '$k_z = 1.0$')
+ semilogy(Ts,plt(Ni(ik15,:)),'DisplayName', '$k_z = 1.5$')
+ semilogy(Ts,plt(Ni(ik20,:)),'DisplayName', '$k_z = 2.0$')
+ xlabel('$t$'); ylabel('$\hat n_i^{00}$'); legend('show');
+ subplot(223); plt = @(x) (abs(x));
+ semilogy(Ts,plt(Ne(ik05,:)),'DisplayName', '$k_z = 0.5$'); hold on
+ semilogy(Ts,plt(Ne(ik10,:)),'DisplayName', '$k_z = 1.0$')
+ semilogy(Ts,plt(Ne(ik15,:)),'DisplayName', '$k_z = 1.5$')
+ semilogy(Ts,plt(Ne(ik20,:)),'DisplayName', '$k_z = 2.0$')
+ xlabel('$t$'); ylabel('$\hat n_e^{00}$'); legend('show');
+FMT = '.fig'; save_figure
+end
+
+if 0
+%% Show frame
+it = min(70,numel(Ts));
+fig = figure; FIGNAME = ['frame',sprintf('_%.2d',JOBNUM)];
+ subplot(221); plt = @(x) (abs(x));
+ plot(kz,plt(PH(:,it)))
+ xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts(it))); legend('$\hat\phi$');
+ subplot(222); plt = @(x) (abs(x));
+ plot(kz,plt(Ni(:,it)))
+ xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts(it))); legend('$\hat n_i^{00}$');
+ subplot(223); plt = @(x) (abs(x));
+ plot(kz,plt(Ne(:,it)))
+ xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts(it))); legend('$\hat n_e^{00}$');
+FMT = '.fig'; save_figure
+end
\ No newline at end of file
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
new file mode 100644
index 00000000..25ad5d17
--- /dev/null
+++ b/wk/analysis_2D.m
@@ -0,0 +1,252 @@
+%% load results
+JOBNUM = 00;
+filename = [BASIC.SIMID,'_','%.2d.h5'];
+filename = sprintf(filename,JOBNUM); disp(['Analysing ',filename])
+[Nipj, p_, j_, kr, kz, Ts, dt] = load_5D_data(filename, 'moments_i');
+Nepj = load_5D_data(filename, 'moments_e');
+Ni = squeeze(Nipj(1,1,:,:,:));
+Ne = squeeze(Nepj(1,1,:,:,:));
+PH = load_2D_data(filename, 'phi');
+Ts = Ts';
+Ns = numel(Ts);
+dt_samp = mean(diff(Ts));
+if strcmp(OUTPUTS.write_non_lin,'.true.')
+ Sipj = load_5D_data(filename, 'Sipj');
+ Sepj = load_5D_data(filename, 'Sepj');
+ SNi = squeeze(Sipj(1,1,:,:,:));
+ SNe = squeeze(Sepj(1,1,:,:,:));
+end
+% Build grids
+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);
+Lk = max(Lkr,Lkz);
+
+dr = 2*pi/Lk; dz = 2*pi/Lk;
+Nr = max(Nkr,Nkz); Nz = Nr;
+r = dr*(-Nr/2:(Nr/2-1)); Lr = max(r)-min(r);
+z = dz*(-Nz/2:(Nz/2-1)); Lz = max(z)-min(z);
+[YY,XX] = meshgrid(z,r);
+% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% IFFT
+ne = zeros(Nr,Nz);
+ni = zeros(Nr,Nz);
+phi = zeros(Nr,Nz);
+
+for it = 1:numel(PH(1,1,:))
+ NE_ = Ne(:,:,it); NN_ = Ni(:,:,it); PH_ = PH(:,:,it);
+ F_ = (ifft2((NE_),Nr,Nz));
+ ne(:,:,it)= real(fftshift(F_));
+ F_ = (ifft2((NN_),Nr,Nz));
+ ni(:,:,it) = real(fftshift(F_));
+ F_ = (ifft2((PH_),Nr,Nz));
+ phi(:,:,it) = real(fftshift(F_));
+end
+
+% Post processing
+Ne_ST_kr = zeros(Nkr,Ns); % Space-Time of max_kz Ne(k)
+Ne_ST_kz = zeros(Nkz,Ns); % '' max_kr Ne(k)
+ne_ST_r = zeros(Nr,Ns); % Space-Time of ne(z==0)
+ne_ST_z = zeros(Nz,Ns); % '' r==0
+ne_00 = zeros(1,Ns); % Time evolution of ne(r,z) at origin
+Ne_gm = zeros(1,Ns); % Time evolution of Ne(k) max gamma (max real)
+Ni_ST_kr = zeros(Nkr,Ns); % same for ions
+Ni_ST_kz = zeros(Nkz,Ns); % .
+ni_ST_r = zeros(Nr,Ns); % .
+ni_ST_z = zeros(Nz,Ns); % .
+ni_00 = zeros(1,Ns); % .
+Ni_gm = zeros(1,Ns); % .
+PH_ST_kr = zeros(Nkr,Ns); % same for ES-potential
+PH_ST_kz = zeros(Nkz,Ns); % .
+phi_ST_r = zeros(Nr,Ns); % .
+phi_ST_z = zeros(Nz,Ns); % .
+phi_00 = zeros(1,Ns); % .
+Sne_norm = zeros(1,Ns); % Time evolution of the amp of e nonlin term
+Sni_norm = zeros(1,Ns); % Time evolution of the amp of i nonlin term
+E_pot = zeros(1,Ns); % Potential energy n^2
+E_kin = zeros(1,Ns); % Kinetic energy grad(phi)^2
+ExB = zeros(1,Ns); % ExB drift intensity \propto |\grad \phi|
+CFL = zeros(1,Ns); % CFL time step
+Ddr = 1i*KR; Ddz = 1i*KZ; lapl = Ddr.^2 + Ddz.^2;
+[~,ir0] = min(abs(r)); % index of r==0
+[~,iz0] = min(abs(z)); % index of z==0
+[~,ikr1] = min(abs(kr-round(1/dkr)*dkr)); % index of kr==1
+[~,ikz1] = min(abs(kz-round(1/dkz)*dkz)); % index of kz==1
+for it = 1:numel(PH(1,1,:))
+ NE_ = Ne(:,:,it); NN_ = Ni(:,:,it); PH_ = PH(:,:,it);
+
+ ne_ST_r (:,it)= ne(:,iz0,it); ne_ST_z (:,it)= ne(ir0,:,it);
+ Ne_ST_kr (:,it)= max(real(Ne(:,:,it)),[],2);
+ Ne_ST_kz (:,it)= max(real(Ne(:,:,it)),[],1);
+ ne_00(it) = ne(ir0,iz0,it);
+ Ne_gm(it) = max(max(real(Ne(:,:,it))));
+
+ ni_ST_r (:,it)= ni(:,iz0,it); ni_ST_z (:,it)= ni(ir0,:,it);
+ Ni_ST_kr (:,it)= max(real(Ni(:,:,it)),[],2);
+ Ni_ST_kz (:,it)= max(real(Ni(:,:,it)),[],1);
+ ni_00(it) = ni(ir0,iz0,it);
+ Ni_gm(it) = max(max(real(Ni(:,:,it))));
+
+ phi_ST_r(:,it) = phi(:,iz0,it); phi_ST_z(:,it) = phi(ir0,:,it);
+ PH_ST_kr(:,it) = max(real(PH(:,:,it)),[],2);
+ PH_ST_kz(:,it) = max(real(PH(:,:,it)),[],1);
+ phi_00(it) = phi(ir0,iz0,it);
+
+if strcmp(OUTPUTS.write_non_lin,'.true.')
+ Sni_norm(it) = sum(sum(abs(SNi(:,:,it))));
+ Sne_norm(it) = sum(sum(abs(SNe(:,:,it))));
+end
+ E_pot(it) = pi/Lr/Lz*sum(sum(abs(NN_).^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))));
+end
+E_kin_KZ = mean(mean(abs(Ddr.*PH(:,:,it)).^2+abs(Ddz.*PH(:,:,it)).^2,3),1);
+E_kin_KR = mean(mean(abs(Ddr.*PH(:,:,it)).^2+abs(Ddz.*PH(:,:,it)).^2,3),2);
+dEdt = diff(E_pot+E_kin)./diff(Ts);
+%% PLOTS
+%% Time evolutions
+fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM)];
+ subplot(221)
+ semilogy(Ts,abs(ne_00),'-','DisplayName','$n_e^{00}$'); hold on;
+ semilogy(Ts,abs(ni_00),'-','DisplayName','$n_i^{00}$');
+ grid on; xlabel('$t$'); ylabel('$|n_a(x=0,y=0)|$');
+ subplot(222)
+ semilogy(Ts,abs(Ni_gm),'-','DisplayName','$\phi$')
+ grid on; xlabel('$t$'); ylabel('$|\tilde n(k_r\approx 1,k_z\approx 1)|$');
+ subplot(223)
+ semilogy(Ts,E_kin+E_pot,'-','DisplayName','$\sum|ik\tilde\phi_i|^2+\sum|\tilde n_i|^2$')
+ hold on;
+ grid on; xlabel('$t$'); ylabel('$E$'); legend('show');
+if strcmp(OUTPUTS.write_non_lin,'.true.')
+ subplot(224)
+ semilogy(Ts,Sne_norm,'-','DisplayName','$\sum|S_e^{00}|$');
+ hold on;
+ semilogy(Ts,Sni_norm,'-','DisplayName','$\sum|S_i^{00}|$');
+ grid on; xlabel('$t$'); ylabel('$S$'); legend('show');
+end
+FMT = '.fig'; save_figure
+
+%% Spectra energy
+% fig = figure; FIGNAME = ['Energy_kin_KZ',sprintf('_%.2d',JOBNUM)];
+% semilogy(kr(floor(end/2)+1:end),E_kin_KR(floor(end/2)+1:end),'o','DisplayName','$\sum_y\langle|ik\tilde\phi_i|^2\rangle_t$')
+% hold on;
+% loglog(kz(floor(end/2)+1:end),E_kin_KZ(floor(end/2)+1:end),'o','DisplayName','$\sum_x\langle|ik\tilde\phi_i|^2\rangle_t$')
+% grid on; xlabel('$k$'); legend('show');
+% FMT = '.fig'; save_figure
+
+%% CFL condition
+fig = figure; FIGNAME = ['CFL',sprintf('_%.2d',JOBNUM)];
+ semilogy(Ts,dz./ExB,'-','DisplayName','$|\nabla \phi|\Delta y$');
+ hold on;
+ plot(Ts,dt*ones(1,numel(Ts)),'--k','DisplayName','$\Delta t$');
+ grid on; xlabel('$t$'); ylabel('$\Delta t$'); legend('show');
+FMT = '.fig'; save_figure
+
+%% Space-Time diagrams at z = 0
+plt = @(x) real(x);
+fig = figure; FIGNAME = ['r_space_time_diag',sprintf('_%.2d',JOBNUM)];
+ [TY,TX] = meshgrid(Ts,z);
+subplot(221)% density
+ pclr = pcolor(TX,TY,(plt(ne_ST_r))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$r\,(z=0)$'); ylabel('$t$'); title('$n_e^{00}$');
+subplot(222)% density
+ pclr = pcolor(TX,TY,(plt(ni_ST_r))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$r\,(z=0)$'); ylabel('$t$'); title('$n_i^{00}$');
+subplot(223)% density
+ pclr = pcolor(TX,TY,(plt(phi_ST_r))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$r\,(z=0)$'); ylabel('$t$'); title('$\phi$');
+FMT = '.fig'; save_figure
+
+%% Space-Time diagrams at r = 0
+plt = @(x) real(x);
+fig = figure; FIGNAME = ['z_space_time_diag',sprintf('_%.2d',JOBNUM)];
+ [TY,TX] = meshgrid(Ts,r);
+subplot(221)% density
+ pclr = pcolor(TX,TY,(plt(ne_ST_z))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$z\,(r=0)$'); ylabel('$t$'); title('$n_e^{00}$');
+subplot(222)% density
+ pclr = pcolor(TX,TY,(plt(ni_ST_z))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$z\,(r=0)$'); ylabel('$t$'); title('$n_i^{00}$');
+subplot(223)% density
+ pclr = pcolor(TX,TY,(plt(phi_ST_z))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$z\,(r=0)$'); ylabel('$t$'); title('$\phi$');
+FMT = '.fig'; save_figure
+
+%% Space-Time diagrams for max_kz(Real)
+plt = @(x) (x);
+fig = figure; FIGNAME = ['kr_space_time_diag',sprintf('_%.2d',JOBNUM)];
+ [TY,TX] = meshgrid(Ts,kr);
+subplot(221)% density
+ pclr = pcolor(TX,TY,(plt(Ne_ST_kr))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$kr$'); ylabel('$t$'); title('$\max_{k_z}(\textrm{Re}(N_e^{00}))$');
+subplot(222)% density
+ pclr = pcolor(TX,TY,(plt(Ni_ST_kr))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$kr$'); ylabel('$t$'); title('$\max_{k_z}(\textrm{Re}(N_i^{00}))$');
+subplot(223)% density
+ pclr = pcolor(TX,TY,(plt(PH_ST_kr))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$kr$'); ylabel('$t$'); title('$\max_{k_z}(\textrm{Re}(\tilde\phi$))');
+FMT = '.fig'; save_figure
+
+%% Space-Time diagrams at max_kr(Real)
+plt = @(x) (x);
+fig = figure; FIGNAME = ['kz_space_time_diag',sprintf('_%.2d',JOBNUM)];
+ [TY,TX] = meshgrid(Ts,kz);
+subplot(221)% density
+ pclr = pcolor(TX,TY,(plt(Ne_ST_kz))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$kz$'); ylabel('$t$'); title('$\max_{k_r}(\textrm{Re}(N_e^{00}))$');
+subplot(222)% density
+ pclr = pcolor(TX,TY,(plt(Ni_ST_kz))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$kz$'); ylabel('$t$'); title('$\max_{k_r}(\textrm{Re}(N_i^{00}))$');
+subplot(223)% density
+ pclr = pcolor(TX,TY,(plt(PH_ST_kz))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$kz$'); ylabel('$t$'); title('$\max_{k_r}(\textrm{Re}(\tilde\phi))$');
+FMT = '.fig'; save_figure
+
+
+if 0
+%% Show frame
+it = min(250,numel(Ts));
+fig = figure; FIGNAME = ['frame',sprintf('_%.2d',JOBNUM)];
+ subplot(221); plt = @(x) fftshift((abs(x)));
+ pclr = pcolor(fftshift(KR),fftshift(KZ),plt(PH(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts(it))); legend('$|\hat\phi|$');
+ subplot(222); plt = @(x) fftshift(abs(x));
+ pclr = pcolor(fftshift(KR),fftshift(KZ),plt(Ni(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts(it))); legend('$|\hat n_i^{00}|$');
+ subplot(223); plt = @(x) fftshift(abs(x));
+ pclr = pcolor(fftshift(KR),fftshift(KZ),plt(Ne(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts(it))); legend('$|\hat n_e^{00}|$');
+if strcmp(OUTPUTS.write_non_lin,'.true.')
+ subplot(224); plt = @(x) fftshift((abs(x)));
+ pclr = pcolor(fftshift(KR),fftshift(KZ),plt(SNi(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$k_r$'); ylabel('$k_z$');legend('$\hat S_i^{00}$');
+end
+FMT = '.fig'; save_figure
+end
+%%
+DELAY = 0.07; skip_ = 10;
+FRAMES = 1:skip_:numel(Ts);
+if 0
+%% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Density electron
+GIFNAME = ['ne',sprintf('_%.2d',JOBNUM)];
+FIELD = real(ne); X = XX; Y = YY; T = Ts;
+FIELDNAME = '$n_e^{00}$'; XNAME = '$r$'; YNAME = '$z$';
+create_gif
+%% Density ion
+GIFNAME = ['ni',sprintf('_%.2d',JOBNUM)];
+FIELD = real(ni); X = XX; Y = YY; T = Ts;
+FIELDNAME = '$n_i^{00}$'; XNAME = '$r$'; YNAME = '$z$';
+create_gif
+%% Phi
+GIFNAME = ['phi',sprintf('_%.2d',JOBNUM)];
+FIELD = real(phi); X = XX; Y = YY; T = Ts;
+FIELDNAME = '$\phi$'; XNAME = '$r$'; YNAME = '$z$';
+create_gif
+%% Density electron frequency
+GIFNAME = ['Ni',sprintf('_%.2d',JOBNUM)];
+FIELD = ifftshift(ifftshift(abs(Ni),2),1); X = fftshift(KR); Y = fftshift(KZ); T = Ts;
+FIELDNAME = '$N_i^{00}$'; XNAME = '$k_r$'; YNAME = '$k_z$';
+create_gif
+end
\ No newline at end of file
diff --git a/wk/setup.m b/wk/setup.m
index 80bdb5c7..3a902cfc 100644
--- a/wk/setup.m
+++ b/wk/setup.m
@@ -5,69 +5,44 @@ default_plots_options
%% Set Up parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 1e-2; % Collision frequency
+NU = 1e-1; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-ETAB = 0.1; % Magnetic gradient
+ETAB = 0.5; % Magnetic gradient
ETAN = 1.0; % Density gradient
ETAT = 0.0; % Temperature gradient
-MU = 1e-2; % Hyper diffusivity coefficient
+MU = 1e-1; % Hyper diffusivity coefficient
LAMBDAD = 0.0;
%% GRID PARAMETERS
N = 64; % Frequency gridpoints (Nkr = N/2)
-L = 10; % Size of the squared frequency domain
-PMAXE = 01; % Highest electron Hermite polynomial degree
-JMAXE = 00; % Highest '' Laguerre ''
-PMAXI = 01; % Highest ion Hermite polynomial degree
-JMAXI = 00; % Highest '' Laguerre ''
+L = 150; % Size of the squared frequency domain
+KREQ0 = 0; % put kr = 0
+PMAXE = 00; % Highest electron Hermite polynomial degree
+JMAXE = 01; % Highest '' Laguerre ''
+PMAXI = 00; % Highest ion Hermite polynomial degree
+JMAXI = 01; % Highest '' Laguerre ''
KPAR = 0.0; % Parellel wave vector component
%% TIME PARAMETERS
-TMAX = 10.0; % Maximal time unit
+TMAX = 30.0; % Maximal time unit
DT = 1e-2; % Time step
SPS = 10; % Sampling per time unit
RESTART = 0; % To restart from last checkpoint
JOB2LOAD= 0;
%% OPTIONS
-SIMID = 'test_kzpos'; % Name of the simulation
-NON_LIN = 0; % activate non-linearity (is cancelled if KREQ0 = 1)
+SIMID = 'Turbulences'; % Name of the simulation
+NON_LIN = 0 *(1-KREQ0); % activate non-linearity (is cancelled if KREQ0 = 1)
CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
DK = 0; % Drift kinetic model (put every kernel to 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ________________________________________________________________________
-% Naming and creating input file
-params = ['P_',num2str(PMAXE),'_J_',num2str(JMAXE),...
- '_nB_',num2str(ETAB),'_nN_',num2str(ETAN),'_mu_%0.0e_'];
-params = sprintf(params,MU);
-if ~NON_LIN; params = ['lin_',params]; end;
-if DK; params = ['DK_',params]; end;
-resolution = ['_',num2str(N/2),'x',num2str(N),'_'];
-gridname = ['L_',num2str(L),'_'];
-BASIC.SIMID = [SIMID,resolution,gridname,params];
-BASIC.nrun = 1e8;
-BASIC.dt = DT;
-BASIC.tmax = TMAX; %time normalized to 1/omega_pe
-if RESTART; BASIC.RESTART = '.true.'; else; BASIC.RESTART = '.false.';end;
-OUTPUTS.nsave_0d = 0;
-OUTPUTS.nsave_1d = 0;
-OUTPUTS.nsave_2d = floor(1.0/SPS/DT);
-OUTPUTS.nsave_5d = floor(1.0/SPS/DT);
-OUTPUTS.write_Ni00 = '.true.';
-OUTPUTS.write_moments = '.true.';
-OUTPUTS.write_phi = '.true.';
-OUTPUTS.write_non_lin = '.true.';
-if NON_LIN == 0; OUTPUTS.write_non_lin = '.false.'; end;
-OUTPUTS.write_doubleprecision = '.true.';
-OUTPUTS.resfile0 = ['''',BASIC.SIMID,''''];
-OUTPUTS.rstfile0 = ['''','../checkpoint/cp_',BASIC.SIMID,''''];
-OUTPUTS.job2load = JOB2LOAD;
% 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.Nr = N; % r grid resolution
-GRID.Lr = L; % r length
+GRID.Nr = N * (1-KREQ0) + KREQ0; % r grid resolution
+GRID.Lr = L * (1-KREQ0); % r length
GRID.Nz = N; % z ''
GRID.Lz = L; % z ''
GRID.kpar = KPAR;
@@ -109,7 +84,34 @@ 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'''];
-
+% Naming and creating input file
+params = ['P_',num2str(PMAXE),'_J_',num2str(JMAXE),...
+ '_nB_',num2str(ETAB),'_nN_',num2str(ETAN),'_nu_%0.0e_',...
+ '_mu_%0.0e_'];
+params = sprintf(params,[NU,MU]);
+if ~NON_LIN; params = ['lin_',params]; end;
+if DK; params = ['DK_',params]; end;
+resolution = ['_',num2str(GRID.Nr),'x',num2str(GRID.Nz/2),'_'];
+gridname = ['L_',num2str(L),'_'];
+BASIC.SIMID = [SIMID,resolution,gridname,params];
+BASIC.nrun = 1e8;
+BASIC.dt = DT;
+BASIC.tmax = TMAX; %time normalized to 1/omega_pe
+% Outputs parameters
+if RESTART; BASIC.RESTART = '.true.'; else; BASIC.RESTART = '.false.';end;
+OUTPUTS.nsave_0d = 0;
+OUTPUTS.nsave_1d = 0;
+OUTPUTS.nsave_2d = floor(1.0/SPS/DT);
+OUTPUTS.nsave_5d = floor(1.0/SPS/DT);
+OUTPUTS.write_Ni00 = '.true.';
+OUTPUTS.write_moments = '.true.';
+OUTPUTS.write_phi = '.true.';
+OUTPUTS.write_non_lin = '.true.';
+if NON_LIN == 0; OUTPUTS.write_non_lin = '.false.'; end;
+OUTPUTS.write_doubleprecision = '.true.';
+OUTPUTS.resfile0 = ['''',BASIC.SIMID,''''];
+OUTPUTS.rstfile0 = ['''','../checkpoint/cp_',BASIC.SIMID,''''];
+OUTPUTS.job2load = JOB2LOAD;
%% Compile and write input file
INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC);
nproc = 1;
--
GitLab