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

use of the new routine growth rate

parent 2c36b823
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
matlab/compute/mode_growth_meter.m
+ 54
− 8
View file @ 6038dc39
......@@ -30,6 +30,8 @@ TW = [OPTIONS.KY_TW; OPTIONS.KX_TW];
[~,ikzf] = max(mean(squeeze(abs(field(1,:,FRAMES))),2));
[wkykx, ekykx] = compute_growth_rates(DATA.PHI(:,:,:,FRAMES),DATA.Ts3D(FRAMES));
FIGURE.fig = figure; %set(gcf, 'Position', [100 100 1200 700])
FIGURE.FIGNAME = 'mode_growth_meter';
for i = 1:2
......@@ -41,14 +43,20 @@ for i = 1:2
if MODES_SELECTOR == 2
switch OPTIONS.ik
case 'sum'
plt_ = @(x,ik) movmean(squeeze(sum(abs((x(:,ik,FRAMES))),1)),Nma);
plt_ = @(f,ik) movmean(squeeze(sum(abs((f(:,ik,FRAMES))),1)),Nma);
MODESTR = '$\sum_{k_y}$';
omega= @(ik) wkykx(1,:,end);
err = @(ik) ekykx(1,:);
case 'max'
plt_ = @(x,ik) movmean(squeeze(max(abs((x(:,ik,FRAMES))),[],1)),Nma);
plt_ = @(f,ik) movmean(squeeze(max(abs((f(:,ik,FRAMES))),[],1)),Nma);
MODESTR = '$\max_{k_y}$';
omega= @(ik) wkykx(1,:,end);
err = @(ik) ekykx(1,:);
otherwise
plt_ = @(x,ik) movmean(abs(squeeze(x(OPTIONS.ik,ik,FRAMES))),Nma);
plt_ = @(f,ik) movmean(abs(squeeze(f(OPTIONS.ik,ik,FRAMES))),Nma);
MODESTR = ['$k_y=$',num2str(DATA.grids.ky(OPTIONS.ik))];
omega= @(ik) wkykx(OPTIONS.ik,:,end);
err = @(ik) ekykx(OPTIONS.ik,:);
end
kstr = 'k_x';
% Number max of modes to plot is kx>0 (1/2) of the non filtered modes (2/3)
......@@ -59,13 +67,19 @@ for i = 1:2
case 'sum'
plt_ = @(x,ik) movmean(squeeze(sum(abs((x(ik,:,FRAMES))),2)),Nma);
MODESTR = '$\sum_{k_x}$';
omega= @(ik) wkykx(:,1,end);
err = @(ik) ekykx(:,1);
case 'max'
plt_ = @(x,ik) movmean(squeeze(max(abs((x(ik,:,FRAMES))),[],2)),Nma);
MODESTR = '$\max_{k_x}$';
omega= @(ik) wkykx(:,1,end);
err = @(ik) ekykx(:,1);
otherwise
plt_ = @(x,ik) movmean(abs(squeeze(x(ik,OPTIONS.ik,FRAMES))),Nma);
MODESTR = ['$k_x=$',num2str(DATA.grids.kx(OPTIONS.ik))];
end
omega= @(ik) wkykx(:,OPTIONS.ik,end);
err = @(ik) ekykx(:,OPTIONS.ik);
end
kstr = 'k_y';
% Number max of modes to plot is ky>0 (1/1) of the non filtered modes (2/3)
Nmax = ceil(DATA.grids.Nky*2/3);
......@@ -81,6 +95,8 @@ for i = 1:2
gamma = MODES;
amp = MODES;
% w_ky = zeros(numel(MODES),numel(FRAMES)-1);
% ce = zeros(numel(MODES),numel(FRAMES));
i_=1;
for ik = MODES
......@@ -88,9 +104,29 @@ for i = 1:2
gr = polyfit(t(it1:it2),to_measure(it1:it2),1);
gamma(i_) = gr(1);
amp(i_) = gr(2);
% % Second method for measuring growth rate (measures frequ. too)
% if MODES_SELECTOR == 1
% to_measure = reshape(DATA.PHI(ik,1,:,FRAMES),DATA.grids.Nz,numel(FRAMES));
% else
% to_measure = reshape(DATA.PHI(1,ik,:,FRAMES),DATA.grids.Nz,numel(FRAMES));
% end
% for it = 2:numel(FRAMES)
% phi_n = to_measure(:,it);
% phi_nm1 = to_measure(:,it-1);
% dt = t(it)-t(it-1);
% ZS = sum(phi_nm1,1); %sum over z
%
% wl = log(phi_n./phi_nm1)/dt;
% w_ky(i_,it) = squeeze(sum(wl.*phi_nm1,1)./ZS);
%
% % for iky = 1:numel(w_ky(:,it))
% % ce(iky,it) = abs(sum(abs(w_ky(iky,it)-wl(iky,:)).^2.*phi_nm1(iky,:),2)./ZS(iky,:));
% % end
% end
i_=i_+1;
end
mod2plot = [2:min(Nmax,OPTIONS.NMODES+1)];
mod2plot = 2:min(Nmax,OPTIONS.NMODES+1);
clr_ = jet(numel(mod2plot));
%plot
% subplot(2+d,3,1+3*(i-1))
......@@ -120,15 +156,25 @@ for i = 1:2
% subplot(2+d,3,3+3*(i-1))
FIGURE.axes(3+3*(i-1)) = subplot(2+d,3,3+3*(i-1),'parent',FIGURE.fig);
plot(k(MODES),gamma,...
'DisplayName',['(',num2str(DATA.inputs.PMAX),',',num2str(DATA.inputs.JMAX),')']); hold on;
% yyaxis("left")
errorbar(k(MODES),real(omega(ik)),real(err(ik)),'-k',...
'DisplayName',...
['$\gamma$, (',num2str(DATA.inputs.PMAX),',',num2str(DATA.inputs.JMAX),')']);
hold on;
for i_ = 1:numel(mod2plot)
plot(k(MODES(mod2plot(i_))),gamma(mod2plot(i_)),'x','color',clr_(i_,:));
end
if MODES_SELECTOR == 2
plot(k(ikzf),gamma(ikzf),'ok');
end
xlabel(['$',kstr,'\rho_s$']); ylabel('$\gamma$');
% ylabel('$\gamma$');
% yyaxis("right")
errorbar(k(MODES),imag(omega(ik)),imag(err(ik)),'--k',...
'DisplayName',...
['$\omega$, (',num2str(DATA.inputs.PMAX),',',num2str(DATA.inputs.JMAX),')']);
hold on;
ylabel('$\gamma,\omega$');
xlabel(['$',kstr,'\rho_s$']);
title('Growth rates')
end
......
wk/lin_run_script.m
+ 13
− 10
View file @ 6038dc39
......@@ -6,8 +6,8 @@
%% Set up the paths for the necessary Matlab modules
gyacomodir = pwd;
gyacomodir = gyacomodir(1:end-2);
% mpirun = 'mpirun';
mpirun = '/opt/homebrew/bin/mpirun'; % for macos
mpirun = 'mpirun';
% mpirun = '/opt/homebrew/bin/mpirun'; % for macos
addpath(genpath([gyacomodir,'matlab'])) % Add matlab folder
addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot folder
addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute folder
......@@ -24,20 +24,23 @@ EXECNAME = 'gyacomo23_sp'; % single precision
% run lin_DTT_HM_rho85
% run lin_DTT_HM_rho98
% run lin_DTT_LM_rho90
run lin_DTT_LM_rho95
% run lin_JET_rho97
% run lin_DTT_LM_rho95
run lin_JET_rho97
% run lin_Entropy
% run lin_ITG
% run lin_KBM
%% Change parameters
NY = 2;
PMAX = 4;
NX = 4;
PMAX = 2;
JMAX = PMAX/2;
ky = 1.0;
ky = 0.5;
LY = 2*pi/ky;
DT = 1e-5;
% SIGMA_E = 0.04;
TMAX = 10;
DTSAVE0D = 200*DT;
DT = 1e-3;
TMAX = 10;
% % SIGMA_E = 0.04;
% TMAX = 10;
% DTSAVE0D = 200*DT;
DTSAVE3D = TMAX/50;
%%-------------------------------------------------------------------------
%% RUN
......
wk/lin_scan_script.m
+ 26
− 32
View file @ 6038dc39
......@@ -16,37 +16,41 @@ addpath(genpath([gyacomodir,'wk/parameters'])) % Add parameters folder
%% Setup run or load an executable
RUN = 1; % To run or just to load
RERUN = 1; % rerun if the data does not exist
RERUN = 0; % rerun if the data does not exist
default_plots_options
EXECNAME = 'gyacomo23_sp'; % single precision
% EXECNAME = 'gyacomo23_dp'; % double precision
%% Setup basic parameters
run lin_DTT_AB_rho85_PT
%% Setup parameters
% run lin_DTT_AB_rho85
% run lin_DTT_AB_rho98
run lin_JET_rho97
% run lin_Entropy
% run lin_ITG
%% Modify parameters
% NZ = 1;
NY = 2;
DT = 0.5e-3;
TMAX = 50;
MU_X = 0.1; MU_Y = 0.1;
%% Change parameters
NY = 2;
% SIGMA_E = 0.023;
%% Scan parameters
SIMID = [SIMID,'_scan'];
P_a = [2 4 6 8];
ky_a = 0.05:0.05:1.5;
ky_a = logspace(-1.5,1.5,30);
%% Scan loop
% arrays for the result
g_ky = zeros(numel(ky_a),numel(P_a),2);
g_avg= g_ky*0;
g_ky = zeros(numel(ky_a),numel(P_a));
g_std= g_ky*0;
w_ky = g_ky*0;
w_std= g_ky*0;
j = 1;
for PMAX = P_a
JMAX = P/2;
i = 1;
for ky = ky_a
LY = 2*pi/ky;
LY = 2*pi/ky;
DT = 1e-4;%min(1e-2,1e-3/ky);
TMAX = 10;%min(10,1.5/ky);
DTSAVE0D = 0.1;
DTSAVE3D = 0.1;
%% RUN
setup
% naming
......@@ -67,8 +71,7 @@ for PMAX = P_a
end
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
[data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
if numel(data_.Ts3D)>10
if numel(data_.Ts3D)>5
if numel(data_.Ts0D)>10
% Load results after trying to run
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
......@@ -83,23 +86,14 @@ for PMAX = P_a
[~,it1] = min(abs(data_.Ts3D-0.5*data_.Ts3D(end))); % start of the measurement time window
[~,it2] = min(abs(data_.Ts3D-1.0*data_.Ts3D(end))); % end of ...
field = 0;
field_t = 0;
for ik = 2:NY/2+1
field = squeeze(sum(abs(data_.PHI),3)); % take the sum over z
field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only
to_measure = log(field_t(it1:it2));
tw = double(data_.Ts3D(it1:it2));
% gr = polyfit(tw,to_measure,1);
gr = fit(tw,to_measure,'poly1');
err= confint(gr);
g_ky(i,j,ik) = gr.p1;
g_std(i,j,ik) = abs(err(2,1)-err(1,1))/2;
end
[wkykx,ekykx] = compute_growth_rates(data_.PHI(:,:,:,it1:it2),data_.Ts3D(it1:it2));
g_ky (i,j) = real(wkykx(2,1,end));
g_std(i,j) = real(ekykx(2,1));
w_ky (i,j) = imag(wkykx(2,1,end));
w_std(i,j) = imag(ekykx(2,1));
[gmax, ikmax] = max(g_ky(i,j,:));
msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
end
end
i = i + 1;
end
......@@ -107,11 +101,11 @@ for PMAX = P_a
end
%% take max growth rate among z coordinate
y_ = g_ky(:,:,2);
e_ = g_std(:,:,2);
y_ = g_ky + 1i*w_ky;
e_ = g_std+ 1i*w_std;
%% Save scan results (gamma)
if(numel(ky_a)>1 && numel(P_a)>1)
if(numel(ky_a)>1 || numel(P_a)>1)
pmin = num2str(min(P_a)); pmax = num2str(max(P_a));
kymin = num2str(min(ky_a)); kymax= num2str(max(ky_a));
filename = [num2str(NX),'x',num2str(NZ),...
......
wk/load_metadata_scan.m
+ 59
− 36
View file @ 6038dc39
......@@ -13,26 +13,32 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0'
% datafname = 'lin_DTT_AB_rho85_PT_scan/16x24_ky_0.1_0.5_P_2_8_DGGK_0.05_be_0.0034.mat';
% datafname = 'lin_DTT_AB_rho85_PT_scan/16x24_ky_0.1_0.5_P_2_8_DGGK_0.05_be_0.0039.mat';
% datafname = 'lin_DTT_AB_rho85_PT_scan/16x24_ky_0.1_1_P_2_8_DGGK_0.05_be_0.0039.mat';
datafname = 'lin_DTT_AB_rho85_PT_scan/16x24_ky_0.1_0.9_P_2_8_kN_1.33_DGGK_0.56901_be_0.0039.mat';
% datafname = 'lin_DTT_AB_rho85_PT_scan/16x24_ky_0.1_0.9_P_2_8_kN_1.33_DGGK_0.56901_be_0.0039.mat';
% datafname = 'lin_DTT_AB_rho85_PT_scan/8x24_ky_0.05_1.5_P_2_8_kN_1.33_DGGK_0.1_be_0.0039.mat';
% datafname = 'lin_JET_rho97_scan/4x32_ky_0.01_10_P_2_8_kN_10_DGGK_0.1_be_0.0031.mat';
% datafname = 'lin_JET_rho97_scan/4x32_ky_0.01_10_P_2_8_kN_10_DGGK_0.05_be_0.0031.mat';
% datafname = 'lin_JET_rho97_scan/8x32_ky_0.01_10_P_2_8_kN_10_DGGK_0.05_be_0.0031.mat';
% datafname = 'lin_JET_rho97_scan/8x32_ky_0.01_10_P_2_2_kN_10_DGGK_0.05_be_0.0031.mat';
datafname = 'lin_JET_rho97_scan/8x32_ky_0.031623_31.6228_P_2_2_kN_10_DGGK_0.1_be_0.0031.mat';
%% Chose if we filter gamma>0.05
FILTERGAMMA = 0;
%% Load data
fname = ['../results/',datafname];
d = load(fname);
gamma = real(d.data); g_err = real(d.err);
omega = imag(d.data); w_err = imag(d.err);
if FILTERGAMMA
d.data = d.data.*(d.data>0.025);
d.err = d.err.*(d.data>0.025);
d.data = real(gamma).*(real(gamma)>0.025) + imag(gamma);
end
if 0
%% Pcolor of the peak
figure;
% [XX_,YY_] = meshgrid(d.s1,d.s2);
[XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2));
pclr=imagesc_custom(XX_,YY_,d.data'.*(d.data>0)');
% pclr=contourf(1:numel(d.s1),1:numel(d.s2),d.data'.*(d.data>0)');
% pclr=surf(1:numel(d.s1),1:numel(d.s2),d.data'.*(d.data>0)');
pclr=imagesc_custom(XX_,YY_,gamma'.*(gamma>0)');
% pclr=contourf(1:numel(d.s1),1:numel(d.s2),gamma'.*(gamma>0)');
% pclr=surf(1:numel(d.s1),1:numel(d.s2),gamma'.*(gamma>0)');
title(d.title);
xlabel(d.s1name); ylabel(d.s2name);
set(gca,'XTick',1:numel(d.s1),'XTicklabel',d.s1)
......@@ -48,23 +54,40 @@ if 1
%% Scan along first dimension
figure
colors_ = jet(numel(d.s2));
subplot(121)
for i = 1:numel(d.s2)
% plot(d.s1,d.data(:,i),'s-',...
plot(d.s1(d.data(:,i)>0),d.data((d.data(:,i)>0),i),'s-',...
'LineWidth',2.0,...
'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
'color',colors_(i,:));
% errorbar(d.s1,d.data(:,i),d.err(:,i),'s-',...
% 'LineWidth',2.0,...
% 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
% 'color',colors_(i,:));
% plot(d.s1,gamma(:,i),'s-',...
% plot(d.s1(gamma(:,i)>0),gamma((gamma(:,i)>0),i),'s-',...
% 'LineWidth',2.0,...
% 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
% 'color',colors_(i,:));
errorbar(d.s1,gamma(:,i),g_err(:,i),'s-',...
'LineWidth',2.0,...
'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
'color',colors_(i,:));
hold on;
end
xlabel(d.s1name); ylabel(d.dname);title(d.title);
xlim([d.s1(1) d.s1(end)]);
subplot(122)
for i = 1:numel(d.s2)
% plot(d.s1,gamma(:,i),'s-',...
% plot(d.s1(gamma(:,i)>0),gamma((gamma(:,i)>0),i),'s-',...
% 'LineWidth',2.0,...
% 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
% 'color',colors_(i,:));
errorbar(d.s1,omega(:,i),w_err(:,i),'s-',...
'LineWidth',2.0,...
'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
'color',colors_(i,:));
hold on;
end
xlabel(d.s1name); ylabel('$\omega c_s/R$');title(d.title);
xlim([d.s1(1) d.s1(end)]);
colormap(colors_);
clb = colorbar;
% caxis([d.s2(1)-0.5,d.s2(end)+0.5]);
clim([1 numel(d.s2)+1]);
clb.Ticks=linspace(d.s2(1),d.s2(end),numel(d.s2));
clb.Ticks =1.5:numel(d.s2)+1.5;
......@@ -78,11 +101,11 @@ if 0
figure
colors_ = jet(numel(d.s1));
for i = 1:numel(d.s1)
plot(d.s2,d.data(i,:),'s-',...
plot(d.s2,gamma(i,:),'s-',...
'LineWidth',2.0,...
'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
'color',colors_(i,:));
% errorbar(d.s2,d.data(i,:),d.err(i,:),'s-',...
% errorbar(d.s2,gamma(i,:),g_err(i,:),'s-',...
% 'LineWidth',2.0,...
% 'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
% 'color',colors_(i,:));
......@@ -104,13 +127,13 @@ end
if 0
%% Convergence analysis
figure
% target_ = 0.25*(d.data(end,end)+d.data(end-i_,end)+d.data(end,end-i_)+d.data(end-i_,end-i_));
% target_ = 0.25*(gamma(end,end)+gamma(end-i_,end)+gamma(end,end-i_)+gamma(end-i_,end-i_));
colors_ = jet(numel(d.s1));
for i = 1:numel(d.s1)
% target_ = d.data(i,end);
% target_ = gamma(i,end);
target_ = 2.79666916212537142172e-01; % Value for nuDGDK = 0.05, kT=6.96, (40,20), Nkx=8
if target_ > 0
eps_ = abs(target_ - d.data(i,1:end-1))/abs(target_);
eps_ = abs(target_ - gamma(i,1:end-1))/abs(target_);
semilogy(d.s2(1:end-1),eps_,'-s',...
'LineWidth',2.0,...
'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
......@@ -140,30 +163,30 @@ figure; i_ = 0;
% target_ = 1.39427e-01; % Value for nuDGDK = 0.001, kT=5.3, (30,16), Nkx=8
target_ = 2.72510405826983714839e-01 % Value for nuDGDK = 0.001, kT=6.96, (50,25), Nkx=8
% target_ = 2.73048910051283844069e-01; % Value for nuDGDK = 0.0, kT=6.96, (40,20), Nkx=8
% target_ = 0.25*(d.data(end,end)+d.data(end-i_,end)+d.data(end,end-i_)+d.data(end-i_,end-i_));
% eps_ = log(abs(target_ - d.data)/abs(target_));
eps_ = max(-10,log(abs(target_ - d.data)/abs(target_)));
sign_ = 1;%sign(d.data - target_);
eps_ = d.data;
% target_ = 0.25*(gamma(end,end)+gamma(end-i_,end)+gamma(end,end-i_)+gamma(end-i_,end-i_));
% eps_ = log(abs(target_ - gamma)/abs(target_));
eps_ = max(-10,log(abs(target_ - gamma)/abs(target_)));
sign_ = 1;%sign(gamma - target_);
eps_ = gamma;
for i = 1:numel(d.s1)
for j = 1:numel(d.s2)
% target_ = d.data(i,end);
% target_ = d.data(i,end);
% target_ = d.data(end,j);
eps_(i,j) = log(abs(target_ - d.data(i,j))/target);
% target_ = gamma(i,end);
% target_ = gamma(i,end);
% target_ = gamma(end,j);
eps_(i,j) = log(abs(target_ - gamma(i,j))/target);
if target_ > 0
% eps_(i,:) = max(-12,log(abs(target_ - d.data(i,1:end))/abs(target_)));
% eps_(i,:) = log(abs(target_ - d.data(i,1:end))/abs(target_));
% eps_(i,:) = min(100,100*abs(target_ - d.data(i,1:end))/abs(target_));
% eps_(i,:) = max(-12,log(abs(target_ - gamma(i,1:end))/abs(target_)));
% eps_(i,:) = log(abs(target_ - gamma(i,1:end))/abs(target_));
% eps_(i,:) = min(100,100*abs(target_ - gamma(i,1:end))/abs(target_));
else
end
end
end
[XX_,YY_] = meshgrid(d.s1,d.s2);
[XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2));
pclr=imagesc_custom(XX_,YY_,eps_'.*(d.data>0)'.*sign_');
% pclr=contourf(1:numel(d.s1),1:numel(d.s2),eps_'.*(d.data>0)'.*sign_',5);
% pclr=surf(1:numel(d.s1),1:numel(d.s2),eps_'.*(d.data>0)'.*sign_');
pclr=imagesc_custom(XX_,YY_,eps_'.*(gamma>0)'.*sign_');
% pclr=contourf(1:numel(d.s1),1:numel(d.s2),eps_'.*(gamma>0)'.*sign_',5);
% pclr=surf(1:numel(d.s1),1:numel(d.s2),eps_'.*(gamma>0)'.*sign_');
title(d.title);
xlabel(d.s1name); ylabel(d.s2name);
set(gca,'XTick',1:numel(d.s1),'XTicklabel',d.s1)
......
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