From e05fe82004e9e955b040d506a48eba53daf051bd Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Wed, 4 Oct 2023 10:09:19 +0200
Subject: [PATCH] update scripts
---
wk/p3_geom_scan_analysis.m | 74 +++++++++++++++++++--------
wk/parameters/lin_RHT.m | 100 +++++++++++++++++++++++++++++++++++++
2 files changed, 154 insertions(+), 20 deletions(-)
create mode 100644 wk/parameters/lin_RHT.m
diff --git a/wk/p3_geom_scan_analysis.m b/wk/p3_geom_scan_analysis.m
index e3a277aa..38d09ed9 100644
--- a/wk/p3_geom_scan_analysis.m
+++ b/wk/p3_geom_scan_analysis.m
@@ -3,10 +3,10 @@ curdir = pwd;
partition= '/misc/gyacomo23_outputs/paper_3/';
% partition= '../results/paper_3/';
% Get the scan directory
-switch 3
+switch 1
case 1 % delta kappa scan
- % scandir = 'DTT_rho85_geom_scan/P2_J1_delta_kappa_scan/';
- scandir = 'DTT_rho85_geom_scan/P4_J2_delta_kappa_scan/';
+ scandir = 'DTT_rho85_geom_scan/P2_J1_delta_kappa_scan/';
+ % scandir = 'DTT_rho85_geom_scan/P4_J2_delta_kappa_scan/';
nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0.23;
nml2 = 'GEOMETRY'; pnam2 = '$\kappa$'; attr2 = 'kappa'; pref2 = 1.53;
case 2 % shear safety factor scan
@@ -15,14 +15,14 @@ switch 3
nml1 = 'GEOMETRY'; pnam1 = '$\hat s$'; attr1 = 'shear'; pref1 = 3.63;
nml2 = 'GEOMETRY'; pnam2 = '$q_0$'; attr2 = 'q0'; pref2 =-2.15;
case 3
- scandir = 'DTT_rho85_geom_scan/P2_J1_delta_nuDGGK_scan/';
+ % scandir = 'DTT_rho85_geom_scan/P2_J1_delta_nuDGGK_scan/';
% scandir = 'DTT_rho85_geom_scan/P4_J2_delta_nuDGGK_scan/';
% scandir = 'DTT_rho85_geom_scan/P2_J1_delta_nuSGGK_scan/';
% scandir = 'DTT_rho85_geom_scan/P4_J2_delta_nuSGGK_scan/';
% scandir = 'DTT_rho85_geom_scan/P2_J1_delta_nuLDGK_scan/';
- % scandir = 'DTT_rho85_geom_scan/P4_J2_delta_nuLDGK_scan/';
+ scandir = 'DTT_rho85_geom_scan/P4_J2_delta_nuLDGK_scan/';
nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0.23;
- nml2 = 'MODEL'; pnam2 = '$\nu$'; attr2 = 'nu'; pref2 = 0.5;
+ nml2 = 'MODEL'; pnam2 = '$\nu$'; attr2 = 'nu'; pref2 = 999;%0.5;
end
scandir = [partition,scandir];
% Get a list of all items in the current directory
@@ -30,6 +30,8 @@ contents = dir(scandir);
% give ref value and param names
REFVAL= 1;
+% normalize to the max
+NORMAL= 0;
% Get and plot the fluxsurface
GETFLUXSURFACE = 0;
@@ -48,7 +50,7 @@ for i = 1:length(contents)
para2 = [para2 param.(nml2).(attr2)];
% Now you are in the subdirectory. You can perform operations here.
[t_all, Pxi_all, Qxi_all, Pxe_all, Qxe_all] = read_flux_out_XX(subdir);
- if(t_all(end) > 200)
+ if(t_all(end) > 100)
[fullAvg,sliceAvg,sliceErr] = sliceAverage(Qxi_all+Qxe_all,3);
Qxavg = [Qxavg fullAvg];
Qxerr = [Qxerr mean(sliceErr)];
@@ -101,14 +103,30 @@ YY = YY(idx1,idx2);
% compute the
if REFVAL
- [~,iref1] = min(abs(XX(:,1)-pref1));
- [~,iref2] = min(abs(YY(1,:)-pref2));
+ if NORMAL
+ Qxname = '$\bar Q_{tot}/\bar Q_{max}[\%]$';
+ [tmp,iref1] = max(Zavg);
+ [~, iref2] = max(tmp);
+ iref1 = iref1(iref2);
+ else
+ Qxname = '$\langle (Q_{tot}-Q_{ref})/Q_{ref} \rangle_t[\%]$';
+ if pref1 ~= 999
+ [~,iref1] = min(abs(XX(:,1)-pref1));
+ else
+ iref1 = 1:N1;
+ end
+ if pref2 ~= 999
+ [~,iref2] = min(abs(YY(1,:)-pref2));
+ else
+ iref2 = 1:N2;
+ end
+ end
+ iref1 = ones(N1,1).*iref1;
+ iref2 = ones(N2,1).*iref2;
xref = XX(iref1,iref2);
yref = YY(iref1,iref2);
Qxref = Zavg(iref1,iref2);
- Zavg = (Zavg-Qxref)/Qxref * 100;
- Qxname = '$\langle (Q_{tot}-Q_{ref})/Q_{ref} \rangle_t[\%]$';
- Qrefname = ['$Q_{ref}=$',num2str(Qxref)];
+ % Qrefname = ['$Q_{ref}=$',num2str(Qxref)];
else
Qxname = '$\langle Q_{tot} \rangle_t$';
end
@@ -118,22 +136,38 @@ figure
subplot(1,2,1)
% contourf(XX(:,1),YY(1,:),Zavg',13); hold on
[xx_,yy_] = meshgrid(XX(:,1),YY(1,:));
-imagesc_custom(xx_,yy_,Zavg'); hold on
if REFVAL
- plot(xref,yref,'xk','MarkerSize',14,'DisplayName',Qrefname)
+ if NORMAL
+ imagesc_custom(xx_,yy_,(Zavg./Qxref)'); hold on
+ CLIM = [0 1];
+ else
+ imagesc_custom(xx_,yy_,((Zavg-Qxref)./Qxref * 100)'); hold on
+ CLIM = 'auto';
+ end
+else
+ imagesc_custom(xx_,yy_,Zavg'); hold on
+ CLIM = 'auto';
end
+% if REFVAL
+% plot(xref,yref,'xk','MarkerSize',14,'DisplayName',Qrefname)
+% end
xlabel(pnam1); ylabel(pnam2);
title(Qxname)
-colormap(bluewhitered); colorbar;
+colormap(bluewhitered); colorbar; clim(CLIM);
+if ~REFVAL
+ colormap(hot);
+end
subplot(1,2,2)
for i = 1:N2
errorbar(XX(:,i),Zavg(:,i),Zerr(:,i),...
'DisplayName',[pnam2,'=',num2str(para2(1,i))]);
hold on;
end
-if REFVAL
- plot(xref,0,'xk','MarkerSize',14,'DisplayName',Qrefname)
-end
+% if REFVAL
+% plot(xref,0,'xk','MarkerSize',14,'DisplayName',Qrefname)
+% end
grid on
-xlabel(pnam1); ylabel(Qxname);
-legend('show');
+xlabel(pnam1); ylabel('$\langle Q_{tot} \rangle_t$');
+legend('show','Location','northwest');
+title([param.COLLISION.collision_model{1}, ...
+ ', $(P,J)=(',num2str(param.GRID.pmax),',',num2str(param.GRID.jmax),')$'])
diff --git a/wk/parameters/lin_RHT.m b/wk/parameters/lin_RHT.m
new file mode 100644
index 00000000..72137f75
--- /dev/null
+++ b/wk/parameters/lin_RHT.m
@@ -0,0 +1,100 @@
+%% Set simulation parameters
+SIMID = 'lin_ITG'; % Name of the simulation
+
+%% Set up physical parameters
+CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
+NU = 0.005; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+K_Ne = 0; % ele Density
+K_Te = 0; % 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)
+NA = 1; % number of kinetic species
+ADIAB_E = (NA==1); % adiabatic electron model
+BETA = 0.0; % electron plasma beta
+EXBRATE = 0.0; % Background ExB shear flow
+%% Set up grid parameters
+P = 4;
+J = 2;%P/2;
+PMAX = P; % Hermite basis size
+JMAX = J; % Laguerre basis size
+NX = 2; % real space x-gridpoints
+NY = 2; % real space y-gridpoints
+LX = 2*pi/0.05; % Size of the squared frequency domain in x direction
+LY = 2*pi/0.01; % Size of the squared frequency domain in y direction
+NZ = 8; % 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.1; % 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'; % 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 = 50; % Maximal time unit
+DT = 1e-3; % Time step
+DTSAVE0D = 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 = 'RK3'; % 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_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
+MU_Z = 0.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.0; % Initial background
+k_gB = 1.0; % Magnetic gradient strength
+k_cB = 1.0; % Magnetic curvature strength
+COLL_KCUT = 1; % Cutoff for collision operator
+S_KAPPA = 0.0;
+S_DELTA = 0.0;
+S_ZETA = 0.0;
+PB_PHASE= 0;
+ADIAB_I = 0;
+MHD_PD = 0;
\ No newline at end of file
--
GitLab