Skip to content
Snippets Groups Projects
Commit 207d8e18 authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann
Browse files

FCGK is available now

parent a58f14c5
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
matlab/setup.m
+ 2
− 1
View file @ 207d8e18
...@@ -54,7 +54,8 @@ if (abs(CO) == 2) %Sugama operator ...@@ -54,7 +54,8 @@ if (abs(CO) == 2) %Sugama operator
elseif (abs(CO) == 3) %pitch angle operator elseif (abs(CO) == 3) %pitch angle operator
INITIAL.mat_file = ['''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5''']; INITIAL.mat_file = ['''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'''];
elseif (CO == 4) % Full Coulomb GK elseif (CO == 4) % Full Coulomb GK
INITIAL.mat_file = ['''../../../iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0.h5''']; INITIAL.mat_file = ['''../../../iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_4.h5'''];
% INITIAL.mat_file = ['''../../../iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_k2.h5'''];
elseif (CO == -1) % DGDK elseif (CO == -1) % DGDK
disp('Warning, DGDK not debugged') disp('Warning, DGDK not debugged')
end end
......
wk/linear_study.m
+ 54
− 45
View file @ 207d8e18
for NU = [0.1]
for ETAN = [2.0]
for CO = [1]
%clear all; %clear all;
addpath(genpath('../matlab')) % ... add addpath(genpath('../matlab')) % ... add
default_plots_options default_plots_options
...@@ -19,7 +16,7 @@ LAMBDAD = 0.0; ...@@ -19,7 +16,7 @@ LAMBDAD = 0.0;
NOISE0 = 1.0e-5; NOISE0 = 1.0e-5;
%% GRID PARAMETERS %% GRID PARAMETERS
N = 50; % Frequency gridpoints (Nkx = N/2) N = 50; % Frequency gridpoints (Nkx = N/2)
L = 100; % Size of the squared frequency domain L = 75; % Size of the squared frequency domain
KXEQ0 = 1; % put kx = 0 KXEQ0 = 1; % put kx = 0
MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
...@@ -67,7 +64,8 @@ q0 = 1.0; % safety factor ...@@ -67,7 +64,8 @@ q0 = 1.0; % safety factor
shear = 0.0; % magnetic shear shear = 0.0; % magnetic shear
eps = 0.0; % inverse aspect ratio eps = 0.0; % inverse aspect ratio
%% PARAMETER SCANS %% PARAMETER SCANS
if 1
if 0
%% Parameter scan over PJ %% Parameter scan over PJ
% PA = [2 4 6 10]; % PA = [2 4 6 10];
% JA = [1 2 3 5]; % JA = [1 2 3 5];
...@@ -141,65 +139,79 @@ saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.png']); ...@@ -141,65 +139,79 @@ saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.png']);
end end
end end
if 0 if 1
%% Parameter scan over eta_B %% Parameter scan over CO
PMAXE = 6; PMAXI = 6; PMAXE = 6; PMAXI = 6;
JMAXE = 3; JMAXI = 3; JMAXE = 3; JMAXI = 3;
TMAX = 400; TMAX = 200;
DT = 0.001; DT = 0.01;
eta_B = [0.45, 0.5, 0.55, 0.6, 0.65]; CO_A = [4];
Nparam = numel(eta_B); CONAME_A = {};
param_name = 'etaB'; Nparam = numel(CO_A);
CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty) param_name = 'CO';
NU = 1e-2; % Collision frequency ETAN = 2.0;
NU = 1e-1; % Collision frequency
Bohm_transport = zeros(Nparam,1); Bohm_transport = zeros(Nparam,1);
gamma_Ni = zeros(Nparam,N); gamma_Ni00 = zeros(Nparam,floor(N/2)+1);
gamma_Nipj = zeros(Nparam,floor(N/2)+1);
for i = 1:Nparam for i = 1:Nparam
% Change scan parameter % Change scan parameter
ETAB = eta_B(i); CO = CO_A(i);
setup setup
CONAME_A{i} = CONAME;
% Run linear simulation % Run linear simulation
system(... system(...
['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ./../../../bin/helaz; cd ../../../wk']... ['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz_3 1 6; cd ../../../wk']...
) )
% Load and process results %% Load an process results
filename = ['../results/',SIMID,'/',PARAMS,'/outputs_00.h5'];
load_results load_results
tend = Ts2D(end); tstart = 0.4*tend; tend = Ts3D(end); tstart = 0.4*tend;
[~,itstart] = min(abs(Ts3D-tstart));
[~,itend] = min(abs(Ts3D-tend));
for ikx = 1:N/2+1 for ikx = 1:N/2+1
gamma_Ni(i,ikx) = LinearFit_s(Ts2D,squeeze(abs(Ni00(ikx,1,:))),tstart,tend); gamma_Ni00(i,ikx) = (LinearFit_s(Ts3D(itstart:itend)',(squeeze(abs(Ni00(ikx,1,itstart:itend))))));
Ni00_ST(i,ikx,1:numel(Ts2D)) = squeeze((Ni00(ikx,1,:))); Ni00_ST(i,ikx,1:numel(Ts3D)) = squeeze((Ni00(ikx,1,:)));
end end
gamma_Ni(i,:) = real(gamma_Ni(i,:) .* (gamma_Ni(i,:)>=0.0)); tend = Ts5D(end); tstart = 0.4*tend;
[gmax,ikymax] = max(gamma_Ni(i,:)); [~,itstart] = min(abs(Ts5D-tstart));
kymax = abs(kx(ikymax)); [~,itend] = min(abs(Ts5D-tend));
Bohm_transport(i) = ETAB/ETAN*gmax/kymax^2; for ikx = 1:N/2+1
gamma_Nipj(i,ikx) = LinearFit_s(Ts5D(itstart:itend)',squeeze(max(max(abs(Nipj(:,:,ikx,1,itstart:itend)),[],1),[],2)));
end
gamma_Ni00(i,:) = real(gamma_Ni00(i,:) .* (gamma_Ni00(i,:)>=0.0));
gamma_Nipj(i,:) = real(gamma_Nipj(i,:) .* (gamma_Nipj(i,:)>=0.0));
% kymax = abs(kx(ikymax));
% Bohm_transport(i) = ETAB/ETAN*gmax/kymax^2;
% Clean output % Clean output
system(['rm -r ',BASIC.RESDIR]) system(['rm -r ',BASIC.RESDIR]);
end end
if 0 if 1
%% Plot %% Plot
SCALE = 1;%sqrt(2);
fig = figure; FIGNAME = 'linear_study'; fig = figure; FIGNAME = 'linear_study';
plt = @(x) circshift(x,N/2-1); plt = @(x) x;
for i = 1:Nparam % subplot(211)
clr = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:); for i = 1:Nparam
linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1); clr = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:);
plot(plt(ky),plt(gamma_Ni(i,:)),... linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1);
'Color',clr,... semilogx(plt(SCALE*kx(2:numel(kx))),plt(gamma_Ni00(i,2:end)),...
'LineStyle',linestyle{1},... 'Color',clr,...
'DisplayName',['$\eta_B=$',num2str(eta_B(i))]); 'LineStyle',linestyle{1},'Marker','^',...
hold on; 'DisplayName',[CONAME_A{i}]);
end hold on;
grid on; xlabel('$k_z\rho_s$'); ylabel('$\gamma(N_i^{00})\rho_2/c_s$'); xlim([0.0,max(ky)]); end
title(['$P_e=',num2str(PMAXE),'$',', $J_e=',num2str(JMAXE),'$',... grid on; xlabel('$k_z\rho_s^{R}$'); ylabel('$\gamma(N_i^{00})L_\perp/c_s$'); xlim([0.0,max(kx)]);
', $P_i=',num2str(PMAXE),'$',', $J_i=',num2str(JMAXI),'$']) title(['$\eta=',num2str(ETAB/ETAN),'$, $\nu=',num2str(NU),'$',', $P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$'])
legend('show') legend('show'); xlim([0.01,10])
saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.fig']); saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.fig']);
saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.png']);
end end
if 1 if 0
%% Plot %% Plot
fig = figure; FIGNAME = 'mixing_length'; fig = figure; FIGNAME = 'mixing_length';
plot(eta_B, Bohm_transport) plot(eta_B, Bohm_transport)
...@@ -209,7 +221,4 @@ title(['$P_e=',num2str(PMAXE),'$',', $J_e=',num2str(JMAXE),'$',... ...@@ -209,7 +221,4 @@ title(['$P_e=',num2str(PMAXE),'$',', $J_e=',num2str(JMAXE),'$',...
saveas(fig,[SIMDIR,FIGNAME,'_vs_',param_name,'_',PARAMS,'.fig']); saveas(fig,[SIMDIR,FIGNAME,'_vs_',param_name,'_',PARAMS,'.fig']);
end end
%% %%
end
end
end
end end
\ No newline at end of file
wk/plot_cosol_mat.m
+ 46
− 9
View file @ 207d8e18
...@@ -63,17 +63,18 @@ subplot(224) ...@@ -63,17 +63,18 @@ subplot(224)
title('imag$<$0'); title('imag$<$0');
%% SGGK %% FCGK
P_ = 6; J_ = 3; P_ = 6; J_ = 3;
mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0.h5'; mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0.h5';
FCGK_self = h5read(mat_file_name,'/00000/Caapj/Ciipj'); kp = 1.2; matidx = round(kp/(8/150));
FCGK_ei = h5read(mat_file_name,'/00000/Ceipj/CeipjF')+h5read(mat_file_name,'/00000/Ceipj/CeipjT'); matidx = sprintf('%5.5i',matidx);
FCGK_ie = h5read(mat_file_name,'/00000/Ciepj/CiepjF')+h5read(mat_file_name,'/00000/Ciepj/CiepjT'); FCGK_self = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']);
FCGK_ei = h5read(mat_file_name,['/',matidx,'/Ceipj/CeipjF'])+h5read(mat_file_name,'/00000/Ceipj/CeipjT');
FCGK_ie = h5read(mat_file_name,['/',matidx,'/Ciepj/CiepjF'])+h5read(mat_file_name,'/00000/Ciepj/CiepjT');
figure figure
MAT = 1i*FCGK_self; suptitle('FCGK ii,P=20,J=10, k=0'); MAT = 1i*FCGK_self; suptitle(['FCGK ii,P=20,J=10, k=',num2str(kp)]);
% MAT = 1i*SGGK_ei; suptitle('FCGK ei,P=20,J=10'); % MAT = 1i*FCGK_ei; suptitle('FCGK ei,P=20,J=10');
% MAT = 1i*SGGK_ie; suptitle('FCGK ie,P=20,J=10'); % MAT = 1i*FCGK_ie; suptitle('FCGK ie,P=20,J=10');
subplot(221) subplot(221)
imagesc(abs(MAT)); imagesc(abs(MAT));
title('log abs') title('log abs')
...@@ -86,4 +87,40 @@ subplot(223) ...@@ -86,4 +87,40 @@ subplot(223)
subplot(224) subplot(224)
imagesc(imag(MAT)<0); imagesc(imag(MAT)<0);
title('imag$<$0'); title('imag$<$0');
\ No newline at end of file %% Eigenvalue analysis
if 0
%%
mfns = {'/home/ahoffman/HeLaZ/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_4.h5',...
'/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_k2.h5',...
};
CONAME_A = {'SG','PA','FC NFLR 4', 'FC NFLR k2'};
kp_a = linspace(0,8,149);
gmax = zeros(size(kp_a));
wmax = zeros(size(kp_a));
figure
for j_ = 1:4
mat_file_name = mfns{j_};
i = 1;
for kp = kp_a
matidx = floor(kp/(8/149));
matidx = sprintf('%5.5i',matidx);disp(matidx)
MAT = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']);
gmax(i) = max(real(eig(MAT)));
wmax(i) = max(imag(eig(MAT)));
i = i + 1;
end
subplot(121)
plot(kp_a,gmax,'DisplayName',CONAME_A{j_}); hold on;
subplot(122)
plot(kp_a,wmax,'DisplayName',CONAME_A{j_}); hold on;
end
subplot(121)
legend('show'); grid on;
xlabel('$k_\perp$'); ylabel('$\gamma_{max}$ from Eig(iCa)')
subplot(122)
legend('show'); grid on;
xlabel('$k_\perp$'); ylabel('$\omega_{max}$ from Eig(iCa)')
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment