From 3f97b0a292012a892eb9726dcd3520aa5788543a Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Fri, 24 Nov 2023 15:19:39 +0100
Subject: [PATCH] scripts update and corrections
---
matlab/compute/compute_coll.m | 9 +++++
matlab/compute/compute_growth.m | 15 ++++++++
matlab/compute/mode_growth_meter.m | 3 +-
matlab/plot/plot_ballooning.m | 1 +
matlab/plot/plot_metric.m | 57 +++++++++++++++++++++++++++---
matlab/process_field.m | 20 +++++------
wk/get_param_from_profiles.m | 3 ++
7 files changed, 91 insertions(+), 17 deletions(-)
create mode 100644 matlab/compute/compute_coll.m
create mode 100644 matlab/compute/compute_growth.m
diff --git a/matlab/compute/compute_coll.m b/matlab/compute/compute_coll.m
new file mode 100644
index 0000000..c0f182e
--- /dev/null
+++ b/matlab/compute/compute_coll.m
@@ -0,0 +1,9 @@
+function [nuGENE, nuGYAC] = compute_coll(Lref,n_e,T_e,T_i)
+
+ tau = T_i/T_e;
+
+ nuGENE = 2.3031E-5*Lref*(n_e)/(T_e)^2*(24.-log(sqrt(n_e*1.0E13)/T_e*0.001));
+
+ nuGYAC = 3/8*sqrt(pi/2)*tau.^(3/2)*nuGENE;
+
+end
\ No newline at end of file
diff --git a/matlab/compute/compute_growth.m b/matlab/compute/compute_growth.m
new file mode 100644
index 0000000..9ef8818
--- /dev/null
+++ b/matlab/compute/compute_growth.m
@@ -0,0 +1,15 @@
+function [gamma,err] = compute_growth(t,y)
+gamma = zeros(size(y));
+ for it = 2:numel(t)
+ y_n = y(it);
+ y_nm1 = y(it-1);
+ dt = t(it)-t(it-1);
+ ZS = sum(y_nm1,1);
+ wl = log(y_n./y_nm1)/dt;
+ gamma(it) = squeeze(sum(wl.*y_nm1,1)./ZS);
+ end
+ % error estimation
+ n = 5;
+ [gamma, ~, err] = sliceAverage(gamma, n);
+ err = mean(err);
+end
\ No newline at end of file
diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index d449552..9c53fe1 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -1,4 +1,4 @@
-function [FIGURE, kykx, wkykx, ekykx] = mode_growth_meter(DATA,OPTIONS)
+function [FIGURE, wkykx, ekykx] = mode_growth_meter(DATA,OPTIONS)
NORMALIZED = OPTIONS.NORMALIZED;
Nma = OPTIONS.NMA; %Number moving average
@@ -41,7 +41,6 @@ TW = [OPTIONS.KY_TW; OPTIONS.KX_TW];
[~,it1] = min(abs(t-TW(1,1)));
[~,it2] = min(abs(t-TW(1,2)));
[wkykx, ekykx] = compute_growth_rates(FIELD(:,:,:,it1:it2),DATA.Ts3D(it1:it2));
-kykx = meshgrid(DATA.grids.ky,DATA.grids.kx)';
FIGURE = struct();
if OPTIONS.SHOWFIG
diff --git a/matlab/plot/plot_ballooning.m b/matlab/plot/plot_ballooning.m
index 5d8f47a..f43866e 100644
--- a/matlab/plot/plot_ballooning.m
+++ b/matlab/plot/plot_ballooning.m
@@ -1,5 +1,6 @@
function [FIG, b_angle, phib, psib, ky_] = plot_ballooning(DATA,OPTIONS)
FIG.fig = figure; iplot = 1;
+ psib = 0;
[~,it0] = min(abs(DATA.Ts3D - OPTIONS.time_2_plot(1)));
[~,it1] = min(abs(DATA.Ts3D - OPTIONS.time_2_plot(end)));
[~,ikyarray] = min(abs(DATA.grids.ky - OPTIONS.kymodes));
diff --git a/matlab/plot/plot_metric.m b/matlab/plot/plot_metric.m
index 39d3855..1267999 100644
--- a/matlab/plot/plot_metric.m
+++ b/matlab/plot/plot_metric.m
@@ -1,4 +1,4 @@
-function [ fig, arrays, R, Z ] = plot_metric( data, options )
+function [ fig, arrays, R, Z] = plot_metric( data, options )
fig = 0;
% names = {'Jacobian','gradxB','gradyB','gradzB','gradz_coeff',...
% 'gxx','gxy','gxz','gyy','gyz','gzz','hatB','hatR','hatZ'};
@@ -11,12 +11,25 @@ for i_ = 1:numel(names)
namae = names{i_};
geo_arrays(:,i_) = h5read(data.outfilenames{end},['/data/metric/',namae])';
end
+gxx = geo_arrays(:,1);
+gxy = geo_arrays(:,2);
+gxz = geo_arrays(:,3);
+gyy = geo_arrays(:,4);
+gyz = geo_arrays(:,5);
+% gzz = geo_arrays(:,6);
+B = geo_arrays(:,7);
+dxB = geo_arrays(:,8);
+dyB = geo_arrays(:,9);
+dzB = geo_arrays(:,10);
+Jac = geo_arrays(:,11);
+
try
- NPLOT = options.SHOW_FLUXSURF + options.SHOW_METRICS;
+ NPLOT = options.SHOW_FLUXSURF + options.SHOW_METRICS + options.SHOW_CURVOP;
catch
NPLOT = 2;
options.SHOW_FLUXSURF = 1;
- options.SHOW_METRICS = 1;
+ options.SHOW_METRICS = 0;
+ options.SHOW_CURVOP = 1;
end
if NPLOT > 0
fig = figure;
@@ -41,12 +54,46 @@ if NPLOT > 0
legend('Interpreter','none')
end
if options.SHOW_FLUXSURF
- subplot(1,2,1)
+ subplot(1,NPLOT,1)
R = [geo_arrays(:,12);geo_arrays(1,12)];
Z = [geo_arrays(:,13);geo_arrays(1,13)];
- plot(R,Z);
+ plot(R./data.inputs.EPS,Z./data.inputs.EPS,'-k');
+ xlabel('R [m]'); ylabel('Z [m]');
+ axis tight
axis equal
end
+ if options.SHOW_CURVOP
+ G1 = gxx.*gyy - gxy.*gxy;
+ G2 = gxx.*gyz - gxy.*gxz;
+ G3 = gxy.*gyz - gyy.*gxz;
+ Ckx0 = -(dyB./B + G2./G1.*dzB./B);
+ C0ky = (dxB./B - G3./G1.*dzB./B);
+ % subplot(1,NPLOT,2);
+ % plot(data.grids.z,Ckx0,'DisplayName','$\mathcal C(k_x,k_y=0)/ik_x$'); hold on
+ % plot(data.grids.z,C0ky,'DisplayName','$\mathcal C(k_x=0,k_y)/ik_y$');
+ % xlabel('$z$'); legend("show");
+ subplot(1,NPLOT,2);
+ % z = [data.grids.z; data.grids.z(1)]+pi;
+ z = data.grids.z;
+ % rho = 1+0.5*(Ckx0)./max(abs(Ckx0));
+ rho = Ckx0;
+ % rho = [rho; rho(1)];
+ minr = min(rho); maxr = max(rho);
+ plot(z,rho,'DisplayName','$\mathcal C(k_x,k_y=0)/ik_x$'); hold on
+ % rho = 1+0.5*(C0ky)./max(abs(C0ky));
+ rho = C0ky;
+ % rho = [rho; rho(1)];
+ minr = min(minr,min(rho)); maxr = max(maxr,max(rho));
+ plot(z,rho,'DisplayName','$\mathcal C(k_x=0,k_y)/ik_y$');
+ % rho = 1+0.5*(B.^(-1)./Jac)./max(abs(B.^(-1)./Jac));
+ rho = B.^(-1)./Jac;
+ % rho = [rho; rho(1)];
+ minr = min(minr,min(rho)); maxr = max(maxr,max(rho));
+ plot(z,rho,'DisplayName','$\hat B^{-1}/J_{xyz}$');
+ % polarplot(z,zeros(size(rho)),'-k','DisplayName','$0$')
+ legend("show",'Location','south','FontSize',14);
+ % rlim(1.1*[minr maxr]);
+ end
end
%outputs
arrays = squeeze(geo_arrays(:,:));
diff --git a/matlab/process_field.m b/matlab/process_field.m
index 06290b4..9390b1b 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -230,31 +230,31 @@ switch OPTIONS.NAME
FLD_ = 0.*DATA.PHI(:,:,:,FRAMES);
OPE_ = 1;
for it = 1:numel(FRAMES)
- tmp = zeros(DATA.Ny,DATA.Nx,Nz);
- for iz = 1:DATA.Nz
+ tmp = zeros(DATA.grids.Ny,DATA.grids.Nx,Nz);
+ for iz = 1:DATA.grids.Nz
vx_ = real((ifourier_GENE(-1i*KY.*(DATA.PHI (:,:,iz,FRAMES(it))))));
ni_ = real((ifourier_GENE( DATA.DENS_I(:,:,iz,FRAMES(it)))));
gx_ = vx_.*ni_;
-% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.Ny,DATA.Nx))),2));
- tmp(:,:,iz) = abs((squeeze(fft2(gx_,DATA.Ny,DATA.Nx))));
+% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.grids.Ny,DATA.grids.Nx))),2));
+ tmp(:,:,iz) = abs((squeeze(fft2(gx_,DATA.grids.Ny,DATA.grids.Nx))));
end
FLD_(:,:,it)= squeeze(compr(tmp(1:DATA.grids.Nky,1:DATA.grids.Nkx,:)));
end
case 'Q_x' % ion heat flux
NAME = 'Qx';
- FLD_ = 0.*DATA.PHI(:,:,:,FRAMES);
+ % FLD_ = 0.*DATA.PHI(:,:,:,FRAMES);
OPE_ = 1;
for it = 1:numel(FRAMES)
- tmp = zeros(DATA.Ny,DATA.Nx,Nz);
- for iz = 1:DATA.Nz
+ tmp = zeros(DATA.grids.Ny,DATA.grids.Nx,Nz);
+ for iz = 1:DATA.grids.Nz
vx_ = real((ifourier_GENE(-1i*KY.*(DATA.PHI (:,:,iz,FRAMES(it))))));
ni_ = real((ifourier_GENE( DATA.DENS_I(:,:,iz,FRAMES(it)))));
Ti_ = real((ifourier_GENE( DATA.TEMP_I(:,:,iz,FRAMES(it)))));
qx_ = vx_.*ni_.*Ti_;
-% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.Ny,DATA.Nx))),2));
- tmp(:,:,iz) = abs((squeeze(fft2(qx_,DATA.Ny,DATA.Nx))));
+% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.grids.Ny,DATA.grids.Nx))),2));
+ tmp(:,:,iz) = abs((squeeze(fft2(qx_,DATA.grids.Ny,DATA.grids.Nx))));
end
- FLD_(:,:,it)= squeeze(compr(tmp(1:DATA.grids.Nky,1:DATA.grids.Nkx,:)));
+ FLD_(:,:,:,it)= squeeze(tmp(1:DATA.grids.Nky,1:DATA.grids.Nkx,:));
end
case 'f_i'
SKIP_COMP = 1;
diff --git a/wk/get_param_from_profiles.m b/wk/get_param_from_profiles.m
index 6376d03..2a8f461 100644
--- a/wk/get_param_from_profiles.m
+++ b/wk/get_param_from_profiles.m
@@ -35,4 +35,7 @@ end
params.BETA = 403.0E-5*n_e*T_e/(Bref*Bref);
cref= sqrt(T_e*1000*1.6e-19/mref);
params.EXBRATE= K_w*wExB*1000*Lref/cref;
+ params.ne = n_e;
+ params.Te = T_e;
+ params.Ti = T_i;
end
\ No newline at end of file
--
GitLab