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

scripts update and corrections

parent 3627a023
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
matlab/compute/compute_coll.m 0 → 100644
+ 9
− 0
View file @ 3f97b0a2
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
matlab/compute/compute_growth.m 0 → 100644
+ 15
− 0
View file @ 3f97b0a2
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
matlab/compute/mode_growth_meter.m
+ 1
− 2
View file @ 3f97b0a2
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
......
matlab/plot/plot_ballooning.m
+ 1
− 0
View file @ 3f97b0a2
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));
......
matlab/plot/plot_metric.m
+ 52
− 5
View file @ 3f97b0a2
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(:,:));
......
matlab/process_field.m
+ 10
− 10
View file @ 3f97b0a2
......@@ -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;
......
wk/get_param_from_profiles.m
+ 3
− 0
View file @ 3f97b0a2
......@@ -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
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