diff --git a/matlab/create_gif.m b/matlab/create_gif.m
index f83f2f2ff66f5094f68c2f767cf727825b7d1924..179a71e591aaea02bcd8b0c61eed05dfb84cf383 100644
--- a/matlab/create_gif.m
+++ b/matlab/create_gif.m
@@ -1,45 +1,34 @@
-function [ ] = create_gif(x, y, t, h, BASIC, GRID, MODEL, delay, GIFNAME, saturate)
title1 = GIFNAME;
-FIGDIR = ['../results/', BASIC.SIMID,'/'];
+FIGDIR = ['../results/',BASIC.SIMID,'/'];
if ~exist(FIGDIR, 'dir')
mkdir(FIGDIR)
end
-if MODEL.CO == -1; CONAME = 'FC';
-elseif MODEL.CO == -2; CONAME = 'DC';
-elseif MODEL.CO == 0; CONAME = 'LB'; end;
-
-GIFNAME = [GIFNAME,'_Pe_',num2str(GRID.pmaxe),'_Je_',num2str(GRID.jmaxe),...
- '_Pi_',num2str(GRID.pmaxi),'_Ji_',num2str(GRID.jmaxi),...
- '_etan_',num2str(MODEL.eta_n),'_etaB_',num2str(MODEL.eta_B),'_etaT_',...
- num2str(MODEL.eta_T),'_nu_',num2str(MODEL.nu),'_',CONAME];
GIFNAME = [FIGDIR, GIFNAME,'.gif'];
% Set colormap boundaries
-hmax = max(max(max(h)));
-hmin = min(min(min(h)));
+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;
- pcolor(x,y,h(:,:,1)); % to set up
+% Setup figure frame
+fig = figure;
+ pcolor(X,Y,FIELD(:,:,1)); % to set up
colormap jet
colorbar
- if not(flag) && saturate
- caxis([hmin hmax])
- end
axis tight manual % this ensures that getframe() returns a consistent size
+ shading interp;
hold on
-
+
n = 0;
- nbytes = fprintf(2,'frame %d/%d',n,numel(h(1,1,:)));
- for n = 1:numel(h(1,1,:)) % time loop
- pclr = pcolor(x,y,h(:,:,n)); % frame plot
- set(pclr, 'edgecolor','none')
- title([title1,', $t \approx$', sprintf('%.3d',ceil(t(n)))])
+ nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,:)));
+ for n = 1:numel(FIELD(1,1,:)) % time loop
+ pclr = pcolor(X,Y,FIELD(:,:,n)); shading interp; % frame plot
+ set(pclr, 'edgecolor','none');
+ title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))]);
drawnow
% Capture the plot as an image
frame = getframe(fig);
@@ -49,17 +38,16 @@ else
if n == 1
imwrite(imind,cm,GIFNAME,'gif', 'Loopcount',inf);
else
- imwrite(imind,cm,GIFNAME,'gif','WriteMode','append', 'DelayTime',delay);
+ 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(h(1,1,:)));
+ nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,:)));
end
-
+ disp(' ')
disp(['Gif saved @ : ',GIFNAME])
end
-end
diff --git a/matlab/default_plots_options.m b/matlab/default_plots_options.m
index 9f4bfb3eace935347857a9a8e87dc086d0f4645e..ed53cc039f90ba647c6f2c417f7ad72e716e7578 100644
--- a/matlab/default_plots_options.m
+++ b/matlab/default_plots_options.m
@@ -6,15 +6,4 @@ set(0,'defaultAxesFontSize',16)
set(0, 'DefaultLineLineWidth', 2.0);
line_colors = [[0, 0.4470, 0.7410];[0.8500, 0.3250, 0.0980];[0.9290, 0.6940, 0.1250];...
[0.4940, 0.1840, 0.5560];[0.4660, 0.6740, 0.1880];[0.3010, 0.7450, 0.9330];[0.6350, 0.0780, 0.1840]];
-line_styles = {'-','-.','--',':'};
-
-%Title management
-if MODEL.CO == -1; CONAME = 'FC';
-elseif MODEL.CO == -2; CONAME = 'DC';
-elseif MODEL.CO == 0; CONAME = 'LB'; end;
-TITLE = [];
-TITLE = [TITLE,'$\eta_n=',num2str(MODEL.eta_n),'$, '];
-TITLE = [TITLE,'$\eta_B=',num2str(MODEL.eta_B),'$, '];
-TITLE = [TITLE,'$\eta_T=',num2str(MODEL.eta_T),'$, '];
-TITLE = [TITLE, '$\nu=',num2str(MODEL.nu),'$, '];
-TITLE = [TITLE, '$(P,J)=(',num2str(GRID.pmaxe),',',num2str(GRID.jmaxe),')$'];
\ No newline at end of file
+line_styles = {'-','-.','--',':'};
\ No newline at end of file
diff --git a/matlab/load_2D_data.m b/matlab/load_2D_data.m
index ebdc15385b2a1fd43689d6f6fd3374e5a45aa6d7..ddb42e1c7e35694fd3209bb4afcdcee065d3f16c 100644
--- a/matlab/load_2D_data.m
+++ b/matlab/load_2D_data.m
@@ -3,9 +3,14 @@ function [ data, kr, kz, time ] = load_2D_data( filename, variablename )
time = h5read(filename,'/data/var2d/time');
kr = h5read(filename,['/data/var2d/',variablename,'/coordkr']);
kz = h5read(filename,['/data/var2d/',variablename,'/coordkz']);
+
+ dt = h5readatt(filename,'/data/input','dt');
+ nsave = h5readatt(filename,'/data/input','nsave_5d');
+ cstart= time(1)/dt/nsave;
+
data = zeros(numel(kr),numel(kz),numel(time));
for it = 1:numel(time)
- tmp = h5read(filename,['/data/var2d/',variablename,'/', num2str(it,'%06d')]);
+ tmp = h5read(filename,['/data/var2d/',variablename,'/', num2str(cstart+it,'%06d')]);
data(:,:,it) = tmp.real + 1i * tmp.imaginary;
end
diff --git a/matlab/load_5D_data.m b/matlab/load_5D_data.m
index 6f51cc1c345abfbb15a239ae285b19ff01382071..34cb449f5b50e0fa88489eb1c47febbe14aef195 100644
--- a/matlab/load_5D_data.m
+++ b/matlab/load_5D_data.m
@@ -5,9 +5,15 @@ function [ data, p, j, kr, kz, time ] = load_5D_data( filename, variablename )
j = h5read(filename,['/data/var5d/',variablename,'/coordj']);
kr = h5read(filename,['/data/var5d/',variablename,'/coordkr']);
kz = h5read(filename,['/data/var5d/',variablename,'/coordkz']);
+
+ dt = h5readatt(filename,'/data/input','dt');
+ nsave = h5readatt(filename,'/data/input','nsave_5d');
+ cstart= time(1)/dt/nsave;
+
data = zeros(numel(p),numel(j),numel(kr),numel(kz),numel(time));
+
for it = 1:numel(time)
- tmp = h5read(filename,['/data/var5d/', variablename,'/', num2str(it,'%06d')]);
+ tmp = h5read(filename,['/data/var5d/', variablename,'/', num2str(cstart+it,'%06d')]);
data(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
end
end
\ No newline at end of file
diff --git a/matlab/save_figure.m b/matlab/save_figure.m
index 4d569a630d44ec754f5830b2d46ebd26ff87b3da..9fae159ba479532926bd5f2df9f5bb50e17f1191 100644
--- a/matlab/save_figure.m
+++ b/matlab/save_figure.m
@@ -5,10 +5,6 @@ if ~exist(FIGDIR, 'dir')
mkdir(FIGDIR)
end
-FIGNAME = [FIGNAME,'_Pe_',num2str(GRID.pmaxe),'_Je_',num2str(GRID.jmaxe),...
- '_Pi_',num2str(GRID.pmaxi),'_Ji_',num2str(GRID.jmaxi),...
- '_etan_',num2str(MODEL.eta_n),'_etaB_',num2str(MODEL.eta_B),'_etaT_',...
- num2str(MODEL.eta_T),'_nu_',num2str(MODEL.nu),'_',CONAME];
-FIGNAME = [FIGDIR, FIGNAME,'.fig'];
-savefig(fig,FIGNAME);
+FIGNAME = [FIGDIR, FIGNAME,FMT];
+saveas(fig,FIGNAME);
disp(['Figure saved @ : ',FIGNAME])
\ No newline at end of file
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index b15369864ea9661be9b7161f40d49e42b42adbc7..2e818e7f0e97f58f2e247cb4c992d33c3990a268 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -2,24 +2,25 @@ function [INPUT] = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASI
% Write the input script "fort.90" with desired parameters
INPUT = '../wk/fort.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,[' RESTART = ', num2str(BASIC.RESTART),'\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,[' nkr = ', num2str(GRID.nkr),'\n']);
-fprintf(fid,[' krmin = ', num2str(GRID.krmin),'\n']);
-fprintf(fid,[' krmax = ', num2str(GRID.krmax),'\n']);
-fprintf(fid,[' nkz = ', num2str(GRID.nkz),'\n']);
-fprintf(fid,[' kzmin = ', num2str(GRID.kzmin),'\n']);
-fprintf(fid,[' kzmax = ', num2str(GRID.kzmax),'\n']);
-fprintf(fid,[' Pad = ', num2str(GRID.Pad),'\n']);
+fprintf(fid,[' Nr = ', num2str(GRID.Nr),'\n']);
+fprintf(fid,[' Lr = ', num2str(GRID.Lr),'\n']);
+fprintf(fid,[' Nz = ', num2str(GRID.Nz),'\n']);
+fprintf(fid,[' Lz = ', num2str(GRID.Lz),'\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']);
@@ -32,11 +33,14 @@ fprintf(fid,[' write_non_lin = ', OUTPUTS.write_non_lin,'\n']);
fprintf(fid,[' write_doubleprecision = ', OUTPUTS.write_doubleprecision,'\n']);
fprintf(fid,[' resfile0 = ', OUTPUTS.resfile0,'\n']);
fprintf(fid,[' rstfile0 = ', OUTPUTS.rstfile0,'\n']);
+fprintf(fid,[' job2load = ', num2str(OUTPUTS.job2load),'\n']);
fprintf(fid,'/\n');
+
fprintf(fid,'&MODEL_PAR\n');
fprintf(fid,' ! Collisionality\n');
fprintf(fid,[' CO = ', num2str(MODEL.CO),'\n']);
fprintf(fid,[' NON_LIN = ', MODEL.NON_LIN,'\n']);
+fprintf(fid,[' mu = ', num2str(MODEL.mu),'\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']);
@@ -49,8 +53,8 @@ fprintf(fid,[' eta_T = ', num2str(MODEL.eta_T),'\n']);
fprintf(fid,[' eta_B = ', num2str(MODEL.eta_B),'\n']);
fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),'\n']);
fprintf(fid,'/\n');
+
fprintf(fid,'&INITIAL_CON\n');
-fprintf(fid,[' RESTART =', INITIAL.RESTART,'\n']);
fprintf(fid,[' only_Na00 =', INITIAL.only_Na00,'\n']);
fprintf(fid,[' initback_moments =', num2str(INITIAL.initback_moments),'\n']);
fprintf(fid,[' initnoise_moments =', num2str(INITIAL.initnoise_moments),'\n']);
@@ -59,8 +63,10 @@ fprintf(fid,[' selfmat_file =', INITIAL.selfmat_file,'\n']);
fprintf(fid,[' eimat_file =', INITIAL.eimat_file,'\n']);
fprintf(fid,[' iemat_file =', INITIAL.iemat_file,'\n']);
fprintf(fid,'/\n');
+
fprintf(fid,'&TIME_INTEGRATION_PAR\n');
fprintf(fid,[' numerical_scheme=', TIME_INTEGRATION.numerical_scheme,'\n']);
fprintf(fid,'/');
+
fclose(fid);
end
diff --git a/src/auxval.F90 b/src/auxval.F90
index 85ef9e9b833e10aaa2592ed654fcbd0b153f50c0..f0c451c91a5bf5d4feafe7ed7f8d3d930699a33d 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -11,6 +11,8 @@ subroutine auxval
INTEGER :: irows,irowe, irow, icol
WRITE(*,*) '=== Set auxiliary values ==='
+ call set_rgrid
+ call set_zgrid
call set_krgrid
call set_kzgrid
call set_pj
diff --git a/src/compute_Sapj.F90 b/src/compute_Sapj.F90
index f6c2fd2f2a8e2c4e3c3823ae0a39a1a50468c0a0..8a38278967623d862f779540085c32f71398b7ff 100644
--- a/src/compute_Sapj.F90
+++ b/src/compute_Sapj.F90
@@ -1,6 +1,6 @@
-SUBROUTINE compute_Sapj(Sepj, Sipj)
+SUBROUTINE compute_Sapj
- USE array, ONLY : dnjs
+ USE array, ONLY : dnjs, Sepj, Sipj
USE basic
USE fourier
USE fields!, ONLY : phi, moments_e, moments_i
@@ -11,8 +11,8 @@ SUBROUTINE compute_Sapj(Sepj, Sipj)
IMPLICIT NONE
COMPLEX(dp), DIMENSION(Nkr,Nkz) :: F_, G_, CONV
- COMPLEX(dp), DIMENSION(ips_e:ipe_e, ijs_e:ije_e, Nkr, Nkz) :: Sepj
- COMPLEX(dp), DIMENSION(ips_i:ipe_i, ijs_i:ije_i, Nkr, Nkz) :: Sipj
+ ! COMPLEX(dp), DIMENSION(ips_e:ipe_e, ijs_e:ije_e, Nkr, Nkz) :: Sepj
+ ! COMPLEX(dp), DIMENSION(ips_i:ipe_i, ijs_i:ije_i, Nkr, Nkz) :: Sipj
INTEGER :: in, is
REAL(dp):: kr, kz, kernel, be_2, bi_2, factn
REAL(dp):: sigmae2_taue_o2, sigmai2_taui_o2
@@ -34,7 +34,7 @@ SUBROUTINE compute_Sapj(Sepj, Sipj)
kernel = be_2**(in-1)/factn * EXP(-be_2)
! First convolution term
- F_(ikr,ikz) = (kz - kr) * phi(ikr,ikz) * kernel
+ F_(ikr,ikz) = (kz - kr) * phi(ikr,ikz) !* kernel
! Second convolution term
G_(ikr,ikz) = 0._dp ! initialization of the sum
DO is = 1, MIN( in+ij-1, jmaxe+1 ) ! sum truncation on number of moments
@@ -78,7 +78,7 @@ SUBROUTINE compute_Sapj(Sepj, Sipj)
bi_2 = sigmai2_taui_o2 * (kr**2 + kz**2)
kernel = bi_2**(in-1)/factn * EXP(-bi_2)
- F_(ikr,ikz) = (kz - kr) * phi(ikr,ikz) * kernel
+ F_(ikr,ikz) = (kz - kr) * phi(ikr,ikz) !* kernel
G_(ikr,ikz) = 0._dp ! initialization of the sum
DO is = 1, MIN( in+ij-1, jmaxi+1 )
diff --git a/src/control.F90 b/src/control.F90
index 9218566ace85a08c12a50bc6b44fbe991107a22d..59abdcf2b463d48c30a37545635f32c7d0128f93 100644
--- a/src/control.F90
+++ b/src/control.F90
@@ -27,7 +27,7 @@ SUBROUTINE control
CALL readinputs
WRITE(*,'(a/)') '...input parameters read'
- ! 1.4 Set auxiliary values (allocate arrays, set laplace operator, ...)
+ ! 1.4 Set auxiliary values (allocate arrays, set grid, ...)
WRITE(*,*) 'Calculate auxval...'
CALL auxval
WRITE(*,'(a/)') '...auxval calculated'
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index 47e9dacf784463330e3a5bc86bbc64c67e33a1a0..0c8a7dcb2cca0f8ad9bcadc2b38922f1fa49ac82 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -173,6 +173,11 @@ SUBROUTINE diagnose(kstep)
!
IF (kstep .GE. 0) THEN
+ ! Terminal info
+ IF (MOD(cstep, INT(1.0/dt)) == 0) THEN
+ WRITE(*,"(F4.0,A,F4.0)") time,"/",tmax
+ ENDIF
+
! 2.1 0d history arrays
IF (nsave_0d .NE. 0) THEN
IF ( MOD(cstep, nsave_0d) == 0 ) THEN
@@ -187,9 +192,6 @@ SUBROUTINE diagnose(kstep)
IF (nsave_2d .NE. 0) THEN
IF (MOD(cstep, nsave_2d) == 0) THEN
CALL diagnose_2d
- IF (MOD(cstep, nsave_2d*10) == 0) THEN
- WRITE(*,"(F4.0,A,F4.0)") time,"/",tmax
- ENDIF
END IF
END IF
@@ -360,7 +362,7 @@ SUBROUTINE checkpoint_save
USE time_integration
IMPLICIT NONE
- WRITE(rstfile,'(a,a3)') TRIM(rstfile0),'.h5'
+ WRITE(rstfile,'(a,a1,i2.2,a3)') TRIM(rstfile0),'_',jobnum,'.h5'
CALL creatf(rstfile, fidrst, real_prec='d')
CALL creatg(fidrst, '/Basic', 'Basic data')
@@ -377,6 +379,6 @@ SUBROUTINE checkpoint_save
CALL putarr(fidrst, '/Basic/moments_i', moments_i(ips_i:ipe_i,ijs_i:ije_i,&
ikrs:ikre,ikzs:ikze,1),ionode=0)
CALL closef(fidrst)
- WRITE(*,'(3x,a)') "Checkpoint file "//TRIM(rstfile)//" saved!"
+ WRITE(*,'(3x,a)') "Checkpoint file "//TRIM(rstfile)//" saved"
END SUBROUTINE checkpoint_save
diff --git a/src/diagnostics_par_mod.F90 b/src/diagnostics_par_mod.F90
index a9005ea86262e395fdfc0b29823d84d19c40f7cb..074aa71a30a6a6c2863d8df1a4e7deb936b1e5a2 100644
--- a/src/diagnostics_par_mod.F90
+++ b/src/diagnostics_par_mod.F90
@@ -12,15 +12,16 @@ MODULE diagnostics_par
LOGICAL, PUBLIC, PROTECTED :: write_doubleprecision=.FALSE.
INTEGER, PUBLIC, PROTECTED :: nsave_0d , nsave_1d , nsave_2d , nsave_5d
- INTEGER, PUBLIC, PROTECTED :: nsave_cp = 1e3
+ INTEGER, PUBLIC, PROTECTED :: nsave_cp = 1e4
! HDF5 file
- CHARACTER(len=64), PUBLIC :: resfile0 = "results" ! Head of main result file name
- CHARACTER(len=64), PUBLIC :: resfile ! Main result file
- INTEGER, PUBLIC :: fidres ! FID for resfile
- CHARACTER(len=64), PUBLIC :: rstfile0 = "restart" ! Head of restart file name
- CHARACTER(len=64), PUBLIC :: rstfile ! Full restart file
- INTEGER, PUBLIC :: fidrst ! FID for restart file
+ CHARACTER(len=128), PUBLIC :: resfile0 = "results" ! Head of main result file name
+ CHARACTER(len=128), PUBLIC :: resfile ! Main result file
+ INTEGER, PUBLIC :: job2load ! jobnum of the checkpoint to load
+ INTEGER, PUBLIC :: fidres ! FID for resfile
+ CHARACTER(len=128), PUBLIC :: rstfile0 = "restart" ! Head of restart file name
+ CHARACTER(len=128), PUBLIC :: rstfile ! Full restart file
+ INTEGER, PUBLIC :: fidrst ! FID for restart file
PUBLIC :: output_par_readinputs, output_par_outputinputs
@@ -36,7 +37,7 @@ CONTAINS
NAMELIST /OUTPUT_PAR/ nsave_0d , nsave_1d , nsave_2d , nsave_5d
NAMELIST /OUTPUT_PAR/ write_Ni00, write_moments, write_phi, write_non_lin, write_doubleprecision
- NAMELIST /OUTPUT_PAR/ resfile0, rstfile0
+ NAMELIST /OUTPUT_PAR/ resfile0, rstfile0, job2load
READ(lu_in,output_par)
diff --git a/src/inital.F90 b/src/inital.F90
index 69312a55cbf39fd49bbcf3bdd7e00c2b2ffc626e..36d14818e3d51129094525272bf9442ef7ea8463 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -7,20 +7,27 @@ SUBROUTINE inital
USE model, ONLY : CO, NON_LIN
USE prec_const
USE coeff
+ USE time_integration
USE array, ONLY : Sepj,Sipj
implicit none
+
+ CALL set_updatetlevel(1)
!!!!!! Set the moments arrays Nepj, Nipj !!!!!!
+ write(*,*) 'Init moments'
IF ( RESTART ) THEN
- CALL init_moments
+ CALL load_cp
ELSE
- CALL load_moments
+ CALL init_moments
ENDIF
!!!!!! Set phi !!!!!!
+ write(*,*) 'Init phi'
CALL poisson
!!!!!! Set Sepj, Sipj and dnjs coeff table !!!!!!
IF ( NON_LIN ) THEN;
- CALL compute_Sapj(Sepj,Sipj)
+ write(*,*) 'Init Sapj'
+ CALL compute_Sapj
+ WRITE(*,*) 'Building Dnjs table'
CALL build_dnjs_table
ENDIF
!!!!!! Load the full coulomb collision operator coefficients !!!!!!
@@ -38,20 +45,19 @@ SUBROUTINE init_moments
USE basic
USE grid
USE fields
- USE initial_par
- USE time_integration
USE prec_const
+ USE utility, ONLY: checkfield
+ USE initial_par
IMPLICIT NONE
- INTEGER, DIMENSION(12) :: iseedarr
REAL(dp) :: noise
+ REAL(dp) :: kr, kz, sigma, gain
+ INTEGER, DIMENSION(12) :: iseedarr
! Seed random number generator
iseedarr(:)=iseed
CALL RANDOM_SEED(PUT=iseedarr)
- CALL set_updatetlevel(1)
-
IF ( only_Na00 ) THEN ! Spike initialization on density only
DO ikr=ikrs,ikre
@@ -61,30 +67,42 @@ SUBROUTINE init_moments
END DO
END DO
- ELSE ! Broad noise initialization
-
- DO ikr=ikrs,ikre
- DO ikz=ikzs,ikze
- DO ip=ips_e,ipe_e
- DO ij=ijs_e,ije_e
- CALL RANDOM_NUMBER(noise)
- moments_e( ip,ij, ikr, ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
- moments_e( ip,ij, Nkz-ikz, ikr, :) = moments_e( ip,ij, ikr, ikz, :) ! Symmetry
- moments_e( ip,ij, ikz,Nkr-ikr, :) = moments_e( ip,ij, ikr, ikz, :) ! Symmetry
- moments_e( ip,ij, Nkz-ikz,Nkr-ikr, :) = moments_e( ip,ij, ikr, ikz, :) ! Symmetry
- END DO
- END DO
- DO ip=ips_i,ipe_i
- DO ij=ijs_i,ije_i
- CALL RANDOM_NUMBER(noise)
- moments_i( ip,ij, ikr, ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
- moments_i( ip,ij, Nkz-ikz, ikr, :) = moments_i( ip,ij, ikr, ikz, :) ! Symmetry
- moments_i( ip,ij, ikz,Nkr-ikr, :) = moments_i( ip,ij, ikr, ikz, :) ! Symmetry
- moments_i( ip,ij, Nkz-ikz,Nkr-ikr, :) = moments_i( ip,ij, ikr, ikz, :) ! Symmetry
- END DO
+ ELSE
+ sigma = 2._dp
+ gain = 0.5_dp
+ !**** Gaussian initialization (Hakim 2017)
+ moments_i = 0; moments_e = 0;
+ DO ikr=ikrs,ikre
+ kr = krarray(ikr)
+ DO ikz=ikzs,ikze
+ kz = kzarray(ikz)
+ moments_i( 1,1, ikr,ikz, :) = gain*sigma/SQRT2 * EXP(-(kr**2+kz**2)*sigma**2/4._dp)
+ moments_e( 1,1, ikr,ikz, :) = gain*sigma/SQRT2 * EXP(-(kr**2+kz**2)*sigma**2/4._dp)
END DO
END DO
- END DO
+ ! Broad noise initialization
+ ! DO ikr=ikrs,ikre
+ ! DO ikz=ikzs,ikze
+ ! DO ip=ips_e,ipe_e
+ ! DO ij=ijs_e,ije_e
+ ! CALL RANDOM_NUMBER(noise)
+ ! moments_e( ip,ij, ikr, ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
+ ! moments_e( ip,ij, Nkz-ikz, ikr, :) = moments_e( ip,ij, ikr, ikz, :) ! Symmetry
+ ! moments_e( ip,ij, ikz,Nkr-ikr, :) = moments_e( ip,ij, ikr, ikz, :) ! Symmetry
+ ! moments_e( ip,ij, Nkz-ikz,Nkr-ikr, :) = moments_e( ip,ij, ikr, ikz, :) ! Symmetry
+ ! END DO
+ ! END DO
+ ! DO ip=ips_i,ipe_i
+ ! DO ij=ijs_i,ije_i
+ ! CALL RANDOM_NUMBER(noise)
+ ! moments_i( ip,ij, ikr, ikz, :) = initback_moments + initnoise_moments*(noise-0.5_dp)
+ ! moments_i( ip,ij, Nkz-ikz, ikr, :) = moments_i( ip,ij, ikr, ikz, :) ! Symmetry
+ ! moments_i( ip,ij, ikz,Nkr-ikr, :) = moments_i( ip,ij, ikr, ikz, :) ! Symmetry
+ ! moments_i( ip,ij, Nkz-ikz,Nkr-ikr, :) = moments_i( ip,ij, ikr, ikz, :) ! Symmetry
+ ! END DO
+ ! END DO
+ ! END DO
+ ! END DO
ENDIF
END SUBROUTINE init_moments
!******************************************************************************!
@@ -92,21 +110,36 @@ END SUBROUTINE init_moments
!******************************************************************************!
!!!!!!! Load moments from a previous save
!******************************************************************************!
-SUBROUTINE load_moments
-! USE basic
-! USE initial_par
-! USE fields, ONLY : moments_e, moments_i
-! USE futils, ONLY : openf, getarr, closef
-! IMPLICIT NONE
-!
-! INTEGER :: fid
-!
-! CALL openf(backup_file,fid, 'r', 'D');
-! CALL getarr(fid, '/moments_e', moments_e) ! Nepj
-! CALL getarr(fid, '/moments_i', moments_i) ! Nipj
-! CALL getarr(fid, '/time', time) ! time
-! CALL closef(fid)
-END SUBROUTINE
+SUBROUTINE load_cp
+ USE basic
+ USE futils, ONLY: openf, closef, getarr, getatt
+ USE grid
+ USE fields
+ USE diagnostics_par
+ USE time_integration
+ IMPLICIT NONE
+
+ WRITE(rstfile,'(a,a1,i2.2,a3)') TRIM(rstfile0),'_',job2load,'.h5'
+
+ WRITE(*,'(3x,a)') "Resume from previous run"
+
+ CALL openf(rstfile, fidrst)
+ CALL getatt(fidrst, '/Basic', 'cstep', cstep)
+ CALL getatt(fidrst, '/Basic', 'time', time)
+ CALL getatt(fidrst, '/Basic', 'jobnum', jobnum)
+ jobnum = jobnum+1
+ CALL getatt(fidrst, '/Basic', 'iframe2d',iframe2d)
+ CALL getatt(fidrst, '/Basic', 'iframe5d',iframe5d)
+ iframe2d = iframe2d-1; iframe5d = iframe5d-1
+
+ ! Read state of system from restart file
+ CALL getarr(fidrst, '/Basic/moments_e', moments_e(ips_e:ipe_e,ijs_e:ije_e,ikrs:ikre,ikzs:ikze,1),ionode=0)
+ CALL getarr(fidrst, '/Basic/moments_i', moments_i(ips_i:ipe_i,ijs_i:ije_i,ikrs:ikre,ikzs:ikze,1),ionode=0)
+ CALL closef(fidrst)
+
+ WRITE(*,'(3x,a)') "Reading from restart file "//TRIM(rstfile)//" completed!"
+
+END SUBROUTINE load_cp
!******************************************************************************!
!******************************************************************************!
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 230afb3db09c53c284f144fffd1971f48129b9b5..9d6f3ff916578f8831d95a34f19b4941fdfde6f9 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -7,6 +7,7 @@ MODULE model
INTEGER, PUBLIC, PROTECTED :: CO = 0 ! Collision Operator
LOGICAL, PUBLIC, PROTECTED :: NON_LIN = .true. ! To turn on non linear bracket term
+ REAL(dp), PUBLIC, PROTECTED :: mu = 0._dp ! Hyperdiffusivity coefficient (for num. stability)
REAL(dp), PUBLIC, PROTECTED :: nu = 1._dp ! Collision frequency
REAL(dp), PUBLIC, PROTECTED :: tau_e = 1._dp ! Temperature
REAL(dp), PUBLIC, PROTECTED :: tau_i = 1._dp !
@@ -30,7 +31,7 @@ CONTAINS
USE prec_const
IMPLICIT NONE
- NAMELIST /MODEL_PAR/ CO, NON_LIN, nu, tau_e, tau_i, sigma_e, sigma_i, &
+ NAMELIST /MODEL_PAR/ CO, NON_LIN, mu, nu, tau_e, tau_i, sigma_e, sigma_i, &
q_e, q_i, eta_n, eta_T, eta_B, lambdaD
READ(lu_in,model_par)
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index 0b85811c9caa211046a9e4aacc3ac040c2b726ea..4ad31556154b5b51710f205dbc09dc37665b9d19 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -7,6 +7,7 @@ SUBROUTINE moments_eq_rhs
USE grid
USE model
USE prec_const
+ USE utility, ONLY : is_nan
IMPLICIT NONE
INTEGER :: ip2, ij2 ! loops indices
@@ -20,8 +21,8 @@ SUBROUTINE moments_eq_rhs
REAL(dp) :: xCapj, xCa20, xCa01, xCa10 ! Coll. factors depending on the pj loop
COMPLEX(dp) :: TNapj, TNapp2j, TNapm2j, TNapjp1, TNapjm1, Tphi
COMPLEX(dp) :: TColl, TColl20, TColl01, TColl10 ! terms of the rhs
+ COMPLEX(dp) :: test_nan
REAL(dp) :: nu_e, nu_i, nu_ee, nu_ie ! Species collisional frequency
-!write(*,*) '----------------------------------------'
!Precompute species dependant factors
taue_qe_etaB = tau_e/q_e * eta_B ! factor of the magnetic moment coupling
@@ -181,7 +182,11 @@ SUBROUTINE moments_eq_rhs
IF ( (ip .EQ. 1) .OR. (ip .EQ. 3) ) THEN ! kronecker p0 or p2
kernelj = b_e2**(ij-1) * exp(-b_e2)/factj
kerneljp1 = kernelj * b_e2 /(ij_dp + 1._dp)
- kerneljm1 = kernelj * ij_dp / b_e2
+ IF ( b_e2 .NE. 0 ) THEN
+ kerneljm1 = kernelj * ij_dp / b_e2
+ ELSE
+ kerneljm1 = 0.5_dp
+ ENDIF
Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz)
ELSE
Tphi = 0._dp
@@ -192,10 +197,11 @@ SUBROUTINE moments_eq_rhs
-imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 - Tphi)&
+ TColl
- ! Adding non linearity
+ ! Adding non linearity and Hyperdiffusivity
IF ( NON_LIN ) THEN
moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) = &
- moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) - Sepj(ip,ij,ikr,ikz)
+ moments_rhs_e(ip,ij,ikr,ikz,updatetlevel) - Sepj(ip,ij,ikr,ikz) &
+ - mu*kperp2**2 * moments_rhs_e(ip,ij,ikr,ikz,updatetlevel)
ENDIF
END DO kzloope
@@ -352,7 +358,11 @@ SUBROUTINE moments_eq_rhs
IF ( (ip .EQ. 1) .OR. (ip .EQ. 3) ) THEN ! kronecker p0 or p2
kernelj = b_i2**(ij-1) * exp(-b_i2)/factj
kerneljp1 = kernelj * b_i2 /(ij_dp + 1._dp)
- kerneljm1 = kernelj * ij_dp / b_i2
+ IF ( b_i2 .NE. 0 ) THEN
+ kerneljm1 = kernelj * ij_dp / b_i2
+ ELSE
+ kerneljm1 = 0.5_dp
+ ENDIF
Tphi = (xphij*kernelj + xphijp1*kerneljp1 + xphijm1*kerneljm1) * phi(ikr,ikz)
ELSE
Tphi = 0._dp
@@ -363,10 +373,11 @@ SUBROUTINE moments_eq_rhs
-imagu * kz * (TNapj + TNapp2j + TNapm2j + TNapjp1 + TNapjm1 - Tphi)&
+ TColl
- ! Adding non linearity
+ ! Adding non linearity and Hyperdiffusivity
IF ( NON_LIN ) THEN
moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) = &
- moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) - Sipj(ip,ij,ikr,ikz)
+ moments_rhs_i(ip,ij,ikr,ikz,updatetlevel) - Sipj(ip,ij,ikr,ikz)&
+ - mu*kperp2**2 * moments_rhs_i(ip,ij,ikr,ikz,updatetlevel)
ENDIF
END DO kzloopi
END DO krloopi
diff --git a/src/poisson.F90 b/src/poisson.F90
index fd7f9030f1e7fe37db1b6c53f3fd6e4282e40842..d91fe1ffbd50a381517a5ab4bbd69bd84a19104d 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -42,8 +42,6 @@ SUBROUTINE poisson
sum_kernel_mom_e = Kne*moments_e(1, 1, ikr, ikz, updatetlevel)
sum_kernel2_e = Kne**2
-!write(*,*) 'K0e = ', Kne
-
! loop over n only if the max polynomial degree is 1 or more
if (jmaxe .GT. 0) then
DO ine=2,jmaxe+1 ! ine = n+1
@@ -53,7 +51,6 @@ SUBROUTINE poisson
sum_kernel_mom_e = sum_kernel_mom_e + Kne * moments_e(1, ine, ikr, ikz, updatetlevel)
sum_kernel2_e = sum_kernel2_e + Kne**2 ! ... sum recursively ...
-!write(*,*) 'K',ine-1,'e = ', Kne
END DO
endif
@@ -65,8 +62,6 @@ SUBROUTINE poisson
sum_kernel_mom_i = Kni*moments_i(1, 1, ikr, ikz, updatetlevel)
sum_kernel2_i = Kni**2
-!write(*,*) 'K0i = ', Kni
-
! loop over n only if the max polynomial degree is 1 or more
if (jmaxi .GT. 0) then
DO ini=2,jmaxi+1
@@ -76,7 +71,6 @@ SUBROUTINE poisson
sum_kernel_mom_i = sum_kernel_mom_i + Kni * moments_i(1, ini, ikr, ikz, updatetlevel)
sum_kernel2_i = sum_kernel2_i + Kni**2 ! ... sum recursively ...
-!write(*,*) 'K',ini-1,'i = ', Kni
END DO
endif
!!! Assembling the poisson equation
@@ -85,12 +79,7 @@ SUBROUTINE poisson
+ qi2_taui * (1._dp - sum_kernel2_i)
gammaD_phi = q_e * sum_kernel_mom_e + q_i * sum_kernel_mom_i
- IF ( (alphaD .EQ. 0._dp) .AND. (gammaD .EQ. 0._dp) ) THEN
- write(*,*) "Warning : 0/0 occuring"
- phi(ikr,ikz) = 0._dp
- ELSE
- phi(ikr, ikz) = gammaD_phi/gammaD
- ENDIF
+ phi(ikr, ikz) = gammaD_phi/gammaD
END DO
END DO
diff --git a/wk/analysis.m b/wk/analysis.m
new file mode 100644
index 0000000000000000000000000000000000000000..229dbd1963c75a2511372d0cee086616720abb2c
--- /dev/null
+++ b/wk/analysis.m
@@ -0,0 +1,201 @@
+%% load results
+JOBNUM = 00;
+filename = [BASIC.SIMID,'_','%.2d.h5'];
+filename = sprintf(filename,JOBNUM); disp(['Analysing ',filename])
+[Nipj, p_, j_, kr, kz, Ts] = 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 = 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
+phi_ST_r = zeros(Nr,Ns); % Space-Time diagram of ES potential
+ne_ST_r = zeros(Nr,Ns); % Space-Time diagram of density
+ni_ST_r = zeros(Nr,Ns); % Space-Time diagram of density
+phi_ST_z = zeros(Nz,Ns); % Space-Time diagram of ES potential
+ne_ST_z = zeros(Nz,Ns); % Space-Time diagram of density
+ni_ST_z = zeros(Nz,Ns); % Space-Time diagram of density
+phi_00 = zeros(1,Ns); % Time evolution of ES potential at origin
+ne_00 = zeros(1,Ns); % Time evolution of density at origin
+ni_00 = zeros(1,Ns); % Time evolution of density at origin
+[~,ir0] = min(abs(r)); [~,iz0] = min(abs(z));
+Ne_11 = zeros(1,Ns); % Time evolution of F density at 1,1
+Ni_11 = zeros(1,Ns); % Time evolution of F density at 1,1
+[~,ikr1] = min(abs(kr-round(1/dkr)*dkr)); [~,ikz1] = min(abs(kz-round(1/dkz)*dkz));
+Sni_norm= zeros(1,Ns); % Time evolution of the amp of density nonlin term
+Sne_norm= zeros(1,Ns); % Time evolution of the amp of vorti. 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;
+for it = 1:numel(PH(1,1,:))
+ NE_ = Ne(:,:,it); NN_ = Ni(:,:,it); PH_ = PH(:,:,it);
+ phi_ST_r(:,it) = phi(:,iz0,it); phi_ST_z(:,it) = phi(ir0,:,it);
+ ne_ST_r (:,it)= ne(:,iz0,it); ne_ST_z (:,it)= ne(ir0,:,it);
+ ni_ST_r (:,it)= ni(:,iz0,it); ni_ST_z (:,it)= ni(ir0,:,it);
+ phi_00(it) = phi(ir0,iz0,it);
+ ne_00(it) = ne(ir0,iz0,it);
+ ni_00(it) = ni(ir0,iz0,it);
+ Ne_11(it) = Ne(ikr1,ikz1,it);
+ Ni_11(it) = Ni(ikr1,ikz1,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_11),'-','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
+
+%% phi
+fig = figure; FIGNAME = ['phi_ST',sprintf('_%.2d',JOBNUM)];
+ [TY,TX] = meshgrid(Ts,z);
+ pclr = pcolor(TX,TY,(plt(phi_ST_r))); set(pclr, 'edgecolor','none'); colorbar;
+ xlabel('$x\,(y=0)$'); ylabel('$t$'); title('$\phi$');
+FMT = '.fig'; save_figure
+
+if 0
+%% Show frame
+it = min(1,numel(Ts));
+fig = figure; FIGNAME = ['frame',sprintf('_%.2d',SID)];
+ subplot(221); plt = @(x) fftshift((real(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(real(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(real(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(log10(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.1; skip_ = 2;
+if 0
+%% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Density electron
+GIFNAME = ['ne',sprintf('_%.2d',JOBNUM)]; FIELDNAME = '$n_e^{00}$';
+FIELD = real(ne(:,:,1:skip_:end)); X = XX; Y = YY; T = Ts(1:skip_:end);
+create_gif
+%% Density ion
+GIFNAME = ['ni',sprintf('_%.2d',JOBNUM)]; FIELDNAME = '$n_i^{00}$';
+FIELD = real(ni(:,:,1:skip_:end)); X = XX; Y = YY; T = Ts(1:skip_:end);
+create_gif
+%% Phi
+GIFNAME = ['phi',sprintf('_%.2d',JOBNUM)]; FIELDNAME = '$\phi$';
+FIELD = real(phi(:,:,1:skip_:end)); X = XX; Y = YY; T = Ts(1:skip_:end);
+create_gif
+end
\ No newline at end of file
diff --git a/wk/fort.90 b/wk/fort.90
index b02202518128b93b13ec1009856e3020e899771f..c11d29567d34becbe94fc861c90d030b6a6efb84 100644
--- a/wk/fort.90
+++ b/wk/fort.90
@@ -1,16 +1,59 @@
&BASIC
- nrun=100000
- dt=0.01
- tmax=200
+ nrun = 100000000
+ dt = 0.001
+ tmax = 200
+ RESTART = .true.
/
&GRID
- pmaxe =5
- jmaxe = 3
- pmaxi = 5
- jmaxi = 3
- nkr = 1
- krmin = 0
- krmax = 0
- nkz = 1
- kzmin = 0.1
- kzmax = 0.1
+ pmaxe =1
+ jmaxe = 0
+ pmaxi = 1
+ jmaxi = 0
+ Nr = 64
+ Lr = 20
+ Nz = 64
+ Lz = 20
+/
+&OUTPUT_PAR
+ nsave_0d = 0
+ nsave_1d = 0
+ nsave_2d = 1000
+ nsave_5d = 1000
+ write_Ni00 = .true.
+ write_moments = .true.
+ write_phi = .true.
+ write_non_lin = .true.
+ write_doubleprecision = .true.
+ resfile0 = 'SwoK_32x64_L_20_P_1_J_0_nB_0.5_nN_1_mu_1e-04_'
+ rstfile0 = '../checkpoint/cp_SwoK_32x64_L_20_P_1_J_0_nB_0.5_nN_1_mu_1e-04_'
+ job2load = 0
+/
+&MODEL_PAR
+ ! Collisionality
+ CO = 0
+ NON_LIN = .true.
+ mu = 0.0001
+ nu = 0
+ tau_e = 1
+ tau_i = 1
+ sigma_e = 0.023338
+ sigma_i = 1
+ q_e = -1
+ q_i = 1
+ eta_n = 1
+ eta_T = 0
+ eta_B = 0.5
+ lambdaD = 0
+/
+&INITIAL_CON
+ only_Na00 =.false.
+ initback_moments =0.0001
+ initnoise_moments =5e-05
+ iseed =42
+ selfmat_file ='../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_1_Jmaxe_0_Pmaxi_1_Jmaxi_0_pamaxx_10.h5'
+ eimat_file ='../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_1_Jmaxe_0_Pmaxi_1_Jmaxi_0_pamaxx_10_tau_1.0000_mu_0.0233.h5'
+ iemat_file ='../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_1_Jmaxe_0_Pmaxi_1_Jmaxi_0_pamaxx_10_tau_1.0000_mu_0.0233.h5'
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme='RK4'
+/
\ No newline at end of file
diff --git a/wk/run.m b/wk/run.m
new file mode 100644
index 0000000000000000000000000000000000000000..47b7d4332aeb8434ccdecbadcc0cf562d4d2ca07
--- /dev/null
+++ b/wk/run.m
@@ -0,0 +1,5 @@
+%% Run HeLaZ
+EXEC = ' ../bin/helaz ';
+RUN = ['mpirun -np ' num2str(nproc)];
+CMD = [RUN, EXEC, INPUT];
+system(CMD);
\ No newline at end of file
diff --git a/wk/setup.m b/wk/setup.m
new file mode 100644
index 0000000000000000000000000000000000000000..c1c45262f8e083a5cede87e93f9f2af168c30b41
--- /dev/null
+++ b/wk/setup.m
@@ -0,0 +1,113 @@
+clear all; close all;
+addpath(genpath('../matlab')) % ... add
+default_plots_options
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Set Up parameters
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% PHYSICAL PARAMETERS
+NU = 0.0; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+ETAB = 0.5; % Magnetic gradient
+ETAN = 1.0; % Density gradient
+ETAT = 0.0; % Temperature gradient
+MU = 1e-4; % Hyper diffusivity coefficient
+%% GRID PARAMETERS
+N = 64; % Frequency gridpoints
+L = 20; % 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 ''
+%% TIME PARAMETERS
+TMAX = 200.0; % Maximal time unit
+DT = 1e-3; % Time step
+SPS = 1; % Sampling per time unit
+RESTART = 1; % To restart from last checkpoint
+JOB2LOAD= 0;
+%% OPTIONS
+SIMID = 'SwoK'; % Name of the simulation
+NON_LIN = 1; % activate non-linearity
+CO = 0;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%% ________________________________________________________________________
+% 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;
+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.Nz = N; % z ''
+GRID.Lz = L; % z ''
+% Model parameters
+MODEL.CO = CO; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+if NON_LIN; MODEL.NON_LIN = '.true.'; else; MODEL.NON_LIN = '.false.';end;
+MODEL.mu = MU;
+MODEL.nu = NU; % hyper diffusive coefficient nu for HW
+% temperature ratio T_a/T_e
+MODEL.tau_e = TAU;
+MODEL.tau_i = TAU;
+% mass ratio sqrt(m_a/m_i)
+MODEL.sigma_e = 0.0233380;
+MODEL.sigma_i = 1.0;
+% charge q_a/e
+MODEL.q_e =-1.0;
+MODEL.q_i = 1.0;
+% gradients L_perp/L_x
+MODEL.eta_n = ETAN; % source term kappa for HW
+MODEL.eta_T = ETAT; % Temperature
+MODEL.eta_B = ETAB; % Magnetic
+MODEL.lambdaD = 0.0;
+% Time integration and intialization parameters
+TIME_INTEGRATION.numerical_scheme = '''RK4''';
+INITIAL.only_Na00 = '.false.';
+INITIAL.initback_moments = 1.0e-4;
+INITIAL.initnoise_moments = 5.0e-5;
+INITIAL.iseed = 42;
+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'''];
+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'''];
+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'''];
+
+%% Compile and write input file
+INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC);
+nproc = 1;
+MAKE = 'cd ..; make; cd wk';
+system(MAKE);
+%%
+disp(['Set up ', BASIC.SIMID]);
\ No newline at end of file