diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index 79df88e435f8f8d08c229955b31bd73272b68942..0dd9e6b5aaac566c1b432ff2a5fefbfe18141b33 100644
--- a/matlab/extract_fig_data.m
+++ b/matlab/extract_fig_data.m
@@ -4,7 +4,7 @@
% tw = [3000 4000];
% tw = [4000 4500];
% tw = [4500 5000];
-tw = [100 1000];
+tw = [200 1000];
fig = gcf;
axObjs = fig.Children;
diff --git a/matlab/process_field.m b/matlab/process_field.m
index 9ead884c631372d4790cf003240b3e5aebbad2df..6e9b66452d3b43d5df2198042cfbc53330ad40ca 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -201,14 +201,22 @@ switch OPTIONS.NAME
OPE_ = 1;
FLD_ = DATA.DENS_I(:,:,:,FRAMES)...
-DATA.DENS_E(:,:,:,FRAMES);
- case 'u_i' % ion density
- NAME = 'ui';
+ case 'upar_i' % ion density
+ NAME = 'upari';
FLD_ = DATA.UPAR_I(:,:,:,FRAMES);
OPE_ = 1;
- case 'u_e' % electron density
- NAME = 'ue';
+ case 'upar_e' % electron density
+ NAME = 'upare';
FLD_ = DATA.UPAR_E(:,:,:,FRAMES);
OPE_ = 1;
+ case 'uper_i' % ion density
+ NAME = 'uperi';
+ FLD_ = DATA.UPER_I(:,:,:,FRAMES);
+ OPE_ = 1;
+ case 'uper_e' % electron density
+ NAME = 'upere';
+ FLD_ = DATA.UPER_E(:,:,:,FRAMES);
+ OPE_ = 1;
case 'T_i' % ion temperature
NAME = 'Ti';
FLD_ = DATA.TEMP_I(:,:,:,FRAMES);
diff --git a/matlab/setup.m b/matlab/setup.m
index 2cdb3cc20077c9614bf5316cf2a133a4dfcb501c..b8850da41972a321e0e46789c475754a62dc4d6a 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -42,6 +42,11 @@ MODEL.ExBrate = EXBRATE;
MODEL.mu_x = MU_X;
MODEL.mu_y = MU_Y;
MODEL.N_HD = N_HD;
+try
+ MODEL.N_HDz = N_HDz;
+catch
+ MODEL.N_HDz = 4;
+end
MODEL.mu_z = MU_Z;
MODEL.HYP_V = ['''',HYP_V,''''];
MODEL.mu_p = MU_P;
diff --git a/src/miller_mod.F90 b/src/miller_mod.F90
index 994022031f6d31b3235411c1e95f30ab9fe8619a..5df1d153fa1c8f9a2b8c79bea474f06f5950c6b7 100644
--- a/src/miller_mod.F90
+++ b/src/miller_mod.F90
@@ -6,7 +6,7 @@ MODULE miller
USE basic
USE parallel, ONLY: my_id, num_procs_z, nbr_U, nbr_D, comm0
! use coordinates,only: gcoor, get_dzprimedz
- USE grid, ONLY: local_Nky, local_Nkx, local_nz, ngz, nzgrid, kyarray, kxarray, zarray, total_nz, local_nz_offset, iodd, ieven
+ USE grid, ONLY: local_Nky, local_Nkx, local_nz, ngz, nzgrid, kyarray, kxarray, zarray, total_nz, local_nz_offset, iodd, ieven, deltaz
! use discretization
USE lagrange_interpolation
@@ -516,6 +516,7 @@ CONTAINS
!(R,Z) derivatives for visualization
dxdR_s = dx_drho/drPsi*psi1_s*cosu_s
dxdZ_s = dx_drho/drPsi*psi1_s*sinu_s
+
! GHOSTS ADAPTED VERSION
if (edge_opt==0._xp) then
!gyacomo z-grid wo ghosts
@@ -524,45 +525,49 @@ CONTAINS
chi_out=linspace(-pi*Npol-4._xp*pi*Npol/total_nz,pi*Npol+2._xp*pi*Npol/total_nz,total_nz+ngz)
else
ERROR STOP '>> ERROR << ghosts not implemented for edge_opt yet'
- !new parallel coordinate chi_out==zprime
- !see also tracer_aux.F90
- if (Npol>1) ERROR STOP '>> ERROR << Npol>1 has not been implemented for edge_opt=\=0._xp'
- do k=1,total_nz
+ !new parallel coordinate chi_out==zprime
+ !see also tracer_aux.F90
+ if (Npol>1) ERROR STOP '>> ERROR << Npol>1 has not been implemented for edge_opt=\=0._xp'
+ do k=1,total_nz
chi_out(k)=sinh((-pi+k*2._xp*pi/total_nz)*log(edge_opt*pi+sqrt(edge_opt**2*pi**2+1))/pi)/edge_opt
- enddo
- !transform metrics according to chain rule
- do k=1,np_s
- !>dz'/dz conversion for edge_opt as function of z
+ enddo
+ !transform metrics according to chain rule
+ do k=1,np_s
+ !>dz'/dz conversion for edge_opt as function of z
if (edge_opt.gt.0) then
- dzprimedz = edge_opt*pi/log(edge_opt*pi+sqrt((edge_opt*pi)**2+1._xp))/&
+ dzprimedz = edge_opt*pi/log(edge_opt*pi+sqrt((edge_opt*pi)**2+1._xp))/&
sqrt((edge_opt*chi_s(k))**2+1)
else
- dzprimedz = 1.0
+ dzprimedz = 1.0
endif
gxz(k)=gxz(k)*dzprimedz
gyz(k)=gyz(k)*dzprimedz
gzz(k)=gzz(k)*dzprimedz**2
jacobian(k)=jacobian(k)/dzprimedz
dBdz(k)=dBdz(k)/dzprimedz
- enddo
+ enddo
endif !edge_opt
- ! interpolate with ghosts
- call lag3interp(gxx,chi_s,np_s,gxx_out,chi_out,total_nz+ngz)
- call lag3interp(gxy,chi_s,np_s,gxy_out,chi_out,total_nz+ngz)
- call lag3interp(gxz,chi_s,np_s,gxz_out,chi_out,total_nz+ngz)
- call lag3interp(gyy,chi_s,np_s,gyy_out,chi_out,total_nz+ngz)
- call lag3interp(gyz,chi_s,np_s,gyz_out,chi_out,total_nz+ngz)
- call lag3interp(gzz,chi_s,np_s,gzz_out,chi_out,total_nz+ngz)
- call lag3interp(B_s,chi_s,np_s,Bfield_out,chi_out,total_nz+ngz)
- call lag3interp(dBdx,chi_s,np_s,dBdx_out,chi_out,total_nz+ngz)
- call lag3interp(dBdz,chi_s,np_s,dBdz_out,chi_out,total_nz+ngz)
- call lag3interp(jacobian,chi_s,np_s,jacobian_out,chi_out,total_nz+ngz)
- call lag3interp(R_s,chi_s,np_s,R_out,chi_out,total_nz+ngz)
- call lag3interp(Z_s,chi_s,np_s,Z_out,chi_out,total_nz+ngz)
- call lag3interp(dxdR_s,chi_s,np_s,dxdR_out,chi_out,total_nz+ngz)
- call lag3interp(dxdZ_s,chi_s,np_s,dxdZ_out,chi_out,total_nz+ngz)
- ! Fill the local geom arrays with the results
- do eo=iodd,ieven
+
+ !! Loop over the z-grids and interpolate the results on it
+ do eo=ieven,iodd
+ ! Shift the grid
+ chi_out = chi_out + REAL(eo-1,xp)*0.5*deltaz
+ ! interpolate with ghosts
+ call lag3interp(gxx,chi_s,np_s,gxx_out,chi_out,total_nz+ngz)
+ call lag3interp(gxy,chi_s,np_s,gxy_out,chi_out,total_nz+ngz)
+ call lag3interp(gxz,chi_s,np_s,gxz_out,chi_out,total_nz+ngz)
+ call lag3interp(gyy,chi_s,np_s,gyy_out,chi_out,total_nz+ngz)
+ call lag3interp(gyz,chi_s,np_s,gyz_out,chi_out,total_nz+ngz)
+ call lag3interp(gzz,chi_s,np_s,gzz_out,chi_out,total_nz+ngz)
+ call lag3interp(B_s,chi_s,np_s,Bfield_out,chi_out,total_nz+ngz)
+ call lag3interp(dBdx,chi_s,np_s,dBdx_out,chi_out,total_nz+ngz)
+ call lag3interp(dBdz,chi_s,np_s,dBdz_out,chi_out,total_nz+ngz)
+ call lag3interp(jacobian,chi_s,np_s,jacobian_out,chi_out,total_nz+ngz)
+ call lag3interp(R_s,chi_s,np_s,R_out,chi_out,total_nz+ngz)
+ call lag3interp(Z_s,chi_s,np_s,Z_out,chi_out,total_nz+ngz)
+ call lag3interp(dxdR_s,chi_s,np_s,dxdR_out,chi_out,total_nz+ngz)
+ call lag3interp(dxdZ_s,chi_s,np_s,dxdZ_out,chi_out,total_nz+ngz)
+ ! Fill the local geom arrays with the results
DO iz = 1,local_nz + ngz
gxx_(iz,eo) = gxx_out(iz+local_nz_offset)
gxy_(iz,eo) = gxy_out(iz+local_nz_offset)
@@ -580,7 +585,7 @@ CONTAINS
dxdR_(iz,eo) = dxdR_out(iz+local_nz_offset)
dxdZ_(iz,eo) = dxdZ_out(iz+local_nz_offset)
ENDDO
- ENDDO
+ ENDDO
contains
!> Generate an equidistant array from min to max with n points
function linspace(min,max,n) result(out)
diff --git a/wk/analysis_gbms.m b/wk/analysis_gbms.m
index 376ae1b2561e698b6a6e336ce22cecbe0fcceb60..936a5c200d7cbe8e75ed0f04604ab2681b45efc2 100644
--- a/wk/analysis_gbms.m
+++ b/wk/analysis_gbms.m
@@ -9,9 +9,9 @@
% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/RH_test_kine/';
% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/KBM/';
% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/TEM/';
-resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/ITG/';
+% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/ITG/';
% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/linear_cyclone/';
-% resdir = '/home/ahoffman/molix/';
+resdir = '/home/ahoffman/Documents/gbms/results/MTM/';
% system(['cd ',resdir,';','./gbms < parameters.in; cd /home/ahoffman/HeLaZ/wk']);
outfile = [resdir,'field.dat.h5'];
diff --git a/wk/fast_analysis.m b/wk/fast_analysis.m
index 490a4f926aa4c094b8c4000d0cb0bdb70adc4bf6..c89503f713eddf82ca4d957a11b1780f514ea75e 100644
--- a/wk/fast_analysis.m
+++ b/wk/fast_analysis.m
@@ -6,20 +6,27 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add
default_plots_options
% Partition of the computer where the data have to be searched
% PARTITION='/Users/ahoffmann/gyacomo/results/';
-PARTITION='/home/ahoffman/gyacomo/results/';
-% PARTITION='/misc/gyacomo23_outputs/paper_3/';
+% PARTITION='/home/ahoffman/gyacomo/results/';
+% PARTITION='/misc/gyacomo23_outputs/';
% resdir = 'AE_3x2x128x32x24/PT';
% resdir = 'AE_3x2x128x32x24/PT';
% resdir = 'AE_5x3x128x32x24/NT';
-resdir = 'IS_5x3x128x32x24/PT';
+% resdir = 'IS_5x3x128x32x24/NT';
+% PARTITION='/home/ahoffman/gyacomo/simulations/';
+% resdir ='cheap_CBC_baseline';
+% resdir ='test_HEL_closure';
+% resdir ='dmax_closure/';
+PARTITION = '/misc/gyacomo23_outputs/reduced_fluid_paper/';
+% resdir = '/Npol_study/AE_CBC_s0_beta0_P4J2';
+resdir = '/Npol_study/RF_CBC_s0_beta0/Npol_11';
% Triangularity paper
-PARTITION = '/misc/gyacomo23_outputs/triangularity_paper/';
+% PARTITION = '/misc/gyacomo23_outputs/triangularity_paper/';
% Nominal parameters
% resdir = 'ion_scale/3x2x256x64x32/0T';
% resdir = 'ion_scale/5x3x256x64x32/0T';
% resdir = 'ion_scale/5x3x192x48x24/0T';
-resdir = 'ion_scale/9x5x256x64x32/0T';
+% resdir = 'ion_scale/9x5x256x64x32/0T';
% resdir = 'ion_scale/restart/5x3x256x64x32/0T';
% resdir = 'ion_scale/restart/9x5x192x48x24/0T';
% resdir = 'adiabatic_electrons/5x2x256x64x32/0T';
@@ -48,7 +55,7 @@ resdir = 'ion_scale/9x5x256x64x32/0T';
DATADIR = [PARTITION,resdir,'/'];
read_flux_out_XX(DATADIR,1,1);
%%
-J0 = 00; J1 = 10;
+J0 = 00; J1 = 00;
% Load basic info (grids and time traces)
data = {};
@@ -60,13 +67,18 @@ data.Ne00 = reshape(data.Na00(2,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grid
catch
end
[data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
+if data.inputs.BETA > 0
+ [data.PSI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'psi');
+end
if 1
%%
[data.TEMP, data.Ts3D] = compile_results_3Da(data.folder,J0,J1,'temp');
-% [data.UPAR, data.Ts3D] = compile_results_3Da(data.folder,J0,J1,'upar');
+[data.UPAR, data.Ts3D] = compile_results_3Da(data.folder,J0,J1,'upar');
+[data.UPER, data.Ts3D] = compile_results_3Da(data.folder,J0,J1,'uper');
[data.DENS, data.Ts3D] = compile_results_3Da(data.folder,J0,J1,'dens');
data.TEMP_I = reshape(data.TEMP(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
-% data.UPAR_I = reshape(data.UPAR(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+data.UPAR_I = reshape(data.UPAR(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
+data.UPER_I = reshape(data.UPER(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
data.DENS_I = reshape(data.DENS(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
if data.inputs.Na > 1
@@ -75,7 +87,7 @@ if data.inputs.Na > 1
data.Ne00 = reshape(data.Na00(2,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
end
end
-if 0
+if 1
%% Plot transport and phi radial profile
% [data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
% [data.PSI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'psi');
@@ -102,33 +114,33 @@ end
if 0
%% 2D field snapshots
% Options
-options.INTERP = 0;
+options.INTERP = 1;
options.POLARPLOT = 0;
options.AXISEQUAL = 0;
options.NORMALIZE = 0;
options.LOGSCALE = 0;
options.CLIMAUTO = 1;
options.TAVG = 1;
-options.NAME = ['N_i^{00}'];
+% options.NAME = ['N_i^{00}'];
% options.NAME = 'n_e';
-% options.NAME = 'u_i';
+% options.NAME = 'upar_i';
% options.NAME = 'n_i';
% options.NAME = 'Q_{xi}';
% options.NAME = 'v_{Ey}';
% options.NAME = 'w_{Ez}';
% options.NAME = '\omega_z';
-% options.NAME = '\phi';
+options.NAME = '\phi';
% options.NAME = 'n_i-n_e';
loc =11;
[~,i_] = min(abs(loc - data.grids.y));
options.COMP =i_;
% options.PLAN = '3D';
% options.PLAN = 'xy'; options.COMP =floor(data.grids.Nz/2)+1;
-options.PLAN = 'xz'; options.COMP ='avg';
+options.PLAN = 'kxz'; options.COMP ='avg';
% options.COMP ='avg';
options.XYZ =[-11 20 0];
% options.TIME = [100 250 350 500]; options.TAVG = 0;
-options.TIME = [100:500]; options.TAVG = 1;
+options.TIME = [50:500]; options.TAVG = 1;
options.RESOLUTION = 256;
fig = photomaton(data,options);
% colormap(gray)
@@ -200,7 +212,7 @@ data.Napjz(1,3,1,:,:) = data.Napjz(1,3,1,:,:)*data.inputs.tau(1);
data.Napjz(1,1,2,:,:) = data.Napjz(1,1,2,:,:)*data.inputs.tau(1);
% [data.Napjz, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'Nipjz');
options.ST = 1;
-options.NORMALIZED = 1;
+options.NORMALIZED = 0;
options.LOGSCALE = 1;
options.FILTER = 0; %filter the 50% time-average of the spectrum from
options.TAVG_2D = 0; %Show a 2D plot of the modes, 50% time averaged
@@ -238,6 +250,20 @@ options.NORMALIZE = 0;
[fig] = plot_spectrum(data,options);
end
+if 0
+%% 1D radial plot
+options.TIME = [100 300]; % averaging time window
+options.NAME = ['N_i^{00}'];
+% options.NAME = 'n_i';
+% options.NAME = 'T_i';
+% options.NAME = 'Q_{xi}';
+% options.NAME = 's_{Ey}';
+% options.NAME = '\phi';
+% options.NAME = '\psi';
+options.NORMALIZE = 0;
+[fig] = plot_spectrum(data,options);
+end
+
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -255,13 +281,13 @@ options.INTERP = 1;
options.POLARPLOT = 0;
options.BWR = 1; % bluewhitered plot or gray
options.CLIMAUTO = 0; % adjust the colormap auto
-% options.NAME = '\phi';
+options.NAME = '\phi';
% options.NAME = 'w_{Ez}';
% options.NAME = '\psi';
-% options.NAME = 'n_i';
+% options.NAME = 'T_i';
% options.NAME = '\phi^{NZ}';
+% options.NAME = ['N_i^{00}'];
% options.NAME = ['N_e^{00}'];
-options.NAME = ['N_i^{00}'];
options.PLAN = 'xy'; options.COMP =floor(data.grids.Nz/2)+1;
% options.PLAN = 'xz'; options.COMP ='avg';
% options.PLAN = '3D';
diff --git a/wk/heat_flux_to_MW_converter b/wk/heat_flux_to_MW_converter
new file mode 100644
index 0000000000000000000000000000000000000000..6933768a1913ba18859c26ae3888e4dc871e2d42
--- /dev/null
+++ b/wk/heat_flux_to_MW_converter
@@ -0,0 +1,21 @@
+n0 = 2.5e19; % density (electrons) [1e19m-3]
+T0 = 0.3; % temperature (electrons) [keV]
+R0 = 1.7; % length (major radius) [m]
+eps= 0.3; % []
+B0 = 2.0; % magnetic field [T]
+m0 = 1.7e-27; % mass [nucleid mass]
+q0 = 1.6e-19; % charge [elementary charge]
+kB = 1.4e-23; % [J/K];
+%
+T0 = q0*(1000*T0)/kB; % [K]
+S0 = 4*pi^2*eps*R0^2;
+%
+p0 = n0*kB*T0; % [Pa]
+c0 = sqrt(kB*T0/m0); % [m/s]
+w0 = q0*B0/m0; %[s]
+r0 = c0/w0; %[m]
+Q0 = p0*c0*r0^2/R0^2; % [W/m^2]
+
+Qxi = 20; Qxe = 10;
+P0 = Q0*S0*(Qxi+Qxe)/1000; % output power in MW
+disp([num2str(P0),' MW'])
\ No newline at end of file
diff --git a/wk/lin_run_script.m b/wk/lin_run_script.m
index 402d2cafe4c628b67c1751864d8e402f85f1d872..61a07f4a59367363ee5bb1154d17d651756a2e6d 100644
--- a/wk/lin_run_script.m
+++ b/wk/lin_run_script.m
@@ -26,12 +26,13 @@ EXECNAME = 'gyacomo23_dp'; % double precision
% run lin_DTT_HM_rho85
% run lin_DTT_HM_rho98
% run lin_DIIID_LM_rho90
-run lin_DIIID_LM_rho95
+% run lin_DIIID_LM_rho95
% run lin_DIIID_LM_rho95_HEL
% run lin_JET_rho97
% run lin_Entropy
% run lin_ITG
% run lin_KBM
+run lin_MTM_Hamed
% run lin_RHT
% run lin_Ivanov
% rho = 0.95; TRIANG = 'NT'; READPROF = 1;
@@ -41,32 +42,32 @@ run lin_DIIID_LM_rho95
% run lin_DIIID_data
% run lin_STEP_EC_HD_psi71
% run lin_STEP_EC_HD_psi49
-if 0
+% if 1
% Plot the profiles
- plot_params_vs_rho(geom,prof_folder,rho,0.5,1.1,Lref,mref,Bref,READPROF);
-end
+% plot_params_vs_rho(geom,prof_folder,rho,0.01,1.1,Lref,mref,Bref,READPROF);
+% end
% SIMID = ['rho_scan_DIIID_AUSTIN_2019/3x2x192x96x32/rho',num2str(rho)];
%% Change parameters
% GEOMETRY = 's-alpha';
-PMAX = 8; JMAX = 4;
+% PMAX = 2; JMAX = 1;
% DELTA = 0.2;
% K_tB = 0; K_mB = 0; K_ldB = 1;
% K_Ni = 0;
% K_Ne = 0;
% K_Ti = 0;
% K_Te = 0;
-DELTA = 0.0;
+% DELTA = 0.0;
% DELTA = 0.2;
-S_DELTA = DELTA/2;
-LY = 2*pi/0.1;
-TMAX = 40;
-NY = 20;
+% S_DELTA = DELTA/2;
+LY = 2*pi/0.7;
+% TMAX = 40;
+% NY = 20;
DT = 0.001;
-CO = 'DG';
+% CO = 'DG';
% TAU = 1; NU = 0.05;
-NA = 1; ADIAB_E = 1; DT = 1e-2; DTSAVE3D = 1e-2; TMAX = 60;
-TAU = 1e-3; K_Ti = K_Ti/2/TAU; K_Ni = 0;
-NU = 3*NU/8/TAU; PMAX = 2; JMAX = 1;
+% NA = 1; ADIAB_E = 1; DT = 1e-2; DTSAVE3D = 1e-2; TMAX = 60;
+% TAU = 1e-3; K_Ti = K_Ti/2/TAU; K_Ni = 0;
+% NU = 3*NU/8/TAU; PMAX = 2; JMAX = 1;
% K_Ti = K_Ti/4;
% MU_X = 1; MU_Y = 1;
%% RUN
@@ -77,8 +78,8 @@ system(['cp ../results/',SIMID,'/',PARAMS,'/fort_00.90 ../results/',SIMID,'/.'])
if RUN
MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
% RUN =['time ',mpirun,' -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
- RUN =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 4 1 0;'];
- % RUN =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 2 2 1 0;'];
+ % RUN =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 4 1 0;'];
+ RUN =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
% RUN =['time ',mpirun,' -np 8 ',gyacomodir,'bin/',EXECNAME,' 2 2 2 0;'];
% RUN =['time ',mpirun,' -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
% RUN = ['./../../../bin/gyacomo23_sp 0;'];
@@ -92,7 +93,7 @@ filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARA
LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
FIGDIR = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
% Load outputs from jobnummin up to jobnummax
-J0 = 0; J1 = 0;
+J0 = 00; J1 = 00;
data = {}; % Initialize data as an empty cell array
% load grids, inputs, and time traces
data = compile_results_low_mem(data,LOCALDIR,J0,J1);
@@ -117,11 +118,11 @@ options.SHOWFIG = 1;
[fig, wkykx, ekykx] = mode_growth_meter(data,options);
% %%
% kx = (1:data.grids.Nx/2)'*2*pi/data.fort_00.GRID.Lx;
-ky = (1:data.grids.Ny/2)'*2*pi/data.fort_00.GRID.Ly;
-gkxky = real(wkykx(2:end,1:data.grids.Nx/2))';
-gkxky(isnan(gkxky)) =0;
-gkxky(isinf(gkxky)) =0;
-figure; plot(ky,gkxky(1,:));
+% ky = (1:data.grids.Ny/2)'*2*pi/data.fort_00.GRID.Ly;
+% gkxky = real(wkykx(2:end,1:data.grids.Nx/2))';
+% gkxky(isnan(gkxky)) =0;
+% gkxky(isinf(gkxky)) =0;
+% figure; plot(ky,gkxky(1,:));
end
if (0 && NZ>4)
@@ -162,52 +163,21 @@ if 0
plot_metric(data);
end
-if 0
- %% Compiled plot for lin EM analysis
- [data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
- if data.inputs.BETA > 0
- [data.PSI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'psi');
- end
- options.time_2_plot = [120];
- options.kymodes = [100];
- options.normalized = 1;
- options.PLOT_KP = 0;
- options.SHOWFIG = 0;
- options.NORMALIZED = 0;
- options.TIME = data.Ts3D;
- % Time window to measure the growth of kx/ky modes
- options.KX_TW = [0.5 1]*data.Ts3D(end);
- options.KY_TW = [0.5 1]*data.Ts3D(end);
- options.NMA = 1; % Set NMA option to 1
- options.NMODES = 999; % Set how much modes we study
- options.iz = 'avg'; % Compressing z
- options.ik = 1; %
- options.GOK2 = 0; % plot gamma/k^2
- options.fftz.flag = 0; % Set fftz.flag option to 0
- options.FIELD = 'phi';
- options.SHOWFIG = 0;
- [~, chi, phib, psib, ~] = plot_ballooning(data,options);
- [fig, wkykx, ekykx] = mode_growth_meter(data,options);
- kx = (1:data.grids.Nx/2)'*2*pi/data.fort_00.GRID.Lx;
- ky = (1:data.grids.Ny/2)'*2*pi/data.fort_00.GRID.Ly;
- [~,~,R,Z] = plot_metric(data,options);
- figure;
- subplot(3,2,1); plot(chi,real(phib),'-b'); hold on;
- plot(chi,imag(phib),'-r'); xlabel('$\chi$'); ylabel('$\phi$')
- subplot(3,2,3); plot(chi,real(psib),'-b'); hold on;
- plot(chi,imag(psib),'-r'); xlabel('$\chi$'); ylabel('$\psi$')
- subplot(3,2,5); errorbar(ky,squeeze(real(wkykx(2:end,1))),...
- squeeze(real(ekykx(2:end,1))),...
- 'ok--','MarkerSize',8,'LineWidth',1.5); hold on;
- errorbar(ky,squeeze(imag(wkykx(2:end,1))),...
- squeeze(imag(ekykx(2:end,1))),...
- '*k--','MarkerSize',8,'LineWidth',1.5);
- xlabel('$k_y\rho_s$'); ylabel('$\gamma (o),\,\omega (*)$');
- R = R*geom.R0; Z = Z*geom.R0;
- subplot(1,2,2); plot(R,Z,'-k'); xlabel('$R$'); ylabel('$Z$');
- axis equal
- title(data.paramshort);
+if (1 && NZ>4)
+%% Ballooning plot
+% [data.PHI, data.Ts3D] = compile_results_3D(DATADIR,J0,J1,'phi');
+if data.inputs.BETA > 0
+[data.PSI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'psi');
end
+options.time_2_plot = [10];
+options.kymodes = 1.5;
+options.normalized = 1;
+options.PLOT_KP = 0;
+% options.field = 'phi';
+options.SHOWFIG = 1;
+[fig, chi, phib, psib, ~] = plot_ballooning(data,options);
+end
+
if 0
%% Hermite-Laguerre spectrum
[data.Napjz, data.Ts3D] = compile_results_3Da(LOCALDIR,J0,J1,'Napjz');
diff --git a/wk/p3_geom_scan_analysis.m b/wk/p3_geom_scan_analysis.m
index d81207835636db7190a19f13bcbfc7d9d95f7581..8136f419b06b03441a60b35055a5fb900c119c48 100644
--- a/wk/p3_geom_scan_analysis.m
+++ b/wk/p3_geom_scan_analysis.m
@@ -14,7 +14,7 @@ GETFLUXSURFACE = 0;
% partition= '../results/paper_3/';
% Get the scan directory
-switch 2
+switch 4
case 1 % delta K_T tau=1
casename = 'DIIID rho95 $\tau=1$';
partition= '../results/paper_3/DIIID_tau_1_rho95_geom_scan/';
@@ -46,6 +46,13 @@ switch 2
nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0;
nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 5; scale2 =500;
t1 = 100; t2 = 150; zfactor = 1;
+ case 4 % HEL CBC
+ casename = 'HEL CBC';
+ partition= '/misc/gyacomo23_outputs/thesis_ch_6/HEL_CBC/128x32x24/';
+ scandir = '.'; scanname= 'CBC HEL';
+ nml1 = 'SPECIES'; pnam1 = '$R_0/L_T\times T_i/T_e$'; attr1 = 'k_T_'; pref1 = 0; scale1 =500;
+ nml2 = 'SPECIES'; pnam2 = '$R_0/L_N$'; attr2 = 'k_N_'; pref2 = 5; scale2 =1.0;
+ t1 = 100; t2 = 150; zfactor = 1;
end
scanname= [casename scanname];
scandir = [partition,scandir,'/'];
@@ -248,7 +255,10 @@ end
subplot(1,2,2)
clrs = jet(N2);
for i = 1:N2
- errorbar(XX(:,i),Zavg(:,i),Zerr(:,i),...
+ % errorbar(XX(:,i),Zavg(:,i),Zerr(:,i),...
+ % 'DisplayName',[pnam2,'=',num2str(p2(i))],...
+ % 'Color',clrs(i,:));
+ plot(XX(:,i),Zavg(:,i),...
'DisplayName',[pnam2,'=',num2str(p2(i))],...
'Color',clrs(i,:));
hold on;
diff --git a/wk/parameters/lin_MTM_Hamed.m b/wk/parameters/lin_MTM_Hamed.m
new file mode 100644
index 0000000000000000000000000000000000000000..113b48fb022155c1a0f0f4e46018b2961174ba05
--- /dev/null
+++ b/wk/parameters/lin_MTM_Hamed.m
@@ -0,0 +1,104 @@
+%% Set simulation parameters
+SIMID = 'lin_MTM_Hamed'; % Name of the simulation
+% These parameters are taken from Frei et al. 2023 MTM illustration
+%% Set up physical parameters
+CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
+NU = 0.05; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+K_Ne = 0; % ele Density '''
+K_Te = 8; % ele Temperature '''
+K_Ni = 0; % ion Density gradient drive
+K_Ti = 0; % ion Temperature '''
+% SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+SIGMA_E = sqrt(0.0027);%0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+NA = 2; % number of kinetic species
+ADIAB_E = (NA==1); % adiabatic electron model
+BETA = 0.0155; % electron plasma beta
+%% Set up grid parameters
+PMAX = 4; % Hermite basis size
+JMAX = 2; % Laguerre basis size
+NX = 24; % real space x-gridpoints
+NY = 2; % real space y-gridpoints
+LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
+LY = 2*pi/0.3; % Size of the squared frequency domain in y direction
+NZ = 64; % number of perpendicular planes (parallel grid)
+SG = 0; % Staggered z grids option
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+
+%% GEOMETRY
+GEOMETRY= 's-alpha';
+% GEOMETRY= 'miller';
+EPS = 0.22; % inverse aspect ratio
+Q0 = 1.3; % safety factor
+SHEAR = 1.1; % magnetic shear
+KAPPA = 1.0; % elongation
+DELTA = 0.0; % triangularity
+ZETA = 0.0; % squareness
+PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
+SHIFT_Y = 0.0; % Shift in the periodic BC in z
+NPOL = 1; % Number of poloidal turns
+
+%% TIME PARAMETERS
+TMAX = 10; % Maximal time unit
+DT = 2e-3; % Time step
+DTSAVE0D = 0.1; % Sampling time for 0D arrays
+DTSAVE2D = -1; % Sampling time for 2D arrays
+DTSAVE3D = 0.1; % Sampling time for 3D arrays
+DTSAVE5D = 100; % Sampling time for 5D arrays
+JOB2LOAD = -1; % Start a new simulation serie
+
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+GKCO = 0; % Gyrokinetic operator
+ABCO = 1; % INTERSPECIES collisions
+INIT_ZF = 0; % Initialize zero-field quantities
+HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+DMAX = -1;
+NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+NMAX = 0;
+KERN = 0; % Kernel model (0 : GK)
+INIT_OPT = 'mom00'; % Start simulation with a noisy mom00/phi/allmom
+NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5)
+
+%% OUTPUTS
+W_DOUBLE = 1; % Output flag for double moments
+W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy)
+W_HF = 1; % Output flag for high-frequency potential energy
+W_PHI = 1; % Output flag for potential
+W_NA00 = 1; % Output flag for nalpha00 (density of species alpha)
+W_DENS = 1; % Output flag for total density
+W_FVEL = 1; % Output flag for fluid velocity
+W_TEMP = 1; % Output flag for temperature
+W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha)
+W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% UNUSED PARAMETERS
+% These parameters are usually not to play with in linear runs
+MU = 0.0; % Hyperdiffusivity coefficient
+MU_X = MU; % Hyperdiffusivity coefficient in x direction
+MU_Y = MU; % Hyperdiffusivity coefficient in y direction
+N_HD = 4; % Degree of spatial-hyperdiffusivity
+N_HDz = 4;
+MU_Z = 2.0; % Hyperdiffusivity coefficient in z direction
+HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model
+MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
+MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
+LAMBDAD = 0.0; % Lambda Debye
+NOISE0 = 0.0e-5; % Initial noise amplitude
+BCKGD0 = 1.0e-5; % Initial background
+S_KAPPA = 0;
+S_DELTA = 0;
+S_ZETA = 0;
+PB_PHASE= 0;
+ADIAB_I = 0;
+MHD_PD = 1;
+EXBRATE = 0;
+K_gB = 1.0; % Magnetic gradient tuner (1 is the real value)
+K_cB = 1.0; % Magnetic curvature tuner (1 is the real value)
+K_mB = 1.0; % mirror force tuner (1 is the real value)
+K_tB = 1.0; % trapping term tuner (1 is the real value)
+K_ldB = 1.0; % Landau damping tuner (1 is the real value)
+COLL_KCUT = 99; % Cutoff for collision operator
diff --git a/wk/triangularitypaper_nonlinear_conv_analysis.m b/wk/triangularitypaper_nonlinear_conv_analysis.m
index cf557d7b2e7676304cd3552902d0932e39598f40..ad145b39b1dcb9663481dd7b4ec24a6a01109586 100644
--- a/wk/triangularitypaper_nonlinear_conv_analysis.m
+++ b/wk/triangularitypaper_nonlinear_conv_analysis.m
@@ -3,19 +3,19 @@ Gx = [...
2.8750,0.8484;...
3.6025,0.8057;...
3.3209,0.3283;...
- 9.4672,12.3946;...
+ 3.1712,0.8761;...
];
Qxe= [...
42.3161,6.4356;...
23.9335,5.4995;...
22.9832,2.3751;...
- 21.2910,3.1200;...
+ 24.4796,8.1372;...
];
Qxi= [...
54.7350,12.8821;...
64.2177,15.8709;...
60.5471,7.2580;...
- 63.4687,10.1195;...
+ 38.5331,10.8445;...
];
@@ -25,6 +25,7 @@ errorbar(PJ,Qxe(:,1),Qxe(:,2)); hold on;
errorbar(PJ,Qxi(:,1),Qxi(:,2));
errorbar(PJ,Gx(:,1),Gx(:,2));
set(gca,'Xscale','log')
+legend({'e heat flux';'i heat flux';'particle flux'})
xlim([5 100])
xticks([6 15 45]);
xticklabels({'(2,1)','(4,2)','(8,4)'});
\ No newline at end of file