Skip to content
Snippets Groups Projects
Commit 7e6cccbc authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann :seedling:
Browse files

Scripts save

parent f34ee449
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 598 additions and 153 deletions
matlab/compute/compute_3D_zpinch_growth_rate.m
+ 7
− 5
View file @ 7e6cccbc
......@@ -13,6 +13,8 @@ else
end
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);
omega = zeros(nky,nkx,nkz);
......@@ -41,16 +43,16 @@ nplots = OPTIONS.kxky + OPTIONS.kzkx + OPTIONS.kzky;
iplot = 1;
if OPTIONS.kxky
[X_XY,Y_XY] = meshgrid(ky(posky),kx(poskx));
[Y_XY,X_XY] = meshgrid(ky(posky),kx(poskx));
subplot(1,nplots,iplot)
if ~OPTIONS.keq0
toplot = squeeze(max(real(omega(:,:,:)),[],3));
pclr= pcolor(X_XY,Y_XY,toplot); set(pclr,'EdgeColor','none');
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');
pclr= pcolor(X_XY,Y_XY,toplot'); set(pclr,'EdgeColor','none');
xlabel('$k_x$'); ylabel('$k_y$');
title('$\gamma(k_z=0)$');
end
......@@ -65,7 +67,7 @@ subplot(1,nplots,iplot)
xlabel('$k_x$'); ylabel('$k_y$');
title('$\max(\gamma)_{kx}$');
else
toplot = squeeze(real(omega(kxeq0,:,:)));
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)$');
......@@ -81,7 +83,7 @@ subplot(1,nplots,iplot)
xlabel('$k_z$'); ylabel('$k_x$');
title('$\max(\gamma)_{ky}$');
else
toplot = squeeze(real(omega(:,kyeq0,:)));
toplot = squeeze(real(omega(kyeq0,:,:)));
pclr= pcolor(Z_ZY,Y_ZY,toplot'); set(pclr,'EdgeColor','none');
xlabel('$k_z$'); ylabel('$k_x$');
title('$\gamma(k_y=0)$');
......
matlab/extract_fig_data.m
+ 2
− 2
View file @ 7e6cccbc
......@@ -4,7 +4,7 @@
% tw = [3000 4000];
% tw = [4000 4500];
% tw = [4500 5000];
tw = [00 5000];
tw = [00 1.3];
fig = gcf;
axObjs = fig.Children;
......@@ -22,7 +22,7 @@ end
mvm = @(x) movmean(x,1);
shift = 0;%X_(n0);
% shift = 0;
skip = 50;
skip = 1;
figure;
plot(mvm(X_(n0:skip:n1)-shift),mvm(Y_(n0:skip:n1))); hold on;
......
matlab/plot/real_plot_1D.m
+ 1
− 1
View file @ 7e6cccbc
......@@ -154,7 +154,7 @@ subplot(1,nplots,3)
legend('show')
otherwise
for it = 1:numel(toplot.FRAMES)
Y = compx(toplot.FIELD(:,:,toplot.FRAMES(it)));
Y = compx(toplot.FIELD(:,:,it));
Y = squeeze(Y);
if options.NORM
Y = Y./max(abs(Y));
......
matlab/plot/statistical_transport_averaging.m
+ 1
− 1
View file @ 7e6cccbc
function [ fig ] = statistical_transport_averaging( data, options )
function [ fig, Gx_infty_avg, Gx_infty_std ] = statistical_transport_averaging( data, options )
scale = data.scale;
Trange = options.T;
[~,it0] = min(abs(Trange(1)-data.Ts0D));
......
matlab/setup.m
+ 2
− 2
View file @ 7e6cccbc
......@@ -66,10 +66,10 @@ switch CO
case 'LR'
COLL.mat_file = '''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5''';
case 'LD'
% COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5''';
COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5''';
% COLL.mat_file = '''../../../iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5''';
% COLL.mat_file = '''../../../iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_30.h5''';
COLL.mat_file = '''../../../iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5''';
% COLL.mat_file = '''../../../iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5''';
end
COLL.coll_kcut = COLL_KCUT;
% Time integration and intialization parameters
......
scripts/scan_kN_pipeline.sh 0 → 100644
+ 59
− 0
View file @ 7e6cccbc
#! /bin/bash
# lastjid=$(sbatch submit_00.cmd)
nu_=0.01
dnu=0.0
kn_=1.7
dkn=0.2
kt_=0.425
dkt=0.05
Lx_=150
dLx=030
Tm_=4000
dTm=1000
# First submit
istart=0
lastji=$(sbatch submit_0$istart.cmd)
lastjid=${lastji##* }
echo $lastji
for i in {1..5}; do
# Setup indices of job id (current and previous one)
im1=$(awk "BEGIN {print $i-1}")
idm1=$(printf '%02d' $im1)
id=$(printf '%02d' $i)
# Create new submit file from older one
awk -v "ID=$id" '{
if (NR == 8) print "#SBATCH --error=err_"ID".txt";
else if (NR == 9) print "#SBATCH --output=out_"ID".txt";
else if (NR == 12) print "srun --cpu-bind=cores ./gyacomo 2 24 1 "ID;
else print $0}' submit_$idm1.cmd > submit_$id.cmd
# Create new fort file from older one
awk -v "NU=$nu_" -v "TM=$Tm_" -v "J2L=$im1" -v "KN=$kn_" -v "KT=$kt_" -v "LX=$Lx_" '{
if (NR == 04) print " tmax = "TM;
else if (NR == 40) print " job2load = "J2L;
else if (NR == 13) print " Lx = "LX;
else if (NR == 54) print " nu = "NU;
else if (NR == 61) print " K_Ne = "KN;
else if (NR == 62) print " K_Ni = "KN;
else if (NR == 63) print " K_Te = "KT;
else if (NR == 64) print " K_Ti = "KT;
else print $0}' fort_$idm1.90 > fort_$id.90
# Retrieve last jobid and launch next job with dep
lastji=$(sbatch --dependency=afterok:$lastjid submit_$id.cmd)
lastjid=${lastji##* }
echo $lastjid
# Increment variables
Lx_=$(awk "BEGIN {print $Lx_+$dLx}")
nu_=$(awk "BEGIN {print $nu_+$dnu}")
kn_=$(awk "BEGIN {print $kn_+$dkn}")
kt_=$(awk "BEGIN {print $kt_+$dkt}")
Tm_=$(awk "BEGIN {print $Tm_+$dTm}")
done
testcases/zpinch_example/fort.90
+ 9
− 9
View file @ 7e6cccbc
&BASIC
nrun = 100000000
dt = 0.005
tmax = 100
dt = 0.01
tmax = 200
maxruntime = 356400
/
&GRID
......@@ -58,20 +58,20 @@
sigma_i = 1
q_e = -1
q_i = 1
K_Ne = 2.22
K_Te = 6.96
K_Ni = 2.22
K_Ti = 6.96
K_Ne = 2.0
K_Te = 0.4
K_Ni = 2.0
K_Ti = 0.4
GradB = 1
CurvB = 1
lambdaD = 0
beta = 1e-4
beta = 0
/
&COLLISION_PAR
collision_model = 'SG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
gyrokin_CO = .false.
interspecies = .true.
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
!mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
&INITIAL_CON
INIT_OPT = 'phi'
......
wk/CBC_nu_CO_scan.m
+ 107
− 65
View file @ 7e6cccbc
% NU_a = [0.05 0.15 0.25 0.35 0.45];
NU_a = [0:0.1:0.5];
g_max= NU_a*0;
g_avg= NU_a*0;
g_std= NU_a*0;
k_max= NU_a*0;
CO = 'DG';
K_T = 7;
DT = 5e-3;
TMAX = 20;
ky_ = 0.3;
SIMID = 'linear_CBC_nu_scan_kT_7_ky_0.3_DGGK'; % Name of the simulation
gyacomodir = '/home/ahoffman/gyacomo/';
% EXECNAME = 'gyacomo_1.0';
% EXECNAME = 'gyacomo_dbg';
EXECNAME = 'gyacomo';
%%
NU_a = [0.005 0.01:0.01:0.1];
% P_a = [2 4 6 8 10 12 16];
% NU_a = 0.05;
P_a = [10];
CO = 'SG';
COLL_KCUT = 1.75;
K_Ti = 6.96;
DT = 1e-2;
TMAX = 100;
kymin = 0.05;
Ny = 40;
SIMID = 'linear_CBC_nu+PJ_scan_kT_6.96_SGGK'; % Name of the simulation
% SIMID = 'linear_CBC_nu_scan_kT_11_ky_0.3_DGGK'; % Name of the simulation
RUN = 1;
figure
for P = [6 8 10]
g_ky = zeros(numel(NU_a),numel(P_a),Ny/2);
g_avg= g_ky*0;
g_std= g_ky*0;
j = 1;
for P = P_a
i=1;
i = 1;
for NU = NU_a
%Set Up parameters
for j = 1
%Set Up parameters
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
TAU = 1.0; % e/i temperature ratio
K_N = 2.22; K_Ne = K_N;
K_Te = K_T; % Temperature '''
K_Ni = 2.22; K_Ne = K_Ni;
K_Te = K_Ti; % 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 = 0e-1; % electron plasma beta
J = P/2;
% J = 2;
PMAXE = P; JMAXE = J;
PMAXI = P; JMAXI = J;
NX = 12; % real space x-gridpoints
NY = 2; % '' y-gridpoints
NX = 2; % real space x-gridpoints
NY = Ny; % '' y-gridpoints
LX = 2*pi/0.1; % Size of the squared frequency domain
LY = 2*pi/ky_;
LY = 2*pi/kymin;
NZ = 16; % number of perpendicular planes (parallel grid)
NPOL = 1; SG = 0;
NPOL = 2; SG = 0;
GEOMETRY= 's-alpha';
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear (Not implemented yet)
EPS = 0.18; % inverse aspect ratio
SHEAR = 0.8; % magnetic shear
KAPPA = 0.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
SPS0D = 1; SPS2D = 0; SPS3D = 1;SPS5D= 1/5; SPSCP = 0;
JOB2LOAD= -1;
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
......@@ -62,56 +76,84 @@ for j = 1
INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
MU_Y = MU; N_HD = 4;
MU_Z = 2.0; MU_P = 0.0; %
MU_Z = 0.2; MU_P = 0.0; %
MU_J = 0.0; LAMBDAD = 0.0;
NOISE0 = 0.0e-5; % Init noise amplitude
BCKGD0 = 1.0; % Init background
GRADB = 1.0;CURVB = 1.0;
end
%%-------------------------------------------------------------------------
% RUN
setup
if RUN
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 2 1 3 0; cd ../../../wk'])
end
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
GRADB = 1.0;CURVB = 1.0;
%%-------------------------------------------------------------------------
% RUN
setup
if RUN
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 2 3 1 0; cd ../../../wk'])
end
% Load results
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
% linear growth rate (adapted for 2D zpinch and fluxtube)
options.TRANGE = [0.5 1]*data.Ts3D(end);
options.NPLOTS = 0; % 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);
% Load results
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [HELAZDIR,'results/',filename,'/'];
data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
%linear growth rate (adapted for 2D zpinch and fluxtube)
trange = [0.5 1]*data.Ts3D(end);
nplots = 0;
lg = compute_fluxtube_growth_rate(data,trange,nplots);
[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);
g_max(i) = gmax;
k_max(i) = kmax;
g_ky(i,j,:) = lg.avg_g;
g_avg(i) = lg.avg_g;
g_std(i) = lg.std_g;
g_avg(i,j,:) = lg.avg_g;
g_std(i,j,:) = lg.std_g;
i = i + 1;
end
%%
j = j + 1;
end
if 1
%% Study of the peak growth rate
figure
% plot(KT_a,max(g_max,0));
y_ = g_avg;
e_ = g_std;
% filter to noisy data
y_ = y_.*(y_-e_>0);
e_ = e_ .* (y_>0);
errorbar(NU_a,y_,e_,...
'LineWidth',1.2,...
'DisplayName',['(',num2str(P),',',num2str(P/2),')']);
hold on;
title(['$\kappa_T=$',num2str(K_T),' $k_y=$',num2str(ky_),' (CLOS = 0)']);
legend('show'); xlabel('$\nu_{DGGK}$'); ylabel('$\gamma$');
drawnow
[y_,idx_] = max(g_avg,[],3);
for i = 1:numel(idx_)
e_ = g_std(:,:,idx_(i));
end
subplot(121)
for i = 1:numel(NU_a)
errorbar(P_a,y_(i,:),e_(i,:),...
'LineWidth',1.2,...
'DisplayName',['$\nu=$',num2str(NU_a(i))]);
hold on;
end
title(['$\kappa_T=$',num2str(K_Ti),' $k_y=k_y^{max}$']);
legend('show'); xlabel('$P$, $J=P/2$'); ylabel('$\gamma$');
subplot(122)
for j = 1:numel(P_a)
errorbar(NU_a,y_(:,j),e_(:,j),...
'LineWidth',1.2,...
'DisplayName',['(',num2str(P_a(j)),',',num2str(P_a(j)/2),')']);
hold on;
end
title(['$\kappa_T=$',num2str(K_Ti),' $k_y=$',num2str(data.ky(idx_))]);
legend('show'); xlabel(['$\nu_{',CO,'}$']); ylabel('$\gamma$');
end
if 0
%% Pcolor of the peak
figure;
[XX_,YY_] = meshgrid(NU_a,P_a);
pclr=pcolor(XX_,YY_,y_'); set(pclr,'EdgeColor','none');
end
\ No newline at end of file
wk/Zpinch_coll_scan_kN_1.7.m
+ 123
− 4
View file @ 7e6cccbc
......@@ -59,25 +59,144 @@ xlabel('$\nu R/c_s$'); ylabel('$\Gamma_x^\infty/\kappa_N$');
set(gca,'Yscale','log');
end
%% KN SCAN (4,2) 200x64, muHD = 0.1 to 1.0, NL = -1, DT=1e-2 nu = 0.1
if 0
%%
figure
nu = 0.1;
% FCGK 4,2
% SGGK 4,2
kn_a = [1.60 1.80 2.00 2.20 2.40];
Gavg_a = [0.0264 0.7829 3.4742 26.3825 73.3567];
Gstd_a = [0.0039 0.1203 0.3892 4.9250 11.5890];
errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Sugama (4,2)'); hold on
% DGGK 4,2
kn_a = [1.60 1.80 2.00 2.20 2.40];
Gavg_a = [0.0080 0.0843 1.6752 28.0470 69.5804];
Gstd_a = [0.0010 0.1303 0.3329 6.0057 11.4145];
errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Dougherty (4,2)'); hold on
% LRGK 4,2
kn_a = [1.60 1.80 2.00 2.20 2.40];
Gavg_a = [1.11e-1 6.86e-1 3.44e-0 1.12e+1 2.87e+1];
Gstd_a = [7.98e-3 1.10e-1 4.03e-1 2.03e+0 7.36e+0];
Gavg_a = [0.8244 2.7805 6.4115 21.1632 62.2215];
Gstd_a = [0.0621 0.2755 0.6950 3.6594 15.0836];
errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Lorentz (4,2)'); hold on
% FCGK 4,2
kn_a = [1.60 1.80 2.00 2.20];
Gavg_a = [0.5071 3.1178 10.1663 29.3928];
Gstd_a = [0.0491 0.3237 1.1579 4.0400];
errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Coulomb (4,2)'); hold on
% LDGKii 4,2
% kn_a = [1.60 1.80 2.00 2.20 2.40];
% Gavg_a = [1.11e-1 6.86e-1 3.44e-0 1.12e+1 2.87e+1];
% Gstd_a = [7.98e-3 1.10e-1 4.03e-1 2.03e+0 7.36e+0];
% errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Landauii (4,2)'); hold on
% % Collisionless
% plot([0 1], 0.02343*[1 1],'--k','DisplayName','$\nu=0$');
%
title('Transport w.r.t. gradient level at $\nu = 0.1$');
xlim([1.6 2.5]);
legend('show');
xlabel('$\nu R/c_s$'); ylabel('$\Gamma_x^\infty/\kappa_N$');
xlabel('$L_n/R$'); ylabel('$\Gamma_x^\infty R/L_n$');
end
%% KN SCAN (4,2) 200x64, muHD = 0.1 to 1.0, NL = -1, DT=1e-2 nu = 0.01
if 0
%%
figure
Mksz = 10; Lnwd = 2.0;
nu = 0.01;
pnlty = 1;
% SGGK 4,2
kn_a = [1.60 1.80 2.00 2.20 2.40 1.70 1.90 2.10 2.30 2.50];
Gavg_a = [0.0104 0.1735 0.6430 3.8566 37.0602 pnlty*0.0433 pnlty*0.5181 pnlty*1.7723 pnlty*14.5713 pnlty*85.4020];
Gstd_a = [0.0045 0.0868 0.3587 1.8440 8.6218 pnlty*0.0214 pnlty*0.2620 pnlty*0.9196 pnlty* 4.1692 pnlty*21.5764];
[~,Idx] = sort(kn_a);
kn_a = kn_a(Idx); Gavg_a = Gavg_a(Idx); Gstd_a = Gstd_a(Idx);
% errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Sugama (4,2)'); hold on
loglog(kn_a, Gavg_a./kn_a,'s','MarkerSize',Mksz,'LineWidth',Lnwd,'DisplayName','Sugama (4,2)'); hold on
% DGGK 4,2
kn_a = [1.60 1.80 2.00 2.20 2.40 1.70 1.90 2.10 2.30 2.50];
Gavg_a = [0.0010 0.1245 0.5219 1.7153 38.6814 pnlty*0.0020 pnlty*0.3208 pnlty*1.0034 pnlty*2.6310 pnlty*79.4978];
Gstd_a = [0.0002 0.1154 0.5585 0.9062 12.814 pnlty*0.0005 pnlty*0.1853 pnlty*0.5899 pnlty*1.5473 pnlty*20.7551];
[~,Idx] = sort(kn_a);
kn_a = kn_a(Idx); Gavg_a = Gavg_a(Idx); Gstd_a = Gstd_a(Idx);
% errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Dougherty (4,2)'); hold on
loglog(kn_a, Gavg_a./kn_a,'v','MarkerSize',Mksz,'LineWidth',Lnwd,'DisplayName','Dougherty (4,2)'); hold on
% LRGK 4,2
kn_a = [1.60 1.80 2.00 2.20 2.40 1.70 1.90 2.10 2.30 2.50];
Gavg_a = [0.1422 0.5389 1.2577 4.4465 50.4417 pnlty*0.2009 pnlty*0.6764 pnlty*1.8266 pnlty*9.0066 pnlty*83.4325];
Gstd_a = [0.0119 0.1684 0.4304 2.8511 4.7634 pnlty*0.0434 pnlty*0.2677 pnlty*0.9563 pnlty*6.9812 pnlty*27.4134];
[~,Idx] = sort(kn_a);
kn_a = kn_a(Idx); Gavg_a = Gavg_a(Idx); Gstd_a = Gstd_a(Idx);
% errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Lorentz (4,2)'); hold on
loglog(kn_a, Gavg_a./kn_a,'^','MarkerSize',Mksz,'LineWidth',Lnwd,'DisplayName','Lorentz (4,2)'); hold on
plot(1,1);
% FCGK 4,2
kn_a = [1.60 1.80 2.00 2.20 2.40 1.70 1.90 2.10 2.30 2.50];
Gavg_a = [0.1914 0.8660 2.2324 5.3947 44.2023 pnlty*0.3803 pnlty*1.3516 pnlty*3.5982 pnlty*16.6212 pnlty*88.5276];
Gstd_a = [0.0097 0.1009 0.4804 1.2491 9.8378 pnlty*0.0202 pnlty*0.3688 pnlty*0.8374 pnlty*3.1974 pnlty*19.0415];
[~,Idx] = sort(kn_a);
kn_a = kn_a(Idx); Gavg_a = Gavg_a(Idx); Gstd_a = Gstd_a(Idx);
% errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Coulomb (4,2)'); hold on
loglog(kn_a, Gavg_a./kn_a,'d','MarkerSize',Mksz,'LineWidth',Lnwd,'DisplayName','Coulomb (4,2)'); hold on
%
% title('Transport w.r.t. gradient level at $\nu = 0.01$');
xlim([1.55 2.55]);
% legend('show');
xlabel('$\kappa_N$'); ylabel('$\Gamma_x^\infty /\kappa_n$');
grid on
end
if 0
%% Mixing length argument
figure;
nu = 0.0;
% collisionless
kn= [1.6 2.0 2.5 ];
k = [0.65 0.60 0.50];
g = [0.19 0.44 0.85];
y = g.^2./k.^3;
plot(kn,y*(9.6/y(3)),'--k'); hold on;
nu = 0.01;
% DG
kn= [1.6 2.0 2.5 ];
k = [0.38 0.60 0.66];
g = [0.13 0.50 0.90];
y = g.^2./k.^3;
plot(kn,y*(9.6/y(3)),'--r');
% SG
kn= [1.6 2.0 2.5 ];
k = [0.41 0.63 0.69];
g = [0.16 0.51 0.89];
y = g.^2./k.^3;
plot(kn,y*(9.6/y(3)),'--b');
% LR
kn= [1.6 2.0 2.5 ];
k = [0.41 0.57 0.60];
g = [0.18 0.50 0.87];
y = g.^2./k.^3;
plot(kn,y*(9.6/y(3)),'--y');
% LD
kn= [1.6 2.0 2.5 ];
k = [0.38 0.53 0.57];
g = [0.15 0.47 0.85];
y = g.^2./k.^3;
plot(kn,y*(9.6/y(3)),'--g');
end
\ No newline at end of file
wk/analysis_gene.m
+ 20
− 17
View file @ 7e6cccbc
helazdir = '/home/ahoffman/HeLaZ/';
addpath(genpath([helazdir,'matlab'])) % ... add
addpath(genpath([helazdir,'matlab/plot'])) % ... add
addpath(genpath([helazdir,'matlab/compute'])) % ... add
addpath(genpath([helazdir,'matlab/load'])) % ... add
gyacomodir = '/home/ahoffman/gyacomo/';
addpath(genpath([gyacomodir,'matlab'])) % ... add
addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
addpath(genpath([gyacomodir,'matlab/load'])) % ... add
% folder = '/misc/gene_results/shearless_cyclone/miller_output_1.0/';
% folder = '/misc/gene_results/shearless_cyclone/miller_output_0.8/';
......@@ -12,14 +12,17 @@ addpath(genpath([helazdir,'matlab/load'])) % ... add
% folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_0.5/';
% folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_1.0/';
% folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_0.8/';
folder = '/misc/gene_results/Z-pinch/HP_fig_2a_mu_1e-2/';
% folder = '/misc/gene_results/Z-pinch/HP_fig_2b_mu_5e-2/';
% folder = '/misc/gene_results/Z-pinch/HP_fig_2a_mu_1e-2/';
folder = '/misc/gene_results/Z-pinch/HP_fig_2b_mu_5e-2/';
% folder = '/misc/gene_results/Z-pinch/HP_fig_2c_mu_5e-2/';
% folder = '/misc/gene_results/LD_zpinch_1.6/';
% folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
% folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
% folder = '/misc/gene_results/Z-pinch/ZP_HP_kn_1.6_HRES/';
% folder = '/misc/gene_results/ZP_kn_2.5_large_box/';
% folder = '/misc/gene_results/Z-pinch/kN_2.0_HD_transport_spectrum_00/';
% folder = '/misc/gene_results/Z-pinch/kN_2.5_HD_transport_spectrum_00/';
% folder = '/misc/gene_results/CBC/128x64x16x24x12/';
% folder = '/misc/gene_results/CBC/196x96x20x32x16_02/';
% folder = '/misc/gene_results/CBC/128x64x16x6x4/';
......@@ -39,6 +42,7 @@ options.TAVG_1 = gene_data.Ts3D(end); % Averaging times duration
options.NMVA = 1; % Moving average for time traces
options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x, Q_x)
options.INTERP = 1;
options.NCUT = 4; % Number of cuts for averaging and error estimation
gene_data.FIGDIR = folder;
fig = plot_radial_transport_and_spacetime(gene_data,options);
save_figure(gene_data,fig,'.png')
......@@ -57,15 +61,15 @@ options.INTERP = 1;
options.POLARPLOT = 0;
options.AXISEQUAL = 0;
% options.NAME = 'Q_x';
options.NAME = '\phi';
% options.NAME = 'n_i';
% options.NAME = '\phi';
options.NAME = 'n_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
options.PLAN = 'kxky';
options.PLAN = 'xy';
% options.NAME ='f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
options.TIME = [325 400];
options.TIME = [1:10];
gene_data.a = data.EPS * 2000;
fig = photomaton(gene_data,options);
save_figure(gene_data,fig,'.png')
......@@ -76,12 +80,11 @@ if 0
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
% options.NAME = '\phi';
options.NAME = '\phi';
% options.NAME = 'v_y';
options.NAME = 'n_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
options.PLAN = 'xy';
options.PLAN = 'kxky';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
......@@ -118,7 +121,7 @@ end
if 0
%% Show f_i(vpar,mu)
options.times = 250:500;
options.times = 600:5000;
options.specie = 'i';
options.PLT_FCT = 'contour';
options.folder = folder;
......@@ -131,7 +134,7 @@ end
if 0
%% Time averaged spectrum
options.TIME = [4000 8000];
options.TIME = [60 10000];
options.NORM =1;
% options.NAME = '\phi';
% options.NAME = 'n_i';
......@@ -150,7 +153,7 @@ if 0
%% Mode evolution
options.NORMALIZED = 0;
options.K2PLOT = 1;
options.TIME = 100:700;
options.TIME = 1:700;
options.NMA = 1;
options.NMODES = 15;
options.iz = 'avg';
......
wk/analysis_gyacomo.m
+ 27
− 24
View file @ 7e6cccbc
%% UNCOMMENT FOR TUTORIAL
% gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory
gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory
% resdir = '.'; %Name of the directory where the results are located
% JOBNUMMIN = 00; JOBNUMMAX = 10;
%%
......@@ -31,7 +31,7 @@ if 1
i_ = 1;
disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))])
disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))])
options.TAVG_0 = data.TJOB_SE(i_);%0.4*data.Ts3D(end);
options.TAVG_0 = data.TJOB_SE(i_)+600;%0.4*data.Ts3D(end);
options.TAVG_1 = data.TJOB_SE(i_+1);%0.9*data.Ts3D(end); % Averaging times duration
options.NCUT = 4; % Number of cuts for averaging and error estimation
options.NMVA = 100; % Moving average for time traces
......@@ -40,46 +40,49 @@ options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag
options.INTERP = 0;
options.RESOLUTION = 256;
fig = plot_radial_transport_and_spacetime(data,options);
save_figure(data,fig,'.png')
% save_figure(data,fig,'.png')
end
if 0
%% statistical transport averaging
for i_ = 1:2:21
gavg =[]; gstd = [];
for i_ = 3:2:numel(data.TJOB_SE)
% i_ = 3;
disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))])
disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))])
options.T = [data.TJOB_SE(i_) data.TJOB_SE(i_+1)];
options.NPLOTS = 0;
fig = statistical_transport_averaging(data,options);
[fig, gavg_, gstd_] = statistical_transport_averaging(data,options);
gavg = [gavg gavg_]; gstd = [gstd gstd_];
end
disp(gavg); disp(gstd);
end
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
% options.NAME = '\phi';
options.NAME = '\phi';
% options.NAME = '\omega_z';
% options.NAME = 'N_i^{00}';
% options.NAME = 'v_y';
% options.NAME = 'v_x';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
options.NAME = 'n_i';
options.PLAN = 'xy';
% options.NAME = 'n_i';
options.PLAN = 'kxky';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
% options.TIME = data.Ts5D(end-30:end);
% options.TIME = data.Ts3D;
options.TIME = [000:0.1:7000];
options.TIME = [000:1:8000];
data.EPS = 0.1;
data.a = data.EPS * 2000;
options.RESOLUTION = 256;
create_film(data,options,'.gif')
end
if 1
if 0
%% 2D snapshots
% Options
options.INTERP = 1;
......@@ -96,7 +99,7 @@ options.PLAN = 'xy';
% options.NAME 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
options.TIME = [1000 1800 2500 3000 4000];
options.TIME = [1000 3000 5000 8000 10000];
data.a = data.EPS * 2e3;
fig = photomaton(data,options);
......@@ -119,16 +122,16 @@ end
if 0
%% Kinetic distribution function sqrt(<f_a^2>xy) (GENE vsp)
options.SPAR = linspace(-3,3,32)+(6/127/2);
options.XPERP = linspace( 0,6,32);
% options.SPAR = gene_data.vp';
% options.XPERP = gene_data.mu';
% options.SPAR = linspace(-3,3,32)+(6/127/2);
% options.XPERP = linspace( 0,6,32);
options.SPAR = gene_data.vp';
options.XPERP = gene_data.mu';
options.iz = 'avg';
options.T = [250 600];
options.PLT_FCT = 'pcolor';
options.ONED = 0;
options.T = [1500 5000];
options.PLT_FCT = 'contour';
options.ONED = 1;
options.non_adiab = 0;
options.SPECIE = 'i';
options.SPECIE = 'e';
options.RMS = 1; % Root mean square i.e. sqrt(sum_k|f_k|^2) as in Gene
fig = plot_fa(data,options);
% save_figure(data,fig,'.png')
......@@ -175,8 +178,8 @@ options.NAME = '\phi';
% options.NAME = 'n_i';
% options.NAME ='\Gamma_x';
% options.NAME ='s_y';
options.COMPX = 'avg';
options.COMPY = 'avg';
options.COMPX = 1;
options.COMPY = 2;
options.COMPZ = 1;
options.COMPT = 1;
options.MOVMT = 1;
......@@ -220,11 +223,11 @@ end
if 0
%% linear growth rate for 3D Zpinch
trange = [5 15];
trange = [5 20];
options.keq0 = 1; % chose to plot planes at k=0 or max
options.kxky = 1;
options.kzkx = 0;
options.kzky = 1;
[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
save_figure(data,fig,'.png')
% save_figure(data,fig,'.png')
end
wk/header_2DZP_results.m
+ 21
− 3
View file @ 7e6cccbc
......@@ -178,7 +178,7 @@ resdir ='';
% resdir ='Zpinch_rerun/Kn_2.0_256x64x9x5_Lx_240_Ly_120';
% resdir ='Zpinch_rerun/Kn_1.6_256x128x7x4';
% resdir ='Zpinch_rerun/Kn_1.6_200x64x11x6';
% resdir ='Zpinch_rerun/Kn_1.6_200x64x11x6_conv';
resdir ='Zpinch_rerun/Kn_1.6_200x64x11x6_conv';
% resdir ='Zpinch_rerun/Kn_1.6_200x64x11x6_mu_0.5';
% resdir ='Zpinch_rerun/Kn_1.6_256x128x21x11';
......@@ -197,16 +197,18 @@ resdir ='';
% % 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';
%% Convergence cases kN = (1.6 2.2) nu = (0.01 1.0)
%% Convergence cases kN = (1.6 2.2) nu = (0.01 1.0), SGGK
% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x3x2_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x5x3_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x7x4_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x9x5_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x11x6_nu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x3x2_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x5x3_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x7x4_mu_0.01';
resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x9x5_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x9x5_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x3x2_nu_0.1';
% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x5x3_nu_0.1';
......@@ -218,6 +220,22 @@ resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x9x5_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x7x4_nu_0.1';
% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x9x5_nu_0.1';
%% kN scans with nu=0.1 and nu=0.01
% resdir = 'Zpinch_rerun/SGGK_kN_scan_200x64x5x3_nu_0.1';
% resdir = 'Zpinch_rerun/DGGK_kN_scan_200x64x5x3_nu_0.1';
% resdir = 'Zpinch_rerun/LRGK_kN_scan_200x64x5x3_nu_0.1';
% resdir = 'Zpinch_rerun/LDGK_kN_scan_200x64x5x3_nu_0.1';
%
% resdir = 'Zpinch_rerun/SGGK_kN_scan_200x64x5x3_nu_0.01';
% 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/DGGK_kN_1.6_200x64x5x3_nu_0.01';
% resdir = 'Zpinch_rerun/DGGK_kN_1.7_200x64x5x3_nu_0.01';
resdir = 'Zpinch_rerun/LRGK_kN_2.4_300x98x5x3_nu_0.01';
% resdir = 'Zpinch_rerun/LDGK_kN_1.9_200x64x5x3_nu_0.01';
JOBNUMMIN = 00; JOBNUMMAX = 03;
resdir = ['results/',resdir];
......
wk/lin_3D_Zpinch.m 0 → 100644
+ 194
− 0
View file @ 7e6cccbc
%% QUICK RUN SCRIPT
% 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"
%
% SIMID = 'lin_3D_Zpinch'; % Name of the simulation
SIMID = 'NL_3D_zpinch'; % Name of the simulation
% 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_1.0';
EXECNAME = 'gyacomo_dbg';
% EXECNAME = 'gyacomo';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
NU = 0.05; % Collision frequency
TAU = 1.0; % e/i temperature ratio
K_Ne = 2.0; % ele Density '''
K_Te = K_Ne/4; % ele Temperature '''
K_Ni = K_Ne; % ion Density gradient drive
K_Ti = K_Ni/4; % 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
BETA = 0.005; % electron plasma beta
%% GRID PARAMETERS
P = 2;
J = P/2;
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
NX = 128; % real space x-gridpoints
NY = 32; % '' y-gridpoints
LX = 2*pi/0.1; % Size of the squared frequency domain
LY = 2*pi/0.1; % Size of the squared frequency domain
NZ = 20; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
%% GEOMETRY
% GEOMETRY= 's-alpha';
% GEOMETRY= 'miller';
GEOMETRY= 'Z-pinch';
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear
KAPPA = 0.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
%% TIME PARMETERS
TMAX = 20; % Maximal time unit
DT = 1e-2; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
SPS3D = 5; % Sampling per time unit for 2D arrays
SPS5D = 1/5; % 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 = 1; % 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
%% 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 = 2.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 = 1.0;
CURVB = 1.0;
COLL_KCUT = 1000;
%%-------------------------------------------------------------------------
%% RUN
setup
% system(['rm fort*.90']);
% Run linear simulation
if RUN
if NZ > 1
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk'])
else
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
end
end
%% Load results
%%
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [HELAZDIR,'results/',filename,'/'];
% Load outputs from jobnummin up to jobnummax
JOBNUMMIN = 00; JOBNUMMAX = 00;
data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
%% Short analysis
if 1
%% 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 1
%% Ballooning plot
options.time_2_plot = [120];
options.kymodes = [0.3];
options.normalized = 1;
% options.field = 'phi';
fig = plot_ballooning(data,options);
end
if 0
%% Hermite-Laguerre spectrum
% options.TIME = 'avg';
options.P2J = 1;
options.ST = 1;
options.PLOT_TYPE = 'space-time';
% options.PLOT_TYPE = 'Tavg-1D';
% options.PLOT_TYPE = 'Tavg-2D';
options.NORMALIZED = 0;
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 0
%% linear growth rate for 3D Zpinch (kz fourier transform)
trange = [0.5 1]*data.Ts3D(end);
options.keq0 = 1; % chose to plot planes at k=0 or max
options.kxky = 1;
options.kzkx = 0;
options.kzky = 0;
[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
save_figure(data,fig)
end
if 0
%% Mode evolution
options.NORMALIZED = 0;
options.K2PLOT = 1;
options.TIME = [0:1000];
options.NMA = 1;
options.NMODES = 1;
options.iz = 'avg';
fig = mode_growth_meter(data,options);
save_figure(data,fig,'.png')
end
if 1
%% 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
wk/lin_EPY.m → wk/lin_ETPY.m
+ 18
− 14
View file @ 7e6cccbc
%% QUICK RUN SCRIPT for linear entropy mode in a Zpinch
%% 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"
%
SIMID = 'lin_EPY'; % Name of the simulation
SIMID = 'lin_ETPY'; % Name of the simulation
% SIMID = 'dbg'; % Name of the simulation
RUN = 0; % To run or just to load
addpath(genpath('../matlab')) % ... add
......@@ -16,9 +16,9 @@ EXECNAME = 'gyacomo';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
NU = 0.1; % Collision frequency
NU = 0.01; % Collision frequency
TAU = 1.0; % e/i temperature ratio
K_Ne = 1.6; % ele Density '''
K_Ne = 2.5; % ele Density '''
K_Te = K_Ne/4; % ele Temperature '''
K_Ni = K_Ne; % ion Density gradient drive
K_Ti = K_Ni/4; % ion Temperature '''
......@@ -26,7 +26,7 @@ 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 = 2;
P = 4;
J = P/2;
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
......@@ -35,7 +35,7 @@ JMAXI = J; % "
NX = 2; % real space x-gridpoints
NY = 100; % '' y-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
LY = 120;%2*pi/0.05; % Size of the squared frequency domain
LY = 200;%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
......@@ -43,14 +43,17 @@ SG = 0; % Staggered z grids option
GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear
KAPPA = 0.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
%% TIME PARMETERS
TMAX = 500; % Maximal time unit
DT = 1e-2; % Time step
DT = 5e-2; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
SPS3D = 1/5; % 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
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
......@@ -58,7 +61,7 @@ JOB2LOAD= -1;
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
CO = 'SG';
CO = 'LD';
GKCO = 1; % gyrokinetic operator
ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
......@@ -85,18 +88,19 @@ MU_Z = 2.0; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 0.0e-5; % Init noise amplitude
BCKGD0 = 1.0; % Init background
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
GRADB = 1.0;
CURVB = 1.0;
COLL_KCUT = 1000;
%%-------------------------------------------------------------------------
%% RUN
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 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'])
end
%% Load results
......@@ -157,7 +161,7 @@ options.kzky = 0;
[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
save_figure(data,fig)
end
if 1
if 0
%% Mode evolution
options.NORMALIZED = 0;
options.K2PLOT = 10;
......
wk/lin_ITG.m
+ 7
− 6
View file @ 7e6cccbc
......@@ -7,7 +7,7 @@ SIMID = 'lin_ITG'; % Name of the simulation
RUN = 0; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
HELAZDIR = '/home/ahoffman/gyacomo/';
gyacomodir = '/home/ahoffman/gyacomo/';
% EXECNAME = 'gyacomo_1.0';
% EXECNAME = 'gyacomo_dbg';
EXECNAME = 'gyacomo';
......@@ -103,14 +103,14 @@ if RUN
% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk'])
system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
end
%% Load results
%%
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [HELAZDIR,'results/',filename,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
% Load outputs from jobnummin up to jobnummax
JOBNUMMIN = 00; JOBNUMMAX = 00;
data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
......@@ -118,9 +118,10 @@ data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from
%% Short analysis
if 1
%% linear growth rate (adapted for 2D zpinch and fluxtube)
trange = [0.5 1]*data.Ts3D(end);
nplots = 3;
lg = compute_fluxtube_growth_rate(data,trange,nplots);
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.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);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment