From 80d6e7cb0cccd384de0f1a0b4a2d2a5d7eeb58b2 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Fri, 10 Feb 2023 11:34:37 +0100
Subject: [PATCH] save the scripts
---
fort_example.90 | 4 +-
matlab/compile_results.m | 10 +-
matlab/extract_fig_data.m | 2 +-
matlab/load/load_3D_data.m | 14 +-
matlab/load/load_gene_data.m | 2 +
matlab/load/quick_gene_load.m | 4 +-
matlab/plot/plot_metric.m | 6 +-
.../plot_radial_transport_and_spacetime.m | 10 +-
matlab/plot/show_moments_spectrum.m | 168 ++++++++++++------
matlab/setup.m | 4 +-
matlab/write_fort90.m | 4 +-
testcases/cyclone_example/fort.90 | Bin 1378 -> 1376 bytes
testcases/fort.90 | 4 +-
testcases/matlab_testscripts/Hallenbert.m | 4 +-
.../matlab_testscripts/cyclone_test_case.m | 4 +-
.../linear_1D_entropy_mode.m | 4 +-
testcases/matlab_testscripts/linear_damping.m | 4 +-
testcases/miller_example_w_salpha/fort_00.90 | 4 +-
testcases/miller_example_w_salpha/fort_01.90 | 4 +-
testcases/zpinch_example/fort.90 | 4 +-
wk/Ajay_scan_CH4_lin_ITG.m | 4 +-
wk/CBC_kT_PJ_scan.m | 2 +-
wk/CBC_kT_scan_salpha.m | 4 +-
wk/CBC_ky_PJ_scan.m | 2 +-
wk/CBC_nu_CO_scan.m | 4 +-
wk/CBC_nu_CO_scan_miller.m | 4 +-
wk/CBC_nu_CO_scan_salpha.m | 4 +-
wk/Dimits_fig3.m | 2 +-
wk/analysis_gene.m | 15 +-
wk/analysis_gyacomo.m | 50 +++---
wk/header_3D_results.m | 37 +++-
wk/lin_3D_Zpinch.m | 4 +-
wk/lin_3Dzpinch.m | 4 +-
wk/lin_ETPY.m | 4 +-
wk/lin_ITG.m | 56 +++---
wk/lin_KBM.m | 4 +-
wk/lin_MTM.m | 4 +-
wk/lin_RHT.m | 4 +-
wk/lin_TEM.m | 4 +-
wk/local_run.m | 4 +-
wk/marconi_run.m | 4 +-
wk/p2_heatflux_salpha_convergence.m | 67 +++++--
wk/quick_run.m | 4 +-
43 files changed, 333 insertions(+), 218 deletions(-)
diff --git a/fort_example.90 b/fort_example.90
index 69c7414d..7607e1e1 100644
--- a/fort_example.90
+++ b/fort_example.90
@@ -68,8 +68,8 @@
K_Te = 6.96 !electron temperature gradient intensity
K_Ni = 2.22 !ion density gradient intensity
K_Ti = 6.96 !ion density temperature intensity
- GradB = 1 !magnetic field gradient strength
- CurvB = 1 !magnetic curvature strength
+ k_gB = 1 !magnetic field gradient strength
+ k_cB = 1 !magnetic curvature strength
lambdaD = 0 !Debye length (not tested when non zero)
/
&COLLISION_PAR
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index fccc7ee8..7d6962c6 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -160,13 +160,15 @@ while(CONTINUE)
Nepjz = load_3D_data(filename, 'Nepjz');
Nepjz_ = cat(4,Nepjz_,Nepjz); clear Nepjz
catch
- end
+ disp('Cannot load Nepjz');
+ end
end
-% try
+ try
[Nipjz, Ts3D, ~] = load_3D_data(filename, 'Nipjz');
Nipjz_ = cat(4,Nipjz_,Nipjz); clear Nipjz
-% catch
-% end
+ catch
+ disp('Cannot load Nipjz');
+ end
end
if W_DENS
if KIN_E
diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index 2e19518e..8c8f0d92 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 = [450 650];
+tw = [100 180];
fig = gcf;
axObjs = fig.Children;
diff --git a/matlab/load/load_3D_data.m b/matlab/load/load_3D_data.m
index 9c9df71d..2d87dfb1 100644
--- a/matlab/load/load_3D_data.m
+++ b/matlab/load/load_3D_data.m
@@ -6,14 +6,24 @@ function [ data, time, dt ] = load_3D_data( filename, variablename )
% Find array size by loading the first output
tmp = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+1,'%06d')]);
- sz = size(tmp.real);
+ try % check if it is complex or real
+ sz = size(tmp.real);
+ cmpx = 1;
+ catch
+ sz = size(tmp);
+ cmpx = 0;
+ end
% add time dimension
sz = [sz numel(time)];
data = zeros(sz);
for it = 1:numel(time)
tmp = h5read(filename,['/data/var3d/',variablename,'/', num2str(cstart+it,'%06d')]);
- data(:,:,:,it) = tmp.real + 1i * tmp.imaginary;
+ if cmpx
+ data(:,:,:,it) = tmp.real + 1i * tmp.imaginary;
+ else
+ data(:,:,:,it) = tmp;
+ end
end
end
diff --git a/matlab/load/load_gene_data.m b/matlab/load/load_gene_data.m
index 2fdff459..8c877d2e 100644
--- a/matlab/load/load_gene_data.m
+++ b/matlab/load/load_gene_data.m
@@ -22,6 +22,8 @@ DATA.Ny = DATA.Nky*2-1;
DATA.z = h5read([folder,coofile],'/coord/z');
DATA.Nz = numel(DATA.z);
+DATA.paramshort = [num2str(DATA.Nkx),'x',num2str(DATA.Nky),'x',num2str(DATA.Nz),...
+ 'x',num2str(DATA.Nvp),'x',num2str(DATA.Nmu)];
if numel(DATA.kx)>1
dkx = DATA.kx(2);
else
diff --git a/matlab/load/quick_gene_load.m b/matlab/load/quick_gene_load.m
index 4d1fc2bd..3edaeae6 100644
--- a/matlab/load/quick_gene_load.m
+++ b/matlab/load/quick_gene_load.m
@@ -73,12 +73,12 @@ path = '/home/ahoffman/gene/linear_CBC_results/';
%----------Convergence nvpar shearless CBC
% fname = 'CBC_salpha_nz_24_nv_scan_nw_16_adiabe.txt';
% fname = 'CBC_miller_nz_24_nv_scan_nw_16_adiabe.txt';
-fname = 'CBC_salpha_nz_24_nv_scan_nw_16_kine.txt';
+% fname = 'CBC_salpha_nz_24_nv_scan_nw_16_kine.txt';
% fname = 'CBC_miller_nz_24_nv_scan_nw_16_kine.txt';
%---------- CBC
% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_kine.txt';
-% fname = 'CBC_miller_nx_20_nz_32_nv_32_nw_12_adiabe.txt';
+fname = 'CBC_miller_nx_20_nz_32_nv_32_nw_12_adiabe.txt';
% fname = 'CBC_miller_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
% fname = 'CBC_miller_nx_20_nz_32_nv_32_nw_12_kine.txt';
% fname = 'CBC_miller_nx_8_nz_24_nv_36_nw_16_kine.txt';
diff --git a/matlab/plot/plot_metric.m b/matlab/plot/plot_metric.m
index 0eb8c096..f03fa6be 100644
--- a/matlab/plot/plot_metric.m
+++ b/matlab/plot/plot_metric.m
@@ -1,9 +1,9 @@
-function [ fig ] = plot_metric( data, options )
+function [ fig, arrays ] = plot_metric( data, options )
% names = {'Jacobian','gradxB','gradyB','gradzB','gradz_coeff',...
% 'gxx','gxy','gxz','gyy','gyz','gzz','hatB','hatR','hatZ'};
names = {'gxx','gxy','gxz','gyy','gyz','gzz',...
- 'hatB', 'gradxB', 'gradyB', 'gradzB',...
+ 'hatB', 'dBdx', 'dBdy', 'dBdz',...
'Jacobian','hatR','hatZ','gradz_coeff'};
geo_arrays = zeros(2,data.Nz,numel(names));
@@ -40,5 +40,7 @@ if NPLOT > 0
axis equal
end
end
+%outputs
+arrays = squeeze(geo_arrays(1,:,:));
end
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index 884a919e..098eae03 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -10,8 +10,6 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS,CODE)
Gx_infty_std = std (DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
Qx_infty_avg = mean(DATA.HFLUX_X(its0D:ite0D))*SCALE;
Qx_infty_std = std (DATA.HFLUX_X(its0D:ite0D))*SCALE;
- disp(['G_x=',sprintf('%2.2e',Gx_infty_avg),'+-',sprintf('%2.2e',Gx_infty_std)]);
- disp(['Q_x=',sprintf('%2.2e',Qx_infty_avg),'+-',sprintf('%2.2e',Qx_infty_std)]);
% disp(['Q_x=',sprintf('%2.2e',Qx_infty_avg),'+-',sprintf('%2.2e',Qx_infty_std)]);
f_avg_z = squeeze(mean(DATA.PHI(:,:,:,:),3));
[~,ikzf] = max(squeeze(mean(abs(f_avg_z(1,:,its3D:ite3D)),3)));
@@ -28,7 +26,7 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS,CODE)
end
Qx_avg = mean(Qx_ee);
Qx_err = std(Qx_ee);
-% disp(['Q_avg=',sprintf('%2.2e',Qx_avg),'+-',sprintf('%2.2e',Qx_err)]);
+ disp(['Q_avg=',sprintf('%2.2e',Qx_avg),'+-',sprintf('%2.2e',Qx_err)]);
%% computations
% Compute zonal and non zonal energies
@@ -66,14 +64,14 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; %set(gcf, 'Position', [500, 1000, 1000, 600])
FIGURE.ax1 = subplot(3,1,1,'parent',FIGURE.fig);
plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI*SCALE),'--',...
- 'color',clr_((DATA.Pmaxi-1)/2-1,:),...
+ 'color',clr_((DATA.Pmaxi-1)/2,:),...
'DisplayName',['$\Gamma_x$ ',DATA.paramshort]); hold on;
plot(mvm(DATA.Ts0D),mvm(DATA.HFLUX_X*SCALE),'-',...
- 'color',clr_((DATA.Pmaxi-1)/2-1,:),...
+ 'color',clr_((DATA.Pmaxi-1)/2,:),...
'DisplayName',['$Q_x$ ',DATA.paramshort]); hold on;
ylabel('Transport')
if(~isnan(Qx_infty_avg))
- plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Qx_infty_avg, '-k',...
+ plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Qx_avg, '-k',...
'DisplayName',['$Q_{avg}=',sprintf('%2.2f',Qx_avg),'\pm',sprintf('%2.2f',Qx_err),'$']); legend('show');
ylim([0,5*abs(Qx_infty_avg)]);
else
diff --git a/matlab/plot/show_moments_spectrum.m b/matlab/plot/show_moments_spectrum.m
index 3f000729..d3839297 100644
--- a/matlab/plot/show_moments_spectrum.m
+++ b/matlab/plot/show_moments_spectrum.m
@@ -2,12 +2,13 @@ function [ FIGURE ] = show_moments_spectrum( DATA, OPTIONS )
Pi = DATA.Pi;
Ji = DATA.Ji;
-if (isempty(DATA.Nipjz))
+if ~(isempty(DATA.Nipjz))
Time_ = DATA.Ts3D;
Nipj = sum(abs(DATA.Nipjz),3);
else
Time_ = DATA.Ts5D;
- Nipj = sum(sum(sum(abs(DATA.Nipj),3),4),5);
+% Nipj = sum(sum(sum(abs(DATA.Nipj).^2,3),4),5);
+ Nipj = sum(sum(sum(conj(DATA.Nipj).*DATA.Nipj,3),4),5);
end
Nipj = squeeze(Nipj);
@@ -89,75 +90,128 @@ if ~OPTIONS.ST
title([TITLE,' He-La elec. spectrum']);
end
else
- [JJ,PP] = meshgrid(Ji,Pi);
- P2Ji = PP + 2*JJ;
- % form an axis of velocity ordered moments
- p2ji = unique(P2Ji);
- % weights to average
- weights = zeros(size(p2ji));
- % space time of moments amplitude wrt to p+2j degree
- Ni_ST = zeros(numel(p2ji),numel(Time_));
- % fill the st diagramm by averaging every p+2j=n moments together
- for ip = 1:numel(Pi)
- for ij = 1:numel(Ji)
- [~,ip2j] = min(abs(p2ji-(Pi(ip)+2*Ji(ij))));
- Ni_ST(ip2j,:) = Ni_ST(ip2j,:) + transpose(squeeze(Nipj(ip,ij,:)));
- weights(ip2j) = weights(ip2j) + 1;
+ if OPTIONS.P2J
+ plotname = '$\langle\sum_k |N_i^{pj}|\rangle_{p+2j=const}$';
+ [JJ,PP] = meshgrid(Ji,Pi);
+ P2Ji = PP + 2*JJ;
+ % form an axis of velocity ordered moments
+ y_ = unique(P2Ji); yname = '$p+2j$';
+ % weights to average
+ weights = zeros(size(y_));
+ % space time of moments amplitude wrt to p+2j degree
+ Ni_ST = zeros(numel(y_),numel(Time_));
+ % fill the st diagramm by averaging every p+2j=n moments together
+ for ip = 1:numel(Pi)
+ for ij = 1:numel(Ji)
+ [~,ip2j] = min(abs(y_-(Pi(ip)+2*Ji(ij))));
+ Ni_ST(ip2j,:) = Ni_ST(ip2j,:) + transpose(squeeze(Nipj(ip,ij,:)));
+ weights(ip2j) = weights(ip2j) + 1;
+ end
end
- end
- % doing the average
- for ip2j = 1:numel(p2ji)
- Ni_ST(ip2j,:) = Ni_ST(ip2j,:)/weights(ip2j);
- end
- if DATA.KIN_E
- % same for electrons!!
- [JJ,PP] = meshgrid(Je,Pe);
- P2Je = PP + 2*JJ;
- % form an axis of velocity ordered moments
- p2je = unique(P2Je);
- % weights to average
- weights = zeros(size(p2je));
- % space time of moments amplitude wrt to p+2j degree
- Ne_ST = zeros(numel(p2je),numel(Time_));
- % fill the st diagramm by averaging every p+2j=n moments together
- for ip = 1:numel(Pe)
- for ij = 1:numel(Je)
- [~,ip2j] = min(abs(p2ji-(Pe(ip)+2*Je(ij))));
- Ne_ST(ip2j,:) = Ne_ST(ip2j,:) + transpose(squeeze(Nepj(ip,ij,:)));
- weights(ip2j) = weights(ip2j) + 1;
+ % doing the average
+ for ip2j = 1:numel(y_)
+ Ni_ST(ip2j,:) = Ni_ST(ip2j,:)/weights(ip2j);
end
- end
- % doing the average
- for ip2j = 1:numel(p2ji)
- Ne_ST(ip2j,:) = Ne_ST(ip2j,:)/weights(ip2j);
- end
+ if DATA.KIN_E
+ % same for electrons!!
+ [JJ,PP] = meshgrid(Je,Pe);
+ P2Je = PP + 2*JJ;
+ % form an axis of velocity ordered moments
+ p2je = unique(P2Je);
+ % weights to average
+ weights = zeros(size(p2je));
+ % space time of moments amplitude wrt to p+2j degree
+ Ne_ST = zeros(numel(p2je),numel(Time_));
+ % fill the st diagramm by averaging every p+2j=n moments together
+ for ip = 1:numel(Pe)
+ for ij = 1:numel(Je)
+ [~,ip2j] = min(abs(y_-(Pe(ip)+2*Je(ij))));
+ Ne_ST(ip2j,:) = Ne_ST(ip2j,:) + transpose(squeeze(Nepj(ip,ij,:)));
+ weights(ip2j) = weights(ip2j) + 1;
+ end
+ end
+ % doing the average
+ for ip2j = 1:numel(y_)
+ Ne_ST(ip2j,:) = Ne_ST(ip2j,:)/weights(ip2j);
+ end
+ end
+ ticks_labels = {};
+ else % We just order our moments w.r.t. to the convention ..
+ % (0,0) (1,0) (2,0) (0,1) (3,0) (1,1) (4,0) (2,1) (0,2) etc.
+ plotname = '$\langle\sum_k |N_i^{pj}|^2\rangle_z$';
+ Nmoments = numel(Nipj(:,:,1)); % total number of moments
+ HL_deg = zeros(Nmoments,2); % list of degrees, first column Hermite, second Laguerre
+ im = 2;
+ deg = 1; % the zero degree is always here first place
+ ticks_labels = cell(10,1);
+ ticks_labels{1} = '(0,0)';
+ while(im<=Nmoments)
+ FOUND = 1;
+ while(FOUND) % As we find a pair matching the degree we retry
+ FOUND = 0;
+ for ij = 1:DATA.Jmaxi
+ for ip = 1:DATA.Pmaxi
+ if((ip-1)+2*(ij-1) == deg)
+ % Check if the pair is already added
+ check_ = HL_deg == [DATA.Pi(ip) DATA.Pi(ij)];
+ check_ = sum(check_(:,1) .* check_(:,2));
+ if ~check_ % if not add it
+ HL_deg(im,1) = DATA.Pi(ip);
+ HL_deg(im,2) = DATA.Pi(ij);
+ ticks_labels{im} = ['(',num2str(DATA.Pi(ip)),',',num2str(DATA.Ji(ij)),')'];
+ im = im + 1; FOUND = 1;
+ end
+ end
+ end
+ end
+ end
+ % No pair found anymore, increase the degree
+ deg = deg + 1;
+ end
+
+ % form an axis of velocity ordered moments
+ y_ = 1:Nmoments; yname = '$(P,J)$';
+ % space time of moments amplitude wrt to p+2j degree
+ Ni_ST = zeros(Nmoments,numel(Time_));
+ for i_ = 1:Nmoments
+ Ni_ST(i_,:) = Nipj(HL_deg(i_,1)+1,HL_deg(i_,2)+1,:);
+ end
+ if DATA.KIN_E
+ % space time of moments amplitude wrt to p+2j degree
+ Ne_ST = zeros(Nmoments,numel(Time_));
+ for i_ = 1:Nmoments
+ Ne_ST = Nepj(HL_deg(i_,1)+1,HL_deg(i_,2)+1,:);
+ end
+ end
+
end
% plots
% scaling
-% plt = @(x,ip2j) x./max(max(x));
if OPTIONS.NORMALIZED
- plt = @(x,ip2j) x(ip2j,:)./max(x(ip2j,:));
+ plt = @(x,i) x(i,:)./max(x(i,:));
else
- plt = @(x,ip2j) x;
+ plt = @(x,i) x;
end
if DATA.KIN_E
subplot(2,1,1)
end
- imagesc(Time_,p2ji,plt(Ni_ST,1:numel(p2ji)));
- set(gca,'YDir','normal')
-% pclr = pcolor(XX,YY,plt(Ni_ST));
-% set(pclr, 'edgecolor','none'); hold on;
- xlabel('$t$'); ylabel('$p+2j$')
- title('$\langle\sum_k |N_i^{pj}|\rangle_{p+2j=const}$')
+
+ imagesc(Time_,y_,plt(Ni_ST,1:numel(y_)));
+ set(gca,'YDir','normal')
+ xlabel('$t$'); ylabel(yname);
+ if ~isempty(ticks_labels)
+ yticks(y_);
+ yticklabels(ticks_labels)
+ end
+ title(plotname)
+
if DATA.KIN_E
subplot(2,1,2)
- imagesc(Time_,p2je,plt(Ne_ST,1:numel(p2ji)));
+ imagesc(Time_,p2je,plt(Ne_ST,1:numel(y_)));
set(gca,'YDir','normal')
-% pclr = pcolor(XX,YY,plt(Ne_ST));
-% set(pclr, 'edgecolor','none'); hold on;
- xlabel('$t$'); ylabel('$p+2j$')
- title('$\langle\sum_k |N_e^{pj}|\rangle_{p+2j=const}$')
- suptitle(DATA.param_title);
+ xlabel('$t$'); ylabel(yname)
+ title(plotname)
+ suptitle(DATA.param_title);
end
end
diff --git a/matlab/setup.m b/matlab/setup.m
index 8f398119..2c66ec48 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -52,8 +52,8 @@ MODEL.K_Ni = K_Ni;
MODEL.K_Ne = K_Ne;
MODEL.K_Ti = K_Ti;
MODEL.K_Te = K_Te;
-MODEL.GradB = GRADB; % Magnetic gradient
-MODEL.CurvB = CURVB; % Magnetic curvature
+MODEL.k_gB = k_gB; % Magnetic gradient
+MODEL.k_cB = k_cB; % Magnetic curvature
MODEL.lambdaD = LAMBDAD;
% Collision parameters
COLL.collision_model = ['''',CO,''''];
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index e708ca0f..92c9d1f4 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -77,8 +77,8 @@ fprintf(fid,[' K_Ne = ', num2str(MODEL.K_Ne),'\n']);
fprintf(fid,[' K_Te = ', num2str(MODEL.K_Te),'\n']);
fprintf(fid,[' K_Ni = ', num2str(MODEL.K_Ni),'\n']);
fprintf(fid,[' K_Ti = ', num2str(MODEL.K_Ti),'\n']);
-fprintf(fid,[' GradB = ', num2str(MODEL.GradB),'\n']);
-fprintf(fid,[' CurvB = ', num2str(MODEL.CurvB),'\n']);
+fprintf(fid,[' k_gB = ', num2str(MODEL.k_gB),'\n']);
+fprintf(fid,[' k_cB = ', num2str(MODEL.k_cB),'\n']);
fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),'\n']);
fprintf(fid,[' beta = ', num2str(MODEL.beta),'\n']);
fprintf(fid,'/\n');
diff --git a/testcases/cyclone_example/fort.90 b/testcases/cyclone_example/fort.90
index d562465a47ad3fcff4858ddad9ea05cd9b707fec..e595b41a0a5c39dee1064a074fc77283af24b09c 100644
GIT binary patch
delta 29
kcmaFF^?+*wCyPLKe7ci@f`YArA(w(ec6{<?X%=@z0Dz1KmH+?%
delta 31
mcmaFB^@wW&CyS7KQDTadf`Wprf+3fJf^%t6*=9)=cSZn{E(lft
diff --git a/testcases/fort.90 b/testcases/fort.90
index abce02cf..fa0b6dea 100644
--- a/testcases/fort.90
+++ b/testcases/fort.90
@@ -63,8 +63,8 @@
K_Ti = 6.92
K_Ne = 0
K_Te = 0
- GradB = 1
- CurvB = 1
+ k_gB = 1
+ k_cB = 1
lambdaD = 0
beta = 0
/
diff --git a/testcases/matlab_testscripts/Hallenbert.m b/testcases/matlab_testscripts/Hallenbert.m
index 768341d0..5ded4b11 100644
--- a/testcases/matlab_testscripts/Hallenbert.m
+++ b/testcases/matlab_testscripts/Hallenbert.m
@@ -24,8 +24,8 @@ J = 2;
Q0 = 1.0; % safety factor
SHEAR = 0.0; % magnetic shear
EPS = 0.0; % inverse aspect ratio
-GRADB = 1.0; % Magnetic gradient
-CURVB = 1.0; % Magnetic curvature
+k_gB = 1.0; % Magnetic gradient
+k_cB = 1.0; % Magnetic curvature
SG = 0; % Staggered z grids option
%% TIME PARAMETERS
TMAX = 2000; % Maximal time unit
diff --git a/testcases/matlab_testscripts/cyclone_test_case.m b/testcases/matlab_testscripts/cyclone_test_case.m
index 0a6b241b..f1c72abc 100644
--- a/testcases/matlab_testscripts/cyclone_test_case.m
+++ b/testcases/matlab_testscripts/cyclone_test_case.m
@@ -23,8 +23,8 @@ J = 2;
Q0 = 1.4; % safety factor
SHEAR = 0.0; % magnetic shear
EPS = 0.18; % inverse aspect ratio
-GRADB = 1.0; % Magnetic gradient
-CURVB = 1.0; % Magnetic curvature
+k_gB = 1.0; % Magnetic gradient
+k_cB = 1.0; % Magnetic curvature
SG = 1; % Staggered z grids option
%% TIME PARAMETERS
TMAX = 500; % Maximal time unit
diff --git a/testcases/matlab_testscripts/linear_1D_entropy_mode.m b/testcases/matlab_testscripts/linear_1D_entropy_mode.m
index e6008157..02de8d1a 100644
--- a/testcases/matlab_testscripts/linear_1D_entropy_mode.m
+++ b/testcases/matlab_testscripts/linear_1D_entropy_mode.m
@@ -75,8 +75,8 @@ MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
-GRADB = 1.0;
-CURVB = 1.0;
+k_gB = 1.0;
+k_cB = 1.0;
%% PARAMETER SCANS
if 1
diff --git a/testcases/matlab_testscripts/linear_damping.m b/testcases/matlab_testscripts/linear_damping.m
index 394e93b5..46b6c4e8 100644
--- a/testcases/matlab_testscripts/linear_damping.m
+++ b/testcases/matlab_testscripts/linear_damping.m
@@ -67,8 +67,8 @@ MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 0*1.0e-5; % Init noise amplitude
BCKGD0 = 1.0; % Init background
-GRADB = 0.0;
-CURVB = 0.0;
+k_gB = 0.0;
+k_cB = 0.0;
%% PARAMETER SCANS
if 1
diff --git a/testcases/miller_example_w_salpha/fort_00.90 b/testcases/miller_example_w_salpha/fort_00.90
index 2701de9d..45bb97e7 100644
--- a/testcases/miller_example_w_salpha/fort_00.90
+++ b/testcases/miller_example_w_salpha/fort_00.90
@@ -63,8 +63,8 @@
K_Ti = 6.92
K_Ne = 0
K_Te = 0
- GradB = 1
- CurvB = 1
+ k_gB = 1
+ k_cB = 1
lambdaD = 0
beta = 0
/
diff --git a/testcases/miller_example_w_salpha/fort_01.90 b/testcases/miller_example_w_salpha/fort_01.90
index 7deb16c3..965375d1 100644
--- a/testcases/miller_example_w_salpha/fort_01.90
+++ b/testcases/miller_example_w_salpha/fort_01.90
@@ -63,8 +63,8 @@
K_Ti = 6.92
K_Ne = 0
K_Te = 0
- GradB = 1
- CurvB = 1
+ k_gB = 1
+ k_cB = 1
lambdaD = 0
beta = 0
/
diff --git a/testcases/zpinch_example/fort.90 b/testcases/zpinch_example/fort.90
index 53a1a749..d11d88e8 100644
--- a/testcases/zpinch_example/fort.90
+++ b/testcases/zpinch_example/fort.90
@@ -62,8 +62,8 @@
K_Te = 0.4
K_Ni = 2.0
K_Ti = 0.4
- GradB = 1
- CurvB = 1
+ k_gB = 1
+ k_cB = 1
lambdaD = 0
beta = 0
/
diff --git a/wk/Ajay_scan_CH4_lin_ITG.m b/wk/Ajay_scan_CH4_lin_ITG.m
index 1c1d2e5a..5836570d 100644
--- a/wk/Ajay_scan_CH4_lin_ITG.m
+++ b/wk/Ajay_scan_CH4_lin_ITG.m
@@ -115,8 +115,8 @@ for NU = NU_a
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
- GRADB = 1.0;
- CURVB = 1.0;
+ k_gB = 1.0;
+ k_cB = 1.0;
%% RUN
setup
% naming
diff --git a/wk/CBC_kT_PJ_scan.m b/wk/CBC_kT_PJ_scan.m
index 22581920..7ff5ff63 100644
--- a/wk/CBC_kT_PJ_scan.m
+++ b/wk/CBC_kT_PJ_scan.m
@@ -77,7 +77,7 @@ P = P_a(i); J = J_a(i);
MU_J = 0.0; LAMBDAD = 0.0;
NOISE0 = 1.0e-4; % Init noise amplitude
BCKGD0 = 0.0; % Init background
- GRADB = 1.0;CURVB = 1.0;
+ k_gB = 1.0;k_cB = 1.0;
%%-------------------------------------------------------------------------
% RUN
diff --git a/wk/CBC_kT_scan_salpha.m b/wk/CBC_kT_scan_salpha.m
index db4aab4b..6229c864 100644
--- a/wk/CBC_kT_scan_salpha.m
+++ b/wk/CBC_kT_scan_salpha.m
@@ -116,8 +116,8 @@ for KT = KT_a
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
- GRADB = 1.0;
- CURVB = 1.0;
+ k_gB = 1.0;
+ k_cB = 1.0;
%% RUN
setup
% naming
diff --git a/wk/CBC_ky_PJ_scan.m b/wk/CBC_ky_PJ_scan.m
index d9bc9e05..b53dde9d 100644
--- a/wk/CBC_ky_PJ_scan.m
+++ b/wk/CBC_ky_PJ_scan.m
@@ -69,7 +69,7 @@ for j = 1
MU_J = 0.0; LAMBDAD = 0.0;
NOISE0 = 0.0e-5; % Init noise amplitude
BCKGD0 = 1.0; % Init background
-GRADB = 1.0;CURVB = 1.0;
+k_gB = 1.0;k_cB = 1.0;
end
%%-------------------------------------------------------------------------
% RUN
diff --git a/wk/CBC_nu_CO_scan.m b/wk/CBC_nu_CO_scan.m
index 019a372c..437e0065 100644
--- a/wk/CBC_nu_CO_scan.m
+++ b/wk/CBC_nu_CO_scan.m
@@ -115,8 +115,8 @@ for NU = NU_a
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
- GRADB = 1.0;
- CURVB = 1.0;
+ k_gB = 1.0;
+ k_cB = 1.0;
%% RUN
setup
% naming
diff --git a/wk/CBC_nu_CO_scan_miller.m b/wk/CBC_nu_CO_scan_miller.m
index fcece647..e9a6384d 100644
--- a/wk/CBC_nu_CO_scan_miller.m
+++ b/wk/CBC_nu_CO_scan_miller.m
@@ -118,8 +118,8 @@ for NU = NU_a
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
- GRADB = 1.0;
- CURVB = 1.0;
+ k_gB = 1.0;
+ k_cB = 1.0;
%% RUN
setup
% naming
diff --git a/wk/CBC_nu_CO_scan_salpha.m b/wk/CBC_nu_CO_scan_salpha.m
index 62830349..c4acfa75 100644
--- a/wk/CBC_nu_CO_scan_salpha.m
+++ b/wk/CBC_nu_CO_scan_salpha.m
@@ -117,8 +117,8 @@ for NU = NU_a
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
- GRADB = 1.0;
- CURVB = 1.0;
+ k_gB = 1.0;
+ k_cB = 1.0;
%% RUN
setup
% naming
diff --git a/wk/Dimits_fig3.m b/wk/Dimits_fig3.m
index 21ad776c..3571b46c 100644
--- a/wk/Dimits_fig3.m
+++ b/wk/Dimits_fig3.m
@@ -1,4 +1,4 @@
-%% Heat flux Qi [R/rhos^2/cs]
+%% Heat flux Qi
kN = 2.22;
%-------------- GM ---------------
%(P,J)=(2,1)
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 11c0b16a..09e9901f 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -47,7 +47,7 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add
%Paper 2
% folder = '/misc/gene_results/CBC/KT_6.96_64x32x32x24x12_Nexc_5/';
% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x8x4_Nexc_5_00/';
-% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x16x8_Nexc_5_00/';
+folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x16x8_Nexc_5_00/';
% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x32x16_Nexc_5_00/';
% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x32x16_Nexc_5_01/';
@@ -55,9 +55,12 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add
% folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x32x16_Nexc_5_01/';
% folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x16x8_Nexc_5/';
% folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x8x4_Nexc_5/';
-folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x8x4_Nexc_5_smallvbox/';
+% folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x8x4_Nexc_5_smallvbox/';
% folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x16x8_Nexc_5_largexbox/';
+% debug ? shearless
+% folder = '/misc/gene_results/CBC/shearless_CBC_128x64x24x24x12_00/';
+% folder = '/misc/gene_results/CBC/shearless_CBC_128x64x24x24x12_01/';
gene_data = load_gene_data(folder);
gene_data.FIGDIR = folder;
gene_data = invert_kxky_to_kykx_gene_results(gene_data);
@@ -95,10 +98,10 @@ options.INTERP = 0;
options.POLARPLOT = 0;
options.AXISEQUAL = 0;
options.NORMALIZE = 0;
-% options.NAME = '\phi';
+options.NAME = '\phi';
% options.NAME = '\psi';
% options.NAME = '\omega_z';
-options.NAME = 'n_i';
+% options.NAME = 'n_i';
% options.NAME = 'n_i-n_e';
% options.NAME = '\phi^{NZ}';
% options.NAME = 'N_i^{00}';
@@ -107,8 +110,8 @@ options.NAME = 'n_i';
% options.NAME = 'Q_x';
% options.NAME = 'k^2n_e';
options.PLAN = 'xy';
-options.COMP = 1;
-options.TIME = [0];
+options.COMP = 'avg';
+options.TIME = [50 200 500];
options.RESOLUTION = 256;
data.a = data.EPS * 2e3;
diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m
index 000f5022..b05f7cb7 100644
--- a/wk/analysis_gyacomo.m
+++ b/wk/analysis_gyacomo.m
@@ -30,12 +30,14 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
% data.scale = 1;%/(data.Nx*data.Ny)^2;
-i_ = 1;
-disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))])
-disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))])
-options.TAVG_0 = data.TJOB_SE(i_)+600;%0.4*data.Ts3D(end);
-options.TAVG_1 = data.TJOB_SE(i_+1);%0.9*data.Ts3D(end); % Averaging times duration
-options.NCUT = 4; % Number of cuts for averaging and error estimation
+% jid_ = 0;
+% disp([num2str(data.TJOB_SE(2*jid_+1)),' ',num2str(data.TJOB_SE(2*(jid_+1)))])
+% disp([num2str(data.NU_EVOL(2*jid_+1)),' ',num2str(data.NU_EVOL(2*(jid_+1)))])
+% options.TAVG_0 = data.TJOB_SE(2*jid_+1);%0.4*data.Ts3D(end);
+% options.TAVG_1 = data.TJOB_SE(2*(jid_+1));%0.9*data.Ts3D(end); % Averaging times duration
+options.TAVG_0 = 100;
+options.TAVG_1 = 1000;
+options.NCUT = 5; % Number of cuts for averaging and error estimation
options.NMVA = 1; % Moving average for time traces
% options.ST_FIELD = '\Gamma_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
@@ -67,23 +69,23 @@ end
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
-options.INTERP = 1;
+options.INTERP = 0;
options.POLARPLOT = 0;
-options.NAME = '\phi';
+% options.NAME = '\phi';
% options.NAME = '\omega_z';
% options.NAME = 'N_i^{00}';
% options.NAME = 's_{Ey}';
% options.NAME = 'n_i^{NZ}';
-% options.NAME = 'G_x';
-% options.NAME = 'n_i';
+% options.NAME = 'Q_x';
+options.NAME = 'n_i';
% options.NAME = 'n_i-n_e';
-options.PLAN = 'xy';
+options.PLAN = 'yz';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
% options.TIME = data.Ts5D(end-30:end);
-options.TIME = data.Ts3D;
-% options.TIME = [0:10000];
+% options.TIME = data.Ts3D;
+options.TIME = [800:1000];
data.EPS = 0.1;
data.a = data.EPS * 2000;
options.RESOLUTION = 256;
@@ -100,7 +102,7 @@ options.NORMALIZE = 0;
options.NAME = '\phi';
% options.NAME = '\psi';
% options.NAME = '\omega_z';
-% options.NAME = 'n_i';
+% options.NAME = 'T_i';
% options.NAME = 'n_i-n_e';
% options.NAME = '\phi^{NZ}';
% options.NAME = 'N_i^{00}';
@@ -109,8 +111,8 @@ options.NAME = '\phi';
% options.NAME = 'Q_x';
% options.NAME = 'k^2n_e';
options.PLAN = 'xy';
-options.COMP = 1;
-options.TIME = [29.5 30 30.5];
+options.COMP = 'avg';
+options.TIME = [50 200 500 1000];
options.RESOLUTION = 256;
data.a = data.EPS * 2e3;
@@ -140,7 +142,7 @@ options.XPERP = linspace( 0,sqrt(6),16).^2;
% options.SPAR = gene_data.vp';
% options.XPERP = gene_data.mu';
options.iz = 'avg';
-options.T = [0.5:0.1:1]*data.Ts3D(end);
+options.T = [0.9:0.1:1]*data.Ts3D(end);
% options.PLT_FCT = 'contour';
% options.PLT_FCT = 'contourf';
options.PLT_FCT = 'surfvv';
@@ -156,13 +158,9 @@ if 0
%% Hermite-Laguerre spectrum
% options.TIME = 'avg';
options.P2J = 0;
-options.ST = 0;
-options.PLOT_TYPE = 'space-time';
+options.ST = 1;
options.NORMALIZED = 0;
-options.JOBNUM = 0;
-options.TIME = [200:500];
-options.specie = 'i';
-options.compz = 'avg';
+options.TIME = [180:10000];
fig = show_moments_spectrum(data,options);
% fig = show_napjz(data,options);
% save_figure(data,fig,'.png');
@@ -170,7 +168,7 @@ end
if 0
%% Time averaged spectrum
-options.TIME = [5000 9000];
+options.TIME = [180 9000];
options.NORM =1;
% options.NAME = '\phi';
% options.NAME = 'N_i^{00}';
@@ -206,7 +204,7 @@ if 0
%% Mode evolution
options.NORMALIZED = 1;
options.TIME = [000:9000];
-options.KX_TW = [25 55]; %kx Growth rate time window
+options.KX_TW = [1 20]; %kx Growth rate time window
options.KY_TW = [0 20]; %ky Growth rate time window
options.NMA = 1;
options.NMODES = 800;
@@ -229,7 +227,7 @@ if 0
%% Metric infos
options.SHOW_FLUXSURF = 0;
options.SHOW_METRICS = 1;
-fig = plot_metric(data,options);
+[fig, geo_arrays] = plot_metric(data,options);
end
if 0
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 9d8b87fe..9a554ec4 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -33,7 +33,7 @@ PARTITION = '/misc/gyacomo_outputs/';
% resdir = 'CBC/old/128x64x16x5x3';
% resdir = 'CBC/96x96x16x3x2_Nexc_6';
% resdir = 'CBC/128x96x16x3x2_Nexc_0';
-% resdir = 'CBC/old/192x96x24x13x7';
+% resdir = 'CBC/192x96x24x13x7';
% resdir = 'CBC/128x96x16x3x2_Nexc_0_periodic_chi';
% resdir = 'CBC/64x32x16x3x2_Nexc_0_periodic_chi';
@@ -77,13 +77,20 @@ PARTITION = '/misc/gyacomo_outputs/';
%% CBC Miller
% resdir = 'GCM_CBC/daint/Miller_GX_comparison';
% resdir = 'GCM_CBC/daint/Salpha_GX_comparison';
-%% Paper 2 simulations
+% resdir = 'Mandell_benchmark/5x3';
+% resdir = 'Mandell_benchmark/new_5x3';
+% resdir = 'Mandell_benchmark/new_5x3_more_diff';
+% resdir = 'Mandell_benchmark/7x5';
+% resdir = 'Mandell_benchmark/new_7x5';
+% resdir = 'Mandell_benchmark/new_7x5_more_diff';
+% % Paper 2 simulations
% convergence CBC and Dimits regime
% resdir = 'paper_2_nonlinear/kT_6.96/CBC_7x4x128x64x24';
% resdir = 'paper_2_nonlinear/kT_6.96/CBC_3x2x128x64x24_Nexc_5';
% resdir = 'paper_2_nonlinear/kT_6.96/CBC_5x3x128x64x24_Nexc_5';
% resdir = 'paper_2_nonlinear/kT_6.96/CBC_7x4x128x64x24_Nexc_5';
% resdir = 'paper_2_nonlinear/kT_6.96/CBC_9x5x128x64x24_Nexc_5';
+% resdir = 'paper_2_nonlinear/kT_6.96/CBC_21x11x128x64x24_Nexc_5';
% resdir = 'paper_2_nonlinear/kT_5.3/CBC_11x6x128x64x24_Nexc_5';
% resdir = 'paper_2_nonlinear/kT_5.3/CBC_kT_5.3_3x2x128x64x24_Nexc_5';
% resdir = 'paper_2_nonlinear/kT_5.3/CBC_kT_5.3_5x3x128x64x24_Nexc_5';
@@ -95,12 +102,28 @@ PARTITION = '/misc/gyacomo_outputs/';
% resdir = 'paper_2_nonlinear/kT_5.3/CBC_kT_5.3_17x9x128x64x24_Nexc_5';
% Scan in kT
% resdir = 'paper_2_nonlinear/kT_scan_DGGK_0.05/9x5x128x64x24';
-
% resdir = 'paper_2_nonlinear/CBC_rerun/rerun_CBC_3x2x128x64x16';
-% resdir = 'dev/init_ppj';
-% resdir = 'dev/hatB_NL';
-resdir = 'dev/CBC_wave_study';
+% resdir = 'paper_2_nonlinear/CBC_rerun/rerun_CBC_5x3x128x64x16';
+% resdir = 'dev/Napjz_spectrum';
+% resdir = 'dev/CBC_test';
+
+% Reruns
+% resdir = 'paper_2_nonlinear/kT_6.96_rerun/5x3x192x96x24';
+% resdir = 'paper_2_nonlinear/kT_6.96_rerun/7x4x192x96x24';
+% resdir = 'paper_2_nonlinear/kT_6.96_rerun/9x5x192x96x24';
+% resdir = 'paper_2_nonlinear/kT_6.96/5x3x128x64x24';
+% resdir = 'paper_2_nonlinear/kT_6.96/7x4x128x64x24';
+% resdir = 'paper_2_nonlinear/2jm1_coeff_kT_6.96/5x3x192x96x24';
+% resdir = 'paper_2_nonlinear/2jm1_coeff_kT_6.96/7x4x192x96x24';
+% resdir = 'paper_2_nonlinear/2jm1_coeff_kT_6.96/9x5x192x96x24';
+
+resdir = 'paper_2_nonlinear/kT_6.96/2jm1_version_5x3x128x64x24';
+% resdir = 'paper_2_nonlinear/dbg/CBC_test_Napjmir_3x2x128x64x24';
+% resdir = 'paper_2_nonlinear/dbg/CBC_bench_3x2x128x64x24';
+
+% debug shearless
+% resdir = 'paper_2_nonlinear/shearless_CBC/7x4x192x96x24';
resdir = ['results/',resdir];
-JOBNUMMIN = 04; JOBNUMMAX = 10;
+JOBNUMMIN = 00; JOBNUMMAX = 10;
run analysis_gyacomo
diff --git a/wk/lin_3D_Zpinch.m b/wk/lin_3D_Zpinch.m
index a1de410d..b56220c7 100644
--- a/wk/lin_3D_Zpinch.m
+++ b/wk/lin_3D_Zpinch.m
@@ -95,8 +95,8 @@ MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
-GRADB = 1.0;
-CURVB = 1.0;
+k_gB = 1.0;
+k_cB = 1.0;
COLL_KCUT = 1000;
%%-------------------------------------------------------------------------
%% RUN
diff --git a/wk/lin_3Dzpinch.m b/wk/lin_3Dzpinch.m
index 95353bc9..be6be1d9 100644
--- a/wk/lin_3Dzpinch.m
+++ b/wk/lin_3Dzpinch.m
@@ -99,8 +99,8 @@ MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
-GRADB = 0.1;
-CURVB = 0.1;
+k_gB = 0.1;
+k_cB = 0.1;
COLL_KCUT = 1.5;
%%-------------------------------------------------------------------------
%% RUN
diff --git a/wk/lin_ETPY.m b/wk/lin_ETPY.m
index 8f25bcd5..fff23bff 100644
--- a/wk/lin_ETPY.m
+++ b/wk/lin_ETPY.m
@@ -100,8 +100,8 @@ MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
-GRADB = 1.0;
-CURVB = 1.0;
+k_gB = 1.0;
+k_cB = 1.0;
COLL_KCUT = 1000;
%%-------------------------------------------------------------------------
%% RUN
diff --git a/wk/lin_ITG.m b/wk/lin_ITG.m
index 2cb5704b..5390f576 100644
--- a/wk/lin_ITG.m
+++ b/wk/lin_ITG.m
@@ -10,12 +10,14 @@ addpath(genpath([gyacomodir,'matlab'])) % ... add
addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
+% SIMID = 'linear_CBC_benchmark_GX'; % Name of the simulation
SIMID = 'dbg'; % Name of the simulation
RUN = 1; % To run or just to load
default_plots_options
% EXECNAME = 'gyacomo_debug';
+% EXECNAME = 'gyacomo_2jm1';
+% EXECNAME = 'gyacomo_dlnBdz';
EXECNAME = 'gyacomo_alpha';
-% EXECNAME = 'gyacomo';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
@@ -33,15 +35,15 @@ KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 1e-4; % electron plasma beta
%% GRID PARAMETERS
P = 4;
-J = P/2;
+J = 2;%P/2;
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
NX = 8; % real space x-gridpoints
-NY = 2; % '' y-gridpoints
+NY = 12; % '' y-gridpoints
LX = 2*pi/0.1; % Size of the squared frequency domain
-LY = 2*pi/0.3; % Size of the squared frequency domain
+LY = 2*pi/0.1; % Size of the squared frequency domain
NZ = 24; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
@@ -59,13 +61,13 @@ PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
SHIFT_Y = 0.0;
%% TIME PARMETERS
-TMAX = 50; % Maximal time unit
-DT = 1e-2; % Time step
+TMAX = 40; % Maximal time unit
+DT = 2e-2; % Time step
% DT = 5e-4; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = -1; % Sampling per time unit for 2D arrays
SPS3D = 1; % Sampling per time unit for 2D arrays
-SPS5D = 1/2; % Sampling per time unit for 5D arrays
+SPS5D = 1/20; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
%% OPTIONS
@@ -96,14 +98,14 @@ INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
MU_Y = MU; %
N_HD = 4;
-MU_Z = 0.2; %
+MU_Z = 2.0; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
-GRADB = 1.0;
-CURVB = 1.0;
+k_gB = 1.0;
+k_cB = 1.0;
COLL_KCUT = 1000;
%%-------------------------------------------------------------------------
%% RUN
@@ -151,24 +153,6 @@ options.normalized = 1;
fig = plot_ballooning(data,options);
end
-if 0
-%% Hermite-Laguerre spectrum
-% options.TIME = 'avg';
-options.P2J = 1;
-options.ST = 1;
-options.PLOT_TYPE = 'space-time';
-% options.PLOT_TYPE = 'Tavg-1D';ls
-
-% options.PLOT_TYPE = 'Tavg-2D';
-options.NORMALIZED = 1;
-options.JOBNUM = 0;
-options.TIME = [0 50];
-options.specie = 'i';
-options.compz = 'avg';
-fig = show_moments_spectrum(data,options);
-% fig = show_napjz(data,options);
-save_figure(data,fig)
-end
if 0
%% linear growth rate for 3D Zpinch (kz fourier transform)
@@ -182,16 +166,16 @@ save_figure(data,fig)
end
if 1
%% Mode evolution
-options.NORMALIZED = 0;
-options.K2PLOT = 1;
-options.TIME = [0:1000];
-% options.TIME = [0.5:0.01:1].*data.Ts3D(end);
+options.NORMALIZED = 1;
+options.TIME = [000:9000];
+options.KX_TW = [0.5 1]*data.Ts3D(end); %kx Growth rate time window
+options.KY_TW = [0.5 1]*data.Ts3D(end); %ky Growth rate time window
options.NMA = 1;
-options.NMODES = 32;
-options.iz = 'avg';
-options.ik = 1;
+options.NMODES = 800;
+options.iz = 'avg'; % avg or index
+options.ik = 1; % sum, max or index
+options.fftz.flag = 0;
fig = mode_growth_meter(data,options);
-save_figure(data,fig,'.png')
end
diff --git a/wk/lin_KBM.m b/wk/lin_KBM.m
index abeafa21..f4f03f62 100644
--- a/wk/lin_KBM.m
+++ b/wk/lin_KBM.m
@@ -91,8 +91,8 @@ MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 0.0e-5; % Init noise amplitude
BCKGD0 = 1.0; % Init background
-GRADB = 1.0;
-CURVB = 1.0;
+k_gB = 1.0;
+k_cB = 1.0;
%%-------------------------------------------------------------------------
%% RUN
setup
diff --git a/wk/lin_MTM.m b/wk/lin_MTM.m
index fb45a8a6..3330f265 100644
--- a/wk/lin_MTM.m
+++ b/wk/lin_MTM.m
@@ -91,8 +91,8 @@ MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 0.0e-5; % Init noise amplitude
BCKGD0 = 1.0; % Init background
-GRADB = 1.0;
-CURVB = 1.0;
+k_gB = 1.0;
+k_cB = 1.0;
%%-------------------------------------------------------------------------
%% RUN
setup
diff --git a/wk/lin_RHT.m b/wk/lin_RHT.m
index d9a23009..7c43f8fd 100644
--- a/wk/lin_RHT.m
+++ b/wk/lin_RHT.m
@@ -92,8 +92,8 @@ MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 0.0e-5; % Init noise amplitude
BCKGD0 = 1.0; % Init background
-GRADB = 1.0;
-CURVB = 1.0;
+k_gB = 1.0;
+k_cB = 1.0;
COLL_KCUT = 1000;
%%-------------------------------------------------------------------------
%% RUN
diff --git a/wk/lin_TEM.m b/wk/lin_TEM.m
index 35eab68d..70e5f26a 100644
--- a/wk/lin_TEM.m
+++ b/wk/lin_TEM.m
@@ -91,8 +91,8 @@ MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 0.0e-5; % Init noise amplitude
BCKGD0 = 1.0; % Init background
-GRADB = 1.0;
-CURVB = 1.0;
+k_gB = 1.0;
+k_cB = 1.0;
%%-------------------------------------------------------------------------
%% RUN
setup
diff --git a/wk/local_run.m b/wk/local_run.m
index f2581ae2..7031cacb 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -23,8 +23,8 @@ J = 2;
Q0 = 1.0; % safety factor
SHEAR = 0.0; % magnetic shear
EPS = 0.0; % inverse aspect ratio
-GRADB = 1.0; % Magnetic gradient
-CURVB = 1.0; % Magnetic curvature
+k_gB = 1.0; % Magnetic gradient
+k_cB = 1.0; % Magnetic curvature
SG = 0; % Staggered z grids option
%% TIME PARAMETERS
TMAX = 120; % Maximal time unit
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 5e321406..4f78af0a 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -66,8 +66,8 @@ W_NAPJ = 1; W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% unused
KIN_E = 1; % Kinetic (1) or adiabatic (2) electron model
-GRADB = 1.0; % Magnetic gradient
-CURVB = 1.0; % Magnetic curvature
+k_gB = 1.0; % Magnetic gradient
+k_cB = 1.0; % Magnetic curvature
SG = 0; % Staggered z grids option
PMAXE = P; % Highest electron Hermite polynomial degree
JMAXE = J; % Highest '' Laguerre ''
diff --git a/wk/p2_heatflux_salpha_convergence.m b/wk/p2_heatflux_salpha_convergence.m
index 49773c32..092f776f 100644
--- a/wk/p2_heatflux_salpha_convergence.m
+++ b/wk/p2_heatflux_salpha_convergence.m
@@ -1,27 +1,27 @@
%% Heat flux convergence for kt=6.96 and 5.3
figure
title('s-$\alpha$ turb. heat flux conv.');
-% KT 6.96, nuDGDK = 0.05, 128x64x16, Nexc 1
-P = [2 4 12];
-Qx = [50.00 46.00 41.00];
-std_= [6.600 2.300 6.600];
- errorbar(P,Qx,std_/2,'o-.r',...
+% KT 6.96, nuDGDK = 0.05, 128x64x24, Nexc 5
+P = [2 4 12];
+Qx = [52.4 46.1 0];
+std_= [12.1 4.98 0];
+ errorbar(P,Qx,std_/2,'s-r',...
'LineWidth',2.0,'MarkerSize',8,...
'DisplayName','KT 6.9, nuDGDK 0.05'); hold on
xlabel('$P$, $J=P/2$'); ylabel('$Q_x$');
-% KT 6.96, nuDGDK = 0.05, 128x64x24, Nexc 5
-P = [4 6 8 10 ];
-Qx = [67.62 67.50 59.21 64.17];
-std_= [15.42 20.32 17.25 16.05];
+% KT 6.96, nuDGDK = 0.05, 192x96x24, Nexc 5
+P = [4 6 8 10 ];
+Qx = [42.9 44.9 00.0 00.0];
+std_= [9.15 7.12 0.00 0.00];
errorbar(P,Qx,std_/2,'o-r',...
'LineWidth',2.0,'MarkerSize',8,...
'DisplayName','KT 6.9, nuDGDK 0.05'); hold on
xlabel('$P$, $J=P/2$'); ylabel('$Q_x$');
% KT 5.3, nuDGDK = 0.05, 128x64x24, Nexc 5
P = [4 6 8 10 12 ];
-Qx = [44.10 21.61 16.04 0.558 0.901];
-std_= [10.61 6.952 4.166 0.025 0.122];
-errorbar(P,Qx,std_/2,'o-b',...
+Qx = [0 0 0 0 0];
+std_= [0 0 0 0 0];
+errorbar(P,Qx,std_/2,'o--r',...
'LineWidth',2.0,'MarkerSize',8,...
'DisplayName','KT 5.3, nuDGDK 0.05');
xlabel('$P$, $J=P/2$'); ylabel('$Q_x$');
@@ -30,7 +30,7 @@ xlabel('$P$, $J=P/2$'); ylabel('$Q_x$');
Nvp = [32 16 8 8];
Qx = [34.53 37.96 1.948 13.0];
std_= [7.830 6.048 0.629 3.09];
-errorbar(Nvp,Qx,std_/2,'s--r',...
+errorbar(Nvp,Qx,std_/2,'o-b',...
'LineWidth',2.0,'MarkerSize',8,...
'DisplayName','KT 6.9 GENE');
xlabel('$P=N_{v\parallel}$, $J=P/2=N_\mu$'); ylabel('$Q_x$');
@@ -42,11 +42,16 @@ xlabel('$P=N_{v\parallel}$, $J=P/2=N_\mu$'); ylabel('$Q_x$');
Nvp = [32 16 8];
Qx = [0.284 0.000 0.370];
std_= [0.177 0.000 0.140];
-errorbar(Nvp,Qx,std_/2,'s--b',...
+errorbar(Nvp,Qx,std_/2,'o--b',...
'LineWidth',2.0,'MarkerSize',8,...
'DisplayName','KT 5.3 GENE');
xlabel('$P=N_{v\parallel}$, $J=P/2=N_\mu$'); ylabel('$Q_x$');
+% DIMITS
+Dimits_result = 7.67*2*2.5;
+plot([0 32], Dimits_result*[1 1],'-.k','DisplayName','Dimits CBC');
+legend('show');
+
%% KT scans 9x5x128x64x24
clrs = lines(10);
figure
@@ -57,3 +62,37 @@ errorbar(kT_,Qx_,std_/2,'s--','color',clrs(4,:),...
'LineWidth',2.0,'MarkerSize',8,...
'DisplayName','(8,4)');
xlabel('$K_T$'); ylabel('$Q_x$');
+
+%% Add Dimits results on current plot
+plot([0 500], Dimits_result*[1 1],'-.k');
+%% Add Mandell Miller results on current plot
+Mandell_result = 7.67*7.5;
+plot([0 500], Mandell_result*[1.3 1.3],'-.k');
+plot([0 500], Mandell_result*[1 1],'--k');
+plot([0 500], Mandell_result*[0.7 0.7],'-.k');
+
+%% Old results (before the 2j-1 factor in mirror term)
+% KT 6.96, nuDGDK = 0.05, 128x64x16, Nexc 1
+% P = [2 4 12];
+% Qx = [50.00 46.00 41.00];
+% std_= [6.600 2.300 6.600];
+% errorbar(P,Qx,std_/2,'s-r',...
+% 'LineWidth',2.0,'MarkerSize',8,...
+% 'DisplayName','KT 6.9, nuDGDK 0.05'); hold on
+% xlabel('$P$, $J=P/2$'); ylabel('$Q_x$');
+% % KT 6.96, nuDGDK = 0.05, 128x64x24, Nexc 5
+% P = [4 6 8 10 ];
+% Qx = [67.62 67.50 59.21 64.17];
+% std_= [15.42 20.32 17.25 16.05];
+% errorbar(P,Qx,std_/2,'o-r',...
+% 'LineWidth',2.0,'MarkerSize',8,...
+% 'DisplayName','KT 6.9, nuDGDK 0.05'); hold on
+% xlabel('$P$, $J=P/2$'); ylabel('$Q_x$');
+% % KT 5.3, nuDGDK = 0.05, 128x64x24, Nexc 5
+% P = [4 6 8 10 12 ];
+% Qx = [44.10 21.61 16.04 0.558 0.901];
+% std_= [10.61 6.952 4.166 0.025 0.122];
+% errorbar(P,Qx,std_/2,'o--r',...
+% 'LineWidth',2.0,'MarkerSize',8,...
+% 'DisplayName','KT 5.3, nuDGDK 0.05');
+% xlabel('$P$, $J=P/2$'); ylabel('$Q_x$');
\ No newline at end of file
diff --git a/wk/quick_run.m b/wk/quick_run.m
index aff73ce3..dffb3b2a 100644
--- a/wk/quick_run.m
+++ b/wk/quick_run.m
@@ -91,8 +91,8 @@ MU_J = 0.0; %
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
-GRADB = 1.0;
-CURVB = 1.0;
+k_gB = 1.0;
+k_cB = 1.0;
%%-------------------------------------------------------------------------
%% RUN
setup
--
GitLab