diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index 6e1afdce7c2ef6762ad096e1cfc94a730dde0083..17b9a01892ae518e96bd9e84c8addadd131b0d49 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -277,7 +277,7 @@ else
% grids
DATA.Pe = Pe; DATA.Pi = Pi;
DATA.Je = Je; DATA.Ji = Ji;
- DATA.kx = kx; DATA.ky = ky; DATA.z = z;
+ DATA.kx = kx; DATA.ky = ky; DATA.z = z; DATA.Npol = -z(1)/pi;
DATA.x = x; DATA.y = y;
DATA.ikx0 = ikx0; DATA.iky0 = iky0;
DATA.Nx = Nx; DATA.Ny = Ny; DATA.Nz = Nz; DATA.Nkx = Nkx; DATA.Nky = Nky;
diff --git a/matlab/compute/compute_3D_zpinch_growth_rate.m b/matlab/compute/compute_3D_zpinch_growth_rate.m
index 521e5a1f0b738935818fcc5551f45dc77a11fe18..c89597b199ae13fc9730969a80577c60710c94ca 100644
--- a/matlab/compute/compute_3D_zpinch_growth_rate.m
+++ b/matlab/compute/compute_3D_zpinch_growth_rate.m
@@ -4,7 +4,7 @@ FIGURE.FIGNAME = ['growth_rate_kx=0_ky=0_planes',DATA.PARAMS];
t = DATA.Ts3D;
[~,its] = min(abs(t-TRANGE(1)));
[~,ite] = min(abs(t-TRANGE(end)));
-nkx = DATA.Nkx; nky = DATA.Nky; nkz = DATA.Nz;
+nkx = DATA.Nkx; nky = DATA.Nky; nkz = DATA.Nz/2;
% Remove AA part
if DATA.Nx > 1
ikxnz = abs(DATA.PHI(1,:,1,1)) > 0;
@@ -15,84 +15,77 @@ ikynz = abs(DATA.PHI(:,1,1,1)) > 0;
ikynz(1) = 1; ikxnz(1) = 1; %put k=0 in the analysis
-phi = fft(DATA.PHI(ikynz,ikxnz,:,:),[],3);
+Y = fft(DATA.PHI(ikynz,ikxnz,:,:),[],3);
+phi = Y(:,:,2:2:end,:);
omega = zeros(nky,nkx,nkz);
-for iz = 1:nkz
- for iy = 1:nky
+for ikz = 1:nkz
+ for iky = 1:nky
for ix = 1:nkx
- omega(iy,ix,iz) = LinearFit_s(t(its:ite),squeeze(abs(phi(iy,ix,iz,its:ite))));
+% omega(iy,ix,iz) = LinearFit_s(t(its:ite),squeeze(abs(phi(iy,ix,iz,its:ite))));
+ to_measure = squeeze(log(phi(iky,ix,ikz,its:ite)));
+ tmp = polyfit(t(its:ite),to_measure(:),1);
+ if ~(isnan(tmp(1)) || isinf(tmp(1)))
+ omega(iky,ix,ikz) = tmp(1);
+ end
end
end
end
%% plot
-kx = DATA.kx; ky = DATA.ky; kz = [(0:nkz/2), (-nkz/2+1):-1];
-poskx = kx>=0;
-posky = ky>=0;
-poskz = kz>=0;
+kx = DATA.kx; ky = DATA.ky; kz = [(0:nkz/2), (-nkz/2+1):-1]/DATA.Npol;
kxeq0 = kx==0;
kyeq0 = ky==0;
kzeq0 = kz==0;
-omega = omega(posky,poskx,poskz);
+% kx = fftshift(kx,1);
+% ky = fftshift(ky,1);
+% kz = fftshift(kz,1);
FIGURE.fig = figure;
nplots = OPTIONS.kxky + OPTIONS.kzkx + OPTIONS.kzky;
iplot = 1;
if OPTIONS.kxky
-[Y_XY,X_XY] = meshgrid(ky(posky),kx(poskx));
+[Y_XY,X_XY] = meshgrid(ky,kx);
subplot(1,nplots,iplot)
- if ~OPTIONS.keq0
- toplot = squeeze(max(real(omega(:,:,:)),[],3));
- pclr= pcolor(X_XY,Y_XY,toplot'); set(pclr,'EdgeColor','none');
- xlabel('$k_x$'); ylabel('$k_y$');
- title('$\max(\gamma)_{kz}$');
- else
- toplot = squeeze(real(omega(:,:,kzeq0)));
- pclr= pcolor(X_XY,Y_XY,toplot'); set(pclr,'EdgeColor','none');
- xlabel('$k_x$'); ylabel('$k_y$');
- title('$\gamma(k_z=0)$');
- end
+ toplot = squeeze(real(omega(:,:,kzeq0)));
+ toplot = fftshift(toplot,2);
+ X_XY = fftshift(X_XY,1);
+ Y_XY = fftshift(Y_XY,1);
+ pclr= pcolor(X_XY,Y_XY,toplot'); set(pclr,'EdgeColor','none');
+ xlabel('$k_x$'); ylabel('$k_y$');
+ title('$\gamma(k_z=0)$');
if OPTIONS.INTERP; shading interp; end
caxis(max(max(abs(toplot))).*[-1,1]);
iplot = iplot + 1;
end
if OPTIONS.kzky
-[Y_ZY,Z_ZY] = meshgrid(ky(posky),kz(poskz));
+[Y_ZY,Z_ZY] = meshgrid(ky,kz);
subplot(1,nplots,iplot)
- if ~OPTIONS.keq0
- toplot = squeeze(max(real(omega(:,:,:)),[],1));
- pclr= pcolor(Z_ZY,Y_ZY,toplot'); set(pclr,'EdgeColor','none');
- xlabel('$k_x$'); ylabel('$k_y$');
- title('$\max(\gamma)_{kx}$');
- else
- toplot = squeeze(real(omega(:,kxeq0,:)));
- pclr= pcolor(Z_ZY,Y_ZY,toplot'); set(pclr,'EdgeColor','none');
- xlabel('$k_z$'); ylabel('$k_y$');
- title('$\gamma(k_x=0)$');
- end
+ toplot = squeeze(real(omega(:,kxeq0,:)));
+ toplot = fftshift(toplot,2);
+ Z_ZY = fftshift(Z_ZY,1);
+ Y_ZY = fftshift(Y_ZY,1);
+ pclr= pcolor(Z_ZY,Y_ZY,toplot'); set(pclr,'EdgeColor','none');
+ xlabel('$k_z$'); ylabel('$k_y$');
+ title('$\gamma(k_x=0)$');
if OPTIONS.INTERP; shading interp; end
caxis(max(max(abs(toplot))).*[-1,1]);
iplot = iplot + 1;
end
if OPTIONS.kzkx
-[X_ZX,Z_ZX] = meshgrid(kx(poskx),kz(poskz));
+[X_ZX,Z_ZX] = meshgrid(kx,kz);
subplot(1,nplots,iplot)
- if ~OPTIONS.keq0
- toplot = squeeze(max(real(omega(:,:,:)),[],2));
- pclr= pcolor(Z_ZX,X_ZX,toplot'); set(pclr,'EdgeColor','none');
- xlabel('$k_z$'); ylabel('$k_x$');
- title('$\max(\gamma)_{ky}$');
- else
- toplot = squeeze(real(omega(kyeq0,:,:)));
- pclr= pcolor(Z_ZX,X_ZX,toplot'); set(pclr,'EdgeColor','none');
- xlabel('$k_z$'); ylabel('$k_x$');
- title('$\gamma(k_y=0)$');
- end
+ toplot = squeeze(real(omega(kyeq0,:,:)));
+ toplot = fftshift(toplot,2);
+ Z_ZX = fftshift(Z_ZX,1);
+ X_ZX = fftshift(X_ZX,1);
+ pclr= pcolor(Z_ZX,X_ZX,toplot'); set(pclr,'EdgeColor','none');
+ xlabel('$k_z$'); ylabel('$k_x$');
+ title('$\gamma(k_y=0)$');
end
if OPTIONS.INTERP; shading interp; end
caxis(max(max(abs(toplot))).*[-1,1]);
diff --git a/matlab/load/load_params.m b/matlab/load/load_params.m
index d3d141fbbe3e4a14d028505805cd34420b86347a..6d657839c7a70c7c3b031ba752eb987945113119 100644
--- a/matlab/load/load_params.m
+++ b/matlab/load/load_params.m
@@ -24,6 +24,12 @@ catch
DATA.K_T = h5readatt(filename,'/data/input','k_Ti');
end
end
+DATA.sigma_e = h5readatt(filename,'/data/input','sigma_e');
+DATA.sigma_i = h5readatt(filename,'/data/input','sigma_i');
+DATA.tau_e = h5readatt(filename,'/data/input','tau_e');
+DATA.tau_i = h5readatt(filename,'/data/input','tau_i');
+DATA.q_e = h5readatt(filename,'/data/input','q_e');
+DATA.q_i = h5readatt(filename,'/data/input','q_i');
DATA.Q0 = h5readatt(filename,'/data/input','q0');
DATA.SHEAR = h5readatt(filename,'/data/input','shear');
DATA.EPS = h5readatt(filename,'/data/input','eps');
diff --git a/matlab/plot/show_moments_spectrum.m b/matlab/plot/show_moments_spectrum.m
index 8c1468aa89901d04f2b3916375b2a98df4a416ee..2fc32fa49077cc0639693fe5f4f6d588d02371ac 100644
--- a/matlab/plot/show_moments_spectrum.m
+++ b/matlab/plot/show_moments_spectrum.m
@@ -22,7 +22,7 @@ if ~OPTIONS.ST
[~,it0] = min(abs(t0-DATA.Ts5D));
[~,it1] = min(abs(t1-DATA.Ts5D));
Nipj = mean(Nipj(:,:,it0:it1),3);
- if DATA.K_E
+ if DATA.KIN_E
Nepj = mean(Nepj(:,:,it0:it1),3);
end
if numel(OPTIONS.TIME) == 1
@@ -34,13 +34,13 @@ if ~OPTIONS.ST
ymini = min(Nipj); ymaxi = max(Nipj);
- if DATA.K_E
+ if DATA.KIN_E
Nepj = squeeze(Nepj);
ymine = min(Nepj); ymaxe = max(Nepj);
ymax = max([ymaxi ymaxe]);
ymin = min([ymini ymine]);
end
- if DATA.K_E
+ if DATA.KIN_E
subplot(121)
end
if ~OPTIONS.P2J
@@ -60,7 +60,7 @@ if ~OPTIONS.ST
legend('show');
title([TITLE,' He-La ion spectrum']);
- if DATA.K_E
+ if DATA.KIN_E
subplot(122)
if ~OPTIONS.P2J
for ij = 1:numel(Je)
diff --git a/matlab/process_field_v2.m b/matlab/process_field_v2.m
index 4525a46a4398524ffcfa416acd4ab4848f34f984..48122855a277532ef8f9aac340e42cd1f56dfa65 100644
--- a/matlab/process_field_v2.m
+++ b/matlab/process_field_v2.m
@@ -65,7 +65,29 @@ 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);
+% --
switch OPTIONS.COMP
case 'avg'
compr = @(x) mean(x,COMPDIM);
@@ -130,6 +152,7 @@ end
switch OPTIONS.NAME
case '\phi' %ES pot
NAME = 'phi';
+% FLD_ = DATA.PHI(:,:,:,FRAMES); % data to plot
FLD_ = DATA.PHI(:,:,:,FRAMES); % data to plot
OPE_ = 1; % Operation on data
case '\psi' %EM pot
@@ -168,13 +191,17 @@ switch OPTIONS.NAME
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);
+ FLD_ = DATA.DENS_I(:,:,:,FRAMES) - adiab_i.* DATA.PHI(:,:,:,FRAMES);
OPE_ = 1;
case 'n_e' % electron density
NAME = 'ne';
- FLD_ = DATA.DENS_E(:,:,:,FRAMES);
+ FLD_ = DATA.DENS_E(:,:,:,FRAMES) - adiab_e.* DATA.PHI(:,:,:,FRAMES);
OPE_ = 1;
case 'k^2n_e' % electron vorticity
NAME = 'k2ne';
@@ -183,7 +210,8 @@ switch OPTIONS.NAME
case 'n_i-n_e' % polarisation
NAME = 'pol';
OPE_ = 1;
- FLD_ = DATA.DENS_I(:,:,:,FRAMES)+DATA.DENS_E(:,:,:,FRAMES);
+ 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);
diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m
index 423c3a09bf6d72c4f9aa12515f84bb8121d2cdb8..88578a5804636c1d1f167a79cd3460bb5c4b6f24 100644
--- a/wk/analysis_gyacomo.m
+++ b/wk/analysis_gyacomo.m
@@ -14,7 +14,7 @@ MISCDIR = ['/misc/gyacomo_outputs/',resdir,'/']; %For long term storage
system(['mkdir -p ',MISCDIR]);
system(['mkdir -p ',LOCALDIR]);
CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
-system(CMD);
+% system(CMD);
% Load outputs from jobnummin up to jobnummax
data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
data.localdir = LOCALDIR;
@@ -68,14 +68,15 @@ options.POLARPLOT = 0;
% options.NAME = 'v_x';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
-options.NAME = 'n_i-n_e';
+options.NAME = 'n_e';
+% options.NAME = 'n_i-n_e';
options.PLAN = 'xy';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
% options.TIME = data.Ts5D(end-30:end);
% options.TIME = data.Ts3D;
-options.TIME = [500:1:9000];
+options.TIME = [00:1:9000];
data.EPS = 0.1;
data.a = data.EPS * 2000;
options.RESOLUTION = 256;
@@ -92,16 +93,17 @@ options.NORMALIZE = 0;
options.NAME = '\phi';
% options.NAME = '\psi';
% options.NAME = '\omega_z';
-% options.NAME = 'n_i';
+% options.NAME = 'n_e';
% options.NAME = 'n_i-n_e';
% options.NAME = '\phi^{NZ}';
% options.NAME = 'N_i^{00}';
+% options.NAME = 'N_i^{00}-N_e^{00}';
% options.NAME = 's_{Ex}';
% options.NAME = 'Q_x';
% options.NAME = 'k^2n_e';
options.PLAN = 'xy';
options.COMP = 'avg';
-options.TIME = [800];
+options.TIME = [500];
options.RESOLUTION = 256;
data.a = data.EPS * 2e3;
diff --git a/wk/header_2DZP_results.m b/wk/header_2DZP_results.m
index 0b25b98b25a49228b6dff11aae87b8f292189d06..24b1929c77514f5835c0a7d0c6dd9b4260c910af 100644
--- a/wk/header_2DZP_results.m
+++ b/wk/header_2DZP_results.m
@@ -194,7 +194,7 @@ resdir ='Zpinch_rerun/Kn_1.6_256x128x21x11';
% resdir ='Zpinch_rerun/nu_0.1_FCGK_200x48x5x3_kN_scan';
% resdir ='Zpinch_rerun/nu_0.1_SGGK_200x48x5x3_kN_scan';
% resdir = 'Zpinch_rerun/kN_1.7_FCGK_200x32x5x3_nu_scan';
-% % resdir = 'Zpinch_rerun/kN_1.7_SGGK_200x32x7x4_nu_scan';
+% resdir = 'Zpinch_rerun/kN_1.7_SGGK_200x32x7x4_nu_scan';
% resdir = 'Zpinch_rerun/kN_2.2_SGGK_200x32x5x3_nu_scan';
% resdir = 'Zpinch_rerun/kN_1.7_SGGK_256x64x5x3_nu_scan';
@@ -230,7 +230,7 @@ resdir ='Zpinch_rerun/Kn_1.6_256x128x21x11';
% resdir = 'Zpinch_rerun/DGGK_kN_scan_200x64x5x3_nu_0.01';
% resdir = 'Zpinch_rerun/even_DGGK_kN_scan_200x64x5x3_nu_0.01';
% resdir = 'Zpinch_rerun/LRGK_kN_scan_200x64x5x3_nu_0.01';
-% resdir = 'Zpinch_rerun/LDGK_kN_scan_200x64x5x3_nu_0.01';
+% resdir = 'Zpinch_rerun/LDGK_kN_scan_202x64x5x3_nu_0.01';
% resdir = 'Zpinch_rerun/DGGK_kN_1.6_200x64x5x3_nu_0.01';
% resdir = 'Zpinch_rerun/DGGK_kN_1.7_200x64x5x3_nu_0.01';
@@ -240,6 +240,6 @@ resdir ='Zpinch_rerun/Kn_1.6_256x128x21x11';
% resdir = 'Zpinch_rerun/COLLESS_kN_1.7_200x64x5x3';
% resdir = 'Zpinch_rerun/COLLESS_kN_1.7_200x64x9x5';
-JOBNUMMIN = 00; JOBNUMMAX = 10;
+JOBNUMMIN = 01; JOBNUMMAX = 03;
resdir = ['results/',resdir];
run analysis_gyacomo
\ No newline at end of file
diff --git a/wk/lin_3Dzpinch.m b/wk/lin_3Dzpinch.m
new file mode 100644
index 0000000000000000000000000000000000000000..66ae45c1423946599fea2bff67d1bdc08154dbdc
--- /dev/null
+++ b/wk/lin_3Dzpinch.m
@@ -0,0 +1,201 @@
+%% QUICK RUN SCRIPT for linear entropy mode (ETPY) in a Zpinch
+% This script create a directory in /results and run a simulation directly
+% from matlab framework. It is meant to run only small problems in linear
+% for benchmark and debugging purpose since it makes matlab "busy"
+%
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % ... add
+addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
+addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
+addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
+SIMID = 'dbg'; % Name of the simulation
+RUN = 0; % To run or just to load
+default_plots_options
+% EXECNAME = 'gyacomo_debug';
+% EXECNAME = 'gyacomo_alpha';
+EXECNAME = 'gyacomo';
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Set Up parameters
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% PHYSICAL PARAMETERS
+NU = 0.1; % Collision frequency
+TAU = 1e-2; % e/i temperature ratio
+K_Ne = 0; % ele Density '''
+K_Te = 0; % ele Temperature '''
+K_Ni = 0; % ion Density gradient drive
+K_Ti = 70; % ion Temperature '''
+SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 0.0; % electron plasma beta
+%% GRID PARAMETERS
+P = 4;
+J = P/2;
+PMAXE = P; % Hermite basis size of electrons
+JMAXE = J; % Laguerre "
+PMAXI = P; % " ions
+JMAXI = J; % "
+NX = 20; % real space x-gridpoints
+NY = 20; % '' y-gridpoints
+LX = 2*pi/0.4; % Size of the squared frequency domain
+LY = 2*pi/0.4; % Size of the squared frequency domain
+NZ = 48; % number of perpendicular planes (parallel grid)
+NPOL = 1;
+SG = 0; % Staggered z grids option
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+%% GEOMETRY
+%% GEOMETRY
+GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
+EPS = 0.18; % inverse aspect ratio
+Q0 = 1.4; % safety factor
+SHEAR = 0.0; % magnetic shear
+KAPPA = 1.0; % elongation
+DELTA = 0.0; % triangularity
+ZETA = 0.0; % squareness
+PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
+% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
+SHIFT_Y = 0.0;
+%% TIME PARMETERS
+TMAX = 50; % Maximal time unit
+DT = 5e-2; % Time step
+SPS0D = 1; % Sampling per time unit for 2D arrays
+SPS2D = -1; % Sampling per time unit for 2D arrays
+SPS3D = 1; % Sampling per time unit for 2D arrays
+SPS5D = 1/2; % Sampling per time unit for 5D arrays
+SPSCP = 0; % Sampling per time unit for checkpoints
+JOB2LOAD= -1;
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+% Collision operator
+% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+CO = 'DG';
+GKCO = 0; % gyrokinetic operator
+ABCO = 1; % interspecies collisions
+INIT_ZF = 0; ZF_AMP = 0.0;
+CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
+NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
+KERN = 0; % Kernel model (0 : GK)
+INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom
+NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
+%% 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
+HD_CO = 0.0; % Hyper diffusivity cutoff ratio
+MU = 0.0; % Hyperdiffusivity coefficient
+INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
+MU_X = MU; %
+MU_Y = MU; %
+N_HD = 4;
+MU_Z = 1.0; %
+MU_P = 0.0; %
+MU_J = 0.0; %
+LAMBDAD = 0.0;
+NOISE0 = 1.0e-5; % Init noise amplitude
+BCKGD0 = 0.0; % Init background
+GRADB = 0.1;
+CURVB = 0.1;
+COLL_KCUT = 1000;
+%%-------------------------------------------------------------------------
+%% RUN
+setup
+% system(['rm fort*.90']);
+% Run linear simulation
+if RUN
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 3 2 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
+end
+
+%% Load results
+%%
+filename = [SIMID,'/',PARAMS,'/'];
+LOCALDIR = [gyacomodir,'results/',filename,'/'];
+FIGDIR = LOCALDIR;
+% Load outputs from jobnummin up to jobnummax
+JOBNUMMIN = 00; JOBNUMMAX = 01;
+data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
+data.FIGDIR = ['../results/',SIMID,'/',PARAMS,'/'];
+
+%% Short analysis
+if 0
+%% linear growth rate (adapted for 2D zpinch and fluxtube)
+options.TRANGE = [0.5 1]*data.Ts3D(end);
+options.NPLOTS = 3; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
+options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3
+lg = compute_fluxtube_growth_rate(data,options);
+[gmax, kmax] = max(lg.g_ky(:,end));
+[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
+msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg);
+msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg);
+end
+
+if 0
+%% Ballooning plot
+options.time_2_plot = [10 50];
+options.kymodes = [0.1 0.2 0.4];
+options.normalized = 1;
+% options.field = 'phi';
+fig = plot_ballooning(data,options);
+end
+
+if 0
+%% Hermite-Laguerre spectrum
+% options.TIME = 'avg';
+options.P2J = 0;
+options.ST = 0;
+options.PLOT_TYPE = 'space-time';
+% options.PLOT_TYPE = 'Tavg-1D';ls
+
+% options.PLOT_TYPE = 'Tavg-2D';
+options.NORMALIZED = 1;
+options.JOBNUM = 0;
+options.TIME = [0 50];
+options.specie = 'i';
+options.compz = 'avg';
+fig = show_moments_spectrum(data,options);
+% fig = show_napjz(data,options);
+% save_figure(data,fig)
+end
+
+if 1
+%% linear growth rate for 3D Zpinch (kz fourier transform)
+trange = [0.5 1]*data.Ts3D(end);
+options.INTERP = 0;
+options.kxky = 0;
+options.kzkx = 0;
+options.kzky = 1;
+[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
+end
+if 0
+%% Mode evolution
+options.NORMALIZED = 0;
+options.K2PLOT = 1;
+options.TIME = [0:1000];
+% options.TIME = [0.5:0.01:1].*data.Ts3D(end);
+options.NMA = 1;
+options.NMODES = 32;
+options.iz = 'avg';
+options.ik = 1;
+fig = mode_growth_meter(data,options);
+save_figure(data,fig,'.png')
+end
+
+
+if 0
+%% RH TEST
+ikx = 2; iky = 2; t0 = 0; t1 = data.Ts3D(end);
+[~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D));
+plt = @(x) squeeze(mean(real(x(iky,ikx,:,it0:it1)),3))./squeeze(mean(real(x(iky,ikx,:,it0)),3));
+figure
+plot(data.Ts3D(it0:it1), plt(data.PHI));
+xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$')
+title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',data.kx(ikx),data.ky(iky)))
+end
\ No newline at end of file
diff --git a/wk/lin_ETPY.m b/wk/lin_ETPY.m
index cdc2a5f583946c0b67b2ce17750125100d3135b7..8f25bcd5c1ac1015624c261853944945d58ce7c3 100644
--- a/wk/lin_ETPY.m
+++ b/wk/lin_ETPY.m
@@ -3,13 +3,17 @@
% from matlab framework. It is meant to run only small problems in linear
% for benchmark and debugging purpose since it makes matlab "busy"
%
-SIMID = 'lin_ETPY'; % Name of the simulation
-% SIMID = 'dbg'; % Name of the simulation
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % ... add
+addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
+addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
+addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
+SIMID = 'dbg'; % Name of the simulation
RUN = 1; % To run or just to load
-addpath(genpath('../matlab')) % ... add
default_plots_options
-PROGDIR = '/home/ahoffman/gyacomo/';
-EXECNAME = 'gyacomo_dbg';
+% EXECNAME = 'gyacomo_debug';
+EXECNAME = 'gyacomo_alpha';
% EXECNAME = 'gyacomo';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
@@ -26,35 +30,40 @@ SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 0.0; % electron plasma beta
%% GRID PARAMETERS
-P = 8;
+P = 4;
J = P/2;
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
NX = 2; % real space x-gridpoints
-NY = 100; % '' y-gridpoints
+NY = 40; % '' y-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
-LY = 200;%2*pi/0.05; % Size of the squared frequency domain
+LY = 2*pi/0.05; % Size of the squared frequency domain
NZ = 1; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+%% GEOMETRY
%% GEOMETRY
GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
+EPS = 0.18; % inverse aspect ratio
Q0 = 1.4; % safety factor
-SHEAR = 0.8; % magnetic shear
-KAPPA = 0.0; % elongation
+SHEAR = 0.0; % magnetic shear
+KAPPA = 1.0; % elongation
DELTA = 0.0; % triangularity
ZETA = 0.0; % squareness
-NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
-EPS = 0.18; % inverse aspect ratio
+PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
+% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
+SHIFT_Y = 0.0;
%% TIME PARMETERS
-TMAX = 500; % Maximal time unit
+TMAX = 50; % Maximal time unit
+% DT = 1e-2; % Time step
DT = 1e-2; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = -1; % Sampling per time unit for 2D arrays
SPS3D = 1; % Sampling per time unit for 2D arrays
-SPS5D = 1/5; % Sampling per time unit for 5D arrays
+SPS5D = 1/2; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
%% OPTIONS
@@ -62,14 +71,14 @@ LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
CO = 'DG';
-GKCO = 1; % gyrokinetic operator
+GKCO = 0; % gyrokinetic operator
ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
KERN = 0; % Kernel model (0 : GK)
-INIT_OPT= 'mom00'; % Start simulation with a noisy mom00/phi/allmom
-NUMERICAL_SCHEME = 'RK2'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
+INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom
+NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
%% OUTPUTS
W_DOUBLE = 1;
W_GAMMA = 1; W_HF = 1;
@@ -85,7 +94,7 @@ INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
MU_Y = MU; %
N_HD = 4;
-MU_Z = 2.0; %
+MU_Z = 0.2; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
@@ -100,23 +109,27 @@ setup
% system(['rm fort*.90']);
% Run linear simulation
if RUN
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',PROGDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',PROGDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
end
%% Load results
%%
filename = [SIMID,'/',PARAMS,'/'];
-LOCALDIR = [PROGDIR,'results/',filename,'/'];
+LOCALDIR = [gyacomodir,'results/',filename,'/'];
+FIGDIR = LOCALDIR;
% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 00;
+JOBNUMMIN = 00; JOBNUMMAX = 01;
data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
+data.FIGDIR = ['../results/',SIMID,'/',PARAMS,'/'];
%% Short analysis
-if 1
+if 0
%% linear growth rate (adapted for 2D zpinch and fluxtube)
options.TRANGE = [0.5 1]*data.Ts3D(end);
-options.NPLOTS = 2; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
+options.NPLOTS = 3; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3
lg = compute_fluxtube_growth_rate(data,options);
[gmax, kmax] = max(lg.g_ky(:,end));
@@ -127,8 +140,8 @@ end
if 0
%% Ballooning plot
-options.time_2_plot = [120];
-options.kymodes = [0.9];
+options.time_2_plot = [10 50];
+options.kymodes = [0.1 0.2 0.4];
options.normalized = 1;
% options.field = 'phi';
fig = plot_ballooning(data,options);
@@ -140,9 +153,10 @@ if 0
options.P2J = 1;
options.ST = 1;
options.PLOT_TYPE = 'space-time';
-% options.PLOT_TYPE = 'Tavg-1D';
+% options.PLOT_TYPE = 'Tavg-1D';ls
+
% options.PLOT_TYPE = 'Tavg-2D';
-options.NORMALIZED = 0;
+options.NORMALIZED = 1;
options.JOBNUM = 0;
options.TIME = [0 50];
options.specie = 'i';
@@ -155,21 +169,24 @@ end
if 0
%% linear growth rate for 3D Zpinch (kz fourier transform)
trange = [0.5 1]*data.Ts3D(end);
+options.INTERP = 0;
options.keq0 = 1; % chose to plot planes at k=0 or max
options.kxky = 1;
-options.kzkx = 0;
-options.kzky = 0;
+options.kzkx = 1;
+options.kzky = 1;
[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
save_figure(data,fig)
end
-if 0
+if 1
%% Mode evolution
options.NORMALIZED = 0;
-options.K2PLOT = 10;
+options.K2PLOT = 1;
options.TIME = [0:1000];
+% options.TIME = [0.5:0.01:1].*data.Ts3D(end);
options.NMA = 1;
-options.NMODES = 10;
+options.NMODES = 32;
options.iz = 'avg';
+options.ik = 1;
fig = mode_growth_meter(data,options);
save_figure(data,fig,'.png')
end
@@ -184,4 +201,4 @@ figure
plot(data.Ts3D(it0:it1), plt(data.PHI));
xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$')
title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',data.kx(ikx),data.ky(iky)))
-end
+end
\ No newline at end of file
diff --git a/wk/lin_ITG.m b/wk/lin_ITG.m
index 39cc2067266fccb80dcd89cd3521216462d031ce..2cb5704ba09b6c5a455be304b3ba08bc9800a856 100644
--- a/wk/lin_ITG.m
+++ b/wk/lin_ITG.m
@@ -11,9 +11,9 @@ addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
SIMID = 'dbg'; % Name of the simulation
-RUN = 0; % To run or just to load
+RUN = 1; % To run or just to load
default_plots_options
-% EXECNAME = 'gyacomo_dbg';
+% EXECNAME = 'gyacomo_debug';
EXECNAME = 'gyacomo_alpha';
% EXECNAME = 'gyacomo';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -29,10 +29,10 @@ K_Ni = 1*2.22; % ion Density gradient drive
K_Ti = 6.96; % ion Temperature '''
% SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
-KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
+KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 1e-4; % electron plasma beta
%% GRID PARAMETERS
-P = 20;
+P = 4;
J = P/2;
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
@@ -60,8 +60,8 @@ PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
SHIFT_Y = 0.0;
%% TIME PARMETERS
TMAX = 50; % Maximal time unit
-% DT = 1e-2; % Time step
-DT = 5e-4; % Time step
+DT = 1e-2; % Time step
+% DT = 5e-4; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = -1; % Sampling per time unit for 2D arrays
SPS3D = 1; % Sampling per time unit for 2D arrays