update scripts

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),')$'])

wk/parameters/lin_RHT.m

+%% 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