diff --git a/matlab/compute/compute_coll.m b/matlab/compute/compute_coll.m new file mode 100644 index 0000000000000000000000000000000000000000..c0f182ed033f1504772240b48f27f5a0fb6ea39b --- /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 0000000000000000000000000000000000000000..9ef8818b18b9c9f2dc4a2c4394459992f4458f52 --- /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 d449552ed0bc09b73e36b63f6241864e318e4abc..9c53fe171a521cd2cf8f4642ec3ae7f9a89128a4 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 5d8f47a4c4d39ce199529ce938d989fb5a81bc0e..f43866e82dd9cc6cf5cb02dad988bbd859b9f083 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 39d3855107422f1215465ced486610a4c5fc94a5..1267999b25d0121511e955858f9b6b53c6509e4f 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 06290b4c8af281f769afe8aad7ff0da7f2b91bd9..9390b1b2cd33e40186f458a292373df32dad64bd 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 6376d03511eb590356d6d43c633f33b129535e5b..2a8f461ecf6d421a7d4e83bf4cd77f3d344c3765 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